Math help... calculation of transfer orbits between non coplanar, non circular orbits

[zillion42]

Hi,

I wrote a small simulation of our solar system using the irrlicht engine. It works very well by now. The planets are positioned by:

1. Finding the julian day since the date of the elements (J2000)
2. Finding the average (or mean) orbital elements of the planet
3. Finding the true anomaly, the angle to the planet from perihelion.
4. Finding the heliocentric radius, distance from planet to sun
5. Finding the heliocentric ecliptic rectangular coordinates of the planet

additionally I have one satellite, being the camera which can be steered around using keyboard inputs. Gravitational pulls influence the satellite of course and a number of realistic results can be achieved, like geo syncronous orbits, hohmann transfers between different altitudes etc. The numbers that run my gravity sim are indeed in a very realistic range, and the approach could be described as a patched conic trajectory simulation. One sim frame equals 100 seconds.
I recently started adding a projection mode or call it planning mode where the user can utilize at any point and in any direction a nifty "burn" button which de- or accelerates the satellite to a certain speed in km/s and projects onward to plot the trajectory.

Green line from earth to mars. I also added a function with the name calcLaunchSolution which attempts to find valid launch windows between planets and prints a valid delta-V or lets just call it velocity to enter a transfer orbit. It also works very well and I just keep on missing mars by a few hundred thousend km. Unfortunately.
I guess most of you here already know why that is so, but to state it again, it stems from the fact that all the math stuff one is able to find for hohmann transfers assume perfect circular and perfectly coplanar orbits, which none of our planets actually has. Well you could argue that since the eccliptic is considered to be on the plane of the coordinate system we use (heliocentric eccliptic) earth does, but no there is inclination and of course eccentricity in all 8 planets orbits.

So, additionally to finding the right day, I need to find for the departure planet the distance from the sun at time of departure and the distance of the arrival planet from sun at the time of arrival. Those two distances are then the perigee and apogee of my transfer. But because in finding the right day I use the mean distances, the transfer period is schoolbook style 258 days(mars), in reality it varies +/- 3/4 days everything is then corrupt.
Secondly even if I get close by trial and error the inclination change between the orbits makes me arrive on a completly different plane, which corrupts even more.
Well, long talk, I hope someone knows what I'm talking about and can provide me with interesting links and readings on the subject so I can get that fixed asap.

Many Thanks in advance

EDIT:
Long day, I found a promising fortran app that I could translate with f2c using a lambert solver, unfortunately the provided source only extends a library which costs 179 $, so I guess I need to work it out using the Java Astrodynamics Toolkit lambert solver, but indeed looks promising...

...and for the record, these are my attempts so far, be warned it is a mess and reflects many sleepless nights...

