package Orbit;

# Copyright (C) 2008  Bernhard W. Adams

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

use strict;
use Math::Trig;
use POSIX; # floor()
use aux;        # aux. routines:
                # epoch_to_date, date_to_epoch, dumpvars, msg, ems

my ($pi,$deg,$AU,$day,$np,$INF,$INVALID,$E_ZZ,$ORBIT_VAL,$ORBIT_INV);
$pi=2*asin(1.0);
$deg=$pi/180.0;
$AU=1.5e8;
$day=3600*24;
$np=1.0e-12; # num. prec.
$INF=1.0e30;
$INVALID=1.0e30;
$E_ZZ=-2;      # flag for e="0/0" in determination of ellipse, below
$ORBIT_INV=-1;
$ORBIT_VAL=1; # orbit is valid 
# scratchpad for given set of end points (space and time)
my (@tr_t,@orb_per,@ecc,@Tfa); # [0 ... (ntheta+2)]

##############################
#
# construct an orbit
# with the orbital parameters 
# $a    km,  s.m. axis 
# $e    rad, eccentricity
# $i    rad, inclination
# $an   rad, long. of ascending node
# $om   rad, arg. of periapsis
# $M0   rad, mean anomaly at epoch
# $tp   s,   time of periapsis closest to epoch
# $ep   s,   standard epoch
# !!!! only one of M0, tp needed, the other will be calculated
# $mu   km^2/s^2, grav. constant
# $I  in rad, not used here
# $L  in rad, mean long at time t_L
# $om
# $no
# $t_L
# $mu in km^3/s^2 
# $ps, $pe, $tra, $orb: filenames for positions, trajectory, orbit
#### $tp in s
#
##############################
sub new {
 my ($class) = shift; 
 my ($a,$dadt,$e,$dedt,$i,$didt,$an,$dandt,$ep,$om,$omb,$domdt,$L0,$M0,$tp,$mu,
     $dir,$name,$ps,$pe,$tra,$orb,$rps,$rpe,$rtra,$rorb,$deltav,$ttlog,$from,
     $to,$nlp,$diag)=@_;
 my ($M,$op,$t);
 
# aux::dumpvars("L0 M0",($L0,$M0));
 my $this = {
   _name  => undef,
   _group => undef,
   _a     => undef,
   _dadt  => undef,
   _e     => undef,
   _dedt  => undef,
   _i     => undef,
   _didt  => undef,
   _an    => undef,
   _dandt => undef,
   _ep    => undef,
   _om    => undef,
   _domdt => undef,
   _M0    => undef,
   _mu    => undef,
   _dir   => undef,
   _op    => undef,
   _PS    => "undef",
   _PE    => "undef",
   _TR    => "undef",
   _OR    => "undef",
   _RPS   => "undef",
   _RPE   => "undef",
   _RTR   => "undef",
   _ROR   => "undef",
   _deltav => "undef",
   _from  => undef,
   _to    => undef,
   _nlp   => undef,
   _complete => undef,   
   _ntheta  => undef,
   _rfx   => undef,
   _rfy   => undef,
   _rfz   => undef,
   _rfa   => undef,
   _rtx   => undef,
   _rty   => undef,
   _rtz   => undef,
   _rta   => undef,
   _DT    => undef,
   _tf    => undef,
   _tt    => undef,
   _ntheta => undef,
   _diag  => undef
 };
 # _M0 is the mean anomaly at ref. epoch _ep
 # M is the mean anomaly last calculated for given time _t
 bless $this, $class;
 $this->{_ntheta}=0;
 $this->{_complete} = $ORBIT_VAL; # set to 0 if not all orbital parameters are init.
                         # set to -1 if invalid orbit
 if(($this->{_e}=$e) eq "undef") {$this->{_complete}=0;}
 if(($this->{_dedt}=$dedt) eq "undef") {$this->{_complete}=0;}
 if(($this->{_a}=$a) eq "undef") {$this->{_complete}=0;}
 if(($this->{_dadt}=$dadt) eq "undef") {$this->{_complete}=0;}
 if(($this->{_i}=$i) eq "undef") {$this->{_complete}=0;} 
 if(($this->{_didt}=$didt) eq "undef") {$this->{_complete}=0;}
 if(($this->{_an}=$an) eq "undef") {$this->{_complete}=0;}
 if(($this->{_dandt}=$dandt) eq "undef") {$this->{_complete}=0;}
 if($om ne "undef") {
   $this->{_om} = $om;
   $this->{_domdt} = $domdt;
 }
 else {
   if($omb eq "undef") {$this->{_complete}=0;}
   $this->{_om} = $omb-$an;
   $this->{_domdt} = $domdt-$dandt;
 }
 if(($this->{_mu}=$mu) eq "undef") {$this->{_complete}=0;}
 $op=2*$pi*sqrt($a*$a*$a)/sqrt($mu);   # orb. period in seconds
 $this->{_dir}=$dir;
 $this->{_op}=$op;
 # http://perldoc.perl.org/perlsyn.html#Loop-Control for switch statement
 if(   (($L0 == "undef") || ($ep == "undef"))
    && (($M0 == "undef") || ($ep == "undef")) 
    && ($tp  == "undef")) {$this->{_complete}=0;}
 if(($L0 != "undef") && ($ep != "undef")) {
   $M0=$L0-$this->{_om}-$this->{_an};
   $tp=$ep-$op*$M0/(2*$pi);
 }
 else {
   if(($M0 != "undef") && ($ep != "undef")) {
     $tp=$ep-$op*$M0/(2*$pi);
   }
   else { # $tp
     print $name." tp\n";
     $M0=0;
     $ep=$tp;
   }
 }
 $this->{_M0} = $M0;
 $this->{_ep} = $ep;
 $this->{_PS} = $ps;
 $this->{_PE} = $pe;
 $this->{_TR} = $tra;
 $this->{_OR} = $orb;
 $this->{_RPS} = $rps;
 $this->{_RPE} = $rpe;
 $this->{_RTR} = $rtra;
 $this->{_ROR} = $rorb;
 $this->{_deltav} = $deltav;
 $this->{_name} = $name;
 $this->{_from} = $from;
 $this->{_ttlog} = $ttlog;
 $this->{_to} = $to;
 $this->{_nlp} = $nlp;
 if($diag eq "yes") {$this->{_diag}=1;} else {$this->{_diag}=0;}
 return $this;
}

