#!/usr/bin/perl

# 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 Math::Trig;
use POSIX; # floor()
use Time::HiRes qw (usleep);
# use Math::Gsl::Sf;
# use Math::Gsl::Sf qw(:Gamma);
# use Math::Complex;

use Orbit;
use aux;
# use array;
# misc. parameters
$pi=2*asin(1.0);
$deg=$pi/180.0;
$rad2deg=180/$pi;
# $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;
$INF=1.0e30;
$ORBIT_VAL=1;
$ORBIT_INV=-1;
# reference date: 13.04.2029
# 

# Sun
$mu=1.32712440018e11;

# JPL:  http://ssd.jpl.nasa.gov/txt/p_elem_t1.txt
# Origin (ex. earth)
# $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
##########################################
#
$mu=       0;          # gravitational parameter mu
$n_orb=    0;          # number of orbits
$n_t_orb=  0;          # number of transit orbits
#                      #
@orbits=   qw();       # array of indices to orbits
#                      #
@O_name=   qw();       # names of orbits
@O_a=      qw();       # a of orbits 
@O_dadt=   qw();       # time derivative of a of orbits 
@O_e=      qw();       # e of orbits
@O_dedt=   qw();       # time derivative of e of orbits
@O_i=      qw();       # i of orbits
@O_didt=   qw();       # time derivative of i of orbits
@O_an=     qw();       # ascending node of orbits
@O_dandt=  qw();       # time derivative of ascending node of orbits
@O_ep=     qw();       # epoch for which L0 or M0 is given !!!!!!!! IN SECONDS
@O_om=     qw();       # argument of periapsis from node
@O_omb=    qw();       # longitude of periapsis
@O_domdt=  qw();       # time der. of argument of periapsis from node
@O_L0=     qw();       # argument of planet from node at epoch
@O_M0=     qw();       # mean anom. of planet from periapsis at epoch
@O_tp=     qw();       # time of periapsis (seconds)

@O_pos_s=  qw();       # start pos. file
@O_pos_e=  qw();       # end pos. file
@O_tra=    qw();       # trajectory file
@O_orb=    qw();       # orbit file
@O_ttlog=  qw();       # log file for travel times

#                      #
@transit=  qw();       # array of indices to transit orbits
#                      #
@O_from=   qw();       # for transfer traj
@O_to=     qw();       # for transfer traj
@O_nlp=    qw();       # for transfer traj
#
#
$t_s=      0;
$t_e=      0;
#
$np_orb=   0;          # no. of points in orbit plots
$np_tra=   0;          # no. of points in trajectory plots
@pl_pos_s= qw();       # files with position at start time
@pl_pos_e= qw();       # files with position at start time
@pl_orb=   qw();       
@pl_tra=   qw();       # trajectories on resp. orbits between start, end times
@pl_orbt=  qw();
@pl_trat=  qw();
$todo=   "";
$o_t=     0;           # unit of time in output
$o_l=     0;           # unit of length in output
$o_v=     0;           # unit of velocity in output
$o_vdt=   0;           # delta t for drawing velocity vectors
$o_dvl="no";           # draw velocity vectors on positions, yes or no