void calcLaunchSolution(u32 start, u32 destination)// f64 *launchday, core::vector3df *velocity
{
	// 1st.
	//calculate mean transfer orbit period
	f64 distPerigee = pMeanDist2[start]*1000.0f;
	f64 distApogee = pMeanDist2[destination]*1000.0f;

	f64 semiMajorAxisM = ((distPerigee+distApogee)/2);
	f64 pConst = Gconst*pMass[9];

	//calculate orbital period in days, for an ellipse that is tangential to start and destination planet
	//the formula for an ellipse is
	//p = 2PI*sqrt(semi-major axis^3/GM) where p is in seconds and semi-major axis in meters
	//then divide it by 2 to give us a transfer orbit period
	f64 TOPeriod = ((((2*core::PI*sqrt(pow(semiMajorAxisM,3)/pConst))/60.0f)/60.0f)/24.0f);
	//printf("Period in days: %f\n",oPeriod);
	//Half of the elipse is the Transfer Orbit, therefore we divide by 2;
	TOPeriod /= 2.0f;
	printf("\n\nTransfer Period = %.9f\n",TOPeriod);
	
	//2nd
	//Find time delta angle for the launch window
	//calculate destination Planet orbital Period DIVIDED 2 MINUS transfer orbit Period
	//if the destination planets orbital period is greater than the start planets orbital period,
	//ie. it moves more slow, then the launch window occurs after the day with matching delta angle,
	//otherwise before... 
	
	f64 angleOffset = 180.0f + (core::PI * (1.0 - TOPeriod / orbitT[destination])*core::RADTODEG64);
	
	if(angleOffset > 360)
		angleOffset - 360.0;
	
	if(orbitT[start] < orbitT[destination])
		angleOffset = angleOffset;
	else
		angleOffset = 360.0 - angleOffset;
	
	printf("Angle Offset = %.9f\n",angleOffset);
	
	if(orbitT[start] < orbitT[destination])
	{
		angleOffset = angleOffset;
	}
	else
	{
		angleOffset = 360.0f - angleOffset;
	}

	//3rd 
	//find a day with this angular offset between start and destination planet
	bool windowNotFound = true;

	f64 testday = Day;
	f64 approxWindowDay;
	core::stringc date;
	
	u32 dates = 0;
	
	while(dates < 1)
	{
		testday += 0.1f;
		date = mySolarSystem->J2000ToGreg(testday).c_str();
		core::array<core::vector3df> startPos = mySolarSystem->calculateAtDay(testday, start);
		core::array<core::vector3df> DestPos = mySolarSystem->calculateAtDay(testday, destination);
		core::vector3df startAngle = startPos[0].getHorizontalAngle();
		core::vector3df DestAngle = DestPos[0].getHorizontalAngle();
		
		f32 delta = DestAngle.Y - startAngle.Y;

		//printf("Delta:%f,- %s\n",delta,date);

		if( core::abs_(delta - angleOffset) <= 0.1f)
		{
			approxWindowDay = testday;
			date = mySolarSystem->J2000ToGreg(approxWindowDay).c_str();
			printf("Approximate Launch Date = %s\n",date);
			windowNotFound = false;
			dates++;
		}
		
	};
	//Day = approxWindowDay;

	//4th
	//calculate precise transfer period based on determined launch date
	core::array<core::vector3df> perigeePos = mySolarSystem->calculateAtDay(approxWindowDay,start);
	core::array<core::vector3df> apogeePos = mySolarSystem->calculateAtDay(approxWindowDay + TOPeriod,destination);
	f32 lDPerigee = perigeePos[0].getLength()*1000.0f;
	f32 lDApogee = apogeePos[0].getLength()*1000.0f;
	f64 semiMajorAxisLD = ((lDPerigee+lDApogee)/2);
	f64 lDTOPeriod = ((((2*core::PI*sqrt(pow(semiMajorAxisLD,3)/pConst))/60.0f)/60.0f)/24.0f)/2.0f;
	printf("precise Transfer Period = %f\n",lDTOPeriod);
	
	//4th
	//Calculate precise launch date
	angleOffset = 180.0f + (core::PI * (1.0 - lDTOPeriod / orbitT[destination])*core::RADTODEG64);
	if(angleOffset > 360)
		angleOffset - 360;
	printf("True Angle Offset = %.9f\n",angleOffset);
	
	if(orbitT[start] < orbitT[destination])
	{
		angleOffset = angleOffset;
	}
	else
	{
		angleOffset = 360.0f - angleOffset;
	}
	
	testday = Day;
	dates = 0;
	
	while(dates < 1)
	{
		testday += 0.1f;
		date = mySolarSystem->J2000ToGreg(testday).c_str();
		core::array<core::vector3df> startPos = mySolarSystem->calculateAtDay(testday, start);
		core::array<core::vector3df> DestPos = mySolarSystem->calculateAtDay(testday, destination);
		core::vector3df startAngle = startPos[0].getHorizontalAngle();
		core::vector3df DestAngle = DestPos[0].getHorizontalAngle();
		f32 delta = DestAngle.Y - startAngle.Y;
		
		//printf("Delta:%f,- %s\n",delta,date);

		if( core::abs_(delta - angleOffset) <= 0.1f)
		{
			approxWindowDay = testday;
			date = mySolarSystem->J2000ToGreg(approxWindowDay).c_str();
			printf("Precise Launch Date = %s\n",date);
			windowNotFound = false;
			dates++;
		}
		
	};
	
	Day = approxWindowDay;
	f64 arrivalDay = Day+lDTOPeriod;
	core::stringc arrivalDate = mySolarSystem->J2000ToGreg(arrivalDay).c_str();
	printf("Precise Arrival Date = %s\n",arrivalDate);
	
	//5th
	//Calculate the Orbital Velocity at perigee
	//Velocity at Perigee  	=      sqrt( GM * ( ( 2 * distApogee) / ( distPerigee * ( distPerigee+distApogee ) ) ) )
	core::array<core::vector3df> nperigeePos = mySolarSystem->calculateAtDay(approxWindowDay,start);
	core::array<core::vector3df> napogeePos = mySolarSystem->calculateAtDay(approxWindowDay + lDTOPeriod,destination);
	f64 nlDPerigee = nperigeePos[0].getLength()*1000.0f;
	f64 nlDApogee = napogeePos[0].getLength()*1000.0f;
	f64 ellipticOrbVelocity = sqrt( pConst * ( ( 2 * nlDApogee) / ( nlDPerigee * ( nlDPerigee + nlDApogee ) ) ) );
	ellipticOrbVelocity /= 1000.0;
	printf("Velocity = %.9f\n",ellipticOrbVelocity);
	//Day = approxWindowDay;
	ellipticOrbVelocity /= scaleF;
	printf("Dist Perigee %.9f\n",lDPerigee/1000.0);
	printf("Dist  Apogee %.9f\n",lDApogee/1000.0);

	f64 yDiff = core::max_(perigeePos[0].Y,apogeePos[0].Y) - core::min_(perigeePos[0].Y,apogeePos[0].Y);
	f64 yDiff2d = (yDiff/scaleF) / (lDTOPeriod/fpsincrement);
	yDiffSpeed = core::vector3df(0,yDiff2d,0);
}

 

