Rework of 2D arrays, implemented vector math library

Vector math library now converts from perifocal frame to ECI frame for propagation

c_code/makefile

@@ -8,4 +8,7 @@ debug:
 	gcc -o main main.c -O3 -Wall -g -Wall -lm
 
 python:
-	gcc -fPIC -shared -o ../python_code/supernova.so main.c -O3 -lm -fopt-info-vec
+	gcc -fPIC -shared -o ../python_code/supernova.so main.c -O3 -lm -fopt-info-vec
+
+test:
+	gcc -o main testvec.c -Wall -lm

c_code/orbit.h

@@ -35,12 +35,28 @@ void orbit(char solver[], double tSpan[], double y0[], char output[], double ATO
 
     for (int i = 0; i < result->n; i++) {
         fprintf(fptr, "%.15f, %.15f, %.15f, %.15f, %.15f, %.15f, %.15f\n", result->t[i], result->y[i][0], result->y[i][1], result->y[i][2], result->y[i][3], result->y[i][4], result->y[i][5]);
-        free(result->y[i]);
     }
     fclose(fptr);
 
     free(result->t);
     free(result);
 }
 
+solution* orbitPTR(char solver[], double tSpan[], double y0[], double ATOL) {
+    /*
+    Propagates orbit using solver of choice and returns pointer to solution
+    solver: name, either "RK810" or "RK1012"
+    tSpan: t0, tf of times.
+    y0: initial x y z vx vy vz state vector
+    ATOL: tolerance
+    */
+    solution* result;
+    if (!strcmp(solver, "RK810")) result = RK810vec(combined_perturbations, tSpan, y0, ATOL);
+    else if (!strcmp(solver, "RK1012")) result = RK1012vec(combined_perturbations, tSpan, y0, ATOL);
+    else {
+        printf("Invalid solver chosen.\n");
+    }
+    return result;
+}
+
 #endif

c_code/solvers.h

@@ -8,7 +8,7 @@
 
 typedef struct AdaptiveSolution {
     double* t; // timesteps
-    double** y; // function value
+    double (*y)[VEC_SIZE]; // function value
     int n; // number of steps taken
 } solution;
 
