#!/usr/bin/perl

use Math::Trig;
# use Math::Gsl::Sf;
# use Math::Gsl::Sf qw(:Gamma);
# use Math::Complex;

use Km;
use af;
# misc. parameters
$pi=2*asin(1.0);
$deg=$pi/180.0;
# $AU=1.49598e8;
# $AU=1.4958e8;
$AU=1.49597870691e8;
$day=3600*24;
$np=1.0e-12; # num. prec.
$npt=100;    # plot points
$LARGE=1.0e30;
$INVALID=1.0e30;
# reference date: 13.04.2029
# 

# Sun
$mu=1.32712440018e11;

# JPL:  http://ssd.jpl.nasa.gov/txt/p_elem_t1.txt
# Origin (ex. earth)
$E_a=1.0000001124*$AU;           # semimajor axis, Wiki
$E_a=1.00000261*$AU;             # s.m., JPL
$E_e=0.01671123;                 # eccentricity, JPL
$E_I=-0.00001531*$deg;           # inclination, not used here, JPL
$E_L=100.46457166*$deg;          # L at J2000, JPL
$E_om=102.93768193*$deg;         # peri long., JPL
$E_no=0.0*$deg;                  # asc. node, JPL
$E_tL=2451545.0;                 # J2000
# $E_tp=2454104.67344;             # time of perihelion
# Besselian
# $E_tp=2454104.22878957;           # time of perihelion
# Julian
# $E_tp=2454104.75-5;           # time of perihelion
# 2007-01-03, 20h (http://aa.usno.navy.mil/data/docs/EarthSeasons.html)
# 1950.0 is 2433282.42344