[zillion42]

:woohoo:

Several days later, and a few grey hair, I have a working lambert solver, ported from the Java Astrodynamical Toolkit... unfortunately precision still wont allow for mars to catch my satellite. I am convinced however, that I have done nearly everything right. I have a few issues with inverted vectors multiplying by -1.0 helps...
For the sake of sharing, here it is:

replace f64 for double, u32 for unsigned int, s32 for signed int and core::vector3df for whatever vector typedef you use, and it becomes standard c++
well...

someVector.getLength() = magnitude of vector;
someVector.dotProduct() = dotProduct of vector;
.crossProduct, normalize = self-explanatory

someVector - someVector2 = someNewVector(x-x²),(y-y²),(z-z²)
same with +
* and / with skalars and so forth...

enjoy

#include <irrlicht.h>

using namespace irr;

f64 s = 0.0;
f64 c = 0.0;
f64 tof = 0.0;
f64 mu = 0.0;
bool aflag = false;
bool bflag = false;
bool debug_print = true;

f64 getalpha(f64 a)
{
	f64 alpha = 2.0 * asin(sqrt(s / (2.0 * a)));
	if (aflag)
	{
		alpha = 2.0 * core::PI64 - alpha;
	}
	return alpha;
}

f64 getbeta(f64 a)
{
	f64 beta = 2.0 * asin(sqrt((s - c) / (2.0 * a)));
	if (bflag)
	{
		beta = -1.0 * beta;
	}
	return beta;
}

f64 getdt(f64 a, f64 alpha, f64 beta)
{
	f64 sa = sin(alpha);
	f64 sb = sin(beta);
	f64 dt = pow(a, 1.5) * (alpha - sa - beta + sb) / sqrt(mu);
	return dt;
}

/** Evaluate the delta-t function f
* @param a semi-major axis.
* @return
*/
f64 evaluate(f64 a)
{
	f64 alpha = getalpha(a);
	f64 beta = getbeta(a);
	f64 out = tof - getdt(a,alpha,beta);
	return out;
}

/** Find the solution using RegulaFalsi.
 * @param x1 lower limit on x.
 * @param x2 upper limit on x.
 * @param dxmin
 * @return
 */
