% Experiment #1 Stall Speed and CLmax Determination
% Ref: NTPS, pages D24-D26.
% This MATLAB script computes stall airspeed and CLmax
% given Vindicated, flap setting, indicated geopotential pressure altitude.
% It assumes standard outside air temperature for the calibrated pressure altitude.
disp(' ')
disp('Start here')
disp('Experiment 1') %Example 1
disp(' ')
clear all
echo on
% List all inputs: Data Inputs
Hi=[4500,4100,4400]       % Item 1 of the Table of Outputs: indicated geopotential feet at stall. Put your data in this vector.
Vi=  [61,60,59]           % Item 2: stall speed knots. Put your data in this vector and the one below.
flap=[ 0, 0, 0]           % Flap setting %flap=0 for no flap %flap=0.5 for half flap %flap=1 for full flap
FuelWeight=[900,800,700]  % Pounds (lbf).
Viwarn=[64,66,63]         % Item 10, indicated stall warning speed knots.
ZeroFuelWeight=4500       % Pounds <<< Check this number. Use the number Céline gives you.
Ws=5500                   % Pounds (lbf), standard weight is maximum takeoff weight (Baron 58 manual p. 1-12)
S=199.2;                  % Wing area, ft squared
% When you do your experiment you will replace this data with your data.
% All three vectors above must be of the same length.

% Compute stall speed and Clmax and build a Table of Outputs 
% as called for in the Data Reduction section of the experiment.

% Item 1 in the Table of Outputs is indicated pressure altitude 
% for the flight experiment and has been given above.

% Item 2 in the Table of Outputs is indicated stall airspeed
% for the flight experiment and has been given above.

% The following step is not listed in the NTPS Data Reduction procedure, 
% but we will do this in this and other experiments.
% We will use the function ALTBaron58 which performs two useful 
% data corrections at once. It computes calibrated pressure altitude
% from indicated pressure altitude, indicated airspeed and flap setting.
% In the process it determines the instrument corrections and 
% static position error corrections. Unfortunately it does not
% return these two errors individually as is required for the Table of Outputs,
% so we will have to do a few awkward steps to get data in the required format
% for the Table of Outputs

help ALTBaron58

% Correct indicated altitude for instrument and position errors.
Hcalibrated=ALTBaron58(Hi,Vi,flap) % Altitude units are feet.
% Notice that we are using array operations to compute Hcalibrated, i.e., 
% Hi and Vi are arrays and returned variable Hcalibrated is an array.
% For most computations done below we will use array operations.


% Items 3 and 4  are the individual corrections applied to 
% airspeed measurements, 
%  dVic is the speed increment to add to indicated airspeed to
% correct for instrument error.
%         Vic=Vi+dVic
% dVpc is the increment to add to airspeed to correct for 
% static position error.
%         Vcalibrated=Vic+dVpc
%                    =Vi+dVic+dVpc

% Item 3 in the Table of Outputs: airspeed instrument correction (dVic)
% The Baron manual assumes that Vinstrumentcorrected (Vic) = Vindicated (Vi)
% so dVic = 0
dVic=zeros(size(Vi)); % for Item 3
Vic=Vi+dVic; % sort of trivial since we assumed dVic=0
% Airspeed units are KNOTS.

% Item 4 in the Table of Outputs: airspeed position correction (dVpc)
% To do this we will use a function called AS1Baron58. This
% function corrects indicated airspeed for instrument error 
% and static position error. Unfortunately it does not
% return these two errors individually, as is required for the Table of Outputs,
% so we will have to do a few awkward steps to get data in the required format.

help AS1Baron58

Vcalibrated=AS1Baron58(Vi,flap) % Airspeed units are KNOTS.
dVpc=Vcalibrated-Vic; % for Item 4. This is the awkward step.
% Airspeed units are KNOTS.

% Items 5 and 6 Compute dVc and equivalent airspeed
% dVc is the Scale Altitude Correction a.k.a. "Compressibility Correction".

