% crank_rocker_design.m 
%
% This program computes the link lengths of crank rocker mechanism.
% The input variables are the oscillation angle, the time ratio, 
% and the length of one of the links.  

%Variables
%   Variables for calculation 
%   r1=frame length (cm)
%   r2=crank length (cm)
%   r3=coupler length (cm)
%   r4=rocker length (cm)
%   r = scaled link lengths [r1 r2 r3 r4]
%   R = actual link lengths [R1 R2 R3 R4]
%   theta=rocker oscillation angle (deg)
%   thetar=rocker oscillation angle (rad)
%   Bx, By=crank pin location - x and y coordinates
%   Cx, Cy=rocker pin location - x and y coordinates
%   Q22=crank angle throughout animation
%   Q44=rocker angle throughout animation
%   mode=assembly mode
%   circleb,circleb=x,y coordinates of the B1 circle
%   radi=distance based on frame and crank link lengths
%   t3,t4=limiting rocker angles (rad)
%   bushx1,bushy1=x,y coordinates of the first bushing
%   bushx2,bushy2=x,y coordinates of the second bushing
%   ipos=number of linkage positions for the animation
%   fact=conversion factor, from degrees to radians
%   ans=string variable controlling animation repetition
%   lknown = link with the known length
%   rknown = known link length

%   Graphics Variables

%   joint1=line for joint at crank/coupler
%   joint2=line for joint at coupler/rocker
%   crank=line for the crank
%   coupler=line for the coupler
%   rocker=line for the rocker
%   bush1,bush2=lines representing the bushings
%   xmin, xmax= x limits for the axes
%   ymin, ymax= y limits for the axes

% program 

clear all;
clf;
ans='y';
ipos=40;  
disp('   ')
disp('      Crank-Rocker Synthesis Program')
disp('  ')

% Input values for linkage variables 

lknown=input('Enter link number corresponding to known length [1]:  ');
if isempty(lknown); lknown=1; end
rknown=input('Enter known link length [2]: ');
if isempty(rknown);rknown = 2; end  
ansinput='y';
while ansinput=='y'
	theta=5;
	while theta < 10
		theta=input('Enter rocker oscillation angle (deg) [80]: ');
		if isempty(theta);theta=80; end  
		if theta < 0.001
	  	  disp ('Theta too small; reinput ')
		end
	end
	
% Enter time ratio or alpha

	which=input('Enter time ratio (Q) or alpha (a)? [a]:  ','s');
	if isempty(which); which ='a'; end
	if which=='Q'
		Q=input('Enter time ratio for mechanism [1.32258065]: ');
		if isempty(Q);Q=1.32258065; end
		if Q < 0; Q=0; end
		alpha=180*(Q-1)/(Q+1);
	end
	if which=='a'
		alpha=input('Enter alpha angle (in deg) for mechanism [10]: ');
		if isempty(alpha); alpha = 10; end
		Q=(180+alpha)/(180-alpha);
	end
	ansi=input('reenter input variables? y/n [n]: ', 's');
	if isempty(ansi);
		ansi='n';
	end

% Print heading for input angles

	disp('   ')
	fprintf('    theta (deg)   Q    alpha (deg)') 

% Print the input values.

	disp('   ')
	fprintf('%10.3f %10.3f %10.3f\n',theta, Q, alpha)
	disp('   ')

% Convert angles to radians

	fact=pi/180;
	thetar=theta*fact;
	alphar=alpha*fact;

% Locate the radius and center of the B2 circle.
	
	denom=sin(pi-thetar+alphar);
	rb=100000;
	sliderflag=1;
	if denom~=0
		sliderflag=0;
		rb=sin(thetar/2)/denom;
	end
	xg=rb*cos(thetar/2 - alphar);
	yg=-rb*sin(thetar/2 - alphar);

% Find the limiting angles for the B2 circle;

	betamax=pi+alphar-thetar/2;
	ym=0;
	xm=2*rb*cos(thetar/2-alphar);
	if alpha >=0
		tt=tan(thetar);
		if abs(tt) > 0.00000001 & abs(tt) < 1000000
			A=1+1/(tan(thetar))^2;
			B=-2*((xg-1)/tan(thetar) + yg);
		end
		if abs(tt) > 1000000;
			A=1;
			B=-2*yg;
		end
		C=(1+yg^2+xg^2-2*xg-rb^2);
		ym=100000;
		arg=B*B-4*A*C;
		if arg > 0
			if A ~= 0
				ym1=(-B+sqrt(arg)/(2*A));
				ym2=(-B-sqrt(arg)/(2*A));
				ym=max([ym1, ym2]);