f64 regulaFalsi(f64 x1, f64 x2) {
	f64 xlow;
	f64 xhigh;
	f64 del;
	f64 out = 0.0;
	f64 f;
	
	u32 err = 0;
	f64 dxmin = 1.0e-15;
	f64 accuracy = 1.0e-6;

	/*NO IDEA !

	f64 fl = this.func.evaluate(x1);
	f64 fh = this.func.evaluate(x2);

	//from JAT ZeroFinder:
	
	* The ZeroFinder Class provides a way to solve a scalar f(x) = 0.
	* These functions have been translated from Numerical Recipes.
	* Currently there are: Regula Falsi, Secant and Ridder's methods.
	* The function f is passed via the ScalarFunction interface.

	//from JAT Scalarfunction:

	* The ScalarFunction interface provides the mechanism for passing a method
	* that evaluates a function to a solver.
	*/
	
	//HOPE THIS WORKS ?!?
	f64 fl = evaluate(x1);
	f64 fh = evaluate(x2);

	f64 test = fl * fh;

	if (test > 0.0) {
		if (fl == 0.0) return x1;
		if (fh == 0.0) return x2;
		err++;
		printf("Root must be bracketed in ZeroFinder\n");
	}

	if (fl < 0.0) 
	{
		xlow = x1;
		xhigh = x2;
	} 
	else 
	{
		xlow = x2;
		xhigh = x1;
		f64 temp = fl;
		fl = fh;
		fh = temp;
	}

	f64 dx = xhigh - xlow;

	for(int i = 1; i < 10000; i++) 
	{
		out = xlow + dx * fl / (fl - fh);
		//Same as above, hope it works...
		//f = this.func.evaluate(out);
		f = evaluate(out);

		if (f < 0.0) {
			del = xlow - out;
			xlow = out;
			fl = f;
		} else {
			del = xhigh - out;
			xhigh = out;
			fh = f;
		}

		dx = xhigh - xlow;

		if ((abs(del) < dxmin) || (abs(f) < accuracy)) {
			return out;
		}
	}

	printf("Regula Falsi exceeded 10000 iterations\n");
	return out;
}

