attachment:sunint.c of LightOrbit

   1 // gcc -o sunint sunint.c -lm ; ./sunint
   2 
   3 #define  NM          10000
   4 #define  N0          (3*NM)
   5 #define  N1          (5*NM)
   6 #define  N2          (7*NM)
   7 #define  N3          (9*NM)
   8 #define  NPOINTS    (12*NM)
   9 
  10 #include  <stdio.h>
  11 #include   <math.h>
  12 
  13 main() {
  14    
  15    double pi2  = 8.0*atan( 1.0 );
  16    double npts = (double) NPOINTS ;
  17    double dpt  = pi2/npts ;
  18 
  19    double f1   = 0.0  ;  // first term integrated
  20    double f2s  = 0.0  ;  // second sin term integrated
  21    double f2c  = 0.0  ;  // second cos term integrated
  22 
  23    int    i           ;
  24    for( i = 0 ; i < NPOINTS ; i++ ) {
  25       double s = dpt * i ;
  26       double p = 0.0     ;
  27       if( ( i > N3 ) || ( i < N0 ) ) { p = 1.0 ; }
  28       else {
  29          double si = sin( s );
  30          double co = cos( s );
  31          if(      i < N1 ) { p = ( si-0.5)/sqrt(1.25-si) ; }
  32 	 else if( i > N2 ) { p = (-si-0.5)/sqrt(1.25+si) ; }
  33       }
  34       if( p > 0.0 ) {
  35          s   *= 2.0 ;
  36          f1  += p   ;
  37          f2c += p*cos(s);
  38          f2s += p*sin(s);
  39       }
  40    }
  41    double dpd = 1.0/npts ;
  42    f1  *= dpd     ;
  43    f2c *= dpd/3.0 ;
  44    f2s *= dpd/3.0 ;
  45 
  46    printf( "%12.8f integrated thrust\n", f1     ); 
  47    printf( "%12.8f integ. sin thrust\n", f2s    ); 
  48    printf( "%12.8f integ. cos thrust\n", f2c    ); 
  49    printf( "%12.8f integ. tot thrust\n", f1-f2c ); 
  50 }