##############################
# number of loops
##############################
sub nlp {my ($ob) = shift; $ob->{_nlp};}
##############################
# time of periapsis
##############################
sub tp {
 my ($ob) = shift;
 my $tp; 
 $tp=$ob->{_ep}-$ob->{_op}*$ob->{_M0}/(2*$pi);
 return $tp;
}
##############################
# orbital period, using the values a and mu from orbit initialization
##############################
sub P {my ($ob) = shift; $ob->{_op};}
##############################
# returns the object's name
##############################
sub name {my ($ob) = shift;  $ob->{_name};}
##############################
# returns the object's "from" name
##############################
sub from {my ($ob) = shift; $ob->{_from};}
##############################
# returns the object's "to" name
##############################
sub to {my ($ob) = shift; $ob->{_to};}
##############################
# returns the object's "to" name
##############################
sub valid {
  my ($ob) = shift; 
  if($ob->{_complete} == $ORBIT_VAL) {return 1;} else {return 0;}
}
################################
# position of orbiting body at time t
# in polar coordinates relative to this orbit
# t: time 
# returns array (r,M,E,T)
################################
sub posn_pol_rel {
 my $ME="posn_pol_rel"; # for error tracking
 my ($ob) = shift;      # this
 my ($t,$CF)=@_;        # time, called from (error tracking)
 my ($tr,$no,$r,$M,$E,$T,@res);
 my ($e,$a,$tp,$om,$mu,$op,$Dt); 

 $Dt=$t-$ob->{_ep};
 $a=$ob->{_a}+$ob->{_dadt}*$Dt;
 $e=$ob->{_e}+$ob->{_dedt}*$Dt;
# $tp=$ob->{_tp};
 $tp=$ob->{_ep}-$ob->{_op}*$ob->{_M0}/(2*$pi); 
 $om=$ob->{_om}+$ob->{_domdt}*$Dt;
 $mu=$ob->{_mu};
 $op=$ob->{_op};
 $tr=$t-$tp; # time since ref periapsis
 $no=floor($tr/$op); # no. of full orbits
 $tr-=$no*$op;      # time since last periapsis
 if($ob->{_dir}==1) {$M=2*$pi*$tr/$op;} # mean anomaly
 else {$M=2*$pi*(1-$tr/$op);}
 if(($M<0) || ($M>2*$pi)) {aux::ems("bug in ".$ME." ".$CF.", M =".$M);}
# $E=$ob->Kepler($M,$ME." ".$CF);
# $T=2*atan(tan($E/2)*sqrt(1+$e)/sqrt(fabs(1-$e))); # check this: fabs()
 ($E,$T)=$ob->Kepler($M,$ME." ".$CF);
  if($e<1) {$r=$a*(1-$e*$e)/(1+$e*cos($T));}
  else {$r=$a*($e*$e-1)/(1+$e*cos($T));}
 @res=($r,$M,$E,$T);
}
################
#
# solves the Kepler equation M=E-e*sin E
# for E with given M, e
# using Danby's method
# http://www.cdeagle.com/ccnum/pdf/demokep1.pdf
# http://www.projectpluto.com/kepler.htm
################
sub Kepler {
 my $ME="Kepler";
 my ($ob)=shift;
 my ($M,$CF)=@_;
#  my ($i,$E,$e,$m,$dm,$dE,$n,@Ea);
 my ($e,$E,$T,$f,$fp,$fpp,$fppp,$D,$Da,$nit,@res);
 my $maxit=300;

 $e=$ob->{_e};
 if($e<0) {die "negative eccentricity not allowed here";}
 $nit=0;
 if($e<1) { # ellipse
   if(sin($M)>0) {$E=$M+0.85*$e;}
   else          {$E=$M-0.85*$e;}
   do {
     $f=$E-$e*sin($E)-$M;
     $fp=1-$e*cos($E);
     $fpp=$e*sin($E);
     $fppp=$e*cos($E);
     $D=-$f/$fp;
     $Da=-$f/($fp+0.5*$D*$fpp);
     $E-=$f/($fp+0.5*$D*$fpp+$Da*$Da*$fppp/6);
     $nit+=1;
     if($nit>$maxit) {die;}
   } while(fabs($f)>$np);
#   $T=atan(sqrt(1-$e*$e)*sin($E))/(cos($E)-$e);
   $T=2*atan(sqrt((1+$e)/(1-$e))*tan(0.5*$E));
 }
 else { # hyperbola
   $E=log(2*$M/$e+1.8); # log is the natural logarithm, not base-10
   do{
     $f=$e*sinh($E)-$E-$M;
     $fp=$e*cosh($E)-1;
     $fpp=$e*sinh($E);
     $fppp=$e*cosh($E);
     $D=-$f/$fp;
     $Da=-$f/($fp+0.5*$D*$fpp);
     $E-=$f/($fp+0.5*$D*$fpp+$Da*$Da*$fppp/6);
     $nit+=1;
     if($nit>$maxit) {die;}     
  } while(fabs($f)>$np);
  $T=atan(sqrt($e*$e-1)*sinh($E))/($e-cosh($E));
   
 }
 return @res=($E,$T);
}
################
#
# solves the Kepler equation M=E-e*sin E
# for E with given M, e
# http://www.cdeagle.com/ccnum/pdf/demokep1.pdf
# http://www.projectpluto.com/kepler.htm
################
sub _Kepler {
 my $ME="Kepler";
 my ($ob)=shift;
 my ($M,$CF)=@_;
 my ($i,$E,$e,$m,$dm,$dE,$n,@Ea,$T,@res);
 my $maxit=300;
 
 $e=$ob->{_e};
 $E=$M;
 $n=0;
 do {
   $m=$E-$e*sin($E);
   $dm=1-$e*cos($E);
#   if(1/$dm>$pi) {$dm=1/$pi;} 
   $dE=($m-$M)/$dm;
   if(fabs($dE)>1) {$dE=$dE/fabs($dE);}
   $E=$E-$dE;
   if($E<0) {$E=0;}
   if($E>2*$pi) {$E=2*$pi;}   
   $Ea[$n]=$E;
   if($n>$maxit) {
     print "M = ".$M."\n";
     for($i=0;$i<scalar(@Ea);++$i) {print $Ea[$i]."\n";}
     die $ME." ".$CF." ".$e." ".$M;
   }
   $n+=1;
 } while(fabs($m-$M)>$M*$np);
# $T=atan(sqrt(1-$e*$e)*sin($E)/(cos($E)-$e));
 $T=2*atan(sqrt((1+$e)/(1-$e))*tan(0.5*$E));
 @res=($E,$T);
}
################################
#
# position of orbiting body at time t
# in cartesian coordinates relative to this orbit
# t: time (seconds, not epoch)
# returns array (x,y)
################################
sub posn_cart_rel {
 my $ME="posn_cart_rel"; # for error tracking
 my ($ob) = shift;      # this
 my ($t,$CF)=@_;        # time, called from (error tracking)
 my (@pol,@xy);
  
 @pol=$ob->posn_pol_rel($t,$ME." ".$CF);
 $xy[0]=$pol[0]*cos($pol[3]);
 $xy[1]=$pol[0]*sin($pol[3]); 
 @xy;
}
################################
#
# position of orbiting body at time t
# in polar coordinates in absolute coordinate system
# t: time (seconds, not epoch)
# returns array (r,phi,Theta)
################################
sub posn_pol_abs {
 my $ME="pos_pol_abs";
 my ($ob)=shift;
 my ($t,$CF)=@_;
 my (@rel,$T,$r,$theta,$phi,@res);
 my ($om,$Om,$i,$Dt);
 
 $Dt=$t-$ob->{_ep}; 
 $om=$ob->{_om}+$ob->{_domdt}*$Dt; # long. of periapsis
 $Om=$ob->{_an}+$ob->{_dandt}*$Dt; # long of ascending node
 $i=$ob->{_i}+$ob->{_didt}*$Dt; # inclination
 @rel=$ob->posn_pol_rel($t,$ME." ".$CF);
 $r=$rel[0]; # radius
 $T=$rel[3]; # true anomaly Theta
 $theta=arcsin(sin($T+$om)*sin($i));
 $phi=$Om+atan(sin($Om+$om)*sin($i)/cos($Om+$om));
 $res[0]=$r; $res[1]=$theta; $res[2]=$phi;
 @res;
}

