attachment:locus01.pl of NearCircularOrbits

   1 #!/usr/local/bin/perl
   2 # draw orbit with semimajor axis 1,  radius eccentricity e
   3 # then rotate by mean anomaly
   4 
   5 use Math::Trig;
   6 
   7 my   $pi           = 4.0*atan(1.0)      ;
   8 my   $radian       = $pi / 180.0        ;
   9 my   @ea = (1e-3, 3e-3, 1e-2,  3e-2, 1e-1, 3e-1 );
  10 ## my   @ea = ( 1e-2 );
  11 my   $da = 0.5 ;
  12 ## my   $da = 15.0 ;
  13 
  14 ## printf "  theta  sin   cos  r      xr    yr    ce    se" ;
  15 ## printf "    EE     M      x         y\n";
  16 
  17 ## foreach my $e ( @ea ) { printf "%11.3e   ", $e } ; printf  "\n" ;
  18 
  19 # compute from apogee to apogee
  20 for( my $an  = -180.0 ; $an < 180.00001 ; $an += $da ) {
  21    my $theta = $radian * $an ;
  22    $s = sin( $theta );
  23    $c = cos( $theta );
  24    foreach my $e ( @ea ) {
  25       my $rr = (1.0-$e*$e) / (1.0 + $e*$c );
  26       my $xr = -$rr * $s ;
  27       my $yr =  $rr * $c ;
  28 
  29       # it is difficult to compute E with the standard acos($ce) 
  30       # formula, because quadrant information is lost.
  31       #
  32       # However, if we compute both the cosine of E ($ce)  and the
  33       # sin ($se), we can use the atan2 function to get $E, since
  34       # $ce and $se together specify the quadrant.
  35 
  36       my $ce = ( $e + $c ) / ( 1.0 + $e*$c );
  37       my $se = sqrt( 1.0 - $e*$e ) * $s / ( 1.0 + $e*$c ) ;
  38 
  39       my $E  = atan2( $se, $ce );
  40 
  41       my $M = $E - $e * sin( $E );
  42 
  43       my $cm = cos( $M );
  44       my $sm = sin( $M );
  45 
  46       # rotate $xr, $yr by angle $M
  47 
  48       my $x = ( $cm * $xr + $sm * $yr ) / $e ;
  49       my $y = ( $cm * $yr - $sm * $xr - 1.0 ) / $e ;
  50 
  51       ## my $x = ( $cm * $xr + $sm * $yr ) ;
  52       ## my $y = ( $cm * $yr - $sm * $xr ) ;
  53 
  54 ##  printf "%7.3f%6.2f%6.2f%6.3f%6.2f%6.2f%6.2f%6.2f%7.3f%7.3f%10.2e%10.2e",
  55 ##          $theta, $s,  $c, $rr, $xr, $yr, $ce, $se,  $E,  $M,   $x,   $y ;
  56       printf "%7.3f%7.3f",  $x, $y ;
  57    }
  58    printf "\n";
  59 }