/** Computes the delta-v's required to go from r0,v0 to rf,vf.
* @return Total delta-v (magnitude) required.
* @param dt Time of flight
* @param r0 Initial position vector.
* @param v0 Initial velocity vector.
* @param rf Desired final position vector.
* @param vf Desired final velocity vector.
* //pointers to return results
* @param deltaV0 computed Initial delta-V.
* @param deltaV1 computed Final delta-V.
* @param totalV0 computed Initial total-V.
* @param totalV0 computed Final total-V.
*/
void LambertCompute(f64 GM, 
					core::vector3df r0,
					core::vector3df v0,
					core::vector3df rf,
					core::vector3df vf, 
					f64 dt,
					core::vector3df* deltaV0,
					core::vector3df* deltaV1,
					core::vector3df* totalV0,
					core::vector3df* totalV1)
{
	s = 0.0;
	c = 0.0;
	aflag = false;
	bflag = false;
	debug_print=true;

	f64 tp = 0.0;
	f64 magr0 = r0.getLength();
	f64 magrf = rf.getLength();
	//GM is expected in km^3/s^2
	mu = GM / 1000000000.0;
	//time of flight expected in seconds
	dt *= 86400;
	tof = dt;

	core::vector3df dr = r0 - (rf);
	c = dr.getLength();
	s = (magr0 + magrf + c) / 2.0;
	f64 amin = s / 2.0;
	
	if(debug_print)
		printf("amin = %.9E\n", amin);
	
	f64 dtheta = acos(r0.dotProduct(rf) / (magr0 * magrf));

	//dtheta = 2.0 * core::PI64 - dtheta;

	if(debug_print)
		printf("dtheta = %.9E\n", dtheta);
	
	if (dtheta < core::PI64)
	{
		tp = sqrt(2.0 / (mu)) * (pow(s, 1.5) - pow(s - c, 1.5)) / 3.0;
	}
	if (dtheta > core::PI64)
	{
		tp = sqrt(2.0 / (mu)) * (pow(s, 1.5) + pow(s - c, 1.5)) / 3.0;
		bflag = true;
	}

	if(debug_print)
		printf("tp = %.9f\n", tp);

	f64 betam = getbeta(amin);
	f64 tm = getdt(amin, core::PI64, betam);

	if(debug_print)
		printf("tm = %.9E\n", tm);

	if (dtheta == core::PI64)
	{
		printf(" dtheta = 180.0. Do a Hohmann\n");
		return;
	}

	f64 ahigh = 1000.0 * amin;
	f64 npts = 3000.0;
	
	if(debug_print)
		printf("dt = %.9E seconds\n", dt);

	if(debug_print)
		printf("************************************************\n");

	if (dt < tp)
	{
		printf(" No elliptical path possible \n");
		return;
	}

	if (dt > tm)
	{
		aflag = true;
	}

	f64 fm = evaluate(amin);
	f64 ftemp = evaluate(ahigh);

	if ((fm * ftemp) >= 0.0)
	{
		printf(" initial guesses do not bound \n");
		return;
	}

	//ZeroFinder regfalsi = new ZeroFinder(this, 10000, 1.0E-6, 1.0E-15);

	f64 sma = regulaFalsi(amin, ahigh);

	f64 alpha = getalpha(sma);
	f64 beta = getbeta(sma);
	
	f64 de = alpha - beta;
	
	f64 f = 1.0 - (sma / magr0) * (1.0 - cos(de));
	f64 g = dt - sqrt(sma * sma * sma / mu) * (de - sin(de));
	
	core::vector3df newv0;
	core::vector3df newvf;

	newv0.X = (rf.X - f * r0.X) / g;
	newv0.Y = (rf.Y - f * r0.Y) / g;
	newv0.Z = (rf.Z - f * r0.Z) / g;
	
	//if it wont work for you, take away the -1.0 multiplications
	//I guess my coordinate system is little screwed...
	*deltaV0 = (newv0 - (v0*-1.0))*-1.0;
	*totalV0 = (newv0*-1.0);

	if(debug_print)
		printf("deltav-0 X=%.9f, Y=%.9f, Z=%.9f\n",deltaV0->X,deltaV0->Y,deltaV0->Z);

	f64 dv0 = deltaV0->getLength();

	f64 fdot = -1.0 * (sqrt(mu * sma) / (magr0 * magrf)) * sin(de);
	f64 gdot = 1.0 - (sma / magrf) * (1.0 - cos(de));

	newvf.X = fdot * r0.X + gdot * newv0.X;
	newvf.Y = fdot * r0.Y + gdot * newv0.Y;
	newvf.Z = fdot * r0.Z + gdot * newv0.Z;
	
	//Same here:
	//if it wont work for you, take away the -1.0 multiplications
	//I guess my coordinate system is little screwed...
	*deltaV1 = ((vf*-1.0) - newvf)*-1.0;
	*totalV1 = (newvf*-1.0);

	if(debug_print)
		printf("deltav-f X=%.9f, Y=%.9f, Z=%.9f\n",deltaV1->X,deltaV1->Y,deltaV1->Z);

	f64 dvf = deltaV1->getLength();
	f64 totaldv = dv0 + dvf;

	if(debug_print)
		printf("\n\nInitial DeltaV dv0 = %.9f\nFinal DeltaV   dvf = %.9f\nTotal DeltaV    dv = %.9f\nSemi Major Axis    = %.9f\n\n",dv0,dvf,totaldv,sma);
		
	/*debug
	printf("************************************************\n");
	printf("alpha = %.9f\n",alpha);
	printf("beta = %.9f\n",beta);
	printf("de = %.9f\n",de);
	printf("f = %.9f\n",f);
	printf("g = %.9f\n",g);
	printf("mu = %.9E\n",mu);

	f64 firstPartOfG = sqrt(sma * sma * sma / mu);
	f64 secondPartOfG = (de - sin(de));

	printf("firstPart = %.9f\n",firstPartOfG);
	printf("secondPart = %.9f\n",secondPartOfG);
	
	f64 firstTimesSecond = firstPartOfG * secondPartOfG;
	f64 timeMinusfirstTimesSecond = dt - test;

	printf("firstTimesSecond = %.9E\n",firstTimesSecond);
	printf("timeMinusfirstTimesSecond = %.9f\n",timeMinusfirstTimesSecond);

	printf("start X=%.9f, Y=%.9f, Z=%.9f\n",r0.X,r0.Y,r0.Z);
	printf("desti X=%.9f, Y=%.9f, Z=%.9f\n\n",rf.X,rf.Y,rf.Z);
	
	printf("start Length: %.9Ef\n",r0.getLength());
	printf("desti Length: %.9Ef\n\n",rf.getLength());
	*/
}

 

[Enjo]

Just a quick answer without doing much of a research: Make sure that you got masses of bodies absolutely right and precise. Orbiter's cfgs are good sources. We've had such problem with precision in TransX back in the old days, when the masses weren't quite right in Orbiter's cfgs, but the motion of planets were hard-coded as if the masses were OK.

Nicely looking app, BTW

 

