Код:

// Only for circular orbits! // Output: // phaseAngle - angle in radians (0 to 2pi) // phaseAngleTime - time of angle value (UT seconds) declare parameter initP,initV,destP,destV. set pi to constant():pi. set deg2rad to pi/180. // method 0 set phaseAngle to vang(initP,destP)*deg2rad. // method 1 //set phaseAngle to arccos(vdot(initP,destP)/(initP:mag*destP:mag))*deg2rad. set phaseAngleTime to time:seconds. if vdot(initV,initP-destP)>0 { set phaseAngle to 2*pi-phaseAngle. }

Код:

// For circular orbits only! // Output: // tpT - time // tpA - angular alignment (rad) // tpDV1 - delta-v to initial burn // tpDV2 - delta-v to final burn // tpAV1 - angular velocity at r1 orbit (rad) // tpAV2 - the same as above but for r2 declare parameter r1,r2,gravParam. set pi to constant():pi. set pi2 to 2*pi. set tpAV1 to sqrt(gravParam/r1^3). set tpAV2 to sqrt(gravParam/r2^3). set tpT to pi*sqrt((r1+r2)^3/(8*gravParam)). // method 0 // http://en.wikipedia.org/wiki/Hohmann_transfer_orbit set tpA to pi-tpAV2*tpT. set tpMethod to 0. // method 1 // http://forum.kerbalspaceprogram.com/threads/16511-Tutorial-Interplanetary-How-To-Guide //set tpA to pi-sqrt(gravParam/r2)*tpT/r2. set tpMethod to 1. // method 2 // https://docs.google.com/document/d/1IX6ykVb0xifBrB4BRFDpqPO6kjYiLvOcEo3zwmZL0sQ/edit //set tpA to pi-pi2*(0.5*(((r1+r2)/(2*r2))^1.5)). set tpMethod to 2. set tpDV1 to sqrt(gravParam/r1)*(sqrt(2*r2/(r1+r2))-1). set tpDV2 to sqrt(gravParam/r2)*(1-sqrt(2*r1/(r1+r2))).

Код:

declare parameter dest. set sps to "             ". print "--------------------------------------"+sps at (0,0). print " Hohmann transfer orbit demo.         "+sps at (0,1). print " Mission: Kerbin to "+(dest:name)+"   "+sps at (0,2). print "======================================"+sps at (0,3). print "     Initial R (km):                  "+sps at (0,4). print " Destination R (km):                  "+sps at (0,5). print "Angular algn.(grad):                  "+sps at (0,6). print "  Transfer time (s):                  "+sps at (0,7). print "    Burn 1 DV (m/s):                  "+sps at (0,8). print "    Burn 2 DV (m/s):                  "+sps at (0,9). print "     Total DV (m/s):                  "+sps at (0,10). print "  Init. AV (grad/s):                  "+sps at (0,11). print "  Dest. AV (grad/s):                  "+sps at (0,12). print "           Reserved:                  "+sps at (0,13). print " Phase angle (grad):                  "+sps at (0,14). print " Start point (grad):                  "+sps at (0,15). print "            ETA (s):                  "+sps at (0,16). print "      Angle formula:                  "+sps at (0,17). print "              DEBUG:                  "+sps at (0,18). print "--------------------------------------"+sps at (0,19). print "AG8 to make node. AG9 to exit.                    " at (0,20). print "If AG not working, try to toggle gears, lights or " at (0,21). print "brakes first. The kOS window must be unfocused.   " at (0,22). print "Seems AG not working in map mode.                 " at (0,23). print "--------------------------------------" at (0,24). set prpos to 21. lock destP to dest:position-body:position. lock destV to dest:velocity. lock initP to body:position*-1. lock initV to velocity:orbit. lock r1 to body:radius+(periapsis+apoapsis)/2. lock r2 to body:radius+dest:altitude. set pi to constant():pi. set pi2 to 2*pi. set rad2deg to 180/pi. set deg2rad to pi/180. if ag8="True" { toggle ag8. } if ag9="True" { toggle ag9. } set exit to 0. until exit=1 {     run getTransferParamsCO(r1,r2,body:mu).     run getPhaseAngleCO(initP,initV,destP,destV).     if r1>r2 {         set phaseAngle to (pi2-phaseAngle)*-1.         set startAngle to phaseAngle-tpA.         if startAngle>0 {             set startAngle to (pi2-startAngle)*-1.         }     }     if r1<r2 {         set startAngle to phaseAngle-tpA.         if startAngle<0 {             set startAngle to pi2+startAngle.         }     }         // method 0     //set t to startAngle/(pi2/obt:period-pi2/(pi2/tpAV2)).         // method 2     set t to startAngle/(tpAV1-tpAV2).         set x to round(r1/1000,3).              print x+sps at (prpos,4).     set x to round(r2/1000,3).              print x+sps at (prpos,5).     set x to round(tpA*rad2deg,8).          print x+sps at (prpos,6).         run fmtTime(tpT).                       print fmtTime+sps at (prpos,7).     set x to round(tpDV1,2).                print x+sps at (prpos,8).     set x to round(tpDV2,2).                print x+sps at (prpos,9).     set x to round(tpDV1+tpDV2,2).          print x+sps at (prpos,10).     set x to round(tpAV1*rad2deg,8).        print x+sps at (prpos,11).     set x to round(tpAV2*rad2deg,8).        print x+sps at (prpos,12).     set x to round(phaseAngle*rad2deg,8).   print x+sps at (prpos,14).     set x to round(startAngle*rad2deg,8).   print x+sps at (prpos,15).     run fmtTime(t).                         print fmtTime+sps at (prpos,16).     print tpMethod at (prpos,17).         if ag8="True" {         set x to node(phaseAngleTime+t, 0, 0, tpDV1).         add x.         toggle ag8.     }     if ag9="True" {         set exit to 1.         toggle ag9.     } }

Код:

// Format time in seconds to human-readable string. // Output: fmtTime declare parameter time2format. // time in seconds // 1 year = ??? days = ???????????????????????????? //            1  day = 6 hours = 360 mins = 21600 s //                     1  hour =  60 mins =  3600 s //                                1  min =     60 s set fmtD to floor(time2format/21600). set x to time2format-fmtD*21600. set fmtH to floor(x/3600). set x to x-fmtH*3600. set fmtM to floor(x/60). set x to x-fmtM*60. set fmtS to floor(x). set fmtTime to "D". if fmtD<100 { set fmtTime to fmtTime+"0". } if fmtD< 10 { set fmtTime to fmtTime+"0". } set fmtTime to fmtTime+fmtD+",". if fmtH< 10 { set fmtTime to fmtTime+"0". } set fmtTime to fmtTime+fmtH+":". if fmtM< 10 { set fmtTime to fmtTime+"0". } set fmtTime to fmtTime+fmtM+":". if fmtS< 10 { set fmtTime to fmtTime+"0". } set fmtTime to fmtTime+fmtS.

Орбитальная скорость

Орбитальная скорость круговой орбиты рассчитывается по формуле

v = sqrt(Mu / r)

где Mu - гравитационный параметр (равно GM, где G - гравитационная постоянная, M - масса планеты), r - расстояние между центрами масс. Например, для выхода на экваториальную орбиту на высоте 300км над поверхностью Кербина потребуется разогнаться до скорости

sqrt(3531600000000 / 900000) = 1980 m/s

Восхождение на орбиту можно разделить на три фазы:
1) выход за пределы плотных слоев атмосферы;
2) гравитационный маневр до выхода на орбиту;
3) корректировка орбиты;


Длительность прожига

First, some basic definitions:

m0 = Mass of spacecraft before burn (in metric tons, i.e. Mg).
m1 = Mass of spacecraft after burn.
dv = Change in velocity that will result from our burn (m/s).
ve = Our engines' effective exaust velocity (m/s) = g0 Isp.
g0 = Standard gravity = 9.81 m/s^2 (you can estimate this as 10, if needed).
Isp = Our engines' specific impulse (s) -- we'll assume a single type of engine.
dm = Mass of fuel expended during burn (Mg)
mdot = Fuel flow rate (Mg/s)
t = Duration of burn (s)
T = Total thrust of all firing engines (kN)