% If ym is negative for both roots, set ym = 0,

				if ym<0; ym=0; end
			end
		end
		if arg >= 0; 
			ansinput='n'; 
		end
		if arg < 0
			disp ('No solution exists; change input data ');
		end
		xm = 1+ym/tan(thetar);
	end
end
betamin=atan2((ym-yg), (xm-xg));
	
% Set up coordinates of bushings and pins.

rpivot=0.03;
ninc=20;  
npoints=ninc+1;
centerx=0;
centery=0;
coord=circle(rpivot,centerx,centery,ninc);

% find coordinates of first pin

for i=1:1:npoints
	xcoord1(i)=coord(i,1);
	xcoord2(i)=xcoord1(i)+1;
	ycoord1(i)=coord(i,2);
end
ycoord2=ycoord1;

% Find coordinates of first bushing

binc=ninc/2;
coord1=bushing(rpivot,centerx,centery,binc,0);
bpoints=ninc/2+16;
for i=1:1:bpoints
	xbush1(i)=coord1(i,1);
	xbush2(i)=xbush1(i)+1;		
	ybush1(i)=coord1(i,2);
end
ybush2=ybush1;

% Set up coordinates for B2 arc if rb is a reasonable size.

ninc=180;
if rb < 20
	inc=(betamax-betamin)/ninc;
	betamean=(betamax+betamin)/2;
	npoints = ninc+1;
	i=0;
	for j=betamin:inc:betamax
		i=i+1;
		coordx(i)=xg+rb*cos(j);
		coordy(i)=yg+rb*sin(j);
	end

% Set up coordinates to end line to B2 arc.

	xline1=[1, xg+rb*cos(betamin)];
	yline1=[0, yg+rb*sin(betamin)];

end

% Define B2 "arc" when rb is large

if rb >= 20
	coordx(1)=0;
	coordy(1)=0;
	angB2=alphar-thetar/2+pi/2;
	coordx(2)=2*cos(angB2);
	coordy(2)=2*sin(angB2);
	xline1=[0 0];
	yline1=[0 0];
end

% Set up the axes for first drawing.

ansdesign='y';
while ansdesign=='y'
	clf;
	width=0.85;
	height=0.75; 
	h1=axes('position',[0.1 0.1 width height],'box','on','xcolor','k',...
	        'ycolor','k'); 
	set(gcf,'color', 'w')

	cx=[xcoord1, xbush1, xcoord2, xbush2,1,0,coordx,xline1,];
	cy=[ycoord1,ybush1,ycoord2, ybush2,0,coordy,yline1,];
	xmina=min(cx);
	xmaxa=max(cx);
	ymina=min(cy);
	if ymina>-.5; ymina=-0.5; end
	ymaxa=max(cy);
	rangex=xmaxa-xmina;
	rangey=ymaxa-ymina;
	xmin=xmina-0.4*rangex;
	xmax=xmaxa+0.4*rangex;
	ymin=ymina-0.1*rangey;
	ymax=ymaxa+0.6*rangey;
	values=axisadjust(xmin, xmax, ymin, ymax, width, height);
	xmin=values(1);
	xmax=values(2);
	ymin=values(3);
	ymax=values(4);
	axis equal
	axis([xmin xmax ymin ymax]);