################################
#
# returns absolute coordinates x,y of orbiting body at time t
# t: time (seconds, not epoch)
# returns array (x,y,z) in km
################################
sub posn_cart_abs {
 my $ME="posn_cart_abs";
 my ($ob) = shift;
 my ($t,$CF)=@_;
 my ($i,$om,$Om); 
 my ($Dt,@R,@Z,@OM);
# my $test=0;
 $Dt=$t-$ob->{_ep};
 $i=$ob->{_i}+$ob->{_didt}*$Dt;
 $om=$ob->{_om}+$ob->{_domdt}*$Dt;
 $Om=$ob->{_an}+$ob->{_dandt}*$Dt; 
 @R=$ob->posn_cart_rel($t,$ME." ".$CF);
# if(($R[0] > 0) && (fabs($R[0]/$R[1]) > 100)) {$test=1;}
 push @R,0; # z component
# if($test==1) {print "Om, om, i = ".($Om/$deg).", ".($om/$deg).", ".($i/$deg)."\n";}
# if($test==1) {print "--R = (".$R[0].", ".$R[1].", ".$R[2].")\n";}
 @Z=(0,0,1);
 @OM=(cos($Om),sin($Om),0);
 @R=aux::rotate(@R,@Z,$om+$Om);
# if($test==1) {print "--R = (".$R[0].", ".$R[1].", ".$R[2].")\n";} 
 @R=aux::rotate(@R,@OM,$i);
# if($test==1) {print "--R = (".$R[0].", ".$R[1].", ".$R[2].")\n";
#   my @tst=(-0.8,-0.6,0); 
#   @tst=aux::rotate(@tst,-0.8,-0.6,0,2);
#   print "tst = (".$tst[0].", ".$tst[1].", ".$tst[2].")\n";
# } 
 @R;
}

