/* capture22.c KHL 2019 Feb 21 gcc -o capture22 capture22.c -lm ; ./capture22 > capture22.out */ #define NAME "capture22" /* ========================================================================= Times are relative to construction orbit apogee, hence negative eccentricity. No correction for apsidal precession. This starts with the [tstart] launch time, and uses an iterative loop to find the best CAN_A extra apogee value, Rather than try to be too clever, the search will start with a lower bound of CAN_A_min = [1 km/s] * [tstart] and a starting upper bound of CAN_A_max = [2 km/s] * ([tstart] + 2 [s] ). The midpoint of the two bounds is computed as the can_a value, which is evaluated for the error time, the difference between the construction and the candidate orbit times at the geometric orbit crossing If the error time is negative, the lower bound is replaced with can_a . If the error time is positive, the upper bound is replaced with can_a . This is a fast calculation, so we will just iterate 60 times and let the bounds converge. Due to numerical noise, the calculation will be non-monotonic for small errors, but we will be close to the "actual" value. In a real system, other perturbations (thruster error, tidal effects from the equatorial bulge and the Moon, Sun, and Jupiter, and general mechanical cussedness) will be much more important, and require small midflight correction thrusts. Below [tstart] of 1.0 seconds, the calculations become erratic; search algorithms should lower-bounds-check at the projected radius of the construction orbit at the [tstart] delay. The curvature of the candidate orbit with that apogee will always be less than the construction orbit, so there cannot be a lower but intersecting candidate orbit. Add a small epsilon just to make sure the calculation is bounded. If an intersection is not found, OR the error time cannot be made to converge to zero, skip that table entry and move on. Set printme = 1 at the start of the loop, and set it to zero if either the binomial root is imaginary or the loop does not converge. [tstart] CAN_A 0.003 3.0646 from 3.05 to 1e9 0.01 3.084670 from 0.19 to 1e9 0.03 3.14398 from 0.00 to 1e9 0.1 3.34698 from 0.00 to 1e9 0.3 3.899743 from 0.00 to 1e9 1.0 5.64861821 from 0.00 to 1e9 3.0 10.02098638 from 0.00 to 1e9 10.0 23.45589910 from 0.00 to 1e9 30.0 58.49035158 from 0.00 to 1e9 45.0 83.802205003 from 0.00 to 1e9 60.0 108.75943134887 from 0.00 to 1e9 90.0 158.08913293178 from 0.00 to 1e9 180.0 304.03736662938 from 0.00 to 1e9 300.0 496.931100566385 from 0.00 to 1e9 450.0 737.308802551446 from 0.00 to 1e9 600.0 977.727815882209 from 0.00 to 1e9 1200.0 1947.92465854122 1800.0 2946.23116529545 2400.0 3992.32125234230 from 0.00 to 1e9 3000.0 5108.68173943019 3600.0 6321.18761716776 from 0.00 to 1e9 ======================================================================*/ #define DIAG "%s.diag" // diagnostic output file #define WIKI "%s.wiki" // main wiki output file #define SUM "%s.sum" // main summary file #define TLSTEP 100.0 // sec launch time step #define TLMIN -900.0 // sec launch time start #define TLMAX 900.0 // sec launch time range, minus to plus #define LOOPA 80.0 // km loop altitude #define LOOPL -8.0 // degrees south loop latitude #define CON_PA 2000.0 // km construction orbit perigee altitude #define CON_SD 1.0 // sidereal day construction orbit period #define MU0 398600.4418 // km3/s2 Earth standard gravitational param #define RE 6378.0 // km/s equatorial radius #define SOL 86400.0 // sec solar day #define YEAR 365.2422 // solar days per year #define NUMCAN 60 // number of convergence iterations #define EPS 1e-5 // crossing window #include #include //======================================================================= //======================================================================= // wiki output subroutines ----------------------------------------- //======================================================================= //======================================================================= // print wiki header ------------ void wikiheader( FILE * fr, double p_constr, double p_prime, double p_sday , double p_secs ) { // first header line fprintf( fr, "||<-7> %s ", NAME ); fprintf( fr, "construction perigee =%5.0f km", p_constr ); fprintf( fr, " launch perigee =%5.0f km", p_prime ); fprintf( fr, " ||\n" ); // second header line fprintf( fr, "||<-7> " ); fprintf( fr, "sidereal period: %2.0f days,", p_sday ); fprintf( fr, "%10.3f seconds ", p_secs ); fprintf( fr, " ||\n" ); // third header line fprintf( fr, "|| " ); // name of parameter fprintf( fr, "|| period " ); // orbit period sec fprintf( fr, "|| arrival " ); // arrival time sec fprintf( fr, "|| apogee " ); // apogee km fprintf( fr, "||<-3> velocity change m/s " ); fprintf( fr, "||\n" ); // end of first header line // fourth header line fprintf( fr, "|| " ); // name of parameter fprintf( fr, "|| sec " ); // orbit period sec fprintf( fr, "|| sec " ); // arrival time sec fprintf( fr, "|| km " ); // apogee km fprintf( fr, "|| tangent" ); fprintf( fr, "|| radial " ); fprintf( fr, "|| plane " ); fprintf( fr, "||\n" ); // end of second header line } // end of wiki header // print blank wiki line ------------- void blankwiki( FILE * fr ) { fprintf( fr, "||<-7> "); fprintf( fr, " ||\n"); } // print line of wiki data --------- void printwiki( FILE * fo, char* p_name, double p_period, double p_launch, double p_arrival, double p_apogee, double p_perigee, double p_vt, double p_vr, double p_vp ) { fprintf( fo, "|| %s" , p_name ); // name of parameter fprintf( fo, "||%10.3f " , p_period ); // orbit period sec fprintf( fo, "||%9.3f " , p_arrival ); // arrival time sec fprintf( fo, "||%10.3f " , p_apogee ); // apogee km fprintf( fo, "||%7.2f " , 1000.0*p_vt ); // tangential m/s fprintf( fo, "||%8.2f " , 1000.0*p_vr ); // radial m/s fprintf( fo, "||%6.2f " , 1000.0*p_vp ); // plane change m/s fprintf( fo, "||\n" ); // end of table line } // end of wiki print line //======================================================================= // print wiki, output or diagnostic header void printhead( FILE * fo, char* p_name, double p_period, double p_launch, double p_arrival, double p_apogee, double p_perigee, double p_vt, double p_vr, double p_vp ) { fprintf( fo, "|| %s" , p_name ); // name of parameter fprintf( fo, "||%10.3f " , p_period ); // orbit period sec fprintf( fo, "||%9.3f " , p_arrival ); // arrival time sec fprintf( fo, "||%10.3f " , p_apogee ); // apogee km fprintf( fo, "||%7.2f " , 1000.0*p_vt ); // tangential m/s fprintf( fo, "||%8.2f " , 1000.0*p_vr ); // radial m/s fprintf( fo, "||%6.2f " , 1000.0*p_vp ); // plane change m/s fprintf( fo, "||\n" ); // end of table line } // end of printhead; output or diagnostic header //======================================================================= // print construction and prime launch orbit, to output or diagnostic void printorbits( FILE * fo, double con_time, double con_v0, double con_semi, double con_v_peri, double con_peri, double con_v_apo, double con_apo, double con_e, // negative double pri_time, double pri_ltime, double pri_v0, double pri_semi, double pri_v_peri, double pri_peri, double pri_v_apo, double pri_apo, double pri_e, // negative double pri_v_loop, double pri_v_ins ) { fprintf( fo, "\n construction orbit ----------" ); fprintf( fo, "--------------------------------------------\n" ); fprintf( fo, "%14.3f km semimajor axis", con_semi ); fprintf( fo, "%14.6f km/s V0 velocity\n", con_v0 ); fprintf( fo, "%14.3f km perigee", con_peri ); fprintf( fo, "%21.6f km/s perigee velocity\n", con_v_peri ); fprintf( fo, "%14.3f km apogee", con_apo ); fprintf( fo, "%22.6f km/s apogee velocity\n", con_v_apo ); fprintf( fo, "%14.3f sec orbit period", con_time ); fprintf( fo, "%16.6f eccentricity\n", con_e ); fprintf( fo, "\n prime launch orbit ----------" ); fprintf( fo, "--------------------------------------------\n" ); fprintf( fo, "%14.3f km semimajor axis", pri_semi ); fprintf( fo, "%14.6f km/s V0 velocity\n", pri_v0 ); fprintf( fo, "%14.3f km perigee", pri_peri ); fprintf( fo, "%21.6f km/s perigee velocity\n", pri_v_peri ); fprintf( fo, "%14.3f km apogee", pri_apo ); fprintf( fo, "%22.6f km/s apogee velocity\n", pri_v_apo ); fprintf( fo, "%14.3f sec orbit period", pri_time ); fprintf( fo, "%16.6f eccentricity\n", pri_e ); fprintf( fo, "%14.3f sec launch time\n", pri_ltime ); fprintf( fo, "%14.6f km/s insertion velocity", pri_v_ins ); fprintf( fo, "%10.5f km/s loop exit velocity\n", pri_v_loop ); } //======================================================================= // print most of data to diagnostic or summary file /// void printmost{ FILE * fo, double tstart, /// double can_angle, double can_v_loop, /// double can_time, double can_ltime, double can_e, /// double can_semi, double can_peri, double can_apo, /// double can_v0, double can_v_peri, double can_v_apo, /// double rcon, double rcan, /// double theta1, double theta2, // double theta, double thetac, /// double omega_con, double omega_can, /// double E_con, double E_can, /// double M_con, double M_can, /// double delay_con, double delay_can, /// double pri_time, double apogee_delay, double can_arrival, /// double v_canr, double v_cant, /// double v_conr, double v_cont, /// double error_time, double can_a_min, double can_a_max } void printmost( FILE * fo, double tstart, double can_angle, double can_v_loop, double can_time, double can_ltime, double can_e, double can_semi, double can_peri, double can_apo, double can_v0, double can_v_peri, double can_v_apo, double rcon, double rcan, double theta1, double theta2, double theta, double thetac, double omega_con, double omega_can, double E_con, double E_can, double M_con, double M_can, double delay_con, double delay_can, double pri_time, double apogee_delay, double can_arrival, double v_canr, double v_cant, double v_conr, double v_cont, double error_time, double can_a_min, double can_a_max ) { double r2d = 45.0/atan(1.0) ; double degree = r2d * theta ; double degree1 = r2d * theta1 ; double degree2 = r2d * theta2 ; double degreec = r2d * thetac ; double canad = r2d * can_angle ; fprintf( fo, "\n candidate launch orbit -------------" ); fprintf( fo, "-------------------------------------\n" ); fprintf( fo, "%14.3f km semimajor axis", can_semi ); fprintf( fo, "%14.6f km/s V0 velocity\n", can_v0 ); fprintf( fo, "%14.3f km perigee", can_peri ); fprintf( fo, "%21.6f km/s perigee velocity\n", can_v_peri ); fprintf( fo, "%14.3f km apogee", can_apo ); fprintf( fo, "%22.6f km/s apogee velocity\n", can_v_apo ); fprintf( fo, "%14.3f sec orbit period", can_time ); fprintf( fo, "%16.6f eccentricity\n", can_e ); fprintf( fo, "%14.3f sec launch time", can_ltime ); fprintf( fo, "%+17.3f sec start time\n", tstart ); fprintf( fo, "%14.3f km/s loop velocity ", can_v_loop ); if( can_angle < 0.0 ) { fprintf( fo, "%+13.6f rad,%7.4f° west of prime\n\n", can_angle, -canad ); } else { fprintf( fo, "%+13.6f rad,%7.4f° east of prime\n\n", can_angle, canad ); } fprintf( fo, "theta1=%12.7f theta2=%12.7f rad\n", theta1, theta2 ); fprintf( fo, " "); fprintf( fo, " construction candidate\n"); fprintf( fo, "theta:%45.7f%16.7f rad\n", theta, thetac ); fprintf( fo, "radius:%44.5f%16.5f km\n", rcon, rcan ); fprintf( fo, "omega:%45.5e%16.5e rad/sec\n", omega_con, omega_can ); fprintf( fo, "eccentric anomaly E:%31.7f%16.7f rad\n", E_con, E_can ); fprintf( fo, "mean anomaly M:%36.7f%16.7f rad\n", M_con, M_can ); fprintf( fo, "time after apogee:" ); fprintf( fo, " %18.6f%16.6f sec\n", delay_con, delay_can ); fprintf( fo, "\n%+8.5f rad,%+9.4f° candidate orbit launch angle\n", thetac, degreec ); fprintf( fo, "%12.6f sec construction delay from apogee to crossing\n", delay_con ); fprintf( fo, "%12.6f sec candidate delay from apogee to crossing\n", delay_can ); fprintf( fo, "%12.6f sec candidate orbit period\n", can_time ); fprintf( fo, "%12.6f sec prime orbit period\n", pri_time ); fprintf( fo, "%12.6f sec candidate apogee difference\n", apogee_delay ); fprintf( fo, "%12.6f sec candidate orbit arrival\n", can_arrival ); fprintf( fo, "v_canr%9.5f v_cant%9.5f v_conr%9.5f v_cont%9.5f km/s\n", v_canr, v_cant, v_conr, v_cont ); fprintf( fo, "error_time=%15.4e sec can_apogee: min=%12.5f km max=%12.5f km\n", error_time, can_a_min, can_a_max ); } // ======================================================================= // print header for diagnostic and summary file void fileheader( FILE * fo, double vloop ) { fprintf(fo, "============ /home/keithl/launchloop/capture/%s.c =============\n", NAME); fprintf(fo, "=================%2.0f Sidereal Day Construction Orbit =================\n", CON_SD ); fprintf(fo, "\n Launch loop at %3.0f km altitude %3.0f°S latitude\n", LOOPA, LOOPL ); fprintf(fo, "%14.6f km/s loop latitude rotation velocity\n", vloop ); } //======================================================================== //======================================================================== // main program --------------------------------------------------- //======================================================================== //======================================================================== int main(void) { double pi = 4.0 * atan(1.0) ; // double pi2 = 8.0 * atan(1.0) ; // double d2r = pi2 / 360.0 ; // degrees to radians double ladv = 0.0 ; // secs launch advance double sid = SOL*YEAR/(YEAR+1.0) ; // secs sidereal day double rid = pi2 / sid ; // rad/sec loop angular velocity double vearth = pi2 * RE / sid ; // km/s equatorial velocity double vloop = vearth*cos(d2r*LOOPL ) ; // km/s loop rotation velocity double trad ; // sec/rad for computing double eps = 1.0/1024.0 ; // comparison "fuzz" double tmp ; double exp13 = 1.0/3.0 ; double eps1 = 1.0/65536.0 ; eps1 = eps1*eps1 ; // 2^(-32), small increment // ------------------------------------------------------------------------- // set up output files char fname[16]; // open the diagnostic output file, very large ... sprintf( fname, DIAG, NAME ); FILE * fd = fopen( fname , "w"); // open the summary output file, smaller ... sprintf( fname, SUM , NAME ); FILE * fs = fopen( fname , "w"); // open the wiki table output file sprintf( fname, WIKI, NAME ); FILE * fw = fopen( fname , "w"); // ------------------------------------------------------------------------- // print headers for diagnostic and output files fileheader( fd, vloop ); fileheader( fs, vloop ); // ------------------------------------------------------------------------- // inner loop parameters, declared here to facilitate end of loop printout double con_peri = RE + CON_PA ; // km construction perigee double pri_peri = RE + LOOPA ; // km prime perigee double con_time = sid * CON_SD ; // secs construction period // most values initialized to NAN to detect failed calculation double can_e = NAN ; // candidate eccentricity double E_con = NAN ; // construction eccentric anomaly double E_can = NAN ; // candidate eccentric anomaly double M_con = NAN ; // construction mean anomaly double M_can = NAN ; // candidate mean anomaly double delay_con = NAN ; // construction delay from apogee double delay_can = NAN ; // candidate delay from apogee double apogee_delay = NAN ; // FIXME describe double can_arrival = NAN ; // candidate arrival time double can_apo = NAN ; // candidate apogee radius double error_time = NAN ; // iterate in inner loop to zero double omega_con = NAN ; // construction angular velocity double omega_can = NAN ; // candidate angular velocity double v_cant = NAN ; // candidate double v_cont = NAN ; // construction double v_canr = NAN ; // candidate double v_conr = NAN ; // construction // four possible pairs of crossing radii, pairs should match // plus or minus angles (always should be plus for post-apogee capture?) double rcan = NAN ; // candidate crossing radius double rcon = NAN ; // construction crossing radius // ------------------------------------------------------------------------- // output header for wiki table wikiheader( fw, con_peri, pri_peri, CON_SD, con_time ); wikiheader( fs, con_peri, pri_peri, CON_SD, con_time ); wikiheader( fd, con_peri, pri_peri, CON_SD, con_time ); // ------------------------------------------------------------------------- // other construction orbit parameters double con_e ; // eccentricity (negative) double con_semi, con_apo ; // km semimajor apogee double con_v_semi, con_v_apo, con_v_peri ; // tangent velocities double con_v0 ; // characteristic velocity trad = con_time / pi2 ; // sec/rad con_semi = pow( MU0 *trad *trad , exp13 ) ; // km constr. semimajor axis con_apo = 2.0*con_semi - con_peri ; // km constr. apogee con_e = 1.0 - ( con_apo / con_semi ) ; // negative construction eccentr. con_v_semi = sqrt( MU0 / con_semi ) ; // km/s constr. semi. velocity con_v0 = sqrt( MU0 * con_semi / // km/s constr. char. velocity ( con_apo * con_peri ) ) ; con_v_apo = sqrt( con_peri / con_apo ) // km/s constr. apogee tangent * con_v_semi ; // velocity con_v_peri = sqrt( con_apo / con_peri ) // km/s constr. perigee tangent * con_v_semi ; // velocity // print information to wiki output file ------------------------- /// void printwiki( FILE * fo, char* p_name, /// double p_period, double p_launch, double p_arrival, /// double p_apogee, double p_perigee, /// double p_vt, double p_vr, double p_vp ) printwiki( fw, "construction ", con_time, NAN, 0.0, con_apo, con_peri, 0.0, 0.0, 0.0 ); printwiki( fd, "construction ", con_time, NAN, 0.0, con_apo, con_peri, 0.0, 0.0, 0.0 ); printwiki( fs, "construction ", con_time, NAN, 0.0, con_apo, con_peri, 0.0, 0.0, 0.0 ); // prime launch orbit parameters ------------------------- double pri_e ; // eccentricity (negative) double pri_time, pri_ltime ; // sec period, launch time double pri_semi, pri_apo ; // km semimajor apogee double pri_v_semi, pri_v_apo, pri_v_peri ; // km/s tangent velocities double pri_v0 ; // km/s characteristic velocity double pri_v_loop ; // km/s loop relative launch double pri_v_ins ; // km/s orbit insertion v. pri_apo = con_apo ; // km prime apogee pri_semi = 0.5*( pri_peri + pri_apo ) ; // km prime semimajor axis pri_e = 1.0 - ( pri_apo / pri_semi ) ; // prime eccentricity (negative) pri_time = pi2* sqrt(pow(pri_semi,3.0)/MU0 ); // sec prime period pri_v0 = sqrt( MU0 * pri_semi / // km/s prime char. velocity ( pri_apo * pri_peri ) ) ; pri_v_semi = sqrt( MU0 / pri_semi ) ; // km/s prime semi velocity pri_v_apo = sqrt( pri_peri / pri_apo ) // km/s prime apogee tangent * pri_v_semi ; // velocity pri_v_peri = sqrt( pri_apo / pri_peri ) // km/s prime perigee tangent * pri_v_semi ; // velocity pri_ltime = -0.5 * pri_time ; // sec minus launch time pri_v_loop = pri_v_peri - vloop ; // km/s prime launch velocity pri_v_ins = ( con_v_apo - pri_v_apo ); // km/s prime insertion /// void printwiki( FILE * fo, char* p_name, /// double p_period, double p_launch, double p_arrival, /// double p_apogee, double p_perigee, /// double p_vt, double p_vr, double p_vp ) printwiki( fw, " prime cargo ", pri_time, 0.0, 0.0, pri_apo, pri_peri, pri_v_ins, 0.0, 0.0 ); blankwiki( fw ); printwiki( fs, " prime cargo ", pri_time, 0.0, 0.0, pri_apo, pri_peri, pri_v_ins, 0.0, 0.0 ); printwiki( fd, " prime cargo ", pri_time, 0.0, 0.0, pri_apo, pri_peri, pri_v_ins, 0.0, 0.0 ); // print information to summary and diagnostic file ------------- /// printorbits( FILE( * fo, /// double con_time, /// double con_v0, double con_semi, /// double con_v_peri, double con_peri, /// double con_v_apo, double con_apo, /// double con_e, /// double pri_time, double pri_ltime; /// double pri_v0, double pri_semi, /// double pri_v_peri, double pri_peri, /// double pri_v_apo, double pri_apo, /// double pri_e, /// double pri_v_loop, double pri_v_ins ) ; printorbits( fs, con_time, con_v0, con_semi, con_v_peri, con_peri, con_v_apo, con_apo, con_e, pri_time, pri_ltime, pri_v0, pri_semi, pri_v_peri, pri_peri, pri_v_apo, pri_apo, pri_e, pri_v_loop, pri_v_ins ) ; printorbits( fd, con_time, con_v0, con_semi, con_v_peri, con_peri, con_v_apo, con_apo, con_e, pri_time, pri_ltime, pri_v0, pri_semi, pri_v_peri, pri_peri, pri_v_apo, pri_apo, pri_e, pri_v_loop, pri_v_ins ) ; // candidate launch orbit parameters --------------------- double can_time, can_ltime ; // sec period, launch time double can_semi, can_peri ; // km semimajor apogee perigee double can_v_semi, can_v_apo, can_v_peri ; // km/s tangent velocities double can_v0 ; // km/s characteristic velocity double can_v_loop ; // km/s loop relative launch double can_angle ; // rad perigee angle double can_v_ins ; // km/s orbit insert double theta_i ; // rad double tstart ; // from outer loop double thetac, degreec ; // angle for candidate orbit double theta, degree ; // angle for construction orbit double theta1, degree1 ; // angle possibility for crossing double theta2, degree2 ; // angle possibility for crossing // simplify calculations // double contop = con_semi * ( 1.0 - con_e * con_e ); //----------------------------------------------------------------------- //----------------------------------------------------------------------- // here is the main outer loop over tlstart //----------------------------------------------------------------------- for( tstart=TLMIN ; tstart <= TLMAX+eps ; tstart += TLSTEP ) {// L1 outer loop // test for near-zero, print prime orbit info and break if( fabs( tstart ) < eps ) { /// void printwiki( FILE * fo, char* p_name, /// double p_period, double p_launch, double p_arrival, /// double p_apogee, double p_perigee, /// double p_vt, double p_vr, double p_vp ) printwiki( fw, " prime cargo ", pri_time, 0.0, 0.0, pri_apo, pri_peri, pri_v_ins, 0.0, 0.0 ); printwiki( fs, " prime cargo ", pri_time, 0.0, 0.0, pri_apo, pri_peri, pri_v_ins, 0.0, 0.0 ); printwiki( fd, " prime cargo ", pri_time, 0.0, 0.0, pri_apo, pri_peri, pri_v_ins, 0.0, 0.0 ); } else { // if not the prime cargo orbit, tstart negative or positive // // iterative binary search for apogee with smallest arrival time error // estimated search bounds double can_a, can_a_min, can_a_max ; // given tstart, find the (clockwise, westward) apogee // angle for the candidate orbit. can_angle = rid * tstart ; // rad candidate launch angle // bounds for apogee guess (con_e is negative) can_a_min = contop / ( 1.0 + con_e * cos( can_angle ) ) ; can_a_max = con_semi*( 1.0 - 2.0*con_e ) ; // Kluge, HUGE fprintf( fd, "TEMP can_angle =%20.6e", can_angle ); fprintf( fd, " || can_a_min =%20.6e", can_a_min ); fprintf( fd, " || can_a_max =%20.6e\n", can_a_max ); //---------------------------------------------------------------------------- // inner loop to iterate can_a, start near center //---------------------------------------------------------------------------- error_time = 65536.0 ; // default if no solution found int numcan; // iterate NUMCAN times for( numcan = NUMCAN ; numcan-- > 0 ; ) { // L2 inner loop, can_a search can_a = 0.5*(can_a_min + can_a_max); // test midpoint // can_apo = pri_apo + can_a ; // WTF? can_apo = can_a ; can_ltime = pri_ltime + tstart ; // sec candidate launch time // if tstart is positive, launch is east of prime, can_angle positive; can_peri = pri_peri ; // km launch perigee can_semi = 0.5*( can_peri+can_apo ) ; // km candidate semimajor axis can_time = pi2 * sqrt( pow( can_semi,3.0 ) / MU0 ) ; // sec candidate period can_e = 1.0 - ( can_apo / can_semi ) ; // negative can. eccentricity can_v0 = sqrt( MU0 * can_semi / // km/s candidate char velocity ( can_apo * can_peri ) ) ; can_v_semi = sqrt( MU0 / can_semi ) ; // km/s candidate semi velocity can_v_apo = sqrt( can_peri / can_apo ) // km/s candidate apogee * can_v_semi ; // tangent velocity can_v_peri = sqrt( can_apo / can_peri ) // km/s candidate perigee * can_v_semi ; // tangent velocity can_v_loop = can_v_peri - vloop ; // km/s prime launch velocity // simplify the calculations: double cantop = can_semi * ( 1.0 - can_e * can_e ); // candidate launch orbit parameters --------------------- // find the crossing point of the construction and candidate launch orbits, // angle and radius, and the time that each object crosses them // // The apogee of the construction orbit is at theta == 0 // The apogee of the candidate launch orbit is at -can_angle, // for early launches, west of the construction station apogee // for later launches, east of the construction orbit apogee // // We will start by assuming the orbits are coplanar. They aren't; there // will be a small "northward" delta V where the planes intersect. // // The radius r at the crossing point of the candidate and construction // orbits is: // //1) r = cantop / ( 1 + can_e * cos(theta - can_angle) ) //2) r = contop / ( 1 + con_e * cos(theta) ) // // we know: con_semi, con_e, and can_angle // we've guessed: can_semi, can_e ... given the guess for: can_apo // // we want to find: theta // //-- First, combine the two equations and remove r (for now): --------------- // //3) cantop / ( 1 + can_e * cos(theta - can_angle) ) // = contop / ( 1 + con_e * cos(theta) ) // //--- Move the constants together: ------------------------------------------- // //4) cantop / contop // = ( 1 + can_e * cos(theta - can_angle) ) / ( 1 + con_e * cos(theta) ) // // Define K = cantop / contop double K = cantop / contop ; //5) K = ( 1 + can_e * cos(theta - can_angle) ) / ( 1 + con_e * cos(theta) ) //6) K * ( 1 + con_e * cos(theta) ) = 1 + can_e * cos(theta - can_angle) //7) K + K * con_e * cos(theta) = 1 + can_e * cos(theta - can_angle) // //-- Rearranging: ----------------------------------------------------------- // //8) can_e * cos(theta - can_angle) = K-1 + K * con_e * cos(theta) // //-- Take apart the compound cos with cos(x-y) = cos(x)cos(y)+sin(x)sin(y) --- // //9) can_e * cos(theta)*cos(can_angle) + can_e * sin(theta)*sin(can_angle) // = (K-1) + K * con_e*cos(theta) // //-- We are solving for theta, can_angle is a given "constant" // move constants to the front: // //10) can_e*cos(can_angle) * cos(theta) + can_e*sin(can_angle)* sin(theta) // =(K-1) + K*con_e* cos(theta) // //-- Collect the cos(theta) terms -------------------------------------------- // //11) ( can_e*cos(can_angle) - K*con_e ) * cos(theta) // = (K-1) - can_e*sin(can_angle)* sin(theta) // //-- Define another "constant" to simplify this equation: double D = can_e*cos(can_angle) - K*con_e ; //12) D * cos(theta) = (K-1) - can_e*sin(can_angle)* sin(theta) // //-- Divide by the constant factor in front: -------------------------------- // //13) cos(theta) = (K-1)/D - ( can_e*sin(can_angle)/D )*sin(theta) // //-- Define two more constants: double M = (K-1.0) / D ; double N = can_e * sin(can_angle) / D ; //-- the equation simplifies to: --------------------------------------------- // //14) cos(theta) = M - N*sin(theta) // //-- square the equation: ---------------------------------------------------- // //15) cos(theta)^2 = M^2 - 2*M*N* sin(theta) + N^2 * sin(theta)^2 // //-- replace cos^2 with 1-sin^2: --------------------------------------------- // //15) 1 - sin(theta)^2 = M^2 - 2*M*N* sin(theta) + N^2 * sin(theta)^2 // //-- Collect terms: ---------------------------------------------------------- // //16) 0 = (N^2 + 1)* sin(theta)^2 - 2*M*N* sin(theta) + (M^2 - 1) // //-- Define (but don't compute) three more constants: ------------------------ // //17) A = N^2 + 1 ; //18) B = -2*M*N ; //19) C = M^2 - 1 ; //-- We now have a simple polynomial: ---------------------------------------- // //20) 0 = A*sin(theta)^2 + B*sin(theta) + C // //-- And we can solve for that with the binomial equation // //21) sin(theta) = ( -B (+/-) sqrt[ B^2 - 4 A C ) ]/ 2 A // //-- "undefine" the constants: ----------------------------------------------- // //22) sin(theta) = ( 2*M*N (+/-) sqrt[ 4*M^2*N^2 - 4*((N^2 + 1)*(M^2 - 1) ] ) // / 2*( N^2 + 1 ) // //-- Factor out the 2 and sqrt(4): ------------------------------------------- // //23) sin(theta) = ( M*N (+/-) sqrt[ M^2*N^2 - (N^2 + 1)*(M^2 - 1) ] ) // / ( N^2 + 1 ) // //-- Replace the radical in the square root: --------------------------------- // //24) R = M^2*N^2 - (N^2 + 1)*(M^2 - 1) // //25) R = M^2*N^2 - N^2*M^2 + N^2 - M^2 + 1 = N^2 + 1 - M^2 // //-- and compute it: --------------------------------------------------------- double RAD = N*N + 1.0 - M*M ; //-- If RAD is negative, big trouble, no solution! // if( RAD < 0.0 ) { fprintf( fd, "\nFAIL\n\n"); fprintf( fd, "tstart=%5.1f can_semi=%9.3f\n", tstart, can_semi ); fprintf( fd, "N=%9.3e M=%9.3e" , N, M ); fprintf( fd, " ... RAD = %9.3e; imaginary roots\n", RAD ); fprintf( fd, "perhaps candidate apogee too low\n\nFAIL\n\n"); break; // fails here, something wrong with math, next loop } // 26) sin(theta) = ( M*N (+/-) sqrt( RAD ) ) / ( N^2 + 1 ) double sin1 = ( -M*N + sqrt( RAD ) ) / ( N*N + 1.0 ) ; double sin2 = ( -M*N - sqrt( RAD ) ) / ( N*N + 1.0 ) ; double theta1 = asin( sin1 ); double theta2 = asin( sin2 ); // DEBUG ============================================================ fprintf( fd, "\nDEB cantop=%11.3f contop=%11.3f K=%12.8f\n", cantop, contop, K ); fprintf( fd, "DEB D=%12.8f M=%12.8f N=%12.8f RAD=%12.8f\n", D, M, N, RAD ); fprintf( fd, "DEB sin1=%13.9f sin2=%13.9f\n", sin1, sin2 ); // DEBUG ============================================================ // RAD is positive! we can take the square root of it and compute two values if( tstart < 0.0 ) { theta = theta2 ; // choose this angle for early launches } else { theta = theta1 ; // choose this angle for late launches } // additional launch angle for candidate orbit // if the launch time is displaced west, then the candidate // orbit apogee is also displaced west, which means thetac must // be more positive to arrive at the same angle in the east thetac = theta + can_angle ; // compute rcon and rcan (should be the same!) rcon = contop / ( 1.0 + con_e * cos( theta ) ); rcan = cantop / ( 1.0 + can_e * cos( thetac ) ); // compute omega angular velocity for orbits omega_con = sqrt( MU0 / ( con_semi * con_semi * con_semi ) ); omega_can = sqrt( MU0 / ( can_semi * can_semi * can_semi ) ); // compute eccentric anomaly from apogee E_con = acos( (con_e + cos(theta )) / (1.0 + con_e * cos(theta )) ); E_can = acos( (can_e + cos(thetac)) / (1.0 + can_e * cos(thetac)) ); // compute mean anomaly from apogee M_con = E_con - con_e*sin( E_con ); M_can = E_can - can_e*sin( E_can ); // Compute delay from apogee to crossing delay_con = M_con / omega_con ; delay_can = M_can / omega_can ; // Compute half of difference between full orbit periods apogee_delay = 0.5*( can_time - pri_time ); // Compute candidate arrival and difference from construction arrival // // The construction orbit apogee is time=0, and the construction orbit crosses // The candidate orbit at time delay_con. // // The candidate orbit apogee time is apogee_delay + tstart // // The candidate orbit crosses the construction orbit later, // after an additional delay from apogee of delay_can can_arrival = apogee_delay + delay_can + tstart ; error_time = can_arrival - delay_con ; ///*** PRINT SUMMARY INFORMATION TO diagnostic file HERE /// void printmost{ FILE * fo, double tstart, /// double can_angle, double can_v_loop, /// double can_time, double can_ltime, double can_e, /// double can_semi, double can_peri, double can_apo, /// double can_v0, double can_v_peri, double can_v_apo, /// double rcon, double rcan, /// double theta1, double theta2, /// double theta, double thetac; /// double omega_con, double omega_can, /// double E_con, double E_can, /// double M_con, double M_can, /// double delay_con, double delay_can, /// double pri_time, double apogee_delay, double can_arrival, /// double v_canr, double v_cant, /// double v_conr, double v_cont, /// double error_time, double can_a_min, double can_a_max } printmost( fd, tstart, can_angle, can_v_loop, can_time, can_ltime, can_e, can_semi, can_peri, can_apo, can_v0, can_v_peri, can_v_apo, rcon, rcan, theta1, theta2, theta, thetac, omega_con, omega_can, E_con, E_can, M_con, M_can, delay_con, delay_can, pri_time, apogee_delay, can_arrival, v_canr, v_cant, v_conr, v_cont, error_time, can_a_min, can_a_max ); fprintf( fd, "\n-------------------------------------------" ); fprintf( fd, "-------------------------------------------\n" ); ///void printwiki( FILE * fo, char* p_name, /// double p_period, double p_launch, double p_arrival, /// double p_apogee, double p_perigee, /// double p_vt, double p_vr, double p_vp ) printwiki( fd, " inner loop ", can_time, tstart, can_arrival, can_apo, can_peri, NAN, NAN, NAN ); fprintf( fd, "--------------------------------------------" ); fprintf( fd, "------------------------------------------\n" ); // Loop bound recalculate. if error_time is negative, make the current // candidate apogee increment can_a into the new can_a_min. If the // error is positive, make can_a into the new can_a_max . fprintf( fd, "can_a=%13.8f OLD min=%13.8f max=%13.8f\n", can_a, can_a_min, can_a_max ); if( error_time < 0.0 ) can_a_min = can_a ; else can_a_max = can_a ; fprintf( fd, "can_a=%13.8f NEW min=%13.8f max=%13.8f\n", can_a, can_a_min, can_a_max ); } ///-------------------------------------------------------------- /// L2 end of inner loop /// ------------------------------------------------------------------------- fprintf( fd, "should be +0 or -0: error_time%+6.0f\n", error_time ); // End of iteration to estimate can_a, the apogee of the candidate orbit. // -------------------------- // Compute the tangential and radial velocities for // the capture and candidate orbits at angle theta // con_e, can_e are negative v_canr = can_e * can_v0 * sin( thetac ); v_cant = ( 1.0 - can_e*can_e ) * can_v0 / ( 1.0 - can_e * cos( thetac ) ); v_conr = con_e * con_v0 * sin( theta ); v_cont = ( 1.0 - con_e*con_e ) * con_v0 / ( 1.0 - con_e * cos( theta ) ); double delta_t = v_cont - v_cant ; double delta_r = v_conr - v_canr ; // Compute the plane crossing angle (half of the launch angle) // and then the "total" velocity vector at the candidate orbit plane crossing. // Then compute the N/S deltaV needed for the plane change. // // The plane change angle is the average of the candidate launch angle and the // prime launch angle (zero), hence half of the tstart angle, PRECEEDING // the prime launch angle. vcr, the total candidate orbit velocity at this // angle (vector sum of tangential and radial velocities) is multiplied by the // 2*sin of the change angle, times sin( latitude ). Note that both thetacr // and LOOPL are negative. // // modified - eccentricity is negative double thetacr = pi * tstart / sid ; // half of launch angle double vcr_r = can_e * can_v0 * sin( thetacr ); double vcr_t = ( 1.0 - can_e*can_e ) * can_v0 / ( 1.0 - can_e * cos( thetacr ) ) ; // Compute total plane crossing velocity // Not sure about sign double vcr = sqrt( vcr_r*vcr_r + vcr_t*vcr_t ); double delta_p = 2.0 * vcr * sin(thetacr) * sin( d2r * LOOPL ); // Print orbit parameters to FILE.wiki, wiki-formatted // first, build label for orbit line // char ssp[16] ; sprintf( ssp, "%5.0fs cargo ", tstart ); fprintf( fd, "\n==========================================="); fprintf( fd, "===========================================\n"); fprintf( fs, "\n==========================================="); fprintf( fs, "===========================================\n"); /// void printwiki( FILE * fo, char* p_name, /// double p_period, double p_launch, double p_arrival, /// double p_apogee, double p_perigee, /// double p_vt, double p_vr, double p_vp ) if( fabs(error_time) > 1e-6 ) { // if no convergence, zero out candidate information can_time = NAN ; can_arrival = NAN ; can_apo = NAN ; delta_t = NAN ; delta_r = NAN ; delta_p = NAN ; } printwiki( fw, ssp, can_time, tstart, can_arrival, can_apo, can_peri, delta_t, delta_r, delta_p ); printwiki( fs, ssp, can_time, tstart, can_arrival, can_apo, can_peri, delta_t, delta_r, delta_p ); printwiki( fd, ssp, can_time, tstart, can_arrival, can_apo, can_peri, delta_t, delta_r, delta_p ); printf( "%6.0f %20.12f\n", tstart, can_apo ); fprintf( fd, "============================================"); fprintf( fd, "==========================================\n"); fprintf( fs, "============================================"); fprintf( fs, "==========================================\n"); ///*** PRINT SUMMARY INFORMATION TO summary file HERE /// void printmost{ FILE * fo, double tstart, /// double can_angle, double can_v_loop, /// double can_time, double can_ltime, double can_e, /// double can_semi, double can_peri, double can_apo, /// double can_v0, double can_v_peri, double can_v_apo, /// double rcon, double rcan, /// double theta1, double theta2, /// double theta, double thetac, /// double omega_con, double omega_can, /// double E_con, double E_can, /// double M_con, double M_can, /// double delay_con, double delay_can, /// double pri_time, double apogee_delay, double can_arrival, /// double v_canr, double v_cant, /// double v_conr, double v_cont, /// double error_time, double can_a_min, double can_a_max } printmost( fs , tstart, can_angle, can_v_loop, can_time, can_ltime, can_e, can_semi, can_peri, can_apo, can_v0, can_v_peri, can_v_apo, rcon, rcan, theta1, theta2, theta, thetac, omega_con, omega_can, E_con, E_can, M_con, M_can, delay_con, delay_can, pri_time, apogee_delay, can_arrival, v_canr, v_cant, v_conr, v_cont, error_time, can_a_min, can_a_max ); } // L2 end of inner loop } // L1 end of outer loop // --------------------------------------------------- // end of outer loop, close files and return fclose( fd ); // diagnostic output file fclose( fs ); // summary output file fclose( fw ); // wiki output file return( 0 ); }