% Set up vectors to be drawn

	joint1=line('xdata', [],'ydata', [], 'marker', 'o','markersize', ...
      6,'erasemode', 'xor', 'color', 'r');
	joint2=line('xdata', [], 'ydata', [], 'marker','o','markersize', ...
      6, 'erasemode', 'xor', 'color','r');
	joint12=line('xdata', [],'ydata', [], 'marker', 'o','markersize', ...
      6,'erasemode', 'xor', 'color', 'k');
	joint22=line('xdata', [], 'ydata', [], 'marker','o','markersize', ...
      6, 'erasemode', 'xor', 'color','k');
	joint13=line('xdata', [],'ydata', [], 'marker', 'o','markersize', ...
      6,'erasemode', 'xor', 'color', 'k');
	joint23=line('xdata', [], 'ydata', [], 'marker','o','markersize', ...
      6, 'erasemode', 'xor', 'color','k');

	line1=line('xdata', [],'ydata', [], 'linestyle', '-','markersize', ...
      8,'erasemode', 'xor', 'color', 'b','linewidth',1);
	line2=line('xdata', [],'ydata', [], 'linestyle', '-','markersize', ...
      8,'erasemode', 'xor', 'color', 'k','linewidth',1);
	line3=line('xdata', [],'ydata', [], 'linestyle', '-','markersize', ...
      8,'erasemode', 'xor', 'color', 'k','linewidth',1);

	crank=line('xdata', [], 'ydata' ,[], 'linewidth' ,2, 'erasemode',...
      'xor','color', 'r');
	coupler=line('xdata', [], 'ydata' ,[], 'linewidth' ,2, 'erasemode',...
      'xor','color', 'r');
	rocker=line('xdata', [], 'ydata' ,[], 'linewidth',2,...
      'erasemode', 'xor','color', 'r');
	crank2=line('xdata', [], 'ydata' ,[], 'linewidth' ,2, 'erasemode',...
      'xor','color', 'k');
	coupler2=line('xdata', [], 'ydata' ,[], 'linewidth' ,2, 'erasemode',...
      'xor','color', 'k');
	rocker2=line('xdata', [], 'ydata' ,[], 'linewidth',2,...
      'erasemode', 'xor','color', 'k');
	crank3=line('xdata', [], 'ydata' ,[], 'linewidth' ,2, 'erasemode',...
      'xor','color', 'k');
	coupler3=line('xdata', [], 'ydata' ,[], 'linewidth' ,2, 'erasemode',...
      'xor','color', 'k');
	rocker3=line('xdata', [], 'ydata' ,[], 'linewidth',2,...
      'erasemode', 'xor','color', 'k');

	bush1=line('xdata',[],'ydata',[], 'erasemode',...
       'none','color', 'r');
	bush2=line('xdata',[],'ydata',[], 'erasemode',...
       'none','color', 'r');
	pin1=line('xdata',[],'ydata',[], 'erasemode',...
       'none','color', 'r');
	pin2=line('xdata',[],'ydata',[], 'erasemode',...
       'none','color', 'r');

	arc1=line('xdata',[],'ydata',[],'linestyle','-','markersize',...
      6,'erasemode','none', 'color','b');
	arc2=line('xdata',[],'ydata',[],'linestyle','-', 'markersize',...
      6,'erasemode','none', 'color','b');
		
% draw the two pins and bushings.

	set(pin1,'xdata', xcoord1, 'ydata',ycoord1);
	set(bush1,'xdata', xbush1,'ydata', ybush1);
	set(pin2,'xdata', xcoord2, 'ydata',ycoord2);
	set(bush2,'xdata', xbush2,'ydata', ybush2);

% Draw the B2 arc.

	set (arc1,'xdata', coordx, 'ydata', coordy)

% Draw line to end of B2 arc

	set (line1, 'xdata', xline1, 'ydata',yline1);

% Draw x at point G.

	set (line2, 'xdata', [xg-0.05, xg+0.05], 'ydata',[yg, yg])
	set (line3, 'xdata', [xg, xg], 'ydata',[yg-0.05, yg+0.05])

% Ask user for mode of input.	
	
	disp('   ')
	disp('Possible modes of inputting data: ')
	disp('  1 - mouse ')
	disp('  2 - keyboard ')
	typei=input('Designate mode of input [1]: ');
	if isempty(typei)
		typei = 1;
	end
	ans1='y';
		
% while loop to input B2.

	while ans1=='y'

% Loop for B2 input.  Input the point using mouse if itype1=1.

		if typei==1
			disp('   ')
			disp('  Use mouse to locate B2')
			[fx1,fy1]=ginput(1);

% Find angle based on coordinates.

			beta=atan2((fy1-yg),(fx1-xg));
		end

