attachment:crosser.pl of Brown2013

   1 #!/usr/local/bin/perl
   2 # crosser.pl    v0.1  2013 Fri 25
   3 #
   4 # This finds the crossing density through the equatorial plane
   5 #
   6 # ./crosser.pl < {two line element file}
   7 
   8 use Astro::SpaceTrack;
   9 use Astro::Coord::ECI;
  10 use Astro::Coord::ECI::TLE;
  11 use Astro::Coord::ECI::TLE::Set;
  12 use Astro::Coord::ECI::Utils qw(deg2rad rad2deg);
  13 
  14 my  $rdel   = 200.0              ;  # km
  15 my  $rmin   = 6400.0             ;  # km
  16 my  $rmax   = 50000.0            ;  # km
  17 my  $r288   = 12789.0            ;  # km
  18 my  $nmin   = int( $rmin/$rdel)  ;  # minimum bin
  19 my  $nmax   = 1 +int($rmax/$rdel);  # number of bins
  20 
  21 my  $cntall = 0                  ;
  22 my  $cnt288 = 0                  ;
  23 my  $pi     = 4.0*atan2(1.0,1.0) ;
  24 my  $mu     = 3.986004418e5      ; # km3/s2 (no oblate effects)
  25 my  $v288   = sqrt($mu/$r288)    ;
  26 
  27 my  @crosstab = (0.0)*$nmax ;
  28 my  $fluxtab  = (0.0)*$nmax ;
  29 my  $r        = (0.0)*2     ;
  30 
  31 foreach $el ( Astro::Coord::ECI::TLE->parse( <> ) ) {
  32   
  33    my $nam = $el->get('name');               # name
  34    my $bod = $el->body_type();               # body type
  35    my $apo = $el->apogee();                  # apogee km
  36    my $per = $el->perigee();                 # perigee km
  37    my $sem = $el->semimajor();               # semimajor axis km
  38    my $prd = $el->period();                  # period sec
  39    my $ecc = $el->get('eccentricity');       # eccentricity
  40    my $inc = $el->get('inclination');        # inclination radians
  41    my $apg = $el->get('argumentofperigee');  # arg of perigee radians
  42    my $v01 = sqrt( $mu*( 1.0 - $ecc**2 )/$sem ); # orbit velocity
  43 
  44    $cntall += 1 ;
  45 
  46    # radius at equatorial crossings (two places)
  47    $r[0] = $sem * ( 1.0 - $ecc**2 ) / ( 1.0 + $ecc*cos($apg) ) ;
  48    $r[1] = $sem * ( 1.0 - $ecc**2 ) / ( 1.0 - $ecc*cos($apg) ) ;
  49 
  50    for( my $cross = 0 ; $cross <= 1 ; $cross += 1 ) {
  51       my $bin = int( $r[$cross] / $rdel ) ;
  52       $bin = $nmax if ( $bin > $nmax )    ;
  53       $crosstab[$bin] += 1.0/$rdel        ; # crossings in that bin
  54       # crossings per square meter per second 
  55       $fluxtab[$bin]  += 5e-7 / ( $prd * $pi * $rdel * $r[$cross] );
  56    }
  57 }
  58 
  59 printf( "# total objects=%8d\n", $cntall );
  60 printf( "# radius  crossers   flux\n"    );
  61 
  62 for( my $oc = $nmin ; $oc <= $nmax ; $oc += 1.0 ) {
  63    my $rr = $oc * $rdel / 1000.0 ;  # thousands of kilometers
  64    printf( "%10.3f%10.3f%10.2e\n", $rr, $crosstab[$oc], $fluxtab[$oc] );
  65 }