# Destination (ex. asteroid)
# http://ssd.jpl.nasa.gov/sbdb.cgi?sstr=99942;orb=1
$A_a=0.92226141513913*$AU;       # 
$A_e=0.191059415070228;          #
$A_I=3.3*$deg;                   # not used here
$A_om=126.3855713119676*$deg;
$A_M=307.3630785264554*$deg;
$A_L=$A_M+$A_om-2*$pi;
# L= M+omega, from M=307.3630785264554
$A_no=204.4591523029436*$deg;
$A_tL=2454200.5;
# $A_tL=2451545.0;
# $A_or=126.3855713119676*$deg;  # orientation of major axis
# $A_or=135.0*$deg;  # orientation of major axis
### $A_tp=2454247.800688934931;    # epoch, 2007-May-27.30068894
# $A_ep=2462239.17344;           # 2029-04-13
##########################################
#
#########################################################
#
$Origin=Km->     new($E_a,$E_e,$E_I,$E_L,$E_om,$E_no,$E_tL,$mu);
$Destination=Km->new($A_a,$A_e,$A_I,$A_L,$A_om,$A_no,$A_tL,$mu);
$transfer[0]=Km->   new(1,$INVALID,  -1,  -1,  -1,   -1,$INVALID,$mu); # dummy parameters
$transfer[1]=Km->   new(1,$INVALID,  -1,  -1,  -1,   -1,$INVALID,$mu); # dummy parameters
# transfer[2] is scratchpad
$transfer[2]=Km->   new(1,$INVALID,  -1,  -1,  -1,   -1,$INVALID,$mu); # dummy parameters
#####################
$pl_traj_of_origin="E_traj";
$pl_traj_of_dest="A_traj";
$pl_transfer_traj[0]="T_traj_0";
$pl_transfer_traj[1]="T_traj_1";
$pl_transfer_orbit[0]="T_orbit_0";
$pl_transfer_orbit[1]="T_orbit_1";
$pl_posn_of_origin="E_posn";
$pl_posn_of_dest="A_posn";
$pl_time_diff="td";
$plot_delta_v="dv";
########################
$task=$ARGV[0];
if(($ARGV[1] ne "-date")&&($ARGV[1] ne "-epoch")) {die "invalid time format";}
if($ARGV[1] eq "-date") {
  @ymd=split(/[,\/]/,$ARGV[2]);
  $epoch=af::date_to_epoch($ymd[0],$ymd[1],$ymd[2]);
  # print "epoch = ".$epoch."\n";
  $interval=$ARGV[3];
}
else {
  $epoch=$ARGV[2];
  $interval=$ARGV[3]-$epoch;
}
###########
$ndvd=0; $ndva=0;
$nloop=0; # multiple loops in transfer orbit
for($m=4;$m<scalar(@ARGV);++$m) {
  if($ARGV[$m] eq "-no_dv") {
    if(scalar(@ARGV)<$m+3) {die "invalid command line parameter";}   
    $ndvd=$ARGV[$m+1]; $ndva=$ARGV[$m+2];
    if(($ndvd<1)||($ndvd>1000)||($ndva<1)||($ndva>1000))
      {die "invalid value ".$no_pt_dv." for no_dv";}
  }
  if($ARGV[$m] eq "-allow_multiple_loops") {
    if(scalar(@ARGV)<$m+2) {die "invalid command line parameter";}    
    $nloop=$ARGV[$m+1];
  }
}
$td=$epoch*86400; # time in seconds
$ta=$td+$interval*86400; # time in seconds
#
if(($ndvd==0)||($ndva==0)) {$ndvd=50; $ndva=50;}
######################
#
#
#######################
if($task eq "-dv") { # delta v
  open F,">".$plot_delta_v;
  print F "1:  m: index\n2:  n: index\n3:  td: time of departure (epoch)\n".
          "4:  ta: time of arrival (epoch)\n".
          "5:  dvda0: delta v at departure, sol.01\n".
          "6:  dvda1: delta v at departure, sol. 1\n".
          "7:  dvaa0: delta v at arrival, sol. 0\n".
          "8:  dvaa1: delta v at arrival, sol. 1\n".
          "9:  dvs0: total delta v, sol. 0\n".
          "10: dvs1: total delta v, sol. 1\n".
          "11: tdiff0: time mismatch, sol. 0\n".
          "12: tdiff1: time mismatch, sol. 1\n";
  $dept_range=$ARGV[4];
  $arrv_range=$ARGV[5];    
  $mindv_i=$LARGE;  $mindv_r=$LARGE; 
  for($m=0;$m<$ndvd;++$m) {
    $tdv=$td+((($m/($ndvd-1))-0.5)*$dept_range)*86400;
    for($n=0;$n<$ndva;++$n) {
      $tav=$tdv+($interval+(($n/($ndva-1))-0.5)*$arrv_range)*86400;      
      if($tav<$tdv) {af::die("tav<tdv","main");}
      @po=$Origin->posn_orb($tdv,"main dept"); # here !!!
      $Thd=$po[2];
      $omd=$po[3];
      $rd=$po[4];                            # dist. from sun at departure
      @po=$Destination->posn_orb($tav,"main arrv");
      $Tha=$po[2];
      $oma=$po[3];
      $ra=$po[4];                            # dist. from sun at arrival
      @res=
        find_matching_trajectories($rd,$ra,$Thd,$Tha,$omd,$oma,$tdv,$tav,$n_or,"");
      # make $transfer object hold the values determined above
      for($o=0;$o<2;++$o) { 
        $beta_o=$res[$o]; $tdiff[$o]=$res[$o+2];
        if($tdiff[$o]!=$INVALID) {
          $Th=$Thd+$omd-$beta_o;
          $dth=$Tha+$oma-$Thd-$omd;
          $tt[$o]=$transfer[$o]->transfer_orbit($rd,$ra,$Th,$dth,$tdv,$beta_o,"main");
	}
      # determine delta v
        if(($tt[$o]!=$INVALID)&&($tdiff[$o]!=$INVALID)) {
          @dvd=delta_v($Origin,$transfer[$o],$tdv,"main_dvd");
	  @dva=delta_v($transfer[$o],$Destination,$tav,"main_dva");          
	  $vdb[0+2*$o]=$dvd[2]; $vdb[1+2*$o]=$dvd[3]; # before departure burn
	  $vda[0+2*$o]=$dvd[4]; $vda[1+2*$o]=$dvd[5]; # after dept. burn
	  $vab[0+2*$o]=$dva[2]; $vab[1+2*$o]=$dva[3]; # before arrv. burn
	  $vaa[0+2*$o]=$dva[4]; $vaa[1+2*$o]=$dva[5]; # after arrv. burn 
	  $dvda[$o]=sqrt($dvd[0]*$dvd[0]+$dvd[1]*$dvd[1]);
          $dvaa[$o]=sqrt($dva[0]*$dva[0]+$dva[1]*$dva[1]);
	  $dvd_[2*$o+0]=$dvd[0]; $dvd_[2*$o+1]=$dvd[1];
	  $dva_[2*$o+0]=$dva[0]; $dva_[2*$o+1]=$dva[1];
	}
        else {$dvda[$o]=20; $dvaa[$o]=20;}
      }
      # print to file
      print F $m."\t".$n."\t".($tdv/86400)."\t".($tav/86400)."\t".
	      $dvda[0]."\t".$dvda[1]."\t".$dvaa[0]."\t".$dvaa[1]."\t".
	      ($dvda[0]+$dvaa[0])."\t".($dvda[1]+$dvaa[1])."\t".
	      $tdiff[0]."\t".$tdiff[1]."\n";
      # low score for intercept (no terminal velocity matching)
      if($dvda[0]<$mindv_i) {
        $mindv_i=$dvda[0];
	$mintd_i=$tdv; $minta_i=$tav;  
        @mindvd_i=($dvd_[0],$dvd_[1]); @mindva_i=($dva_[0],$dva[1]);
	$mintdiff_i=$tdiff[0];
	$minot_i=0;
        @vdbn=@vdb; @vdan=@vda;	
        @vabn=@vdb; @vaan=@vda;
      }
      if($dvda[1]<$mindv_i) {
        $mindv_i=$dvda[1];
	$mintd_i=$tdv; $minta_i=$tav;  
        @mindvd_i=($dvd_[2],$dvd_[3]); @mindva_i=($dva_[2],$dva[3]);
	$mintdiff_i=$tdiff[1];
        $minot_i=1;
        @vdbn=@vdb; @vdan=@vda;	
        @vabn=@vdb; @vaan=@vda;	      
      }      
      # low score for rendezvous
      if($dvda[0]+$dvaa[0]<$mindv_r) {
        $mindv_r=$dvda[0]+$dvaa[0];
	$mintd_r=$tdv; $minta_r=$tav;
        @mindvd_r=($dvd_[0],$dvd_[1]); @mindva_r=($dva_[0],$dva[1]);
	$mintdiff_r=$tdiff[0];
        $minot_r=0;
        @vdbm=@vdb; @vdam=@vda;	
        @vabm=@vdb; @vaam=@vda;	      
      }    
      if($dvda[1]+$dvaa[1]<$mindv_r) {
        $mindv_r=$dvda[1]+$dvaa[1];
	$mintd_r=$tdv; $minta_r=$tav;  
        @mindvd_r=($dvd_[2],$dvd_[3]); @mindva_r=($dva_[2],$dva[3]);
	$mintdiff_r=$tdiff[0];
        $minot_r=1;
        @vdbm=@vdb; @vdam=@vda;	
        @vabm=@vdb; @vaam=@vda;	
      }    
    }
    print F "\n";
  }
  close F;
  $minepd_i=$mintd_i/86400; $mindtd_i=af::epoch_to_date($minepd_i);
  $minepa_i=$minta_i/86400; $mindta_i=af::epoch_to_date($minepa_i);
  $minepd_r=$mintd_r/86400; $mindtd_r=af::epoch_to_date($minepd_r);
  $minepa_r=$minta_r/86400; $mindta_r=af::epoch_to_date($minepa_r);
  print "min. delta v (intercept) =    ".$mindv_i."\n".
        "departure delta v = (".$mindvd_i[0].", ".$mindvd_i[1].")\n".
        " (".$vdbn[0+2*0].", ".$vdbn[1+2*0].
        ") - (".$vdan[0+2*0].", ".$vdan[1+2*0]."), ".
        sqrt(sq($vdbn[0+2*0]) + sq($vdbn[1+2*0])).", ".
        sqrt(sq($vdan[0+2*0]) + sq($vdan[1+2*0])).        
        "\n".
        "arrival delta v   = (".$mindva_i[0].", ".$mindva_i[1].")\n".
        " (".$vabn[0+2*0].", ".$vabn[1+2*0].
        ") - (".$vaan[0+2*0].", ".$vaan[1+2*0]."), ".
        sqrt(sq($vabn[0+2*0]) + sq($vabn[1+2*0])).", ".
        sqrt(sq($vaan[0+2*0]) + sq($vaan[1+2*0]))."\n".
        "departure at epoch   ".$minepd_i."\t(".$mindtd_i.")\n".
        "arrival at epoch     ".$minepa_i."\t(".$mindta_i.")\n".
        "time mismatch =      ".$mintdiff_i."\n";
  print "min. delta v (rendezvous) =    ".$mindv_r."\n".
        "departure delta v = (".$mindvd_r[0].", ".$mindvd_r[1].")\n".
        " (".$vdbm[0+2*0].", ".$vdbm[1+2*0].
        ") - (".$vdam[0+2*0].", ".$vdam[1+2*0]."), ".
        sqrt(sq($vdbm[0+2*0]) + sq($vdbm[1+2*0])).", ".
        sqrt(sq($vdam[0+2*0]) + sq($vdam[1+2*0]))."\n".
        "arrival delta v   = (".$mindva_r[0].", ".$mindva_r[1].")\n".
        " (".$vabm[0+2*0].", ".$vabm[1+2*0].
        ") - (".$vaam[0+2*0].", ".$vaam[1+2*0]."), ".
        sqrt(sq($vabm[0+2*0]) + sq($vabm[1+2*0])).", ".
        sqrt(sq($vaam[0+2*0]) + sq($vaam[1+2*0]))."\n".
        "departure at epoch   ".$minepd_r."\t(".$mindtd_r.")\n".
        "arrival at epoch     ".$minepa_r."\t(".$mindta_r.")\n".
        "time mismatch =      ".$mintdiff_r."\n";
  $td=$mintd_r; $ta=$minta_r;
} # end if
#
@po=$Origin->posn_orb($td,"main");
$Thd=$po[2];
$omd=$po[3];
$rd=$po[4];                            # dist. from sun at departure
@po=$Destination->posn_orb($ta,"main");
$Tha=$po[2];
$oma=$po[3];
$ra=$po[4];                            # dist. from sun at arrival
#
if($task eq "-tran") {
  $beta_o=$ARGV[4];
  $Th=$Thd+$omd-$beta_o;
  $dth=$Tha+$oma-$Thd-$omd;
  $tt[0]=$transfer[0]->transfer_orbit($rd,$ra,$Th,$dth,$td,$beta_o,"main");
}
######################
#
#
######################
 if(($task eq "-scan") || ($task eq "-dv")) {
# if($task eq "-scan") {
# af::dumpvars("rd ra Th dth td beta",($rd,$ra,$Th,$dth,$td,$beta_o));
 @res=
   find_matching_trajectories($rd,$ra,$Thd,$Tha,$omd,$oma,$td,$ta,$n_or,$pl_time_diff);
 for($o=0;$o<2;++$o) {
   $beta_o=$res[$o]; $tdiff[$o]=$res[$o+2];
   if($tdiff[$o]!=$INVALID) {
     $Th=$Thd+$omd-$beta_o;
     $dth=$Tha+$oma-$Thd-$omd;
     $tt[$o]=$transfer[$o]->transfer_orbit($rd,$ra,$Th,$dth,$td,$beta_o,"main");
     print "o       tt ".$o."\t".$tt[$o]."\n";
   }
   else {$tt[$o]=$INVALID;}
   if($tt[$o]!=$INVALID) {
     @dvd=delta_v($Origin,$transfer[$o],$td,"main_dvd");
     @dva=delta_v($transfer[$o],$Destination,$ta,"main_dva");
     $dvda=sqrt($dvd[0]*$dvd[0]+$dvd[1]*$dvd[1]);
     $dvaa=sqrt($dva[0]*$dva[0]+$dva[1]*$dva[1]);
   }
   else {$dvda=20; $dvaa=20;}      
   print "delta v:\ndeparture: ".$dvda."\narrival:   ".$dvaa."\ntotal:     ".
         ($dvda+$dvaa)."\ntime mismatch: ".$tdiff[$o]."\n";
 }
}
#################
if(($task eq "-tran") || ($task eq "-scan") || ($task eq "-dv")) {
  for($o=0;$o<2;++$o) {
    if($tt[$o]<$LARGE) {
      plot_trajectory($td,$ta,$npt,$transfer[$o],$pl_transfer_traj[$o],1,"c");
      $op=$transfer[$o]->P();
      plot_trajectory($td,$td+$op,$npt,$transfer[$o],$pl_transfer_orbit[$o],0,"d");
    }
    else {
      plot_trajectory($td,$ta,-1,$transfer[$o],$pl_transfer_traj[$o],0,"d");    
      plot_trajectory($td,$ta,-1,$transfer[$o],$pl_transfer_orbit[$o],0,"d");    
    }
  }
}
# plot trajectories within interval of given departure and arrival times
plot_trajectory($td,$ta,$npt,$Origin,$pl_traj_of_origin,0,"a");
plot_trajectory($td,$ta,$npt,$Destination,$pl_traj_of_dest,0,"b");
plot_position($td,$Origin,$pl_posn_of_origin,"a");
plot_position($ta,$Destination,$pl_posn_of_dest,"b");
#
##################
#
# times from - to in days (epoch)
sub plot_trajectory {
 my $ME="plot trajectory";
 my ($tf,$tt,$npt,$object,$file,$pars,$CF)=@_;
 my ($m,$t,@pa,@po,$a,$e,$om);

 if($pars==1) {
   $a=$object->{_a};
   $e=$object->{_e};
   $om=$object->{_om};
   af::dumpvars("a e om",($a,$e,$om));
 }
 open F,">".$file;
 for($m=0;$m<$npt;++$m) {
   $t=$tf+$m*($tt-$tf)/($npt-1); # time in seconds
   @pa=$object->posn_abs($t,$ME." ".$CF);
   @po=$object->posn_orb($t,"main, a");
   print F $pa[0]."\t".$pa[1]."\t".$po[4]."\n";
 }
 close F;
}
#######################
#
#
#######################
sub plot_position {
 my $ME="plot_position";
 my ($t,$object,$file,$CF)=@_;   

 open F,">".$file;
 @pa=$object->posn_abs($t,$ME." ".$CF);
 @po=$object->posn_orb($t,"main, a");
 print F $pa[0]."\t".$pa[1]."\t".$po[4]."\n";
 close F; 
}
#######################
#
#
#######################
sub delta_v {
 my $ME="delta_v";   
 my ($orbit_a,$orbit_b,$t,$CF)=@_;
 my (@va,@vb,@dv);
    
 @va=$orbit_a->V($t,$ME." ".$CF);
 @vb=$orbit_b->V($t,$ME." ".$CF);
 @dv=($vb[0]-$va[0],$vb[1]-$va[1],$va[0],$va[1],$vb[0],$vb[1]);
 # af::dumpvars("d d b b a a",(@dv));
 @dv;
}
#######################
#
#
#######################
sub find_matching_trajectories {
 my ($rd,$ra,$Thd,$Tha,$omd,$oma,$td,$ta,$n_or,$pl_time_diff)=@_;
 my (@r,@res);
 my $n_or=360;
    
 @r=scan_orbit_orient($rd,$ra,$Thd,$Tha,$omd,$oma,$td,$ta,$n_or,$pl_time_diff); 
 @r=iterate_orbit_orient($rd,$ra,$Thd,$Tha,$omd,$oma,$td,$ta,@r);
 @res=($r[0],$r[0],$r[1],$r[1]); # ($beta_s,$beta_l,$tdiff,...)
# @res=Lambert_transfer($rd,$ra,$Thd,$Tha,$omd,$oma,$td,$ta,$n_or);
}
#######################
#
# using Lambert
#######################
sub Lambert_transfer {
 my $ME="Lambert_transfer";
 my ($rd,$ra,$Thd,$Tha,$omd,$oma,$td,$ta,$n_a)=@_;
# my ($m,$n,$a_from,$a_to,@a,$dth,$op,$c,$s,$al,$be,$ts,$tl,@tdiffs,@nlps);
# my (@tdiffl,@nlpl,$min,$mms,$mml,@res);
 my ($dth,@a,@ts,@tl,$ah,$tsh,$tlh,@mms,@mml,@as,@al,$bs,$bl,$a_h,$ts_h,$tl_h);
 my ($m); 
 
# af::dumpvars("dt",($ta-$td));
 $dth=$Tha+$oma-$Thd-$omd; 
 if($dth>2*$pi) {$dth-=2*$pi;}
 if($dth<0) {$dth+=2*$pi;}

# $a_h=$ah; $ts_h=$tsh; $tl_h=$tlh;
# @a=@$a_h; @ts=@$ts_h; @tl=@$tl_h;
#
# af::dumpvars("ta td",($ta,$td));
 ($ah,$tsh)=gen_times($dth,$rd,$ra,$n_a,"s");
 @a=@$ah; @ts=@$tsh;
 @mms=find_min_td($ta-$td,\@a,\@ts); # go through array
 af::dumpvars("al ah nlp",@mms);
 if($mms[0]!=$INVALID) {
   @as=iterate_min_td($ta-$td,$rd,$ra,$dth,$mms[0],$mms[1],$mms[2],'s'); 
   $bs=$Thd+$omd-find_perihelion($dth,$rd,$ra,$as[0]);
   if($bs<0) {$bs+=2*$pi;} if($bs>2*$pi) {$bs-=2*$pi;} 
#   af::dumpvars("_a tdiff al ah nlp",@as);
  # @res=($aa,$tdiff,$a_l,$a_h,$nlp);
 }
 else {$bs=$INVALID;}
 ($ah,$tlh)=gen_times($dth,$rd,$ra,$n_a,"l");
 @a=@$ah; @tl=@$tlh;    
 @mml=find_min_td($ta-$td,\@a,\@tl); # $a[low], $a[high],$nlp
 # af::dumpvars("al ah nlp",@mml); 
 if($mml[0]!=$INVALID) {
   @al=iterate_min_td($ta-$td,$rd,$ra,$dth,$mml[0],$mml[1],$mml[2],'l'); 
   $bl=$Thd+$omd-find_perihelion($dth,$rd,$ra,$al[0]);
   if($bl<0) {$bl+=2*$pi;} if($bl>2*$pi) {$bl-=2*$pi;}
#   af::dumpvars("__a tdiff al ah nlp bl",(@al,$bl));
   # af::dumpvars("_a",($al[0])) 
 }
 else {$bl=$INVALID;}
 open L,">alog";
 for($m=0;$m<scalar(@a);++$m) {print L $a[$m]."\t".$ts[$m]."\t".$tl[$m]."\n"}
 close L;
 @res=($bs,$bl,$as[1],$al[1]);
}

