// t26.c
// v0.1 KHL 2009-04-21
// gcc -o t26 t26.c -lm -lgd ; ./t26
// static drawing, though slow to compute
// earth and toroid around orbit

#define TITLE    "M288, 2500 orbits"
#define GIFOUT   "t26.gif"
#define NORB    -24
#define DELTANG 0.2
#define ORBWIDTH  9
#define ANG0      0
#define ANG1     90
#define ANG2    142
#define NELIP    32
#define RADIUS    0.44367

#include "tor00.hc"

int main () {
   double  t_ang ;
   int     e_cnt ;
   int     orb ;

   gdPoint ellipse0[NELIP+2];
   gdPoint ellipse1[NELIP+2];

   displaystart();

   torstart();

   gdImageGifAnimBegin( im1, gifout, 1, -1 ) ; // no repeat
   // gdImageGifAnimBegin( im1, gifout, 1,  4 ) ; // repeat 4 times
   // gdImageGifAnimBegin( im1, gifout, 1,  0 ) ; // continuous repeat

   framestart(    90, TITLE );

   // draw back portion of the orbit toroid
  
   double acolor =  0.0 ;

   for( t_ang=ANG0 ; t_ang < ANG2 ; t_ang += DELTANG ) {

      double ecolor ;
      if( acolor < 0.5 ) { ecolor = gray ; }
      else {
         ecolor = lgray ;
         if( acolor > 14.5 ) { acolor -= 15.0 ; }
      }
      acolor += DELTANG ;

      double angle_rad=DEG2RAD( t_ang );

      // ellipse disk - since it is at an angle, we will have to draw it
      // as a filled polygon
    
      // ellipse disk on left
      for( e_cnt=0 ; e_cnt < NELIP ; e_cnt++ ) {
         double toa = angle_rad + (2.0*PI*e_cnt)/NELIP  ;
         double tos = RADIUS          ;
         double tor = tos * cos( toa );
         double toy = tos * sin( toa );
         double tox = -(major+tor)*sin( angle_rad );
	 double toz = -(major+tor)*cos( angle_rad );
	 ellipse0[e_cnt] = sr0( tox, toy, toz );
      }
      gdImagePolygon( im, ellipse0, NELIP, ecolor );

      // ellipse disk on right
      for( e_cnt=0 ; e_cnt < NELIP ; e_cnt++ ) {
         double toa = -angle_rad + (2.0*PI*e_cnt)/NELIP  ;
         double tos = RADIUS          ;
         double tor = tos * cos( toa );
         double toy = tos * sin( toa );
         double tox = -(major+tor)*sin( -angle_rad );
         double toz = -(major+tor)*cos( -angle_rad );
         ellipse1[e_cnt] = sr0( tox, toy, toz );
      }
      gdImagePolygon( im, ellipse1, NELIP, ecolor );
   }

   gdImageFilledPolygon( im, ellipse0, NELIP, black );
   gdImagePolygon( im, ellipse0, NELIP, lgray  );

   gdImageFilledPolygon( im, ellipse1, NELIP, black );
   gdImagePolygon( im, ellipse1, NELIP, lgray  );

   int    ptx ;
   int    pty ;
   gdPoint ptxy ;
   double ptd = 0.1*step ;  
   double ang2_rad =DEG2RAD( ANG2 );
   int    ptmy = (int) sqrt(800.0);

   // orbit points on left and right
   for( pty = -ptmy ; pty <= ptmy ; pty++ ) {
      int ptmx = (int) sqrt( (double)( 800 - pty*pty ) );
      for( ptx = -ptmx ; ptx <= ptmx ; ptx++ ) {
         double ptx0 = ptd * ptx ;
         double pty0 = ptd * pty ;
         double ptr1 = ptx0 * cos( ang2_rad ) - pty0 * sin( ang2_rad ) + major ;
         double pty1 = pty0 * cos( ang2_rad ) + ptx0 * sin( ang2_rad ) ;
         double ptz1 = -ptr1 * cos( ang2_rad ) ;
         double ptx1 = -ptr1 * sin( ang2_rad ) ;

         ptxy = sr0( ptx1, pty1, ptz1 );
         gdImageSetPixel( im, ptxy.x, ptxy.y, white );
    
         ptxy = sr0( -ptx1, pty1, ptz1 );
         gdImageSetPixel( im, ptxy.x, ptxy.y, white );
      }
   }
    
   drawearth();
   sprintf( bottom, "%5.0fkm spacing between orbits", 637.80*step/re );

   frameend();
   displayend();
   return 0;
}