% Input the angle beta using keyboard if itype1=2.

		if typei==2
			if rb < 20
				disp('   ')
				beta=input(' Input beta (in deg) to locate B2: ');
				if isempty(beta); beta = (betamax-betamin)/(2*fact);end
				beta=beta*fact;
				if beta > betamax; beta = betamax-0.0001; end
				if beta < betamin; beta = betamin+0.0001; end
			end
			if rb >= 20
				disp('   ')
				fy1=input(' Input y coordinate of B2: ');
				if isempty(fy1); fy1 = 1;end
				rfy=fy1/sin(angB2);
				fx1=rfy*cos(angB2);
				beta=atan2((fy1-yg),(fx1-xg));
			end			
		end

% Compute the coordinates of B1 and B2.

		xb2=xg+rb*cos(beta);
		yb2=yg+rb*sin(beta);
		theta22=180+atan2(yb2,xb2)/fact;
		phi=atan2(yb2, (xb2-1));
		r4=sqrt((xb2-1)^2+yb2^2);
		r1=1;
		angb1=phi-thetar;
		xb1=1+r4*cos(angb1);
		yb1=r4*sin(angb1);
		theta21=atan2(yb1,xb1)/fact;

% Compute r2 and r3.

		leng1=sqrt(xb1^2+yb1^2);
		leng2=sqrt(xb2^2+yb2^2);
		r2=(leng1-leng2)/2;
		r3=(leng1+leng2)/2;

% Draw the linkage in an arbitrary position.

		flag=1;
		theta2=(theta21+theta22)/2;
		mode=-1;
		theta1=0;
		points = fourbar(r1,r2,r3,r4,theta2,0,0,mode,theta1,flag);
		Bxx=points(17);
		Byy=points(18);
		Cxx=points(19);
		Cyy=points(20);
		set(joint1,'xdata', Bxx,'ydata',Byy);
		set(joint2,'xdata', Cxx, 'ydata',Cyy);
		set(crank,'xdata',[0 Bxx],'ydata', [0 Byy]);
		set(coupler,'xdata', [Bxx Cxx],'ydata', [Byy Cyy]);
		set(rocker,'xdata', [r1 Cxx],'ydata', [0 Cyy]);

% Draw linkage in extreme (B1) position

		theta2=theta21;
		points = fourbar(r1,r2,r3,r4,theta2,0,0,mode,theta1,flag);
		Bxx=points(17);
		Byy=points(18);
		Cxx=points(19);
		Cyy=points(20);
		set(joint12,'xdata', Bxx,'ydata',Byy);
		set(joint22,'xdata', Cxx, 'ydata',Cyy);
		set(crank2,'xdata',[0 Bxx],'ydata', [0 Byy]);
		set(coupler2,'xdata', [Bxx Cxx],'ydata', [Byy Cyy]);
		set(rocker2,'xdata', [r1 Cxx],'ydata', [0 Cyy]);

% Draw linkage in extreme (B2) position

		theta2=theta22;
		points = fourbar(r1,r2,r3,r4,theta2,0,0,mode,theta1,flag);
		Bxx=points(17);
		Byy=points(18);
		Cxx=points(19);
		Cyy=points(20);
		set(joint13,'xdata', Bxx,'ydata',Byy);
		set(joint23,'xdata', Cxx, 'ydata',Cyy);
		set(crank3,'xdata',[0 Bxx],'ydata', [0 Byy]);
		set(coupler3,'xdata', [Bxx Cxx],'ydata', [Byy Cyy]);
		set(rocker3,'xdata', [r1 Cxx],'ydata', [0 Cyy]);
		ans1=input('reenter point? y/n [n]: ', 's');
		if isempty(ans1);
			ans1='n';
		end		
	end
	
% Unscale the link lengths

	ra=[r1 r2 r3 r4];
	scale = rknown/ra(lknown);
	R=ra*scale;

% Print out the link Lengths

	disp('   ')
	fprintf('     Scaled values for link lengths ') 
	disp('   ')
	fprintf('       r1         r2         r3         r4 ') 
	disp('  ')
	fprintf('%10.3f %10.3f %10.3f %10.3f\n',ra(1),ra(2),ra(3),ra(4))
	disp('   ')
	fprintf('       Actual values for link lengths ') 
	disp('   ')
	fprintf('       R1         R2         R3         R4 ') 
	disp('  ')
	fprintf('%10.3f %10.3f %10.3f %10.3f\n',R(1),R(2),R(3),R(4))

% **********************************************************

% The linkage is designed, now animate the results.

% **********************************************************

	pause;
	clf;
	tt=3;
	w2=1;
	Tout=1;
	Q11=0;
	Q1=0;

