// gcc -o sunint sunint.c -lm ; ./sunint

#define  NM          10000
#define  N0          (3*NM)
#define  N1          (5*NM)
#define  N2          (7*NM)
#define  N3          (9*NM)
#define  NPOINTS    (12*NM)

#include  <stdio.h>
#include   <math.h>

main() {
   
   double pi2  = 8.0*atan( 1.0 );
   double npts = (double) NPOINTS ;
   double dpt  = pi2/npts ;

   double f1   = 0.0  ;  // first term integrated
   double f2s  = 0.0  ;  // second sin term integrated
   double f2c  = 0.0  ;  // second cos term integrated

   int    i           ;
   for( i = 0 ; i < NPOINTS ; i++ ) {
      double s = dpt * i ;
      double p = 0.0     ;
      if( ( i > N3 ) || ( i < N0 ) ) { p = 1.0 ; }
      else {
         double si = sin( s );
         double co = cos( s );
         if(      i < N1 ) { p = ( si-0.5)/sqrt(1.25-si) ; }
	 else if( i > N2 ) { p = (-si-0.5)/sqrt(1.25+si) ; }
      }
      if( p > 0.0 ) {
         s   *= 2.0 ;
         f1  += p   ;
         f2c += p*cos(s);
         f2s += p*sin(s);
      }
   }
   double dpd = 1.0/npts ;
   f1  *= dpd     ;
   f2c *= dpd/3.0 ;
   f2s *= dpd/3.0 ;

   printf( "%12.8f integrated thrust\n", f1     ); 
   printf( "%12.8f integ. sin thrust\n", f2s    ); 
   printf( "%12.8f integ. cos thrust\n", f2c    ); 
   printf( "%12.8f integ. tot thrust\n", f1-f2c ); 
}