##############
#
# hash of all parameters to be read from parameter file
# keys are keys in par. file, values are scalar or array handles
#
#############
%parameters=(
	   "orbits[]"                       => \@orbits,
	   "transit[]"                      => \@transit,
	   "grav. param. mu"                => \$mu,
   	   "todo"                           => \$todo,
	   "output unit of time"            => \$o_t,
	   "output unit of length"          => \$o_l,
	   "output unit of velocity"        => \$o_v,  
	   "output dt of velocity vectors"  => \$o_vdt,  
	   "output, draw velocity vectors"  => \$o_dvl,
	   "orbit _O_, name"                => \@O_name,
	   "orbit _O_, a"                   => \@O_a,
	   "orbit _O_, da/dt"               => \@O_dadt,
	   "orbit _O_, e"                   => \@O_e,
	   "orbit _O_, de/dt"               => \@O_dedt,
	   "orbit _O_, i"                   => \@O_i,
	   "orbit _O_, di/dt"               => \@O_didt,
	   "orbit _O_, asc. node"           => \@O_an,
	   "orbit _O_, d_an/dt"             => \@O_dandt,
	   "orbit _O_, epoch"               => \@O_ep,
	   "orbit _O_, arg. of periapsis"   => \@O_om,
	   "orbit _O_, long. of periapsis"  => \@O_omb,
	   "orbit _O_, d_om/dt"             => \@O_domdt,
	   "orbit _O_, L_0 at ref. epoch"   => \@O_L0,
	   "orbit _O_, M at ref. epoch"     => \@O_M0,
	   "orbit _O_, time of periapsis"   => \@O_tp,
	   "orbit _O_, start posn file"     => \@O_pos_s,
	   "orbit _O_, end posn file"       => \@O_pos_e,
	   "orbit _O_, traj file"           => \@O_tra,
	   "orbit _O_, orbit file"          => \@O_orb,	   
	   "orbit _O_, rel start posn file" => \@O_rpos_s,
	   "orbit _O_, rel end posn file"   => \@O_rpos_e,
	   "orbit _O_, rel traj file"       => \@O_rtra,
	   "orbit _O_, rel orbit file"      => \@O_rorb,
	   "orbit _O_, Delta v file"        => \@O_deltav,
	   "orbit _O_, tt log file"         => \@O_ttlog,	     
           "start time"                     => \$t_s,
           "end time"                       => \$t_e,
	   "scan, start time from"          => \$tsf,
	   "scan, start time to"            => \$tst,
	   "scan, start time step"          => \$tss,
	   "scan, end time from"            => \$tef,
	   "scan, end time to"              => \$tet,
	   "scan, end time step"            => \$tes,	     
	   "no. of points for orbits"       => \$np_orb,	   
           "no. of points for trajectories" => \$np_tra,
	   "orbit _O_, from"                => \@O_from,
	   "orbit _O_, to"                  => \@O_to,	     
	   "orbit _O_, no. of loops"        => \@O_nlp,
	   "maximum Delta v"                => \$maxdv,  
	   "diagnostics"                    => \$diagnostics
	   );
#########################################################
#

# $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";
########################