################################
# velocity of orbiting body at time t
# in rotated rectilinear coordinates relative to this orbit
# t: time 
# returns array (r,v_r,v_t,T,v,vp)
# r: radius
# v_r: radial vel.
# v_t: tang. vel.
# v: abs. velocity
# T: true anomaly
################################
sub vel_rrc_rel {
 my $ME="vel_rrc_rel"; # for error tracking
 my ($ob) = shift;      # this
 my ($t,$CF)=@_;        # time, called from (error tracking)
 my ($tr,$no,$r,$M,$E,$T,@res);
 my ($Dt,$mu,$a,$e,$rp,$vps,$vs,$v,$vr,$vt); 

 if($ob->{_complete}!=$ORBIT_VAL) {return ($ob->{_a},$INF,$INF,$INF,0);}
 $Dt=$t-$ob->{_ep};
 $mu=$ob->{_mu}; 
 $a=$ob->{_a}+$ob->{_dadt}*$Dt;
 $e=$ob->{_e}+$ob->{_dedt}*$Dt;
 if($e==1.0) {$e=1.0-1.0e-12;}
 if($e<1.0) {$rp=(1-$e)*$a;} else {$rp=($e-1)*$a;}
 $vps=(1+$e)*$mu/$rp; # square of periapsis vel.
 ($r,$M,$E,$T)=$ob->posn_pol_rel($t,$ME." ".$CF);
 $vs=$vps-2*$mu/$rp+2*$mu/$r; # velocity squared
 $v=sqrt($vs);                # velocity
 $vt=$rp*sqrt($vps)/$r;       # constant angular momentum 
 $vr=sqrt(fabs($vs-$vt*$vt)); # fabs to be tolerant to small neg. values due
 if(($T>$pi) || ($T<0)) {$vr*=-1;} # to roundoff errors
 if($ob->{_dir}==-1) {$vt*=-1; $vr*=-1;} 
 @res=($r,$vr,$vt,$v,$T);
}

################################
# velocity of orbiting body at time t
# in polar coordinates relative to this orbit
# t: time
# returns array (v,vp)
# v: abs. vel.
# vp: phi of vel.
################################
sub vel_pol_rel {
 my $ME="vel_pol_rel"; # for error tracking
 my ($ob) = shift;      # this
 my ($t,$CF)=@_;        # time, called from (error tracking)
 my ($r,$vr,$vt,$v,$T,@res); 
 
 ($r,$vr,$vt,$v,$T)=$ob->vel_rrc_rel($t,$ME." ".$CF);
 @res=($v,$T+0.5*$pi-atan($vr/$vt));
}

################################
# velocity of orbiting body at time t
# in absolute cartesian coordinates
# t: time 
# returns array (r,v_r,v_t)
################################
sub vel_cart_abs {
 my $ME="vel_cart_abs"; # for error tracking
 my ($ob) = shift;      # this
 my ($t,$CF)=@_;        # time, called from (error tracking)
 my ($Dt,$i,$om,$Om,$r,$vr,$vt,$v,$T,@V,@OM);
 
 $Dt=$t-$ob->{_ep};
 $i=$ob->{_i}+$ob->{_didt}*$Dt;
 $om=$ob->{_om}+$ob->{_domdt}*$Dt;
 $Om=$ob->{_an}+$ob->{_dandt}*$Dt;     
 ($r,$vr,$vt,$v,$T)=$ob->vel_rrc_rel($t,$ME." ".$CF);
 @V=($vr,$vt,0);
 @V=aux::rotate(@V,0,0,1,$T);
 @V=aux::rotate(@V,0,0,1,$om+$Om);
 @OM=(cos($Om),sin($Om),0);
 @V=aux::rotate(@V,@OM,$i);
}

