% Investigation of lateral-directional trim in a crosswind
% See page 223 of Roskam for aircraft data
%
clear
disp(' ')
disp(' Start here')
disp(' ')
%
% INPUTS
%
W=8750 % pounds
rho=.002378 % slugs/ft^3
U1=270 % ft/sec
S=253  % ft^2
LT1=0  % ft-lbf
NT1=0  % ft-lbf
DND1=0 % ft-lbf
b=38   % ft
Cyb=-.0105  % per degree
Cyda=0      % per degree
Cydr= .0021 % per degree
Clb=-.0029  % per degree
Clda=.0024  % per degree
Cldr=.0002  % per degree
Cnb=.0018   % per degree
Cnda=.0005  % per degree
Cndr=-.001  % per degree
gamma=0     % degree
deltabeta=1 % degree
nbeta=12
%
% Problem statement
% Given beta, find phi(deg), da(deg) and dr(deg)
% x=sin(phi),da,dr
% Ax=B*beta
%
% Set up dimensions
%
BETA=zeros(1,nbeta);
DA=zeros(1,nbeta);
DR=zeros(1,nbeta);
PHI=zeros(1,nbeta);
xwindkn=zeros(1,nbeta);
%
% Compute intermediate items
%
d2r=pi/180
qbar=.5*rho*U1*U1
U1knots=U1*.5925 % knots
a1=[W*cos(gamma*d2r)/(qbar*S), 0, 0]'
a2=[Cyda, Clda, Cnda]'
a3=[Cydr, Cldr, Cndr]'
A=[a1 a2 a3]
B=[-Cyb; -Clb; -Cnb]
AIB=inv(A)*B;
%
for i=0:nbeta-1
       beta=i*deltabeta;
       x=AIB*beta;
       PHI(i+1)=asin(x(1))/d2r;
       DA(i+1)=x(2);
       DR(i+1)=x(3);
       BETA(i+1)=beta;
end
windkn=sin(BETA*d2r)*U1*.5925
%
% Plot results vs beta
%
figure(1)
subplot(3,1,1)
plot(BETA,DR)
heading=['Landing in a crosswind, U1= ',num2str(U1knots),' knots, right sideslip (beta+)']
title(heading)
ylabel2('rudder (deg)')
text2(.05,.8,'Dr+ means trailing edge left, push left pedal, nose left of wind, wind hits on the right side')
%
subplot(3,1,2)
plot(BETA,DA)
ylabel2('aileron (deg)')
text2(.05,.8,'Da+ means right aileron up, push stick right')
%
subplot(3,1,3)
plot(BETA,PHI)
text2(.05,.6,'Right wing down, right bank angle')
ylabel2('phi (deg)',-.10,.9)
xlabel('beta (deg)')
%
% Plot results vs crosswind velocity in knots
%
figure(2)
subplot(3,1,1)
plot(windkn,DR)
title(heading)
xlabel('crosswind (knots)')
ylabel2('rudder (deg)')
text2(.05,.8,'Dr+ means trailing edge left, push left pedal, nose left of wind, wind hits on the right side')
text2(.05,.65,'Crosswind=U1*sin(beta)')
%
subplot(3,1,2)
plot(windkn,DA)
xlabel('crosswind (knots)')
ylabel2('aileron (deg)')
text2(.05,.8,'Da+ means right aileron up, push stick right')
%
subplot(3,1,3)
plot(windkn,PHI)
xlabel('crosswind (knots)')
text2(.05,.6,'Right wing down, right bank angle')
ylabel2('phi (deg)',-.10,.9)




