Page Comparison

%SI (m, kg, m/s, N)

%General
MASS_TOTAL = 5.137;
rho = 1.225;
V = 25;
Ixx = 0.06; %moment of inertia about the roll axis (kg*m^2)
Sv = 0.0423;
Sh = 0.0666;

%Wing geometry
b = 1.5;
fuselage_width = 0.132;
C = 0.23;
Sw = 0.345;
taper = 1;
MAC = 0.23;
NACA = 4412;
Cr = C; %

%Control surface effectiveness param data points
CA_C_ratio = [0.004918032786885199,0.03442622950819667,0.06967213114754095, 0.09918032786885245, 0.12213114754098359, 0.19344262295081968, 0.28196721311475414, 0.4024590163934427, 0.5016393442622952, 0.5434426229508198, 0.6778688524590166
];
TAU = [-0.002272727272726982, 0.10000000000000009, 0.20000000000000007, 0.26363636363636367, 0.3045454545454547, 0.4, 0.5, 0.6068181818181819, 0.6886363636363637, 0.7181818181818183, 0.8136363636363637
];

%Aileron inputs 
span_ratio = 0.4;
chord_ratio = 0.25;
bi = 1.5*0.5-span_ratio*1.5*0.5;
bo = 1.5*0.5;
Ca = chord_ratio*C;
tau = interpolate_tau(CA_C_ratio, TAU, chord_ratio); %Aileron efficiency paramter 
delta_max = [25, 25]; % [UP, DOWN]
CLw_slope = 6.03;%cl per rad

C_DR = 0.9; %avg drag from wing, horizontal tail, vertical tail 
yD = (b-fuselage_width)*0.25;
phi_req = 90;
t_req = 1.7;



%EQUATIONS
eq_Cl_deltaA = @(yi, yo) ((0.5*yo^2 + (2/3)*yo^3*(taper-1)/b)-(0.5*yi^2 + ...
    (2/3)*yi^3*(taper-1)/b))*(2*CLw_slope*tau*Cr)/(Sw*b); %Cl derivative
eq_Cl = @(x) x*delta_max(1,1)*pi/180; %Cl - rolling moment coefficient 

eq_LA = @(Cl) 0.5*rho*V^2*Sw*Cl*b; %Lift contribution from aileron 

eq_Pss = @(LA) sqrt((2*LA)/(rho*(Sw+Sv+Sh)*C_DR*yD^3));

eq_phi1 = @(Pss) (Ixx)*(log(Pss^2))/(rho*yD^3 * (Sw+Sv+Sh)*C_DR);

eq_roll_rate_derivative = @(phi1, Pss) 0.5*Pss^2/phi1;

%Assume that phi1 is going to be larger than phi_req for GA
eq_t2 = @(phi, P_dot) sqrt(2*phi/P_dot);

%Testing the equation 
Cl_deltaA = eq_Cl_deltaA(bi, bo);
Cl = eq_Cl(Cl_deltaA);
LA = eq_LA(Cl);
Pss = eq_Pss(LA);
phi1 = eq_phi1(Pss);
P_dot = eq_roll_rate_derivative(phi1, Pss);
tss = sqrt(2*phi1/P_dot);

if (phi1 > phi_req)
    t2 = eq_t2(phi_req*pi/180, P_dot);
else
    delta_t = (phi_req - phi1)/Pss;
    t2 = tss + delta_t;
end
err = (t2-t_req)/t_req *100;

%Plotting some points
time_v_phi(0.01, 5, tss,P_dot,Pss,phi1);

function graph_t_v_phi = time_v_phi(time_step,max_time,tss,P_dot,Pss,phi1)
        %Time vs phi
    t_nonlinear = 0:time_step:tss;
    phi = 0.5.*t_nonlinear.^2 .* P_dot;
    plot(t_nonlinear, phi);
    title("Time vs. Bank Position");
    xlabel('Time (s)');
    ylabel('Bank angle (deg)');
    hold on
    %Linear region
    t_linear = tss:time_step:max_time;
    b = phi1 - Pss*tss;
    phi_lin = b + Pss.*t_linear;
    plot(t_linear, phi_lin);
    hold off
end

function tau = interpolate_tau(x_points, y_points, x)
    data_x = x_points(:);
    data_y = y_points(:);

    if x < min(data_x) || x > max(data_x)
        error('Outside the range!');
    end

    %Find indices of 2 points, 1st point less than x, 2nd point greater
    indx1 = find(data_x <= x, 1, 'last');
    indx2 = find(data_x >= x, 1, 'first');

    %If x matches a data point, return corresponding y
    if data_x(indx1) == x_points
        tau = data_y(indx1);
    end

    %Linear interpolation 
    x1 = data_x(indx1);
    x2 = data_x(indx2);
    y1 = data_y(indx1);
    y2 = data_y(indx2);

    tau = y1 + (y2 - y1) * (x - x1) / (x2 - x1);
end