We can use the rocket equation to find the mass of the fuel that we burn (dm):

dv/ve = ln(m0/m1)
m0/m1 = e^(dv/ve)
dm = m0-m1 = m0 (1 - m1/m0) = m0 (1 - e^(-dv/ve))

Then, given the fuel flow rate (mdot), we can find out how long it takes us to burn that fuel (t).

dm = mdot t, so t = dm / mdot

We can find the fuel flow rate from our thrust (T):

T = mdot ve, so mdot = T / ve

Substituting the formulas for dm and mdot into the formula for t, we get the following formula for burn time:

t = (m0 ve / T) (1 - e^(-dv/ve))

That's the exact equation for the burn time. Don't let it seem too intimidating, because we're going to simplify it in a bit. One thing to notice at this point though, is that burn time is proportional to your starting mass, and inversely proportional to thrust. If you double your mass, the burn time doubles. If you double the thrust, the burn time is cut in half. This may seem obvious, but rocket science isn't always obvious -- for instance, notice that the burn time doesn't scale linearly with dv. That's because our ship is getting lighter as we burn, so our ship gets more bang for the buck (actually, it bucks more for the same amount of bang). That's the nature of the rocket equation.


Орбитальный период

Орбитальный период эллиптической орбиты определяется по формуле

t = 2PI*sqrt(a^3/mu)

где a - длинна большой полуоси, mu - гравитационный параметр планеты.


Фазовый угол

Для определения фазового угла можно использовать доступные в рамках kOS данные о положении объектов. Делать это следует только в статусе ORBITING.

В соответствии с документацией, для значений body:position за центр координатного пространства принимается положение корабля, на котором находится фокус. А расчет фазового угла необходимо выполнять относительно центра масс. В связи с этим, перед началом расчетов необходимо преобразовать позиции объектов таким образом, что бы получить координаты корабля и луны относительно центрального тела. Для этого, из координат луны следует вычесть координаты центрального тела, а координаты центрального тела инвертировать, что бы они стали соответствовать координатам корабля относительно общего центра. После этого, нужно нормализовать полученные векторы (привести к единичной длине) и найти угол между ними. Угол между двумя векторами находится по формуле

phi = arccos( p1*p2 / ( |p1| * |p2| ) )

где p1*p2 - скалярное произведение (dot product) векторов, представляющих положение корабля и луны относительно центра масс, |p1|*|p2| - произведение длин тех же самых векторов.

Для kOS угол между векторами можно получить посредством вызова функции VECTORANGLE или сокращенно VANG.

Решение проблемы потери знака:

// Getting the angle from a dot product tends to hide the sign, and I need the sign.
// (I need to know the difference between being 30 degrees ahead of the target vs
// 30 degrees behind it.)
// If my velocity is towards-ish the target then it's ahead of me.
// if my velocity is away-ish from the target, then it's behind me.
// I figure out which it is by looking at whether the dot product of my
// velocity and the target position relative to me is postive:
set chk to (lowVel:x*(highPos:x-lowPos:x) ) + (lowVel:y*(highPos:y-lowPos:y) ) + (lowVel:z*(highPos:z-lowPos:z) ) .
if chk < 0 {
set haPhi to 360-haPhi.
}.

Здесь lowVel это вектор скорости медленной орбиты, то есть более высокой.

Для объекта на более высокой орбите выражение фазового угла по отношению к объекту на более низкой орбите должно выражаться в интервале (-2PI, 0]. Для таких случаев, полученное выше значение фазового угла следует конвертировать следующим образом

phi = -2 * PI + phi

В результате фактически получаем угол опережения нагоняющего объекта. Чем ближе значение этого угла к нулю, тем ближе нагоняющий нас объект. Дополнительные разъяснения см. ниже в разделе объясняющем способ определения времени начала трансорбитального перехода.