@@ -28,15 +28,14 @@ solution* RK810vec(void (*f)(double, double[], double*), double tSpan[], double
 
     if (tSpan[1] <= tSpan[0]) return NULL; // error case
 
-    ////// Initial estimate for n
+    ////// Initial estimate for n = number of rows in solution
     int n = (int)((tSpan[1] - tSpan[0])/H0810); // 1 step every 10 seconds
     int step = 0;
 
     ////// Allocate result struct
     solution* result = (solution*) malloc(sizeof(solution));
-    result->y = malloc(n * sizeof(double*));
+    result->y = malloc(n * sizeof(double[VEC_SIZE]));
     result->t = malloc(n * sizeof(double));
-    for (int i = 0; i < n; i++) result->y[i] = malloc(VEC_SIZE * sizeof(double));
 
     ////// Prepare Step variables
     double h = H0810; // initial stepsize
@@ -110,10 +109,9 @@ solution* RK810vec(void (*f)(double, double[], double*), double tSpan[], double
             /// step within tolerance, append solution
             // check array size and increase size if needed
             if ((step + 1) >= n) {
-                n *= 2;
-                result->y = realloc(result->y, n * sizeof(double*));
+                n *= 2; // double size of array and reallocate memory
+                result->y = realloc(result->y, n * sizeof(double[VEC_SIZE]));
                 result->t = realloc(result->t, n * sizeof(double));
-                for (int i = step + 1; i < n; i++) result->y[i] = malloc(VEC_SIZE * sizeof(double));
             }
 
             // Append to result
@@ -156,9 +154,9 @@ solution* RK1012vec(void (*f)(double, double[], double*), double tSpan[], double
 
     ////// Allocate result struct
     solution* result = (solution*) malloc(sizeof(solution));
-    result->y = malloc(n * sizeof(double*));
+    result->y = malloc(n * sizeof(double[VEC_SIZE]));
     result->t = malloc(n * sizeof(double));
-    for (int i = 0; i < n; i++) result->y[i] = malloc(VEC_SIZE * sizeof(double));
+
 
     ////// Prepare Step variables
     double h = H01012; // initial stepsize
@@ -239,10 +237,9 @@ solution* RK1012vec(void (*f)(double, double[], double*), double tSpan[], double
             /// step within tolerance, append solution
             // check array size and increase size if needed
             if ((step + 1) >= n) {
-                n *= 2;
-                result->y = realloc(result->y, n * sizeof(double*));
+                n *= 2; // double size of array and reallocate memory
+                result->y = realloc(result->y, n * sizeof(double[VEC_SIZE]));
                 result->t = realloc(result->t, n * sizeof(double));
-                for (int i = step + 1; i < n; i++) result->y[i] = malloc(VEC_SIZE * sizeof(double));
             }
 
             // Append to result

c_code/tables/2D_array_format.md

@@ -0,0 +1,36 @@
+# A Note on Dynamic 2D Array Formatting
+
+Source: https://stackoverflow.com/questions/40847172/best-way-to-allocate-memory-to-a-two-dimensional-array-in-c
+
+2D arrays are most efficiently stored as a contiguous block, i.e. as a single linear array.
+
+Therefore, the array
+```c
+[1, 2, 3]
+[4, 5, 6]
+[7, 8, 9]
+```
+has memory allocated like the following:
+`[1, 2, 3, 4, 5, 6, 7, 8, 9]`
+
+## Array Declaration
+Arrays of doubles should be declared like this:
+```c
+double (*name)[COLUMN_SIZE];
+```
+Wherein `COLUMN_SIZE` is some constant denoting the column length of the array. For Supernova, we use a column length of 6 (x/y/z/vx/vy/vz)
+
+This way of formatting the array makes it so that the first index of the array yields chunks spaced out by `COLUMN_SIZE`.
+
+## Array Memory Allocation
+Memory allocation can be done like the following:
+```c
+name = malloc(n * sizeof(double[COLUMN_SIZE]));
+```
+Wherein `n` is the number of rows in the array. Similar procedures are done for reallocating memory, if you need to make the array bigger.
+
+## Accessing Array indices
+When the array is initialized like this, you can access it just as you would expect with Python.
+```c
+printf("%f", name[1][2]);
+>>> 6

c_code/testvec.c

@@ -0,0 +1,58 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "vecmath.h"
+#include <math.h>
+
+int main() {
+    // Validation from Mathworks (MATLAB documentation)
+
+    // 1. CROSS PRODUCT
+    double A1[3] = {4, -2, 1};
+    double B1[3] = {1, -1, 3};
+    double C1[3];
+    cross(A1, B1, C1);
+    printVec(C1);
+    // Result: [-5, -11, -2]
+
+    // 2. DOT PRODUCT
+    printf("%f\n", dot(A1, C1));
+    // Result: 0
+    printf("%f\n", dot(B1, C1));
+    // Result: 0
+
+    double A2[3] = {4, -1, 2};
+    double B2[3] = {2, -2, -1};
+    printf("%f\n", dot(A2, B2));
+    // Result: 8
+
+    // 3. Matrices
+    double A3[3][3] = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
+    double B3[3][3] = {{0, 1, 0}, {1, 0, 0}, {0, 0, 1}};
+    double (*C3)[3] = malloc(3 * sizeof(double[3]));
+
+    matXmat(A3, B3, C3);
+    printMatrix(C3);
+
+    double (*q)[3] = QpX(0, 97.49588373 * PI/180.0, 227.05566156 * PI/180.0);
+    printMatrix(q);
+
+    // 4. Matrix times Vec
+    double A4[3][3] = {{-0.681288, -0.095495, -0.725760}, {-0.732016, 0.088877, 0.675465}, {0.000000, 0.991454, -0.130455}};
+    double B4[3] = {6877999.998558, 0.000000, 0.000000};
+    double C4[3];
+    matXvec(A4, B4, C4);
+    printVec(C4);
+
+    // 5. Orbit alteration
+    double* Perifocal_pos_vel = PFE(6903e3, 0.003621614, 97.49588373 * PI/180.0, 0, 0);
+    // Perifocal coordinates
+
+    double ECI_position[3];
+    double ECI_vel[3];
+    matXvec(q, Perifocal_pos_vel, ECI_position);
+    matXvec(q, Perifocal_pos_vel+3, ECI_vel);
+
+    printVec(ECI_position);
+    printVec(ECI_vel);
+
+}

c_code/vecmath.h

@@ -3,51 +3,38 @@
 
 // vector math library
 #include <math.h>
+#include <stdio.h>
+#define PI 	3.14159265358979323846
+#define GM 3.986004405e14
 
-void cross(double a[3], double b[3], double u[3]) {
-    // Crosses a, b and stores result in u
-    // 231 312
-    // 120 201
-    u[0] = a[1]*b[2] - b[1]*a[2];
-    u[1] = a[2]*b[0] - b[2]*a[0];
-    u[2] = a[0]*b[1] - b[0]*a[1];
+void cross(double A[3], double B[3], double U[3]) {
+    // Crosses vectors A and B then stores result in vector U
+    U[0] = A[1]*B[2] - B[1]*A[2];
+    U[1] = A[2]*B[0] - B[2]*A[0];
+    U[2] = A[0]*B[1] - B[0]*A[1];
 }
 
-double dot(double a[3], double b[3]) {
-    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
+double dot(double A[3], double B[3]) {
+    // Returns dot product of Vectors A and B as scalar
+    return A[0]*B[0] + A[1]*B[1] + A[2]*B[2];
 }
 
-double** QMatrix(double omega, double inc, double OMEGA) {
-    // Perifocal to ECI matrix
-    double** matrix = malloc(3 * sizeof(double*));
-    for (int i=0; i<3; i++) matrix[i] = malloc(3 * sizeof(double));
-
-    matrix[0][0] = -sin(OMEGA) * cos(inc) * sin(omega) + cos(OMEGA) * cos(omega);
-    matrix[0][1] = cos(OMEGA) * cos(inc) * sin(omega) + sin(OMEGA) * cos(omega);
-    matrix[0][2] = sin(inc) * sin(omega);
-
-    matrix[1][0] = -sin(OMEGA) * cos(inc) * cos(omega) - cos(OMEGA) * sin(omega);
-    matrix[1][1] = cos(OMEGA) * cos(inc) * cos(omega) - sin(OMEGA) * sin(omega);
-    matrix[1][2] = sin(inc) * cos(omega);
-
-    matrix[2][0] = sin(OMEGA) * sin(inc);
-    matrix[2][1] = -cos(OMEGA) * sin(inc);
-    matrix[2][2] = cos(inc);
-    return matrix;
+void matXvec(double (*A)[3], double V[3], double U[3]) {
+    // Multiplies matrix A by vector V, storing result in vector U
+    for (int i=0; i<3; i++) U[i] = A[i][0] * V[0] + A[i][1] * V[1] + A[i][2] * V[2];
 }
 
-void vecmatmul(double** mat, double vec[3], double u[3]) {
-    // vector multiply 3d matrix to 3d vector
-    for (int i=0; i<3; i++) u[i] = mat[i][0] * vec[0] + mat[i][1] * vec[1] + mat[i][2] * vec[2];
-}
-
-void matmatmul(double** A, double** B, double** U) {
-    // vector multiply 3d matrix to 3d matrix
-
+void matXmat(double (*A)[3], double (*B)[3], double (*C)[3]) {
+    // Multiplies matrix A by matrix B, storing result in matrix C
+    for (int i=0; i < 3; i++) {
+        for (int j = 0; j < 3; j++) {
+            C[i][j] = A[i][0] * B[0][j] + A[i][1] * B[1][j] + A[i][2] * B[2][j];
+        }
+    }
 }
 
-void matT(double** A) {
-    // in place matrix transpose
+void transpose(double (*A)[3]) {
+    // in place matrix transpose of A
     double tmp;
     for (int i = 0; i < 2; i++) {
         for (int j = i+1; j < 3; j++) {
@@ -59,4 +46,83 @@ void matT(double** A) {
 }
 
 
+double (*QXp(double aop, double inc, double raan))[3] {
+    // ECI to Perifocal matrix
+    // aop: AOP
+    // raan: RAAN
+    // inc: inclination
+    // all angles in radians please
+    double (*matrix)[3] = malloc(3 * sizeof(double[3]));
+
+    matrix[0][0] = -sin(raan) * cos(inc) * sin(aop) + cos(raan) * cos(aop);
+    matrix[0][1] = cos(raan) * cos(inc) * sin(aop) + sin(raan) * cos(aop);
+    matrix[0][2] = sin(inc) * sin(aop);
+
+    matrix[1][0] = -sin(raan) * cos(inc) * cos(aop) - cos(raan) * sin(aop);
+    matrix[1][1] = cos(raan) * cos(inc) * cos(aop) - sin(raan) * sin(aop);
+    matrix[1][2] = sin(inc) * cos(aop);
+
+    matrix[2][0] = sin(raan) * sin(inc);
+    matrix[2][1] = -cos(raan) * sin(inc);
+    matrix[2][2] = cos(inc);
+    return matrix;
+}
+
+double (*QpX(double aop, double inc, double raan))[3] {
+    // Perifocal to ECI matrix    
+    // aop: AOP
+    // raan: RAAN
+    // inc: inclination
+    // all angles in radians please
+    double (*matrix)[3] = malloc(3 * sizeof(double[3]));
+
+    matrix[0][0] = -sin(raan) * cos(inc) * sin(aop) + cos(raan) * cos(aop);
+    matrix[1][0] = cos(raan) * cos(inc) * sin(aop) + sin(raan) * cos(aop);
+    matrix[2][0] = sin(inc) * sin(aop);
+
+    matrix[0][1] = -sin(raan) * cos(inc) * cos(aop) - cos(raan) * sin(aop);
+    matrix[1][1] = cos(raan) * cos(inc) * cos(aop) - sin(raan) * sin(aop);
+    matrix[2][1] = sin(inc) * cos(aop);
+
+    matrix[0][2] = sin(raan) * sin(inc);
+    matrix[1][2] = -cos(raan) * sin(inc);
+    matrix[2][2] = cos(inc);
+    return matrix;
+}
+
+double* PFE(double a, double e, double i, double aop, double m) {
+    // Returns perifocal position and velocity as 6-vector
+
+    double* perifocals = malloc(6 * sizeof(double));
+
+    double E = m + ((e*sin(m))/(1-sin(m + e) + sin(m))); // Eccentric anomaly 
+    double theta = 2 * atan(sqrt((1+e)/(1-e))*tan(E/2)); // total anomaly
+    double r = a * (1-e*e) / (1+ e*cos(theta)); // semi latus rectum
+    double v = sqrt(GM/(a*(1 - e*e))); // velocity magnitude
+
+    // PERIFOCAL R/V
+    perifocals[0] = r * cos(theta);
+    perifocals[1] = r * sin(theta);
+    perifocals[2] = 0;
+
+    perifocals[3] = v * -sin(theta);
+    perifocals[4] = v * (e + cos(theta));
+    perifocals[5] = 0;
+
+    return perifocals;
+}
+
+/*
+FORMATTERS
+*/
+void printVec(double vec[3]) {
+    printf("[%f, %f, %f]\n", vec[0], vec[1], vec[2]);
+}
+
+void printMatrix(double (*matrix)[3]) {
+    for (int i = 0; i < 3; i++) {
+        printf("[%f, %f, %f]\n", matrix[0][i], matrix[1][i], matrix[2][i]);
+    }
+}
+
 #endif