% Item 6a Get standard pressure (p, lbf/ft^2) for calibrated pressure altitude.
GeometricFlag=0 %Hvector is Geopotential altitude
[temp,p,rho,Hgeopvector]=atmosphere4(Hcalibrated,GeometricFlag);
p
% Item 6b Get equivalent airspeed (Vequivalent)
% To do this we will use the function AS2 to compute
% equivalent velocity from calibrated velocity and static pressure at the test
% condition. Unfortunately, this function does not return dVc, so we
% will have to perform an awkward step to compute it.

help AS2

[Vequivalent,qcOverp0,qcOverp]=AS2(Vcalibrated,p); % for line 6, Knots
dVc=Vequivalent-Vcalibrated; % for Item 5. This is the awkward step.
Vequivalent                  % Airspeed units are KNOTS.

% Item 7 Compute test weights.
Wt=ZeroFuelWeight+FuelWeight % for Item 7 of the Table of Outputs


% Compute Item 8: CLmax
% To do this we will use two different methods.

% Item 8 done using equivalent airspeed (Vequivalent, knots)
rho0=0.00237691267925741;  % Sea level density, slugs per cubic foot
Vefps=1.687808*Vequivalent % convert knots to ft/sec
Clmax=2*Wt./(rho0*Vefps.*Vefps*S) % for Item 9 of the Table of Outputs

% Item 8 of the Table of Outputs done using true airspeed (Vtrue, knots)
% In this calculation we assume standard temperature for the calibrated pressure altitude.
% It would be better if we had measured outside air temperature, but we did not in this experiment.
temp     %degree Rankine
% To do this we use the function AS3. This function computes true airspeed 
% from equivalent airspeed, outside air temperature and calibrated altitude.

help AS3

[Vtrue,M,aknots,rho,sigma]=AS3(Vequivalent,temp,Hcalibrated); % knots
Vtfps=1.687808*Vtrue % convert knots to ft/sec
Clmax2=2*Wt./(rho.*Vtfps.*Vtfps*S) % Note that this is dene with array operations.

% Item 9 of the Table of Outputs: Stall speed corrected to the standard weight
Vs1=Vequivalent.*sqrt(Ws./Wt) % for Item 9 of the Table of Outputs, knots


% *********
% Compute warning speed. To do this we reepeat all the airspeed calibration steps.

% Item 10 Indicated airspeed at stall warning.
% This was input above in the Data Inputs section of this script.
% The following step is not listed in the NTPS Data Reduction procedure, 
% but we will do this in this and other experiments.
% Correct indicated altitude for instrument and position errors.
Hcalibratedwarn=ALTBaron58(Hi,Viwarn,flap)
% The Baron manual assumes that Vinstrumentcorrected (Vic) = Vindicated (Vi)
% so dVic = 0
dVicwarn=zeros(size(Viwarn)); % forItem 11
Vicwarn=Viwarn+dVicwarn;

% Items 11 and 12 Correct for indicated airspeed for instrument, position errors.
Vcalibratedwarn=AS1Baron58(Viwarn,flap) % Airspeed units are KNOTS.
dVpcwarn=Vcalibratedwarn-Vicwarn; % for Item 12

% Items 13 and 14 Compute equivalent airspeed
% Item 14a Get standard pressure (p, lbf/ft^2) for calibrated pressure altitude.
GeometricFlag=0 %Hcalibratedwarn is Geopotential altitude
[temp,pwarn,rho,Hgeopvector]=atmosphere4(Hcalibratedwarn,GeometricFlag);
pwarn
% Item 14b Get equivalent airspeed (Vequivalent)
[Vequivalentwarn,qcOverp0,qcOverp]=AS2(Vcalibratedwarn,pwarn); %Knots
Vequivalentwarn % for Item 14
dVcwarn=Vequivalentwarn-Vcalibratedwarn; % for Item 13

% Item 15 Warning speed corrected for weight
Vwarn=Vequivalentwarn.*sqrt(Ws./Wt) %for Item 15,  knots

% Item 16 Stall warning
StallWarning=Vwarn-Vs1 % for Item 16, knots

% Now let's call a script that writes a formatted table
%to the command window in a pretty format. We use a script here simply to
% cut down on the size of this main script (which is already too long).
echo off
MakeOutputTable