[zillion42]

:stirpot:
Hmm, I'm cooking something and it's not always to my taste...

I figured that the lambert solver works just fine, no inverted vectors at all. The reason why it presensts me a total delta-V 108km/s solution sometimes, is that it chooses to take the other way around (always two ways to go on an ellipse) to arrive at rf(destination Pos) with rv(destination Speed), except for that, it works on the spot.

If I randomly, or lets say pseudo-system-randomly try it with different transfer periods it also works the right way around sometimes, also on the spot, unfortunately it's never really close to total 6-7 km/s delta-V (~3 departure & ~3 arrival) as it should, more like high energy transfers totaling in around 14 km/s. I could extend this system and calculate my way through until I finally find, and only by chance, the right configuration, but that wouldn't be the intended way things should work.

To make it more complicated:
There is one way I can get those total 6-7 km/s transfers (tricked in the right direction), which is when I use one of the inverted solutions and invert the velocity vector again (as I did before by *-1.0), which then unfortunately doesn't quite fit anymore, meaning that of course since it's an ellipse, I pass through rf, but either to late or to early...

Maybe I should add that it only works on the spot after I use a new cpu costly kepler function from David Vallado with double precision to plot the actual trajectory after burning to calculated velocity.
http://www.google.com/codesearch/p?...hrM/ast2body.cpp&q=vallado lambert lang:c&d=3

Now my question would be to one you transfer MFD programmers, is there a way to force the lambert solver to always chose the same direction. Or even better, could someone maybe sum up the steps required to find valid r0,v0,rf,vf with a valid transfer time so the solver keeps it's promises...

Your help would be greatly appreciated.

Thanks in advance

@Enjo
I use:
<Mercury mass="330.2e21" />
<Venus mass="4868.5e21" />
<Earth mass="5973.6e21" />
<Mars mass="641.85e21" />
<Jupiter mass="1898600e21" />
<Saturn mass="568460e21" />
<Uranus mass="86832e21" />
<Neptune mass="102430e21" />
<Pluto mass="13.105e21" />
<Sun mass="1989100000e21" />

 

[Piper]

You have it right by using Lambert's solution, that what I used in my Trajectory Planner. As for forcing it into the same direction every time, it shouldn't be doing that anyways. Lambert's solution works on the premise that for every Launch Date and every Arrival Date, there is only one unique solution. In other words, just because you can go one way around an orbit, it doesn't mean the other way will get you there.

As for the satellite reaching Mars, how do you calculate its position? Are you using Kepler elements for it, or using something like VSOP87? If you're not using VSOP87, I recommend you do. I'm not sure where they are on the top of my head, but if you search around you should be able to find the routines/functions for it online.

If you feel like making your own function for solving Lambert's problem, I recommend [ame="http://www.amazon.ca/Orbital-Mechanics-John-E-Prussing/dp/0195078349/ref=sr_1_2?ie=UTF8&s=books&qid=1266169907&sr=8-2"]Orbital Mechanics[/ame] by John E. Prussing and Bruce A. Conway. It has everything you could ever want to know about this stuff, and a lot more too. You need a solid understanding of Algebra and Calculus for it, but it sounds like you already do.

 

[zillion42]

wow and thanks a lot for your lengthy answer... I guess I get it figured correctly eventually, right now though, I'm still struggling with it. Please excuse my (maybe) lengthy answer concerning stuff you may already know, it's just that I found over the years that formulating my problem precisely helps me tremendously in actually understanding it.

To answer your questions:
The ephemeris calculation I do is based on the work of Stephen R. Schmitt and as much as I just tried to double check I cant find what kind of model he uses. From what I can find the calculation of the planetary ephemerides is based on the linear rate of change per century using the elements of the Epoch J2000. It varies from reality by a few minutes of arc at max, over a few millenia. I could easily switch to VSOP87, there is heap loads of source code around (probably will do that sooner or later) but it wouldn't really change a lot for the problem I have right now, unfortunately.

The lambert solver takes as input:

• departure and destination position vectors aswell as
• departure and destination velocity vectors aswell as a
• transfer period, all in SI units (km, km/s, s)

