testArduino/fsfw/coordinates/Jgm3Model.h

#ifndef FRAMEWORK_COORDINATES_JGM3MODEL_H_
#define FRAMEWORK_COORDINATES_JGM3MODEL_H_

#include <stdint.h>
#include "CoordinateTransformations.h"
#include "../globalfunctions/math/VectorOperations.h"
#include "../globalfunctions/timevalOperations.h"
#include "../globalfunctions/constants.h"
#include <memory.h>


template<uint8_t DEGREE,uint8_t ORDER>
class Jgm3Model {
public:
	static const uint32_t factorialLookupTable[DEGREE+3]; //This table is used instead of factorial calculation, must be increased if order or degree is higher

	Jgm3Model() {
		y0[0] = 0;
		y0[1] = 0;
		y0[2] = 0;
		y0[3] = 0;
		y0[4] = 0;
		y0[5] = 0;

		lastExecutionTime.tv_sec = 0;
		lastExecutionTime.tv_usec = 0;
	}
	virtual ~Jgm3Model(){};

	//double acsNavOrbit(double posECF[3],double velECF[3],timeval gpsTime);

	double y0[6]; //position and velocity at beginning of RK step in EC
	timeval lastExecutionTime; //Time of last execution


	void accelDegOrd(const double pos[3],const double S[ORDER+1][DEGREE+1],const double C[ORDER+1][DEGREE+1],double* accel){
		//Get radius of this position
		double r = VectorOperations<double>::norm(pos,3);



		//Initialize the V and W matrix
		double V[DEGREE+2][ORDER+2] = {{0}};
		double W[DEGREE+2][ORDER+2] = {{0}};

		for(uint8_t m=0;m<(ORDER+2);m++){
			for(uint8_t n=m;n<(DEGREE+2);n++){
				if((n==0) && (m==0)){
					//Montenbruck "Satellite Orbits Eq.3.31"
					V[0][0] = Earth::MEAN_RADIUS / r;
					W[0][0] = 0;
				}else{
					if(n==m){
						//Montenbruck "Satellite Orbits Eq.3.29"
						V[m][m] = (2*m-1)* (pos[0]*Earth::MEAN_RADIUS/pow(r,2)*V[m-1][m-1] - pos[1]*Earth::MEAN_RADIUS/pow(r,2)*W[m-1][m-1]);
						W[m][m] = (2*m-1)* (pos[0]*Earth::MEAN_RADIUS/pow(r,2)*W[m-1][m-1] + pos[1]*Earth::MEAN_RADIUS/pow(r,2)*V[m-1][m-1]);
					}else{
						//Montenbruck "Satellite Orbits Eq.3.30"
						V[n][m] = ((2*n-1)/(double)(n-m))*pos[2]*Earth::MEAN_RADIUS / pow(r,2)*V[n-1][m];
						W[n][m] = ((2*n-1)/(double)(n-m))*pos[2]*Earth::MEAN_RADIUS / pow(r,2)*W[n-1][m];
						if(n!=(m+1)){
							V[n][m] = V[n][m] - (((n+m-1)/(double)(n-m)) * (pow(Earth::MEAN_RADIUS,2) / pow(r,2)) * V[n-2][m]);
							W[n][m] = W[n][m] - (((n+m-1)/(double)(n-m)) * (pow(Earth::MEAN_RADIUS,2) / pow(r,2)) * W[n-2][m]);
						}//End of if(n!=(m+1))
					}//End of if(n==m){
				}//End of if(n==0 and m==0)
			}//End of for(uint8_t n=0;n<(DEGREE+1);n++)
		}//End of for(uint8_t m=0;m<(ORDER+1);m++)


		//overwrite accel if not properly initialized
		accel[0] = 0;
		accel[1] = 0;
		accel[2] = 0;

		for(uint8_t m=0;m<(ORDER+1);m++){
			for(uint8_t n=m;n<(DEGREE+1);n++){
				//Use table lookup to get factorial
				double partAccel[3] = {0};
				double factor = Earth::STANDARD_GRAVITATIONAL_PARAMETER/pow(Earth::MEAN_RADIUS,2);
				if(m==0){
					//Montenbruck "Satellite Orbits Eq.3.33"
					partAccel[0] = factor * (-C[n][0]*V[n+1][1]);
					partAccel[1] = factor * (-C[n][0]*W[n+1][1]);
				}else{
					double factMN = static_cast<double>(factorialLookupTable[n-m+2]) / static_cast<double>(factorialLookupTable[n-m]);
					partAccel[0] = factor * 0.5 * ((-C[n][m]*V[n+1][m+1]-S[n][m]*W[n+1][m+1])+factMN*(C[n][m]*V[n+1][m-1]+S[n][m]*W[n+1][m-1]));
					partAccel[1] = factor * 0.5 * ((-C[n][m]*W[n+1][m+1]+S[n][m]*V[n+1][m+1])+factMN*(-C[n][m]*W[n+1][m-1]+S[n][m]*V[n+1][m-1]));
				}

				partAccel[2] = factor * ((n-m+1)*(-C[n][m]*V[n+1][m]-S[n][m]*W[n+1][m]));


				accel[0] += partAccel[0];
				accel[1] += partAccel[1];
				accel[2] += partAccel[2];
			}//End of for(uint8_t n=0;n<DEGREE;n++)
		}//End of uint8_t m=0;m<ORDER;m++
	}

