Removed atmosisa dependence. Added automation scripts for trajectory runs.

.gitignore

@@ -0,0 +1,4 @@
+*.asv
+Code/ICANSAT/V0/Result
+Code/ICANSAT/V1/Result
+Code/MAD/Result

Code/AllClean.m

@@ -0,0 +1,17 @@
+mfile_name          = mfilename('fullpath');
+[pathstr,name,ext]  = fileparts(mfile_name);
+cd(pathstr);
+
+cd ICANSAT/V0/Result
+delete *.txt
+fprintf("ICANSAT V0 result files deleted.\n");
+
+cd ../../V1/Result
+delete *.txt
+fprintf("ICANSAT V1 result files deleted.\n");
+
+cd ../../../MAD/Result
+delete *.txt
+fprintf("MAD result files deleted.\n");
+
+cd(pathstr);

Code/AllRun.m

@@ -0,0 +1,14 @@
+mfile_name          = mfilename('fullpath');
+[pathstr,name,ext]  = fileparts(mfile_name);
+cd(pathstr);
+
+cd ICANSAT/V0
+ICANSAT
+
+cd ../V1
+ICANSAT_V1
+
+cd ../../MAD
+MAD_Validation
+
+cd(pathstr);

Code/CalcAtmosQuantities.m

@@ -0,0 +1,70 @@
+function [t,p,rho] = CalcAtmosQuantities(h)
+    hg = h;               % geometric alt in meters
+    R = 287;              % Gas Constant
+    re = 6.356766e6;      % radius of earth
+    g0 = 9.8066;          % gravity
+    h = re/(re+hg)*hg;
+
+    % Altitude in meters
+    h0 = 0.0;
+    h1 = 11000.0;
+    h2 = 25000.0;
+    h3 = 47000;
+    h4 = 53000;
+    h5 = 79000;
+    h6 = 90000;
+    h7 = 105000;
+
+    % Temperature in Kelvin.
+    t0 = 288.16;
+    t1 = 216.66;
+    t2 = 216.66;
+    t3 = 282.66;
+    t4 = 282.66;
+    t5 = 165.66;
+    t6 = 165.66;
+    t7 = 225.66;
+
+    %%equations
+    a01 = (t1-t0)/(h1-h0);  % slope of T/h in K/m for 0 to 11 km
+    a23 = (t3-t2)/(h3-h2);  % slope of T/h in K/m for 25 to 47 km
+    a45 = (t5-t4)/(h5-h4);  % slope of T/h in K/m for 53 to 79 km
+    a67 = (t7-t6)/(h7-h6);  % slope of T/h in K/m for 90 to 105 km
+
+    % Pressure in N/m^2
+    p0 = 101325.0;
+    p1 = p0*(t1/t0)^(-g0/(a01*R));
+    p2 = p1*exp(-(g0/(R*t1))*(h2-h1));
+    p3 = p2*(t3/t2)^(-g0/(a23*R));
+    p4 = p3*exp(-(g0/(R*t3))*(h4-h3));
+    p5 = p4*(t5/t4)^(-g0/(a45*R));
+    p6 = p5*exp(-(g0/(R*t5))*(h6-h5));
+    p7 = p6*(t7/t6)^(-g0/(a67*R));
+
+    if h < h1
+        t = t0+a01*h;
+        p = p0*(t/t0)^(-g0/(a01*R));
+    elseif h < h2
+        t = t1;
+        p= p1*exp(-(g0/(R*t))*(h-h1));
+    elseif h < h3
+        t = t2+a23*(h-h2);
+        p = p2*(t/t2)^(-g0/(a23*R));
+    elseif h < h4
+        t = t3;
+        p = p3*exp(-(g0/(R*t))*(h-h3));
+    elseif h < h5
+        t = t4+a45*(h-h4);
+        p = p4*(t/t4)^(-g0/(a45*R));
+    elseif h < h6
+        t = t5;
+        p = p5*exp(-(g0/(R*t))*(h-h5));
+    elseif h < h7
+        t = t6+a67*(h-h6);
+        p = p6*(t/t6)^(-g0/(a67*R));
+    else
+        t= t7;
+        p = p7*exp(-(g0/(R*t))*(h-h7));
+    end
+    rho=p/(R*t); % Density in kg/m^3
+end

Code/DragInterpolations_ICANSAT.m

