implemented improved thrust curve

.vscode/c_cpp_properties.json

@@ -0,0 +1,22 @@
+{
+    "configurations": [
+        {
+            "name": "Win32",
+            "includePath": [
+                "${workspaceFolder}/**",
+                "${workspaceFolder}/include"
+            ],
+            "defines": [
+                "_DEBUG",
+                "UNICODE",
+                "_UNICODE"
+            ],
+            "windowsSdkVersion": "10.0.18362.0",
+            "compilerPath": "C:\\Program Files (x86)\\Microsoft Visual Studio\\2019\\Community\\VC\\Tools\\MSVC\\14.25.28610\\bin\\Hostx64\\x64\\cl.exe",
+            "cStandard": "c17",
+            "cppStandard": "c++17",
+            "intelliSenseMode": "windows-msvc-x64"
+        }
+    ],
+    "version": 4
+}

build/debug/simPlot.m

@@ -59,18 +59,26 @@ figure(2)
 
 % Euler Angles
 subplot(2, 1, 1)
-plot(t, yaw, pitch, roll)
+hold on;
+plot(t, yaw)
+plot(t, pitch)
+plot(t, roll)
 title('Euler Angles vs Time')
 xlabel('Time (ms)')
 ylabel('Euler Angles (deg)')
+legend("yaw", "pitch", "roll")
 
 % Angular Velocity
 subplot(2, 1, 2)
-plot(t, yawdot, pitchdot, rolldot)
+hold on;
+plot(t, yawdot)
+plot(t, pitchdot)
+plot(t, rolldot)
 title('Angular Velocity vs Time')
 xlabel('Time (ms)')
 ylabel('Angular Velocity (deg/s)')
 saveas(gcf,'outputs/Euler Angles vs Time.png')
+legend("yawdot", "pitchdot", "rolldot")
 
 figure(3)
 

include/sVars.h

@@ -38,33 +38,6 @@ struct sVars
 
     double simTime;
     int stepSize;
-
-    std::array<std::array<double, 2>, 26> thrustCurve = {{{0.148, 7.638},
-                                                          {0.228, 12.253},
-                                                          {0.294, 16.391},
-                                                          {0.353, 20.210},
-                                                          {0.382, 22.756},
-                                                          {0.419, 25.260},
-                                                          {0.477, 23.074},
-                                                          {0.520, 20.845},
-                                                          {0.593, 19.093},
-                                                          {0.688, 17.500},
-                                                          {0.855, 16.225},
-                                                          {1.037, 15.427},
-                                                          {1.205, 14.948},
-                                                          {1.423, 14.627},
-                                                          {1.452, 15.741},
-                                                          {1.503, 14.785},
-                                                          {1.736, 14.623},
-                                                          {1.955, 14.303},
-                                                          {2.210, 14.141},
-                                                          {2.494, 13.819},
-                                                          {2.763, 13.338},
-                                                          {3.120, 13.334},
-                                                          {3.382, 13.013},
-                                                          {3.404, 9.352},
-                                                          {3.418, 4.895},
-                                                          {3.450, 0.000}}};
 };
 
 #endif

src/main.cpp

@@ -14,7 +14,7 @@ int main()
     Vars.vz = 0; // [m/s]
 
     // Initial YPR
-    Vars.yaw = 15 * 3.1416 / 180;  // [rad]
+    Vars.yaw = 0 * 3.1416 / 180;   // [rad]
     Vars.pitch = 0 * 3.1416 / 180; // [rad]
     Vars.roll = 0 * 3.1416 / 180;  // [rad]
 
@@ -36,10 +36,10 @@ int main()
     Vars.stepSize = 1; // [ms]
 
     // Other Properties
-    Vars.mp = 0.06;                                             // [kg]
+    Vars.mp = 0.06;                // [kg]
-    Vars.mb = Vars.m0 - Vars.mp;                                // [kg]
+    Vars.mb = Vars.m0 - Vars.mp;   // [kg]
-    Vars.tb = Vars.thrustCurve[25][0] - Vars.thrustCurve[0][0]; // [s]
+    Vars.tb = 3.45 - 0.148;        // [s]
-    Vars.mdot = Vars.mp / Vars.tb;                              // [kg/s]
+    Vars.mdot = Vars.mp / Vars.tb; // [kg/s]
 
     sim(Vars);
 