##########################
# scratchpad for given set of end points (space and time)
# my (@from,@to,@obpl); # [0..2]
# my ($DT,$tf,$tt,$rfa,$rta);
# my (@tr_t,@orb_per,@Tfa); # [0 ... (ntheta+2)]
#
##########################
sub set_endpoints {
 my $ME="set_endpoints";
 my ($ob) = shift;
 my ($tfr,$tto,$rfx,$rfy,$rfz,$rtx,$rty,$rtz,$CF)=@_;
 my ($m,$Dt,$rfa,$rta,$DT,$Tf,$dTf,$ay);
 my ($e,$a,$tt,$op);
#  my ($A0,$A1,$x0,$x1,$y0,$y1,$T0,$T1,$r0,$r1,$cT0,$cT1,$x10,$x11,$x20,$x21);
 
 $ob->{_rfx}=$rfx;
 $ob->{_rfy}=$rfy;
 $ob->{_rfz}=$rfz;
 $ob->{_rfa}=$rfa=aux::vabs($rfx,$rfy,$rfz);
 $ob->{_rtx}=$rtx;
 $ob->{_rty}=$rty;
 $ob->{_rtz}=$rtz;
 $ob->{_rta}=$rta=aux::vabs($rtx,$rty,$rtz);
 $ob->{_tf}=$tfr;
 $ob->{_tt}=$tto;
 $Dt=$tto-$tfr;
 $ob->{_DT}=$DT=acos(aux::sc_prod($rfx,$rfy,$rfz,$rtx,$rty,$rtz)/($rfa*$rta));
 if($ob->{_diag}==1) { 
   open F,">".$ob->{_ttlog};           # for given endpoints in space & time
   print F "#Dt = ".($tto-$tfr)."\n";
   print F "#DT = ".$DT."\n";
   print F "# rfa,rta = ".$rfa.", ".$rta."\n";
   print F "#1:m, 2:Tf, 3:Tf+DT, 4:e, 5:a, 6:tt, 7:op-tt, 8:tt+op, 9:op\n";
 }
 $dTf=2*$pi/$ob->{_ntheta};
 for($m=0;$m<$ob->{_ntheta};++$m) {
   $Tf=$m*$dTf;
   ($a,$e,$tt,$op)=$ob->ell_pars($Tf);
   $Tfa[$m+1]=$Tf;
   $tr_t[$m+1]=$tt; # +1 to leave space at index 0
   $orb_per[$m+1]=$op; # for closure of cycle below
   $ecc[$m+1]=$e;
   # 1-7:  m, Tf, Tt, e, a, tt, tt-Dt
   # 8-10: op-tt-Dt, tt+op-Dt, op
   # 10-17: diag A0, A1, x0, x1, y0, y1, T0, T1,
   # 18-25       r0, r1, cosT0, cosT1, x10, x11, x20, x21
   if($ob->{_diag}==1) {
     print F $m."\t".$Tf."\t".($Tf+$DT)."\t".$e."\t".$a."\t".$tt."\t".($tt-$Dt).
           "\t".($op-$tt-$Dt)."\t".($tt+$op-$Dt)."\t".$op."\n";
#  	   "\t".$A0."\t".$A1."\t".$x0."\t".$x1."\t".$y0."\t".$y1.
#      	   "\t".$T0."\t".$T1."\t".$r0."\t".$r1."\t".$cT0."\t".$cT1.
#	   "\t".$x10."\t".$x11."\t".$x20."\t".$x21."\n";
   }
 }
 $Tfa[0]=$Tfa[$ob->{_ntheta}];
 $Tfa[$ob->{_ntheta}+1]=$Tfa[1];
 $tr_t[0]=$tr_t[$ob->{_ntheta}];
 $tr_t[$ob->{_ntheta}+1]=$tr_t[1]; # close cycle on both ends
 $orb_per[0]=$orb_per[$ob->{_ntheta}];
 $orb_per[$ob->{_ntheta}]=$orb_per[1];   # close cycle on both ends
 if($ob->{_diag}==1) {close F;}
}
##################
#
# times from - to in days (epoch)
# by definition of DT (acos above) the ellipses are all
# calculated so that the arc does not sweep across the sun
#
##############################
sub transfer_trajectory {
 my $ME="transfer_trajectory";
 my ($ob)=shift;
 my ($nl)=@_;
 my (@from,@to,@obpl,@pa,@an,@po); # all [0... 2]
 my ($Dt);
 my ($m,$a,$e,$tt,$op,$ta,$ea,$ay,$mm);
 my ($time,$Tl,$Tr,$Tm,$Tmn,$oTm,$ttl,$ttr,$ttm,$ottm,$expand);
 my ($dT,$trt,$hyp,$nit,$min_dt);
 my (@mt); # [1...ntheta]               # ntheta is def'd in new()
 my (@res); # [0..3]
 my $maxit=100;
    
 @from=($ob->{_rfx},$ob->{_rfy},$ob->{_rfz});
 @to=($ob->{_rtx},$ob->{_rty},$ob->{_rtz});
 @obpl=aux::cross(@from,@to);
 @obpl=aux::v_prod_s(@obpl,1/aux::vabs(@obpl));
 if($obpl[2]<0) {$nl=-$nl-1;} 
 if($nl>-1) {$ob->{_dir}=1;} else {$ob->{_dir}=-1;}
 @an=aux::cross(0,0,1,@obpl); # vector towards ascending node
 if(aux::vabs(@an)==0) {@an=(1,0,0);} # if origin, dest. trajectories have i=0
 @an=aux::v_prod_s(@an,1/aux::vabs(@an)); # unit vector ...
 $ay=aux::sc_prod(0,1,0,@an);
 if($ay>0) {$ob->{_an}=acos(aux::sc_prod(@an,1,0,0));}
 else      {$ob->{_an}=2*$pi-acos(aux::sc_prod(@an,1,0,0));} 
 $ob->{_dandt}=0;
 if($obpl[2]>0) {$ob->{_i}=asin(aux::vabs(aux::cross(0,0,1,@obpl)));}
 else {$ob->{_i}=$pi-asin(aux::vabs(aux::cross(0,0,1,@obpl)));}
 $ob->{_didt}=0;
 $Dt=$ob->{_tt}-$ob->{_tf};
 $min_dt=$INF; 
 # $nl=0... : add $nl full orb. periods
 # $nl=0    : less than 1 orb. period forward
 # $nl=-1   : less than 1 orb. period backward
 # $nl<-1   : less than 1 orb. backw. PLUS -$nl-1 full orb backw.
 # find closest match over $m
 for($m=1;$m<=$ob->{_ntheta};++$m) {   #######
   if($ecc[$m]<1) { # ellipse
     if($nl>=0) {$time= $tr_t[$m]+$nl*$orb_per[$m];}
     else       {$time=-$nl*$orb_per[$m]-$tr_t[$m];}
   }
   else { # hyperbola
     if(($nl==0)||($nl==-1)){if(($time=2*($nl+0.5)*$tr_t[$m])<0) {$time=$INF;}}
     else {$time=$INF;}
   }
   $mt[$m]=$time-$Dt;
 }
 for($m=1;$m<=$ob->{_ntheta};++$m) {   ########
   if((fabs($mt[$m])<fabs($min_dt))&&
      ($mt[$m-1]<0.5*$INF)&&($mt[$m+1]<0.5*$INF)&&
      ($ecc[$m-1]>=0)&&($ecc[$m]>=0)&&($ecc[$m+1]>=0)) {
     $min_dt=$mt[$m];
     $mm=$m;
   } # closest transit time if
     # tr. times left and right are both valid
 } # end for($m)
 $Tl=$Tfa[$mm-1]; $ttl=$mt[$mm-1];
 $Tm=$Tfa[$mm+0]; $ttm=$mt[$mm+0];
 $Tr=$Tfa[$mm+1]; $ttr=$mt[$mm+1];
 if($ob->{_diag}==1) {
   open F,">exp.log";
   print F "#nl =   ".$nl."\n";
   print F "#Dt =   ".$Dt."\n";
   print F "#Tl =   ".$Tl."\n";
   print F "#Tm =   ".$Tm."\n";     
   print F "#Tr =   ".$Tr."\n";   
   print F "#ttl =  ".$ttl."\n";
   print F "#ttm =  ".$ttm."\n";
   print F "#ttr =  ".$ttr."\n";
 } 
 $expand=0;
 if($ttl*$ttr>0) {# no zero crossing between Tl and Tr ? --> expand interval
   if($ob->{_diag}==1) {
     print F "#min_dt = ".$min_dt."\n";
     print F "#Dt =    ".$Dt."\n";   
     print F "#ecc_l =    ".$ecc[$mm-1]."\n";   
     print F "#ecc_m =    ".$ecc[$mm]."\n";   
     print F "#ecc_r =    ".$ecc[$mm+1]."\n";      
   }
   if(($ttl+$ttr)*($ttr-$ttl)>0) {$expand=-1;} # mean * 'slope'   
   else {$expand=1;} # expand Theta interval to left
   if($ob->{_diag}==1) {print F "#expand = ".$expand."\n";}
   ($a,$e,$trt,$op)=$ob->ell_pars($Tl);
   if($e<0) {die;}   ##############?????
   ($a,$e,$trt,$op)=$ob->ell_pars($Tr);
   if($e<0) {die;}
   $nit=0;
   $dT=($Tr-$Tl); # 2* Tf resolution
   do { # extrapolation ######### here
     $dT/=2; # 1* Tf resolution
     if($expand == 1) { # expand interval towards right
#       print ">";
       ($a,$e,$trt,$op)=$ob->ell_pars($Tr+$dT);
       if($e==$E_ZZ) {
         $ob->{_complete}=$ORBIT_INV;
       	 if($ob->{_diag}==1) {
       	   print F "# ZZ>\n";
           print F $Tl."\t".$Tr."\t".$a."\t".$e."\t".$trt."\t".$op."\n";
       	 }
       	 return @res=(0,-1,-1,$INF);
       }
       if($e>=0) { # if valid ellipse/hyp. 
         $Tr+=$dT;
         if($e<1) { # ellipse
           if($nl>=0) {$ttr= $trt+$nl*$op-$Dt;}
           else       {$ttr=-$nl*$op-$trt-$Dt;}
         }
         else { # hyperbola
           if(($nl==0)||($nl==-1)) {if(($ttr=2*($nl+0.5)*$trt)<0) {$ttr=$INF;}}
           else {$ttr=$INF;}
         }
       }
     }
     else {
#       print "<";
       ($a,$e,$trt,$op)=$ob->ell_pars($Tl-$dT);
       if($e==$E_ZZ) {
         $ob->{_complete}=$ORBIT_INV;
       	 if($ob->{_diag}==1) {
       	   print F "# ZZ<\n";
       	   print F "# rf, rt = ".$ob->{_rfa}."\t".$ob->{_rta}."\n";
           print F $Tl."\t".$Tr."\t".$a."\t".$e."\t".$trt."\t".$op."\n";
       	 }
       	 return @res=(0,-1,-1,$INF);
       }
       if($e>=0) {
         $Tl-=$dT;
         if($e<1) { # ellipse
           if($nl>=0) {$ttl= $trt+$nl*$op-$Dt;}
           else       {$ttl=-$nl*$op-$trt-$Dt;}
         }
         else { # hyperbola
           if(($nl==0)||($nl==-1)) {if(($ttl=2*($nl+0.5)*$trt)<0) {$ttl=$INF;}}
           else {$ttl=$INF;}
         } 
       }
     }
     # 1-6
     if($ob->{_diag}==1) {
       print F $Tl."\t".$Tr."\t".$a."\t".$e."\t".$trt."\t".$op.
               "\t".$ttl."\t".$ttr."\n";
     }
     # 7-8
     $nit+=1;
     if($nit>$maxit) {close F; return @res=(0,-1,-1,$INF);}
#   } while(($ttr+$ttl)*($ttr-$ttl)>0);
   } while(($ttl*$ttr>0) || ($ttl == $INF) || ($ttr == $INF));
   if($ob->{_diag}==1) {close F;}
 } # end if($ttl*$ttr>0) 
# print "----------";
 # now we need only to interpolate, not extrapolate
 $nit=0;
# print "tt_{l,m,r} = ".$ttl.", ".$ttm.", ".$ttr."\n";
# print "T_{l,m,r}, ttm = ".$Tl.", ".$Tm.", ".$Tr." ".$ttm."\n"; 
 if($ob->{_diag}==1)
   {open F, ">int.log"; print F "# Tl, Tr = ".$Tl.", ".$Tr."\n";}
 if($Tm<$Tl) {$Tm+=2*$pi;} # for wraparound 2*pi - \epsilon --> + \epsilon
 if($Tr<$Tl) {$Tr+=2*$pi;} # for wraparound 2*pi - \epsilon --> + \epsilon
 if(($Tl>$Tr) || ($Tm>$Tr) || ($Tm<$Tl)) {die $Tl.", ".$Tm.", ".$Tr;} 
 if($ob->{_diag}==1) {
   print F "#expand = ".$expand."\n";
   print F "#nl = ".$nl."\n";
   print F (-1).", ".$Tl.", ".$Tm.", ".$Tr.", ".
            $ttl.", ".$ttm.", ".$ttr.", ".$ecc[$mm-1].", ".
	    $ecc[$mm].", ".$ecc[$mm+1]."\n";} 
 do { # we are entering this loop if $ttl*$ttr < 0
   $Tm=0.5*($Tl+$Tr);
   ($a,$e,$trt,$op)=$ob->ell_pars($Tm);   
   if($e<0) {
     $ob->{_complete}=$ORBIT_INV;
     return @res=(0,-1,-1,$INF); # (a, e, Tm, op)
   }
   if($a==-1) {die $e.", ".$trt;}
   if($e<1) { # ellipse
     if($nl>=0) {$ttm= $trt+$nl*$op-$Dt;}
     else       {$ttm=-$nl*$op-$trt-$Dt;}
   }
   else { # hyperbola
     if(($nl==0)||($nl==-1)) {if(($ttm=2*($nl+0.5)*$trt)<0) {$ttm=$INF;}}
     else {$ttm=$INF;}
   }
   if($ttm*$ttr>0) {$Tr=$Tm; $ttr=$ttm;}
   else            {$Tl=$Tm; $ttl=$ttm;}
   if($ob->{_diag}==1) 
     {print F $nit.", ".$Tl.", ".$Tm.", ".$Tr.", ".
              $ttl.", ".$ttm.", ".$ttr."\n";}     
   $nit+=1;
  }
 while(($nit<$maxit) && (fabs($ttm/$Dt) > 1.0e-12) && ($Tl != $Tr));
 if($ob->{_diag}==1) {close F;}
 if($nit==$maxit) {print F "max no. of iterations ".$maxit."\n";}
 $ob->{_a}=$a;
 $ob->{_dadt}=0;
 $ob->{_e}=$e;
 $ob->{_dedt}=0;
 $ob->{_op}=$op; 
 # calculate orbital parameters
 @pa=aux::rotate(@from,@obpl,-$Tm);
 @pa=aux::v_prod_s(@pa,1/aux::vabs(@pa)); # unit vector ...
 my $fa=aux::vabs(@from);
 @po=aux::cross(@obpl,@an);
 if(aux::sc_prod(@po,@pa)>0) {$ob->{_om}=acos(aux::sc_prod(@an,@pa));}
 else {$ob->{_om}=2*$pi-acos(aux::sc_prod(@an,@pa));}
 $ob->{_domdt}=0;
 $ob->{_ep}=$ob->{_tf};
 $ta=$Tm; # true anomaly
 if($e<1) { # ellipse
   $ea=2*atan(tan($ta/2)*sqrt(1-$e)/sqrt(1+$e)); # eccentric anomaly
   if($ob->{_dir}==1) {$ob->{_M0}=$ea -$e*sin($ea);}
   else {$ob->{_M0}=2*$pi-($ea -$e*sin($ea));}
   $ob->{_op}=2*$pi*sqrt($a*$a*$a)/sqrt($ob->{_mu});
   $ob->{_complete}=1; 
 }
 @res=($a,$e,$Tm,$op);
}

