attachment:t26.c of ToroidalOrbits

   1 // t26.c
   2 // v0.1 KHL 2009-04-21
   3 // gcc -o t26 t26.c -lm -lgd ; ./t26
   4 // static drawing, though slow to compute
   5 // earth and toroid around orbit
   6 
   7 #define TITLE    "M288, 2500 orbits"
   8 #define GIFOUT   "t26.gif"
   9 #define NORB    -24
  10 #define DELTANG 0.2
  11 #define ORBWIDTH  9
  12 #define ANG0      0
  13 #define ANG1     90
  14 #define ANG2    142
  15 #define NELIP    32
  16 #define RADIUS    0.44367
  17 
  18 #include "tor00.hc"
  19 
  20 int main () {
  21    double  t_ang ;
  22    int     e_cnt ;
  23    int     orb ;
  24 
  25    gdPoint ellipse0[NELIP+2];
  26    gdPoint ellipse1[NELIP+2];
  27 
  28    displaystart();
  29 
  30    torstart();
  31 
  32    gdImageGifAnimBegin( im1, gifout, 1, -1 ) ; // no repeat
  33    // gdImageGifAnimBegin( im1, gifout, 1,  4 ) ; // repeat 4 times
  34    // gdImageGifAnimBegin( im1, gifout, 1,  0 ) ; // continuous repeat
  35 
  36    framestart(    90, TITLE );
  37 
  38    // draw back portion of the orbit toroid
  39   
  40    double acolor =  0.0 ;
  41 
  42    for( t_ang=ANG0 ; t_ang < ANG2 ; t_ang += DELTANG ) {
  43 
  44       double ecolor ;
  45       if( acolor < 0.5 ) { ecolor = gray ; }
  46       else {
  47          ecolor = lgray ;
  48          if( acolor > 14.5 ) { acolor -= 15.0 ; }
  49       }
  50       acolor += DELTANG ;
  51 
  52       double angle_rad=DEG2RAD( t_ang );
  53 
  54       // ellipse disk - since it is at an angle, we will have to draw it
  55       // as a filled polygon
  56     
  57       // ellipse disk on left
  58       for( e_cnt=0 ; e_cnt < NELIP ; e_cnt++ ) {
  59          double toa = angle_rad + (2.0*PI*e_cnt)/NELIP  ;
  60          double tos = RADIUS          ;
  61          double tor = tos * cos( toa );
  62          double toy = tos * sin( toa );
  63          double tox = -(major+tor)*sin( angle_rad );
  64 	 double toz = -(major+tor)*cos( angle_rad );
  65 	 ellipse0[e_cnt] = sr0( tox, toy, toz );
  66       }
  67       gdImagePolygon( im, ellipse0, NELIP, ecolor );
  68 
  69       // ellipse disk on right
  70       for( e_cnt=0 ; e_cnt < NELIP ; e_cnt++ ) {
  71          double toa = -angle_rad + (2.0*PI*e_cnt)/NELIP  ;
  72          double tos = RADIUS          ;
  73          double tor = tos * cos( toa );
  74          double toy = tos * sin( toa );
  75          double tox = -(major+tor)*sin( -angle_rad );
  76          double toz = -(major+tor)*cos( -angle_rad );
  77          ellipse1[e_cnt] = sr0( tox, toy, toz );
  78       }
  79       gdImagePolygon( im, ellipse1, NELIP, ecolor );
  80    }
  81 
  82    gdImageFilledPolygon( im, ellipse0, NELIP, black );
  83    gdImagePolygon( im, ellipse0, NELIP, lgray  );
  84 
  85    gdImageFilledPolygon( im, ellipse1, NELIP, black );
  86    gdImagePolygon( im, ellipse1, NELIP, lgray  );
  87 
  88    int    ptx ;
  89    int    pty ;
  90    gdPoint ptxy ;
  91    double ptd = 0.1*step ;  
  92    double ang2_rad =DEG2RAD( ANG2 );
  93    int    ptmy = (int) sqrt(800.0);
  94 
  95    // orbit points on left and right
  96    for( pty = -ptmy ; pty <= ptmy ; pty++ ) {
  97       int ptmx = (int) sqrt( (double)( 800 - pty*pty ) );
  98       for( ptx = -ptmx ; ptx <= ptmx ; ptx++ ) {
  99          double ptx0 = ptd * ptx ;
 100          double pty0 = ptd * pty ;
 101          double ptr1 = ptx0 * cos( ang2_rad ) - pty0 * sin( ang2_rad ) + major ;
 102          double pty1 = pty0 * cos( ang2_rad ) + ptx0 * sin( ang2_rad ) ;
 103          double ptz1 = -ptr1 * cos( ang2_rad ) ;
 104          double ptx1 = -ptr1 * sin( ang2_rad ) ;
 105 
 106          ptxy = sr0( ptx1, pty1, ptz1 );
 107          gdImageSetPixel( im, ptxy.x, ptxy.y, white );
 108     
 109          ptxy = sr0( -ptx1, pty1, ptz1 );
 110          gdImageSetPixel( im, ptxy.x, ptxy.y, white );
 111       }
 112    }
 113     
 114    drawearth();
 115    sprintf( bottom, "%5.0fkm spacing between orbits", 637.80*step/re );
 116 
 117    frameend();
 118    displayend();
 119    return 0;
 120 }