// cap11.c Capture rail GEO // ------------------------------------------------------------------------ // GEO, capture rail 24 hours of simulation // ------------------------------------------------------------------------ // cc -o cap11 cap11.c -lgd -lpng -lm // ------------------------------------------------------------------------ // Uses the libgd library, documentation: /usr/share/doc/gd-*/index.html // // The image is drawn large, then resampled down with ImageMagick "convert" // Slow, but resulting in better looking output. // // png2swf segmentation faults if given too many large frames #define NAME "cap11" #define TITLETEXT "Capture Rail" #define DEBUG 1 // earth and orbits ----------------------------------------------- // we will use a 24 hour, not a 23.9346 hour, sidereal day. // #define CAPTURE 30411.507 // climbing capture #define SCALE 3400 // bigger for less magnification #define OTIME 24.0 // orbit time #define RAILHIGH 45000.0 #define RAILLOW 28000.0 #define NSTAR 100 #define RE 6378.0 #define ETIME 24.0 // hours per turn #define MU 398600.0 // grav parameter km3/s2 #define ALTITUDE 80.0 // perigee altitude #define GMAX 0.003 // 3 m/s^2, 0.003km/s^2 // plot rate ------------------------------------------------------ #define RATE 5 #define NPLOT 120 #define MAXHOUR 24.0 #define STEPS 1000 // orbit simulation steps per plot // plot sizing ----------------------------------------------------- #define PSCALE 32 // global scaling factor for drawing #define LINSCALE 1 #define XSIZE (96*PSCALE) #define YSIZE (32*PSCALE) #define XOUT 960 // size of output image after convert #define YOUT 320 #define XCENT1 (16*PSCALE) // #1 fixed frame orbit #define YCENT1 (18*PSCALE) #define MAG1 1.0 #define XCENT2 (48*PSCALE) // #2 rotating frame orbit #define YCENT2 (18*PSCALE) #define MAG2 1.0 #define XCENT3 (80*PSCALE) // #3 magnified rotating frame orbit #define YCENT3 (58*PSCALE) #define MAG3 4.0 #define VSIZE (PSCALE/2) // Vehicle circle size #define FNT "DejaVuMonoSans" #define FSZ (1.25*PSCALE) #include #include #include #include #include //----------------------------------------------------------------------- int white, lgray, gray, dgray ; // color indexes int black, blue, cyan, dcyan ; int red, dred, green, dgreen ; int yellow, dyellow ; // elapsed time long start_time = 0L ; // epoch start time int end_time = 0 ; // estimated end time int elapsed_time = 0 ; // elapsed time int remaining_time = 0 ; // remaining time double speedup ; // simulation to display time double pi2 = 6.28318530717958646 ; double third = 0.33333333333333 ; // files and graphics arrays gdImagePtr im ; // Declare the image FILE *pngout ; // Declare output file char command[200] ; // char filename[80] ; // // simulation time parameters int iplot ; // number of plot double ip1 ; // fraction of complete plot double hour ; // time of plot in hours from start double eangle ; // earth turn angle double oangle ; // destination orbit turn angle //----------------------------------------------------------------------- // orbit parameters // distances in km // speeds in km/sec // target orbit parameters double omega ; // rad/sec angular velocity double orbsize ; // km orbit radius // transfer orbit parameters double rp ; // km perigee radius double ra ; // km apogee radius double rc ; // km capture radius, apogee double vc ; // km/s capture velocity double sem ; // km semimajor axis double e ; // eccentricity double h ; // km^2/s angular momentum double vp ; // km/s perigee velocity double vct ; // km/s transverse capture velocity double v0 ; // km/s orbit velocity double r0 ; // km/s orbit radius double costh ; // cos of orbit angle at capture double sinth ; // sin of orbit angle at capture double theta ; // orbit angle radians double dtheta ; // orbit angle degrees double period ; // orbit period, seconds double tr_time ; // orbit transfer time, seconds double arg_p ; // argument of perigee double energy ; // orbit energy, MJ double enc ; // energy at capture, MJ double av ; // double rv ; // double vv ; // double tstep ; // simulation timestep // background stars ; double starx[NSTAR], stary[NSTAR] ; // star field position // transfer orbit simulation radius (km) and angle (radians), fixed frame double cr[NPLOT], ca[NPLOT] ; // transfer orbit, with capture double mr[NPLOT], ma[NPLOT] ; // transfer orbit, missed capture double rr[NPLOT] ; // transfer orbit, rotation int nc = NPLOT ; // capture step int nv = 0 ; // number of transfer orbit points int scale = SCALE/PSCALE ; // pixels per kilometer // ======================================================================== void initialize() { int i ; double tmp ; double E ; double M ; sprintf( command, "rm -rf %sdir ; mkdir %sdir", NAME, NAME ); system( command ) ; // make directory start_time = (long) time( NULL ); speedup = 3600.0 * MAXHOUR * RATE / NPLOT ; // sim/display time ratio tstep = 3600.0*MAXHOUR/(NPLOT*STEPS) ; // simulation timestep // orbit initialize // fake: sidereal 24 hour day omega = pi2 / ( 3600.0 * OTIME ) ; // angular velocity orbsize = pow( ( MU /(omega*omega)), third ) ; // orbit radius rp = RE + ALTITUDE ; // perigee radius // three different possibilities - #ifdef CAPTURE // ---------------- new capture rail rc = CAPTURE ; // see http://www.nature1st.net/bogan/orbits/kepler/orbteqtn.html // though be careful of the redefinition of q, both perigee and // true anomaly vct = rc*omega ; // transverse capture velocity h = rc*vct ; // orbit angular momentum vp = h/rp ; // perigee velocity km/sec v0 = MU/h ; // orbit velocity r0 = h/v0 ; // orbit radius e = (r0/rp)-1.0 ; // eccentricity ra = r0/(1.0-e) ; // apogee radius sem = r0/(1.0-e*e) ; // semimajor axis energy = 0.5*vp*vp-MU/rp ; // orbit energy MJ enc = 0.5*vct*vct-MU/rc ; // vc = sqrt(2.0*(energy-enc)); // tmp = fabs(sem*sem*sem)/MU ; // period = pi2*sqrt( tmp ) ; // orbit period costh = ((r0/rc)-1.0)/e ; sinth = sqrt(1.0-costh*costh) ; theta = atan2(sinth,costh) ; // true anomaly dtheta = (360.0/pi2)*theta ; av = vct*vct/rc-MU/(rc*rc) ; // acceleration at capture if( sem > 0.0 ) { // ----- elliptical double cosE = (e+costh)/(1.0+e*costh); // possibly two solutions for E, M, and tr_time ! E = acos( cosE ) ; // eccentric anomaly M = E - e*sin( E ) ; // mean anomaly } else { // ----- hyperbolic double coshE = (e+costh)/(1.0+e*costh); // possibly two solutions for E, M, and tr_time ! E = acosh( coshE ) ; // eccentric anomaly M = e*sinh( E ) - E ; // mean anomaly } #else // --------------------- classical capture rc = 0.7*orbsize ; tmp = 2.0 * MU / (omega*omega); for( i = 0 ; i < 77 ; i++ ) { rc = pow( tmp*rp/(rc+rp), third);// } sem = 0.5*(rp+rc) ; // semimajor axis e = 1.0 - rp/sem ; // eccentricity vct = rc*omega ; // apogee (capture) velocity h = rc*vct ; // orbit angular momentum v0 = MU/h ; // orbit velocity r0 = h/v0 ; // orbit radius ra = rc ; // apogee vp = h/rp ; // perigee velocity km/sec energy = 0.5*vp*vp-MU/rp ; // orbit energy MJ enc = 0.5*vct*vct-MU/rc ; // vc = sqrt(2.0*(energy-enc)); // radial capture velocity tmp = fabs(sem*sem*sem)/MU ; // period = pi2*sqrt( tmp ) ; // orbit period costh = -1.0 ; // sinth = 0.0 ; // theta = 0.5*pi2 ; // dtheta = 180.0 ; // E = 0.5*pi2 ; // M = E ; // av = vct*vct/rc-MU/(rc*rc) ; // acceleration at capture #endif tr_time = M*period / pi2 ; // time from perigee, seconds arg_p = omega*tr_time-theta ; // argument of perigee #ifdef DEBUG printf( "|omega %11.3e rd/s ", omega ); printf( "|orbsize%11.2f km ", orbsize ); printf( "|tstep %11.5f sec ", tstep ); printf( "|\n" ); printf( "|rp %11.2f km ", rp ); printf( "|ra %11.2f km ", ra ); printf( "|sem %11.2f km ", sem ); printf( "|\n" ); printf( "|e %11.6f ", e ); printf( "|vp %11.4f km/s ", vp ); printf( "|vct %11.4f km/s ", vct ); printf( "|\n" ); printf( "|h %11.2fkm2/s ", h ); printf( "|r0 %11.2f km ", r0 ); printf( "|v0 %11.4f km/s ", v0 ); printf( "|\n" ); printf( "|mu %11.2fm3/s2 ", MU ); printf( "|energy %11.4fMJ/kg ", energy ); printf( "|vc %11.4f km/s ", vc ); printf( "|\n" ); printf( "|rc %11.2f km ", rc ); printf( "|E %11.4f rad ", E ); printf( "|M %11.4f rad ", M ); printf( "|\n" ); printf( "|tr_time%11.2f sec ", tr_time ); printf( "|period %11.4f sec ", period ); printf( "|p.hours%11.4f hr ", period/3600.0 ); printf( "|\n" ); printf( "|theta %11.4f rad ", theta ); printf( "|arg_p %11.6f rad ", arg_p ); printf( "|av_c %11.4f m/s2 ", 1000.0*av ); printf( "|\n" ); #endif // DEBUG // start simulations at perigee cr[0] = rp ; // captured radius mr[0] = rp ; // missed capture radius ca[0] = 0.0 ; // captured angle ma[0] = 0.0 ; // missed capture angle rr[0] = 0.0 ; // rotation at capture // --------------------------------------------- initialize stars // generate a disk of stars with // radius < 1.414*orbsize double stmult = orbsize*sqrt(8.0); for( i=0 ; i < NSTAR ; ) { double x = drand48()-0.5; double y = drand48()-0.5; if( x*x+y*y < 0.25 ) { starx[i]=stmult*x ; stary[i]=stmult*y ; i++ ; } } } //----------------------------------------------------------------------- void openplot() { im = gdImageCreateTrueColor( (int)XSIZE, (int)YSIZE ); black = gdImageColorAllocate(im, 0, 0, 0); dgray = gdImageColorAllocate(im, 96, 96, 96); gray = gdImageColorAllocate(im, 128, 128, 128); lgray = gdImageColorAllocate(im, 160, 160, 160); white = gdImageColorAllocate(im, 255, 255, 255); red = gdImageColorAllocate(im, 255, 0, 0); dred = gdImageColorAllocate(im, 160, 0, 0); green = gdImageColorAllocate(im, 0, 255, 0); dgreen = gdImageColorAllocate(im, 0, 160, 0); cyan = gdImageColorAllocate(im, 0, 160, 160); dcyan = gdImageColorAllocate(im, 0, 80, 80); yellow = gdImageColorAllocate(im, 255, 255, 0); dyellow = gdImageColorAllocate(im, 160, 160, 0); blue = gdImageColorAllocate(im, 0, 0, 255); gdImageSetThickness( im, LINSCALE ); #if DEBUG > 1 printf("\n---------------------------------------------------\n"); #endif // DEBUG } //----------------------------------------------------------------------- // draw orbiting GEO tether // ex, ey are center of drawing in pixels // ang and rot in radians void tether( double ang, double rot, double mag, int ex, int ey ) { int x0, y0 ; int x1, y1 ; double rang = ang+rot; double rlow = mag*RAILLOW /scale ; double rhigh = mag*RAILHIGH/scale ; int orb = 2*(int)(mag*orbsize/scale) ; int wide = (7*orb)/10 ; // Draw a black box over background. This covers up overlap // from other views. gdImageFilledRectangle( im, ex-wide, 0, ex+wide, YSIZE, black ); // draw orbit circle gdImageSetThickness( im, (int)5.0*mag*LINSCALE ); gdImageArc( im, ex, ey, orb, orb, 0, 360, dgray ); // draw tether gdImageSetThickness( im, (int)3.0*mag*LINSCALE ); x0 = ex + (int)( rlow * sin(rang)); y0 = ey - (int)( rlow * cos(rang)); x1 = ex + (int)( rhigh * sin(rang)); y1 = ey - (int)( rhigh * cos(rang)); gdImageLine( im, x0, y0, x1, y1, white ); } //----------------------------------------------------------------------- // draw stars void stars( double rot, double mag, int ex, int ey ) { int x, y, i ; double s = sin( rot ); double c = cos( rot ); double ms = mag/scale ; int size = 4*(int)sqrt(mag) ; int wide = (int)(1.3*ms*orbsize) ; gdImageSetThickness( im, 1); for( i=0 ; i < NSTAR ; i++ ) { int x = (int)(-ms*(c*starx[i] + s*stary[i])); if( abs( x ) < wide ) { x += ex ; int y = ey + (int)(ms*(c*stary[i] - s*starx[i])); gdImageFilledArc(im, x, y, size, size, 0, 360, white, gdArc); } } } //----------------------------------------------------------------------- // draw rotating earth // ex, ey are center of drawing in pixels // ang and rot in radians void earth( double ang, double rot, double mag, int ex, int ey ) { int x0, y0 ; int x1, y1 ; double an0 ; int i ; double r2d = 360.0/pi2 ; double drot = r2d*rot + 360000.0 ; // make sure it is positive double rang = ang+rot; double re = mag*RE/scale ; int esize = 2*((int)re) ; double ea1 = fmod( drot + 180.0, 360.0 ) ; double ea2 = fmod( drot , 360.0 ) ; gdImageSetThickness( im, 1); gdImageFilledArc(im, ex, ey, esize, esize, 0, 360, cyan, gdArc); gdImageFilledArc(im, ex, ey, esize, esize, ea2, ea1, dcyan, gdArc); // draw longitude lines ( 6 = 180/30 degrees ) gdImageSetThickness( im, LINSCALE ); for( i=5 ; i >= 0 ; i-- ) { an0 = rang + pi2*( (double) i ) / 12.0 ; x0 = ex + (int)( re * sin(an0 )); y0 = ey - (int)( re * cos(an0 )); x1 = ex - (int)( re * sin(an0 )); y1 = ey + (int)( re * cos(an0 )); gdImageLine( im, x0, y0, x1, y1, black ); } gdImageSetThickness( im, 3*LINSCALE ); gdImageLine( im, x0, y0, ex, ey, black ); // draw launch site pointer x0 = ex + (int)( re * sin( rang + arg_p )); y0 = ey - (int)( re * cos( rang + arg_p )); x1 = ex + (int)( 0.5*re * sin( rang + arg_p )); y1 = ey - (int)( 0.5*re * cos( rang + arg_p )); gdImageSetThickness( im, 6*LINSCALE ); gdImageLine( im, x0, y0, x1, y1, yellow ); // draw latitude lines gdImageSetThickness( im, LINSCALE ); for( i=1 ; i < 3 ; i++ ) { an0 = pi2*( (double) i ) / 12.0 ; x0 = (int)( 2.0 * ((double)re) * cos(an0) ); gdImageArc( im, ex, ey, x0, x0, 0, 360, black ); } } //----------------------------------------------------------------------- // draw vehicle transfer orbit void trans( int rotflag, double rot, double mag, int ex, int ey ) { int i ; int x, y ; if( nv < iplot ) { // ------------ compute new transfer orbit point double th0 = ma[nv] ; #if DEBUG > 1 printf( "%5dnv%6dnc%6diplot%16.8fth0\n", nv, nc, iplot, th0 ); printf( "theta%9.4f rot%9.4f\n", theta, rot ); #endif // DEBUG for( i = 0; i < STEPS ; i++ ) { double rx = r0/(1.0+e*cos(th0)) ; // double rx = r0/(1.0+e*cos(th0)) + 0.5*v0*e*tstep*sin(th0); th0 += h*tstep/(rx*rx); #if DEBUG > 1 printf( "rx%12.4e th0%10.6f\n", rx, th0 ); #endif // DEBUG if( nc != NPLOT ) { #ifdef CAPTURE double atest = 0.5*vv*vv/fabs(orbsize-rv) ; if( atest < GMAX ) { av = omega*omega*rv-MU/(rv*rv) ; } else { av = -GMAX ; } vv += tstep*av ; if( ( vv < 0.0 ) || ( rv > orbsize ) ) { vv=0.0 ; } rv += tstep*vv ; #endif // CAPTURE } else if( th0 > theta ) { nc = iplot ; rv = r0/(1.0+e*cos(th0)) ; vv = v0*e*sin(th0) ; #if DEBUG > 0 printf("\nCapture nc=%04d rv=%9.2f vv=%7.4f\n", nc, rv, vv ); #endif } } nv++ ; rr[nv] = rot ; mr[nv] = r0/(1.0+e*cos(th0)); ma[nv] = th0 ; if( nc == NPLOT ) { cr[nv] = mr[nv]; ca[nv] = ma[nv]; } else { cr[nv] = rv ; ca[nv] = oangle-arg_p; #if DEBUG > 0 printf("%8.4f%7.3f%9.1f", 1000.0*av, vv, rv ); #endif } } // plot transfer orbit points double man0 ; int mx ; int my ; double can0 ; int cx ; int cy ; double ms = mag/scale ; int vsize = 2*((int)(0.5*mag*VSIZE)) ; double ar ; for( i = 0 ; i <= nv ; i++ ) { if( rotflag == 1 ) { ar = arg_p + rr[i]; } else { ar = arg_p + rot; } man0 = ma[i] + ar ; mx = ex + (int)(ms*mr[i]*sin(man0)); my = ey - (int)(ms*mr[i]*cos(man0)); #if DEBUG > 1 printf( "%4dnv%4dnc%4diplot%4di%5dmx%5dmy%8.1f_r%9.4f_an0\n", nv, nc, iplot, i, mx, my , mr[i], man0 ); #endif // DEBUG if( i <= nc ) { gdImageFilledArc( im, mx, my, vsize, vsize, 0, 360, dyellow, gdArc ); } else { gdImageFilledArc( im, mx, my, vsize, vsize, 0, 360, dred , gdArc ); can0 = ca[i] + ar ; cx = ex + (int)(ms*cr[i]*sin(can0)); cy = ey - (int)(ms*cr[i]*cos(can0)); gdImageFilledArc( im, cx, cy, vsize, vsize, 0, 360, dgreen , gdArc ); } } // current points if( iplot <= nc ) { gdImageFilledArc( im, mx, my, vsize, vsize, 0, 360, yellow, gdArc ); } else { gdImageFilledArc( im, mx, my, vsize, vsize, 0, 360, red , gdArc ); gdImageFilledArc( im, cx, cy, vsize, vsize, 0, 360, green , gdArc ); } } //----------------------------------------------------------------------- void closeplot() { sprintf( command, "%s.png", NAME ); pngout = fopen( command, "wb"); gdImagePngEx(im, pngout, 1); fclose(pngout); // write out image gdImageDestroy(im); sprintf( command, "convert %s.png -resize %dx%d\\> %sdir/b%04d.png", NAME, XOUT, YOUT, NAME, iplot ); system( command ); } // ======================================================================== int main() { initialize() ; for( iplot = 0 ; iplot < NPLOT ; iplot++ ) { openplot() ; ip1 = ((double)iplot)/((double) NPLOT) ; hour = ip1*MAXHOUR ; eangle= pi2 * hour/ETIME ; // radians oangle= pi2 * hour/OTIME ; // radians tether( oangle, -oangle, MAG3, XCENT3, YCENT3 ) ; stars( -oangle, MAG3, XCENT3, YCENT3 ) ; earth( eangle, -oangle, MAG3, XCENT3, YCENT3 ) ; trans( 1, -oangle, MAG3, XCENT3, YCENT3 ) ; // note - the next two drawings overlap and obscure the // first magnified plot tether( oangle, -oangle, MAG2, XCENT2, YCENT2 ) ; stars( -oangle, MAG2, XCENT2, YCENT2 ) ; earth( eangle, -oangle, MAG2, XCENT2, YCENT2 ) ; trans( 1, -oangle, MAG2, XCENT2, YCENT2 ) ; tether( oangle, 0.0, MAG1, XCENT1, YCENT1 ) ; stars( 0.0, MAG1, XCENT1, YCENT1 ) ; earth( eangle, 0.0, MAG1, XCENT1, YCENT1 ) ; trans( 0, 0.0, MAG1, XCENT1, YCENT1 ) ; // labels ----------------------------------------------------------- gdImageStringFT( im, NULL, // imagespace, bounding white, FNT, 1.6*FSZ, // color, font, fontsize 0.0, 1*PSCALE, 3*PSCALE, // angle, x, y TITLETEXT ); // text sprintf( command, "View from South %04.1f/%04.1f hours, %4.0fx speedup", hour, MAXHOUR, speedup ); gdImageStringFT( im, NULL, // imagespace, bounding white, FNT, FSZ, // color, font, fontsize 0.0, 44*PSCALE, 3*PSCALE, // angle, x, y command ); // text gdImageStringFT( im, NULL, // imagespace, bounding white, FNT, 0.75*FSZ, // color, font, fontsize 0.0, 1*PSCALE, 31*PSCALE, // angle, x, y "Fixed Frame" ); // text gdImageStringFT( im, NULL, // imagespace, bounding white, FNT, 0.75*FSZ, // color, font, fontsize 0.0, 32*PSCALE, 31*PSCALE, // angle, x, y "Rotating Frame" ); // text gdImageStringFT( im, NULL, // imagespace, bounding white, FNT, 0.75*FSZ, // color, font, fontsize 0.0, 64*PSCALE, 31*PSCALE, // angle, x, y "Rotating Frame, 4x larger" );// text closeplot(); elapsed_time = (int)(((long) time( NULL )) - start_time ) ; if( elapsed_time > 1 ) { end_time = (long)(((double) elapsed_time ) / ip1 ); remaining_time = end_time - elapsed_time ; } #if DEBUG < 2 printf( "%4d/%d-plot%4d-nc%5.1f/%04.1f-hour%4ld/%4ld/%4ldsec\r", iplot+1, NPLOT, nc, hour, ((double)MAXHOUR), elapsed_time, remaining_time, end_time ); fflush(stdout); #endif // DEBUG } printf( "\n" ); fflush(stdout); sprintf( command, "png2swf -o %s.swf -r%3d %sdir/*.png", NAME, RATE, NAME ); system( command ); return(0); }