#######################
#
# quick-and-dirty
#
#######################
sub find_perihelion {
 my($dth,$ra,$rb,$a)=@_;
 my ($m,$n,$Th,$sa,$sb,@ea1,@ea2,@eb1,@eb2,$min,$el,$eh,$nit);
 my $n_th=200;
 my $max_it=1000;
    
 for($m=0;$m<$n_th;++$m) {
   $Th=2*$pi*$m/($n_th-1);   
   $sa=$ra*$ra*cos($Th)*cos($Th)          -4*$a*($ra-$a);
   $sb=$rb*$rb*cos($Th+$dth)*cos($Th+$dth)-4*$a*($rb-$a);
#   af::dumpvars("sa sb Th dth",($sa,$sb,$Th,$dth));
   if(($sa<0)||($sb<0)) 
     {$ea1[$m]=$ea2[$m]=$eb1[$m]=$eb2[$m]=$INVALID;} 
   else {
     $ea1[$m]=(-$ra*cos($Th     )+sqrt($sa))/(2*$a);
     $ea2[$m]=(-$ra*cos($Th     )-sqrt($sa))/(2*$a);
     $eb1[$m]=(-$rb*cos($Th+$dth)+sqrt($sb))/(2*$a);
     $eb2[$m]=(-$rb*cos($Th+$dth)-sqrt($sb))/(2*$a);
   }
 }
 af::dumpvars("ARRAY","ea1",@ea1);
 af::dumpvars("ARRAY","ea2",@ea2);
 af::dumpvars("ARRAY","eb1",@eb1);
 af::dumpvars("ARRAY","eb2",@eb2);
 for($min=$LARGE,$mm=-1,$m=0;$m<$n_th;++$m) {
   if(($ea1[$m]!=$INVALID)&&($ea1[$m]>=0)&&($ea1[$m]<=1)) {
     if((abs($ea1[$m]-$eb1[$m])<$min)&&($eb1[$m]>=0)&&($eb1[$m]<=1))
       {$min=abs($ea1[$m]-$eb1[$m]); $mm=$m;}
     if((abs($ea1[$m]-$eb2[$m])<$min)&&($eb2[$m]>=0)&&($eb2[$m]<=1))
       {$min=abs($ea1[$m]-$eb2[$m]); $mm=$m;}
   }
   if(($ea2[$m]!=$INVALID)&&($ea2[$m]>=0)&&($ea2[$m]<=1)) {
     if((abs($ea2[$m]-$eb1[$m])<$min)&&($eb1[$m]>=0)&&($eb1[$m]<=1))
       {$min=abs($ea2[$m]-$eb1[$m]); $mm=$m;}
     if((abs($ea2[$m]-$eb2[$m])<$min)&&($eb2[$m]>=0)&&($eb2[$m]<=1))
       {$min=abs($ea2[$m]-$eb2[$m]); $mm=$m;}
   }
 }
# af::dumpvars("mm ea1 ea2 eb1 eb2",($mm,$ea1[$mm],$ea2[$mm],$eb1[$mm],$eb2[$mm]));
# if(($mm!=0)&&($mm!=$n_th-1)) {$ml=$mm-1; $mh=$mm+1;}
# else {if($mm==0) {$ml=$n_th-1;$mh=1;} else {$ml=$n_th-2; $mh=0;}}     
# do {
#   $ee=($el+$eh)/2;
#   $Th=acos(($a*(1-$ee*$ee)-$ra)/($ra*$ee));
#   $er=$a*(1-$ee*$ee)/(1+$ee*cos($Th+$dth))-$rb;
#   if($er<0) {$el=$ee;} else {$eh=$ee;} 
#   ++$nit;
# } while(($er>1.0e-9)&&($nit<$max_it));
# if($nit==$max_it) {
#     af::dumpvars("ARRAY","mm",@mm);
#     af::dumpvars("dth ra rb a el eh",($dth,$ra,$rb,$a,$el,$eh));     
#     die "bisection does not converge";}
# $th=acos(($a*(1-$ee*$ee)-$ra)/($ra*$ee)); 
 $Th=2*$pi*$mm/($n_th-1);
}
#######################
sub comp {
 my ($ea1,$ea2,$eb1,$eb2,$min)=@_;
 my ($a,$b,$c,$d);
 
 $a=$b=$c=$d=$LARGE;
 if($ea1 != $INVALID) {$a=abs($ea1-$eb1); $b=abs($ea1-$eb2);}
 if($ea2 != $INVALID) {$c=abs($ea2-$eb1); $d=abs($ea2-$eb2);}
 
 $min; 
}

