%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This Matlab script was created by:
%Kristopher Seluga, USA www.rclandsailing.com copyright 2002
%To run this program you must have Matlab by MathWorks Inc.
%Place this file in Matlab/bin directory and type "landsail" at the prompt to execute
%If you have comments or questions please write to kseluga@rclandsailing.com
%Also, if you use this program to simulate an actual model or full sized landyacht
%please send feedback on its usefulness and how the predictions compare to observations
%so that the algorithm may be validated and improved.  
%The author makes no claims to the accuracy the results of this program.
%Thank you and enjoy.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This script finds the maximum velocity vs heading of a landyacht based on:
%1. True wind speed
%2. Symmetric wing Lift/Drag correlations
%			a. Drag coefficient at zero angle of attack
%			b. Wing span, area & aspect ratio
%			c. Lift reduction after stall angle is reached
%3. Parasitic drag
%			a. Fuselage Drag
%			b. Crossbeam drag
%			c. Wheel rolling resistance
%4. Flipping criteria
%			a. Vehicle weight and geometry
%			b. Down force generation
%5. Side slipping
%			a. Available friction force
%			b. Down force generation
%6. Number of wings/sails
%7. Number of wheels (3 or 4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MAIN VARIABLE DEFINITIONS
%all units are inches, pounds, seconds and radians unless otherwise stated
%phi = angle between true wind velocity and vehicle velocity
%theta = angle of wing(s) with respect to vehicle velocity (wing/sail trim)
%beta = angle between apparent wind and vehicle velocity
%alpha = angle of attack for wing(s) (alpha = beta - theta)
%V_TRUE = velocity of true wind
%vy = vehicle velocity relative to ground
%AR_WING = wing aspect ratio
%F_LIFT_WING = lift force on wing (perpendicular to apparent wind)
%F_DRAG_WING = drag force on wing (parallel to apparent wind)
%F_THRUST = net force "thrust" on vehicle in direction of vehicle velocity
%F_SIDE = sideways force on vehicle that causes slipping
%FLIP = a test value for flipping criteria, 0=no, 1=yes
%SLIP = a test value for slipping criteria, 0=no, 1=yes
%V_LIMIT = limiting yacht velocity when slipping and flipping criteria are ignored
%V_MAX = actual predicted maximum yacht velocity based on flip/slip criteria
%"CG" refers to center of gravity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;	%clear memory before running program

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SINGLE VALUE PARAMETERS THAT REMAIN CONSTANT FOR ENTIRE SIMULATION
%Change these values to match your vehicle's geometry and configuration
%air properties and other constants
RHO_AIR = 1.175*10^-7;					%density of air (lb*sec^2/in^4)
MU_AIR = 2.7326*10^-8;					%viscosity of air (lb*sec/in^2)
g = 32.2;									%gravity (ft/sec^2)
%true wind speed
V_TRUE_MPH = 20;							%true wind velocity (mph)
V_TRUE = V_TRUE_MPH*(528/360)*12;	%true wind velocity (in/sec)
%wing properties (assuming symmetric section)
N_WING = 1;									%number of wings/sails
CD0_WING = 0.2;							%wing drag coefficient at zero angle of attack
S_WING = 48;								%wing span (in)
A_WING = 250;								%wing area (in^2)
HCE = 26;									%height of wing center of effort from ground (in)
psi_deg = 0;								%lean of mast from side stay slack (deg)
stall_deg = 18;							%wing stall angle- normally 16-20 deg (deg)
CL_DROP_FACTOR = 1.0;					%factor for exponetial lift dropoff after stall
psi = psi_deg*pi/180;					%lean of mast from side stay slack (rad)
stall = stall_deg*pi/180;				%stall angle (rad)
ARM_WING = HCE*cos(psi);				%moment arm for rig flipping moment (in)
C_WING = A_WING/S_WING;					%average wing chord length (in)
AR_WING = (S_WING^2)/A_WING;			%aspect ratio
CL_WING_MAX = 2*pi*sin(stall)/(1+2/AR_WING);		%maximum coefficient of lift before stall
%total yacht weight
WEIGHT = 4.0;								%total vehicle weight including pilot (lb)
%fuselage/body properties
AP_BODY = 10;								%total frontal fuselage projected area (in^2)
CD_BODY = 0.5;								%vehicle fuselage drag coefficient 
TRACK = 36;									%vehicle width/track (wheel center to center) (in)
WHEELBASE = 48;							%vehicle wheelbase (in)- only necessary for 3 wheel vehicle
DIST_CG = 9;								%horizontal distance from CG to rear wheels (in)- only necessary for 3 wheels
%wheel properties
N_WHEEL = 3;								%number of wheels
FRCT = 0.75;								%coefficient of friction between wheels and road
DIAM_WHEEL_FT = 3.0;						%front wheel diameter (in)
DIAM_WHEEL_RE = 4.0;						%rear wheel diameter (in)
B_WHEEL = 0.007;							%drag coefficient of rolling wheel, per wheel (lb/(rev/sec))
F_DRAG_WHEEL_MIN = 0.5;					%total rolling resistance at zero velocity (lbs)
%for small models on wheels with bearings B_WHEEL is on the order of 0.010 lb/(rev/sec)
%for full sized landyachts B_WHEEL is on the order of 0.10 lb/(rev/sec)
%crossbeam properties (lift-drag relations are for a symmetric naca section)
W_BEAM = 36;								%width of crossbeam available for generating down force (in)
T_BEAM = 0.5;								%maximum thickness of crossbeam wing section (in)
C_BEAM = 3.0;								%average crossbeam chord length (in)
A_BEAM = C_BEAM*W_BEAM;					%foil surface area per crossbeam (in^2)
AP_BEAM = T_BEAM*W_BEAM;				%projected frontal beam area (in^2)
AR_BEAM = (W_BEAM^2)/A_BEAM;			%aspect ratio of crossbeam
%relations for symmetric airfoil crossbeam with constant angle of attack
gamma_deg = 0;												%crossbeam angle of attack for generating down force (deg)
gamma = gamma_deg*pi/180;								%crossbeam angle of attack for generating down force (rad)
CD0_BEAM = 0.10;											%crossbeam drag coefficient at zero angle of attack
CL_BEAM = 2*pi*sin(gamma)/(1+2/AR_BEAM);			%crossbeam coefficient of lift
CD_BEAM = CD0_BEAM+((CL_BEAM^2)/(pi*AR_BEAM));	%crossbeam coefficient of drag
%END OF YACHT DATA ENTRY

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PRE-LOOP CALCULATIONS
%CG/BEAM moment arm for righting moment calculation
if N_WHEEL == 4										%for a 4 wheeled symmetric yacht (track front = track rear)
   N_WHEEL_FT = 2;									%two front wheels
   N_WHEEL_RE = 2;									%two rear wheels
   zeta = 0;											%angle between flipping force and F_SIDE is zero (rad)
   ARM_CG = TRACK/2;									%cg moment arm for righting moment (in)
   ARM_BEAM = TRACK/2;								%beam moment arm for righting moment (in)
   N_BEAM = 2;											%number of crossbeams for down force generation (probably has 2 crossbeams)
elseif N_WHEEL == 3									%for a 3 wheeled symmetric yacht (2 wheels in rear)
   N_WHEEL_FT = 1;									%one front wheel
   N_WHEEL_RE = 2;									%two rear wheels
   zeta = atan((TRACK/2)/WHEELBASE);			%angle made between front wheel and each rear wheel (rad)
   ARM_CG = (WHEELBASE-DIST_CG)*sin(zeta);	%cg moment arm for righting moment (in)
   ARM_BEAM = (TRACK/2)*cos(zeta);				%beam moment arm for righting moment (in)
	N_BEAM = 1;											%number of crossbeams for down force generation (rear crossbeam only)
end
%main program loop parameters/starting values
e = 0.050;						%error sensitivity for thrust equilibrium check (lbs)
dv = 0;							%initial value for velocity change (in/sec)
vy0 = 100;						%initial value for vehicle velocity (in/sec)
vy = vy0;						%assign initial value to vy (in/sec)
F_THRUST = 1;					%initial non-zero value for thrust force (lbs)
%define matrices PHI and THETA for main program loops
%for smoother plot make N_phi and N_theta larger (this will make the program take longer to run)
N_phi = 39;											%number of divisions in PHI
phi_max = pi;										%maximum phi = 180 degrees (pi rad)
PHI = [0:phi_max/N_phi:phi_max];				%vector phi in radians from head-to-wind to 180 degrees (rad)
P = N_phi+1;										%number of loops to be made for varying phi
N_theta = 120;										%number of divisions in THETA
theta_max = pi/3;									%maximum theta to test during loop (rad)
THETA = [0:theta_max/N_theta:theta_max];	%define vector of wing angles from zero to 90 deg (rad)
T = N_theta+1;										%number of loops to be made for varying theta

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%BEGIN MAIN PROGRAM LOOP THROUGH ALL VALUES OF PHI
for p = 1:1:P						%loop for entire range of phi (all possible headings)
	Steps_Remaining = P-p+1		%display number of steps remaining in loop to track progress
	phi = PHI(p);					%lookup value phi for this loop (rad)
   t = 1;							%initialize t counting variable

	while t <= T													%find VY_max for loop through for all values of theta
   %setup theta and vy for this value of t  
   theta=THETA(t);												%lookup current theta (rad)
   SIGN_F_THRUST = F_THRUST/abs(F_THRUST);				%find the sign of F_THRUST
   vy = vy + SIGN_F_THRUST*dv;								%increment/decrement guess for vehicle velocity (in/sec)
	   %do calculations necessary to find net thrust
		if V_TRUE*cos(phi) + vy > 0																	%if apparent wind is from bow
			beta = atan((V_TRUE*sin(phi))/(vy+V_TRUE*cos(phi)));								%calc beta for beta < pi/2 (rad)
   	else																									%if apparent wind is from stern
	   	beta = atan((V_TRUE*sin(phi-(pi/2))-vy)/(V_TRUE*cos(phi-(pi/2))))+pi/2;		%calc beta for beta > pi/2 (rad)
   	end																									%end if loop
   alpha = beta-theta;															%calc alpha - wing angle of attack (rad)
   V_APPARENT = sqrt((vy+V_TRUE*cos(phi))^2+(V_TRUE*sin(phi))^2);	%calc magnitude of apparent wind velocity (in/sec)
	V_FRONT = V_APPARENT*cos(beta);											%calc apparent wind velocity in direction of vehicle velocity (in/sec)
  	SIGN_V_FRONT = cos(beta)/abs(cos(beta));								%direction of front component of apparent wind (+1 from bow, -1 from stern)
   	%wing lift and drag with rough approximation to lift drop off after stall
   	if alpha > stall					 												%if angle of attack is less than stall angle
      	CL_WING = (CL_WING_MAX*exp(stall))*exp(-CL_DROP_FACTOR*alpha);	%calc reduced wing lift coefficient
   		CD_WING = CD0_WING+((CL_WING^2)/(pi*AR_WING));						%calc wing drag coefficient for this angle of attack
   	else																					%if angle of attack is greter than stall angle
      	CL_WING = 2*pi*sin(alpha)/(1+2/AR_WING);								%calc wing lift coefficient for this angle of attack
   		CD_WING = CD0_WING+((CL_WING^2)/(pi*AR_WING));						%calc wing drag coefficient for this angle of attack
   	end																					%end if loop
   AP_WING = A_WING*(1-((1-cos(psi))*alpha)/(pi/2));						%linear approximation to adjust effective wing area based on alpha and psi (in^2)
   %NOTE: when alpha is zero, AP_WING = A_WING, if alpha is 90 deg, AP_WING = A_WING*cos(psi)
   %this leads to the linear approximation: AP_WING = A_WING*(1-((1-cos(psi))*alpha)/(pi/2));
   F_LIFT_WING = N_WING*0.5*CL_WING*RHO_AIR*(V_APPARENT^2)*AP_WING;	%calc total wing lift force (lbs)
	F_DRAG_WING = N_WING*0.5*CD_WING*RHO_AIR*(V_APPARENT^2)*AP_WING;	%calc total wing drag force (lbs)
   %yacht/fuselage drag (not including crossbeams)
   F_DRAG_BODY = SIGN_V_FRONT*0.5*CD_BODY*RHO_AIR*(V_FRONT^2)*AP_BODY;	%calc fuselage drag force (lbs)
   %crossbeam lift (down force) and drag
   	if SIGN_V_FRONT == 1																	%for apparent wind from bow
		   F_LIFT_BEAM = N_BEAM*0.5*CL_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;		%calc total crossbeam lift force (lbs)
   		F_DRAG_BEAM = N_BEAM*0.5*CD_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;		%calc total crossbeam drag force (lbs)
   	else																						%for apparent wind from stern
      	F_LIFT_BEAM = 0;																	%crossbeam generates zero lift
   		F_DRAG_BEAM = -2*N_BEAM*0.5*CD_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;	%calc total crossbeam drag force (lbs)
   	end																						%end if loop
   %wheel drag (rolling viscous friction)
   RPS_WHEEL_FT = vy/(pi*DIAM_WHEEL_FT);														%calc front wheel rev per sec (rps)
   RPS_WHEEL_RE = vy/(pi*DIAM_WHEEL_RE);														%calc rear wheel rev per sec (rps)
   F_DRAG_WHEEL = B_WHEEL*(N_WHEEL_FT*RPS_WHEEL_FT + N_WHEEL_RE*RPS_WHEEL_RE);	%calc total wheel drag force (lbs)
   F_DRAG_WHEEL = max(F_DRAG_WHEEL,F_DRAG_WHEEL_MIN);										%use minimum wheel rolling resistance for low speeds (lbs)
   %sum up all forward thrust forces
	F_SIDE = F_LIFT_WING*cos(beta)*cos(psi)+F_DRAG_WING*sin(beta);		%calc sideways sliding force (lbs)
   F_THRUST = F_LIFT_WING*sin(beta)*cos(psi)-F_DRAG_WING*cos(beta)-F_DRAG_BODY-F_DRAG_BEAM-F_DRAG_WHEEL;	%calc net thrust force (lbs)
   %SLIP/FLIP required quantities (righting/flipping moment and side friction)
	F_DOWN = WEIGHT + F_LIFT_BEAM + F_LIFT_WING*sin(psi);					%total down force from crossbeam and wing (lb)
	M_RIGHT = WEIGHT*ARM_CG + F_LIFT_BEAM*ARM_BEAM;							%available righting moment (in*lb)
   F_FLIP = F_SIDE*cos(zeta) + F_THRUST*sin(zeta);							%component of rig force perpendicular to supporing wheels (lbs)
   M_FLIP = F_FLIP*ARM_WING;														%flipping moment from rig (in*lb)
	F_FRICT = F_DOWN*FRCT;															%total available sideways friction force (lbs)
   %END OF ALL FORCE AND MOMENT CALCULATIONS FOR THIS LOOP
   %scale velocity step change for next loop
   X = 10;													%scale factor for calculating dv
   dv = X*abs(F_THRUST);								%scale velocity change by F_THRUST (if X is too large loop becomes unstable)
		%perform slip criteria calculations
  		if F_SIDE > F_FRICT			%if side sliding force from rig exceeds available friction force
			SLIP(t) = 1;				%set SLIP to 1 (yes)
		else								%if friction is greater than sliding force
			SLIP(t) = 0;				%set SLIP to 0 (no)
		end								%end SLIP if loop
		%perform flip criteria calculations
		if M_FLIP > M_RIGHT			%if flipping moment from rig exceeds available righting moment
			FLIP(t) = 1;				%set FLIP to 1 (yes)
		else								%if righting moment is greater than flipping moment
			FLIP(t) = 0;				%set FLIP to 0 (no)
		end								%end FLIP if loop
		%check for positive thrust force near zero, if so go to next theta value
		if vy > 0 & alpha > 0		%only continue for loops that lead to positive velocity and angle of attack
			if abs(F_THRUST) < e		%if thrust is close to zero (lbs)
      	  	VY(t) = vy;				%store final value of vehicle speed from this loop (in/sec) 
        		t = t+1;					%increment t to calculate vehicle speed for next theta
     		end 	  						%end thrust check if loop
   	else								%if yacht velocity or angle of attack is not positive
     		VY(t) = 0;					%assign this velocity to zero (do not keep track of negative velocities)
      	t = t+1;						%increment t to calculate vehicle speed for next theta
      	vy = vy0;					%reset initial velocity guess
   	end								%end else loop
   end	 								%end of while  t <= T loop
	%truncate VY matrix to remove infinite element
   V_LIMIT(p) = max(VY);			%store V_LIMIT vector as maximum velocity before SLIP and FLIP are considered (in/sec)
	%change all velocity entries to zero if SLIP or FLIP = 1 (yes)
   for z = 1:1:T						%loop through entire matrix for each theta
 	  	if (SLIP(z) + FLIP(z) > 0)	%if FLIP or SLIP = 1 (yes)
  	  		VY(z) = 0;					%set this velocity element to zero
     	end								%end if loop
     end									%end for loop
   %store theta that corresponds to this VMAX
   V_MAX(p) = max(VY);				%store maximum velocity found for this angle phi to VMAX(p)
	q = 1;								%initialize counting var q = 1
   while VY(q) ~= V_MAX(p)			%while VY_TRUNK(q) not equal VMAX
      q = q+1;							%increment q until VY_TRUNK(q) = VMAX(p)
	end									%end while loop
   THETA_MAX(p) = THETA(q);		%store theta that corresponds to VMAX(p) (rad)
end   									%end of "for 1:1:P" loop
%END MAIN PROGRAM LOOP

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ORGANIZE DATA INTO MAXIMUM SPEED INFORMATION
TOP_SPEED_IN_SEC = max(V_MAX);							%store top speed in in/sec (in/sec)
TOP_SPEED_FT_SEC = TOP_SPEED_IN_SEC/12;				%convert top speed to ft/sec (ft/sec)
TOP_SPEED_MPH = TOP_SPEED_FT_SEC*(360/528);			%convert top speed to mph (mph)
V_LIMIT_MPH = (V_LIMIT/12)*(360/528);					%convert limit speed to mph (mph)
V_MAX_FT_SEC = V_MAX/12;									%convert maximum speed matrix to mph (ft/sec)
V_MAX_MPH = (V_MAX_FT_SEC)*(360/528);					%convert maximum speed matrix to mph (mph)
for j= 1:1:P-1													%loop through all of matrix V_MAX
   if V_MAX(j) == TOP_SPEED_IN_SEC						%look for entry of maximum velocity
      J = j;													%if this is that entry, save entry number
   end															%end if loop
end																%end for loop
THETA_TOP_SPEED_RAD = THETA_MAX(J);						%store theta that corresponds to TOP_SPEED_MPH (rad)
PHI_TOP_SPEED_RAD = PHI(J);								%store phi that corresponds to TOP_SPEED_MPH (rad)
THETA_TOP_SPEED_DEG = THETA_TOP_SPEED_RAD*180/pi;	%convert to degrees (deg)
PHI_TOP_SPEED_DEG = PHI_TOP_SPEED_RAD*180/pi;		%convert to degrees (deg)
%plot velocity results in polar coordinates
RHO_MAX = 3*ceil(max(V_LIMIT_MPH)/3);					%get maximum radius in whole numbers divisible by 3
MP = 1;															%subplot dimension
NP = 1;															%subplot dimension
subplot(MP,NP,1),polar(0,RHO_MAX);						%create polar plot with whole number scale
subplot(MP,NP,1),hold;										%hold plot so data can be plotted
subplot(MP,NP,1),polar(PHI,V_LIMIT_MPH,'ro--');		%polar plot limit velocity vs phi in red
subplot(MP,NP,1),polar(-PHI,V_LIMIT_MPH,'ro--');	%polar plot limit velocity vs -phi (mirror image)
subplot(MP,NP,1),polar(PHI,V_MAX_MPH,'bx-'); 		%polar plot maximum velocity vs phi in blue
subplot(MP,NP,1),polar(-PHI,V_MAX_MPH,'bx-');		%polar plot maximum velocity vs -phi (mirror image)
%subplot(MP,NP,1),title('Polar Vel. Plot (mph)');		%title plot
subplot(MP,NP,1),hold;										%release plot to return to default condition
%display desired values
TOP_SPEED_MPH = round(TOP_SPEED_MPH)					%round and print maximum velocity (mph)
PHI_TOP_SPEED_DEG = round(PHI_TOP_SPEED_DEG)			%round and print phi for TOP SPEED (deg)
THETA_TOP_SPEED_DEG = round(THETA_TOP_SPEED_DEG)	%round and print theta for TOP SPEED (deg)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%END MAIN PROGRAM
%BEGIN POST LOOP CALCULATIONS FOR STORING & DISPLAYING OTHER RELEVANT INFO

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CALCULATE APPROXIMATE TIME REQUIRED TO REACH 95% of TOP SPEED IGNORING SLIP/FLIP CRITERIA
phi = PHI_TOP_SPEED_RAD;												%heading for top speed
theta = THETA_TOP_SPEED_RAD;											%wing angle at top speed
vy = 88;																		%begin at 5 mph yacht velocity
dt = 0.2;																	%time increment (sec)
TIME_MAX = 60*10;															%maximum time before loop is cut off (sec)
time = 0;																	%start at time = zero
while vy < 0.95*TOP_SPEED_IN_SEC & time < TIME_MAX											%keep looping until vy is 95% of top speed
	if V_TRUE*cos(phi) + vy > 0																	%if apparent wind is from bow
		beta = atan((V_TRUE*sin(phi))/(vy+V_TRUE*cos(phi)));								%calc beta for beta < pi/2 (rad)
   else																									%if apparent wind is from stern
	   beta = atan((V_TRUE*sin(phi-(pi/2))-vy)/(V_TRUE*cos(phi-(pi/2))))+pi/2;		%calc beta for beta > pi/2 (rad)
   end																									%end if loop
   alpha = beta-theta;															%calc alpha - wing angle of attack (rad)
   V_APPARENT = sqrt((vy+V_TRUE*cos(phi))^2+(V_TRUE*sin(phi))^2);	%calc magnitude of apparent wind velocity (in/sec)
	V_FRONT = V_APPARENT*cos(beta);											%calc apparent wind velocity in direction of vehicle velocity (in/sec)
  	SIGN_V_FRONT = cos(beta)/abs(cos(beta));								%direction of front component of apparent wind (+1 from bow, -1 from stern)
   %wing lift and drag with rough approximation to lift drop off after stall
   if alpha > stall					 												%if angle of attack is less than stall angle
    	CL_WING = (CL_WING_MAX*exp(stall))*exp(-CL_DROP_FACTOR*alpha);	%calc reduced wing lift coefficient
   	CD_WING = CD0_WING+((CL_WING^2)/(pi*AR_WING));						%calc wing drag coefficient for this angle of attack
   else																					%if angle of attack is greter than stall angle
     	CL_WING = 2*pi*sin(alpha)/(1+2/AR_WING);								%calc wing lift coefficient for this angle of attack
   	CD_WING = CD0_WING+((CL_WING^2)/(pi*AR_WING));						%calc wing drag coefficient for this angle of attack
   end																					%end if loop
	AP_WING = A_WING*(1-((1-cos(psi))*alpha)/(pi/2));						%linear approximation to adjust projected wing area based on alpha and psi (in^2)
   %NOTE: when alpha is zero, AP_WING = A_WING, if alpha is 90 deg, AP_WING = A_WING*cos(psi)
   %this leads to the linear approximation: AP_WING = A_WING*(1-((1-cos(psi))*alpha)/(pi/2));
   F_LIFT_WING = N_WING*0.5*CL_WING*RHO_AIR*(V_APPARENT^2)*AP_WING;	%calc total wing lift force (lbs)
	F_DRAG_WING = N_WING*0.5*CD_WING*RHO_AIR*(V_APPARENT^2)*AP_WING;	%calc total wing drag force (lbs)
   %yacht/fuselage drag (not including crossbeams)
   F_DRAG_BODY = SIGN_V_FRONT*0.5*CD_BODY*RHO_AIR*(V_FRONT^2)*AP_BODY;	%calc fuselage drag force (lbs)
   %crossbeam lift (down force) and drag
   if SIGN_V_FRONT == 1																	%for apparent wind from bow
	   F_LIFT_BEAM = N_BEAM*0.5*CL_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;		%calc total crossbeam lift force (lbs)
   	F_DRAG_BEAM = N_BEAM*0.5*CD_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;		%calc total crossbeam drag force (lbs)
   else
      F_LIFT_BEAM = 0;
   	F_DRAG_BEAM = -2*N_BEAM*0.5*CD_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;	%calc total crossbeam drag force (lbs)
   end
   %wheel drag (rolling viscous friction)
   RPS_WHEEL_FT = vy/(pi*DIAM_WHEEL_FT);														%calc front wheel rev per sec (rps)
   RPS_WHEEL_RE = vy/(pi*DIAM_WHEEL_RE);														%calc rear wheel rev per sec (rps)
   F_DRAG_WHEEL = B_WHEEL*(N_WHEEL_FT*RPS_WHEEL_FT + N_WHEEL_RE*RPS_WHEEL_RE);	%calc wheel drag force (lbs)
   F_DRAG_WHEEL = max(F_DRAG_WHEEL,F_DRAG_WHEEL_MIN);										%use minimum wheel rolling resistance for low speeds (lbs)
	%sum up all forward thrust forces
   F_THRUST = F_LIFT_WING*sin(beta)*cos(psi)-F_DRAG_WING*cos(beta)-F_DRAG_BODY-F_DRAG_BEAM-F_DRAG_WHEEL;	%calc net thrust force (lbs)
	ACCEL_G = F_THRUST/WEIGHT;														%calculate acceleration (g's)
	ACCEL_IN_SEC = ACCEL_G*g*12;													%acceleration in (in/sec^2)
   vy = vy + dt*ACCEL_IN_SEC;														%increase velocity based on acceleration and time step (in/sec)
   time = time + dt;																	%increment time by one time step (sec)
end
if time < TIME_MAX
   TIME_TO_TOP_SPEED_SEC = time;													%time to top speed (sec)
	TIME_TO_TOP_SPEED_MIN = TIME_TO_TOP_SPEED_SEC/60;						%time to top speed (min)
	TIME_TO_TOP_SPEED_MIN = round(10*TIME_TO_TOP_SPEED_MIN)/10			%display time (min)
else
   disp('Accleration calc failed')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CALCULTE & STORE RELEVANT PARAMETERS AT FULL SPEED FOR ALL HEADINGS
for p = 1:1:P													%loop for entire range of phi (all possible headings)
	vy = V_MAX(p);												%top speed (in/sec)
   phi = PHI(p);												%heading (rad)
   %do calculations necessary to find net thrust
   if vy > 0																								%only proceed for non zero velocities
		if V_TRUE*cos(phi) + vy > 0																	%if apparent wind is from bow
			BETA(p) = atan((V_TRUE*sin(phi))/(vy+V_TRUE*cos(phi)));							%calc beta for beta < pi/2 (rad)
   	else																									%if apparent wind is from stern
		   BETA(p) = atan((V_TRUE*sin(phi-(pi/2))-vy)/(V_TRUE*cos(phi-(pi/2))))+pi/2;	%calc beta for beta > pi/2 (rad)
	   end																									%end if loop
      ALPHA(p) = BETA(p)-THETA_MAX(p);												%calc alpha (rad)
		V_APPARENT(p) = sqrt((vy+V_TRUE*cos(phi))^2+(V_TRUE*sin(phi))^2);	%calc magnitude of apparent wind velocity (in/sec)
		V_FRONT = V_APPARENT(p)*cos(BETA(p));										%calc apparent wind velocity along vehicle velocity direction (in/sec)
	  	SIGN_V_FRONT = cos(BETA(p))/abs(cos(BETA(p)));							%direction of front component of apparent wind (+1 from bow, -1 from stern)
   	%wing lift and drag with rough approximation to lift drop off after stall
   	if ALPHA(p) > stall				 													%if angle of attack is less than stall angle
      	CL_WING = (CL_WING_MAX*exp(stall))*exp(-CL_DROP_FACTOR*ALPHA(p));	%calc reduced wing lift coefficient
   		CD_WING = CD0_WING+((CL_WING^2)/(pi*AR_WING));							%calc wing drag coefficient for this angle of attack
   	else																						%if angle of attack is greter than stall angle
      	CL_WING = 2*pi*sin(ALPHA(p))/(1+2/AR_WING);								%calc wing lift coefficient for this angle of attack
   		CD_WING = CD0_WING+((CL_WING^2)/(pi*AR_WING));							%calc wing drag coefficient for this angle of attack
   	end																						%end if loop
		AP_WING = A_WING*(1-((1-cos(psi))*ALPHA(p))/(pi/2));							%linear approximation to adjust projected wing area based on alpha and psi (in^2)
		F_LIFT_WING(p) = N_WING*0.5*CL_WING*RHO_AIR*(V_APPARENT(p)^2)*AP_WING;	%calc total wing lift force (lbs)
		F_DRAG_WING(p) = N_WING*0.5*CD_WING*RHO_AIR*(V_APPARENT(p)^2)*AP_WING;	%calc total wing drag force (lbs)
   	%yacht/fuselage drag (not including crossbeams)
		F_DRAG_BODY(p) = SIGN_V_FRONT*0.5*CD_BODY*RHO_AIR*(V_FRONT^2)*AP_BODY;	%calc fuselage drag force (lbs)
		%crossbeam lift (down force) and drag
      if SIGN_V_FRONT == 1																		%for apparent wind from bow
		   F_LIFT_BEAM(p) = N_BEAM*0.5*CL_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;		%calc total crossbeam lift force (lbs)
   		F_DRAG_BEAM(p) = N_BEAM*0.5*CD_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;		%calc total crossbeam drag force (lbs)
   	else
      	F_LIFT_BEAM(p) = 0;
   		F_DRAG_BEAM(p) = -2*N_BEAM*0.5*CD_BEAM*RHO_AIR*(V_FRONT^2)*A_BEAM;	%calc total crossbeam drag force (lbs)
   	end
   	%wheel drag (rolling viscous friction)
   	RPS_WHEEL_FT = vy/(pi*DIAM_WHEEL_FT);															%calc front wheel rev per sec (rps)
   	RPS_WHEEL_RE = vy/(pi*DIAM_WHEEL_RE);															%calc rear wheel rev per sec (rps)
   	F_DRAG_WHEEL(p) = B_WHEEL*(N_WHEEL_FT*RPS_WHEEL_FT + N_WHEEL_RE*RPS_WHEEL_RE);	%calc wheel drag force (lbs)
	   F_DRAG_WHEEL(p) = max(F_DRAG_WHEEL(p),F_DRAG_WHEEL_MIN);									%use minimum wheel rolling resistance for low speeds (lbs)
		%sum up all forward (thrust) and sideways (slide) forces
		F_THRUST(p) = F_LIFT_WING(p)*sin(BETA(p))*cos(psi)-F_DRAG_WING(p)*cos(BETA(p))-F_DRAG_BODY(p)-F_DRAG_BEAM(p)-F_DRAG_WHEEL(p);		%calc net thrust force (lbs)
		F_SIDE(p) = F_LIFT_WING(p)*cos(BETA(p))*cos(psi)+F_DRAG_WING(p)*sin(BETA(p));						%calc sideways sliding force (lbs)
   else																					%if vy is zero, then set all stored values to zero
      BETA(p) = 0;
      ALPHA(p) = 0;
		V_APPARENT(p) = 0;
		F_LIFT_WING(p) = 0;
		F_DRAG_WING(p) = 0;
		F_DRAG_BODY(p) = 0;
		F_LIFT_BEAM(p) = 0;
		F_DRAG_BEAM(p) = 0;
		F_DRAG_WHEEL(p) = F_DRAG_WHEEL_MIN;
		F_THRUST(p) = 0;
		F_SIDE(p) = 0;
   end																				%end if loop
end																					%end for loop
%print out desired values
REYNOLDS = RHO_AIR*C_WING*V_APPARENT/MU_AIR;								%Reynolds number
APPARENT_WIND_SPEED_MPH = round(V_APPARENT*(360/528)/12);			%apparent wind speed (mph)
DOWN_FORCE_LB = round(10*(F_LIFT_BEAM))/10;								%down force from crossbeam (lbs)
SIDE_FORCE_LB = round(10*F_SIDE)/10;										%sliding force (lbs)
BODY_AERO_DRAG_LB = round(10*(F_DRAG_BODY + F_DRAG_BEAM))/10;		%total body aerodynamic drag (mph)
WHEEL_DRAG_LB = round(10*F_DRAG_WHEEL)/10;								%wheel drag (lbs)
RIG_LIFT_LB = round(10*F_LIFT_WING)/10;									%total rig lift (lbs)
RIG_DRAG_LB = round(10*F_DRAG_WING)/10;									%total rig drag (lbs)
PHI_DEG = PHI*180/pi;															%convert PHI to degrees
ALPHA_DEG = ALPHA*180/pi;														%convert ALPHA to degrees
BETA_DEG = BETA*180/pi;															%convert BETA to degrees
THETA_MAX_DEG = THETA_MAX*180/pi;											%convert THETA_MAX to degrees
%plot angular info (change NP to 2 to use this plot)
%subplot(MP,NP,2),plot(PHI_DEG,BETA_DEG,'k');							%plot apparent wind angle v. phi
%subplot(MP,NP,2),hold;															%hold plot
%subplot(MP,NP,2),plot(PHI_DEG,THETA_MAX_DEG,'b');						%plot wing trim angle angle v. phi
%subplot(MP,NP,2),plot(PHI_DEG,ALPHA_DEG,'r');							%plot angle of attack wind angle v. phi
%subplot(MP,NP,2),title('BETA, THETA MAX & ALPHA');					%title plot
%subplot(MP,NP,2),hold;															%hold plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%END OF PROGRAM