#!/usr/local/bin/perl
# insertrb.pl    v0.2  2013 Sat 26
#
# This computes the two pairs of delta V burns to transfer a NORAD
# object to m288 orbit for rocket bodies
#
# ./insertrb.pl < {three line element file} > {output file}

use Astro::SpaceTrack;
use Astro::Coord::ECI;
use Astro::Coord::ECI::TLE;
use Astro::Coord::ECI::TLE::Set;
use Astro::Coord::ECI::Utils qw(deg2rad rad2deg);

my $r288   = 12789.0            ;  # km

my $pi     = 4.0*atan2(1.0,1.0) ;
my $mu     = 3.986004418e5      ; # km3/s2 (no oblate effects)
my $v288   = sqrt($mu/$r288)    ;
my $leo    = 2000+6370          ; # low earth orbit semimajor radius

my @cislist = ("BLOCK","BREEZE","DNEPR","FREGAT","STRELA","SL-","SS-","SSN-");

my $cntall  = 0                 ;
my $cntdeb  = 0                 ;
my $cntrb   = 0                 ;
my $cntleo  = 0                 ;
my $cnthigh = 0                 ;
my $cntcis  = 0                 ;
my $cntsl6  = 0                 ;

# array for delta V histograms
my $rbhist   = (0) x 30         ; # all rocket bodies
my @leohist  = (0) x 30         ; # rocket bodies in LEO
my @highhist = (0) x 30         ; # rocket bodies above LEO
my $cishist  = (0) x 30         ; # CIS/Soviet rockets
my @sl6hist  = (0) x 30         ; # SL-6 above LEO (Molniya)

my $dvstep   = 0.5              ; # array bin size 500m/s
my $dvmax    = 11.0             ; # max delta V binned, 10 km/s
my $dvinmax  = 22               ;

# printf("\n%7.0f=r288%7.3f=v288\n", $r288, $v288 );

open(  F , ">rb.deltav");

