# определение необходимых delta-v, прилагаемых по курсу движения.
sub hohmann_dv {
  my ($rc, $rt) = @_; # $rc - Radius of Current orbit; $rt - Radius of Target orbit;

  my $dv1 = sqrt($u / $rc) * (sqrt(2 * $rt / ($rc + $rt)) - 1);
  my $dv2 = sqrt($u / $rt) * (1 - sqrt(2 * $rc / ($rc + $rt)));

  return ($dv1, $dv2);
}

# приложение импульса силой $dv из точки $x, $y в направлении $cw перпендикулярно радиус-вектору
sub direct_dv {
  my ($dv, $x, $y, $cw) = @_;
  return map { $_ * $dv } vec_norm(vec_vp($x, $y, 0, 0, 0, $cw))
}

...

  if (1 == $stage && !defined $x1) {
    # first impulse, circular -> elliptic
    ($dvx, $dvy) = direct_dv($dv[0], $x, $y, $cw);
    vm_input(0x2 => $dvx, 0x3 => $dvy);
    ($x1, $y1) = ($x, $y);
    $stage++;
  }

  if (2 == $stage && 0 > $cw * vec_sp($x1, $y1, $x, $y)) {
    # vec_sp - скалярное произведение векторов, $x1, $y1 - координаты первого импульса, $cw - clockwiseness
    # second impulse, elliptic -> circular
    ($dvx, $dvy) = direct_dv($dv[1], $x, $y, $cw);
    vm_input(0x2 => $dvx, 0x3 => $dvy);
    $stage++;
  }

...

# а это переход с произвольной эллиптической на круговую.
# то есть первый импульс рассчитывается по той же формуле гомана,
# только для «якобы целевой» орбиты радиусом сильно больше нужного,
# а второй импульс рассчитывается непосредственно в месте пересечения текущей и целевой орбит

  if (2 == $stage && $v > (($tgt - $alt))) {
  # $v - модуль скорости в данный момент,
  # $alt - наше текущее расстояние от начала координат,
  # $tgt - желаемое расстояние от начала координат

    my ($vx, $vy) = ($x0 - $x, $y0 - $y);
    # $x0, $y0 - наши координаты на прошлом тике,
    # $x, $y - на текущем

    my ($vdx, $vdy) = map { $_ * $vf } vec_norm(vec_vp($x, $y, 0, 0, 0, $cw));
    # $vdx = velocity desired x, $vdy = velocity desired y
    # vec_vp - векторное произведение, vec_norm - нормализация вектора по длине

    ($dvx, $dvy) = ($vdx - $vx, $vdy - $vy);
    vm_input(0x2 => $dvx, 0x3 => $dvy);
    $stage++;
  }