#######################
#
# quick-and-dirty
#
#######################
sub find_perihelion_e {
 my($dth,$ra,$rb,$a)=@_;
 my ($m,$n,$e,$Th,$cTh,@mm,@zc,$el,$eh,$nit);
 my $n_e=1000;
 my $max_it=1000;
    
 for($m=0;$m<$n_e;++$m) {
   $e=1.0e-6+$m/($n_e-1);
   $cTh=($a*(1-$e*$e)-$ra)/($ra*$e);
#   $Th=acos($cTh);
   if(abs($cTh)<1) {
     $mm[$m]=$a*(1-$e*$e)/(1+$e*cos($cTh+$dth))-$rb;
   }
   else {$mm[$m]=$INVALID;}
 }
 for($m=1,$n=0;$m<$n_e-1;++$m) {
   if(($mm[$m]*$mm[$m-1]<0)&&($mm[$m]!=$INVALID)&&($mm[$m-1]!=$INVALID))
     {$zc[$n]=$m; $n++;} 
 } 
 if($n==0) {
     af::dumpvars("ARRAY","mm",@mm);
     af::dumpvars("dth ra rb a",($dth,$ra,$rb,$a));
     die "no zero crossing";
 }
 if($n>1) {
     af::dumpvars("ARRAY","mm",@mm);
     af::dumpvars("n dth ra rb a",($n,$dth,$ra,$rb,$a));     
     af::dumpvars("zc0 zc1",($zc[0],$zc[1]));
     die "more than 1 zero crossing";}
 if($mm[$zc[0]-1]<$mm[$zc[0]]) {$el=($zc[0]-1)/($n_e-1); $eh=$zc[0]/($n_e-1);}
 else {$eh=($zc[0]-1)/($n_e-1); $el=$zc[0]/($n_e-1);}
 $nit=0;
 do {
   $ee=($el+$eh)/2;
   $Th=acos(($a*(1-$ee*$ee)-$ra)/($ra*$ee));
   $er=$a*(1-$ee*$ee)/(1+$ee*cos($Th+$dth))-$rb;
   if($er<0) {$el=$ee;} else {$eh=$ee;} 
   ++$nit;
 } while(($er>1.0e-9)&&($nit<$max_it));
 if($nit==$max_it) {
     af::dumpvars("ARRAY","mm",@mm);
     af::dumpvars("dth ra rb a el eh",($dth,$ra,$rb,$a,$el,$eh));     
     die "bisection does not converge";}
 $th=acos(($a*(1-$ee*$ee)-$ra)/($ra*$ee));
}