######################
#
#
#
######################
sub _transfer_trajectory {
 my $ME="transfer_trajectory";
 my ($ob)=shift;
 my ($tfr,$tto,$rfx,$rfy,$rfz,$rtx,$rty,$rtz,$nl,$CF)=@_;
 my ($rfa,$rta,$Dt,$DT);
 my ($c,$e,$a,$ay);
 my (@from,@to,@obpl,@an);

#
 $ob->{_rfx}=$rfx;
 $ob->{_rfy}=$rfy;
 $ob->{_rfz}=$rfz;
 $ob->{_rfa}=$rfa=aux::vabs($rfx,$rfy,$rfz);
 $ob->{_rtx}=$rtx;
 $ob->{_rty}=$rty;
 $ob->{_rtz}=$rtz;
 $ob->{_rta}=$rta=aux::vabs($rtx,$rty,$rtz);
 $ob->{_tf}=$tfr;
 $ob->{_tt}=$tto;
 $Dt=$tto-$tfr;
 $ob->{_DT}=$DT=acos(aux::sc_prod($rfx,$rfy,$rfz,$rtx,$rty,$rtz)/($rfa*$rta));
#
 @from=($rfx,$rfy,$rfz);
 @to=($rtx,$rty,$rtz);
 @obpl=aux::cross(@from,@to);
 @obpl=aux::v_prod_s(@obpl,1/aux::vabs(@obpl));
 if($obpl[2]<0) {$nl=-$nl-1;} 
 if($nl>-1) {$ob->{_dir}=1;} else {$ob->{_dir}=-1;}
 @an=aux::cross(0,0,1,@obpl); # vector towards ascending node
 @an=aux::v_prod_s(@an,1/aux::vabs(@an)); # unit vector ...
 $ay=aux::sc_prod(0,1,0,@an);
 if($ay>0) {$ob->{_an}=acos(aux::sc_prod(@an,1,0,0));}
 else      {$ob->{_an}=2*$pi-acos(aux::sc_prod(@an,1,0,0));} 
 $ob->{_dandt}=0;
 if($obpl[2]>0) {$ob->{_i}=asin(aux::vabs(aux::cross(0,0,1,@obpl)));}
 else {$ob->{_i}=$pi-asin(aux::vabs(aux::cross(0,0,1,@obpl)));}
 $ob->{_didt}=0;
 $Dt=$ob->{_tt}-$ob->{_tf};
 
}
##################
#
#
###################
sub ell_pars {
  my $ME="ell_pars";
  my ($ob)=shift;
  my ($Tf)=@_;
  my ($m,$a,$b,$tt,$e,$op,$den,@res);
  my (@r,@T,@A,@x,@y); # [0 ... 1]
             	     
 $T[0]=$Tf; $T[1]=$Tf+$ob->{_DT}; # --here
 $r[0]=$ob->{_rfa}; $r[1]=$ob->{_rta};    
 if($T[1]>2*$pi) {$T[1]-=2*$pi;}
 $den=$r[1]*cos($T[1])-$r[0]*cos($T[0]);
 # return invalid result if denominator is too large, especially for cases
 # of e ~ "0/0"
 if($den*$den/($r[0]*$r[1])<1.0e-25) {return @res=(0,$E_ZZ,-1,$INF);}
 $e=($r[0]-$r[1])/$den;
 if($e==1.0) {$e=1.0-1.0e-12;}
 $a=$r[0]*(1+$e*cos($T[0]))/(1-$e*$e); 
 if($e>=1) {$a*=-1;} # hyperbola
 else {if($a<=1.0e-12) {$a=1;}} # kludge 
 if(($e>=0) && ($a>0)) { # if $a<=0 Theta too large, no point on hyp.
# if($e>=0) {
   if($e*$e<1) { # ellipse
     $b=$a*sqrt(1-$e*$e);
     for($m=0;$m<2;++$m) {
       $x[$m]=$r[$m]*cos($T[$m])+$a*$e; # this produces a glitch
#       if(($x[$m]>0) && (fabs($T[$m]-$pi)<0.1)) {
#         die $r[$m].", ".$T[$m].", ".cos($T[$m]).", ".$a.", ".$e;
#       }
       $y[$m]=$r[$m]*sin($T[$m]);
       if(fabs($x[$m]/$a)>1.0+1.0e-15) {
#         print "x, a, e, T = ".$x[$m].", ".$a.", ".$e.", ".$T[$m]."\n";
       }
       if(fabs($x[$m]/$a)>1.0-1.0e-12) 
         {if($x[$m]>0) {$x[$m]=$a;} else {$x[$m]=-$a;}}
       if(fabs($x[$m]/$a)>1) {
       	 die "x, a, e, r = ".$x[$m].", ".$a.", ".$e."".$r[$m].
	     ", x/a = ".($x[$m]/$a);
       }
       $A[$m]=0.5*(fabs($a*$b*acos($x[$m]/$a))-fabs($a*$e*$y[$m]));
       if($y[$m]<0) {$A[$m]=$pi*$a*$b-$A[$m];} # area of ellipse sector
     } # end for
     if($A[1]<$A[0]) {$A[1]+=$pi*$a*$b;}
     $op=2*$pi*sqrt($a*$a*$a/$ob->{_mu}); # orbital period
   }
   else { # hyperbola
     $b=$a*sqrt($e*$e-1);
     for($m=0;$m<2;++$m) {     
       $x[$m]=$r[$m]*cos($T[$m])-$a*$e; if($x[$m]>0) {die;}
       $y[$m]=$r[$m]*sin($T[$m]);
       if((($x[$m]==0)) || (fabs($y[$m]/$x[$m])>fabs($b/$a))) 
         {return @res=(0,-1,-1,$INF);}
       else {$A[$m]=0.5*(fabs($a*$e*$y[$m])-fabs($a*$b*acosh(-$x[$m]/$a)));}
       if($y[$m]<0) {$A[$m]*=-1;}
     } # end for
     $op=$INF;
   } # end else
   $tt=2*($A[1]-$A[0])/sqrt(fabs($a*$ob->{_mu}*(1-$e*$e)));
 } # end if(($e>=0) && ($a>=0)) {
 else {return @res=(-1,-1,-1,$INF);}
 if($tt>$INF/2) {die;}
 @res=($a,$e,$tt,$op);
}

return 1;