Тред №135884

Исходный текст:

#!/usr/bin/perl -w
use Math::Trig ":pi";

use strict;
use vars qw($Arx $Ary $Akx $Aky $Ar $Ak); # компоненты радиус-вектора и касательного вектора и их модули
use vars qw($Kx $Ky); # координаты корабля
use vars qw($Vx $Vy $V); # компонеты вектора скорости и его модуль
use vars qw($Gx $Gy $G); # компоненты вектора силы тяжестии его модули
use vars qw($FDA $FDK $FDAx $FDKx $FDAy $FDKy $FDx $FDy $FD $FDt1 $FDt2); # компоненты ветора тяги двигателя и его модуль в локальной
                                            # и общей координатах и временной значение для расчетов.
use vars qw($Fy $Fx $F); # сумма всех сил и ее компонеты.
use vars qw($Ay $Ax $A); # ускорение и его компонеты.
use vars qw($mass0 $mass $dmass $rash_topl_sek); # начальная масса, текущая масса,
                                                # изменение массы за цикл, секундный расход топлива
use vars qw($dT $Td $Tr $T); # шаг по времени, время работы двигателя, время расчета, время от начала
use vars qw($R $H $Hkm); # Расстояние от центра Земли и высота.
use vars qw($V0x $V0y);
use vars qw($Vras); # Идеальная расчетная первая космическая скрость для высоты
use vars qw($Vmin $Vmax $Hmin $Hmax); # Максимальные и минимальные скорости и высоты
use vars qw($Urys $Utang); # Углы тангажа и рыскания

use constant GCONST => 6.67259E-11; # Гравитационная постоянная
use constant MASSZ => 5.97371E24; # Масса Земли
use constant RZ => 6366250; # Радиус Земли

# Начальные значения

$H = 180000;
$Kx = RZ + $H;
$Ky = 0;
$mass0 = 124000;
$mass = 124000;
$Tr = 36000;
$dT = 0.001;
$Td = 305;
$T = 0;
$rash_topl_sek = 210.88;
$dmass = $rash_topl_sek * $dT;
$FDx = 0;
$FDy = 0;
$FD = 0;
$Urys = 16;
$Utang = 67;

$mass = $mass0 - ($dmass / 2); # Средняя масса за период расчета одной итерации

sub VECT # Вычисление опорных весторов и силы тяжести (средней)
{
   $R = sqrt($Kx**2 + $Ky**2); # Вычисление расстояния от центра земли и высоты от ее поверхности.
   $H = $R - RZ;

   $Arx = $Kx/$R; # Вычисление радиус-вектора, единичного
   $Ary = $Ky/$R;
   $Ar = sqrt($Arx**2 + $Ary**2);

   $Akx = $Ary; # Вычисление касательного вектора в направлении полета, единичного
   $Aky = -$Arx;
   $Ak = sqrt($Akx**2 + $Aky**2);

   $G = (GCONST * $mass * MASSZ) / ($R**2); # Вычисление вектора силы тяжести
   $Gx = -($Arx * $G);
   $Gy = -($Ary * $G);
}

$Vras = sqrt (GCONST * MASSZ / (RZ + $H)); # Расчет 1-й космической для данной высоты

print "Расчетная V1 = $Vras\n";

$Vy = -$Vras; # Задаем начальную скорость
$Vx = 0;
$V = $Vras;

&VECT;

print "R = $R, H = $H\n";
print "Ar = $Ar, Arx = $Arx, Ary = $Ary\n";
print "Ak = $Ak, Akx = $Akx, Aky = $Aky\n";
print "G = $G, Gx = $Gx, Gy = $Gy\n";
print "успешно\n";

use vars qw($i);

$FD = 891351 * cos(pi2/360*$Urys); # величина проекции тяги двигателя на плоскость орбиты.
$FDA = $FD * cos(pi2/360*$Utang); # Составляющая тяги на радиус-вектор, пересчет их градусов в радианы
$FDK = $FD * cos(pi2/360*(90-$Utang)); # Составляющая тяги на касательный вектор, пересчет их градусов в радианы
$FDt1 = sqrt($FDA**2 + $FDK**2); # Контроль совпадения
$FDx = $FDA * $Arx + $FDK * $Akx; # Пересчет составляющей тяги на ось X
$FDy = $FDA * $Ary + $FDK * $Aky; # Пересчет составляющей тяги на ось Y
$FDt2 = sqrt($FDx**2 + $FDy**2); # Контроль совпадения

print "FD = $FD, FDt1 = $FDt1, FDx = $FDx, FDy = $FDy, FDt2 = $FDt2\n";

$i = 0; # Итератор

$Vmin = $V; # Начальные для минимальных и максимальных скоростей и высот
$Vmax = $V;
$Hmin = $H;
$Hmax = $H;

while ($T < $Tr) # начало главного цикла
{
$T = $dT*$i; # Текущее время

if ($T < $Td) # Если текущее время меньше времени работы двигателя, учитываем
             # измение массы за счет расхода топлива и берем среднюю за время итерации.
{
$mass = $mass0 - ($dmass * $i) - ($dmass / 2);
}

$Fx = $Gx; # Сила, действующая на корабль, равна силе гравитации
$Fy = $Gy;

if ($T < $Td) # Если двигатель работатет, добавляем силу тяги двигателя.
{
   $Fx = $Gx + $FDA * $Arx + $FDK * $Akx;
   $Fy = $Gy + $FDA * $Ary + $FDK * $Aky;
}

$Ax = $Fx / $mass; # Ускорение по X
$Ay = $Fy / $mass; # Ускорение по Y
$V0x = $Vx; # Начальная скорость равна вычисленной в пред. цикле
$V0y = $Vy;
$Vx = $V0x + $Ax * $dT; #Скорость в конце цикла
$Vy = $V0y + $Ay * $dT;
$Kx = $Kx + ($V0x + $Vx) * $dT / 2; # Координаты в конце цикла
$Ky = $Ky + ($V0y + $Vy) * $dT / 2;

$mass = $mass0 - ($dmass * ($i+1)) - ($dmass / 2); # Средняя масса для следующего цикла
&VECT; # Вычисление опорных весторов и силы тяжести (средней) для конца цикла

if (($i % 10000) == 0) # Для каждой 10000-й итерации (10 сек) вычисляем и печататем
{
$V = sqrt($Vx**2 + $Vy**2);
$Hkm = $H / 1000;

my $Vkas;
my $Vvert;
my $cosf;
my $Akasmod;
my $Armod;
$Akasmod = sqrt($Akx**2 + $Aky**2); # Определяем через скалярное умножение векторов выртикальну и горизонтальную составляющие скорости
$cosf = ($Akx * $Vx + $Aky * $Vy) / ($V * $Akasmod);
$Vkas = $V * $cosf;
$Armod = sqrt($Arx**2 + $Ary**2);
$cosf = ($Arx * $Vx + $Ary * $Vy) / ($V * $Armod);
$Vvert = $V * $cosf;

print "T = $T, H = $Hkm, V = $V, Vk = $Vkas, Vv = $Vvert\n";

}

if ($V < $Vmin) {$Vmin = $V} # Опрелеяем макс и мин. значения высоты и скорости
if ($V > $Vmax) {$Vmax = $V}
if ($H < $Hmin) {$Hmin = $H}
if ($H > $Hmax) {$Hmax = $H}


$i = $i + 1; # Увеличиваем итератор

} # Конец цикла

print "Vmin = $Vmin, Vmax = $Vmax, Hmin = $Hmin, Hmax = $Hmax\n"