#######################
#
#
#######################
sub gen_times {
 my($dth,$ra,$rb,$n_a,$sl)=@_;
 my ($m,@a,$op,@al,@be,@t,@res);
 my ($a_from,$a_to,$c,$s);
    
 $c=sqrt($ra*$ra+$rb*$rb-2*$ra*$rb*cos($dth));
 $s=0.5*($ra+$rb+$c);
 # a must be > s/2 and > (s-c)/2 to make acos, below, real-valued
 $a_from=1.000000001*$s/2; # (s-c)/2 is less, anyway
 $a_to=10*$AU;
 for($m=0;$m<$n_a;++$m) {
   $a[$m]=$a_from+($a_to-$a_from)*$m/($n_a-1);
   $t[$m]=t_a($a[$m],$ra,$rb,$dth,$sl);
 }
# af::dumpvars("aa ra rb dth sl",($aa,$ra,$rb,$dth,$sl)); 
 
 @res=(\@a,\@t);
}
#######################
#
#
#######################
sub iterate_min_td {
 my ($dt,$ra,$rb,$dth,$a_l,$a_h,$nlp,$sl)=@_;
 my ($aa,$nit,$t,$tdiff);
 my $max_it=1000;

   $aa=$a_l;
   $t=t_a($aa,$ra,$rb,$dth,$sl);
#   print "-- ".$aa."\t".($dt-$t)."\n";
   $aa=$a_h;
   $t=t_a($aa,$ra,$rb,$dth,$sl);  
#   print "-- ".$aa."\t".($dt-$t)."\n";
# af::dumpvars("aa ra rb dt dth sl",($aa,$ra,$rb,$dt,$dth,$sl)); 
 $nit=0;    
 do {
   $aa=($a_l+$a_h)/2;
   $t=t_a($aa,$ra,$rb,$dth,$sl);
   $tdiff=$dt-($t+$nlp*$op);
#   af::dumpvars("a al ah tdiff",($aa,$a_l,$a_h,$tdiff));
   if($tdiff<0) {$a_l=$aa} else {$a_h=$aa;}
   $nit+=1;
 } while((abs($tdiff)>1.0e-12*abs($dt))&&($nit<$max_it));
 # af::dumpvars("_aa _al _ah nit tdiff",($aa,$a_l,$a_h,$nit,$tdiff)); 
 if($nit==$max_it) {$tdiff=$INVALID;}
 @res=($aa,$tdiff,$a_l,$a_h,$nlp);
}