% Define the axes for the first graph (input/output angle graph)

	width=0.38;
	height=0.35;
	h1=axes('position',[0.6 .55 width height],'box','on','xcolor','k','ycolor',...
      'k');   %axis location
	set(gcf,'color', 'w')
	xlabel('crank angle','color', 'k')
	ylabel('rocker angle, deg','color', 'k')
	Displ=line('xdata', [], 'ydata', [], 'erasemode', 'none',...
      'color', 'r', 'linestyle', '-');
	dbead=line('xdata', [], 'ydata', [],'marker','o','markersize', 6, 'erasemode',...
      'xor');

% Define the axes for the second graph (rocker velocity graph)

	h2=axes('position', [0.6 .1 width height],'box', 'on','xcolor', 'k','ycolor',...
      'k');
	set(gcf,'color', 'w')
	xlabel('crank angle','color', 'k')
	ylabel('transmission angle, deg','color', 'k')
	Tangle=line('xdata', [], 'ydata', [], 'erasemode', 'none', ...
      'color', 'r', 'linestyle', '-');
	vbead=line('xdata', [], 'ydata', [],'marker','o','markersize', 6, ...
      'erasemode','xor');

% Define the axes for the third graph (torque graph)

	h3=axes('position', [0.1 .1 width height],'box', 'on','xcolor', 'k','ycolor',...
      'k');
	set(gcf,'color', 'w')
	xlabel('crank angle','color', 'k')
	ylabel('input torque/output torque','color', 'k')
	Torque=line('xdata', [], 'ydata', [], 'erasemode','none', ...
      'color', 'r', 'linestyle', '-');
	tbead=line('xdata', [], 'ydata', [],'marker','o','markersize', 6, 'erasemode','xor');

  
% Set up the vectors needed to draw the linkage

	h4=axes('position', [0.1 .55 width height],'box', 'on', 'xcolor', 'k','ycolor',...
       'k');
	set(gcf,'color', 'w')
	axis equal;
	axis([xmin xmax ymin ymax]);
	joint1=line('xdata', [], 'ydata', [], 'marker', 'o','markersize', ...
       6, 'erasemode', 'xor', 'color', 'k');
	joint2=line('xdata', [], 'ydata', [], 'marker', 'o','markersize', ...
       6, 'erasemode', 'xor', 'color', 'k');
	crank=line('xdata', [0 0], 'ydata' ,[], 'linewidth' ,2, 'erasemode', ...
       'xor','color', 'k');
	coupler=line('xdata', [], 'ydata' ,[], 'linewidth' ,2, 'erasemode', ...
       'xor','color', 'k');
	rocker=line('xdata', [], 'ydata' ,[], 'linewidth',2,...
       'erasemode', 'xor','color', 'k');
	bush1=line('xdata',[],'ydata',[], 'erasemode',...
       'none','color', 'r');
	bush2=line('xdata',[],'ydata',[], 'erasemode','none','color', 'r');
	pin1=line('xdata',[],'ydata',[], 'erasemode','none','color', 'r');
	pin2=line('xdata',[],'ydata',[], 'erasemode','none','color', 'r');

% Define constants

	i=0;
	r2d=180/pi;
	mode=-1;
	dt=5;
 
% Loop to analyze the linkage for tt full cycles

	for Q21 = 0:dt:360*tt;
		i=i+1;
		Q2 = Q21*fact;

% Call function to compute the position

		para = fourbar(r1,r2,r3,r4,Q21,w2,0,mode,Q1,2); 
		Q3=para(7);
		Q33(i)=Q3;
		Q4=para(8);
		Q44(i)=Q4;
		w33(i)=para(11);
		w44=para(12);
		tang(i)=Q44(i)-Q33(i);
		Tin1(i)=(Tout*w44)/w2;
		assem=para(29);

% Store position information in matrices.

		Bx(i)= para(17);
		By(i)= para(18);
		Cx(i)= para(19);
		Cy(i)= para(20);
		Q22(i) = Q21;
		w22(i) = w2; 
	end	

