theta_beta_mach_relation

%% NIR KALUSH
%% Implementation of Equation 9.23 from Anderson's Fundamentals of Aerodynamics 
%% ALL RIGHTS RESERVED to Mr. Kalush, University of Maryland, College Park

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NOTE!!!!!! YOU WILL NEED THE SYMBOLIC TOOLBOX TO RUN THIS PROGRAM. Check to see you have it installed  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function theta_beta_mach_relation
%% This program accepts 2 of the 3 values: sigma, beta, Mach values, and returns the third
fprintf('please enter the unknown value to solve for:\n\n 1. Sigma \t 2. beta \t 3. Mach number\n\n');
variable = input('');

if variable == 1
    B = input('Enter Beta value in degrees: '); fprintf('\n');
    M = input('Enter Mach number: '); fprintf('\n');
    T = (180/pi)*atan(2*cot(B*pi/180)*((M*sin(B*pi/180))^2-1)/(M^2*(1.4+cos(2*B*pi/180))+2));
    if (T < 0)
        fprintf('The theta for the given Mach = %2g and Beta = %2g vals falls below zero: %f degrees - unreasonable',M,B,T);
    else
        fprintf('The theta for the given Mach = %2g and Beta = %2g vals is: \t %f degrees',M,B,T);
    end
    
elseif variable == 2
    T = input('Enter Theta value in degrees: '); fprintf('\n');
    M = input('Enter Mach number: '); fprintf('\n');
    LHS = tan(T*pi/180);
    B = 1;
    RHS = inline('2*cot(B*pi/180)*((M*sin(B*pi/180))^2-1)/(M^2*(1.4+cos(2*B*pi/180))+2)','B','M')
    rhs_beta_1_mach_given = RHS(1,M)
    counter = 0;
    while (RHS(B,M) < LHS*0.98 | RHS(B,M) > LHS*1.02)
        if (B > 90 | B < 0 | counter >200)
            break;
        end
        if (RHS(B,M) < LHS*0.70)
            B=B+1;
        elseif(RHS(B,M) < LHS*0.90)
            B=B+0.5;
        elseif(RHS(B,M) < LHS*0.95)
            B=B+0.3;
        elseif(RHS(B,M) < LHS*0.98)
            B=B+0.1;
        end
        
        if(RHS(B,M) > LHS*1.03)
            B=B-0.1;
        elseif(RHS(B,M) > LHS*1.10)
            B = B-0.5;
        end
        counter = counter+1;
        lhs_curr = LHS;
        rhs_curr = RHS(B,M);
    
    end
    
    counter = 0;
    
    while (RHS(B,M) < LHS*0.99 | RHS(B,M) > LHS*1.01)
        if (B > 90 | B < 0 | counter > 100)
            break;
        end
        if (RHS(B,M) < LHS*0.98)
            B=B+0.05;
        elseif(RHS(B,M) < LHS*0.985)
            B=B+0.03;
        elseif(RHS(B,M) < LHS*0.99)
            B=B+0.01;
        end
        
        if(RHS(B,M) > LHS*1.01)
            B=B-0.01;
        elseif(RHS(B,M) > LHS*1.02)
            B = B-0.1;
        end
    
        counter = counter+1;
        lhs_curr = LHS;
        rhs_curr = RHS(B,M);
    end
    
    counter = 0;
    
    while (RHS(B,M) ~= LHS)
        if (B > 90 | B < 0 | counter > 40)
            break;
        end
        
        if (RHS(B,M) < LHS)
            B=B+0.005;
        end
        if(RHS(B,M) > LHS)
            B=B-0.005;
        end
        counter = counter+1;
        lhs_curr = LHS;
        rhs_curr = RHS(B,M);
    end
    
    if (B < 0)
        fprintf('The Beta angle for the given Mach = %2g and Theta = %2g vals is not in bounds (0,90): B = %2g degrees',M,T,B);
    elseif (B > 90)
        fprintf('The Shock is detached from the body \n\n for the given Mach = %2g and Theta = %2g \t B = %2g degrees - a useless quantity',M,T,B);
    else
        fprintf('The LHS of equation 9.23 yields %f. The RHS approximation yields %f \n\n',lhs_curr,rhs_curr);
        fprintf('The approximation was within %2g %% of LHS \n\n',100*abs(lhs_curr-rhs_curr)/lhs_curr);
        fprintf('The Beta angle for the given Mach = %2g and Theta = %2g vals is: B = %2g degrees',M,T,B);
    end
elseif variable == 3
    B = input('Enter Beta value in degrees: '); fprintf('\n');
    T = input('Enter Theta value in degrees: '); fprintf('\n');
    
    LHS = tan(T*pi/180);
    M=0;
    RHS = inline('2*cot(B*pi/180)*((M*sin(B*pi/180))^2-1)/(M^2*(1.4+cos(2*B*pi/180))+2)','M','B')
    rhs_mach_0_beta_given = RHS(0,B)
    counter = 0;
    while (RHS(M,B) < LHS*0.98 | RHS(M,B) > LHS*1.02)
        if (M > 20 | counter >200)
            break;
        end
        if (RHS(M,B) < LHS*0.85)
            M=M+0.1;
        elseif(RHS(M,B) < LHS*0.90)
            M=M+0.05;
        elseif(RHS(M,B) < LHS*0.95)
            M=M+0.01;
        elseif(RHS(M,B) < LHS*0.98)
            M=M+0.005;
        end
        
        if(RHS(M,B) > LHS*1.02)
            M=M-0.005;
        elseif(RHS(M,B) > LHS*1.10)
            M = M-0.1;
        end
        counter = counter+1;
        lhs_curr = LHS;
        rhs_curr = RHS(M,B);
    
    end
    
    counter = 0;
    
    while (RHS(M,B) < LHS*0.99 | RHS(M,B) > LHS*1.01)
        if (M > 20 | counter > 100)
            break;
        end
        if (RHS(M,B) < LHS*0.98)
            M=M+0.005;
        elseif(RHS(M,B) < LHS*0.985)
            M=M+0.001;
        elseif(RHS(M,B) < LHS*0.99)
            M=M+0.0001;
        end
        
        if(RHS(M,B) > LHS*1.01)
            M=M-0.0005;
        elseif(RHS(M,B) > LHS*1.02)
            M = M-0.001;
        end
    
        counter = counter+1;
        lhs_curr = LHS;
        rhs_curr = RHS(M,B);
    end
    
    counter = 0;
    
    while (RHS(M,B) ~= LHS)
        if (M>20 | counter > 40)
            break;
        end
        
        if (RHS(M,B) < LHS)
            M=M+0.0001;
        end
        if(RHS(M,B) > LHS)
            M=M-0.0001;
        end
        counter = counter+1;
        lhs_curr = LHS;
        rhs_curr = RHS(M,B);
    end
    
    fprintf('The LHS of equation 9.23 yields %f. The RHS approximation yields %f \n\n',lhs_curr,rhs_curr);
    fprintf('The approximation was within %2g %% of LHS \n\n',100*abs(lhs_curr-rhs_curr)/lhs_curr);
    fprintf('The Mach number for the given Beta = %2g and Theta = %2g vals is: M = %2g',B,T,M);
end