$parfile=$ARGV[0];
get_pars($parfile,\%parameters);
dump_pars("pardump",%parameters);
process_pars(\%parameters);
$ph=$parameters{"start time"};
print "departure: epoch ".($t_s/$day)."\t".aux::epoch_to_date($t_s/$day)."\n";
$ph=$parameters{"end time"};
print "arrival:   epoch ".($t_e/$day)."\t".aux::epoch_to_date($t_e/$day)."\n";
# $n_orb_h=$parameters{"number of orbits"};
# $n_orb=$$n_orb_h;
# die "mu = ".$mu;
# die "n_orb = ".$n_orb;
# $nlpmin=$INF;
# $nlpmax=-$INF;
for($m=0;$m<$n_orb;++$m) {
#  print "initializing orbit ".$O_name[$m]."\n";
#  print "i : ".$O_i[$m].", ".$O_didt[$m]."\n";
  $Orb[$m]=Orbit->new($O_a[$m], $O_dadt[$m],
		      $O_e[$m], $O_dedt[$m],
		      $O_i[$m], $O_didt[$m],
		      $O_an[$m],$O_dandt[$m],
		      $O_ep[$m],
		      $O_om[$m],$O_omb[$m],$O_domdt[$m],
		      $O_L0[$m],$O_M0[$m],$O_tp[$m],
		      $mu,1,
		      $O_name[$m],
		      $O_pos_s[$m],$O_pos_e[$m],$O_tra[$m],$O_orb[$m],
		      $O_rpos_s[$m],$O_rpos_e[$m],$O_rtra[$m],$O_rorb[$m],
		      $O_deltav[$m],
		      $O_ttlog[$m],
		      $O_from[$m],$O_to[$m],$O_nlp[$m],
		      $diagnostics);
  if(($Orb[$m]->{_from} ne "undef") && ($Orb[$m]->{_to} ne "undef") &&
     ($Orb[$m]->{_nlp} ne "undef")) { # this is a transfer orbit
    $group_name=$Orb[$m]->{_from}."--".$Orb[$m]->{_to};
  }
  else {$group_name=$Orb[$m]->{_name};}
#  if(exists $groups{$group_name}) {$groups{$group_name}->add($m);}
#  else {$groups{$group_name}=array->new($m);}
#  $Orb[$m]->{_group} = $groups{$group_name};

# print $O_name[$m].": ".($Orb[$m]->{_om}/$deg).", ".($Orb[$m]->{_M0}/$deg).
#        ", ".($Orb[$m]->{_ep}/$day).", ".($Orb[$m]->{_tp}/$day)."\n"; 
#
}
if($todo =~ m/orbital data/) {
  $day=3600*24;
  $AU=1.49597870691e8;
  for($m=0;$m<$n_orb;++$m) {
    print "perihelion of ".$O_name[$m]." at epoch ".($Orb[$m]->tp()/$day).
        " (".aux::sec_to_date($Orb[$m]->tp()).")\n";  
  }
  for($m=0;$m<$n_orb;++$m) {
    print "At epoch ".($t_s/$day)." (".aux::sec_to_date($t_s)."):\n";
    @rMET=$Orb[$m]->posn_pol_rel($t_s,"main");  
    print $O_name[$m].": M0 = ".($rMET[1]*$rad2deg)." deg, ".
          "r = ".($rMET[0]/$AU)." AU\n";
  }
#  plot_rel_trajectory($pl_pos1,1,$Origin,$t_s,$t_s,"main");           # position
#  plot_rel_trajectory($pl_tra1,$np_tra,$Origin,$t_s,$t_e,"main");     # trajectory
#  plot_rel_trajectory($pl_orb1,$np_orb,$Origin,$t_s,$t_s+$Origin->P(),"main"); # orbit
}
#####################################
if($todo =~ m/scan Delta v/) {
#  $tss=($tst-$tsf)/$nts;
#  $tes=($tet-$tef)/$nte;
  for($o=0,$m=0;$m<$n_orb;++$m) { # loop over all orbits, see which ones are transfer
    if(($Orb[$m]->from() ne "undef") && ($Orb[$m]->to() ne "undef") && 
       ($Orb[$m]->nlp() ne "undef")) {
      for($fi=-1,$ti=-1,$n=0;$n<$n_orb;++$n) {
        if($Orb[$m]->from() eq $Orb[$n]->name()) {$fia[$o]=$n;}
        if($Orb[$m]->to() eq $Orb[$n]->name()) {$tia[$o]=$n;}
      } # end for
      $Oa[$o]=$m;
#      print "----- ".$m.": ".$Orb[$m]->{_deltav}."\n";
      open $fh[$o],">".$Orb[$m]->{_deltav};
      $o+=1;
    } # end if
  } # end for
#  print "tst-tsf, tss = ".($tst-$tsf).", ".$tss."\n";
  $mindv=1.0e30;
  open L,">scan.log";
  $lp=0;
  for($ts=$tsf;$ts<=$tst;$ts+=$tss) {
    print "ts = ".($ts/$day).", ".aux::epoch_to_date($ts/$day)."\n";
    for($te=$tef;$te<=$tet;$te+=$tes,$lp+=1) {
  #    print "te = ".$te."\n";
      if($lp%4 == 0) {print "\r-";}
      if($lp%4 == 1) {print "\r/";}
      if($lp%4 == 2) {print "\r|";}
      if($lp%4 == 3) {print "\r\\";}
      for($o=0;$o<scalar(@Oa);++$o) { # loop over all orbits, see which ones are transfer 
        $fhc=$fh[$o];
	@from=$Orb[$fia[$o]]->posn_cart_abs($ts,"main");
        @to=  $Orb[$tia[$o]]->posn_cart_abs($te,"main");
        @v_from=$Orb[$fia[$o]]->vel_cart_abs($ts,"main");
        @v_to=  $Orb[$tia[$o]]->vel_cart_abs($te,"main");      
        $Orb[$Oa[$o]]->{_ntheta}=60;
	do {
  	  $Orb[$Oa[$o]]->set_endpoints($ts,$te,@from,@to);
          ($a,$e,$Tm,$op)=$Orb[$Oa[$o]]->transfer_trajectory($O_nlp[$Oa[$o]]);
          if($a==0) {
	    $Orb[$Oa[$o]]->{_ntheta}*=5;
	  #  print "N*5, ";
	  }
	} while(($a==0) && ($Orb[$Oa[$o]]->{_ntheta} < 10000));
	if($Orb[$Oa[$o]]->{_ntheta}>=10000) 
	  {$Orb[$Oa[$o]]->{_complete}= $ORBIT_INV;}
        if($Orb[$Oa[$o]]->{_complete}==$ORBIT_VAL) {
          @vt_from=$Orb[$Oa[$o]]->vel_cart_abs($ts,"main");
          @vt_to=  $Orb[$Oa[$o]]->vel_cart_abs($te,"main");      
          @dvf=aux::vdiff(@v_from,@vt_from);
          @dvt=aux::vdiff(@v_to,@vt_to);
	  $advf=aux::vabs(@dvf);
	  $advt=aux::vabs(@dvt);
	  if(($advf+$advt)<$mindv) {
	    $mindv=$advf+$advt;
	    $mindvf=$advf;
	    $mindvt=$advt;
	    $t_s=$ts;
	    $t_e=$te;
	    @mv_from=@v_from;
	    @mvt_from=@vt_from;
	    @mv_to=@vt_to;
	    @mvt_to=@vt_to;
	    print L "mindv, t_s, t_e = ".$mindv.", ".($t_s/$day).
		    ", ".($t_e/$day)."\n";
	  }
	  $sdv=$advf+$advt;
	  if($sdv>$maxdv) {$sdv=$maxdv;}
	  # ts, te, Dvs+Dve, Dvs, Dve, a, e
          if(0) {
	    print $fhc ($ts/$o_t)."\t".($te/$o_t)."\t".($sdv/$o_v)."\t".
	          $advf."\t".$advt."\t".($Orb[$Oa[$o]]->{_a}/$o_l)."\t".
	          $Orb[$Oa[$o]]->{_e}."\n";
	  }
	  else {
	    print $fhc int($ts/$o_t+0.5)."\t".int($te/$o_t+0.5)."\t".
	               ($sdv/$o_v)."\n";
	  }
	}
	else {
	  print $fhc int($ts/$o_t+0.5)."\t".int($te/$o_t+0.5)."\t".(1000)."\n";
	}
#	if($Orb[$Oa[$o]]->valid() == 0) {print "invalid transfer orbit\n";}
      } # end for  
    }
    for($o=0;$o<scalar(@fh);++$o) {$fhc=$fh[$o]; print $fhc "\n";}
  }
  for($o=0;$o<scalar(@fh);++$o) {close $fh[$o];}
  print "v_from = (".$mv_from[0].", ".$mv_from[1].", ".$mv_from[2].")\n";
  print "vt_from = (".$mvt_from[0].", ".$mvt_from[1].", ".$mvt_from[2].")\n";
  print "v_to = (".$mv_to[0].", ".$mv_to[1].", ".$mv_to[2].")\n";
  print "vt_to = (".$mvt_to[0].", ".$mvt_to[1].", ".$mvt_to[2].")\n";
#  
  print "minimum Delta v = ".$mindv." for departure at ".($t_s/$day)." (".
        aux::epoch_to_date($t_s/$day).") with Delta v = ".
        $mindvf.", arrival at ".($t_e/$day)." (".aux::epoch_to_date($t_e/$day).
        ") with Delta v = ".$mindvt."\n";
  close L;
} # end if(todo)
##################################
if($todo =~ m/transfer trajectory/) {
  for($m=0;$m<$n_orb;++$m) { # loop over all orbits, see which ones are transfer
    if(($Orb[$m]->from() ne "undef") && ($Orb[$m]->to() ne "undef") && 
       ($Orb[$m]->nlp() ne "undef")) {
      for($fi=-1,$ti=-1,$n=0;$n<$n_orb;++$n) {
        if($Orb[$m]->from() eq $Orb[$n]->name()) {$fi=$n;}
        if($Orb[$m]->to() eq $Orb[$n]->name()) {$ti=$n;}
      } # end for
      if(($fi==-1)||($ti==-1)) {die "from ".$Orb[$m]->from()." to ".$Orb[$m]->to();}
      @from=$Orb[$fi]->posn_cart_abs($t_s,"main");
      @to=  $Orb[$ti]->posn_cart_abs($t_e,"main");
      @v_from=$Orb[$fi]->vel_cart_abs($t_s,"main");
      @v_to=  $Orb[$ti]->vel_cart_abs($t_e,"main");      
      $Orb[$m]->{_ntheta}=3000; ########### here     
      $Orb[$m]->set_endpoints($t_s,$t_e,@from,@to);
      ($a,$e,$Tm,$op)=$Orb[$m]->transfer_trajectory($O_nlp[$m]);
      if($Orb[$m]->{_complete}==$ORBIT_VAL) {
	@tfrom=$Orb[$m]->posn_cart_abs($t_s,"main");
	@tto=$Orb[$m]->posn_cart_abs($t_e,"main");
	@vt_from=$Orb[$m]->vel_cart_abs($t_s,"main");
        @vt_to=  $Orb[$m]->vel_cart_abs($t_e,"main");      
        @dvf=aux::vdiff(@v_from,@vt_from);
        @dvt=aux::vdiff(@v_to,@vt_to);
	$advf=aux::vabs(@dvf);
	$advt=aux::vabs(@dvt);
	@dposf=aux::vdiff(@from,@tfrom);
	@dpost=aux::vdiff(@to,@tto);
	print "a, e = ".$a.", ".$e."\n";
        print "dpos at dept.: (".$dposf[0].", ".$dposf[0].", ".$dposf[0].")\n";
	print "dpos at arrv.: (".$dpost[0].", ".$dpost[0].", ".$dpost[0].")\n";
        if(aux::vabs(@dposf)>0.1) {die "> 100 m difference in dept. positions";}
        if(aux::vabs(@dpost)>0.1) {die "> 100 m difference in arrv. positions";}
	print "v_from = (".$v_from[0].", ".$v_from[1].", ".$v_from[2].")\n";
        print "vt_from = (".$vt_from[0].", ".$vt_from[1].", ".$vt_from[2].")\n";
        print "v_to = (".$v_to[0].", ".$v_to[1].", ".$v_to[2].")\n";
        print "vt_to = (".$vt_to[0].", ".$vt_to[1].", ".$vt_to[2].")\n";
        print "Delta v for ".$Orb[$m]->{_name}.":\n".
              "departure: ".$advf."\n".
              "arrival:   ".$advt."\n".
              "total:     ".($advf+$advt)."\n".
	      "a, e, Tm, op = ".$a.", ".$e.", ".$Tm.", ".$op."\n";
      }
      if($Orb[$m]->valid() == 0) {print "invalid transfer orbit\n";}
    } # end if(from, to, nlp)
  } # end for
} # end if(todo)
##################################
if($todo =~ m/plot trajectories/) {
#  print "start time = ".$$ph."\t".$t_s."\n";  
  for($m=0;$m<$n_orb;++$m) {    
#    print $Orb[$m]->{_name}."complete ".$Orb[$m]->{_complete}."\n";
    if($Orb[$m]->{_complete} == $ORBIT_VAL) { # only if orbital parameters are defined
      plot_abs_trajectory("_PS",$o_dvl,1,$Orb[$m],$t_s,$t_s,"main"); # position
      plot_abs_trajectory("_PE",$o_dvl,1,$Orb[$m],$t_e,$t_e,"main"); # position
      plot_abs_trajectory("_TR","no",$np_tra,$Orb[$m],$t_s,$t_e,"main"); # trajectory
      plot_abs_trajectory("_OR","no",$np_orb,$Orb[$m],$t_s,$t_s+$Orb[$m]->P(),"main"); # orbit
    }
  }
}
if($todo =~ m/plot relative trajectories/) {
  for($m=0;$m<$n_orb;++$m) {    
    if($Orb[$m]->{_complete} == $ORBIT_VAL) { # only if orbital parameters are defined
#      print "--- ".$Orb[$m]->{_RTR}."\n";	
      plot_rel_trajectory("_RPS",1,$Orb[$m],$t_s,$t_s,"main"); # position
      plot_rel_trajectory("_RPE",1,$Orb[$m],$t_e,$t_e,"main"); # position
      plot_rel_trajectory("_RTR",$np_tra,$Orb[$m],$t_s,$t_e,"main"); # trajectory
      plot_rel_trajectory("_ROR",$np_orb,$Orb[$m],$t_s,$t_s+$Orb[$m]->P(),"main"); # orbit
    }
  }
}
# for(;;) {print "."; usleep(-10);}
exit;
##################
#
# times from - to in days (epoch)
sub dump_pars {
 my ($file,%pars)=@_;
 my ($h,$k,$m,$lk);   
 
 open F,">".$file;
 foreach $key (keys %pars) {
   $h=$pars{$key};
   $k=$key;
   $lk=length($key);
   for($m=$lk;$m<40;++$m) {$k=$k." ";}
   if(($k !~ m/_O_/) && ($k !~ m/_T_/) && ($k !~ m/\[\]/)) {
     print F $k.$$h."\n";
   }
   if($k =~ m/_O_/) {
     @ar=@$h;
     for($n=0;$n<$n_orb;++$n) {       
       print F $k.$ar[$n]."\n";
     }
   }
   if($k =~ m/_T_/) {
     @ar=@$h;
     for($n=0;$n<$n_t_orb;++$n) {
       print F $k.$ar[$n]."\n";       
     }
   }
 }
 close F;
}
##################
#
# times from - to in days (epoch)
sub plot_abs_trajectory {
 my $ME="plot_abs_trajectory";
 my ($type,$vel,$npt,$obj,$tf,$tt,@CF)=@_;
 my ($m,$t,@pa,$vl,$pd,$day,$file);

 $day=24*3600;   
 $file=$obj->{$type};
 open F,">".$file;
 for($m=0;$m<$npt;++$m) {
   if($m==0) {$t=$tf;}
   else      {$t=$tf+$m*($tt-$tf)/($npt-1);} # time in seconds
   @pa=$obj->posn_cart_abs($t,$ME." ".$CF);
   @vl=$obj->vel_cart_abs($t,$ME." ".$CF);
   @pd=($pa[0]+$vl[0]*$o_vdt,$pa[1]+$vl[1]*$o_vdt,$pa[2]+$vl[2]*$o_vdt);    
   print F ($t/$o_t)."\t".
           ($pa[0]/$o_l)."\t".($pa[1]/$o_l)."\t".($pa[2]/$o_l)."\t".$pa[3]."\t".
           ($vl[0]/$o_v)."\t".($vl[1]/$o_v)."\t".($vl[2]/$o_v)."\n";
   if($vel eq "yes") {
   print F ($t/$o_t)."\t".
           ($pd[0]/$o_l)."\t".($pd[1]/$o_l)."\t".($pd[2]/$o_l)."\t".
           ($vl[0]/$o_v)."\t".($vl[1]/$o_v)."\t".($vl[2]/$o_v)."\n";       
   }
 }
 close F;
}
##################
#
# times from - to in days (epoch)
sub plot_rel_trajectory {
 my $ME="plot_rel_trajectory";
 my ($type,$npt,$obj,$tf,$tt,@CF)=@_;
 my ($m,$t,@pa,$vl,$pd,$day,$file);

 $day=24*3600;   
 $file=$obj->{$type};
 open F,">".$file;
 print F "#1:t, 2-5:r,M,E,T, 6-10:r,vr,vt,v,T";
 for($m=0;$m<$npt;++$m) {
   if($m==0) {$t=$tf;}
   else      {$t=$tf+$m*($tt-$tf)/($npt-1);} # time in seconds
   @pa=$obj->posn_pol_rel($t,$ME." ".$CF);
   @vl=$obj->vel_rrc_rel($t,$ME." ".$CF);
   print F ($t/$o_t)."\t".
           ($pa[0]/$o_l)."\t".$pa[1]."\t".$pa[2]."\t".$pa[3]."\t".
           ($vl[0]/$o_l)."\t".($vl[1]/$o_v)."\t".($vl[2]/$o_v)."\t".
           ($vl[3]/$o_v)."\t".$vl[4]."\n";
 }
 close F;
}

