attachment:lj2.pl of Brown2013

   1 #!/usr/local/bin/perl
   2 # lj2.pl   light pressure modified orbits graph
   3 #   lj2.pl > lj2.dat
   4 #   KHL  2013 Feb 16
   5 #
   6 # FIXME:  the search loop for eccentricity Really Sucks
   7 #         and blows up near the asymptote
   8 
   9 
  10 my    $amin =  6600.0           ; # km, bottom of LEO
  11 my    $amax =  384400.0         ; # km, lunar distance
  12 my   $steps =  500.0            ; # logarithmic step
  13 
  14 # grav parameters --------------------------------------
  15 my      $re =  6378.137         ; # earth radius km
  16 my      $J2 = -1.082626683e-03  ; # spherical harmonic
  17 my       $Y =  31556926.0       ; # seconds per year
  18 my      $mu =  398600.4418      ; # grav parameter  km3/s2
  19 my      $pi =  3.14159265358979 ; 
  20 
  21 # optical parameters -----------------------------------
  22 my      @th = (50,100)          ; # thickness μm
  23 my    $dens = 2.5e-3            ; # density kg/μ-m3
  24 my  $ipower = 1367.0            ; # light power W/m2
  25 my  $ifudge = 0.9               ; # for shadow and tilt
  26 my      $cl = 2.998e8           ; # speed of light m/s
  27 
  28 
  29 # derived ----------------------------------------------
  30 my    $adel = ($amax/$amin)**(1.0/($steps-1.0))  ;
  31 my    $npr0 = (-1.5*$J2*$Y/$pi)*sqrt($mu/$re**3) ;
  32 my     $ac0 = 1e-3*$ipower*$ifudge/($cl*$dens)   ; # acceleration
  33 my     $el0 = 3.0*$ac0*$Y / ( 4.0*$pi*sqrt($mu)) ;
  34 my     @el  = ( $el0/$th[0] , $el0/$th[1] )      ; # km/s^2
  35 my     @ee  = ( 0.0 , 0.0 )                      ;
  36 
  37 my     $npr ;
  38 
  39 # printf( "#npr0 = %12.6f\n", $npr0 );
  40 printf( "#s.major.axs     e(%4.0f)     e(%4.0f)\n", $th[0], $th[1] );
  41 
  42 # loop for the x range of semimajor axis
  43 for( my $a = $amin ; $a <= $amax+10 ; $a *= $adel ) {
  44    my     $are = $a/$re               ;
  45    my    $npr1 = $npr0 * $are**(-3.5) ;
  46 
  47    # loop for each acceleration
  48    foreach my $idx ( 0..1 ) {
  49       my  $el2 = $el[$idx] ;
  50       my   $e  = 0.0       ;
  51       my   $e0 = 0.5       ;
  52       my $elam = $el2 * sqrt( $a ) ;
  53    
  54       # loop to converge on $e eccentricity
  55       while( ( abs($e-$e0) > 1e-8 ) && ( abs($e) < 0.9 ) ) {
  56          my $nprm = ( $npr1 / (1.0-$e*$e)**2 ) - 1.0 ;
  57          my $enew  = $elam / $nprm ;
  58             $e0    = $e ;
  59             $e     = 0.99*$e + 0.01*$enew ;  ## KLUDGE
  60      }
  61      $ee[ $idx ] = abs( $e ) ;
  62   }   
  63   if( ( $ee[0]< 0.4 ) && ( $ee[1] < 0.4 ) ) {
  64      # 1000 km
  65      printf( "%12.4f%12.6f%12.6f", 0.001*$a, $ee[0], $ee[1] );
  66   }
  67   printf("\n");
  68 }