#!/usr/local/bin/perl
# crosser.pl    v0.1  2013 Fri 25
#
# This finds the crossing density through the equatorial plane
#
# ./crosser.pl < {two line element 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  $rdel   = 200.0              ;  # km
my  $rmin   = 6400.0             ;  # km
my  $rmax   = 50000.0            ;  # km
my  $r288   = 12789.0            ;  # km
my  $nmin   = int( $rmin/$rdel)  ;  # minimum bin
my  $nmax   = 1 +int($rmax/$rdel);  # number of bins

my  $cntall = 0                  ;
my  $cnt288 = 0                  ;
my  $pi     = 4.0*atan2(1.0,1.0) ;
my  $mu     = 3.986004418e5      ; # km3/s2 (no oblate effects)
my  $v288   = sqrt($mu/$r288)    ;

my  @crosstab = (0.0)*$nmax ;
my  $fluxtab  = (0.0)*$nmax ;
my  $r        = (0.0)*2     ;

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 $prd = $el->period();                  # period sec
   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 ;

   # radius at equatorial crossings (two places)
   $r[0] = $sem * ( 1.0 - $ecc**2 ) / ( 1.0 + $ecc*cos($apg) ) ;
   $r[1] = $sem * ( 1.0 - $ecc**2 ) / ( 1.0 - $ecc*cos($apg) ) ;

   for( my $cross = 0 ; $cross <= 1 ; $cross += 1 ) {
      my $bin = int( $r[$cross] / $rdel ) ;
      $bin = $nmax if ( $bin > $nmax )    ;
      $crosstab[$bin] += 1.0/$rdel        ; # crossings in that bin
      # crossings per square meter per second 
      $fluxtab[$bin]  += 5e-7 / ( $prd * $pi * $rdel * $r[$cross] );
   }
}

printf( "# total objects=%8d\n", $cntall );
printf( "# radius  crossers   flux\n"    );

for( my $oc = $nmin ; $oc <= $nmax ; $oc += 1.0 ) {
   my $rr = $oc * $rdel / 1000.0 ;  # thousands of kilometers
   printf( "%10.3f%10.3f%10.2e\n", $rr, $crosstab[$oc], $fluxtab[$oc] );
}