	void initializeNavOrbit(const double position[3],const double velocity[3], timeval timeUTC){
		CoordinateTransformations::positionEcfToEci(position,&y0[0],&timeUTC);
		CoordinateTransformations::velocityEcfToEci(velocity,position,&y0[3],&timeUTC);
		lastExecutionTime = timeUTC;
	}


	void acsNavOrbit(timeval timeUTC, const double S[ORDER+1][DEGREE+1],const double C[ORDER+1][DEGREE+1], double outputPos[3],double outputVel[3]){

		//RK4 Integration for this timestamp
		double deltaT = timevalOperations::toDouble(timeUTC-lastExecutionTime);

		double y0dot[6] = {0,0,0,0,0,0};
		double yA[6] = {0,0,0,0,0,0};
		double yAdot[6] = {0,0,0,0,0,0};
		double yB[6] = {0,0,0,0,0,0};
		double yBdot[6] = {0,0,0,0,0,0};
		double yC[6] = {0,0,0,0,0,0};
		double yCdot[6] = {0,0,0,0,0,0};

		//Step One
		rungeKuttaStep(y0,y0dot,lastExecutionTime,S,C);

		//Step Two
		VectorOperations<double>::mulScalar(y0dot,deltaT/2,yA,6);
		VectorOperations<double>::add(y0,yA,yA,6);
		rungeKuttaStep(yA,yAdot,lastExecutionTime,S,C);

		//Step Three
		VectorOperations<double>::mulScalar(yAdot,deltaT/2,yB,6);
		VectorOperations<double>::add(y0,yB,yB,6);
		rungeKuttaStep(yB,yBdot,lastExecutionTime,S,C);

		//Step Four
		VectorOperations<double>::mulScalar(yBdot,deltaT,yC,6);
		VectorOperations<double>::add(y0,yC,yC,6);
		rungeKuttaStep(yC,yCdot,lastExecutionTime,S,C);

		//Calc new State
		VectorOperations<double>::mulScalar(yAdot,2,yAdot,6);
		VectorOperations<double>::mulScalar(yBdot,2,yBdot,6);
		VectorOperations<double>::add(y0dot,yAdot,y0dot,6);
		VectorOperations<double>::add(y0dot,yBdot,y0dot,6);
		VectorOperations<double>::add(y0dot,yCdot,y0dot,6);
		VectorOperations<double>::mulScalar(y0dot,1./6.*deltaT,y0dot,6);
		VectorOperations<double>::add(y0,y0dot,y0,6);

		CoordinateTransformations::positionEciToEcf(&y0[0],outputPos,&timeUTC);
		CoordinateTransformations::velocityEciToEcf(&y0[3],&y0[0],outputVel,&timeUTC);

		lastExecutionTime = timeUTC;


	}


	void rungeKuttaStep(const double* yIn,double* yOut,timeval time, const double S[ORDER+1][DEGREE+1],const double C[ORDER+1][DEGREE+1]){
		double rECF[3] = {0,0,0};
		double rDotECF[3] = {0,0,0};
		double accelECF[3] = {0,0,0};
		double accelECI[3] = {0,0,0};


		CoordinateTransformations::positionEciToEcf(&yIn[0],rECF,&time);
		CoordinateTransformations::velocityEciToEcf(&yIn[3],&yIn[0],rDotECF,&time);
		accelDegOrd(rECF,S,C,accelECF);
		//This is not correct, as the acceleration would have derived terms but we don't know the velocity and position at that time
		//Tests showed that a wrong velocity does make the equation worse than neglecting it
		CoordinateTransformations::positionEcfToEci(accelECF,accelECI,&time);
		memcpy(&yOut[0],&yIn[3],sizeof(yOut[0])*3);
		memcpy(&yOut[3],accelECI,sizeof(yOut[0])*3);
	}
};






#endif /* FRAMEWORK_COORDINATES_JGM3MODEL_H_ */