diff --git a/ascent.py b/ascent.py index 23da85f..b209806 100644 --- a/ascent.py +++ b/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 diff --git a/learn-ascent.py b/learn-ascent.py new file mode 100644 index 0000000..10e0b38 --- /dev/null +++ b/learn-ascent.py @@ -0,0 +1,297 @@ +import argparse +import math +import os +import random +import sys + +import ascent +import planet + +STABLE_ITERATIONS = 1000 + +# Using an end angle other than 0 is actually not very helpful. +VARY_END_ANGLE = False + +class Profile: + + # Class variables. + planet = None + alt0 = None + alt1 = None + accel = None + drag = None + + cache = {} + MAX_CACHE_SIZE = 10000 + + def __init__(self, gt0, gt1, curve, endAngle): + # Limit precision of values. + (gt0, gt1, endAngle) = [round(x, 2) for x in (gt0, gt1, endAngle)] + curve = round(curve, 3) + + self.gt0 = min(max(gt0, 0), self.alt1) + self.gt1 = min(max(gt1, 0), self.alt1) + self.curve = min(max(curve, 0), 1) + self.endAngle = min(max(endAngle, -10), 90) if VARY_END_ANGLE else 0 + + self.ascent = None + + self.score = self.cache.get((self.gt0, self.gt1, self.curve, self.endAngle), None) + if self.score is None: + if len(self.cache) == self.MAX_CACHE_SIZE: + del self.cache[random.choice(self.cache.keys())] + self._calc_score() + self.cache[(self.gt0, self.gt1, self.curve, self.endAngle)] = self.score + + self.generation = 1 + + @classmethod + def init(cls, planet, alt0, alt1, accel, drag): + cls.planet = planet + cls.alt0 = alt0 + cls.alt1 = alt1 + cls.accel = accel + cls.drag = drag + + @classmethod + def clear_cache(cls): + cls.cache = {} + + @classmethod + def from_string(cls, str): + tokens = [float(f) for f in str.strip().split(" ")] + if len(tokens) == 5: + tokens.append(-1) + (gt0, gt1, curve, endAngle, score) = tokens[:5] + return Profile(gt0, gt1, curve, endAngle) + + #9.72 23.87 0.3381 1119.067026 + @classmethod + def random(cls): + gt0 = random.random() * (cls.alt1 / 2) + gt1 = gt0 + random.random() * (cls.alt1 - gt0) + curve = random.random() * 2 + endAngle = random.random() * 90 + return Profile(gt0, gt1, curve, endAngle) + + def mutated(self): + vals = [self.gt0, self.gt1, self.curve, self.endAngle] + i = random.randint(0, len(vals) - 1) + vals[i] = self._mutate(vals[i]) + return Profile(vals[0], vals[1], vals[2], vals[3]) + + def _mutate(self, value): + amount = 0.1 + + # Mutate more slowly in later generations. + #amount /= math.log(self.generation + 1) + + return value + (random.random() - 0.5) * (math.sqrt(math.fabs(value)) if value else 1) * amount + + def _calc_score(self): + try: + self.ascent = ascent.climbSlope(self.planet, + orbitAltitude = self.alt1 * 1000, + gravityTurnStart = self.gt0 * 1000, + gravityTurnEnd = self.gt1 * 1000, + gravityTurnCurve = self.curve, + acceleration = self.planet.gravity() * self.accel, + initialAltitude = self.alt0 * 1000, + dragCoefficient = self.drag, + endAngleDeg = self.endAngle + ) + self.score = self.ascent.deltaV() + except ascent.BadFlightPlanException as bfpe: + self.score = -1 + self.ascent = None + + def _combine(self, a, b): + average = (a ** 2 + b ** 2) ** 0.5 + return self._mutate(average) + + def __lt__(self, other): + return self.score != -1 and self.score < other.score + + def combine(self, other): + return Profile(self._combine(self.gt0, other.gt0), + self._combine(self.gt1, other.gt1), + self._combine(self.curve, other.curve), + self._combine(self.endAngle, other.endAngle)) + + def better_than(self, other): + return (not other) or other.score < 0 or (self.score > 0 and self.score < other.score) + + def worse_than(self, other): + return (not other) or self.score > other.score + + def __str__(self): + return "%.2f %.2f %.3f %.2f %f" % (self.gt0, self.gt1, self.curve, self.endAngle, self.score) + + def guide(self): + angleStep = 15 + guide = "" + if self.curve: + guide = "angle altitude" + for angle in range(0, 90 + 1, angleStep): + dy = self.gt1 - self.gt0 + alt = self.gt0 + ((float(angle) / 90) ** (1.0 / self.curve)) * dy + guide += "\n%5d %8.1f" % (angle, alt) + return guide + + @classmethod + def desc_header(cls): + header = "start end shape" + if VARY_END_ANGLE: + header += " endAng" + header += " deltaV" + header += " loss_g atm steer" + return header + + def desc(self): + desc = "" + desc += "%5.2f %6.2f %5.1f " % (self.gt0, self.gt1, self.curve * 100) + if VARY_END_ANGLE: + desc += "%6.2f " % self.endAngle + desc += "%8.2f " % self.score + desc += "%7.2f %7.2f %6.2f" % (self.ascent.loss_gravity, self.ascent.loss_drag, self.ascent.loss_steering) + return desc + +def select(pool, s): + for profile in pool: + if random.random() > s: + return profile + return pool[-1] + +SILENT = True + +def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, poolSize = 20, fileIn = None, genLimit = None): + + p = planet.planets[planetName.lower()] + + if endAlt is None: + endAlt = math.ceil(p.topOfAtmosphere() / 5000) * 5 + + if not SILENT: + print("ascending on %s from %d to %d km" % (p, startAlt, endAlt)) + print("max acceleration: %.2f x surface gravity = %.2f m/s^2" % (accel, accel * p.gravity())) + if drag != 0.2: + print("drag coefficient: %.2f" % drag) + Profile.init(p, startAlt, endAlt, accel, drag) + + fileOut = "%s_%d_%d_%.2f_%.2f.txt" % (p.name, startAlt, endAlt, accel, drag) + + if fileIn is None: + fileIn = fileOut + + pool = [] + + if fileIn and os.path.exists(fileIn) and not SILENT: + with open(fileIn) as f: + for line in f.readlines(): + line = line.strip() + if line and not line.startswith("#"): + pool.append(Profile.from_string(line)) + + while len(pool) < poolSize: + profile = Profile.random() + pool.append(profile) + + bestEver = None + gen = 1 + lastChange = 0 + needNewline = False + + bestThisRound = None + if not SILENT: + print("%6s %s" % ("iter", Profile.desc_header())) + try: + while True: + best = None + worst = None + total = 0 + successes = 0 + candidates = [] + for profile in pool: + profile.generation = gen + if profile.better_than(best): + if profile.better_than(bestEver): + bestEver = profile + best = profile + if profile.better_than(bestThisRound): + lastChange = gen + bestThisRound = profile + + if profile.worse_than(worst): + worst = profile + + if profile.score > 0: + total += profile.score + successes += 1 + candidates.append(profile) + candidates.sort() + if successes and not SILENT: + lineOut = "\r%6d %s" % (gen, bestEver.desc()) + sys.stdout.write(lineOut) + sys.stdout.flush() + newPool = [] + SELECT_P = 0.5 + + # Automatically select the top candidate and a mutant of it. + if candidates: + newPool.append(candidates[0]) + newPool.append(candidates[0].mutated()) + + while len(newPool) < min(poolSize / 2, successes): + a = select(candidates, SELECT_P) + b = select(candidates, SELECT_P) + newPool.append(a.combine(b)) + + if gen >= lastChange + STABLE_ITERATIONS: + #print("\n%6d stable iterations; resetting..." % STABLE_ITERATIONS) + Profile.clear_cache() + newPool = [] + lastChange = gen + bestThisRound = None + + while len(newPool) < poolSize: + newPool.append(Profile.random()) + pool = newPool + + gen += 1 + if genLimit is not None and gen >= genLimit: + break + + except KeyboardInterrupt: + print("") + pool.append(bestEver) + pool.sort() + with open(fileOut, "w") as f: + for profile in pool: + f.write("%s\n" % profile) + print("") + print(bestEver.guide()) + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description = "Learn an ascent.") + args = [ + ('planet', 'planet', str, None, 'planet or moon'), + ('alt0', 'alt0', int, None, 'initial altitude (km)'), + ('alt1', 'alt1', int, None, 'final altitude (km)'), + ('-a', 'accel', float, 2.2, 'ship acceleration as a multiple of planet surface gravity (default: %(default)s)'), + ('-d', 'drag', float, 0.2, 'drag coefficient (default: %(default)s)'), + ('-p', 'poolSize', int, 2, 'pool size (default: %(default)s)'), + ('--profile', 'filename', str, None, 'profile one generation of execution and save results'), + ] + for (name, metavar, type, default, help) in args: + parser.add_argument(name, metavar=metavar, type=type, help=help, default=default) + + args = parser.parse_args(sys.argv[1:]) + + if args.profile: + random.seed(0) + import cProfile + SILENT = True + cProfile.run('learnAscent(args.planet, args.alt0, args.alt1, args.a, args.d, args.p, genLimit = 5)', args.profile) + else: + SILENT = False + learnAscent(args.planet, args.alt0, args.alt1, args.a, args.d, args.p) diff --git a/physics.py b/physics.py index abccefa..20c4464 100644 --- a/physics.py +++ b/physics.py @@ -34,4 +34,7 @@ def quadratic(a, b, c): return ( (-b + sqrtdiscriminant) / a2, (-b - sqrtdiscriminant) / a2) def L2(vector): - return math.sqrt(sum(x*x for x in vector)) + a = 0 + for x in vector: + a += x ** 2 + return math.sqrt(a) diff --git a/planet.py b/planet.py index ba486ae..3700f52 100644 --- a/planet.py +++ b/planet.py @@ -25,21 +25,18 @@ """ # The drag model in KSP 0.18.1 is -# F_{drag} = gamma D v^2 m -# where D is atmospheric density, v is the speed, m the mass (as a standin -# for cross-section). gamma combines a bunch of coefficients in one -# (including the 1/2 that would normally be there). It is not actually -# a constant, but it varies little over reported terminal velocities and over -# the parts we actually use (though aero components have lower drag). -coefficientOfDrag = 0.2 # Assumed constant, not really the case +# F_{drag} = gamma D P v^2 m +# where D is the coefficient of drag, P is atmospheric pressure, v is the speed, +# m the mass (as a standin for cross-section). gamma combines a bunch of +# coefficients in one (including the 1/2 that would normally be there). dragMultiplier = 0.008 kerbinSurfaceDensity = 1.2230948554874 -gamma = 0.5 * coefficientOfDrag * dragMultiplier * kerbinSurfaceDensity +gamma = 0.5 * dragMultiplier * kerbinSurfaceDensity class planet(object): - def __init__(self, name, gravityParam, SOI, radiusKm, siderealPeriod, datumPressure, scale): + def __init__(self, name, gravityParam, SOI, radiusKm, siderealPeriod, datumPressure, scale, ap, pe): self.name = name self.mu = gravityParam # m^3/s^2 self.SOI = SOI # m @@ -49,6 +46,10 @@ def __init__(self, name, gravityParam, SOI, radiusKm, siderealPeriod, datumPress self.siderealRotationSpeed = 2 * pi * self.radius / siderealPeriod # m/s self.datumPressure = datumPressure # atm self.scale = scale # m (pressure falls by e every scale altitude) + self.ap = ap + self.pe = pe + self.sma = (ap + pe) / 2 + self.defaultDragCoefficient = 0.2 def __str__(self): return self.name @@ -62,7 +63,18 @@ def gravity(self, altitude = 0): r = self.radius + altitude return self.mu / (r*r) - def orbitalVelocity(self, altitude): + def recircularize(self, alt1, alt2): + (alt1, alt2) = sorted([alt1, alt2]) + v1 = self.orbitalVelocity(alt1) + v2 = self.orbitalVelocity(alt1, ap = alt2) + v3 = self.orbitalVelocity(alt2, pe = alt1) + v4 = self.orbitalVelocity(alt2) + return (v2 - v1) + (v4 - v3) + + def impactVelocity(self, altitude): + return self.orbitalVelocity(0, ap = altitude, pe = -self.radius) + + def orbitalVelocity(self, altitude, ap = None, pe = None): """ Return the velocity required to maintain an orbit with the given altitude at the semimajor axis (e.g. a circular orbit at that @@ -70,7 +82,11 @@ def orbitalVelocity(self, altitude): altitude is in meters. """ - return sqrt(self.mu / (self.radius + altitude)) + if ap is None: ap = altitude + if pe is None: pe = altitude + r = self.radius + altitude + a = self.radius + (ap + pe) / 2 + return sqrt(self.mu * (2 / r - 1 / a)) def escapeVelocity(self, altitude): """ @@ -121,28 +137,31 @@ def determineOrbit(self, h, velocity): p = h * v * cos(theta) return quadratic(ainv, -2, p*p / self.mu) - def determineOrbit2(self, position, velocity): - """ - Given a position and velocity in some coordinate frame centered about - the center of the planet (doesn't matter, as long as 1 unit is 1 - meter), return the apsides, in meters above the core. - If one of them is negative, the orbit is hyperbolic. - """ - + def determineOrbit3(self, h, theta, velocity): + # h is distance from the planet's core. # theta is the angle of the position above the horizontal. # phi is the angle of the velocity above the horizontal. # psi is the angle of the velocity above the surface tangent. # psi = (90 - theta) + phi # thetabar = 90-theta - theta = math.atan2(position[1], position[0]) phi = math.atan2(velocity[1], velocity[0]) thetabar = (math.pi/2) - theta psi = thetabar + phi - h = L2(position) v = L2(velocity) return self.determineOrbit(h, (v, psi)) + def determineOrbit2(self, position, velocity): + """ + Given a position and velocity in some coordinate frame centered about + the center of the planet (doesn't matter, as long as 1 unit is 1 + meter), return the apsides, in meters above the core. + If one of them is negative, the orbit is hyperbolic. + """ + theta = math.atan2(position[1], position[0]) + h = L2(position) + return self.determineOrbit3(h, theta, velocity) + def hohmann(self, a1, a2): """ Given a circular orbit at altitude a1 and a target orbit at altitude @@ -226,7 +245,7 @@ def soiBurn(self, altitude, soiSpeed): return v - self.orbitalVelocity(altitude) - def terminalVelocity(self, altitude): + def terminalVelocity(self, altitude, dragCoefficient = None): """ Return the terminal velocity at a given altitude. This is the speed at which drag is as strong as gravity. @@ -234,14 +253,15 @@ def terminalVelocity(self, altitude): altitude is in meters. """ - # v = sqrt(g/ (datumPressure * gamma * e^{-altitude/scale})) + # v = sqrt(g/ (datumPressure * gamma * D * e^{-altitude/scale})) # Then pull the e term out, and remember that gravity changes # (slightly) with altitude. + if dragCoefficient is None: dragCoefficient = self.defaultDragCoefficient if altitude is None or altitude >= self.topOfAtmosphere(): return float("inf") return exp(0.5 * altitude / self.scale) * sqrt( - self.gravity(altitude) / (gamma * self.datumPressure) ) + self.gravity(altitude) / (gamma * dragCoefficient * self.datumPressure) ) - def drag(self, altitude, velocity): + def drag(self, altitude, velocity, dragCoefficient = None): """ Return the acceleration due to drag, in m/s^2. @@ -250,13 +270,15 @@ def drag(self, altitude, velocity): altitude is in meters, velocity in m/s """ - # F_{drag} = D v^2 m gamma, where D is atmospheric pressure, v is - # the speed, m the mass (as a standin for cross-section), and gamma - # combines a bunch of variables that are close enough to constant. - # To get the acceleration, divide by mass, which cancels out m: - # a_{drag} = gamma P v^2 + # F_{drag} = P v^2 m D gamma, where P is atmospheric pressure, v is + # the speed, m the mass (as a standin for cross-section), D is the + # coefficient of drag, and gamma combines a bunch of variables that are + # close enough to constant. To get the acceleration, divide by mass, + # which cancels out m: + # a_{drag} = gamma P v^2 D + if dragCoefficient is None: dragCoefficient = self.defaultDragCoefficient if altitude is None or altitude >= self.topOfAtmosphere(): return 0 - return gamma * self.pressure(altitude) * velocity * velocity + return gamma * dragCoefficient * self.pressure(altitude) * (velocity ** 2) def pressure(self, altitude): """ @@ -289,37 +311,40 @@ def topOfAtmosphere(self): """ # self.altitude with pressure = self.datumPressure * 1e-6, inlined: # log(1e6) ~ 13.81551... - return self.scale * 13.81551 + return self.scale * 13.81551 if self.scale else -self.radius # the sun has an infinite SOI; arbitrarily set it to 500x the orbit of Jool sunSOI = 500 * 71950638386 +# the sun doesn't orbit anything +sunAp = sunPe = 0 + planets = dict([ (p.name.lower(), p) for p in ( - # name, mu, SOI, radiusKm, sidereal, atm, scale - planet("Kerbol", 1.172332794E+18, sunSOI, 261600, 432000, 0, 0), - planet("Moho", 245250003655, 11206449 , 250, 1210000, 0, 0), - planet("Eve", 8171730229211, 85109365 , 700, 80500, 5, 7000), - planet("Gilly", 8289450, 126123.27, 13, 28255, 0, 0), - planet("Kerbin", 3.5316E+12, 84159286 , 600, 21600, 1, 5000), - planet("Mun", 65138397521, 2429559.1 , 200, 138984, 0, 0), - planet("Minmus", 1765800026, 2247428.4 , 60, 40400, 0, 0), - planet("Duna", 301363211975, 47921949 , 320, 65518, 0.2, 3000), - planet("Ike", 18568368573, 1049598.9 , 130, 65518, 0, 0), - planet("Dres", 21484488600, 32700000 , 138, 34800, 0, 0), - planet("Jool", 282528004209995, 2455985200 , 600, 36000, 15, 9000), - planet("Laythe", 1962000029236, 3723645.8 , 500, 52981, 0.8, 4000), - planet("Vall", 207481499474, 2406401.4 , 300, 105962, 0, 0), - planet("Tylo", 2825280042100, 10856518 , 600, 211926, 0, 0), - planet("Bop", 2486834944, 993002.8 , 65, 544507, 0, 0), - planet("Pol", 227905920, 1041613 , 44, 901903, 0, 0), - planet("Eeloo", 74410814527, 119087000 , 210, 19460, 0, 0), + # name, mu, SOI, radiusKm, sidereal, atm, scale, apoapsis, periapsis + planet("Kerbol", 1.172332794E+18, sunSOI, 261600, 432000, 0, 0, sunAp, sunPe), + planet("Moho", 245250003655, 11206449 , 250, 1210000, 0, 0, 6315765980, 4210510628), + planet("Eve", 8171730229211, 85109365 , 700, 80500, 5, 7000, 9931011387, 9734357701), + planet("Gilly", 8289450, 126123.27, 13, 28255, 0, 0, 48825000, 14175000), + planet("Kerbin", 3.5316E+12, 84159286 , 600, 21600, 1, 5000, 13599840256, 13599840256), + planet("Mun", 65138397521, 2429559.1 , 200, 138984, 0, 0, 12000000, 12000000), + planet("Minmus", 1765800026, 2247428.4 , 60, 40400, 0, 0, 47000000, 47000000), + planet("Duna", 301363211975, 47921949 , 320, 65518, 0.2, 3000, 21783189163, 19669121365), + planet("Ike", 18568368573, 1049598.9 , 130, 65518, 0, 0, 3296000, 3104000), + planet("Dres", 21484488600, 32700000 , 138, 34800, 0, 0, 46761053522, 34917642884), + planet("Jool", 282528004209995, 2455985200 , 600, 36000, 15, 9000, 71950638386, 65073282253), + planet("Laythe", 1962000029236, 3723645.8 , 500, 52981, 0.8, 4000, 27184000, 27184000), + planet("Vall", 207481499474, 2406401.4 , 300, 105962, 0, 0, 43152000, 43152000), + planet("Tylo", 2825280042100, 10856518 , 600, 211926, 0, 0, 68500000, 68500000), + planet("Bop", 2486834944, 993002.8 , 65, 544507, 0, 0, 129057500, 79942500), + planet("Pol", 720792000, 1041613 , 44, 901903, 0, 0, 210624206, 149155794), + planet("Eeloo", 74410814527, 119087000 , 210, 19460, 0, 0, 113549713200, 66687926800), ) ]) -del sunSOI +del sunSOI, sunAp, sunPe def getPlanet(name): return planets[name.lower()] # register the planets with the module, so users can write planet.kerbin or planet.getPlanet("kerbin") equivalently -for (name, p) in planets.iteritems(): +for (name, p) in planets.items(): globals()[name] = p