src/sim.cpp

@@ -47,38 +47,30 @@ void sim(struct sVars &Vars)
 
 void burnStartTimeCalc(struct sVars &Vars, double g)
 {
-    double v_prev = Vars.vz;
+    double v = Vars.vz;
-    double h_prev = 0;
+    double h = 0;
-    double dt, a, v, h;
+    double dt = 0.001;
+    double a, j, m, thrust;
 
-    for (int i = 0; i < 25; i++)
+    for (double i = 0.148; i < 3.450; i = i + dt)
     {
-        // Integral of thrust curve
+        m = Vars.m0 - i * Vars.mdot;
-        if (i < 25)
-        {
-            dt = Vars.thrustCurve[i + 1][0] - Vars.thrustCurve[i][0];
-            a = (Vars.thrustCurve[i][1] /
-                     (Vars.m0 - Vars.mdot * Vars.thrustCurve[i + 1][0]) +
-                 g);
-        }
-        else
-        {
-            dt = 0;
-            a = (Vars.thrustCurve[i][1] / Vars.mb) + g;
-        }
 
-        v = a * dt + v_prev;
+        if ((i > 0.147) & (i < 0.420))
-        v_prev = v;
+            thrust = 65.165 * i - 2.3921;
-        h = v * dt + h_prev;
+
-        h_prev = h;
+        else if ((i > 0.419) & (i < 3.383))
+            thrust = 0.8932 * pow(i, 6) - 11.609 * pow(i, 5) + 60.739 * pow(i, 4) - 162.99 * pow(i, 3) + 235.6 * pow(i, 2) - 174.43 * i + 67.17;
+
+        else if ((i > 3.382) & (i < 3.46))
+            thrust = -195.78 * i + 675.11;
+
+        v = (((thrust / m) + g) * dt) + v;
+        h = v * dt + h;
     }
 
-    double hb = h; // height that we want to start our burn at
+    Vars.z = h + (pow(v, 2) / (2 * -g)); // starting height
-    double hf =
+    Vars.vb = v;                         // terminal velocity
-        (v * v) / (2 * -g); // amount of height that we need to be in freefall in
-                            // order to reach terminal velocity
-    Vars.z = hb + hf;       // starting height
-    Vars.vb = v;            // terminal velocity
 
     double burnStartTime = Vars.vb / -g;
     Vars.simTime = (Vars.tb + burnStartTime) * 1000;
@@ -110,26 +102,14 @@ void thrustSelection(struct sVars &Vars, int t)
     else
         Vars.burnElapsed = 2000; // arbitrary number to ensure we don't burn
 
-    bool b_thrustSelect;
+    if ((Vars.burnElapsed > 0.147) & (Vars.burnElapsed < 0.420))
+        Vars.thrust = 65.165 * Vars.burnElapsed - 2.3921;
 
-    for (int i = 0; i < 25; i++)
+    else if ((Vars.burnElapsed > 0.419) & (Vars.burnElapsed < 3.383))
-    {
+        Vars.thrust = 0.8932 * pow(Vars.burnElapsed, 6) - 11.609 * pow(Vars.burnElapsed, 5) + 60.739 * pow(Vars.burnElapsed, 4) - 162.99 * pow(Vars.burnElapsed, 3) + 235.6 * pow(Vars.burnElapsed, 2) - 174.43 * Vars.burnElapsed + 67.17;
-        // Compare elapsed burn time to each time entry in the thrust curve
-        b_thrustSelect = (Vars.burnElapsed < (1 + tol) * Vars.thrustCurve[i][0]) &
-                         (Vars.burnElapsed > (1 - tol) * Vars.thrustCurve[i][0]);
 
-        if (b_thrustSelect)
+    else if ((Vars.burnElapsed > 3.382) & (Vars.burnElapsed < 3.46))
-        {
+        Vars.thrust = -195.78 * Vars.burnElapsed + 675.11;
-            // Choose thrust associated with time entry that elapsed burn time was
-            // compared to above
-            Vars.thrust = Vars.thrustCurve[i][1];
-            Vars.thrust_prev = Vars.thrustCurve[i][1];
-        }
-
-        else
-            // latch
-            Vars.thrust = Vars.thrust_prev;
-    }
 }
 
 void lqrCalc(struct sVars &Vars)