to find 2 velocity vectors one for departure and one for arrival. I prefer to refer to what I call total-V, since delta-V is merely the total-V minus the input velocity vector, or in other words the speed of the satellite minus the speed of earth for departure. Long talk, the point is that the planetary model precision is (a nice thing to have) but for this problem negligable since however precise it may be, I already know (and it wont change) where earth is at departure and mars at arrival. Or again in other words the precision of the ephemeris calculation doesnt concern the lambert solver in the slightest as long as the planets are where they are supposed to be at the time they are supposed to be there.

I hope I'm not wrong with that, but from my understanding part of the lambert problem is:

The geometrical problem to find all ellipses that go through the points P1 and and P2 and have a focus at the point F1

wikipedia
and depending on the the transfer angle alpha which can either be

0° < alpha < 180° or 180° < alpha < 360°

wikipedia
the orbital pole can either be:

r1 cross r2 or -r1 cross r2

in the 2nd case resulting in the same ellipse with the opposite direction of motion. So, yes there is only one solution for the shortest way depending on which way is shorter. I guess that's where I have to search for my error.

Which is another good example of talking about it and actually understanding my problem, and oh, how much I hate wikipedias mathematical descriptions.

btw. I have at home Planetary Sciences by Imke De Pater & Jack Jonathan Lissauer which is so so, and I am really really thinking about buying Fundamentals of astrodynamics and applications by David Vallado which comes with the Companion code for Fundamentals of Astrodyanmics and Applications which available in fortran, c++, pascal and matlab

Thx again and please feel free to add whatever comes to your mind.

EDIT:
so it looks like I'm going to have a nice printf debug evening with beer and constantly looking on these few lines of code, testing with different configurations

	if (dtheta < core::PI64)
	{
		tp = sqrt(2.0 / (mu)) * (pow(s, 1.5) - pow(s - c, 1.5)) / 3.0;
	}
	if (dtheta > core::PI64)
	{
		tp = sqrt(2.0 / (mu)) * (pow(s, 1.5) + pow(s - c, 1.5)) / 3.0;
		bflag = true;
	}
	if (dtheta == core::PI64)
	{
		printf(" dtheta = 180.0. Do a Hohmann\n");
		return;
	}

EDIT2:
on the other hand I'm not entirely sure if that concerns the inclination change or the transfer angle... I'll see
EDIT3:
nerve fodder acquired ready to go !

 

[Enjo]

I use:
<Mercury mass="330.2e21" />
<Venus mass="4868.5e21" />
<Earth mass="5973.6e21" />
<Mars mass="641.85e21" />
<Jupiter mass="1898600e21" />
<Saturn mass="568460e21" />
<Uranus mass="86832e21" />
<Neptune mass="102430e21" />
<Pluto mass="13.105e21" />
<Sun mass="1989100000e21" />

Well, that's imprecise. For example for Earth you'd need to use 5.973698968e+24.
Make the effort and do check those Orbiter configs

 

[zillion42]

I'll do that, but to be honest I dont have Orbiter installed atm., and it wont help my problem because right now I'm not concerned with gravity, except for that of the sun and there I use the value that came with the kepler function that is in use only for interplanetary transfers. Once I'll reach the spehere of incluence of mars, and especiallly after the lambert solver does what it's supposed to do, I'll switch back to my old newtonian vector addition gravity model which is a lot faster.
I'll be delighted by the way if you could copy those values for me, spares me from installing orbiter on my cranky old populated hard disk.

And for some updates, I went to bed at 2 yesterday frustrated since I couldn't really find a clear connection between the inverted transfers and the the transfer angle, so I came back to step 1 finding a launch day. Now in finding that, I replaced my old method of finding the headstart angle between earth and mars at launch time with the one I used before, frankly I'm little confused now, can someone maybe tell me which makes more sense, because they yield very similar but noticably different results.

EDIT:
Seems I've been so consumed with math I dont understand, that I even got the most basic things wrong. Both of them below are wrong, right would be:

marsHeadStartAngle = 180.0 - (360.0*(TransOrbitPeriod / marsFullOrbPeriod));

marsHeadStartAngle = TransOrbitPeriod * (360.0/(marsFullOrbPeriod/2.0));
or
marsHeadStartAngle = 180.0 + (PI * (1.0 - TransOrbitPeriod / marsFullOrbPeriod) * RADTODEG);

 

[Enjo]

Here you are, just 90 KB

View attachment OrbiterConfig.zip