Orbital mechanics of circularization

[RizzoTheRat]

I’m fiddling with a KOS script for circularization and can’t seems to get the physics quite spot on.

My intention is to do the circularization burn in a direction that means the Ap doesn’t move, which surely must be possible.

I’m calculating the required velocity for a circular orbit at the altitude of the Apoapsis, and subtracting my current velocity to get a dV vector.  I then figure there’s 3 accelerations acting on the ship, Thrust, Gravity, and Centripetal Acceleration.  So I’m calculating the thrust angle to be that required so that when you add on the vertical deceleration (gravity – centripetal) you get the direction of the dV vector.

It’s close but clearly not quite right as the Ap initially climbs and then drops back down, ending up around 120m below my target altitude, but a lot better than just following the dV vector which ends up around 500-600m low.  I'm not taking in to account increasing acceleration due to reduced mass on the burn time, so I actually start the burn a bit early, but as the Ap is dropping towards the end starting the burn a bit later would presumably be worse.

Any thoughts on what I’m missing?

 

Maths looks like this:

Required speed = Sqrt (mu / (kerbin radius + altitude)

dVel = Required tangential velocity – orbital velocity

Acc Gravity = mu / (radius + alttude)^2

Acc Centripetal = horizontal velocity ^2 / (radius +altitude)

Acc Vertical = Acc Gravity – Acc Centripetal

Thrust Acc = =max thrust / mass

I then use the sine rule to calculate the thrust angle from the thrust magnitude, vertical acceleration magnitude and angle of resultant thrust from the vertical.

 

[Superfluous J]

Your calculations assume thrust is done in one 0-second burst, but in fact it occurs over multiple seconds.

I have no idea how to correct for it, but suspect a lot of the hard work put into MechJeb is exactly for that.

[bewing]

Your horizontal velocity changes during the burn, so are you doing your centripetal calculation as an integral?

[RizzoTheRat]

I'm recalculating it quite a few times a second so it takes in to account changes in velocity during the burn, and updates the dV vector and the heading each time, so I think that should get around both those issues.

How accurate does Mechjeb get the circularization to compared to the target altitude?  It's a while since I've used it but I seem to remember it being pretty accurate.  But while Mechjeb can do the launch it's nowhere near as much fun as working it out myself, and it can't do an entire mission itself (I'm currently earning money flying a suborbital tourist mission while I type here, who needs to play computer games for fun when you can program them to play themselves  )

[Jacke]

The instantaneous dV is the zeroth approximation.   Using the vis-viva equation, get the coast orbit's velocity vector at the circularization altitude.  Then vector subtract the circular orbit velocity to get the dV vector (some amount slightly pointed down).  You should already have that.

The next approximation is to start the dV burn half the burn time before getting to the circularization altitude.  Given the dV in the zeroth approx and the effective exhaust velocity (Isp x g/zero), calculate its burn time.  Then you have to calculate the time you would normally get to the circularization altitude.  To do that to high precision for an elliptical orbit using the formulas for true and mean anomaly is complex but doable.  In this case, you want to calculate that time.  Given that time and the coast orbit velocity determined before, you can then make the linear approximation of the time and altitude of the burn start.

The next step I think would be to dynamically track the burn by tracking the current orbit velocity.  Once it becomes circular (no relative vertical velocity component), cut engines.

[Starman4308]

The simple solution: once you're close to target altitude, run a PID loop using pitch or throttle to control vertical velocity, keeping it at zero.

It's probably not the most efficient way to make a totally circular orbit, but it'll work.

[MrSystems]

I'll share my script if you can tell me how you're running KOS.  Steam automatically downloaded KSP version 1.6.1 and broke my KOS and I didn't have a backup.

This just creates the maneuver node. I have a different script (over 300 lines long) to execute the burn. The KOS github page has a much shorter version that works okay.

// circularizes at the next apsis
set burn_at_periapsis to TRUE.
if APOAPSIS > 0 { // if apoapsis is negative, time to apoapsis is infinite
	if ETA:APOAPSIS < ETA:PERIAPSIS { SET burn_at_periapsis to FALSE. }
}
if burn_at_periapsis { set which_apsis_text to "Periapsis". } else { set which_apsis_text to "Apoapsis". }

if burn_at_periapsis {
	set node_time to TIME:SECONDS + ETA:PERIAPSIS.
	set otherapsis to APOAPSIS.
	set burnapsis to PERIAPSIS.
} ELSE {
	set node_time to TIME:SECONDS + ETA:APOAPSIS.
	set otherapsis to PERIAPSIS.
	set burnapsis to APOAPSIS.
}

set v_old to sqrt(BODY:MU * (2/(burnapsis+BODY:RADIUS) -
                             1/SHIP:OBT:SEMIMAJORAXIS)).
set v_new to sqrt(BODY:MU * (2/(burnapsis+BODY:RADIUS) -
			     1/(BODY:RADIUS+burnapsis))).
set dv to v_new - v_old.

set MyNode to NODE(node_time, 0,0,dv).
ADD MyNode.

print "Setting up node at " + which_apsis_text.
print "ETA=" + ROUND(node_time - TIME:SECONDS,1) + ", dV=" + ROUND(dv,1).
print "Old eccentricity was " + ROUND(SHIP:ORBIT:ECCENTRICITY,4) + ", new eccentricity is " + ROUND(MyNode:ORBIT:ECCENTRICITY,4).

 

[RizzoTheRat]

12 hours ago, Jacke said:

The next approximation is to start the dV burn half the burn time before getting to the circularization altitude.  Given the dV in the zeroth approx and the effective exhaust velocity (Isp x g/zero), calculate its burn time.  Then you have to calculate the time you would normally get to the circularization altitude.  To do that to high precision for an elliptical orbit using the formulas for true and mean anomaly is complex but doable.  In this case, you want to calculate that time.  Given that time and the coast orbit velocity determined before, you can then make the linear approximation of the time and altitude of the burn start.

Bugger, I'm starting the burn at the wrong time having failed to take in to account the increased speed getting me there sooner!  KOS gives me a reading of time to Apoapsis and the velocity at a given time along my path.  I was starting the burn at BurnTime =(RequiredSpeedAtApoapsis-PredictedSpeedAtApoapsis)/(Thrust/Mass) before the Apoapsis.

I've modified that to StartTime = BurnTime * 2 * CurrentSpeed / (CurrentSpeed + RequiredSpeed) , so the burn should finish about when it reaches the current Apoapsis location...

..and now I end up 500m above the target altitude and the engine shuts off about 10 seconds after the original Apoapsis time, from a 25 second burn, whereas prevuiusly it was about 3 seconds before the Apo time.   So clearly that's still not right!  I'm only thrusting a degree or two away from the dV vector so there shouldn't be that much difference.

Ignoring the timing issue for a minute, there's definitely something wrong with my heading calculation  as the Ap climbs by 500-600m before dropping back down again so my initial thrust vector must need to point down a bit more.

 

 

I'm making life hard for myself by trying to work out the physics for everything rather than the use KSP's nodes or PID loops, although I am using basic feedback loops for things like Q limiting the climb.  I should be able to do that properly but can't really be bothered at the moment as I'm more interesting in the orbital mechanics than the aerodynamics.  I was using PID loops a bit on a previous game but found in some cases they're actually more complicated than working out the physics

 

Re versions, I've not played in a while, think I played on 1.4 but missed 1.5.  Downloaded 1.61 the other day, installed KOS, and all worked fine despite being the 1.41 compatibility version of KOS.

 

Code for those interested.

Function Circularize2 {
	Parameter ApOrPe.
	Set Thrott to 0.
	Set kuniverse:timewarp:mode to "PHYSICS".
	if ApOrPe ="Ap" {
		set nodetime to time:seconds+eta:apoapsis.
		set OrbSpd to sqrt(body:mu/(body:radius+ship:apoapsis)).
		Lock steering to velocityat(ship,nodetime):orbit. 
		set TargetAlt to ship:apoapsis.
		wait until vang(ship:facing:vector,velocityat(ship,nodetime):orbit)<1.
	}else if ApOrPe = "Pe"{
		set nodetime to time:seconds+eta:periapsis.
		set OrbSpd to sqrt(body:mu/(body:radius+ship:periapsis)).
		Lock steering to -velocityat(ship,nodetime):orbit.  
		Set TargetAlt to Ship:periapsis.
		wait until vang(ship:facing:vector,-velocityat(ship,nodetime):orbit)<1.			
	}else{
		print "Error, not selected Ap or Pe for circularization".
	}
	Print "Circularizing at : " + TargetAlt + " meters".	
	set nodespd to velocityat(ship,nodetime):orbit:mag.
	Set BurnTime to abs((OrbSpd-nodespd)/(Ship:Maxthrust/Ship:mass)).
	if nodetime-time:seconds>300{
		Set Warp to 4.
	}else{
		set warp to 1.
	}
	wait until nodetime-time:seconds <60 + Burntime.
	Set Warp to 1.
	wait Until nodetime-time:seconds<Burntime+10. 
	set warp to 0.	
	wait 1.	
	set dVel to ship:facing:vector.
	
	//Set StartTime to Burntime * 2 * vxcl(up:vector,velocity:orbit):mag / (sqrt(body:mu/(body:radius+Ship:altitude))+vxcl(up:vector,velocity:orbit):mag).
	
	Print "Start circularization burn at node minus "+burntime+ " seconds".
	Set shipheading to ship:facing:vector.
	lock steering to shipheading. 

	until dVel:mag<0.1 {
		if dVel:mag<50{
			Set Reqspeed to sqrt(body:mu/(body:radius+Ship:altitude)).
		} else {
			Set Reqspeed to sqrt(body:mu/(body:radius+TargetAlt)).
		}		
		Set HorizVec to vxcl(up:vector,velocity:orbit):normalized.       
		Set dVel to ReqSpeed*HorizVec-velocity:orbit. 
		Set GravAcc to (body:mu/(altitude + body:radius)^2).
		Set CentripetalAcc to (vxcl(up:vector,velocity:orbit):mag^2)/(altitude + body:radius).
		set VertAcc to GravAcc-CentripetalAcc.
		set dVangle to 90-vang(dVel,HorizVec).
		set ThrustAcc to ship:maxthrust/ship:mass.
		set OffAngle to arcsin(VertAcc*sin(dVangle)/ThrustAcc).
		Set ThrustAngle to 180-offangle-dVangle.
		
		Set ShipHeading to Heading (LaunchBearing,(90-ThrustAngle)):forevector.

		
		
		Print "   ReqSpeed: " + ReqSpeed + "        " at (0,20).
		print " HorizSpeed: " + vxcl(up:vector,velocity:orbit):mag + "        " at (0,21).
		print "      Delta: " + dVel:mag + "        " at (0,22).
		print "   Altitude: " + ship:altitude + "        " at (0,23).
		print "   Apoapsis: " + ship:apoapsis + "        " at (0,24).
		print "   Apo-Peri: " + (ship:apoapsis-ship:periapsis) + "       " at (0,25).
		print "    ETA Apo: " + eta:apoapsis + "         " at (0,26). 
		print " dVel Angle: " + vang(up:vector,dVel) + "     " at (0,7).
		print "    GravAcc: " + GravAcc + "     " at (0,28).
		print "Centripetal: " + CentripetalAcc + "    " at (0,29).
		print "    dVangle: " + dVangle + "    " at (0,30).		
		print "  ThrustAcc: " + ThrustAcc + "    " at (0,31).		
		print "   OffAngle: " + OffAngle + "    " at (0,32).		
		print "ThrustAngle: " + ThrustAngle + "    " at (0,33).	
		print "   NodeTime: " + (nodetime-time:seconds) at (0,34).
		
		if nodetime-time:seconds<BurnTime{
			Set Thrott to ThrottleDir(shipheading,0.25,dVel:mag/50,0.01).
		}
	}	
	Set Thrott to 0.
	Wait 1.
	print "Stable Orbit Achieved".

	
}


Function ThrottleDir{
	Parameter ReqVector.
	Parameter Accuracy.
	Parameter ThrottSet.
	Parameter MinThrott.
	if vang(ship:facing:vector,ReqVector)<Accuracy{
		Set Thrott to max(minThrott,ThrottSet).
	}else{
		Set Thrott to 0.
	}
	Return Thrott.
}

 

Edited February 2, 2019 by RizzoTheRat

[Jacke]

17 minutes ago, RizzoTheRat said:

I was starting the burn at BurnTime =(RequiredSpeedAtApoapsis-PredictedSpeedAtApoapsis)/(Thrust/Mass) before the Apoapsis.

I've modified that to StartTime = BurnTime * 2 * CurrentSpeed / (CurrentSpeed + RequiredSpeed) , so the burn should finish about when it reaches the current Apoapsis location...

Unless your circularization dv (in your terms "RequiredSpeedAtApoapsis - PredictedSpeedAtApoapsis") is small leading to dm (change in the mass of the spacecraft) being small, the BurnTime can't be accurately figured from the linear equation.  Use dV = c ln (m0/m1), the Rocket Equation, and change it to determine m1, the spacecraft mass after the burn, with c being the effective exhaust velocity (= Isp g0, specific impulse times standard acceleration due to gravity) and m0 the spacecraft mass before the burn.  Or

m1 = m0  exp(dv/c)

dm = m0 - m1

dm = mrate dt

mrate = T / c

T is the thrust and dt is the burn time.  Solving for dt

dt = ( m0 c / T) (1 - exp (- dv / c)

Replacing c with Isp and g0 which you have,

dt = (m0 Isp g0 / T) (1 - exp ( - dv / (Isp g0)))

And that's your burn time, from constants g0 (9.81 m/s/s) and spacecraft initial mass (m0), engine specific impulse (Isp), engine Thrust (T), and the circularization dv.

But that's the linear size of the vector dv.  You have to do the vector math to get the angle to face down.

[MrSystems]

9 hours ago, RizzoTheRat said:

Bugger, I'm starting the burn at the wrong time having failed to take in to account the increased speed getting me there sooner!  KOS gives me a reading of time to Apoapsis and the velocity at a given time along my path.  I was starting the burn at BurnTime =(RequiredSpeedAtApoapsis-PredictedSpeedAtApoapsis)/(Thrust/Mass) before the Apoapsis.

There are two things wrong with that approach. First, you want to execute half your delta-V before the apoapsis and half after. Because the burn takes some time, the part of the burn that happens before apoapsis changes your apoapsis, but the part that happens after, assuming you're pointing along the correct vector the entire time, changes it back to where it was when you started.

Second, your thrust/mass will change as you burn rocket fuel. You'll need to take that into account in deciding a) how long the burn will take (it's less than what thrust_now / mass_now would predict) and b) when to start the burn (more than half of the burn goes before the apoapsis).

Relevant code from my burn_node script, which I use every time I need to execute a node and not just for circularizing at the end of a launch, is:

// calculate the weighted ISP of all active engines (not strictly necessary if you only ever use one engine)

List ENGINES in MyEngines.
set current_stage to stage:number.
set thrust_multiplier to 1 // if my burn takes less than 10 seconds, I have code elsewhere to lower the throttle
set total_thrust to 0.
set total_flow to 0.
set loop_complete to FALSE.  // gotta leave time to complete the FOR loop before proceeding? takes time to count the engines?

set engine_count to 0.
for eng in MyEngines {
	if (eng:stage = current_stage) {
	  set total_thrust TO total_thrust + eng:availablethrust.
	  set total_flow TO total_flow + eng:availablethrust / eng:isp.
	  Print "Thrust=" + round(total_thrust) + ", mass flow=" + 1000*round(total_flow) + " for " + eng:name.
	}
	set engine_count to engine_count + 1.
}.

// total_flow is the rate at which mass flows through the engines
set loop_complete to TRUE.
wait until loop_complete AND engine_count >= MyEngines:LENGTH. // multi-threaded code...

// need to multiply by 1000 (because thrust is kN) and divide by g to get mass flow rate.
set total_flow to thrust_multiplier * total_flow * 1000 / 9.80665.
SET current_mass TO MASS.
SET current_accel TO thrust_multiplier * total_thrust / current_mass.
set weighted_ISP TO 1000*(total_thrust / (total_flow / thrust_multiplier)) / 9.80665.

set v_e  TO weighted_ISP * 9.80665.
set future_mass TO MASS / (constant:e)^(NextNode:DeltaV:mag / v_e). // mass after burn is executed
// mass and acceleration both change linearly with time, so average equals midpoint.
// also, time to burn is equal to change in mass divided by flow rate.

set burn_duration TO (1000*(MASS - future_mass) /  (total_flow)).
If thrust_multiplier < 1 {
  Print "Throttle-adjusted burn time is " + round(burn_duration,1).
} ELSE {
  Print "Burn time: " + (round(burn_duration,1) + TIME - TIME):CLOCK + " (" + round(burn_duration,1) + "s)".
}
                        
// for very long burns significantly more acceleration happens on the back half of the burn due to change in TWR
// the first half of dV is achieved when remaining mass = EXP((ln(initial mass) - ln(final mass))/2)
// that happens at time (initial mass - remaining mass) / flow rate
set first_half_burn_mass TO (constant:e)^((ln(MASS) + ln(future_mass))/2).
set first_half_burn_time TO 1000*(MASS - first_half_burn_mass)/(total_flow). // total flow is in kg, MASS is in metric tons.

Print "Half-burn: " + ((round(first_half_burn_time,1) + TIME - TIME)):CLOCK + " (" + round(first_half_burn_time,1) + "s)".
set target_burn_time to time:seconds + nd:eta - (first_half_burn_time). //nd:eta is the time until the burn node
  // if you're not using a node, it would be set target_burn_time to ETA:APOAPSIS - first_half_burn_time
Print "Burning in " + (target_burn_time - time):clock + " at " + (target_burn_time + TIME - TIME):clock.

That's just for calculating when and how long to burn. My code to execute the burn wouldn't help you with what you're trying to do because it relies on a maneuver node that you don't wish to create.

My process for putting ships into orbit, though, is:

• open a terminal and type 'run launch'
• wait about 4 minutes

My 'launch' script is about 440 lines. At the end it calls my 'circularize' script (pasted in full in my previous reply) and then my 'burn_node' script. 

[RizzoTheRat]

Cheers both, I'll have a play with getting the burn timing right and see if that improves things.

 

47 minutes ago, MrSystems said:

My process for putting ships into orbit, though, is:

• open a terminal and type 'run launch'
• wait about 4 minutes

My 'launch' script is about 440 lines. At the end it calls my 'circularize' script (pasted in full in my previous reply) and then my 'burn_node' script. 

That seems like a lot of effort, I just click the launch button from the VAB 

Boot script loads my functions to the onboard processor and tuns them, and then loads a script file of the same name as the vessel which has the parameters like initial turn velocity and angle, and then calls the relevant functions for the flight.

 

I find KOS to be something that's alternately as frustrating as hell and very rewarding...then again I suppose I could say the same about KSP

 

[MrSystems]

11 minutes ago, RizzoTheRat said:

Cheers both, I'll have a play with getting the burn timing right and see if that improves things.

Let us know if it works! I tried the vector addition method before I switched to just using orbital mechanics methods, and when I failed I assumed it was impossible.

[RizzoTheRat]

On 2/2/2019 at 11:25 AM, Jacke said:

dt = (m0 Isp g0 / T) (1 - exp ( - dv / (Isp g0)))

Still not got this quite right, I've tried to do in steps a bit like @MrSystems code, as I like to be able to follow code fairly easily, and then also done it by @Jacke's single equation.

Slow method:

VExhaust = ISP * 9.81

MassFinal = Mass * e^(-dV/VExhaust)

MeanAcc = 2 * Thrust / (Mass + MassFinal)

BurnTime = dV / MeanAcc

 

One step method:

dT = (Mass * ISP * 9.81/Thrust) * (1-e^(dv / (ISP*9.81)))

 

A test ship the pad is reporting 29.68 seconds for the slow method and 29.28 for the one step.  I can't see a mistake in Jacke's workings, so is my assumption that the rate of change of change of mass, and therefore acceleration, is linear wrong?

I'm also not 100% sure about units here.  MrSystems is multiplying by 1000 as thrust is in KN, but it appears mass is in tonnes, so that shouldn't be needed should it?

Edited February 3, 2019 by RizzoTheRat

[Jacke]

1 hour ago, RizzoTheRat said:

One step method:

dT = (Mass * ISP * 9.81/Thrust) * (1-e^(dv / (ISP*9.81)))

dT = (Mass * ISP * 9.81/Thrust) * (1-e^(- dv / (ISP*9.81)))

That minus sign in the exponent is vital.  Else the formula will return too low a value.

That's the scalar burn time.  You have to do the vector math of the velocity subtraction to get the vector dv and thus the angle down for the burn.

Edited February 3, 2019 by Jacke

[RizzoTheRat]

Oops, that was a typo on here, the code has it negative, so I'm still assuming I've got an assumption wrong on the slow method.

 

The vector bit is what I originally started the thread for, but got sidetracked when we discovered I'd messed up the burn time calculation

 

I'm subtracting the current velocity from the orbital velocity to get the required dV, and the angle of that velocity below the horizon (alpha)

The ship is still suborbital at this point so I assumed its accelerating straight down at a rate of Gravity - Centripetal Acceleration, in which case I need to offset the thrust angle so that the sum of thrust, gravity and centripetal ends up at the angle alpha.  Excuse the dodgy PowerPoint diagrams...

With my current small test ship and the updated timing this gets me to within around 100m...but thrusting on the dV vector gets me to within 10m.  I've just tried it with a bigger ship and dV vector go me about 800m low while the Gravity-Centripetal method got me about 1500m high.

So now I have no idea what's going on.

[Jacke]

Assuming you're trying to circularize at 100km, that's 700km semi-major axis (SMA), ie. the altitude of the desired orbit.  So your +800m high and 1500m low are about 0.1% to 0.2% error.

I think to get closer, use the best estimates to get the start of burn, but then split the problem into components: the vertical SMA and velocity and the horizontal velocity.  You want to zero the vertical velocity component at the desired SMA and change the horizontal velocity to the circular velocity at that SMA.  That's why the spacecraft angles down, to decrease the vertical velocity to zero at the SMA altitude and increase the horizontal velocity to the circular velocity.

Figure out the component accelerations to do that.

'0' indicating start of burn, '1' indicating end of burn, and 'd' indicating desired, given vertical velocity v0, altitude A0 and Ad,  horizontal velocities h0 and hc, and burn time dt, find vertical acceleration av and horizontal acceleration ah

av = v0^2 / (2 ( Ad - A0)) - g

ah = (hc - h0 ) / dt

And

av^2 + ah^2 = ( (2 Thrust ) / ( m0 + m1) )^2

The last is a linear approximation, but I think it'll get low error as long as the burn's not too long and the change in mass too much.  Take all those and solve for A0, the altitude at which to start the burn.

A0 = Ad - ( (v0^2) / 2 ) / sqrt ( (2 T / (m0 + m1) )^2 -  ( (hc - h0) / dt)^2 )

Calculate av and ah and then get the angle from

angle = arctan ( av / ah )

You perhaps should put in an override to cut throttle when vertical velocity is zero.  The orbit is circular at that point.

 

Edited February 3, 2019 by Jacke

[adsii1970]

Question has been moved from Gameplay Questions and Tutorials to Science and Spaceflight.

[RizzoTheRat]

26 minutes ago, Jacke said:

I think to get closer, use the best estimates to get the start of burn, but then split the problem into components: the vertical SMA and velocity and the horizontal velocity.  You want to zero the vertical velocity component at the desired SMA and change the horizontal velocity to the circular velocity at that SMA.  That's why the spacecraft angles down, to decrease the vertical velocity to zero at the SMA altitude and increase the horizontal velocity to the circular velocity.

This is the bit I'm struggling to get my head around. Left to it's own devices the vertical component will decrease to 0 at the Ap.  By adding more horizontal velocity I'm increasing the Centripetal acceleration, so I need to thrust below the horizontal to take that in to account.  I just can't quite make up my mind if I'm double counting it with the dV already pointing down, and the adding in the effect of gravity and centripetal acceleration. 

 

 

[Jacke]

I had to correct my approximate formula for A0 above.  Forgot to square one term and square-root another.

All of this is very approximate.  To solve this exactly, it's a power track--an accelerated trajectory--and to solve that to a high degree is best done with a Lagrangian method.  Which I was introduced to reading Von Braun's The Mars Project.  (And even he only did it for the long burn in the Trans Mars Injection, to get the angle around the Earth of the burn to set the lead to get the right ejection angle.) That was a long time ago and I'm quite rusty with the math.  I think more approximate methods will work here.

 

1 hour ago, RizzoTheRat said:

This is the bit I'm struggling to get my head around. Left to it's own devices the vertical component will decrease to 0 at the Ap.  By adding more horizontal velocity I'm increasing the Centripetal acceleration, so I need to thrust below the horizontal to take that in to account.  I just can't quite make up my mind if I'm double counting it with the dV already pointing down, and the adding in the effect of gravity and centripetal acceleration.

Since we're working in the frame of the rocket, it's a rotating frame, so it means it's centrifugal acceleration.  And you're right, it will be a factor.

Think of the vertical acceleration component.  We want that to change the vertical velocity to zero over the burn.  In the rotating frame, there's the local acceleration of gravity down, g.  And the centrifugal acceleration, ac, is up, h^2/(r0 + A), where h is the instantaneous horizontal velocity, r0 is the radius of the planet, and A is the instantaneous altitude.

g over the period of the burn will be roughly constant.  (r0 + A) will change little.  The centrifugal acceleration will increase as h^2 increases.

The thrust of the craft will change, but hopefully not much.  That means the total acceleration doesn't change much.  And because the burn is at mostly shallow angles, that means the horizontal acceleration doesn't change much.

What will change is the vertical acceleration, as the centrifugal acceleration will increase.  The plan for the burn is for that vertical acceleration to be about constant, which means as the centrifugal acceleration goes up the rocket will have to pitch down.

I think the autopilot mods set this up and then use a PID to drive the numbers at a steady rate towards the goal values at the right time.

For your case I think you can calculate the angle at burn start and the angle at burn end and just pitch the spacecraft through that angle.

angle start = arctan ( (av - g + (h0^2 / (r0 + A0))) / ah)

angle end = arctan ( (av - g + (hc^2 / (r0 + Ad))) / ah)

where av is the roughly constant desired vertical acceleration, g is the local acceleration due to gravity, h0 is the initial horizontal velocity, hc is the desired circular velocity, r0 is the radius of the planet, A0 is the initial altitude, Ad is the altitude of the desired orbit, and ah is the roughly constant desired horizontal velocity.

And here's how to calculate the local g, where g0 is the g at the surface of the planet.

g = g0 ( r0 / (r0 + A) )^2

And I hope this is a little less murky than muddy water.

 

 

Edited February 3, 2019 by Jacke

[RizzoTheRat]

Helpfully KOS can extract the Standard Gravitational Parameter mu directly from KSP so local acceleration due to gravity is really easy at mu/(r0+A)^2, and I can get the centripetal acceleration from HorizontalVelocity^2/(r0+A).  I think your angle calcs are a condensed version of current approach but I'll need to sit down and work it through to be sure.  The beauty of KOS is it can calculate on every physics tick so you can update these things on the fly rather than need to calculate the start and end and interpolate in between.

I was thinking last night about trying to burn in a direction to ensure the vertical component of acceleration is such that it hits 0m/s vertically at the end of the burn, but then it occurred to me there might be a simpler way.  If I start the burn at Ap, I could angle up slightly to counteract the effect of gravity-centripetal and keep the vertical acceleration to 0, this should mean the Ap stays at the target height as the ship accelerates.  I think this would be a bit less efficient* though, and still doesn't solve the issue of working out exactly what the hell is going on.

 

*efficiency is a pretty useful thing, I've managed to design an early craft that can get a pilot and 2 tourists to an 80km orbit with something like 100m/s remaining dv left when launched by my current script.  I tried flying it manually and couldn't make orbit.  It appears the most efficient launch trajectory is one where the fins overheat and explode just before main engine cutoff. 

 

 

Edited February 4, 2019 by RizzoTheRat

[Jacke]

Looks like you've got a good approach.

12 minutes ago, RizzoTheRat said:

..centripetal acceleration....

It's important to remember that centripetal acceleration is in the unrotated frame of reference, as in a frame of reference rigidly attached to the planet.  There, the centripetal acceleration is due to gravity and the rate of angular rotation of the rockets velocity vector due to that is reduced by moving faster.

This is case where I think with these methods, things are better understood in the rotating frame of reference, one attached to the spacecraft.  There you have the centrifugal acceleration caused by the horizontal velocity of the spacecraft as well as the acceleration of gravity down opposing it.  There's also the Coriolis force that provides sideways acceleration and the Euler force due to the change in angular velocity and in the opposite direction, but hopefully both of those will be small.

Moving the viewpoint between the 2 frames during any investigation has to be done carefully, but this approach can help get simpler formulae to manipulate.  Unfortunately, the Wikipedia articles are rather thick with vector math and thin on diagrams.

 

[RizzoTheRat]

I've got about 6 approaches but I've not managed to get any of them right yet

Not sure on the correct terminology but when I say Centripetal acceleration I just mean the acceleration away from Kerbin due to the ships velocity.  I'm subtracting that from gravity at the given altitude to get the total vertical acceleration

The rotating frame of reference thing is one cause of headaches.  For example I'm using the ships actual velocity vector prior to Apoapsis, but for the target velocity I'm using the orbital speed tangential to the planet at the point where the space craft currently is.  In reality this is going to be tilted down a bit as it's really the tangential velocity at a point further around the arc of the planet.  Again I'm not sure if I'm double counting as the change in this velocity is caused by gravity and the centripetal acceleration.  I'm thinking I should have a go at estimating the angle at apoapsis.  It's going to be pretty small unless it's a long burn time, but in a previous game I did have the my original script doing multi pass burns for low TWR ships going interplanetary, I wasn't too bothered about the Periapsis getting higher if I was a bit out, but then again if I can keep it lower for more Oberth effect would make for a more efficient burn.

I've not played KSP for a while but when I first started playing with KOS a couple of versions ago I found I'd forgotten so much maths I had to ask my wife, a maths teacher, for a crash course on vectors. She thought I was after it for work to start with rather than for a computer game

[Jacke]

2 hours ago, RizzoTheRat said:

Not sure on the correct terminology but when I say Centripetal acceleration I just mean the acceleration away from Kerbin due to the ships velocity.

CentriFugal forces and accelerations Fling away from the centre.  CentriPetal forces and accelerations Pull towards the centre.

Centripetal forces in astrodynamics are gravitational forces, so they're in the non-rotating frame of reference, that of the planet if it wasn't rotating (or the distant stars).

Centrifugal forces don't appear to have any source (let's no bring in Mach's Principle, things are complicated enough already), so they're in the rotating frame of reference, that of the rocket itself.

I'm going to have to think on this one.  Will have to go back to the basic problem and physics first principles.

What parameters of the coasting ellipse are available from KOS ?

[Ozymandias_the_Goat]

52 minutes ago, Jacke said:

CentriFugal forces and accelerations Fling away from the centre.  CentriPetal forces and accelerations Pull towards the centre.

Keep in mind that centrifugal forces are not a real force, they are merely a fictitious force that appears to justify the motion of an object in a circle, but in reality does not. 

[RizzoTheRat]

 

1 hour ago, Jacke said:

What parameters of the coasting ellipse are available from KOS ?

You can get the current velocity and altitude, Apo/Periapsis time and altitude, and the ships velocity vector at a given time on it's current path.  Although the latter is tricky as the frame of reference is fixed in relation to thee sun, so rotates relative to the ship and the planet  so you if you get a vector at a future point it won't be correct by the time you get there.

I've just discovered you can get ground speed and vertical speed which I usually calculate from the velocity vector, and it also seems to have angular velocity and angular momentum.