Продолжительность трансорбитального перехода

Время перехода с одной круговой орбиты на другую определяется по формуле

tH = PI * sqrt( (r1 + r2)^3 / (8 * mu) )

где r1 и r2 - радиусы орбит, mu - гравитационный параметр центра орбит.


Угловая скорость

Угловая скорость в радианах рассчитывается по формуле

w = sqrt( mu / r^3)

где mu - гравитационный параметр центрального тела, r - радиус круговой орбиты.


Угловое выравнивание


Угловое выравнивание необходимо для расчета времени начала маневра трансорбитального перехода. Угловое расстояние перехода всегда равно PI. Таким образом, что бы найти угловую разницу, необходимо из PI вычесть расстояние, которое пройдет целевая орбита

a = PI - w * tH

где w - угловая скорость целевой орбиты, tH - продолжительность трансорбитального перехода.



Время начала трансорбитального перехода

Корректное определение момента первого прожига очень важно, с точки зрения точности попадания в нужную точку целевой орбиты. Прежде всего, нужно рассчитать текущее положение точки старта. Для перехода с нижней орбиты на верхнюю, точка старта находится раньше точки противостояния на расстоянии углового выравнивания

startA = phi - a

где phi - фазовый угол, a - величина углового выравнивания. Полученное значение нужно откорректировать, в случае, если оно вышло за область значений [0, 2PI). Если этого не сделать, то придется совершить лишний виток.

Формула, преобразующая угловое расстояние до точки старта во временное смещение выглядит следующим образом

startT = startA / (w1 - w2)

где w1 - угловая скорость начальной орбиты, w2 - скорость конечно орбиты. Таким образом, мы получили секунды, которые остаются от деления расстояния на скорость. Если мы просто разделим угол (угловое расстояние) на угловую скорость начальной орбиты, то мы получим время прохождения указанного угла, если бы цель стояла на одном месте. Но цель движется в том же направлении с постоянной скоростью. Скорость этого сближения рассчитывается как разность угловых скоростей. Эта разница и есть фактическая скорость прохождения нужного нам угла. Таким образом, разделив угловое расстояние на скорость сближения мы узнаем время, за которое достигнем точки начала маневра.

Что бы рассчитать время перехода с верхней орбиты на нижнюю посредством приведенных выше формул, требуются дополнительные действия. Прежде всего, следует обратить внимание на то, что разница угловых скоростей для перехода сверху-вниз будет отрицательная. И в результате расчета «как есть», мы получим отрицательное время, что при сложении с текущим укажет на время в прошлом. Так же следует обратить внимание на формулу определения углового выравнивания. Известно, что наш путь составит PI радиан в масштабе исходной орбиты. Исходя из этого у нас рассчитывается продолжительность перехода, которая в произведении на угловую скорость конечной орбиты дает нам угловое выравнивание. Если конечная орбита имеет более высокую угловую скорость (что характерно, если она ниже исходной), то за соответствующее время на ней будет пройден путь больший, чем PI радиан и в результате получится отрицательное значение углового выравнивания.

Что бы все вышеприведенные формулы без изменений работали для возвращения с верхней орбиты на нижнюю, необходимо изменить отношение к фазовому углу с «насколько мы отстаем от более медленного объекта», на «насколько мы опережаем более быстрый объект». Если в первом случае мы высчитываем угол, который нам осталось пройти впереди, то во втором нам необходимо знать какой угол у нас есть в запасе позади. А выражая угол «позади» в отрицательных величинах мы получаем тот самый минус, который даст нам положительный результат при делении угла до начала маневра на отрицательную разницу угловых скоростей. Проще всего организовать подобную инверсию непосредственно в момент определения фазового угла, воспользовавшись координатами расположения объектов для определения высшей орбиты. На языке kOS это может выглядеть примерно так:

if initP:mag>destP:mag { set phaseAngle to -2*pi+phaseAngle. }

где initP - позиция на исходной орбите, destP - позиция на целевой орбите, phaseAngle - текущий фазовый угол в радианах в пределах [0, 2PI).