################
#
################
sub t_a {
 my ($a,$ra,$rb,$dth,$sl)=@_;
 my ($op,$c,$s,$al,$be,$t,$tdiff);
  
# af::dumpvars("a ra rb dth sl",($a,$ra,$rb,$dth,$sl));
 $op=2*$pi*sqrt($a*$a*$a)/sqrt($mu);
 $c=sqrt($ra*$ra+$rb*$rb-2*$ra*$rb*cos($dth));
 $s=0.5*($ra+$rb+$c);
 $al=acos(1-$s/$a); if(($al<0)||($al>$pi)) {die "";}
 $be=acos(1-($s-$c)/$a); if(($be<0)||($be>$pi)) {die "";}
# $al=0.5;$be=0.5;
 if($sl eq "s") {$t=($op/(2*$pi))*(($al-sin($al))-($be-sin($be)));}
 else           {$t=$op-($op/(2*$pi))*(($al-sin($al))+($be-sin($be)));}
 $t; # diff=$dt-($t+$nlp*$op);
# af::dumpvars("a al ah tdiff",($aa,$a_l,$a_h,$tdiff));
}
#######################
#
#
#######################
sub find_min_td {
 my ($dt,$a_h,$t_h)=@_;
 my ($m,$n,@a,@t,$n_a,$op,@tdiff,@nlp,@res);

 @a=@$a_h;
 @t=@$t_h;    
 $n_a=scalar(@a);
 for($m=0;$m<$n_a;++$m) {
   $op=2*$pi*sqrt($a[$m]*$a[$m]*$a[$m])/sqrt($mu);   
   # find no. of loops
   for($tdiff[$m]=$LARGE,$n=0;$n<=$nloop;++$n) {
     if(abs($dt-($t[$m]+$n*$op))<abs($tdiff[$m])) {
       $tdiff[$m]=$dt-($t[$m]+$n*$op);
       # af::dumpvars("dt t tdiff",($dt,$t[$m],$tdiff[$m]));
       $nlp[$m]=$n;
     }
   }
 }
 #
# af::dumpvars("ARRAY","t",@t); 
 for($mm=-1,$min=$LARGE,$m=0;$m<$n_a;++$m) {
   if(abs($tdiff[$m])<$min) {$min=abs($tdiff[$m]); $mm=$m;}
 }
 if(($mm<1)||($mm>=$n_a-1)||($tdiff[$mm-1]*$tdiff[$mm+1]>=0)) {
   $a[$mm-1]=$a[$mm]=$a[$mm+1]=$INVALID; 
 }
 if(($nlp[$mm-1]!=$nlp[$mm+1])||($nlp[$mm-1]!=$nlp[$mm]))
   {die; "nlp does not agree";}  
 if($tdiff[$mm-1]<$tdiff[$mm+1]) {@res=($a[$mm-1],$a[$mm+1],$nlp[$mm]);}
 else                            {@res=($a[$mm+1],$a[$mm-1],$nlp[$mm]);}
# af::dumpvars("mm tdiff-1 tdiff+1 a-1 a+1",
#	      ($mm,$tdiff[$mm-1],$tdiff[$mm+1],$a[$mm-1],$a[$mm+1]));
 
 @res;
}