% Set the axis limits for the graphs.  Start with the fourth
% graph (linkage position)
		 
	tempx= [xcoord1,xbush1,xcoord2,xbush2,Bx, Cx];
	xmina=min(tempx);
	xmaxa=max(tempx);
	tempy=[ycoord1,ybush1,ycoord2,ybush2,By, Cy];
	ymina=min(tempy);
	ymaxa=max(tempy);
	rangex=xmaxa-xmina;
	rangey=ymaxa-ymina;
	xmin=xmina-0.05*rangex;
	xmax=xmaxa+0.05*rangex;
	ymin=ymina-0.05*rangey;
	ymax=ymaxa+0.05*rangey;
	values=axisadjust(xmin, xmax, ymin, ymax, width, height);
	xmin=values(1);
	xmax=values(2);
	ymin=values(3);
	ymax=values(4);
	axis equal
	axis([xmin xmax ymin ymax]);

	axes(h3);
	axis([0 360 min(Tin1) max(Tin1)]);  %set axis limits for the torque plot
	axes(h1);
	axis([0 360 min(Q44)-3 max(Q44)+3]);  %set axis limits for the rocker angle plot
	axes(h2);
	axis([0 360 min(tang) max(tang)]);  %set axis limits for the velocity plot

% Plot the results.

	itotal=i;
	ans='y';
	
% Draw the fixed curves first.

	set(Torque,'xdata', Q22,'ydata', Tin1);
	set(Displ,'xdata', Q22,'ydata', Q44);
	set(Tangle,'xdata', Q22,'ydata', tang);
	set(pin1,'xdata', xcoord1, 'ydata',ycoord1);
	set(bush1,'xdata', xbush1,'ydata', ybush1);
	set(pin2,'xdata', xcoord2, 'ydata',ycoord2);
	set(bush2,'xdata', xbush2,'ydata', ybush2);

% Animate the results as long as ans = 'y'
	
	timesink=1;
	while ans=='y'
		for i=1:1:itotal
			cycle=fix((i-0.0000001)/(360/dt)); %determine the cycle as an integer
			set(joint1,'xdata', Bx(i), 'ydata',By(i));
			set(joint2,'xdata', Cx(i), 'ydata',Cy(i));
			set(crank,'xdata',[0 Bx(i)],'ydata', [0 By(i)]);
			set(coupler,'xdata', [Bx(i) Cx(i)],'ydata', [By(i) Cy(i)]);
			set(rocker,'xdata', [(r1*cos(Q11)) Cx(i)],'ydata', ...
				[(r1*sin(Q11)) Cy(i)]);
			set(tbead, 'xdata', Q22(i)-cycle*360, 'ydata', Tin1(i));
			set(dbead, 'xdata', Q22(i)-cycle*360, 'ydata', Q44(i));
			set(vbead, 'xdata', Q22(i)-cycle*360, 'ydata', tang(i));
			drawnow;   %flush the draw buffer

% Do meaningless calculation to slow down the linkage

			for k=1:1:timesink
				zzz=log10(k);
			end
		end

% Plot figures in final position.  Plot position first.

		i = 6;
		set(joint1,'xdata', Bx(i), 'ydata',By(i));
		set(joint2,'xdata', Cx(i), 'ydata',Cy(i));
		set(crank,'xdata',[0 Bx(i)],'ydata', [0 By(i)]);
		set(coupler,'xdata', [Bx(i) Cx(i)],'ydata', [By(i) Cy(i)]);
		set(rocker,'xdata', [(r1*cos(Q11)) Cx(i)],'ydata', ...
	  		 [(r1*sin(Q11)) Cy(i)]);
   
% Plot Torque information

		set(tbead, 'xdata', Q22(i), 'ydata', Tin1(i));
 
% Plot angle information

		set(dbead, 'xdata', Q22(i), 'ydata', Q44(i));

% Plot transmission angle information

		set(vbead, 'xdata', Q22(i), 'ydata', tang(i));

%flush the draw buffer

		drawnow; 

% Ask if the animation should be repeated

		ans=input('Repeat animation? y/n [y]: ', 's');
		if isempty(ans);
			ans='y';
		end
		if ans=='y'
			rep=input('Change animation speed? (+/-/0) [-]:','s');
			if isempty(rep);rep='-';end
			if rep=='+'; timesink=timesink/10; end
			if rep=='-'; timesink=timesink*10; end
			if timesink<1; timesink=1; end
		end
	end

	ansdesign=input('Select another design? y/n [y]: ', 's');
	if isempty(ansdesign);
		ansdesign='y';
	end
end
