Minor ascent/planet tweaks

ascent.py

@@ -63,7 +63,13 @@ def __init__(self,
         initialVelocity = None,
         launchInclination = 0,
         acceleration = None,
-        timestep = 1):
+        timestep = 1,
+        gravityTurnCurve = 1,
+        calculateExtraThrust = False,
+        dragCoefficient = None,
+        specificImpulse = None,
+        shipThrust = None,
+        endAngleDeg = 0):
         """
         Compute a climb slope for exiting the atmosphere and achieving orbit.
 
@@ -81,11 +87,14 @@ def __init__(self,
         with an angle that varies linearly with altitude.  The start and end points
         of the turn can be changed, the curve cannot.
 
-        There's no model of variable acceleration.  It would be easy
-        to add in since we are using numerical integration; you'd want to
-        take in a function from (deltaV, altitude, velocity) ->
-        acceleration to model jets.  By default, acceleration is assumed to be 2.2g,
-        which is good low down (except the first ~100m) but too little at altitude.
+        Variable acceleration is modeled if shipThrust (in kN) and specificImpulse
+        (in s) are specified. Initial ship mass is calculated by assuming that the given
+        acceleration is the maximum acceleration achievable at launch. At each step of
+        the ascent, the total delta-V achieved and the given specificImpulse are used to
+        reduce the ship's mass and thus increase its maximum acceleration.
+
+        By default, acceleration is assumed to be 2.2g, which is good low down (except
+        the first ~100m) but too little at altitude.
 
         Specify the timestep to change the accuracy; by default we
         use a timestep of 1s.  It is a fixed timestep.  Smaller timesteps
@@ -128,8 +137,10 @@ def __str__(self):
                         self.deltaV,
                         self.dragLoss,
                         self.thrust,
-                        "" if self.thrustLimited == 0 else
-                                (" needs %g m/s^2 more thrust" % self.thrustLimited)
+                        "" if not self.thrustLimited else
+                                (" needs %smore thrust" % (
+                                    "" if self.thrustLimited is True else
+                                        ("%g m/s^2" % self.thrustLimited)))
                     ))
 
         self.planet = planet
@@ -140,6 +151,12 @@ def __str__(self):
         else:
             v = list(initialVelocity)
 
+        if dragCoefficient is None:
+            dragCoefficient = planet.defaultDragCoefficient
+
+        # These angles are both measured from horizontal.
+        endAngleRad = endAngleDeg * math.pi / 180
+
         # p and v are the current position and velocity; we update these.
         # We use a coordinate system centered at the core of the planet,
         # aligned so that the initial altitude forms the y axis.
@@ -150,9 +167,9 @@ def __str__(self):
         t = 0                                    # s
         dragLoss = 0                             # m/s
         if acceleration is None:
-            a_thrust = 2.2 * planet.gravity(initialAltitude) # m/s^2
-        else:
-            a_thrust = acceleration
+            acceleration = 2.2 * planet.gravity(initialAltitude) # m/s^2
+
+        variable_accel = not (specificImpulse is None and shipThrust is None)
 
         TOA = planet.topOfAtmosphere() # m
         if orbitAltitude is None:
@@ -176,6 +193,14 @@ def __str__(self):
         # decide how much to burn and update our position.
         targetApoapsis = orbitAltitude + planet.radius
         climbSlope = []             # ClimbPoint list
+        loss_steering = 0
+        loss_drag     = 0
+        loss_gravity  = 0
+
+        # Keep track of approximately how much thrust would be needed to reach
+        # the desired apoapsis.
+        thrust_apo = None
+
         while h < targetApoapsis and t < 1000:
             if h < planet.radius:
                 # We crashed!  Bring more thrust next time.
@@ -199,15 +224,16 @@ def __str__(self):
             else:
                 # What fraction of the turn have we done?
                 ratio = (alt - gravityTurnStart) / (gravityTurnEnd - gravityTurnStart)
+                ratio **= gravityTurnCurve
                 # The more we did, the closer we want to be to pi/2 from
                 # the straight-up orientation.
-                phiSurf = - (ratio * math.pi / 2)
+                phiSurf = - (ratio * (math.pi / 2 - endAngleRad))
             # phi is the angle of thrust in the original coordinate frame.
             phi = phiSurf + psi
 
             # Compute the gravity and drag.
             a_grav = planet.gravity(alt)
-            a_drag = planet.drag(alt, L2(v))
+            a_drag = planet.drag(alt, L2(v), dragCoefficient)
 
             # Terminal velocity for our altitude:
             #print ("drag %g m/s^2, grav %g m/s^2, term %g m/s" % (a_drag, a_grav, v_term))
@@ -219,12 +245,24 @@ def __str__(self):
                 - (a_drag * sin(theta) + a_grav * sin(psi))
             )
             v_nothrust = [ v[i] + loss[i] * timestep for i in range(2) ]
-            def thrustResult(thrust):
+
+            def total_accel(pos, vel, thr):
+                altitude = L2(pos) - planet.radius
+                ag = planet.gravity(altitude)
+                ad = planet.drag(altitude, L2(vel), dragCoefficient)
                 return (
-                    v_nothrust[0] + thrust * cos(phi) * timestep,
-                    v_nothrust[1] + thrust * sin(phi) * timestep,
+                    thr * cos(phi) - (ad * cos(theta) + ag * cos(psi)),
+                    thr * sin(phi) - (ad * sin(theta) + ag * sin(psi))
                 )
 
+            def thrustResult2(thrust, tx, ty):
+                return (v_nothrust[0] + thrust * tx, v_nothrust[1] + thrust * ty)
+
+            def thrustResult(thrust):
+                tx = cos(phi) * timestep
+                ty = sin(phi) * timestep
+                return thrustResult2(thrust, tx, ty)
+
             # Compute the acceleration that gets us to terminal velocity in the
             # direction of thrust.
             # Let X be (cos phi, sin phi) => transform to the direction of thrust.
@@ -248,7 +286,7 @@ def thrustResult(thrust):
             def findTerminalThrust():
                 # above the top of the atmosphere, there is no limit!
                 if alt >= TOA: return 1e30
-                v_term = planet.terminalVelocity(alt)
+                v_term = planet.terminalVelocity(alt, dragCoefficient)
                 a = timestep * timestep
                 v_noTX = v_nothrust[0] * cos(phi) + v_nothrust[1] * sin(phi)
                 b = 2 * timestep * v_noTX
@@ -259,7 +297,7 @@ def findTerminalThrust():
             # Binary search the thrust to achieve the apoapsis.
             # I'm sure this could be worked out analytically.
             # Use the terminal velocity thrust to reduce our search space.
-            def findApoapsisThrust(thrust_term):
+            def findApoapsisThrust(thrust_term, guess = None):
 
                 # The most we could want to speed up is enough to immediately
                 # get our momentum to match what it will when we have a
@@ -280,14 +318,23 @@ def findApoapsisThrust(thrust_term):
                                (p[0]*cos(thetaSurf) + p[1]*sin(thetaSurf)))
                 thrust_targetMax = (v_targetMax - L2(v_nothrust)) / timestep
 
+                tx = cos(phi) * timestep
+                ty = sin(phi) * timestep
+
+                theta = math.atan2(p[1], p[0])
+
                 thrust_lo = 0
                 thrust_hi = min(thrust_targetMax, thrust_term)
                 # while we are off by more than 1mm/s, binary search
                 # we waste a bit of time searching for really big solutions.
                 while thrust_hi - thrust_lo > 1e-3:
-                    thrust = (thrust_hi + thrust_lo) / 2
-                    v_next = thrustResult(thrust)
-                    (apoapsis, _) = planet.determineOrbit2(p, v_next)
+                    if guess and guess < thrust_hi:
+                        thrust = guess
+                        guess = None
+                    else:
+                        thrust = (thrust_hi + thrust_lo) / 2
+                    v_next = thrustResult2(thrust, tx, ty)
+                    (apoapsis, _) = planet.determineOrbit3(h, theta, v_next)
 
                     # if apoapsis is too high, or negative (i.e. hyperbolic), reduce thrust
                     if apoapsis > targetApoapsis or apoapsis < 0:
@@ -296,23 +343,89 @@ def findApoapsisThrust(thrust_term):
                         thrust_lo = thrust
                 return thrust_lo # or hi, it's only 1mm/s difference
 
+            def getAcceleration(dv, alt):
+                acc = acceleration
+                if variable_accel:
+                    # dv = isp * g0 * ln(m0/m1)
+                    # m0 / m1 = math.exp(dv / (isp * g0))
+                    # 1 + dm / m0 = math.exp(dv / (isp * g0))
+                    # dm = m0 * (math.exp(dv / (isp * g0)) - 1)
+                    # m = m0 - dm
+                    # m = m0 * (1 - (math.exp(dv / (isp * g0)) - 1))
+                    # m = m0 * (2 - math.exp(dv / (isp * g0)))
+                    # m0 = shipThrust / acceleration
+                    # m = (2 - math.exp(dv / (isp * g0))) * shipThrust / acceleration
+                    # acc = shipThrust / m
+                    # acc = acceleration / (2 - math.exp(dv / (isp * g0)))
+                    acc /= (2 - math.exp(dv / (specificImpulse * 9.81)))
+                return acc
+
+            a_thrust = getAcceleration(dV, alt)
+
             thrust_term = findTerminalThrust()
-            thrust_apo = findApoapsisThrust(thrust_term)
+            guess = thrust_apo
+            thrust_limit = thrust_term
+            if not calculateExtraThrust:
+                # It's nice to know if thrust_apo exceeds possible thrust
+                # (a_thrust). However, determining thrust_apo precisely is very
+                # time-consuming, so let's limit it to 1 m/s^2 greater than
+                # a_thrust. This will still tell us if we need more thrust, it
+                # just won't tell us how much extra thrust we need.
+                thrust_limit = min(thrust_term, a_thrust + 1)
+            thrust_apo = findApoapsisThrust(thrust_limit, guess)
 
             thrust = min(thrust_term, thrust_apo)
+
             #print ("thrust: %g to reach term, %g to reach apo" % (thrust_term, thrust_apo))
             if thrust < a_thrust:
                 thrustLimit = 0
                 if thrust < 0:
                     thrust = 0
             else:
-                thrustLimit = thrust - a_thrust
+                thrustLimit = (thrust - a_thrust) if calculateExtraThrust else True
                 thrust = a_thrust
 
+            vmag = math.sqrt(v[0] ** 2 + v[1] ** 2)
+            vdot = 0 if not vmag else (v[0] * cos(phi) + v[1] * sin(phi)) / vmag
+            loss_steering += timestep * thrust * (1 - vdot)
+            loss_drag     += timestep * a_drag
+            loss_gravity  += timestep * math.fabs(sin(thetaSurf)) * planet.gravity(alt)
+
+            def vmult(val, vector):
+                return [val * x for x in vector]
+
+            def vadd(v1, v2):
+                return [v1[i] + v2[i] for i in range(len(v1))]
+
+            def update_rk4(p, v):
+                dv1 = vmult(timestep, total_accel(p, v, thrust))
+                dx1 = vmult(timestep, v)
+
+                dv2 = vmult(timestep, total_accel(vadd(p, vmult(0.5, dx1)), vadd(v, vmult(0.5, dv1)), thrust))
+                dx2 = vmult(timestep, vadd(v, vmult(0.5, dv1)))
+
+                dv3 = vmult(timestep, total_accel(vadd(p, vmult(0.5, dx2)), vadd(v, vmult(0.5, dv2)), thrust))
+                dx3 = vmult(timestep, vadd(v, vmult(0.5, dv2)))
+
+                dv4 = vmult(timestep, total_accel(vadd(p, dx3), vadd(v, dv3), thrust))
+                dx4 = vmult(timestep, vadd(v, dv3))
+
+                dx = vmult(1.0 / 6.0, vadd(vadd(vadd(dx1, vmult(2, dx2)), vmult(2, dx3)), dx4))
+                dv = vmult(1.0 / 6.0, vadd(vadd(vadd(dv1, vmult(2, dv2)), vmult(2, dv3)), dv4))
+
+                p = vadd(p, dx)
+                v = vadd(v, dv)
+                return (p, v)
+
+            def update_euler(p, v):
+                for i in (0,1): p[i] += v[i] * timestep
+                v = thrustResult(thrust)
+                return (p, v)
+
             # Update everything!
-            for i in (0,1): p[i] += v[i] * timestep
+            (p, v) = update_rk4(p, v)
             h = L2(p)
-            v = thrustResult(thrust)
+
             dV += thrust * timestep
             dragLoss += a_drag * timestep
             t += timestep
@@ -322,6 +435,9 @@ def findApoapsisThrust(thrust_term):
             # Timed out...
             raise BadFlightPlanException
 
+        self.loss_steering = loss_steering
+        self.loss_drag     = loss_drag
+        self.loss_gravity  = loss_gravity
         self._climbSlope = climbSlope
         self.planet = planet
         self.orbit = orbitAltitude