home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_300 / 341_01 / sk.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-02-25  |  2.8 KB  |  93 lines
  1. /*
  2. HEADER:       sk.c        
  3. TITLE:        Solve Kepler's Equation         
  4. VERSION:      1.0    
  5. DESCRIPTION:  This routine solves Kepler's equation for eccentric
  6.               anomaly.  Iterative solution by Newton's method is used.
  7.               Required inputs are initial guess at eccentric anomaly,
  8.               value of mean anomaly, stopping tolerance, and maximum
  9.               number of iterations.  Input and output anomalies are
  10.               in degrees.
  11. KEYWORDS:     Astrodynamics, orbital mechanics, Kepler's equation
  12. SYSTEM:       MS-DOS, PC-DOS (Coded with Ver. 3.3, but should be version
  13.               independent)         
  14. FILENAME:     sk.c    
  15. WARNINGS:      This routine was coded for educational purposes, using
  16.               a standard iteration technique.  No attempt has been made
  17.               to insure optimal numerical stability or optimal convergence
  18.               rate.
  19. SEE-ALSO:     ---    
  20. AUTHORS:      Rodney Long
  21.               19003 Swan Drive
  22.               Gaithersburg, MD 20879     
  23. COMPILERS:    Microsoft C 5.1    
  24. REFERENCE:    Bate, Mueller, White: Fundamentals of Astrodynamics
  25.               Atkinson: An Introduction to Numerical Analysis
  26. */
  27.  
  28.  
  29. #include <stdio.h>
  30. #include <float.h>
  31. #include <math.h>
  32. #include "orbcons.h"
  33.  
  34. double sk();
  35.  
  36. double e[3] = {.1,.9,.999999999999}; //Input orbit eccen
  37.  
  38. main()
  39.  
  40. {
  41.   double eout;
  42.   double m, e0, t;
  43.   int i,mi;
  44.  
  45.  
  46.    e0   =  280;        // Input guess at eccen anom
  47.    m    =  286.879298; // Input mean anom
  48.    t    = .000001;      // Input stopping tolerance
  49.    mi   =  15;         // Input maximum iterations
  50.  
  51.    printf("Kepler's equation solved by Newton method \n");
  52.    printf(" --Beware of e very close to 1!-- \n");
  53.    printf("Initial eccentric anom           = %25.16lf \n",e0);
  54.    printf("Input mean anom                  = %25.16lf \n",m);
  55.    printf("Input maximum its                =      %d \n",mi);
  56.    printf("Input stop tolerance             = %25.16lf \n",t);
  57.  
  58.    for (i=0;i<1;i++) {
  59.      printf("Input eccentricity               = %25.16lf \n",e[i]);
  60.      eout = sk(m,e0,e[i],mi,t);
  61.      printf("Final eccentric anom             = %25.16lf \n",eout);
  62.      printf("Mean anom computed from solution = %25.16lf \n",
  63.             ((eout*DTR - e[i]*sin(eout*DTR))*RTD));
  64.    }
  65. }
  66.   
  67. double sk(double m,double e0,double e,int mi,double t)
  68. {
  69.   /* Solve Kepler's equation by the Newton method. */
  70.   double en,ep;
  71.   int i;
  72.  
  73.   e0 = e0 * DTR;
  74.   m  = m  * DTR;
  75.  
  76.   i  = 0;
  77.   en = e0;
  78.   do {
  79.     i++;
  80.     ep = en;
  81.  
  82.     // The next line of code is the Newton iteration:
  83.     //   x(n+1) = x(n) - (1/(dy/dx)) * y(n). 
  84.  
  85.     en = en - (en - e*sin(en) - m)/(1 - e*cos(en))   ;
  86.     en = fmod(en,TWOPI);
  87.     if ( fabs(en-ep) < t) break;
  88.   } while (i <= mi );
  89.   return(en*RTD);
  90. }
  91.     
  92.    
  93.