#!/usr/local/bin/perl
# lj2.pl   light pressure modified orbits graph
#   lj2.pl > lj2.dat
#   KHL  2013 Feb 16
#
# FIXME:  the search loop for eccentricity Really Sucks
#         and blows up near the asymptote


my    $amin =  6600.0           ; # km, bottom of LEO
my    $amax =  384400.0         ; # km, lunar distance
my   $steps =  500.0            ; # logarithmic step

# grav parameters --------------------------------------
my      $re =  6378.137         ; # earth radius km
my      $J2 = -1.082626683e-03  ; # spherical harmonic
my       $Y =  31556926.0       ; # seconds per year
my      $mu =  398600.4418      ; # grav parameter  km3/s2
my      $pi =  3.14159265358979 ; 

# optical parameters -----------------------------------
my      @th = (50,100)          ; # thickness μm
my    $dens = 2.5e-3            ; # density kg/μ-m3
my  $ipower = 1367.0            ; # light power W/m2
my  $ifudge = 0.9               ; # for shadow and tilt
my      $cl = 2.998e8           ; # speed of light m/s


# derived ----------------------------------------------
my    $adel = ($amax/$amin)**(1.0/($steps-1.0))  ;
my    $npr0 = (-1.5*$J2*$Y/$pi)*sqrt($mu/$re**3) ;
my     $ac0 = 1e-3*$ipower*$ifudge/($cl*$dens)   ; # acceleration
my     $el0 = 3.0*$ac0*$Y / ( 4.0*$pi*sqrt($mu)) ;
my     @el  = ( $el0/$th[0] , $el0/$th[1] )      ; # km/s^2
my     @ee  = ( 0.0 , 0.0 )                      ;

my     $npr ;

# printf( "#npr0 = %12.6f\n", $npr0 );
printf( "#s.major.axs     e(%4.0f)     e(%4.0f)\n", $th[0], $th[1] );

# loop for the x range of semimajor axis
for( my $a = $amin ; $a <= $amax+10 ; $a *= $adel ) {
   my     $are = $a/$re               ;
   my    $npr1 = $npr0 * $are**(-3.5) ;

   # loop for each acceleration
   foreach my $idx ( 0..1 ) {
      my  $el2 = $el[$idx] ;
      my   $e  = 0.0       ;
      my   $e0 = 0.5       ;
      my $elam = $el2 * sqrt( $a ) ;
   
      # loop to converge on $e eccentricity
      while( ( abs($e-$e0) > 1e-8 ) && ( abs($e) < 0.9 ) ) {
         my $nprm = ( $npr1 / (1.0-$e*$e)**2 ) - 1.0 ;
         my $enew  = $elam / $nprm ;
            $e0    = $e ;
            $e     = 0.99*$e + 0.01*$enew ;  ## KLUDGE
     }
     $ee[ $idx ] = abs( $e ) ;
  }   
  if( ( $ee[0]< 0.4 ) && ( $ee[1] < 0.4 ) ) {
     # 1000 km
     printf( "%12.4f%12.6f%12.6f", 0.001*$a, $ee[0], $ee[1] );
  }
  printf("\n");
}