@@ -1,7 +1,5 @@
 function [coeffdrag] = DragInterpolations_ICANSAT(velocity, altitude, time, tb)
 
-% global tb
-
 drag_data = dlmread("ICANSAT/ICANSAT_DragData.txt", "\t",1,0);
 
 % Mach No
@@ -13,7 +11,13 @@ cd_full = drag_data(:,3);
 % No base drag
 cd_zero = drag_data(:,2);
 
-[temp, sound_speed, P, rho] = atmosisa(altitude);
+R = 8.3145;
+gamma = 1.4;
+M_air = 0.0289645;
+
+[T,p,rho] = CalcAtmosQuantities(altitude);
+sound_speed = sqrt(gamma*R*T/M_air);
+
 mach_speed = velocity/sound_speed;
 
 if (time > tb)

Code/DragInterpolations_MAD.m

@@ -27,7 +27,13 @@ cdzero_10000 = drag_data(:,3);
 
 % Interpolate for Cd values at two stated altitudes, then interpolate for
 % Cd value at a specific altitude
-[temp, sound_speed, P, rho] = atmosisa(altitude);
+
+R = 8.3145;
+gamma = 1.4;
+M_air = 0.0289645;
+
+[T,p,rho] = CalcAtmosQuantities(altitude);
+sound_speed = sqrt(gamma*R*T/M_air);
 mach_speed = velocity/sound_speed;
 
 % If propellant is still available

Code/DragInterpolations_Spark.m

@@ -2,7 +2,13 @@ function [coeffdrag] = DragInterpolations_Spark(velocity, altitude)
 
 spark_drag_data = dlmread("Spark/Spark_DragData.txt", "\t",1,0); 
 
-[temp, sound_speed, P, rho] = atmosisa(altitude);
+R = 8.3145;
+gamma = 1.4;
+M_air = 0.0289645;
+
+[T,p,rho] = CalcAtmosQuantities(altitude);
+sound_speed = sqrt(gamma*R*T/M_air);
+
 mach_speed = velocity/sound_speed;
 
 mach = spark_drag_data(:,1);

Code/ICANSAT/V1/ICANSAT_V1.m

@@ -5,15 +5,15 @@ addpath('../..')
 % Df, Lr, Ls, ms, mp, tb
 icansat = SingleStageRocket(0.1567, 2.391, 4, 14, 1.160534, 4.0);
 
-% RunTrajectorySim(icansat, 88);
-% RunTrajectorySim(icansat, 86);
-% RunTrajectorySim(icansat, 84);
-% RunTrajectorySim(icansat, 82);
-% RunTrajectorySim(icansat, 80);
-% RunTrajectorySim(icansat, 78);
-% RunTrajectorySim(icansat, 76);
-% RunTrajectorySim(icansat, 74);
-% RunTrajectorySim(icansat, 72);
+RunTrajectorySim(icansat, 88);
+RunTrajectorySim(icansat, 86);
+RunTrajectorySim(icansat, 84);
+RunTrajectorySim(icansat, 82);
+RunTrajectorySim(icansat, 80);
+RunTrajectorySim(icansat, 78);
+RunTrajectorySim(icansat, 76);
+RunTrajectorySim(icansat, 74);
+RunTrajectorySim(icansat, 72);
 RunTrajectorySim(icansat, 70);
 
 function RunTrajectorySim(rocket, launch_angle)

Code/SingleStageRocket.m

@@ -165,81 +165,9 @@ classdef SingleStageRocket < handle
            obj.CalculateMassFlow(length(resampled_time));
        end
        
