; "Лунолёт-103-12" версия 0.12.15, 5 декабря 2011 года ; 12-разрядная ветка, оптимизированная на скорость ; ; Математика и тестирование Станислава Володарского, Ленинград ; Программирование Ильи Васильева, Москва (оптимизация в процессе) ; Железо и компилятор Михаила Степанищева, Новосибирск ; С уважением к писателю Михаилу Пухову, инженеру Сергею Волкову, космонавту Юрию Глазкову ; и всем любителям ПМК ; ; Ra,10 gc - гравитационная постоянная небесного тела ; Rb,11 dt - шаг по времени (для блока "МЕТОД" вход) ; Rc,12 ra - модуль реактивного ускорения (функция времени) ; Rd,13 da - направление реактивного ускорения в полярных координатах ; Re,14 - адрес блока выбранного метода интегрирования "МАТЕМАТИКА" ; ; R16 m — масса топлива ; R17 mass — полный расход топлива на манёвр ; R18 DT — полное время манёвра ; R19 t — время, прошедшее от начала манёвра "МАТЕМАТИКА" ; ; R20 mp.x - горизонтальная координата (угол места в радианах) "МЕТОД" вход ; R21 mp.y - радиус (расстояние до центра планеты) "МЕТОД" вход ; R22 mp.vx - проекция скорости на местную горизонталь (трансверсальная скорость) "МЕТОД" вход ; R23 mp.vy - проекция скорости корабля на местную вертикаль (радиальная скорость) "МЕТОД" вход ; ; R24 rax - трансверсальное реактивное ускорение "ФИЗИКА" ; R25 ray - радиальное реактивное ускорение "ФИЗИКА" ; ; R26 p0.x — вектор координаты корабля "МАТЕМАТИКА" ; R27 p0.y ; R28 p0.vx — вектор скорости корабля "МАТЕМАТИКА" ; R29 p0.vy ; ; R30 q.x — "МЕТОД" выход, результат LinEx h(), аргумент Fizika g() ; R31 q.y ; R32 q.vx ; R33 q.vy ; ; R34 exy — желаемая ошибка координат на секунду полёта "МАТЕМАТИКА" ; R35 ev — желаемая ошибка скорости на секунду полёта "МАТЕМАТИКА" ; ; R36 p1.xy.x — вектор координаты xy корабля на следующем шаге "МАТЕМАТИКА" ; R37 p1.xy.y ; R38 p1.v.vx — вектор скорости v корабля на следующем шаге "МАТЕМАТИКА" ; R39 p1.v.vy ; ; R40 M — сухая масса корабля (без топлива) ; R41 C — скорость истечения горючего ; R42 R — радиус планеты ; R43 dt_min — минимальный шаг по времени "МАТЕМАТИКА" ; ; R44 rerxy — относительная ошибка координат корабля (оценка) "МАТЕМАТИКА" ; R45 rerv — относительная ошибка скорости корабля (оценка) "МАТЕМАТИКА" ; ; R46 p21.x — вектор координаты корабля на половинном шаге "МАТЕМАТИКА" ; R47 p21.y ; R48 p21.vx — вектор скорости корабля на половинном шаге "МАТЕМАТИКА" ; R49 p21.vy ; ; R50 p2.xy.x — вектор координаты xy корабля через два половинных шага "МАТЕМАТИКА" ; R51 p2.xy.y ; R52 p2.v.vx — вектор скорости v корабля через два половинных шага "МАТЕМАТИКА" ; R53 p2.v.vy ; ; R54 ∆t — шаг по времени "МАТЕМАТИКА" ; R55 rer — временная переменная "МАТЕМАТИКА" ; R56 mt — текущее время внутри манёвра "МЕТОД" вход ; R57 tt — текущее время внутри манёвра "ФИЗИКА" вход ; ; R58 coeff — коэффициент, зависящий от метода "МАТЕМАТИКА" ; R59 mdt — шаг по времени, аргумент функции h(mp, mdt, p'), подпрограмма LinEx ; ; R60 p5'.vx' — проекция полного ускорения p5' на местную горизонталь "МЕТОД ДОРМАНД-ПРИНС" ; R61 p5'.vy' — проекция полного ускорения p5' на местную вертикаль "МЕТОД ДОРМАНД-ПРИНС" ; R62 p6'.vx' ; R63 p6'.vy' ; R64 p7'.vx' // убрано при оптимизации ; R65 p7'.vy' // убрано при оптимизации ; ; R70 p1'.vx' — проекция полного ускорения p1' на местную горизонталь "МЕТОД РУНГЕ-КУТТА" ; R71 p1'.vy' — проекция полного ускорения p1' на местную вертикаль "МЕТОД РУНГЕ-КУТТА" ; R72 p2'.vx' — проекция полного ускорения p2' на местную горизонталь "МЕТОД РУНГЕ-КУТТА" ; R73 p2'.vy' — проекция полного ускорения p2' на местную вертикаль "МЕТОД РУНГЕ-КУТТА" ; R74 p3'.vx' — проекция полного ускорения p3' на местную горизонталь "МЕТОД РУНГЕ-КУТТА" ; R75 p3'.vy' — проекция полного ускорения p3' на местную вертикаль "МЕТОД РУНГЕ-КУТТА" ; R76 p4'.vx' — проекция полного ускорения p4' на местную горизонталь "МЕТОД ДОРМАНД-ПРИНС" ; R77 p4'.vy' — проекция полного ускорения p4' на местную вертикаль "МЕТОД ДОРМАНД-ПРИНС" ; ; R78 p'.vx' — результат Fizika g(), аргумент LinEx h() ; R79 p'.vy' ; ; R80 p1'.x' — полная угловая скорость p1' "МЕТОД РУНГЕ-КУТТА" ; R81 p1'.y' — проекция полной скорости p1' на местную вертикаль "МЕТОД РУНГЕ-КУТТА" ; R82 p2'.x' ; R83 p2'.y' ; R84 p3'.x' ; R85 p3'.y' ; R86 p4'.x' ; R87 p4'.y' ; ; R88 p'.x' — результат Fizika g(), аргумент LinEx h() ; R89 p'.y' ; ; R90 p5'.x' — проекция полной скорости p5' на местную горизонталь "МЕТОД ДОРМАНД-ПРИНС" ; R91 p5'.y' — проекция полной скорости p5' на местную вертикаль "МЕТОД ДОРМАНД-ПРИНС" ; R92 p6'.x' ; R93 p6'.y' ; R94 p7'.x' // убрано при оптимизации ; R95 p7'.y' // убрано при оптимизации ; ; R99 таймер, измерят скорость работы ЭКВМ ; ; R201…226 — коэффициенты таблицы Дорманд-Принс ; ; R0300 o.p ; R0301 ecos ; R0302 esin ; R0303 o.e ; R0304 o.xperi — перицентр ; R0305 apo — апоцентр ; R0306 theta — угловая кооррдината на кеплеровой орбите, измеряется от перицентра ; ; R0307 x — функция Arg ; R0308 y — функция Arg ; R0309 sx — знак x, функция Arg ; R0310 sy — знак y, фукнкция Arg ; ; R0311 l — Root ; R0312 r — Root ; R0313 M — Root, Root1, Root2 ; R0314 m — Root ; R0315 M1 — Root1 ; R0316 E — equation() в Root, Root2; E2theta, E2t ; R0317 C — E2t ; R0318 t — T2E ; R0319 E — FreeFlight ; R0320 E — LandingQ ; R0321 E1 — LandingQ ; R0322 peri — LandingQ ; ; R0325 theta1, E1 — по аналогии с M1 ; ; Байтовые регистры: ; R1000 iMetod — 0 Эйлер, 1 Рунге-Кутт, 2 Дорманд-Принс "МАТЕМАТИКА" ; R1001 k — порядок метода интегрирования (1 Эйлер, 4 Рунге-Кутт, 4 Дорманд-Принс) "МАТЕМАТИКА" ; .ORG 0 PGOTO Init PGOTO LandingQ PGOTO FreeFlight ;0. Начальная инициализация p0, exy, ev, dt_min, gc, C, R, M, m Init: PGSB TajmerStart ; Засечь время Cx PM26 ; p0.x 0 PM29 ; p0.vy 0 PM28 ; p0.vx ; 1738000 ; (Луна), значение из 1980-х 1737400 ; Средний экваториальный радиус Луны по Куликовскому PM42 ; R PM27 ; p0.y ; В 1980-х квадрат радиуса умножался на ускорение свободного падения на поверхности ; PRM42 Fx^2 1,62 * ; (Луна) ; Значительно более точное значение гравитационной постоянной Луны можно получить через ; геоцентрическую гравитационную постоянную GE и отношение масс Луны и Земли μ 3,986005 EE 14 ENT ; GE = 3,986005·10^14 м³/с² 81,30068 ; μ = 1/81,30068 / ; K = μGE Ma ; gc 3660 PM41 ; C 2150 PM40 ; M 3500 PM16 ; m 7 +/- F10^x PM34 PM35 ; exy, ev 0,000001 PM43 ; dt_min 2 PPM1000 ; iMetod := 2 (Дорманд-Принс) 200 M4 ; Коэффициенты таблицы Дорманд-Принс ; 5 F1/x KM4 ; R200 := c21 = 1/5 /* вшит в код */ 3 ENT 40 / KM4 ; R201 := c31 = 3/40 9 ENT 40 / KM4 ; R202 := c32 = 9/40 44 ENT 45 / KM4 ; R203 := c41 = 44/45 56 ENT 15 / +/- KM4 ; R204 := c42 = -56/15 32 ENT 9 / KM4 ; R205 := c43 = 32/9 19372 ENT 6561 / KM4 ; R206 := c51 = 19372/6561 2536 ENT 218,7 / +/- KM4 ; R207 := c52 = -25360/2187 64448 ENT 6561 / KM4 ; R208 := c53 = 64448/6561 212 ENT 729 / +/- KM4 ; R209 := c54 = -212/729 9017 ENT 3168 / KM4 ; R210 := c61 = 9017/3168 355 ENT 33 / +/- KM4 ; R211 := c62 = -355/33 46732 ENT 5247 / KM4 ; R212 := c63 = 46732/5247 49 ENT 176 / KM4 ; R213 := c64 = 49/176 5103 ENT 18656 / +/- KM4 ; R214 := c65 = -5103/18656 35 ENT 384 / KM4 ; R215 := c71 = 35/384 0 KM4 ; R216 := c72 = 0 5 ENT 11,13 / KM4 ; R217 := c73 = 500/1113 125 ENT 192 / KM4 ; R218 := c74 = 125/192 2187 ENT 6784 / +/- KM4 ; R219 := c75 = -2187/6784 11 ENT 84 / KM4 ; R220 := c76 = 11/84 51,79 ENT 576 / KM4 ; R221 := c81 = 5179/57600 0 KM4 ; R222 := c82 = 0 7571 ENT 16695 / KM4 ; R223 := c83 = 7571/16695 39,3 ENT 64 / KM4 ; R224 := c84 = 393/640 920,97 ENT 3392 / +/- KM4 ; R225 := c85 = -92097/339200 1,87 ENT 21 / KM4 ; R226 := c86 = 187/2100 ; 40 F1/x KM4 ; R227 := c87 = 1/40 /* вшит в код */ ;------------------------------------ ; БЛОК ВВОДА-ВЫВОДА ; Развитие "Лунолёта-3", добавлены индикации скоростей и приглашений ; (угол отклонения вектора тяги от вертикали, градусы) С/П (расход топлива, кг) С/П (время, с) С/П Manevr: PGSB TajmerStop ; Измерили время работы Cx PPM9045 ; Градусы .NUM msgUgol PPM9026 ; Вывести приглашение RMa PRM27 / FSQRT ; Y: первая космическая скорость на данной высоте PRM27 PRM42 - ; X: p0.y - R текущая высота полёта PRM29 ; T: p0.vy вертикальная скорость PRM28 ; Z: p0.vx горизонтальная скорость FR FR R/S ENT FPI * 180 / Md ; da := RX * pi / 180 .NUM msgMass PPM9026 ; Вывести приглашение Cx FR R/S PM17 ; Ввести mass .NUM msgDT PPM9026 ; Вывести приглашение Cx FR R/S PM18 ; Ввести DT PGSB TajmerStart ; Начать отсчёт времени работы ЭКВМ ;------------------------------------ ; РАСЧЁТНЫЙ БЛОК "МАТЕМАТИКА" ; совершает один манёвр, заданный пользователем ;2. Начальная инициализация Cx PM19 ; t := 0 PRM18 PM54 ; ∆t := DT 1 PPM9045 ; Радианы ; Заранее установим k, coeff и Re для метода Рунге-Кутта ; Коэффициент для РК будет 16/15. ; Общая формула: coeff = 2^k/(2^k - 1), где k - порядок метода (4 для РК, 1 для Э) ; 4 PPM1001 ; k := 4 16 ENT 15 / PM58 ; coeff := 16/15 .NUM MetRK Me ; Re := адрес блока "МЕТОД РУНГЕ-КУТТА" Shag: PPRM1000 ; iMetod Px=0 NeEjler ; Установка метода Эйлера (iMetod = 0) .NUM MetEjler Me ; Re := адрес блока "МЕТОД ЭЙЛЕРА" 1 PPM1001 ; k := 1 2 PM58 ; coeff := 2 PGOTO Shag1 NeEjler: 1 - ; Если iMetod = 1, уже всё настроено под "МЕТОД РУНГЕ-КУТТА" Px!=0 Shag1 ; Вызов метода Дорманд-Принс (iMetod = 2) ; ; ДП возвращает два вектора p1 - решение, p2 нужен только для расчёта ошибки: ; ;stepAndError(p0, t, ∆t) { ; (p1, p2) = dormandPrince(p0, t, ∆t) ; return (p1, |p1 - p2|) ;} PRM26 PM20 PRM27 PM21 PRM28 PM22 PRM29 PM23 ; mp := p0 PRM19 PM56 ; mt := t PRM54 Mb ; dt := ∆t PGSB DorPrinc ; (p1, p2) := m(mp, mt, dt) PGOTO ShagProv Shag1: ; Пункты 3…5 реализуют следующую inline-функцию для методов Эйлера и Рунге-Кутта: ; ;stepAndError(p0, t, ∆t) { ; p1 := m(p0, t, ∆t) ; p21 := m(p0, t, ∆t/2) ; p2 := m(p21, t + ∆t/2, ∆t/2) ; return (p1, |p1 - p2|) ;} ; ; ;3. Вычисление p1 := m(p0, t, ∆t) // шаг PRM26 PM20 PRM27 PM21 PRM28 PM22 PRM29 PM23 ; mp := p0 PRM19 PM56 ; mt := t PRM54 Mb ; dt := ∆t KGSBe ; q := m(mp, mt, dt) PRM30 PM36 PRM31 PM37 PRM32 PM38 PRM33 PM39 ; p1 := q ;4. Вычисление p21 := f(p0, t, ∆t/2) // первые пол-шага ; PRM26 PM20 PRM27 PM21 PRM28 PM22 PRM29 PM23 ; mp := p0 (уже записано) ; PRM19 PM56 ; mt := t (уже записано) PRM54 2 / Mb ; dt := ∆t/2 KGSBe ; q := m(mp, mt, dt) PRM30 PM46 PRM31 PM47 PRM32 PM48 PRM33 PM49 ; p21 := q ;5. Вычисление p2 := f(p21, t+∆t/2, ∆t/2) // вторые пол-шага PRM46 PM20 PRM47 PM21 PRM48 PM22 PRM49 PM23 ; mp := p21 PRM19 PRM54 2 / + PM56 ; mt := t + ∆t/2 ; PRM54 2 / Mb ; dt := ∆t/2 (уже записано) KGSBe ; q := m(mp, mt, dt) PRM30 PM50 PRM31 PM51 PRM32 PM52 PRM33 PM53 ; p2 := q ; // rerxy, rerv - относительные ошибки, безразмерные ;6. Вычисление rerxy := |p2.xy-p1.xy| / (exy * ∆t) ; // оцениваем ошибку координат erxy как sqrt((p1.x*p1.y - p2.x*p2.y)^2 + (p1.y - p2.y)^2) ShagProv: PRM36 PRM37 * PRM50 PRM51 * - Fx^2 ; (p1.x * p1.y - p2.x * p2.y)^2 PRM37 PRM51 - Fx^2 ; (p1.y - p2.y)^2 + FSQRT PRM34 PRM54 * / PM44 ; rerxy := sqrt(..) / (exy * ∆t) ;7. Вычисление rerv := |p2.v-p1.v| / (ev * ∆t) ; // оцениваем ошибку скорости erv как sqrt((p1.vx - p2.vx)^2 + (p1.vy - p2.vy)^2) PRM38 PRM52 - Fx^2 ; (p1.vx - p2.vx)^2 PRM39 PRM53 - Fx^2 ; (p1.vy - p2.vy)^2 + FSQRT PRM35 PRM54 * / PM45 ; rerv := sqrt(..) / (ev * ∆t) ;8. Вычисляем максимум, коэффициент зависит от метода PRM44 PRM45 Kmax PRM58 * PM55 ; rer := coeff * max(rerxy, rerv) ;9. Если ∆t ≤ dt_min, перейти к пункту 11 // такой маленький шаг объявим точным PRM43 PRM54 - ; dt_min - ∆t Px<0 Tochno ;10. Если ошибка слишком грубая, перейти к пункту 14 для уточнения шага ;// делаем шаг, если можно ; if rer < 1 ; then t := t + ∆t, p0 := p1 PRM55 1 - ; rer - 1 Px<0 KorShag ;*** Точность нас устроила, шаг верен *** ;11. Если p1.y ≤ R, перейти к пункту 16 // врезались! Tochno: PRM42 PRM37 - ; R - p1.y Px<0 Posadka Vzlet: ;12. Полетели! PRM36 PM26 PRM37 PM27 PRM38 PM28 PRM39 PM29 ; p0 := p1 PRM19 PRM54 + PM19 ; t := t + ∆t ;13. Если t + dt_min/2 > DT, остановиться // прилетели! PRM18 PRM19 PRM43 2 / + - ; DT - (t + dt_min/2) Px<0 KorShag PRM16 PRM17 - PM16 ; m := m - mass PGOTO Manevr ; Манёвр завершён! ;*** Коррекция шага *** ;14. Коррекция шага KorShag: ; if (rer != 0) ; then new_dt := 0.9 * ∆t / (rer^(1/k)) ; else new_dt := DT PRM18 ; DT PRM55 ; rer Px!=0 LblZZ ; Перерасчёт dt при ошибке: ; new_dt := 0.9 * ∆t / (rer^(1/k)) PPRM1001 F1/x PRM55 Fx^y ; rer ^ (1/k) 0,9 PRM54 * <-> / ; 0.9 * ∆t / .. <-> LblZZ: <-> PM54 ; ∆t ;15. Ограничить ∆t снизу (dt_min) и сверху (длительностью манёвра) PRM54 PRM43 Kmax PM54 ; ∆t := max (∆t, dt_min) PRM54 PRM18 PRM19 - Kmax <-> PM54 ; ∆t := min (∆t, DT - t) PGOTO Shag ; Перейти к пункту 3 ;16. Здесь будет блок посадки ;------------------------------------ ; БЛОК ПОСАДКИ ; Posadka: Px=0 Avost PRM39 ; p1.vy1 Px<0 Vzlet ; Кораблю - взлёт! Avost: 2 R/S ; "2": Врезались ;------------------------------------ ; БЛОК "МЕТОД ЭЙЛЕРА" ; вычисляет позицию q = R30-33 для начальной точки mp = R20-R23, mt = R56 и dt = Rb ; MetEjler: ;euler(mp, mt, dt) { ; return h(mp, dt, g(mp, mt)) ; p' := g(mp, mt) PRM20 PM30 PRM21 PM31 PRM22 PM32 PRM23 PM33 ; q := mp PRM56 PM57 ; tt := mt PGSB Fizika ; p' := g(q, tt) RMb PM59 ; mdt := dt PGOTO LinEx ; q := h(mp, mdt, p') ;} ;------------------------------------ ; БЛОК "МЕТОД РУНГЕ-КУТТА" ; вычисляет позицию R30-33 для начальной точки mp = R20-R23, mt = R56 и dt = Rb ; методом Рунге-Кутта 4-го порядка для решения системы дифференциальных уравнений ; ; Сложная схема гарантирует, что если решением задачи является многочлен 4-ой степени, ; то задача будет решена точно. Кроме того, если шаг уменьшить в два раза, и проделать ; эти операции два раза, то точность вырастет примерно в 16 раз (2^4). ; MetRK: ;rungeKutta(mp, mt, dt) { ; ; Первый прогноз соответствует методу Эйлера. ; p1' := g(mp , mt ) // начало PRM20 PM30 PRM21 PM31 PRM22 PM32 PRM23 PM33 ; q := mp PRM56 PM57 ; tt := mt PGSB Fizika ; p' := g(q, tt) PRM88 PM80 PRM89 PM81 PRM78 PM70 PRM79 PM71 ; p1' := p' ; Второй прогноз считается для средней точки на отрезке интегрирования ; (пользуемся первым прогнозом). ; p2' := g(h(mp, dt/2, p1'), mt + dt/2) // середина ; PRM80 PM88 PRM81 PM89 PRM70 PM78 PRM71 PM79 ; p' := p1' (уже записано) RMb 2 / PM59 ; mdt := dt/2 PGSB LinEx ; q := h(mp, mdt, p') PRM56 RMb 2 / + PM57 ; tt := mt + dt/2 PGSB Fizika ; p' := g(q, tt) PRM88 PM82 PRM89 PM83 PRM78 PM72 PRM79 PM73 ; p2' := p' ; Третий прогноз также для средней точки, но учитывает результаты второго прогноза. ; p3' := g(h(mp, dt/2, p2'), mt + dt/2) // ещё раз ; PRM82 PM88 PRM83 PM89 PRM72 PM78 PRM73 PM79 ; p' := p2' (уже записано) ; RMb 2 / PM59 ; mdt := dt/2 (уже записано) PGSB LinEx ; q := h(mp, mdt, p') ; PRM56 RMb 2 / + PM57 ; tt := mt + dt/2 (уже записано) PGSB Fizika ; p' := g(q, tt) PRM88 PM84 PRM89 PM85 PRM78 PM74 PRM79 PM75 ; p3' := p' ; Четвёртый прогноз считается для конца шага интегрирования. ; p4' := g(h(mp, dt , p3'), mt + dt ) // конец шага! ; PRM84 PM88 PRM85 PM89 PRM74 PM78 PRM75 PM79 ; p' := p3' (уже записано) RMb PM59 ; mdt := dt PGSB LinEx ; q := h(mp, mdt, p') PRM56 RMb + PM57 ; tt := mt + dt PGSB Fizika ; p' := g(q, tt) ; PRM88 PM86 PRM89 PM87 PRM78 PM76 PRM79 PM77 ; p4' := p' // Оптимизация: p4' не храним ; Затем считаем взвешенную сумму и "передвигаем" корабль. ; // взвешенная сумма производных ; return h(mp, dt, (p1' + 2 * (p2' + p3') + p4') / 6) ; Оптимизация: p4' считываем из p' ; p'.x' := (p1'.x' + 2 * (p2'.x' + p3'.x') + p4'.x') / 6 PRM80 PRM82 PRM84 + 2 * + PRM88 + 6 / PM88 ; p'.y' := (p1'.y' + 2 * (p2'.y' + p3'.y') + p4'.y') / 6 PRM81 PRM83 PRM85 + 2 * + PRM89 + 6 / PM89 ; p'.vx' := (p1'.vx' + 2 * (p2'.vx' + p3'.vx') + p4'.vx') / 6 PRM70 PRM72 PRM74 + 2 * + PRM78 + 6 / PM78 ; p'.vy' := (p1'.vy' + 2 * (p2'.vy' + p3'.vy') + p4'.vy') / 6 PRM71 PRM73 PRM75 + 2 * + PRM79 + 6 / PM79 ; RMb PM59 ; mdt := dt (уже записано) ;} ; PGOTO LinEx ; Идти не надо, мы уже здесь LinEx: ; Линейная экстраполяция q := h(mp, mdt, p') = mp + mdt * p' ; q (R30…33), mp (R20…23), время mdt (R59), производная p' (R88,89,78,79 = x', y', vx', vy') ; Каждый прогноз вычисляет смещение по координатам и скоростям. ;h(mp, mdt, p’) { ; вычисляем новые полярные координаты q.x, q.y PRM20 PRM88 PRM59 * + PM30 ; q.x = mp.x + mdt * p’.x’ PRM21 PRM89 PRM59 * + PM31 ; q.y = mp.y + mdt * p’.y’ ; вычисляем новую скорость q.vx, q.vy PRM22 PRM78 PRM59 * + PM32 ; q.vx = mp.vx + mdt * p’.vx' PRM23 PRM79 PRM59 * + PM33 ; q.vy = mp.vy + mdt * p’.vy' ; return q ;} RTN ;------------------------------------ ; БЛОК "МЕТОД ДОРМАНД-ПРИНСА" ; ; Метод Дорманд-Принса для начальной точки mp = R20-R23, mt = R56 и dt = Rb ; возвращает два вектора p1 = R36-39 - решение, p2 = R50-53 нужен только для расчёта ошибки: ; DorPrinc: ;dormandPrince(mp, mt, dt) { ; p1' := g( mp , mt ) PRM20 PM30 PRM21 PM31 PRM22 PM32 PRM23 PM33 ; q := mp PRM56 PM57 ; tt := mt PGSB Fizika ; p' := g(q, tt) PRM88 PM80 PRM89 PM81 PRM78 PM70 PRM79 PM71 ; p1' := p' ; p2' := g(h2(mp, dt, c21, p1' ), mt + c2 * dt) ; q := h2(mp, dt, c21, p1') = mp + dt * (c21 * p1') ; Оптимизация: c21 = 1/5 PRM80 5 / RMb * PRM20 + PM30 ; q.x = mp.x + dt * c21 * p1’.x’ PRM81 5 / RMb * PRM21 + PM31 ; q.y = mp.y + dt * c21 * p1’.y’ PRM70 5 / RMb * PRM22 + PM32 ; q.vx = mp.vx + dt * c21 * p1’.vx' PRM71 5 / RMb * PRM23 + PM33 ; q.vy = mp.vy + dt * c21 * p1’.vy' PRM56 RMb ; mt, dt 0,2 * ; c2 = 1/5 + PM57 ; tt := mt + c2 * dt PGSB Fizika ; p' := g(q, tt) PRM88 PM82 PRM89 PM83 PRM78 PM72 PRM79 PM73 ; p2' := p' ; p3' := g(h3(mp, dt, c31, p1', c32, p2' ), mt + c3 * dt) ; q := h3(mp, dt, c31, p1', c32, p2') = mp + dt * (c31 * p1' + c32 * p2') PRM80 PPRM201 * PRM82 PPRM202 * + RMb * PRM20 + PM30 ; q.x = mp.x + dt*( c31 * p1’.x’ + c32 * p2'.x' ) PRM81 PPRM201 * PRM83 PPRM202 * + RMb * PRM21 + PM31 ; q.y = mp.y + dt*( c31 * p1’.y’ + c32 * p2'.y' ) PRM70 PPRM201 * PRM72 PPRM202 * + RMb * PRM22 + PM32 ; q.vx = mp.vx + dt*( c31 * p1’.vx' + c32 * p2'.vx') PRM71 PPRM201 * PRM73 PPRM202 * + RMb * PRM23 + PM33 ; q.vy = mp.vy + dt*( c31 * p1’.vy' + c32 * p2'.vy') PRM56 RMb ; mt, dt 0,3 * ; c3 = 3/10 + PM57 ; tt := mt + c3 * dt PGSB Fizika ; p' := g(q, tt) PRM88 PM84 PRM89 PM85 PRM78 PM74 PRM79 PM75 ; p3' := p' ; p4' := g(h4(mp, dt, c41, p1', c42, p2', c43, p3' ), mt + c4 * dt) ; q := h4(mp, dt, c41, p1', c42, p2', c43, p3') = mp + dt * (c41 * p1' + c42 * p2' + c43 * p3') PRM80 PPRM203 * PRM82 PPRM204 * + PRM84 PPRM205 * + RMb * PRM20 + PM30 ; …c41 * p1’.x’ + c42 * p2'.x' + c43 * p3'.x' ) PRM81 PPRM203 * PRM83 PPRM204 * + PRM85 PPRM205 * + RMb * PRM21 + PM31 ; …c41 * p1’.y’ + c42 * p2'.y' + c43 * p3'.y' ) PRM70 PPRM203 * PRM72 PPRM204 * + PRM74 PPRM205 * + RMb * PRM22 + PM32 ; …c41 * p1’.vx'+ c42 * p2'.vx'+ c43 * p3'.vx') PRM71 PPRM203 * PRM73 PPRM204 * + PRM75 PPRM205 * + RMb * PRM23 + PM33 ; …c41 * p1’.vy'+ c42 * p2'.vy'+ c43 * p3'.vy') PRM56 RMb ; mt, dt 0,8 * ; c4 = 4/5 + PM57 ; tt := mt + c4 * dt PGSB Fizika ; p' := g(q, tt) PRM88 PM86 PRM89 PM87 PRM78 PM76 PRM79 PM77 ; p4' := p' ; p5' := g(h5(mp, dt, c51, p1', c52, p2', c53, p3', c54, p4' ), mt + c5 * dt) ; q := h5(mp, dt, c51, p1', c52, p2', c53, p3', c54, p4') = ; = mp + dt * (c51 * p1' + c52 * p2' + c53 * p3' + c54 * p4') PRM80 PPRM206 * PRM82 PPRM207 * + PRM84 PPRM208 * + PRM86 PPRM209 * + RMb * PRM20 + PM30 ; … + c54 * p4'.x' ) PRM81 PPRM206 * PRM83 PPRM207 * + PRM85 PPRM208 * + PRM87 PPRM209 * + RMb * PRM21 + PM31 ; … + c54 * p4'.y' ) PRM70 PPRM206 * PRM72 PPRM207 * + PRM74 PPRM208 * + PRM76 PPRM209 * + RMb * PRM22 + PM32 ; … + c54 * p4'.vx') PRM71 PPRM206 * PRM73 PPRM207 * + PRM75 PPRM208 * + PRM77 PPRM209 * + RMb * PRM23 + PM33 ; … + c54 * p4'.vy') PRM56 RMb ; mt, dt 8 * 9 / ; c5 = 8/9 + PM57 ; tt := mt + c5 * dt PGSB Fizika ; p' := g(q, tt) PRM88 PM90 PRM89 PM91 PRM78 PM60 PRM79 PM61 ; p5' := p' ; p6' := g(h6(mp, dt, c61, p1', c62, p2', c63, p3', c64, p4', c65, p5' ), mt + c6 * dt) ; q := h6(mp, dt, c61, p1', c62, p2', c63, p3', c64, p4', c65, p5') = ; = mp + dt * (c61 * p1' + c62 * p2' + c63 * p3' + c64 * p4' + c65 * p5') PRM80 PPRM210 * PRM82 PPRM211 * + PRM84 PPRM212 * + PRM86 PPRM213 * + PRM90 PPRM214 * + RMb * PRM20 + PM30 ;…+c65*p5'.x' ) PRM81 PPRM210 * PRM83 PPRM211 * + PRM85 PPRM212 * + PRM87 PPRM213 * + PRM91 PPRM214 * + RMb * PRM21 + PM31 ;…+c65*p5'.y' ) PRM70 PPRM210 * PRM72 PPRM211 * + PRM74 PPRM212 * + PRM76 PPRM213 * + PRM60 PPRM214 * + RMb * PRM22 + PM32 ;…+c65*p5'.vx') PRM71 PPRM210 * PRM73 PPRM211 * + PRM75 PPRM212 * + PRM77 PPRM213 * + PRM61 PPRM214 * + RMb * PRM23 + PM33 ;…+c65*p5'.vy') PRM56 RMb ; mt, dt ; c6 = 1 + PM57 ; tt := mt + c6 * dt PGSB Fizika ; p' := g(q, tt) PRM88 PM92 PRM89 PM93 PRM78 PM62 PRM79 PM63 ; p6' := p' ; p7' := g(h7(mp, dt, c71, p1', c72, p2', c73, p3', c74, p4', c75, p5', c76, p6'), mt + c7 * dt) = ; q := h7(mp, dt, c71, p1', c72, p2', c73, p3', c74, p4', c75, p5', c76, p6') = ; = mp + dt * (c71 * p1' + c72 * p2' + c73 * p3' + c74 * p4' + c75 * p5' + c76 * p6') ; Оптимизация: R216 = c72 = 0 ; Оптимизация: p2 := q PRM80 PPRM215 * PRM84 PPRM217 * + PRM86 PPRM218 * + PRM90 PPRM219 * + PRM92 PPRM220 * + RMb * PRM20 + PM30 PM50 ; p2.x := q.x := mp.y + dt * … PRM81 PPRM215 * PRM85 PPRM217 * + PRM87 PPRM218 * + PRM91 PPRM219 * + PRM93 PPRM220 * + RMb * PRM21 + PM31 PM51 ; p2.y := q.y := mp.y + dt * … PRM70 PPRM215 * PRM74 PPRM217 * + PRM76 PPRM218 * + PRM60 PPRM219 * + PRM62 PPRM220 * + RMb * PRM22 + PM32 PM52 ; p2.vx := q.vx := mp.vx + dt * … PRM71 PPRM215 * PRM75 PPRM217 * + PRM77 PPRM218 * + PRM61 PPRM219 * + PRM63 PPRM220 * + RMb * PRM23 + PM33 PM53 ; p2.vy := q.vy := mp.vy + dt * … PRM56 RMb ; mt, dt ; c7 = 1 + PM57 ; tt := mt + c7 * dt PGSB Fizika ; p' := g(q, tt) ; PRM88 PM94 PRM89 PM95 PRM78 PM64 PRM79 PM65 ; p7' := p' /* Оптимизация: p7' не храним */ ; //считаем два решения ; p1 := h8(mp, dt, c81, p1', c82, p2', c83, p3', c84, p4', c85, p5', c86, p6', c87, p7') = ; = mp + dt * (c81 * p1' + c82 * p2' + c83 * p3' + c84 * p4' + c85 * p5' + c86 * p6' + c87 * p7') ; Оптимизация: p7' считываем из p' ; Оптимизация: R222 = c82 = 0 ; Оптимизация: c87 = 1/40 PRM80 PPRM221 * PRM84 PPRM223 * + PRM86 PPRM224 * + PRM90 PPRM225 * + PRM92 PPRM226 * + PRM88 40 / + RMb * PRM20 + PM36 ; p1.x := mp.x + dt * … PRM81 PPRM221 * PRM85 PPRM223 * + PRM87 PPRM224 * + PRM91 PPRM225 * + PRM93 PPRM226 * + PRM89 40 / + RMb * PRM21 + PM37 ; p1.y := mp.y + dt * … PRM70 PPRM221 * PRM74 PPRM223 * + PRM76 PPRM224 * + PRM60 PPRM225 * + PRM62 PPRM226 * + PRM78 40 / + RMb * PRM22 + PM38 ; p1.vx := mp.vx + dt * … PRM71 PPRM221 * PRM75 PPRM223 * + PRM77 PPRM224 * + PRM61 PPRM225 * + PRM63 PPRM226 * + PRM79 40 / + RMb * PRM23 + PM39 ; p1.vy := mp.vy + dt * … ; p2 := h8(mp, dt, d81, p1', d82, p2', d83, p3', d84, p4', d85, p5', d86, p6', d87, p7') ; мы уже рассчитали p2 заранее, при вычислении p7' /* оптимизация */ ; return (p1, p2); ;} RTN ;------------------------------------ ; РАСЧЁТНЫЙ БЛОК "ФИЗИКА" ; вычисляет производную p' = g(q, tt) в точке q = R30-R33, в момент времени tt = R57 ; Fizika: ; g(q, tt) { ;1. вычисляем модуль реактивного ускорения ra для постоянной тяги mass/DT ; // полная масса корабля = сухая + топливо - то что израсходовано ; fm(tt) := M + m - (tt / DT) * mass PRM40 PRM16 + PRM57 PRM18 / PRM17 * - ; fm оставим в стеке ; мгновенное реактивное ускорение PRM18 * PRM41 PRM17 * <-> / Mc ; ra(tt) := C * mass / (DT * fm(tt)) ;2. вычисляем проекции реактивного ускорения rax, ray ; // проекция реактивного ускорения на местную горизонталь (трансверсальное ускорение) RMc RMd Fsin * PM24 ; rax(tt) := ra(tt) * sin(da) ; // проекция реактивного ускорения на местную вертикаль (радиальное ускорение) RMc RMd Fcos * PM25 ; ray(tt) := ra(tt) * cos(da) ;3. вычисляем производную p' PRM32 PRM31 / PM88 ; p'.x' := q.vx / q.y PRM33 PM89 ; p'.y' := q.vy PRM24 PRM33 PRM32 * PRM31 / - PM78 ; p'.vx' := rax(tt) - q.vy * q.vx / q.y ; p'.vy' := ray(tt) + q.vx * q.vx / q.y - gc / q.y^2 PRM25 PRM32 FX^2 PRM31 / + RMa PRM31 FX^2 / - PM79 ; return p' ;} RTN ;------------------------------------ ; РАСЧЁТНЫЙ БЛОК "ОРБИТА" ; Вычисляет параметры орбиты корабля ; ; Ra,10 gc - гравитационная постоянная небесного тела ; ;// положение в полярной системе координат (координаты и скорости) ;p = (x, y, vx, vy) ; ; R26 p0.x — вектор координаты корабля "МАТЕМАТИКА" ; R27 p0.y ; R28 p0.vx — вектор скорости корабля "МАТЕМАТИКА" ; R29 p0.vy ; ;// кеплерова орбита в полярной системе координат ;// (некий важный параметр, ;// эксцентриситет, ;// угловая координата перицентра) ;o = (p, e, xperi) ; ; ; R0300 o.p — некий важный параметр (кеплерова орбита в полярной системе координат) ; R0301 ecos ; R0302 esin ; R0303 o.e — эксцентриситет (кеплерова орбита в полярной системе координат) ; R0304 o.xperi — угловая координата перицентра (кеплерова орбита в полярной системе координат) ; R0305 apo — апоцентр ; R0306 theta — угловая кооррдината на кеплеровой орбите, измеряется от перицентра ; ;// вычисляем кеплерову орбиту и угол по полярным координатам ;p2o (p) { ; o.p = (vx * y)^2 / gc ; ecos = (vx * vx * y) / gc - 1 ; esin = (vx * vy * y) / gc ; o.e = sqrt(ecos^2 + esin^2) ; if (o.e > 0) ; theta = arg(ecos, esin) ; else ; theta = 0 ; o.xperi = x - theta ; return o, theta ;} ; P2O: PRM28 PRM27 * Fx^2 RMa / PPM0300 ; o.p := (p0.vx * p0.y)^2 / gc PRM28 Fx^2 PRM27 * RMa / 1 - PPM0301 ; ecos := (p0.vx * p0.vx * p0.y) / gc - 1 PRM28 PRM29 * PRM27 * RMa / PPM0302 ; esin := (p0.vx * p0.vy * p0.y) / gc PPRM0301 Fx^2 PPRM0302 Fx^2 + FSQRT PPM0303 ; o.e := sqrt(ecos^2 + esin^2) PX!=0 P2Ozero PPRM0301 PPRM0302 PGSB Arg ; RX := arg(ecos, esin) P2Ozero: PPM0306 ; theta := RX PRM26 PPRM0306 - PPM0304 ; o.xperi := p0.x - theta RTN ;// математика полностью не работает для вертикального полёта (vx == 0) ; ;// arg(x, y) - полярный аргумент точки (x, y). Например arg(x, y) = atan2(y, x). ; ; R0307 x — аргумент функции Arg ; R0308 y — аргумент функции Arg ; R0309 sx — знак x, функция Arg ; R0310 sy — знак y, фукнкция Arg ; ;arg(x, y) ; ; (x y -- arg) ; Arg: PPM0308 <-> PPM0307 ; Сохранить x, y PPRM0307 KSGN PPM0309 ; sx := sign(x) PPRM0307 KABS PPM0307 ; x := abs(x) PPRM0308 KSGN PPM0310 ; sy := sign(y) PPRM0308 KABS PPM0308 ; y := abs(y) ; if (y <= x) ; arg = arctg(y / x) ; else ; arg = pi / 2 - arctg(x / y) PPRM0307 PPRM0308 - ; RX := x - y Px>=0 ArgElse PPRM0308 PPRM0307 / Farctg ; RX := arctg(y / x) PGOTO ArgEnd ArgElse: FPI 2 / PPRM0307 PPRM0308 / Farctg - ; RX := pi / 2 - arctg(x / y) ArgEnd: ; if (sx < 0) ; arg = pi - arg PPRM0309 ; sx Px<0 ArgSX <-> FPI <-> - <-> ArgSX: <-> ; if (sy < 0) ; arg = -arg PPRM0310 ; sy Px<0 ArgSY <-> +/- <-> ArgSY: ; return arg <-> RTN ; ;// tgh - гиперболический тангенс ;// sinh - гиперболический синус ;// arctgh - гиперболический арктангенс ;// ;// tgh(x) = (exp(2 * x) - 1) / (exp(2 * x) + 1) Tgh: 2 * FEXP 1 <-> + FANS 1 - <-> / RTN ;// sinh(x) = (exp(2 * x) - 1) / (2 * exp(x)) Sinh: ENT FEXP 2 * <-> 2 * FEXP 1 - <-> / RTN ;// arctgh(x) = log((1 + x) / (1 - x)) / 2 ; ; Делаю в предположении, что log это натуральный логарифм ; Arctgh: 1 <-> + FANS 1 <-> - / Fln 2 / RTN ; ; R0311 l — Root ; R0312 r — Root ; R0313 M — Root, Root1, Root2 ; R0314 m — Root ; R0315 M1 — Root1 ; R0316 E — equation() in Root, Root2; E2theta, E2t ; R9 equation — адрес решаемого уравнения ; ; ;root(l, r, equation, M) { ; while (true) { ; m = (l + r) / 2 ; if (m <= l || r <= m) ; return m ; if (equation(m) <= M) ; l = m ; else ; r = m ; } ;} Root: PPRM0311 PPRM0312 + 2 / PPM0314 ; m := (l + r) / 2 PPRM0311 PPRM0314 - ; RX := l - m Px<0 RootRTN PPRM0314 PPRM0312 - ; RX := m - r Px>=0 RootQQ RootRTN: PPRM0314 RTN RootQQ: PPRM0314 KGSB9 ; RX := equation(m) PPRM0313 <-> - ; RX := M - equation(m) Px>=0 RootR PPRM0314 PPM0311 ; l := m PGOTO Root RootR: PPRM0314 PPM0312 ; r := m PGOTO Root ; ;root1(equation, M) { ; M1 = 2 * pi * floor(M / (2 * pi)) ; return M1 + root(0, 2 * pi, equation, M - M1) ;} ; Root1: PPRM0313 2 FPI * / KINT 2 * FPI * PPM0315 ; M1 := 2 * pi * [ M / (2 * pi) ] PPRM0313 PPRM0315 - PPM0313 ; M := M - M1 Cx PPM0311 ; l := 0 2 FPI * PPM0312 ; r := 2 * pi PGSB Root PPRM0315 + RTN ; return M1 + root(…) ; ; R0311 l — Root ; R0312 r — Root ; R0313 M — Root, Root1, Root2 ; R0314 m — Root ; R0315 M1 — Root1 ; ;ooot2(equation, M) { ; if (M < 0) ; r = 0 ; l = -1 ; while (equation(l) >= M) ; r = l ; l = l * 2 ; else ; l = 0 ; r = 1 ; while (equation(r) < M) ; l = r ; r = r * 2 ; ; return root(l, r, equation, M) ;} ; Root2: PPRM0313 Px<0 Root2p Cx PPM0312 ; r := 0 1 +/- PPM0311 ; l := -1 Root2w1: PPRM0311 KGSB9 ; RX := equation(l) PPRM0313 - ; RX := equation(l) - M Px>=0 Root PPRM0311 PPM0312 2 * PPM0311 ; r := l; l := l * 2 PGOTO Root2w1 Root2p: Cx PPM0311 ; l := 0 1 PPM0312 ; r := 1 Root2w2: PPRM0312 KGSB9 ; RX := equation(r) PPRM0313 - ; RX := equation(r) - M Px<0 Root PPRM0312 PPM0311 2 * PPM0312 ; l := r; r := r * 2 PGOTO Root2w2 ; ; ; R26 p0.x — вектор координаты корабля "МАТЕМАТИКА" ; R27 p0.y ; R28 p0.vx — вектор скорости корабля "МАТЕМАТИКА" ; R29 p0.vy ; ; R0300 o.p — некий важный параметр (кеплерова орбита в полярной системе координат) ; R0301 ecos ; R0302 esin ; R0303 o.e — эксцентриситет (кеплерова орбита в полярной системе координат) ; R0304 o.xperi — угловая координата перицентра (кеплерова орбита в полярной системе координат) ; R0305 apo — апоцентр ; R0306 theta — угловая кооррдината на кеплеровой орбите, измеряется от перицентра ; ; R0325 theta1, E1 — по аналогии с M1 ; ;// вычисляем полярные кординаты по кеплеровой орбите и углу на ней ;o2p (o, theta) O2P: PPRM0304 PPRM0306 + PM26 ; p0.x := o.xperi + theta PPRM0306 Fcos PPRM0303 * 1 + PPRM0300 <-> / PM27 ; p0.y := o.p / (1 + o.e * cos(theta)) RMa PPRM0300 / FSQRT ; RX := sqrt(gc / o.p) PPRM0306 Fcos PPRM0303 * 1 + <-> * PM28 ; p0.vx := sqrt(gc / o.p) * (1 + o.e * cos(theta)) ; используется RX1 !!! FANS PPRM0306 Fsin PPRM0303 * * PM29 ; p0.vy := sqrt(gc / o.p) * o.e * sin(theta) RTN ; ;// переводим угловой параметр theta в E (ещё один угол, по которому будем считать время) ;theta2E(e, theta) Theta2E: PPRM0303 1 - Px=0 Th2En1 PPRM0306 2 / Ftg RTN ; E := tg(theta / 2) /* o.e == 1 */ Th2En1: Px<0 ThEb1 ; if (o.e < 1) ; E := 2 * arg(sqrt(1 + o.e) * cos(theta / 2), sqrt(1 - o.e) * sin(theta / 2)) PPRM0306 2 FPI * / KINT 2 * FPI * PPM0325 ; theta1 := 2 * pi * [ theta / (2 * pi) ] PPRM0303 1 + FSQRT ; sqrt(1 + o.e) PPRM0306 PPRM0325 - 2 / Fcos * ; * cos((theta - theta1) / 2), 1 PPRM0303 - FSQRT ; sqrt(1 - o.e) PPRM0306 PPRM0325 - 2 / Fsin * ; * sin((theta - theta1) / 2) PGSB Arg 2 * PPRM0325 + ; E += theta1 RTN ThEb1: ; if (o.e > 1) ; E := 2 * arctgh(sqrt((o.e - 1) / (o.e + 1)) * tg(theta / 2)) PPRM0303 1 - PPRM0303 1 + / FSQRT PPRM0306 2 / Ftg * PGSB Arctgh 2 * RTN ; R0303 o.e — эксцентриситет (кеплерова орбита в полярной системе координат) ; R0306 theta — угловая кооррдината на кеплеровой орбите, измеряется от перицентра ; R0316 E — equation() in Root, Root2; E2theta, E2t ; ;// переводим угловой параметр E в theta ;E2theta(e, E) E2theta: PPM0316 ; E := RX PPRM0303 1 - Px=0 E2thN1 PPRM0316 Farctg 2 * RTN ; theta = 2 * arctg(E) /* o.e == 1 */ E2thN1: Px<0 E2thB1 ; if (o.e < 1) ; theta := 2 * arg(sqrt(1 - o.e) * cos(E / 2), sqrt(1 + o.e) * sin(E / 2)) PPRM0316 2 FPI * / KINT 2 * FPI * PPM0325 ; E1 := 2 * pi * [ E / (2 * pi) ] 1 PPRM0303 - FSQRT ; sqrt(1 - o.e) PPRM0316 PPRM0325 - 2 / Fcos * ; * cos((E - E1) / 2), 1 PPRM0303 + FSQRT ; sqrt(1 + o.e) PPRM0316 PPRM0325 - 2 / Fsin * ; * sin((E - E1) / 2) PGSB Arg 2 * PPRM0325 + ; theta += E1 RTN E2thB1: ; if (o.e > 1) ; theta := 2 * arctg(sqrt((o.e + 1) / (o.e - 1)) * tgh(E / 2)) PPRM0316 2 / PGSB Tgh PPRM0303 1 + PPRM0303 1 - / FSQRT * Farctg 2 * RTN ; ; ;// перевод E во время полёта от перицентра (E = 0) до указанного параметра (E) ;E2t(o, E) ; E2t: PPM0316 ; E := RX PPRM0303 1 - Px=0 E2tN1 PPRM0300 ENT RMa / Fsqrt * 2 / PPM0317 ; C := o.p * sqrt(o.p / gc) / 2 /* o.e == 1 */ PPRM0316 PGSB Eq2 ; M := E + E^3 / 3 PGOTO E2tMulRTN E2tN1: Px<0 E2tB1 ; if (o.e < 1) PPRM0300 ENT ENT 1 PPRM0303 Fx^2 - / ; a := o.p / (1 - o.e^2), экономим регистр ENT RMa / FSQRT * PPM0317 ; C := a * sqrt(a / gc), используем пред. расчёт PPRM0316 PGSB Eq1 ; M := E - o.e * sin(E) PGOTO E2tMulRTN E2tB1: ; if (o.e > 1) PPRM0300 PPRM0303 Fx^2 1 - / ; a := o.p / (o.e^2 - 1), экономим регистр ENT RMa / FSQRT * PPM0317 ; C := a * sqrt(a / gc), используем пред. расчёт PPRM0316 PGSB Eq3 ; M := o.e * sinh(E) - E E2tMulRTN: PPRM0317 * RTN ; return C * M Eq1: ENT Fsin PPRM0303 * - RTN ; "E - o.e * sin(E)" Eq2: ENT ENT Fx^2 * 3 / + RTN ; "E + E^3 / 3" Eq3: PPM0316 PGSB Sinh PPRM0303 * PPRM0316 - RTN ; "o.e * sinh(E) - E" ;// перевод времени полёта от перицентра (E = 0) в параметр E ;t2E(o, t) ; ; R0300 o.p — некий важный параметр (кеплерова орбита в полярной системе координат) ; R0303 o.e — эксцентриситет (кеплерова орбита в полярной системе координат) ; R0304 o.xperi — угловая координата перицентра (кеплерова орбита в полярной системе координат) ; R0313 M — Root, Root1, Root2 ; R0316 E — equation() в Root, Root2; E2theta, E2t ; R0317 C — E2t ; R0318 t — T2E ; R0319 E — FreeFlight T2E: PPM0318 ; t := RX PPRM0303 1 - Px=0 T2En1 .NUM Eq2 M9 ; equation := Eq2 PPRM0300 ENT RMa / FSQRT * 2 / PPM0317 ; C := o.p * sqrt(o.p / gc) / 2 PGOTO T2Eqq T2En1: Px<0 T2Eb1 .NUM Eq1 M9 ; equation := Eq1 PPRM0300 1 PPRM0303 Fx^2 - / ; RX := o.p / (1 - o.e^2) ENT RMa / FSQRT * PPM0317 ; C := RX * sqrt(RX / gc) PPRM0318 PPRM0317 / PPM0313 ; M := t / C PGOTO Root1 T2Eb1: .NUM Eq3 M9 ; equation := Eq3 PPRM0300 PPRM0303 Fx^2 1 - / ; RX := o.p / (o.e^2 - 1) ENT RMa / FSQRT * PPM0317 ; C := RX * sqrt(RX / gc) T2Eqq: PPRM0318 PPRM0317 / PPM0313 ; M := t / C PGOTO Root2 ;---------------------------------------------- ; ; R42 R — радиус планеты ; R0300 o.p — некий важный параметр (кеплерова орбита в полярной системе координат) ; R0303 o.e — эксцентриситет (кеплерова орбита в полярной системе координат) ; R0306 theta — угловая кооррдината на кеплеровой орбите, измеряется от перицентра ; R0320 E — LandingQ ; R0321 E1 — LandingQ ; R0322 peri — LandingQ ; ;// определение момента и точки посадки по теущему положению ;landing?(p) LandingQ: PGSB TajmerStart 1 PPM9045 ; Радианы PGSB P2O ; o, theta := p2o(p) PPRM0300 1 PPRM0303 + / PPM0322 ; peri := o.p / (1 + o.e) ; if peri >= R // перицентр над поверхностью, посадки не будет ; return "no_landing" PPRM0322 PRM42 - Px>=0 Land1 NoLand: .NUM msgNoLanding PPM9026 ; Вывести приглашение R/S PGOTO Manevr Land1: PGSB Theta2E PPM0320 ; E := Theta2E(o.e, theta) ; // угловая координата момента посадки относительно перицентра ; // посадка всегда происходит до перицентра (в прошлом) ; theta = -arccos((o.p / R - 1) / o.e) PPRM0300 PRM42 / 1 - PPRM0303 / Farccos +/- PPM0306 PGSB Theta2E PPM0321 ; E1 := Theta2E(o.e, theta) ; // текущий момент и момент посадки могут быть переставлены во времени ; if (E1 < E) ; if (o.e < 1) // корректируем "угол" посадки для эллипса и круга ; E1 = E1 + 2 * pi ; else // мы на параболе или гиперболе и поднимаемся ; return "no_landing" PPRM0321 PPRM0320 - Px<0 Land2 PPRM0303 1 - Px<0 NoLand PPRM0321 FPI 2 * + PPM0321 Land2: ; // время до посадки ; dt = E2t(o, E1) - E2t(o, E) PPRM0321 PGSB E2t PM18 PPRM0320 PGSB E2t PRM18 <-> - PM18 PGSB O2P ; p1 = o2p(o, theta) PGOTO Manevr ; return dt, p1 ;------------------------------------------ ;// полёт без реактивного ускорения ;free_flight(p, DT) FreeFlight: PM18 ; DT := RX PGSB TajmerStart ; Засечь время работы ЭКВМ 1 PPM9045 ; Радианы PGSB P2O ; o, theta := p2o(p) ; ; // текущая координата PGSB Theta2E PPM319 ; E := Theta2E(o.e, theta) ; ; // "углы" через dt секунд ; E1 = t2E(o, E2t(o, E) + dt) PGSB E2t PRM18 + PGSB T2E ; theta = E2theta(o.e, E1) PGSB E2theta PPM0306 PGSB O2P ; return o2p(o, theta) PGOTO Manevr ;------------------------------------ ; БЛОК ТАЙМЕР ; Измеряет время работы ЭКВМ ; Записывает результат в R99 ; TajmerStart: 65535 PPM 9051 RTN TajmerStop: 65535 PPRM 9051 - EE 2 +/- PM99 RTN ;------------------------------------ ; БЛОК ДАННЫХ ; Текстовые сообщения ; msgUgol: .DB 93H,0a3H,0aeH,0abH,20H ;Угол .DB 0e2H,0efH,0a3H,0a8H,2cH,20H ;тяги, .DB 0a3H,0e0H,0a0H,0a4H,0e3H,0e1H,0ebH,20H,3fH,0 ;градусы ? msgMass: .DB 090H,0a0H,0e1H,0e5H,0aeH,0a4H,20H ;Расход .DB 0e2H,0aeH,0afH,0abH,0a8H,0a2H,0a0H,2cH,20H ;топлива, .DB 0aaH,0a3H,20H,3fH,0 ;кг ? msgDT: .DB 82H,0e0H,0a5H,0acH,0efH,2cH,20H,0e1H,20H,3fH,0 ;Время, с ? msgNoLanding: .TEXT "No Landing!" .DB 0 ; ; Коэффициенты метода Дорманд-Принс ; ; 0 | ; c2 | c21 ; c3 | c31 c32 ; c4 | c41 c42 c43 ; c5 | c51 c52 c53 c54 ; c6 | c61 c62 c63 c64 c65 ; c7 | c71 c72 c73 c74 c75 c76 ; ----------------------------------- ; | c81 c82 c83 c84 c85 c86 c87 ; | d81 d82 d83 d84 d85 d86 d87 ; ; 0 | ; 1/5 | 1/5 ; 3/10 | 3/40 9/40 ; 4/5 | 44/45 −56/15 32/9 ; 8/ 9 | 19372/6561 −25360/2187 64448/6561 −212/729 ; 1 | 9017/3168 −355/33 46732/5247 49/176 −5103/18656 ; 1 | 35/384 0 500/1113 125/192 −2187/6784 11/84 ; -------------------------------------------------------------------------------------- ; | 5179/57600 0 7571/16695 393/640 −92097/339200 187/2100 1/40 ; | 35/ 384 0 500/1113 125/192 −2187/6784 11/84 0 ; .END