foreach $el ( Astro::Coord::ECI::TLE->parse( <> ) ) {
  
   my $nam = $el->get('name');               # name
   my $bod = $el->body_type();               # body type
   my $apo = $el->apogee();                  # apogee km
   my $per = $el->perigee();                 # perigee km
   my $sem = $el->semimajor();               # semimajor axis km
   my $ecc = $el->get('eccentricity');       # eccentricity
   my $inc = $el->get('inclination');        # inclination radians
   my $apg = $el->get('argumentofperigee');  # arg of perigee radians
   my $v01 = sqrt( $mu*(1.0-$ecc**2)/$sem ); # orbit velocity

   $cntall += 1 ;
   $cntdeb += 1 if ( $nam =~ / DEB/ ) ;

   if( $nam =~ /R\/B/ ) {

      # printf( "\n%7.3f=i%7.3f=ap%7.3f=v0", $inc, $apg, $v01 );
      # printf( "\n%7.0f=sem%7.3f=ecc", $sem, $ecc );
   
      # radius at equatorial crossings (two places)
      my $r11 = $sem * ( 1.0 - $ecc**2 ) / ( 1.0 + $ecc*cos($apg) ) ;
      my $r12 = $sem * ( 1.0 - $ecc**2 ) / ( 1.0 - $ecc*cos($apg) ) ;
   
      # printf( "%7.0f=ra%7.0f=rp%7.0f=r11%7.0f=r12\n",$apo,$per,$r11,$r12 );
   
      #  tangential orbit velocity
      my $vt11 = $v01 * ( 1.0 + $ecc*cos($apg) ) ;
      my $vt12 = $v01 * ( 1.0 - $ecc*cos($apg) ) ;
      
      #  tangential orbit velocity in plane
      my $vh11 = $vt11 * cos($inc);
      my $vh12 = $vt12 * cos($inc);
   
      #  tangential orbit velocity out of plane
      my $vv11 = $vt11 * sin($inc);
      my $vv12 = $vt12 * sin($inc);
   
      #  radial orbit velocity
      my $vr11 =  $ecc * $v01 * sin($apg) ;
      my $vr12 = -$ecc * $v01 * sin($apg) ;
   
      # printf( "%7.3f=vt1%7.3f=vt2%7.3f=vr1%7.3f=vr2\n",$vt11,$vt12,$vr11,$vr12);
   
      #  transfer orbit semimajor axis
      my $sem21 = 0.5 * ( $r11 + $r288 ) ;
      my $sem22 = 0.5 * ( $r12 + $r288 ) ;
   
      #  transfer orbit eccentricity
      my $ecc21 = ( $r11/$sem21 ) - 1.0  ;
      my $ecc22 = ( $r12/$sem22 ) - 1.0  ;
   
      # transfer orbit velocity
      my $v021 = sqrt( $mu * ( 1.0 - $ecc21**2 ) / $sem21 );
      my $v022 = sqrt( $mu * ( 1.0 - $ecc22**2 ) / $sem22 );
   
      # transfer orbit r1 velocity
      my $vh21 = ( 1.0 - $ecc21 ) * $v021 ;
      my $vh22 = ( 1.0 - $ecc22 ) * $v022 ;
   
      # transfer orbit r288 velocity
      my $vh31 = ( 1.0 + $ecc21 ) * $v021 ;
      my $vh32 = ( 1.0 + $ecc22 ) * $v022 ;
   
      # deltaV to transfer orbit
      my $dv11 = sqrt( $vr11**2 + $vv11**2 + ($vh21-$v288)**2 );
      my $dv12 = sqrt( $vr12**2 + $vv12**2 + ($vh22-$v288)**2 );
   
      # deltaV transfer orbit to m288
      my $dv21 = abs( $vh21-$v288 ) ;
      my $dv22 = abs( $vh22-$v288 ) ;
   
      # sort the values, lowest first
      if ( ( $dv12 < $dv11 ) && ( $dv22 < $dv21 ) ) {
         my $tmp = $dv11 ; $dv11 = $dv12 ; $dv12 = $tmp ;
            $tmp = $dv21 ; $dv21 = $dv22 ; $dv22 = $tmp ;
      }

      my $x = " ";    # extra character
      $x = "+" if($dv21 > $dv22) ;
   
      my $dvt = $dv11 + $dv21 ;
      my $dvindex = int ($dvt/$dvstep) ;

          $rbhist[$dvindex]+=1;   $cntrb+=1;  $x="O"; 

      foreach my $cis ( @cislist ) {
        if( $nam =~ $cis )  
        { $cishist[$dvindex]+=1;  $cntcis+=1; $x="S"; }
      }

      if( $nam =~ /SL-6/ )
        { $sl6hist[$dvindex]+=1;  $cntsl6+=1; $x="M"; }

      if( $sem < $leo    )
        { $leohist[$dvindex]+=1;  $cntleo+=1;  }
      else                   
        { $highhist[$dvindex]+=1; $cnthigh+=1; }

      printf F ("%7.3f%7.3f%8.3f%7.3f %s %s\n",
              $dv11, $dv21,  $dv12, $dv22,  $x, $nam);
   }
}
close F ;

my $cntncis = $cntrb - $cntcis ;

open(  F , ">rb.stat");
printf F ("\ntotal=%5d",               $cntall  );
printf F (  "   debris=%5d  ",            $cntdeb  );
printf F (  "   rocket bodies=%5d\n",     $cntrb   );
printf F (  "LEO R/B=%4d   ",          $cntleo  );
printf F (  "   High R/B=%4d\n",          $cnthigh );
printf F (  "CIS/Soviet=%4d",          $cntcis  );
printf F (  "   Soyuz/Molniya=%4d",       $cntsl6  );
printf F (  "   NonCIS=%4d\n",            $cntncis );
close  F ;

open(  F, ">rb.dat" );
printf F (  "delV   ALL  HIGH   CIS  SL-6\n"    );

for( my $i = 0 ; $i <= $dvinmax ; $i++ ) {
   printf F ("%5.1f%6d%6d%6d%6d\n", $i*$dvstep,
      $rbhist[$i], $highhist[$i], $cishist[$i], $sl6hist[$i] );
} 
close F ;