Merge pull request #1264 from indy91/EnckeIntegratorFixes

RTCC Encke free-flight integrator corrections

Orbitersdk/samples/ProjectApollo/src_launch/rtcc.cpp

@@ -21088,6 +21088,10 @@ void RTCC::PMMDAN(AEGHeader Header, AEGDataBlock aeg, int IND, int &ERR, double
 
 void RTCC::EMDCHECK(int veh, int opt, double param, double THTime, int ref, bool feet)
 {
+	//INPUTS:
+	//opt: 1 = GMT, 2 = GET, 3 = MVI, 4 = MVE, 5 = RAD, 6 = ALT, 7 = FPA, 8 = LAT, 9 = LNG
+	//ref: 0 = ECI, 1 = ECT, 2 = MCI, 3 = MCT
+
 	EphemerisData sv_out, sv_conv, sv_inert;
 	MissionPlanTable *table;
 	OrbitEphemerisTable *ephem;
@@ -21160,10 +21164,16 @@ void RTCC::EMDCHECK(int veh, int opt, double param, double THTime, int ref, bool
 	case 7:
 		sprintf(EZCHECKDIS.Option, "FPA");
 		break;
+	case 8:
+		sprintf(EZCHECKDIS.Option, "LAT");
+		break;
+	case 9:
+		sprintf(EZCHECKDIS.Option, "LNG");
+		break;
 	}
 
 	//All options except MVI and MVE
-	if (opt == 1 || opt == 2 || opt == 5 || opt == 6 || opt == 7)
+	if (opt == 1 || opt == 2 || (opt >= 5 && opt <= 9))
 	{
 		EphemerisData svtemp;
 		double THGMT;
@@ -21342,6 +21352,14 @@ void RTCC::EMDCHECK(int veh, int opt, double param, double THTime, int ref, bool
 			{
 				intab.CutoffIndicator = 4;
 			}
+			else if (opt == 8)
+			{
+				intab.CutoffIndicator = 8;
+			}
+			else if (opt == 9)
+			{
+				intab.CutoffIndicator = 7;
+			}
 			else
 			{
 				intab.CutoffIndicator = 1;
@@ -21564,11 +21582,13 @@ void RTCC::EMDCHECK(int veh, int opt, double param, double THTime, int ref, bool
 
 	if (ref == 0)
 	{
+		//ECI: lambda is mean inertial, lambda_D = is mean geographic
 		EZCHECKDIS.lambda_D -= SystemParameters.MCLAMD + sv_out.GMT*OrbMech::w_Earth;
 	}
 	else if (ref == 1)
 	{
-		EZCHECKDIS.lambda_D -= sv_out.GMT*OrbMech::w_Earth;
+		//ECT: lambda is true geographic, lambda_D is true inertial
+		EZCHECKDIS.lambda -= sv_out.GMT*OrbMech::w_Earth;
 	}
 	OrbMech::normalizeAngle(EZCHECKDIS.lambda);
 	OrbMech::normalizeAngle(EZCHECKDIS.lambda_D);
@@ -32034,6 +32054,14 @@ int RTCC::EMGTVMED(std::string med, std::vector<std::string> data)
 		{
 			opt = 7;
 		}
+		else if (data[1] == "LAT")
+		{
+			opt = 8;
+		}
+		else if (data[1] == "LNG")
+		{
+			opt = 9;
+		}
 		else
 		{
 			return 2;
@@ -32062,6 +32090,10 @@ int RTCC::EMGTVMED(std::string med, std::vector<std::string> data)
 			{
 				param = sin(param*RAD);
 			}
+			else if (opt == 8 || opt == 9)
+			{
+				param = param*RAD;
+			}
 		}
 		bool hasTHT;
 		double THTime = -1.0;
@@ -32080,7 +32112,7 @@ int RTCC::EMGTVMED(std::string med, std::vector<std::string> data)
 		}
 		if (opt >= 5)
 		{
-			//Mandatory for RAD, ALT, FPA
+			//Mandatory for RAD, ALT, FPA, LAT, LNG
 			if (hasTHT == false)
 			{
 				return 2;

Orbitersdk/samples/ProjectApollo/src_rtccmfd/ApolloRTCCMFD.cpp

@@ -4931,7 +4931,7 @@ void ApolloRTCCMFD::menuMPTInitM50M55Vehicle()
 void ApolloRTCCMFD::CheckoutMonitorCalc()
 {
 	bool CheckoutMonitorCalcInput(void* id, char *str, void *data);
-	oapiOpenInputBox("Format: U02, CSM or LEM, Indicator (GMT,GET,MVI,MVE,RAD,ALT,FPA), Parameter, Threshold Time (opt.), Reference (ECI,ECT,MCI,MCT) (opt.), FT (opt.);", CheckoutMonitorCalcInput, 0, 50, (void*)this);
+	oapiOpenInputBox("Format: U02, CSM or LEM, Indicator (GMT,GET,MVI,MVE,RAD,ALT,FPA,LAT,LNG), Parameter, Threshold Time (opt.), Reference (ECI,ECT,MCI,MCT) (opt.), FT (opt.);", CheckoutMonitorCalcInput, 0, 50, (void*)this);
 }
 
 bool CheckoutMonitorCalcInput(void *id, char *str, void *data)

Orbitersdk/samples/ProjectApollo/src_rtccmfd/CoastNumericalIntegrator.cpp

@@ -438,7 +438,7 @@ VECTOR3 CoastIntegrator2::adfunc(VECTOR3 R)
 			//Get Earth rotation matrix only during initialization. For the Moon the libration matrix is updated by the PLEFEM call below
 			if (P == BODY_EARTH)
 			{
-				pRTCC->ELVCNV(CurrentTime(), RTCC_COORDINATES_ECT, RTCC_COORDINATES_ECI, Rot);
+				pRTCC->ELVCNV(CurrentTime(), RTCC_COORDINATES_ECI, RTCC_COORDINATES_ECT, Rot);
 			}
 		}
 

Orbitersdk/samples/ProjectApollo/src_rtccmfd/EnckeIntegrator.cpp

@@ -130,7 +130,7 @@ void EnckeFreeFlightIntegrator::Propagate(EMMENIInputTable &in)
 		//1 meter tolerance for radius and height
 		DEV = 1.0;
 	}
-	else if (ISTOPS == 4 || ISTOPS == 6)
+	else if (ISTOPS == 4 || ISTOPS == 6 || ISTOPS == 7 || ISTOPS == 8)
 	{
 		//0.0001� tolerance
 		DEV = 0.0001*RAD;
@@ -221,7 +221,15 @@ void EnckeFreeFlightIntegrator::Edit()
 	}
 	else
 	{
-		dt_max = min(dt_lim, K*OrbMech::power(rr, 1.5) / sqrt(mu));
+		if (DRAG > 0.0 && P == BODY_EARTH && rr < 6499920.0)
+		{
+			//Special atmospheric integration step logic
+			dt_max = 10.0;
+		}
+		else
+		{
+			dt_max = min(dt_lim, K*OrbMech::power(rr, 1.5) / sqrt(mu));
+		}
 	}
 
 	//Should we even check?
@@ -237,6 +245,7 @@ void EnckeFreeFlightIntegrator::Edit()
 	{
 		if (ISTOPS == 2 || ISTOPS == 5)
 		{
+			//Radius
 			FUNCT = length(R);
 			if (P == BODY_EARTH)
 			{
@@ -249,6 +258,7 @@ void EnckeFreeFlightIntegrator::Edit()
 		}
 		else if (ISTOPS == 3)
 		{
+			//Altitude
 			if (P == BODY_EARTH)
 			{
 				FUNCT = length(R) - OrbMech::R_Earth;
@@ -262,6 +272,7 @@ void EnckeFreeFlightIntegrator::Edit()
 		}
 		else if (ISTOPS == 4)
 		{
+			//Sine of flight path angle
 			FUNCT = dotp(unit(R), unit(V));
 			if (P == BODY_EARTH)
 			{
@@ -274,6 +285,7 @@ void EnckeFreeFlightIntegrator::Edit()
 		}
 		else if (ISTOPS == 6)
 		{
+			//Ascending node
 			VECTOR3 NVEC, H_ECT;
 			EphemerisData2 sv_temp, sv_ECT;
 
@@ -287,6 +299,47 @@ void EnckeFreeFlightIntegrator::Edit()
 			NVEC = unit(crossp(_V(0, 0, 1), H_ECT));
 			RCALC = OrbMech::PHSANG(sv_ECT.R, sv_ECT.V, NVEC);
 		}
+		else if (ISTOPS == 7)
+		{
+			//Longitude
+			double lat;
+			OrbMech::latlong_from_r(R_EF, lat, FUNCT);
+
+			if (P == BODY_EARTH)
+			{
+				RCALC = FUNCT - STOPVAE - OrbMech::w_Earth*CurrentTime();
+			}
+			else
+			{
+				RCALC = FUNCT - STOPVAM;
+			}
+
+			//modulo between -PI and PI
+			if (RCALC > 0)
+			{
+				RCALC = fmod(RCALC + PI, PI2) - PI;
+			}
+			else
+			{
+				RCALC = fmod(RCALC - PI, PI2) + PI;
+			}
+		}
+		else if (ISTOPS == 8)
+		{
+			//Latitude
+
+			double lng;
+			OrbMech::latlong_from_r(R_EF, FUNCT, lng);
+
+			if (P == BODY_EARTH)
+			{
+				RCALC = FUNCT - STOPVAE;
+			}
+			else
+			{
+				RCALC = FUNCT - STOPVAM;
+			}
+		}
 	}
 
 	IEND = ISTOPS;
@@ -320,8 +373,20 @@ void EnckeFreeFlightIntegrator::Edit()
 	else if (INITE == -1)
 	{
 		//Not bounded
-		if (RCALC*RES1 >= 0)
+		bool found = false;
+
+		if (RCALC*RES1 < 0.0)
+		{
+			//For longitude search, prevent ambiguity
+			if (!(ISTOPS == 7 && abs(RCALC) > PI05))
+			{
+				found = true;
+			}
+		}
+
+		if (!found)
 		{
+			//No sign switch, not bounded
 			goto EMMENI_Edit_4A;
 		}
 
@@ -502,7 +567,7 @@ void EnckeFreeFlightIntegrator::adfunc()
 			//Get Earth rotation matrix only during initialization. For the Moon the libration matrix is updated by the PLEFEM call below
 			if (P == BODY_EARTH)
 			{
-				pRTCC->ELVCNV(CurrentTime(), RTCC_COORDINATES_ECT, RTCC_COORDINATES_ECI, Rot);
+				pRTCC->ELVCNV(CurrentTime(), RTCC_COORDINATES_ECI, RTCC_COORDINATES_ECT, Rot);
 				U_Z = tmul(Rot, _V(0, 0, 1));
 			}
 		}

Orbitersdk/samples/ProjectApollo/src_rtccmfd/EnckeIntegrator.h

@@ -144,7 +144,7 @@ class EnckeFreeFlightIntegrator : public RTCCModule
 	double VAR;
 	//Maximum time to integrate
 	double TMAX;
-	//Stop condition (1 = time, 2 = radial distance, 3 = altitude above Earth or moon, 4 = flight-path angle, 5 = first reference switch, 6 = first ascending node relative to the Earth)
+	//Stop condition (1 = time, 2 = radial distance, 3 = altitude above Earth or moon, 4 = flight-path angle, 5 = first reference switch, 6 = first ascending node relative to the Earth, 7 = longitude, 8 = latitude)
 	int ISTOPS;
 	//Reference frame of desired stopping parameter (0 = Earth, 1 = Moon, 2 = both)
 	int StopParamRefFrame;

Orbitersdk/samples/ProjectApollo/src_rtccmfd/RTCCTables.h

@@ -376,7 +376,7 @@ struct EMSMISSInputTable
 	int ManCutoffIndicator;
 	//Descent burn indicator
 	bool DescentBurnIndicator = false;
-	//Cut-off indicator (1 = Time, 2 = radial distance, 3 = altitude above Earth or moon, 4 = flight-path angle, 5 = first reference switch)
+	//Cut-off indicator (1 = Time, 2 = radial distance, 3 = altitude above Earth or moon, 4 = flight-path angle, 5 = first reference switch, 6 = first ascending node relative to the Earth, 7 = longitude, 8 = latitude)
 	int CutoffIndicator = 1;
 	//Integration direction indicator (+X-forward, -X-backward)
 	double IsForwardIntegration = 1.0;