#######################
# Scanning orbit orientations
#
#######################
sub scan_orbit_orient {
 my $ME="scan_orbit_orient";
 my ($rd,$ra,$Thd,$Tha,$omd,$oma,$td,$ta,$n_or,$pl_time_diff)=@_;
 my ($m,$n,@beta,$Th,$dth,$a,$mu,$op,@tdiff,@nlp,$max,$min,$mm,@res);
  
# open L, ">logfile";
 for($m=0;$m<$n_or;++$m) {
   $beta[$m]=2*$pi*$m/($n_or-1);
   $Th=$Thd+$omd-$beta[$m];
   $dth=$Tha+$oma-$Thd-$omd;
   $tt=$transfer[2]->transfer_orbit($rd,$ra,$Th,$dth,$td,$beta[$m],$ME);
   # try direct trajectory, and ones looping around the sun a few times ...
   if(abs($tt)<$LARGE) { # $tt=$LARGE indicates subelliptic or hyperbolic traj.
     $a=$transfer[2]->{_a};
     $mu=$transfer[2]->{_mu}; # mu is global, but using object's own doesn't hurt
     $op=2.0*$pi*sqrt($a*$a*$a)/sqrt($mu);
#     print L $tt."\t".$op."\t"."\n";
     # multiple loops on transfer orbit
     for($tdiff[$m]=$LARGE,$n=0;$n<=$nloop;++$n) {
       if(abs(($ta-$td)-($tt+$n*$op))<abs($tdiff[$m])) {
	 $tdiff[$m]=($ta-$td)-($tt+$n*$op);
         $nlp[$m]=$n;
       }
     }
#     $tdiff[$m]=($ta-$td)-$tt;     #  here!
   }
   else {$tdiff[$m]=$LARGE; $nlp[$m]=0;}
 }
 $max=-$LARGE;
 $min=$LARGE;
 # find min and max of all points that are not marked as invalid
 for($m=0;$m<$n_or;++$m) {
   if(abs($tdiff[$m])<$LARGE/2) {
     if($tdiff[$m]<$min) {$min=$tdiff[$m];}
     if($tdiff[$m]>$max) {$max=$tdiff[$m];}
   }
 }
 # limit invalid points to min-max range (to better plot it)
 for($m=0;$m<$n_or;++$m) {
   if($tdiff[$m]<$min) {$tdiff[$m]=$min;}
   if($tdiff[$m]>$max) {$tdiff[$m]=$max;}
 }
 # seek minimum time diff
 for($min=$LARGE,$m=0;$m<$n_or;++$m) {
   if(abs($tdiff[$m])<$min) {$min=abs($tdiff[$m]); $mm=$m;}
 }
 # generate plot of time mismatch 
 if($pl_time_diff ne "") {
   open T,">".$pl_time_diff;
   for($m=0;$m<$n_or;++$m) {
     print T $beta[$m]."\t".$tdiff[$m]."\n";
   }  
   close T;
 }
 # results array
 if(($mm>0)&&($mm<$n_or-1))
   {@res=($beta[$mm],$tdiff[$mm],$beta[$mm-1],$beta[$mm+1],$nlp[$mm]);}
 if($mm==0) {@res=($beta[0],$tdiff[$mm],$beta[$n_or],$beta[1],$nlp[$mm]);} 
 if($mm==$n_or-1)
   {@res=($beta[$mm],$tdiff[$mm],$beta[$n_or-2],$beta[0],$nlp[$mm]);}

# $res[3]=$tdiff[$m]; # =($ta-$td)-$tt;
# close L;
 @res; # res must have 5 el., that is expected by iterate_...
}

