Merge pull request #1297 from indy91/GPMFix

GPM height and circularization burns converge better on highly elliptical orbits; RTCC routine that finds altitude crossing as well

Orbitersdk/samples/ProjectApollo/src_launch/rtcc.cpp

@@ -6245,6 +6245,9 @@ bool RTCC::GeneralManeuverProcessor(GMPOpt *opt, VECTOR3 &dV_i, double &P30TIG)
 	case 8: //PMMAPD error for resultant apo/peri
 		PMXSPT("PMMGPM", 136);
 		break;
+	case 9: //Resultant orbit non-elliptical
+		PMXSPT("PMMGPM", 137);
+		break;
 	}
 
 	P30TIG = PZGPMDIS.GET_TIG;

Orbitersdk/samples/ProjectApollo/src_rtccmfd/GeneralPurposeManeuver.cpp

@@ -717,6 +717,9 @@ void RTCCGeneralPurposeManeuverProcessor::FlightControllerInput()
 void RTCCGeneralPurposeManeuverProcessor::HeightManeuver(bool circ)
 {
 	double DH, dv_h, V_a2, x_a, E_a, V_a, r_b_dot, V_H_a, r_D;
+	int IMAX;
+
+	IMAX = 10;
 
 	if (circ)
 	{
@@ -775,17 +778,37 @@ void RTCCGeneralPurposeManeuverProcessor::HeightManeuver(bool circ)
 		{
 			DH = r_D - sv_temp.R;
 		}
-		if (abs(DH) < eps1 || I > 5)
+		if (abs(DH) < eps1 || I > IMAX)
 		{
 			break;
 		}
 		dv_h = mu * DH / (4.0*pow(sv_a.coe_osc.a, 2)*V_a);
+		//Limit change in DV. Is equation above even correct?
+		if (abs(dv_h) > 1000.0)
+		{
+			if (dv_h > 0.0)
+			{
+				dv_h = 1000.0;
+			}
+			else
+			{
+				dv_h = -1000.0;
+			}
+		}
+
 		V_H_a = V_H_a + dv_h;
 		V_a2 = V_H_a * V_H_a + r_b_dot * r_b_dot;
 		V_a = sqrt(V_a2);
 		sv_a.coe_osc.a = mu * sv_b.R / (2.0*mu - sv_b.R*V_a2);
 		x_a = sv_a.R*sv_a.R*V_H_a*V_H_a;
 		sv_a.coe_osc.e = sqrt(1.0 - x_a / mu / sv_a.coe_osc.a);
+		//Orbit too eccentric
+		if (sv_a.coe_osc.e >= 1.0)
+		{
+			ErrorIndicator = 9;
+			return;
+		}
+
 		E_a = pRTCC->GLQATN(sv_b.R*r_b_dot*sqrt(sv_a.coe_osc.a), sqrt(mu)*(sv_a.coe_osc.a - sv_b.R));
 		sv_a.coe_osc.l = E_a - sv_a.coe_osc.e*sin(E_a);
 		sv_a.f = pRTCC->GLQATN(sv_b.R*r_b_dot*sqrt(x_a), x_a - mu * sv_b.R);
@@ -798,7 +821,7 @@ void RTCCGeneralPurposeManeuverProcessor::HeightManeuver(bool circ)
 		{
 			sv_a.coe_osc.g += PI2;
 		}
-	} while (I < 5);
+	} while (I < IMAX);
 }
 
 int RTCCGeneralPurposeManeuverProcessor::OptimumApsidesChange()

Orbitersdk/samples/ProjectApollo/src_rtccmfd/GeneralPurposeManeuver.h

@@ -198,7 +198,8 @@ class RTCCGeneralPurposeManeuverProcessor : public RTCCModule
 	VECTOR3 DV_Vector;
 	double Pitch_Man, Yaw_Man, Node_Ang, Del_G, DV, H_man, lat_man, lng_man;
 	//0 = good maneuver, 1 = Maneuver cannot be performed in this orbit, 2 = Maneuver cannot be performed at this point in the orbit, 3 = AEG failed to converge, 4 = Unrecoverable AEG error
-	//5 = Unable to obtain libration matrix, 6 = Unable to advance to selenographic argument of latitude, 7 = PMMAPD error for current apo/peri, 8 = PMMAPD error for resultant apo/peri
+	//5 = Unable to obtain libration matrix, 6 = Unable to advance to selenographic argument of latitude, 7 = PMMAPD error for current apo/peri, 8 = PMMAPD error for resultant apo/peri, 
+	//9 = Resultant orbit non-elliptical
 	int ErrorIndicator;
 
 	//CONSTANTS

Orbitersdk/samples/ProjectApollo/src_rtccmfd/rtcc_intermediate_library_programs.cpp

@@ -758,10 +758,11 @@ int RTCC::PITCIR(AEGHeader header, AEGDataBlock in, double R_CIR, AEGDataBlock &
 {
 	//Output: 0 = no error, 1 = essentially circular orbit, 2 = unrecoverable AEG error, 3 = requested height not in orbit, 4 = failed to converge on radius
 
-	double cos_f_CI, f_CI, dt, sgn, ddt, eps_t;
-	int I;
+	double cos_f_CI, f_CI, dt, sgn, ddt, eps_t, l2, dl;
+	int I, IMAX;
 	bool fail;
 
+	IMAX = 10;
 	eps_t = 0.01;
 
 	if (in.ENTRY == 0)
@@ -849,26 +850,28 @@ int RTCC::PITCIR(AEGHeader header, AEGDataBlock in, double R_CIR, AEGDataBlock &
 			f_CI += PI2;
 		}
 		//Calculate angle difference
-		ddt = f_CI - out.f;
-		if (ddt > PI)
+		l2 = OrbMech::TrueToMeanAnomaly(f_CI, out.coe_osc.e);
+		dl = l2 - out.coe_osc.l;
+		if (dl > PI)
 		{
-			ddt -= PI2;
+			dl -= PI2;
 		}
-		else if (ddt < -PI)
+		else if (dl < -PI)
 		{
-			ddt += PI2;
+			dl += PI2;
 		}
+
 		//Calculate ddt
-		ddt = ddt / (out.g_dot + out.l_dot);
+		ddt = dl / (out.g_dot + out.l_dot);
 		if (abs(ddt) > eps_t)
 		{
 			dt = dt + ddt;
 		}
 		I++;
-	} while (abs(ddt) > eps_t && I < 10);
+	} while (abs(ddt) > eps_t && I < IMAX);
 
 	//Failed to converge
-	if (I == 10)
+	if (I >= IMAX)
 	{
 		PMXSPT("PITCIR", 19);
 		return 4;