Attachment 'capture22.c'
Download 1 /* capture22.c KHL 2019 Feb 21
2 gcc -o capture22 capture22.c -lm ; ./capture22 > capture22.out
3 */
4 #define NAME "capture22"
5 /* =========================================================================
6 Times are relative to construction orbit apogee, hence negative eccentricity.
7 No correction for apsidal precession.
8
9 This starts with the [tstart] launch time, and uses an iterative loop to
10 find the best CAN_A extra apogee value, Rather than try to be too clever,
11 the search will start with a lower bound of CAN_A_min = [1 km/s] * [tstart]
12 and a starting upper bound of CAN_A_max = [2 km/s] * ([tstart] + 2 [s] ).
13
14 The midpoint of the two bounds is computed as the can_a value, which is
15 evaluated for the error time, the difference between the construction
16 and the candidate orbit times at the geometric orbit crossing
17
18 If the error time is negative, the lower bound is replaced with can_a .
19 If the error time is positive, the upper bound is replaced with can_a .
20
21 This is a fast calculation, so we will just iterate 60 times and let
22 the bounds converge. Due to numerical noise, the calculation will
23 be non-monotonic for small errors, but we will be close to the "actual"
24 value. In a real system, other perturbations (thruster error, tidal
25 effects from the equatorial bulge and the Moon, Sun, and Jupiter, and
26 general mechanical cussedness) will be much more important, and require
27 small midflight correction thrusts.
28
29 Below [tstart] of 1.0 seconds, the calculations become erratic;
30 search algorithms should lower-bounds-check at the projected
31 radius of the construction orbit at the [tstart] delay. The
32 curvature of the candidate orbit with that apogee will always
33 be less than the construction orbit, so there cannot be a
34 lower but intersecting candidate orbit. Add a small epsilon
35 just to make sure the calculation is bounded.
36
37 If an intersection is not found, OR the error time cannot be
38 made to converge to zero, skip that table entry and move on.
39 Set printme = 1 at the start of the loop, and set it to zero
40 if either the binomial root is imaginary or the loop does not
41 converge.
42
43 [tstart] CAN_A
44 0.003 3.0646 from 3.05 to 1e9
45 0.01 3.084670 from 0.19 to 1e9
46 0.03 3.14398 from 0.00 to 1e9
47 0.1 3.34698 from 0.00 to 1e9
48 0.3 3.899743 from 0.00 to 1e9
49 1.0 5.64861821 from 0.00 to 1e9
50 3.0 10.02098638 from 0.00 to 1e9
51 10.0 23.45589910 from 0.00 to 1e9
52 30.0 58.49035158 from 0.00 to 1e9
53 45.0 83.802205003 from 0.00 to 1e9
54 60.0 108.75943134887 from 0.00 to 1e9
55 90.0 158.08913293178 from 0.00 to 1e9
56 180.0 304.03736662938 from 0.00 to 1e9
57 300.0 496.931100566385 from 0.00 to 1e9
58 450.0 737.308802551446 from 0.00 to 1e9
59 600.0 977.727815882209 from 0.00 to 1e9
60 1200.0 1947.92465854122
61 1800.0 2946.23116529545
62 2400.0 3992.32125234230 from 0.00 to 1e9
63 3000.0 5108.68173943019
64 3600.0 6321.18761716776 from 0.00 to 1e9
65
66 ======================================================================*/
67
68 #define DIAG "%s.diag" // diagnostic output file
69 #define WIKI "%s.wiki" // main wiki output file
70 #define SUM "%s.sum" // main summary file
71
72 #define TLSTEP 100.0 // sec launch time step
73 #define TLMIN -900.0 // sec launch time start
74 #define TLMAX 900.0 // sec launch time range, minus to plus
75
76 #define LOOPA 80.0 // km loop altitude
77 #define LOOPL -8.0 // degrees south loop latitude
78 #define CON_PA 2000.0 // km construction orbit perigee altitude
79 #define CON_SD 1.0 // sidereal day construction orbit period
80
81 #define MU0 398600.4418 // km3/s2 Earth standard gravitational param
82 #define RE 6378.0 // km/s equatorial radius
83 #define SOL 86400.0 // sec solar day
84 #define YEAR 365.2422 // solar days per year
85
86 #define NUMCAN 60 // number of convergence iterations
87 #define EPS 1e-5 // crossing window
88
89 #include <stdio.h>
90 #include <math.h>
91
92 //=======================================================================
93 //=======================================================================
94 // wiki output subroutines -----------------------------------------
95 //=======================================================================
96
97
98 //=======================================================================
99 // print wiki header ------------
100 void wikiheader( FILE * fr, double p_constr, double p_prime,
101 double p_sday , double p_secs ) {
102
103 // first header line
104
105 fprintf( fr, "||<-7> %s ", NAME );
106 fprintf( fr, "construction perigee =%5.0f km", p_constr );
107 fprintf( fr, " launch perigee =%5.0f km", p_prime );
108 fprintf( fr, " ||\n" );
109
110 // second header line
111
112 fprintf( fr, "||<-7> " );
113 fprintf( fr, "sidereal period: %2.0f days,", p_sday );
114 fprintf( fr, "%10.3f seconds ", p_secs );
115 fprintf( fr, " ||\n" );
116
117 // third header line
118
119 fprintf( fr, "|| " ); // name of parameter
120 fprintf( fr, "|| period " ); // orbit period sec
121 fprintf( fr, "|| arrival " ); // arrival time sec
122 fprintf( fr, "|| apogee " ); // apogee km
123 fprintf( fr, "||<-3> velocity change m/s " );
124 fprintf( fr, "||\n" ); // end of first header line
125
126 // fourth header line
127
128 fprintf( fr, "|| " ); // name of parameter
129 fprintf( fr, "|| sec " ); // orbit period sec
130 fprintf( fr, "|| sec " ); // arrival time sec
131 fprintf( fr, "|| km " ); // apogee km
132 fprintf( fr, "|| tangent" );
133 fprintf( fr, "|| radial " );
134 fprintf( fr, "|| plane " );
135 fprintf( fr, "||\n" ); // end of second header line
136 } // end of wiki header
137
138 // print blank wiki line -------------
139
140 void blankwiki( FILE * fr ) {
141 fprintf( fr, "||<-7> ");
142 fprintf( fr, " ||\n");
143 }
144
145 // print line of wiki data ---------
146
147 void printwiki( FILE * fo, char* p_name,
148 double p_period, double p_launch, double p_arrival,
149 double p_apogee, double p_perigee,
150 double p_vt, double p_vr, double p_vp )
151 {
152 fprintf( fo, "|| %s" , p_name ); // name of parameter
153 fprintf( fo, "||%10.3f " , p_period ); // orbit period sec
154 fprintf( fo, "||%9.3f " , p_arrival ); // arrival time sec
155 fprintf( fo, "||%10.3f " , p_apogee ); // apogee km
156 fprintf( fo, "||%7.2f " , 1000.0*p_vt ); // tangential m/s
157 fprintf( fo, "||%8.2f " , 1000.0*p_vr ); // radial m/s
158 fprintf( fo, "||%6.2f " , 1000.0*p_vp ); // plane change m/s
159 fprintf( fo, "||\n" ); // end of table line
160 } // end of wiki print line
161
162 //=======================================================================
163 // print wiki, output or diagnostic header
164
165 void printhead( FILE * fo,
166 char* p_name,
167 double p_period,
168 double p_launch,
169 double p_arrival,
170 double p_apogee,
171 double p_perigee,
172 double p_vt,
173 double p_vr,
174 double p_vp
175 )
176 {
177 fprintf( fo, "|| %s" , p_name ); // name of parameter
178 fprintf( fo, "||%10.3f " , p_period ); // orbit period sec
179 fprintf( fo, "||%9.3f " , p_arrival ); // arrival time sec
180 fprintf( fo, "||%10.3f " , p_apogee ); // apogee km
181 fprintf( fo, "||%7.2f " , 1000.0*p_vt ); // tangential m/s
182 fprintf( fo, "||%8.2f " , 1000.0*p_vr ); // radial m/s
183 fprintf( fo, "||%6.2f " , 1000.0*p_vp ); // plane change m/s
184 fprintf( fo, "||\n" ); // end of table line
185
186 } // end of printhead; output or diagnostic header
187
188 //=======================================================================
189 // print construction and prime launch orbit, to output or diagnostic
190
191 void printorbits( FILE * fo,
192 double con_time,
193 double con_v0, double con_semi,
194 double con_v_peri, double con_peri,
195 double con_v_apo, double con_apo,
196 double con_e, // negative
197 double pri_time, double pri_ltime,
198 double pri_v0, double pri_semi,
199 double pri_v_peri, double pri_peri,
200 double pri_v_apo, double pri_apo,
201 double pri_e, // negative
202 double pri_v_loop, double pri_v_ins ) {
203
204 fprintf( fo, "\n construction orbit ----------" );
205 fprintf( fo, "--------------------------------------------\n" );
206 fprintf( fo, "%14.3f km semimajor axis", con_semi );
207 fprintf( fo, "%14.6f km/s V0 velocity\n", con_v0 );
208 fprintf( fo, "%14.3f km perigee", con_peri );
209 fprintf( fo, "%21.6f km/s perigee velocity\n", con_v_peri );
210 fprintf( fo, "%14.3f km apogee", con_apo );
211 fprintf( fo, "%22.6f km/s apogee velocity\n", con_v_apo );
212 fprintf( fo, "%14.3f sec orbit period", con_time );
213 fprintf( fo, "%16.6f eccentricity\n", con_e );
214 fprintf( fo, "\n prime launch orbit ----------" );
215 fprintf( fo, "--------------------------------------------\n" );
216 fprintf( fo, "%14.3f km semimajor axis", pri_semi );
217 fprintf( fo, "%14.6f km/s V0 velocity\n", pri_v0 );
218 fprintf( fo, "%14.3f km perigee", pri_peri );
219 fprintf( fo, "%21.6f km/s perigee velocity\n", pri_v_peri );
220 fprintf( fo, "%14.3f km apogee", pri_apo );
221 fprintf( fo, "%22.6f km/s apogee velocity\n", pri_v_apo );
222 fprintf( fo, "%14.3f sec orbit period", pri_time );
223 fprintf( fo, "%16.6f eccentricity\n", pri_e );
224 fprintf( fo, "%14.3f sec launch time\n", pri_ltime );
225 fprintf( fo, "%14.6f km/s insertion velocity", pri_v_ins );
226 fprintf( fo, "%10.5f km/s loop exit velocity\n", pri_v_loop );
227 }
228
229 //=======================================================================
230 // print most of data to diagnostic or summary file
231
232 /// void printmost{ FILE * fo, double tstart,
233 /// double can_angle, double can_v_loop,
234 /// double can_time, double can_ltime, double can_e,
235 /// double can_semi, double can_peri, double can_apo,
236 /// double can_v0, double can_v_peri, double can_v_apo,
237 /// double rcon, double rcan,
238 /// double theta1, double theta2,
239 // double theta, double thetac,
240 /// double omega_con, double omega_can,
241 /// double E_con, double E_can,
242 /// double M_con, double M_can,
243 /// double delay_con, double delay_can,
244 /// double pri_time, double apogee_delay, double can_arrival,
245 /// double v_canr, double v_cant,
246 /// double v_conr, double v_cont,
247 /// double error_time, double can_a_min, double can_a_max }
248
249 void printmost( FILE * fo, double tstart,
250 double can_angle, double can_v_loop,
251 double can_time, double can_ltime, double can_e,
252 double can_semi, double can_peri, double can_apo,
253 double can_v0, double can_v_peri, double can_v_apo,
254 double rcon, double rcan,
255 double theta1, double theta2,
256 double theta, double thetac,
257 double omega_con, double omega_can,
258 double E_con, double E_can,
259 double M_con, double M_can,
260 double delay_con, double delay_can,
261 double pri_time, double apogee_delay, double can_arrival,
262 double v_canr, double v_cant,
263 double v_conr, double v_cont,
264 double error_time, double can_a_min, double can_a_max )
265 {
266 double r2d = 45.0/atan(1.0) ;
267 double degree = r2d * theta ;
268 double degree1 = r2d * theta1 ;
269 double degree2 = r2d * theta2 ;
270 double degreec = r2d * thetac ;
271 double canad = r2d * can_angle ;
272
273 fprintf( fo, "\n candidate launch orbit -------------" );
274 fprintf( fo, "-------------------------------------\n" );
275 fprintf( fo, "%14.3f km semimajor axis", can_semi );
276 fprintf( fo, "%14.6f km/s V0 velocity\n", can_v0 );
277 fprintf( fo, "%14.3f km perigee", can_peri );
278 fprintf( fo, "%21.6f km/s perigee velocity\n", can_v_peri );
279 fprintf( fo, "%14.3f km apogee", can_apo );
280 fprintf( fo, "%22.6f km/s apogee velocity\n", can_v_apo );
281 fprintf( fo, "%14.3f sec orbit period", can_time );
282 fprintf( fo, "%16.6f eccentricity\n", can_e );
283 fprintf( fo, "%14.3f sec launch time", can_ltime );
284 fprintf( fo, "%+17.3f sec start time\n", tstart );
285 fprintf( fo, "%14.3f km/s loop velocity ", can_v_loop );
286
287 if( can_angle < 0.0 ) {
288 fprintf( fo, "%+13.6f rad,%7.4f° west of prime\n\n", can_angle, -canad );
289 } else {
290 fprintf( fo, "%+13.6f rad,%7.4f° east of prime\n\n", can_angle, canad );
291 }
292
293 fprintf( fo, "theta1=%12.7f theta2=%12.7f rad\n", theta1, theta2 );
294
295 fprintf( fo, " ");
296 fprintf( fo, " construction candidate\n");
297
298
299 fprintf( fo, "theta:%45.7f%16.7f rad\n", theta, thetac );
300
301 fprintf( fo, "radius:%44.5f%16.5f km\n", rcon, rcan );
302
303 fprintf( fo, "omega:%45.5e%16.5e rad/sec\n", omega_con, omega_can );
304
305 fprintf( fo, "eccentric anomaly E:%31.7f%16.7f rad\n", E_con, E_can );
306
307 fprintf( fo, "mean anomaly M:%36.7f%16.7f rad\n", M_con, M_can );
308
309 fprintf( fo, "time after apogee:" );
310 fprintf( fo, " %18.6f%16.6f sec\n", delay_con, delay_can );
311
312 fprintf( fo, "\n%+8.5f rad,%+9.4f° candidate orbit launch angle\n",
313 thetac, degreec );
314 fprintf( fo, "%12.6f sec construction delay from apogee to crossing\n",
315 delay_con );
316 fprintf( fo, "%12.6f sec candidate delay from apogee to crossing\n",
317 delay_can );
318 fprintf( fo, "%12.6f sec candidate orbit period\n", can_time );
319 fprintf( fo, "%12.6f sec prime orbit period\n", pri_time );
320 fprintf( fo, "%12.6f sec candidate apogee difference\n", apogee_delay );
321 fprintf( fo, "%12.6f sec candidate orbit arrival\n", can_arrival );
322 fprintf( fo,
323 "v_canr%9.5f v_cant%9.5f v_conr%9.5f v_cont%9.5f km/s\n",
324 v_canr, v_cant, v_conr, v_cont );
325 fprintf( fo,
326 "error_time=%15.4e sec can_apogee: min=%12.5f km max=%12.5f km\n",
327 error_time, can_a_min, can_a_max );
328 }
329
330 // =======================================================================
331 // print header for diagnostic and summary file
332
333 void fileheader( FILE * fo, double vloop ) {
334 fprintf(fo,
335 "============ /home/keithl/launchloop/capture/%s.c =============\n", NAME);
336 fprintf(fo,
337 "=================%2.0f Sidereal Day Construction Orbit =================\n",
338 CON_SD );
339 fprintf(fo, "\n Launch loop at %3.0f km altitude %3.0f°S latitude\n",
340 LOOPA, LOOPL );
341 fprintf(fo, "%14.6f km/s loop latitude rotation velocity\n", vloop );
342 }
343
344 //========================================================================
345 //========================================================================
346 // main program ---------------------------------------------------
347 //========================================================================
348 //========================================================================
349
350 int main(void) {
351
352 double pi = 4.0 * atan(1.0) ; //
353 double pi2 = 8.0 * atan(1.0) ; //
354 double d2r = pi2 / 360.0 ; // degrees to radians
355
356 double ladv = 0.0 ; // secs launch advance
357 double sid = SOL*YEAR/(YEAR+1.0) ; // secs sidereal day
358 double rid = pi2 / sid ; // rad/sec loop angular velocity
359 double vearth = pi2 * RE / sid ; // km/s equatorial velocity
360 double vloop = vearth*cos(d2r*LOOPL ) ; // km/s loop rotation velocity
361
362 double trad ; // sec/rad for computing
363
364 double eps = 1.0/1024.0 ; // comparison "fuzz"
365
366 double tmp ;
367 double exp13 = 1.0/3.0 ;
368
369 double eps1 = 1.0/65536.0 ;
370 eps1 = eps1*eps1 ; // 2^(-32), small increment
371
372 // -------------------------------------------------------------------------
373 // set up output files
374
375 char fname[16];
376
377 // open the diagnostic output file, very large ...
378 sprintf( fname, DIAG, NAME );
379 FILE * fd = fopen( fname , "w");
380
381 // open the summary output file, smaller ...
382 sprintf( fname, SUM , NAME );
383 FILE * fs = fopen( fname , "w");
384
385 // open the wiki table output file
386 sprintf( fname, WIKI, NAME );
387 FILE * fw = fopen( fname , "w");
388
389 // -------------------------------------------------------------------------
390 // print headers for diagnostic and output files
391
392 fileheader( fd, vloop );
393 fileheader( fs, vloop );
394
395 // -------------------------------------------------------------------------
396 // inner loop parameters, declared here to facilitate end of loop printout
397
398 double con_peri = RE + CON_PA ; // km construction perigee
399 double pri_peri = RE + LOOPA ; // km prime perigee
400 double con_time = sid * CON_SD ; // secs construction period
401
402 // most values initialized to NAN to detect failed calculation
403
404 double can_e = NAN ; // candidate eccentricity
405 double E_con = NAN ; // construction eccentric anomaly
406 double E_can = NAN ; // candidate eccentric anomaly
407 double M_con = NAN ; // construction mean anomaly
408 double M_can = NAN ; // candidate mean anomaly
409 double delay_con = NAN ; // construction delay from apogee
410 double delay_can = NAN ; // candidate delay from apogee
411 double apogee_delay = NAN ; // FIXME describe
412 double can_arrival = NAN ; // candidate arrival time
413 double can_apo = NAN ; // candidate apogee radius
414 double error_time = NAN ; // iterate in inner loop to zero
415 double omega_con = NAN ; // construction angular velocity
416 double omega_can = NAN ; // candidate angular velocity
417 double v_cant = NAN ; // candidate
418 double v_cont = NAN ; // construction
419 double v_canr = NAN ; // candidate
420 double v_conr = NAN ; // construction
421
422 // four possible pairs of crossing radii, pairs should match
423 // plus or minus angles (always should be plus for post-apogee capture?)
424
425 double rcan = NAN ; // candidate crossing radius
426 double rcon = NAN ; // construction crossing radius
427
428 // -------------------------------------------------------------------------
429 // output header for wiki table
430
431 wikiheader( fw, con_peri, pri_peri, CON_SD, con_time );
432 wikiheader( fs, con_peri, pri_peri, CON_SD, con_time );
433 wikiheader( fd, con_peri, pri_peri, CON_SD, con_time );
434
435 // -------------------------------------------------------------------------
436 // other construction orbit parameters
437
438 double con_e ; // eccentricity (negative)
439 double con_semi, con_apo ; // km semimajor apogee
440 double con_v_semi, con_v_apo, con_v_peri ; // tangent velocities
441 double con_v0 ; // characteristic velocity
442
443 trad = con_time / pi2 ; // sec/rad
444 con_semi = pow( MU0 *trad *trad , exp13 ) ; // km constr. semimajor axis
445 con_apo = 2.0*con_semi - con_peri ; // km constr. apogee
446 con_e = 1.0 - ( con_apo / con_semi ) ; // negative construction eccentr.
447 con_v_semi = sqrt( MU0 / con_semi ) ; // km/s constr. semi. velocity
448
449 con_v0 = sqrt( MU0 * con_semi / // km/s constr. char. velocity
450 ( con_apo * con_peri ) ) ;
451
452 con_v_apo = sqrt( con_peri / con_apo ) // km/s constr. apogee tangent
453 * con_v_semi ; // velocity
454
455 con_v_peri = sqrt( con_apo / con_peri ) // km/s constr. perigee tangent
456 * con_v_semi ; // velocity
457
458 // print information to wiki output file -------------------------
459
460 /// void printwiki( FILE * fo, char* p_name,
461 /// double p_period, double p_launch, double p_arrival,
462 /// double p_apogee, double p_perigee,
463 /// double p_vt, double p_vr, double p_vp )
464
465 printwiki( fw, "construction ",
466 con_time, NAN, 0.0,
467 con_apo, con_peri,
468 0.0, 0.0, 0.0 );
469
470 printwiki( fd, "construction ",
471 con_time, NAN, 0.0,
472 con_apo, con_peri,
473 0.0, 0.0, 0.0 );
474
475 printwiki( fs, "construction ",
476 con_time, NAN, 0.0,
477 con_apo, con_peri,
478 0.0, 0.0, 0.0 );
479
480
481 // prime launch orbit parameters -------------------------
482
483 double pri_e ; // eccentricity (negative)
484 double pri_time, pri_ltime ; // sec period, launch time
485 double pri_semi, pri_apo ; // km semimajor apogee
486 double pri_v_semi, pri_v_apo, pri_v_peri ; // km/s tangent velocities
487 double pri_v0 ; // km/s characteristic velocity
488 double pri_v_loop ; // km/s loop relative launch
489 double pri_v_ins ; // km/s orbit insertion v.
490
491 pri_apo = con_apo ; // km prime apogee
492
493 pri_semi = 0.5*( pri_peri + pri_apo ) ; // km prime semimajor axis
494
495 pri_e = 1.0 - ( pri_apo / pri_semi ) ; // prime eccentricity (negative)
496
497 pri_time = pi2* sqrt(pow(pri_semi,3.0)/MU0 ); // sec prime period
498
499 pri_v0 = sqrt( MU0 * pri_semi / // km/s prime char. velocity
500 ( pri_apo * pri_peri ) ) ;
501
502 pri_v_semi = sqrt( MU0 / pri_semi ) ; // km/s prime semi velocity
503
504 pri_v_apo = sqrt( pri_peri / pri_apo ) // km/s prime apogee tangent
505 * pri_v_semi ; // velocity
506
507 pri_v_peri = sqrt( pri_apo / pri_peri ) // km/s prime perigee tangent
508 * pri_v_semi ; // velocity
509
510 pri_ltime = -0.5 * pri_time ; // sec minus launch time
511
512 pri_v_loop = pri_v_peri - vloop ; // km/s prime launch velocity
513
514 pri_v_ins = ( con_v_apo - pri_v_apo ); // km/s prime insertion
515
516 /// void printwiki( FILE * fo, char* p_name,
517 /// double p_period, double p_launch, double p_arrival,
518 /// double p_apogee, double p_perigee,
519 /// double p_vt, double p_vr, double p_vp )
520
521 printwiki( fw, " prime cargo ",
522 pri_time, 0.0, 0.0,
523 pri_apo, pri_peri,
524 pri_v_ins, 0.0, 0.0 );
525
526 blankwiki( fw );
527
528 printwiki( fs, " prime cargo ",
529 pri_time, 0.0, 0.0,
530 pri_apo, pri_peri,
531 pri_v_ins, 0.0, 0.0 );
532
533 printwiki( fd, " prime cargo ",
534 pri_time, 0.0, 0.0,
535 pri_apo, pri_peri,
536 pri_v_ins, 0.0, 0.0 );
537
538 // print information to summary and diagnostic file -------------
539
540 /// printorbits( FILE( * fo,
541 /// double con_time,
542 /// double con_v0, double con_semi,
543 /// double con_v_peri, double con_peri,
544 /// double con_v_apo, double con_apo,
545 /// double con_e,
546 /// double pri_time, double pri_ltime;
547 /// double pri_v0, double pri_semi,
548 /// double pri_v_peri, double pri_peri,
549 /// double pri_v_apo, double pri_apo,
550 /// double pri_e,
551 /// double pri_v_loop, double pri_v_ins ) ;
552
553 printorbits( fs, con_time,
554 con_v0, con_semi,
555 con_v_peri, con_peri,
556 con_v_apo, con_apo,
557 con_e,
558 pri_time, pri_ltime,
559 pri_v0, pri_semi,
560 pri_v_peri, pri_peri,
561 pri_v_apo, pri_apo,
562 pri_e,
563 pri_v_loop, pri_v_ins ) ;
564
565 printorbits( fd, con_time,
566 con_v0, con_semi,
567 con_v_peri, con_peri,
568 con_v_apo, con_apo,
569 con_e,
570 pri_time, pri_ltime,
571 pri_v0, pri_semi,
572 pri_v_peri, pri_peri,
573 pri_v_apo, pri_apo,
574 pri_e,
575 pri_v_loop, pri_v_ins ) ;
576
577 // candidate launch orbit parameters ---------------------
578
579 double can_time, can_ltime ; // sec period, launch time
580 double can_semi, can_peri ; // km semimajor apogee perigee
581 double can_v_semi, can_v_apo, can_v_peri ; // km/s tangent velocities
582 double can_v0 ; // km/s characteristic velocity
583 double can_v_loop ; // km/s loop relative launch
584 double can_angle ; // rad perigee angle
585 double can_v_ins ; // km/s orbit insert
586 double theta_i ; // rad
587 double tstart ; // from outer loop
588
589 double thetac, degreec ; // angle for candidate orbit
590 double theta, degree ; // angle for construction orbit
591 double theta1, degree1 ; // angle possibility for crossing
592 double theta2, degree2 ; // angle possibility for crossing
593
594 // simplify calculations
595 //
596 double contop = con_semi * ( 1.0 - con_e * con_e );
597
598 //-----------------------------------------------------------------------
599 //-----------------------------------------------------------------------
600 // here is the main outer loop over tlstart
601 //-----------------------------------------------------------------------
602
603 for( tstart=TLMIN ; tstart <= TLMAX+eps ; tstart += TLSTEP ) {// L1 outer loop
604
605 // test for near-zero, print prime orbit info and break
606
607 if( fabs( tstart ) < eps ) {
608
609 /// void printwiki( FILE * fo, char* p_name,
610 /// double p_period, double p_launch, double p_arrival,
611 /// double p_apogee, double p_perigee,
612 /// double p_vt, double p_vr, double p_vp )
613
614 printwiki( fw, " prime cargo ",
615 pri_time, 0.0, 0.0,
616 pri_apo, pri_peri,
617 pri_v_ins, 0.0, 0.0 );
618
619 printwiki( fs, " prime cargo ",
620 pri_time, 0.0, 0.0,
621 pri_apo, pri_peri,
622 pri_v_ins, 0.0, 0.0 );
623
624 printwiki( fd, " prime cargo ",
625 pri_time, 0.0, 0.0,
626 pri_apo, pri_peri,
627 pri_v_ins, 0.0, 0.0 );
628 } else
629 {
630
631 // if not the prime cargo orbit, tstart negative or positive
632 //
633 // iterative binary search for apogee with smallest arrival time error
634 // estimated search bounds
635
636 double can_a, can_a_min, can_a_max ;
637
638 // given tstart, find the (clockwise, westward) apogee
639 // angle for the candidate orbit.
640
641 can_angle = rid * tstart ; // rad candidate launch angle
642
643 // bounds for apogee guess (con_e is negative)
644
645 can_a_min = contop / ( 1.0 + con_e * cos( can_angle ) ) ;
646 can_a_max = con_semi*( 1.0 - 2.0*con_e ) ; // Kluge, HUGE
647
648 fprintf( fd, "TEMP can_angle =%20.6e", can_angle );
649 fprintf( fd, " || can_a_min =%20.6e", can_a_min );
650 fprintf( fd, " || can_a_max =%20.6e\n", can_a_max );
651
652 //----------------------------------------------------------------------------
653 // inner loop to iterate can_a, start near center
654 //----------------------------------------------------------------------------
655
656 error_time = 65536.0 ; // default if no solution found
657 int numcan; // iterate NUMCAN times
658 for( numcan = NUMCAN ; numcan-- > 0 ; ) { // L2 inner loop, can_a search
659
660 can_a = 0.5*(can_a_min + can_a_max); // test midpoint
661 // can_apo = pri_apo + can_a ; // WTF?
662 can_apo = can_a ;
663 can_ltime = pri_ltime + tstart ; // sec candidate launch time
664
665 // if tstart is positive, launch is east of prime, can_angle positive;
666
667 can_peri = pri_peri ; // km launch perigee
668 can_semi = 0.5*( can_peri+can_apo ) ; // km candidate semimajor axis
669
670 can_time = pi2 *
671 sqrt( pow( can_semi,3.0 ) / MU0 ) ; // sec candidate period
672
673 can_e = 1.0 - ( can_apo / can_semi ) ; // negative can. eccentricity
674
675 can_v0 = sqrt( MU0 * can_semi / // km/s candidate char velocity
676 ( can_apo * can_peri ) ) ;
677
678 can_v_semi = sqrt( MU0 / can_semi ) ; // km/s candidate semi velocity
679
680 can_v_apo = sqrt( can_peri / can_apo ) // km/s candidate apogee
681 * can_v_semi ; // tangent velocity
682
683 can_v_peri = sqrt( can_apo / can_peri ) // km/s candidate perigee
684 * can_v_semi ; // tangent velocity
685
686 can_v_loop = can_v_peri - vloop ; // km/s prime launch velocity
687
688 // simplify the calculations:
689
690 double cantop = can_semi * ( 1.0 - can_e * can_e );
691
692 // candidate launch orbit parameters ---------------------
693
694 // find the crossing point of the construction and candidate launch orbits,
695 // angle and radius, and the time that each object crosses them
696 //
697 // The apogee of the construction orbit is at theta == 0
698 // The apogee of the candidate launch orbit is at -can_angle,
699 // for early launches, west of the construction station apogee
700 // for later launches, east of the construction orbit apogee
701 //
702 // We will start by assuming the orbits are coplanar. They aren't; there
703 // will be a small "northward" delta V where the planes intersect.
704 //
705 // The radius r at the crossing point of the candidate and construction
706 // orbits is:
707 //
708 //1) r = cantop / ( 1 + can_e * cos(theta - can_angle) )
709 //2) r = contop / ( 1 + con_e * cos(theta) )
710 //
711 // we know: con_semi, con_e, and can_angle
712 // we've guessed: can_semi, can_e ... given the guess for: can_apo
713 //
714 // we want to find: theta
715 //
716 //-- First, combine the two equations and remove r (for now): ---------------
717 //
718 //3) cantop / ( 1 + can_e * cos(theta - can_angle) )
719 // = contop / ( 1 + con_e * cos(theta) )
720 //
721 //--- Move the constants together: -------------------------------------------
722 //
723 //4) cantop / contop
724 // = ( 1 + can_e * cos(theta - can_angle) ) / ( 1 + con_e * cos(theta) )
725 //
726 // Define K = cantop / contop
727
728 double K = cantop / contop ;
729
730 //5) K = ( 1 + can_e * cos(theta - can_angle) ) / ( 1 + con_e * cos(theta) )
731 //6) K * ( 1 + con_e * cos(theta) ) = 1 + can_e * cos(theta - can_angle)
732 //7) K + K * con_e * cos(theta) = 1 + can_e * cos(theta - can_angle)
733 //
734 //-- Rearranging: -----------------------------------------------------------
735 //
736 //8) can_e * cos(theta - can_angle) = K-1 + K * con_e * cos(theta)
737 //
738 //-- Take apart the compound cos with cos(x-y) = cos(x)cos(y)+sin(x)sin(y) ---
739 //
740 //9) can_e * cos(theta)*cos(can_angle) + can_e * sin(theta)*sin(can_angle)
741 // = (K-1) + K * con_e*cos(theta)
742 //
743 //-- We are solving for theta, can_angle is a given "constant"
744 // move constants to the front:
745 //
746 //10) can_e*cos(can_angle) * cos(theta) + can_e*sin(can_angle)* sin(theta)
747 // =(K-1) + K*con_e* cos(theta)
748 //
749 //-- Collect the cos(theta) terms --------------------------------------------
750 //
751 //11) ( can_e*cos(can_angle) - K*con_e ) * cos(theta)
752 // = (K-1) - can_e*sin(can_angle)* sin(theta)
753 //
754 //-- Define another "constant" to simplify this equation:
755
756 double D = can_e*cos(can_angle) - K*con_e ;
757
758 //12) D * cos(theta) = (K-1) - can_e*sin(can_angle)* sin(theta)
759 //
760 //-- Divide by the constant factor in front: --------------------------------
761 //
762 //13) cos(theta) = (K-1)/D - ( can_e*sin(can_angle)/D )*sin(theta)
763 //
764 //-- Define two more constants:
765
766 double M = (K-1.0) / D ;
767 double N = can_e * sin(can_angle) / D ;
768
769 //-- the equation simplifies to: ---------------------------------------------
770 //
771 //14) cos(theta) = M - N*sin(theta)
772 //
773 //-- square the equation: ----------------------------------------------------
774 //
775 //15) cos(theta)^2 = M^2 - 2*M*N* sin(theta) + N^2 * sin(theta)^2
776 //
777 //-- replace cos^2 with 1-sin^2: ---------------------------------------------
778 //
779 //15) 1 - sin(theta)^2 = M^2 - 2*M*N* sin(theta) + N^2 * sin(theta)^2
780 //
781 //-- Collect terms: ----------------------------------------------------------
782 //
783 //16) 0 = (N^2 + 1)* sin(theta)^2 - 2*M*N* sin(theta) + (M^2 - 1)
784 //
785 //-- Define (but don't compute) three more constants: ------------------------
786 //
787 //17) A = N^2 + 1 ;
788 //18) B = -2*M*N ;
789 //19) C = M^2 - 1 ;
790
791 //-- We now have a simple polynomial: ----------------------------------------
792 //
793 //20) 0 = A*sin(theta)^2 + B*sin(theta) + C
794 //
795 //-- And we can solve for that with the binomial equation
796 //
797 //21) sin(theta) = ( -B (+/-) sqrt[ B^2 - 4 A C ) ]/ 2 A
798 //
799 //-- "undefine" the constants: -----------------------------------------------
800 //
801 //22) sin(theta) = ( 2*M*N (+/-) sqrt[ 4*M^2*N^2 - 4*((N^2 + 1)*(M^2 - 1) ] )
802 // / 2*( N^2 + 1 )
803 //
804 //-- Factor out the 2 and sqrt(4): -------------------------------------------
805 //
806 //23) sin(theta) = ( M*N (+/-) sqrt[ M^2*N^2 - (N^2 + 1)*(M^2 - 1) ] )
807 // / ( N^2 + 1 )
808 //
809 //-- Replace the radical in the square root: ---------------------------------
810 //
811 //24) R = M^2*N^2 - (N^2 + 1)*(M^2 - 1)
812 //
813 //25) R = M^2*N^2 - N^2*M^2 + N^2 - M^2 + 1 = N^2 + 1 - M^2
814 //
815 //-- and compute it: ---------------------------------------------------------
816
817 double RAD = N*N + 1.0 - M*M ;
818
819 //-- If RAD is negative, big trouble, no solution!
820 //
821 if( RAD < 0.0 ) {
822 fprintf( fd, "\nFAIL\n\n");
823 fprintf( fd, "tstart=%5.1f can_semi=%9.3f\n", tstart, can_semi );
824 fprintf( fd, "N=%9.3e M=%9.3e" , N, M );
825 fprintf( fd, " ... RAD = %9.3e; imaginary roots\n", RAD );
826 fprintf( fd, "perhaps candidate apogee too low\n\nFAIL\n\n");
827 break; // fails here, something wrong with math, next loop
828 }
829
830 // 26) sin(theta) = ( M*N (+/-) sqrt( RAD ) ) / ( N^2 + 1 )
831
832 double sin1 = ( -M*N + sqrt( RAD ) ) / ( N*N + 1.0 ) ;
833 double sin2 = ( -M*N - sqrt( RAD ) ) / ( N*N + 1.0 ) ;
834
835 double theta1 = asin( sin1 );
836 double theta2 = asin( sin2 );
837
838 // DEBUG ============================================================
839
840 fprintf( fd, "\nDEB cantop=%11.3f contop=%11.3f K=%12.8f\n",
841 cantop, contop, K );
842 fprintf( fd, "DEB D=%12.8f M=%12.8f N=%12.8f RAD=%12.8f\n",
843 D, M, N, RAD );
844 fprintf( fd, "DEB sin1=%13.9f sin2=%13.9f\n",
845 sin1, sin2 );
846
847 // DEBUG ============================================================
848
849
850 // RAD is positive! we can take the square root of it and compute two values
851
852 if( tstart < 0.0 ) {
853 theta = theta2 ; // choose this angle for early launches
854 } else {
855 theta = theta1 ; // choose this angle for late launches
856 }
857
858 // additional launch angle for candidate orbit
859 // if the launch time is displaced west, then the candidate
860 // orbit apogee is also displaced west, which means thetac must
861 // be more positive to arrive at the same angle in the east
862
863 thetac = theta + can_angle ;
864
865 // compute rcon and rcan (should be the same!)
866
867 rcon = contop / ( 1.0 + con_e * cos( theta ) );
868 rcan = cantop / ( 1.0 + can_e * cos( thetac ) );
869
870 // compute omega angular velocity for orbits
871
872 omega_con = sqrt( MU0 / ( con_semi * con_semi * con_semi ) );
873 omega_can = sqrt( MU0 / ( can_semi * can_semi * can_semi ) );
874
875 // compute eccentric anomaly from apogee
876
877 E_con = acos( (con_e + cos(theta )) / (1.0 + con_e * cos(theta )) );
878 E_can = acos( (can_e + cos(thetac)) / (1.0 + can_e * cos(thetac)) );
879
880 // compute mean anomaly from apogee
881
882 M_con = E_con - con_e*sin( E_con );
883 M_can = E_can - can_e*sin( E_can );
884
885 // Compute delay from apogee to crossing
886
887 delay_con = M_con / omega_con ;
888 delay_can = M_can / omega_can ;
889
890 // Compute half of difference between full orbit periods
891 apogee_delay = 0.5*( can_time - pri_time );
892
893 // Compute candidate arrival and difference from construction arrival
894 //
895 // The construction orbit apogee is time=0, and the construction orbit crosses
896 // The candidate orbit at time delay_con.
897 //
898 // The candidate orbit apogee time is apogee_delay + tstart
899 //
900 // The candidate orbit crosses the construction orbit later,
901 // after an additional delay from apogee of delay_can
902
903 can_arrival = apogee_delay + delay_can + tstart ;
904 error_time = can_arrival - delay_con ;
905
906 ///*** PRINT SUMMARY INFORMATION TO diagnostic file HERE
907
908 /// void printmost{ FILE * fo, double tstart,
909 /// double can_angle, double can_v_loop,
910 /// double can_time, double can_ltime, double can_e,
911 /// double can_semi, double can_peri, double can_apo,
912 /// double can_v0, double can_v_peri, double can_v_apo,
913 /// double rcon, double rcan,
914 /// double theta1, double theta2,
915 /// double theta, double thetac;
916 /// double omega_con, double omega_can,
917 /// double E_con, double E_can,
918 /// double M_con, double M_can,
919 /// double delay_con, double delay_can,
920 /// double pri_time, double apogee_delay, double can_arrival,
921 /// double v_canr, double v_cant,
922 /// double v_conr, double v_cont,
923 /// double error_time, double can_a_min, double can_a_max }
924
925 printmost( fd, tstart,
926 can_angle, can_v_loop,
927 can_time, can_ltime, can_e,
928 can_semi, can_peri, can_apo,
929 can_v0, can_v_peri, can_v_apo,
930 rcon, rcan,
931 theta1, theta2,
932 theta, thetac,
933 omega_con, omega_can,
934 E_con, E_can,
935 M_con, M_can,
936 delay_con, delay_can,
937 pri_time, apogee_delay, can_arrival,
938 v_canr, v_cant,
939 v_conr, v_cont,
940 error_time, can_a_min, can_a_max );
941
942 fprintf( fd, "\n-------------------------------------------" );
943 fprintf( fd, "-------------------------------------------\n" );
944
945 ///void printwiki( FILE * fo, char* p_name,
946 /// double p_period, double p_launch, double p_arrival,
947 /// double p_apogee, double p_perigee,
948 /// double p_vt, double p_vr, double p_vp )
949
950 printwiki( fd, " inner loop ",
951 can_time, tstart, can_arrival,
952 can_apo, can_peri,
953 NAN, NAN, NAN );
954
955 fprintf( fd, "--------------------------------------------" );
956 fprintf( fd, "------------------------------------------\n" );
957
958 // Loop bound recalculate. if error_time is negative, make the current
959 // candidate apogee increment can_a into the new can_a_min. If the
960 // error is positive, make can_a into the new can_a_max .
961
962 fprintf( fd, "can_a=%13.8f OLD min=%13.8f max=%13.8f\n",
963 can_a, can_a_min, can_a_max );
964
965 if( error_time < 0.0 ) can_a_min = can_a ;
966 else can_a_max = can_a ;
967
968 fprintf( fd, "can_a=%13.8f NEW min=%13.8f max=%13.8f\n",
969 can_a, can_a_min, can_a_max );
970
971 } ///--------------------------------------------------------------
972 /// L2 end of inner loop
973 /// -------------------------------------------------------------------------
974
975 fprintf( fd, "should be +0 or -0: error_time%+6.0f\n", error_time );
976
977
978 // End of iteration to estimate can_a, the apogee of the candidate orbit.
979
980 // --------------------------
981 // Compute the tangential and radial velocities for
982 // the capture and candidate orbits at angle theta
983 // con_e, can_e are negative
984
985 v_canr = can_e * can_v0 * sin( thetac );
986 v_cant = ( 1.0 - can_e*can_e ) * can_v0 / ( 1.0 - can_e * cos( thetac ) );
987
988 v_conr = con_e * con_v0 * sin( theta );
989 v_cont = ( 1.0 - con_e*con_e ) * con_v0 / ( 1.0 - con_e * cos( theta ) );
990
991 double delta_t = v_cont - v_cant ;
992 double delta_r = v_conr - v_canr ;
993
994 // Compute the plane crossing angle (half of the launch angle)
995 // and then the "total" velocity vector at the candidate orbit plane crossing.
996 // Then compute the N/S deltaV needed for the plane change.
997 //
998 // The plane change angle is the average of the candidate launch angle and the
999 // prime launch angle (zero), hence half of the tstart angle, PRECEEDING
1000 // the prime launch angle. vcr, the total candidate orbit velocity at this
1001 // angle (vector sum of tangential and radial velocities) is multiplied by the
1002 // 2*sin of the change angle, times sin( latitude ). Note that both thetacr
1003 // and LOOPL are negative.
1004 //
1005 // modified - eccentricity is negative
1006
1007 double thetacr = pi * tstart / sid ; // half of launch angle
1008
1009 double vcr_r = can_e * can_v0 * sin( thetacr );
1010
1011 double vcr_t = ( 1.0 - can_e*can_e ) * can_v0
1012 / ( 1.0 - can_e * cos( thetacr ) ) ;
1013
1014 // Compute total plane crossing velocity
1015 // Not sure about sign
1016
1017 double vcr = sqrt( vcr_r*vcr_r + vcr_t*vcr_t );
1018
1019 double delta_p = 2.0 * vcr * sin(thetacr) * sin( d2r * LOOPL );
1020
1021 // Print orbit parameters to FILE.wiki, wiki-formatted
1022 // first, build label for orbit line
1023 //
1024 char ssp[16] ;
1025 sprintf( ssp, "%5.0fs cargo ", tstart );
1026
1027 fprintf( fd, "\n===========================================");
1028 fprintf( fd, "===========================================\n");
1029 fprintf( fs, "\n===========================================");
1030 fprintf( fs, "===========================================\n");
1031
1032 /// void printwiki( FILE * fo, char* p_name,
1033 /// double p_period, double p_launch, double p_arrival,
1034 /// double p_apogee, double p_perigee,
1035 /// double p_vt, double p_vr, double p_vp )
1036
1037
1038 if( fabs(error_time) > 1e-6 ) {
1039 // if no convergence, zero out candidate information
1040 can_time = NAN ;
1041 can_arrival = NAN ;
1042 can_apo = NAN ;
1043 delta_t = NAN ;
1044 delta_r = NAN ;
1045 delta_p = NAN ;
1046 }
1047
1048 printwiki( fw, ssp,
1049 can_time, tstart, can_arrival,
1050 can_apo, can_peri,
1051 delta_t, delta_r, delta_p );
1052
1053 printwiki( fs, ssp,
1054 can_time, tstart, can_arrival,
1055 can_apo, can_peri,
1056 delta_t, delta_r, delta_p );
1057
1058 printwiki( fd, ssp,
1059 can_time, tstart, can_arrival,
1060 can_apo, can_peri,
1061 delta_t, delta_r, delta_p );
1062
1063 printf( "%6.0f %20.12f\n", tstart, can_apo );
1064
1065 fprintf( fd, "============================================");
1066 fprintf( fd, "==========================================\n");
1067 fprintf( fs, "============================================");
1068 fprintf( fs, "==========================================\n");
1069
1070 ///*** PRINT SUMMARY INFORMATION TO summary file HERE
1071
1072 /// void printmost{ FILE * fo, double tstart,
1073 /// double can_angle, double can_v_loop,
1074 /// double can_time, double can_ltime, double can_e,
1075 /// double can_semi, double can_peri, double can_apo,
1076 /// double can_v0, double can_v_peri, double can_v_apo,
1077 /// double rcon, double rcan,
1078 /// double theta1, double theta2,
1079 /// double theta, double thetac,
1080 /// double omega_con, double omega_can,
1081 /// double E_con, double E_can,
1082 /// double M_con, double M_can,
1083 /// double delay_con, double delay_can,
1084 /// double pri_time, double apogee_delay, double can_arrival,
1085 /// double v_canr, double v_cant,
1086 /// double v_conr, double v_cont,
1087 /// double error_time, double can_a_min, double can_a_max }
1088
1089 printmost( fs , tstart,
1090 can_angle, can_v_loop,
1091 can_time, can_ltime, can_e,
1092 can_semi, can_peri, can_apo,
1093 can_v0, can_v_peri, can_v_apo,
1094 rcon, rcan,
1095 theta1, theta2,
1096 theta, thetac,
1097 omega_con, omega_can,
1098 E_con, E_can,
1099 M_con, M_can,
1100 delay_con, delay_can,
1101 pri_time, apogee_delay, can_arrival,
1102 v_canr, v_cant,
1103 v_conr, v_cont,
1104 error_time, can_a_min, can_a_max );
1105
1106 } // L2 end of inner loop
1107 } // L1 end of outer loop
1108
1109 // ---------------------------------------------------
1110 // end of outer loop, close files and return
1111
1112 fclose( fd ); // diagnostic output file
1113 fclose( fs ); // summary output file
1114 fclose( fw ); // wiki output file
1115 return( 0 );
1116 }
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.