#######################
#
#
#######################
sub iterate_orbit_orient {
 my $ME="iterate_orbit_orient";
 my ($rd,$ra,$Thd,$Tha,$omd,$oma,$td,$ta,@res)=@_;
 my ($beta,$Th,$dth,$tt,$tdiff,$tda,$tdb,$bl,$bh,$nit,$a,$mu,$op,$nlp);
 my $max_it=1000;
   
 $bl=$res[2]; $bh=$res[3]; $nlp=$res[4];
 $Th=$Thd+$omd-$bl;
 $dth=$Tha+$oma-$Thd-$omd;
 $tt=$transfer[2]->transfer_orbit($rd,$ra,$Th,$dth,$td,$bl,$ME);
 $tda=($ta-$td)-$tt;
 $Th=$Thd+$omd-$bh;
 $tt=$transfer[2]->transfer_orbit($rd,$ra,$Th,$dth,$td,$bh,$ME);
 $tdb=($ta-$td)-$tt;
 if($tda>$tdb) {$bl=$res[3]; $bh=$res[2];}
 $nit=0;
# print "nlp = ".$nlp."\n";
 do {
   $beta=($bl+$bh)/2;
   $Th=$Thd+$omd-$beta;
   $dth=$Tha+$oma-$Thd-$omd;
   $tt=$transfer[2]->transfer_orbit($rd,$ra,$Th,$dth,$td,$beta,$ME);
   if($tt==$INVALID) {$tdiff=$INVALID; $nit=$max_it;} # last does not work
   else { # apparently terminates loops in main, as well
     $a=$transfer[2]->{_a};
     $mu=$transfer[2]->{_mu};
     $op=2*$pi*sqrt($a*$a*$a)/sqrt($mu);
     $tdiff=($ta-$td)-($tt+$nlp*$op);
#   $tdiff=($ta-$td)-$tt;
     if($tdiff<0) {$bl=$beta} else {$bh=$beta;}
     $nit+=1;
   }
 } while((abs($tdiff)>1.0e-12*abs($ta-$td))&&($nit<$max_it));
 if($nit==$max_it) {$tdiff=$INVALID;}
 @res=($beta,$tdiff,$bl,$bh,$nlp);
 @res;
}

###############
#
#
#
##############
sub sq {
  my $x=$_[0];
  
  my $y=$x*$x;
}