##################
#
#
##################
sub get_pars {
 my $parfile=shift;
 my $pars_h=shift;  # handle to hash of parameter names and var. handles
 my %pars=%$pars_h; # copy into local hash
 my (@orb,@tra,$Nname,@pf,$vl,$vh,@ar,$n);
 
 open F, $parfile;
 @pf=<F>;
 close F;
# @pf=include_args(@pf);
 @pf=include_files(@pf);
# @pf=multilines(@pf);
 @pf=expand_defines(@pf);

 @orb=get_array("orbits[]",@pf); # this value will be put into %pars, below
 $n_orb=scalar(@orb); # global var. n_orb
 @tra=get_array("transit[]",@pf); # this value will be put into %pars, below
 $n_t_orb=scalar(@tra); # global var. n_t_orb
 foreach $name (keys %pars) {
   $vh=$pars{$name}; # handle to scalar or array
   if($name =~ m/_O_/) {
     for($n=0;$n<scalar(@orb);++$n) {
       $Nname=$name;
       $Nname=~s/_O_/$orb[$n]/;
#       print $Nname;
       $ar[$n]=get_scalar($Nname,@pf);
     }
     @$vh=@ar;
   }
   if($name =~ m/\[\]/) {
     @ar=get_array($name,@pf);
     @$vh=@ar;
   }
   if(($name !~ m/\[\]/) && ($name !~ m/_O_/)) {
     $vl=get_scalar($name,@pf);
     if($vl ne "undef") {
       # $vh=$pars{$name};
       $$vh=$vl;
     }
   }
 }
}
##################
#
#
##################
sub expand_defines {
  my @pf=@_;
  my ($line,$nested,$n,%define);
  my $nd=1; # default for define nesting depth
    
  $nested="#depth of nested defines:";
  foreach $line (@pf) {
    if(substr($line,0,length($nested)) eq $nested) {
      $nd=substr($line,length($nested));
      $nd =~ s/^\s*(.*?)\s*$/$1/; # trim whitespace
    }
  }
  foreach $line (@pf) {
    if(substr($line,0,7) eq "#define") {
      $df=$line; $df=~s/^#define\s+([^\s\r\n]+)\s+.+[\r\n]+$/$1/; 
      $dt=$line; $dt=~s/^#define\s+([^\s\r\n]+)\s+([^\s\r\n]+)\s*[\r\n]+$/$2/g;
      $define{$df}=$dt;    
    }
    else { # only expand those strings that have been #define'd already
      for($n=0;$n<$nd;++$n) {
	foreach $key (keys %define) {$line =~ s/$key/$define{$key}/e;}
      }
    }      
  }
  @pf;
}
##################
#
#
##################
sub include_files {
  my @pf=@_;
  my ($line,$file,@in);
  my $nd=1; # default for include nesting depth
  
  foreach $line (@pf) {
    if(substr($line,0,8) eq "#include") {
      $file=substr($line,8);
      $file =~ s/^\s+([^\s\r\n]+)\s+[\r\n]+/$1/;
      open F,$file;
      @in=<F>;
      
      close F;
    }
  }     
  @pf;
}
##################
#
# if a parameter name appears multiple times, the las value will be used
##################
sub get_scalar {
  my ($name,@pf)=@_;
  my ($n,$val);
    
  $val="undef";
  foreach $line (@pf) {
    if(substr($line,0,length($name)+1) eq $name.":") {
      $val=substr($line,length($name)+1);
      $val =~ s/^\s*(.*?)\s*$/$1/; # trim whitespace
    }
  }
  if($val ne "undef") {
    if($val =~ m/^\".+\"[\r\n]*$/g) {$val =~ s/^\"(.+)\"[\r\n]*$/$1/;}
    else { $val=eval($val);}
  }
  $val;
}
##################
#
#
##################
sub get_array {
  my ($name,@pf)=@_;
  my ($n,$vals,@val);
  
  $val[0]="undef";
  for($n=0;$n<scalar(@pf);++$n) {
    if(substr($pf[$n],0,length($name)+1) eq $name.":") {
      $vals=substr($pf[$n],length($name)+1);
      $vals =~ s/^\s*(.*?)\s*$/$1/; # trim whitespace
      @val=split(/\s*,\s*/,$vals);
      for($m=0;$m<scalar(@val);++$m) {
        if($val[$m] =~ m/^\".+\"[\r\n]*$/g) 
	  {$val[$m] =~ s/^\"(.+)\"[\r\n]*$/$1/;}
        else {$val[$m]=eval($val[$m]);}
      }
#      print $val[0]."\t".$val[1]."\n";
      last;
    }
  }
  @val;
}

#######################
#
#
#######################
sub process_pars {
 my $ph=$_[0];
 my (%pars,$vh,@ln,$m,@ar);
 my @pr_sc=("start time","end time","scan, start time from",
	    "scan, start time to","scan, end time from","scan, end time to");
 my @pr_ar=("orbit _O_, time of periapsis");
    
 %pars=%$ph;
 for($n=0;$n<scalar(@pr_sc);++$n) {
   $vh=$pars{$pr_sc[$n]};
   if($$vh =~ m/\//) {
     @ln=split("/",$$vh);
     $$vh=aux::date_to_sec($ln[0],$ln[1],$ln[2]);   
   } 
 }
 for($n=0;$n<scalar(@pr_ar);++$n) {
   $vh=$pars{$pr_ar[$n]};
   @ar=@$vh;
   for($m=0;$m<scalar(@ar);++$m) {
     if($ar[$m] =~ m/\//) {
       @ln=split("/",$ar[$m]);
       $ar[$m]=aux::date_to_sec($ln[0],$ln[1],$ln[2]);      
     }   
   }
   @$vh=@ar; 
 }
}