%%
%% Wing profile modifications for Short S.23.
%% Anders Gidenstam  2010

%% Airfoil for the main wing:
%% Göttingen 436 airfoil. The S.23 had a modified Gö 436 airfoil.
%% "The modified Gö 436 aerofoil used had a thickness to chord ratio of
%% 18.78 % at the junction with the hull, diminishing to an
%% hypothetical 9% thickness to chord ratio at the tip." [Cassidy:2004:FE].

%% Göttingen 436:
Go436upper = [
   0.00000   0.00000
   0.01250   0.02231
   0.02500   0.03263
   0.05000   0.04625
   0.07500   0.05787
   0.10000   0.06650
   0.15000   0.07925
   0.20000   0.08750
   0.30000   0.09250
   0.40000   0.08950
   0.50000   0.08300
   0.60000   0.07200
   0.70000   0.05850
   0.80000   0.04100
   0.90000   0.02200
   0.95000   0.01125
   1.00000   0.00000
];
Go436lower = [
   0.00000   0.00000
   0.01250  -0.01469
   0.02500  -0.01987
   0.05000  -0.02275
   0.07500  -0.02283
   0.10000  -0.02250
   0.15000  -0.02125
   0.20000  -0.02000
   0.30000  -0.01750
   0.40000  -0.01500
   0.50000  -0.01250
   0.60000  -0.01000
   0.70000  -0.00750
   0.80000  -0.00500
   0.90000  -0.00250
   0.95000  -0.00125
   1.00000   0.00000
];

Go436xcoord = [0:0.005:0.1 0.15:0.05:1];
Go436yupper = interp1(Go436upper(:,1), 1.3*Go436upper(:,2), Go436xcoord, 'spline');
Go436ylower = interp1(Go436lower(:,1), 1.3*Go436lower(:,2), Go436xcoord, 'spline');

%% Plot the profile.
figure();
plot(%upper(:,1), upper(:,2), upper(:,1), lower(:,2),
     Go436xcoord(:,1), 1.3*Go436upper(:,2),
     Go436xcoord(:,1), 1.3*Go436lower(:,2),
     Go436xcoord, Go436yupper,
     Go436xcoord, Go436ylower
    );
axis("equal");

%% PHETE and PHETEP based on the airfoil for DATCOM
%%   PHETE  = 1/2 * (Y_90 - Y_99) / 9
%%   PHETEP = 1/2 * (Y_95 - Y_99) / 4

Y_90wng = interp1(Go436upper(:,1), 1.3*Go436upper(:,2), 0.90, 'spline') - \
          interp1(Go436lower(:,1), 1.3*Go436lower(:,2), 0.90, 'spline');
Y_95wng = interp1(Go436upper(:,1), 1.3*Go436upper(:,2), 0.95, 'spline') - \
          interp1(Go436lower(:,1), 1.3*Go436lower(:,2), 0.95, 'spline');
Y_99wng = interp1(Go436upper(:,1), 1.3*Go436upper(:,2), 0.99, 'spline') - \
          interp1(Go436lower(:,1), 1.3*Go436lower(:,2), 0.99, 'spline');

PHETEwng = 0.5 * (Y_90wng - Y_99wng) / 9
PHETEPwng = 0.5 * (Y_95wng - Y_99wng) / 4


%% Airfoil for the horizontal tail:
%% "The tail plane aerofoils used at Rochester were normally of RAF
%% 30 symmetrical section. For the Empire 'boats the thickness to
%% chord ratio was increased from the standard 12.64% to 13.75%.",
%% [Cassidy:2004:FE].

raf30upper = [
 0.0000000 0.0000000
 0.0125000 0.0180000
 0.0250000 0.0248000
 0.0500000 0.0346000
 0.1000000 0.0468000
 0.1500000 0.0544000
 0.2000000 0.0594000
 0.3000000 0.0632000
 0.4000000 0.0620000
 0.5000000 0.0566000
 0.6000000 0.0478000
 0.7000000 0.0370000
 0.8000000 0.0250000
 0.9000000 0.0130000
 0.9500000 0.0070000
 1.0000000 0.0000000
];
raf30lower = [
 0.0000000 0.0000000
 0.0125000 -.0180000
 0.0250000 -.0248000
 0.0500000 -.0346000
 0.1000000 -.0468000
 0.1500000 -.0544000
 0.2000000 -.0594000
 0.3000000 -.0632000
 0.4000000 -.0620000
 0.5000000 -.0566000
 0.6000000 -.0478000
 0.7000000 -.0370000
 0.8000000 -.0250000
 0.9000000 -.0130000
 0.9500000 -.0070000
 1.0000000 0.0000000
];

raf30xcoord = [0:0.005:0.1 0.15:0.05:1];
raf30yupper = interp1(raf30upper(:,1), 1.09*raf30upper(:,2), raf30xcoord, 'spline');
raf30ylower = interp1(raf30lower(:,1), 1.09*raf30lower(:,2), raf30xcoord, 'spline');
% Plain raf30 for the vertical tail (guess).
%raf30yupper = interp1(raf30upper(:,1), raf30upper(:,2), raf30xcoord, 'spline');
%raf30ylower = interp1(raf30lower(:,1), raf30lower(:,2), raf30xcoord, 'spline');

%% Plot the profile.
figure();
plot(%upper(:,1), upper(:,2), upper(:,1), lower(:,2),
     raf30upper(:,1), raf30upper(:,2), raf30upper(:,1), raf30lower(:,2),
     raf30xcoord, raf30yupper,
     raf30xcoord, raf30ylower
    );
axis("equal");

%% PHETE and PHETEP based on the airfoil for DATCOM
%%   PHETE  = 1/2 * (Y_90 - Y_99) / 9
%%   PHETEP = 1/2 * (Y_95 - Y_99) / 4

Y_90ht = interp1(raf30upper(:,1), 1.09*raf30upper(:,2), 0.90, 'spline') - \
         interp1(raf30lower(:,1), 1.09*raf30lower(:,2), 0.90, 'spline');
Y_95ht = interp1(raf30upper(:,1), 1.09*raf30upper(:,2), 0.95, 'spline') - \
          interp1(raf30lower(:,1), 1.09*raf30lower(:,2), 0.95, 'spline');
Y_99ht = interp1(raf30upper(:,1), 1.09*raf30upper(:,2), 0.99, 'spline') - \
         interp1(raf30lower(:,1), 1.09*raf30lower(:,2), 0.99, 'spline');

PHETEht = 0.5 * (Y_90ht - Y_99ht) / 9
PHETEPht = 0.5 * (Y_95ht - Y_99ht) / 4