-       function [t,p,rho] = CalcAtmosQuantities(obj, h)
-            hg = h;               % geometric alt in meters
-            R = 287;              % Gas Constant
-            re = 6.356766e6;      % radius of earth
-            g0 = 9.8066;          % gravity
-            h = re/(re+hg)*hg;
-
-            % Altitude in meters
-            h0 = 0.0;
-            h1 = 11000.0;
-            h2 = 25000.0;
-            h3 = 47000;
-            h4 = 53000;
-            h5 = 79000;
-            h6 = 90000;
-            h7 = 105000;
-
-            % Temperature in Kelvin.
-            t0 = 288.16;
-            t1 = 216.66;
-            t2 = 216.66;
-            t3 = 282.66;
-            t4 = 282.66;
-            t5 = 165.66;
-            t6 = 165.66;
-            t7 = 225.66;
-
-            %%equations
-
-            a01 = (t1-t0)/(h1-h0);  % slope of T/h in K/m for 0 to 11 km
-            a23 = (t3-t2)/(h3-h2);  % slope of T/h in K/m for 25 to 47 km
-            a45 = (t5-t4)/(h5-h4);  % slope of T/h in K/m for 53 to 79 km
-            a67 = (t7-t6)/(h7-h6);  % slope of T/h in K/m for 90 to 105 km
-
-            % Pressure in N/m^2
-            p0 = 101325.0;
-            p1 = p0*(t1/t0)^(-g0/(a01*R));
-            p2 = p1*exp(-(g0/(R*t1))*(h2-h1));
-            p3 = p2*(t3/t2)^(-g0/(a23*R));
-            p4 = p3*exp(-(g0/(R*t3))*(h4-h3));
-            p5 = p4*(t5/t4)^(-g0/(a45*R));
-            p6 = p5*exp(-(g0/(R*t5))*(h6-h5));
-            p7 = p6*(t7/t6)^(-g0/(a67*R));
-
-            if h < h1
-                t = t0+a01*h;
-                p = p0*(t/t0)^(-g0/(a01*R));
-            elseif h < h2
-                t = t1;
-                p= p1*exp(-(g0/(R*t))*(h-h1));
-            elseif h < h3
-                t = t2+a23*(h-h2);
-                p = p2*(t/t2)^(-g0/(a23*R));
-            elseif h < h4
-                t = t3;
-                p = p3*exp(-(g0/(R*t))*(h-h3));
-            elseif h < h5
-                t = t4+a45*(h-h4);
-                p = p4*(t/t4)^(-g0/(a45*R));
-            elseif h < h6
-                t = t5;
-                p = p5*exp(-(g0/(R*t))*(h-h5));
-            elseif h < h7
-                t = t6+a67*(h-h6);
-                p = p6*(t/t6)^(-g0/(a67*R));
-            else
-                t= t7;
-                p = p7*exp(-(g0/(R*t))*(h-h7));
-            end
-            rho=p/(R*t); % Density in kg/m^3
-       end
-       
        % Thrust correction based on ASTOS manual
        function T = ThrustCorrection(obj, h, thrust)
-           [t,p,rho] = obj.CalcAtmosQuantities(h);
+           [t,p,rho] = CalcAtmosQuantities(h);
            T = thrust*1.095 - 0.014 * p;
        end
        
@@ -277,7 +205,7 @@ classdef SingleStageRocket < handle
        
        % Density calculation
        function [T, rho] = CalculateRhoTemp(obj, h)
-          [T, p, rho] = obj.CalcAtmosQuantities(h);
+          [T,p,rho] = CalcAtmosQuantities(h);
           return;
        end
        
@@ -427,15 +355,18 @@ classdef SingleStageRocket < handle
        % Write data file
        function WriteDataFile(obj)
            if obj.mdot_constant == true
-               result_dat_file = obj.ID + "_Result_" + obj.launch_angle + "deg.txt";
+               result_dat_file = 'Result/' + obj.ID + "_Result_" + obj.launch_angle + "deg.txt";
            else
                text = "VariableMassFlow";
-               result_dat_file = obj.ID + "_Result_" + obj.launch_angle + "_deg_" + text + ".txt";
+               result_dat_file = 'Result/' + obj.ID + "_Result_" + obj.launch_angle + "_deg_" + text + ".txt";
            end
            result_matrix = [obj.t obj.x/1000 obj.y/1000 obj.M obj.thrust obj.DvDt obj.m obj.mdot obj.isp];
            result_data = reshape(result_matrix, [], 9);
            dlmwrite(result_dat_file, 'Time(s),x(km),y(km),Mach,Thrust(N),Acceleration(m/s2),Mass(kg),mdot(kg/s),Isp (s)','delimiter','')
            dlmwrite(result_dat_file, result_data, '-append', 'delimiter', ',');
+           notify_text = 'Results written for ' + obj.ID + ' at ' + obj.launch_angle + " degrees.\n ";
+           fprintf(notify_text);
+           fprintf("-----------------------------------------------------------------------------------\n");
        end
    end
 end