// tor00.hc // v0.1 KHL 2009 Apr 20 // combined includes and C code, because I am lazy // This needs to be a library - or something #define FNT "DejaVuMonoSans" #define FSZ 24 // this scales the font and the whole drawing #define ASTEP0 400 // starting and ending hold times #define ASTEP1 80 // image rate (1.25 images/second ) #define RE 1.0 // earth width (arbitrary units ) #define MAJOR 2.0 // orbit width (arbitrary units ) #define MINOR 0.4 // #define STEP 0.156862 // distance between orbits #define XVIEW 5.0 // display width ( arbitrary units ) #define ATILT 23.432 // Axial Tilt - full summer #define XSIZE 40*FSZ #define YSIZE 26*FSZ #define XCENTER 20*FSZ #define YCENTER 14*FSZ #define EPTS 400 // points in circle to draw earth and lat/long #define OPTS 800 // points in circle to draw orbit #include #include #include #include #define PI 3.14159265359 #define DEG2RAD(x) ((x)*PI/180.) // degrees -> radians FILE *gifout; ; // output file gdImagePtr im0, im1, im; ; // working images int black, white, gray, red, green, blue ; // colors int dgray, cyan, yellow, magenta, trans ; // colors int lgray ; // colors int col00, col01, col02, col03, col04 ; // colors int col05, col06, col07, col08, col09 ; // colors int col10, col11, col12, col13, col14 ; // colors int col15, col16, col17, col18, col19 ; // colors double scale = ((double)XSIZE)/((double)XVIEW) ; // pixels per unit double atilt = DEG2RAD( ATILT ) ; // axial tilt in rads double ax, az ; // axial tilt rotations int astep = ASTEP1 ; // frame display time x10ms double re=RE ; // earth radius double major=MAJOR ; double minor=MINOR ; double xcenter=XCENTER ; double ycenter=YCENTER ; double step=STEP ; double long0=0.0 ; // zero longitude double zz=0.0 ; //-------------------------------------------------------------------- // rotation, transformation, scaling -------- // gdPoint s0( double inx, double iny ) { gdPoint retvar ; retvar.x = (int) ( xcenter + scale * inx ) ; retvar.y = (int) ( ycenter - scale * iny ) ; return retvar; } // without clipping gdPoint sr0( double inx, double iny, double inz ) { gdPoint retvar ; // rotate, IN TWO ONE PLANES - z is out of the plane, toward the observer // y is vertical for this program (inconsistent ) // first rotation on x axis (front to back ) double fx = inx ; double fy = iny * cos( ax ) - inz * sin( ax ) ; double fz = iny * sin( ax ) + inz * cos( ax ) ; // second rotation on z axis double rx = fx * cos( az ) - fy * sin( az ) ; double ry = fx * sin( az ) + fy * cos( az ) ; double rz = fz ; zz = rz ; // copy to testable global retvar.x = (int) ( xcenter + scale * rx ) ; retvar.y = (int) ( ycenter - scale * ry ) ; return retvar; } gdPoint sr( double inx, double iny, double inz ) { gdPoint retvar ; // rotate, IN TWO ONE PLANES - z is out of the plane, toward the observer // y is vertical for this program (inconsistent ) // first rotation on x axis (front to back ) double fx = inx ; double fy = iny * cos( ax ) - inz * sin( ax ) ; double fz = iny * sin( ax ) + inz * cos( ax ) ; // second rotation on z axis double rx = fx * cos( az ) - fy * sin( az ) ; double ry = fx * sin( az ) + fy * cos( az ) ; double rz = fz ; double rr = sqrt( rx*rx + ry*ry) ; zz = rz ; // copy to testable global if( ( rz < 0 ) && ( rr < re ) ) { retvar.x = -1 ; // behind earth retvar.y = -1 ; } else { retvar.x = (int) ( xcenter + scale * rx ) ; retvar.y = (int) ( ycenter - scale * ry ) ; } return retvar; } // ----------------------------------------------------------------- // draw a picture of the earth at a particular "date", 0 to 360 void drawearth() { gdPoint epoint[EPTS]; gdPoint epoint1; double dangle = (2.0 * PI)/( EPTS-1 ); double rr, rx, ry, rz ; double off ; // latitude/longitude starting drawing offset double offt ; // test offset double ox, oz ; int pcount ; int cntr ; double l ; double minz ; gdImageSetThickness( im, 1 ); // draw background color and edge for( cntr=0 ; cntr < EPTS ; cntr++ ) { rx = re*cos( dangle * cntr ) ; ry = re*sin( dangle * cntr ) ; epoint[cntr] = s0( rx, ry ) ; } gdImageFilledPolygon( im, epoint, EPTS, cyan ); gdImagePolygon( im, epoint, EPTS, dgray ); // draw latitude lines for( l = -75.0 ; l < 89.9 ; l += 15.0 ) { rr = re*cos( DEG2RAD(l) ) ; ry = re*sin( DEG2RAD(l) ) ; pcount = 0 ; // find a starting computation offset, behind the earth minz = 1.0 ; for( offt=0.0 ; offt < 2.0*PI ; offt += 0.25*PI ) { rx = rr*sin( offt ); rz = -rr*cos( offt ); epoint1 = sr( rx, ry, rz ); if( minz > zz ) { minz = zz ; off = offt ; } } // draw curve if transformed z > 0.0 for( cntr=0 ; cntr < EPTS ; cntr++ ) { rx = rr*sin( dangle*cntr + off ) ; rz = -rr*cos( dangle*cntr + off ) ; epoint1 = sr( rx, ry, rz ); if( epoint1.x >= 0 ) { epoint[ pcount++ ] = epoint1 ; } } if ( l*l < 1.0 ) { gdImageSetThickness( im, 3 ); // fat equator } else { gdImageSetThickness( im, 1 ); } gdImageOpenPolygon( im, epoint, pcount, dgray ); } // draw longitude lines; start at farthest back point for( l = -90.000 ; l < 89.999 ; l += 15.0 ) { pcount = 0 ; ox = re * sin( DEG2RAD(l-long0) ) ; oz = -re * cos( DEG2RAD(l-long0) ) ; // find a starting computation offset, behind the earth minz = 1.0 ; for( offt=0.0 ; offt < 2.0*PI ; offt += 0.25*PI ) { rx = ox*cos( offt ) ; ry = re*sin( offt ) ; rz = oz*cos( offt ) ; epoint1 = sr( rx, ry, rz ); if( minz > zz ) { minz = zz ; off = offt ; } } // draw curve if transformed z > 0.0 for( cntr=0 ; cntr < EPTS ; cntr++ ) { rx = ox*cos( dangle*cntr + off ) ; ry = re*sin( dangle*cntr + off ) ; rz = oz*cos( dangle*cntr + off ) ; epoint1 = sr( rx, ry, rz ); if( epoint1.x >= 0 ) { epoint[ pcount++ ] = epoint1 ; } } gdImageOpenPolygon( im, epoint, pcount, dgray ); } } // ----------------------------------------------------------------- // draw main orbit void drawmain() { gdPoint opoint[OPTS]; gdPoint opoint1; double dangle = (2.0 * PI)/( OPTS-1 ); double rx, ry, rz ; int pcount ; int cntr ; gdImageSetThickness( im, 3 ); pcount = 0 ; for( cntr=0 ; cntr < OPTS ; cntr++ ) { rx = major*sin( dangle * cntr ) ; rz = -major*cos( dangle * cntr ) ; opoint1 = sr( rx, 0.0, rz ); if( opoint1.x >= 0 ) { opoint[ pcount++ ] = opoint1 ; } } gdImageOpenPolygon( im, opoint, pcount, white ); } // ----------------------------------------------------------------- // draw toroid orbit void drawtor( int color, int dstep, int an_deg) { gdPoint opoint[OPTS]; gdPoint opoint1; double dangle = (2.0 * PI)/( OPTS-1 ); double rx, ry, rz ; int pcount ; int cntr ; double an_rad = DEG2RAD( an_deg ) ; double da ; // ascending angle double dr ; // radius modification double dy ; // vertical modification gdImageSetThickness( im, 3 ); pcount = 0 ; for( cntr=0 ; cntr < OPTS ; cntr++ ) { // toroid modifications da = dangle * cntr + an_rad ; dr = dstep * step * cos( da ); // cheesy toroid rather than real orbit dy = dstep * step * sin( da ); // no apogee skew // main orbit rx = (major+dr)*sin( dangle * cntr ) ; rz = -(major+dr)*cos( dangle * cntr ) ; opoint1 = sr( rx, dy, rz ); if( opoint1.x >= 0 ) { opoint[ pcount++ ] = opoint1 ; } } gdImageOpenPolygon( im, opoint, pcount, color ); } //------------------------------------------------------------------ // draw one body in main orbit void drawbody1(double angle_deg, int scolor, int size) { gdPoint opoint1; double angle_rad = DEG2RAD( angle_deg) ; double rx, ry, rz ; rx = major*sin( angle_rad ) ; rz = -major*cos( angle_rad ) ; opoint1 = sr( rx, 0.0, rz ); if( opoint1.x >= 0 ) { gdImageFilledEllipse( im, opoint1.x, opoint1.y, size, size, white ); } } // --------------------------------------- double fr_top_siz = 1.3*FSZ ; double fr_bot_siz = FSZ ; char bottom[80]; void framestart( int date, char* comment ) { sprintf( bottom, "viewed from the sun" ) ; // default bottom comment double rangle = DEG2RAD( date ); im = gdImageCreate( XSIZE, YSIZE ); gdImagePaletteCopy( im, im1 ); ax = atilt * sin( rangle ); az = -atilt * cos( rangle ); gdImageStringFT( im, NULL, white, FNT, fr_top_siz, 0.0, FSZ, 2*FSZ, comment ); } // --------------------------------------- void frameend() { gdImageStringFT( im, NULL, white, FNT, fr_bot_siz, 0.0, FSZ, YSIZE-FSZ, bottom ); // may be changed in main gdImageGifAnimAdd ( im, gifout, 0, 0, 2, astep, gdDisposalRestoreBackground, im1 ); im0 = im1 ; im1 = im ; gdImageDestroy (im0); } // ----------------------------------------------------------------- // draw all, one frame void drawall( int date, char* comment ) { framestart(date, comment); drawearth(); drawmain(); frameend(); } // ----------------------------------------------------------------- void displaystart() { im1 = gdImageCreate (XSIZE, YSIZE); white = gdImageColorAllocate (im1, 255, 255, 255); black = gdImageColorAllocate (im1, 0, 0, 0); gray = gdImageColorAllocate (im1, 128, 128, 128); dgray = gdImageColorAllocate (im1, 64, 64, 64); lgray = gdImageColorAllocate (im1, 192, 192, 192); cyan = gdImageColorAllocate (im1, 0, 192, 192); red = gdImageColorAllocate (im1, 255, 0, 0); yellow = gdImageColorAllocate (im1, 255, 255, 0); green = gdImageColorAllocate (im1, 0, 255, 0); blue = gdImageColorAllocate (im1, 0, 0, 255); magenta = gdImageColorAllocate (im1, 192, 0, 192); trans = gdImageColorAllocate (im1, 1, 1, 1); col01 = gdImageColorAllocate (im1, 255, 192, 0); col02 = gdImageColorAllocate (im1, 192, 255, 0); col03 = gdImageColorAllocate (im1, 128, 255, 128); col04 = gdImageColorAllocate (im1, 0, 255, 192); col05 = gdImageColorAllocate (im1, 0, 192, 255); col06 = gdImageColorAllocate (im1, 128, 128, 255); col07 = gdImageColorAllocate (im1, 192, 0, 255); col08 = gdImageColorAllocate (im1, 255, 0, 192); col09 = gdImageColorAllocate (im1, 255, 128, 128); gifout = fopen ( GIFOUT , "wb"); } // ----------------------------------------------------------------- // describe toroidal orbits double t_ord[25]; double tz[25]; int tc[25]; double td[25]; double ta[25]; void torstart() { tc[ 0] = white; td[ 0] = 0.0000; ta[ 0] = 0.000; tc[ 1] = yellow; td[ 1] = 1.0000; ta[ 1] = 0.000; tc[ 2] = blue; td[ 2] = 1.0000; ta[ 2] = 180.000; tc[ 3] = red; td[ 3] = 1.0000; ta[ 3] = 90.000; tc[ 4] = green; td[ 4] = 1.0000; ta[ 4] = 270.000; tc[ 5] = col01; td[ 5] = 1.4142; ta[ 5] = 45.000; tc[ 6] = col02; td[ 6] = 1.4142; ta[ 6] = 135.000; tc[ 7] = col03; td[ 7] = 1.4142; ta[ 7] = 225.000; tc[ 8] = col04; td[ 8] = 1.4142; ta[ 8] = 315.000; tc[ 9] = col05; td[ 9] = 2.0000; ta[ 9] = 0.000; tc[10] = col06; td[10] = 2.0000; ta[10] = 180.000; tc[11] = col07; td[11] = 2.0000; ta[11] = 90.000; tc[12] = col08; td[12] = 2.0000; ta[12] = 270.000; tc[13] = col09; td[13] = 2.2361; ta[13] = 26.565; tc[14] = magenta; td[14] = 2.2361; ta[14] = 206.565; tc[15] = cyan; td[15] = 2.2361; ta[15] = 116.565; tc[16] = gray; td[16] = 2.2361; ta[16] = 296.565; tc[17] = dgray; td[17] = 2.2361; ta[17] = 333.435; tc[18] = col01; td[18] = 2.2361; ta[18] = 153.435; tc[19] = col02; td[19] = 2.2361; ta[19] = 63.349; tc[20] = col03; td[20] = 2.2361; ta[20] = 243.435; tc[21] = col04; td[21] = 2.8284; ta[21] = 45.000; tc[22] = col05; td[22] = 2.8284; ta[22] = 135.000; tc[23] = col06; td[23] = 2.8284; ta[23] = 225.000; tc[24] = col07; td[24] = 2.8284; ta[24] = 315.000; } // ----------------------------------------------------------------- void displayend() { gdImageGifAnimEnd(gifout); fclose (gifout); gdImageDestroy (im1); }