source: Sophya/trunk/SophyaExt/XephemAstroLib/deep.c@ 3632

Last change on this file since 3632 was 3477, checked in by cmv, 18 years ago

mise a jour Xephem 3.7.3 , cmv 25/03/2008
File size: 21.1 KB
RevLine 
[1457]1#include <stdlib.h>
2#include <math.h>
3
4#include "deepconst.h"
5#include "satspec.h"
6
7/* * DEEP SPACE 31 OCT 80 */
8/* SUBROUTINE DEEP */
9/* COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O, */
10/* 1 XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50 */
11/* COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD, */
12/* 1 XJ3,XKE,XKMPER,XMNPDA,AE */
13/* COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2 */
14/* DOUBLE PRECISION EPOCH, DS50 */
15/* DOUBLE PRECISION */
16/* * DAY,PREEP,XNODCE,ATIME,DELT,SAVTSN,STEP2,STEPN,STEPP */
17/* DATA ZNS, C1SS, ZES/ */
18/* A 1.19459E-5, 2.9864797E-6, .01675/ */
19/* DATA ZNL, C1L, ZEL/ */
20/* A 1.5835218E-4, 4.7968065E-7, .05490/ */
21/* DATA ZCOSIS, ZSINIS, ZSINGS/ */
22/* A .91744867, .39785416, -.98088458/ */
23/* DATA ZCOSGS, ZCOSHS, ZSINHS/ */
24/* A .1945905, 1.0, 0.0/ */
25/* DATA Q22,Q31,Q33/1.7891679E-6,2.1460748E-6,2.2123015E-7/ */
26/* DATA G22,G32/5.7686396,0.95240898/ */
27/* DATA G44,G52/1.8014998,1.0508330/ */
28/* DATA G54/4.4108898/ */
29/* DATA ROOT22,ROOT32/1.7891679E-6,3.7393792E-7/ */
30/* DATA ROOT44,ROOT52/7.3636953E-9,1.1428639E-7/ */
31/* DATA ROOT54/2.1765803E-9/ */
32/* DATA THDT/4.3752691E-3/ */
33
34#define XMO (sat->elem->se_XMO)
35#define XNODEO (sat->elem->se_XNODEO)
36#define OMEGAO (sat->elem->se_OMEGAO)
37#define EO (sat->elem->se_EO)
38#define XINCL (sat->elem->se_XINCL)
39#define XNO (sat->elem->se_XNO)
40#define XNDT20 (sat->elem->se_XNDT20)
41#define XNDD60 (sat->elem->se_XNDD60)
42#define BSTAR (sat->elem->se_BSTAR)
43#define EPOCH (sat->elem->se_EPOCH)
44
45#define ZNS (1.19459E-5)
46#define C1SS (2.9864797E-6)
47#define ZES (.01675)
48#define ZNL (1.5835218E-4)
49#define C1L (4.7968065E-7)
50#define ZEL (.05490)
51#define ZCOSIS (.91744867)
52#define ZSINIS (.39785416)
53#define ZSINGS (-.98088458)
54#define ZCOSGS (.1945905)
55#define ZCOSHS (1.0)
56#define ZSINHS (0.0)
57
58#define Q22 (1.7891679E-6)
59#define Q31 (2.1460748E-6)
60#define Q33 (2.2123015E-7)
61#define G22 (5.7686396)
62#define G32 (0.95240898)
63#define G44 (1.8014998)
64#define G52 (1.0508330)
65#define G54 (4.4108898)
66#define ROOT22 (1.7891679E-6)
67#define ROOT32 (3.7393792E-7)
68#define ROOT44 (7.3636953E-9)
69#define ROOT52 (1.1428639E-7)
70#define ROOT54 (2.1765803E-9)
71#define THDT (4.3752691E-3)
72
73#define IRESFL (sat->deep->deep_flags.IRESFL)
74#define ISYNFL (sat->deep->deep_flags.ISYNFL)
75
76#define s_SINIQ (sat->deep->deep_s_SINIQ)
77#define s_COSIQ (sat->deep->deep_s_COSIQ)
78#define s_OMGDT (sat->deep->deep_s_OMGDT)
79#define ATIME (sat->deep->deep_ATIME)
80#define D2201 (sat->deep->deep_D2201)
81#define D2211 (sat->deep->deep_D2211)
82#define D3210 (sat->deep->deep_D3210)
83#define D3222 (sat->deep->deep_D3222)
84#define D4410 (sat->deep->deep_D4410)
85#define D4422 (sat->deep->deep_D4422)
86#define D5220 (sat->deep->deep_D5220)
87#define D5232 (sat->deep->deep_D5232)
88#define D5421 (sat->deep->deep_D5421)
89#define D5433 (sat->deep->deep_D5433)
90#define DEL1 (sat->deep->deep_DEL1)
91#define DEL2 (sat->deep->deep_DEL2)
92#define DEL3 (sat->deep->deep_DEL3)
93#define E3 (sat->deep->deep_E3)
94#define EE2 (sat->deep->deep_EE2)
95#define FASX2 (sat->deep->deep_FASX2)
96#define FASX4 (sat->deep->deep_FASX4)
97#define FASX6 (sat->deep->deep_FASX6)
98#define OMEGAQ (sat->deep->deep_OMEGAQ)
99#define PE (sat->deep->deep_PE)
100#define PINC (sat->deep->deep_PINC)
101#define PL (sat->deep->deep_PL)
102#define SAVTSN (sat->deep->deep_SAVTSN)
103#define SE2 (sat->deep->deep_SE2)
104#define SE3 (sat->deep->deep_SE3)
105#define SGH2 (sat->deep->deep_SGH2)
106#define SGH3 (sat->deep->deep_SGH3)
107#define SGH4 (sat->deep->deep_SGH4)
108#define SGHL (sat->deep->deep_SGHL)
109#define SGHS (sat->deep->deep_SGHS)
110#define SH2 (sat->deep->deep_SH2)
111#define SH3 (sat->deep->deep_SH3)
112#define SHS (sat->deep->deep_SHS)
113#define SHL (sat->deep->deep_SHL)
114#define SI2 (sat->deep->deep_SI2)
115#define SI3 (sat->deep->deep_SI3)
116#define SL2 (sat->deep->deep_SL2)
117#define SL3 (sat->deep->deep_SL3)
118#define SL4 (sat->deep->deep_SL4)
119#define SSE (sat->deep->deep_SSE)
120#define SSG (sat->deep->deep_SSG)
121#define SSH (sat->deep->deep_SSH)
122#define SSI (sat->deep->deep_SSI)
123#define SSL (sat->deep->deep_SSL)
124#define STEP2 (sat->deep->deep_STEP2)
125#define STEPN (sat->deep->deep_STEPN)
126#define STEPP (sat->deep->deep_STEPP)
127#define THGR (sat->deep->deep_THGR)
128#define XFACT (sat->deep->deep_XFACT)
129#define XGH2 (sat->deep->deep_XGH2)
130#define XGH3 (sat->deep->deep_XGH3)
131#define XGH4 (sat->deep->deep_XGH4)
132#define XH2 (sat->deep->deep_XH2)
133#define XH3 (sat->deep->deep_XH3)
134#define XI2 (sat->deep->deep_XI2)
135#define XI3 (sat->deep->deep_XI3)
136#define XL2 (sat->deep->deep_XL2)
137#define XL3 (sat->deep->deep_XL3)
138#define XL4 (sat->deep->deep_XL4)
139#define XLAMO (sat->deep->deep_XLAMO)
140#define XLI (sat->deep->deep_XLI)
141#define XNI (sat->deep->deep_XNI)
142#define XNQ (sat->deep->deep_XNQ)
143#define XQNCL (sat->deep->deep_XQNCL)
144#define ZMOL (sat->deep->deep_ZMOL)
145#define ZMOS (sat->deep->deep_ZMOS)
146
147/* * ENTRANCE FOR DEEP SPACE INITIALIZATION */
148
149/* ENTRY DPINIT(EQSQ,SINIQ,COSIQ,RTEQSQ,AO,COSQ2,SINOMO,COSOMO, */
150/* 1 BSQ,XLLDOT,OMGDT,XNODOT,XNODP) */
151
152void
153dpinit(SatData *sat, double EQSQ, double SINIQ, double COSIQ,
154 double RTEQSQ, double AO, double COSQ2, double SINOMO,
155 double COSOMO, double BSQ, double XLLDOT, double OMGDT,
156 double XNODOT, double XNODP)
157{
158 double A1, A10, A2, A3, A4, A5, A6, A7, A8, A9, AINV2, AQNV, BFACT,
159 C, CC, COSQ, CTEM, DAY, DS50, EOC, EQ, F220, F221, F311, F321, F322,
160 F330, F441, F442, F522, F523, F542, F543, G200, G201, G211, G300,
161 G310, G322, G410, G422, G520, G521, G532, G533, GAM, PREEP, S1, S2,
162 S3, S4, S5, S6, S7, SE, SGH, SH, SI, SINI2, SINQ, SL, STEM, TEMP,
163 TEMP1, X1, X2, X3, X4, X5, X6, X7, X8, XMAO, XNO2, XNODCE, XNOI,
164 XPIDOT, Z1, Z11, Z12, Z13, Z2, Z21, Z22, Z23, Z3, Z31, Z32, Z33,
165 ZCOSG, ZCOSGL, ZCOSH, ZCOSHL, ZCOSI, ZCOSIL, ZE, ZMO, ZN, ZSING,
166 ZSINGL, ZSINH, ZSINHL, ZSINI, ZSINIL, ZX, ZY;
167
168 int c;
169#if 0
170 A1=A10=A2=A3=A4=A5=A6=A7=A8=A9=AINV2=AQNV=BFACT = signaling_nan();
171 C=CC=COSQ=CTEM=DAY=DS50=EOC=EQ=F220=F221=F311=F321=F322 = signaling_nan();
172 F330=F441=F442=F522=F523=F542=F543=G200=G201=G211=G300 = signaling_nan();
173 G310=G322=G410=G422=G520=G521=G532=G533=GAM=PREEP=S1=S2 = signaling_nan();
174 S3=S4=S5=S6=S7=SE=SGH=SH=SI=SINI2=SINQ=SL=STEM=TEMP = signaling_nan();
175 TEMP1=X1=X2=X3=X4=X5=X6=X7=X8=XMAO=XNO2=XNODCE=XNOI = signaling_nan();
176 XPIDOT=Z1=Z11=Z12=Z13=Z2=Z21=Z22=Z23=Z3=Z31=Z32=Z33 = signaling_nan();
177 ZCOSG=ZCOSGL=ZCOSH=ZCOSHL=ZCOSI=ZCOSIL=ZE=ZMO=ZN=ZSING = signaling_nan();
178 ZSINGL=ZSINH=ZSINHL=ZSINI=ZSINIL=ZX=ZY = signaling_nan();
179#endif
180 if(!sat->deep)
181 sat->deep = (struct deep_data *) malloc(sizeof(struct deep_data));
182 else
183 return;
184
185 /* init_deep(sat->deep); */
186 PREEP = 0.0;
187
188 ZCOSGL = ZCOSHL = ZCOSIL = ZSINGL = ZSINHL = ZSINIL = 0.0;
189
190 /* Save some of the arguments, for use by dpsec() and dpper() */
191 s_SINIQ = SINIQ;
192 s_COSIQ = COSIQ;
193 s_OMGDT = OMGDT;
194
195 THGR = thetag(EPOCH, &DS50);
196
197 EQ = EO;
198 XNQ = XNODP;
199 AQNV = 1.0/AO;
200 XQNCL = XINCL;
201 XMAO = XMO;
202 XPIDOT = OMGDT + XNODOT;
203 SINQ = sin(XNODEO);
204 COSQ = cos(XNODEO);
205 OMEGAQ = OMEGAO;
206
207 /* INITIALIZE LUNAR SOLAR TERMS */
208
209 DAY = DS50 + 18261.5;
210
211 if(DAY != PREEP) {
212 PREEP = DAY;
213 XNODCE = 4.5236020 - 9.2422029E-4 * DAY;
214 STEM = sin(XNODCE);
215 CTEM = cos(XNODCE);
216 ZCOSIL = .91375164 - .03568096 * CTEM;
217 ZSINIL = sqrt(1.0 - ZCOSIL * ZCOSIL);
218 ZSINHL = .089683511 * STEM / ZSINIL;
219 ZCOSHL = sqrt(1.0 - ZSINHL * ZSINHL);
220 C = 4.7199672 + .22997150 * DAY;
221 GAM = 5.8351514 + .0019443680 * DAY;
222 ZMOL = fmod(C-GAM, TWOPI);
223 ZX = .39785416 * STEM / ZSINIL;
224 ZY = ZCOSHL * CTEM + 0.91744867 * ZSINHL * STEM;
225 ZX = actan(ZX, ZY);
226 ZX = GAM + ZX - XNODCE;
227 ZCOSGL = cos(ZX);
228 ZSINGL = sin(ZX);
229 ZMOS = 6.2565837 + .017201977 * DAY;
230 ZMOS = fmod(ZMOS, TWOPI);
231 }
232
233 /* DO SOLAR TERMS */
234
235 SAVTSN = 1.0E20;
236 ZCOSG = ZCOSGS;
237 ZSING = ZSINGS;
238 ZCOSI = ZCOSIS;
239 ZSINI = ZSINIS;
240 ZCOSH = COSQ;
241 ZSINH = SINQ;
242 CC = C1SS;
243 ZN = ZNS;
244 ZE = ZES;
245 ZMO = ZMOS;
246 XNOI = 1.0 / XNQ;
247
248 for(c = 0; c < 2; c++) {
249 A1 = ZCOSG * ZCOSH + ZSING * ZCOSI * ZSINH;
250 A3 = -ZSING * ZCOSH + ZCOSG * ZCOSI * ZSINH;
251 A7 = -ZCOSG * ZSINH + ZSING * ZCOSI * ZCOSH;
252 A8 = ZSING * ZSINI;
253 A9 = ZSING * ZSINH + ZCOSG * ZCOSI * ZCOSH;
254 A10 = ZCOSG * ZSINI;
255 A2 = COSIQ * A7 + SINIQ * A8;
256 A4 = COSIQ * A9 + SINIQ * A10;
257 A5 = - SINIQ * A7 + COSIQ * A8;
258 A6 = - SINIQ * A9 + COSIQ * A10;
259
260 X1 = A1 * COSOMO + A2 * SINOMO;
261 X2 = A3 * COSOMO + A4 * SINOMO;
262 X3 = - A1 * SINOMO + A2 * COSOMO;
263 X4 = - A3 * SINOMO + A4 * COSOMO;
264 X5 = A5 * SINOMO;
265 X6 = A6 * SINOMO;
266 X7 = A5 * COSOMO;
267 X8 = A6 * COSOMO;
268
269 Z31 = 12.0 * X1 * X1 -3.0 * X3 * X3;
270 Z32 = 24.0 * X1 * X2 -6.0 * X3 * X4;
271 Z33 = 12.0 * X2 * X2 -3.0 * X4 * X4;
272 Z1 = 3.0 * (A1 * A1 + A2 * A2) + Z31 * EQSQ;
273 Z2 = 6.0 * (A1 * A3 + A2 * A4) + Z32 * EQSQ;
274 Z3 = 3.0 * (A3 * A3 + A4 * A4) + Z33 * EQSQ;
275 Z11 = -6.0 * A1 * A5 + EQSQ * (-24.0 * X1 * X7 - 6.0 * X3 * X5);
276
277 Z12 = -6.0 * (A1 * A6 + A3 * A5) +
278 EQSQ * (-24.0 * (X2 * X7 + X1 * X8) - 6.0 * (X3 * X6 + X4 * X5));
279
280 Z13 = -6.0 * A3 * A6 + EQSQ * (-24.0 * X2 * X8 - 6.0 * X4 * X6);
281 Z21 = 6.0 * A2 * A5 + EQSQ * (24.0 * X1 * X5 - 6.0 * X3 * X7);
282
283 Z22 = 6.0 * (A4 * A5 + A2 * A6) +
284 EQSQ * (24.0 * (X2 * X5 + X1 * X6) - 6.0 * (X4 * X7 + X3 * X8));
285
286 Z23 = 6.0 * A4 * A6 + EQSQ * (24.0 * X2 * X6 - 6.0 * X4 * X8);
287 Z1 = Z1 + Z1 + BSQ * Z31;
288 Z2 = Z2 + Z2 + BSQ * Z32;
289 Z3 = Z3 + Z3 + BSQ * Z33;
290 S3 = CC * XNOI;
291 S2 = -.5 * S3 / RTEQSQ;
292 S4 = S3 * RTEQSQ;
293 S1 = -15.0 * EQ * S4;
294 S5 = X1 * X3 + X2 * X4;
295 S6 = X2 * X3 + X1 * X4;
296 S7 = X2 * X4 - X1 * X3;
297 SE = S1 * ZN * S5;
298 SI = S2 * ZN * (Z11 + Z13);
299 SL = -ZN * S3 * (Z1 + Z3 - 14.0 - 6.0 * EQSQ);
300 SGH = S4 * ZN * (Z31 + Z33 - 6.0);
301 SH = -ZN * S2 * (Z21 + Z23);
302
303 if(XQNCL < 5.2359877E-2)
304 SH = 0.0;
305
306 EE2 = 2.0 * S1 * S6;
307 E3 = 2.0 * S1 * S7;
308 XI2 = 2.0 * S2 * Z12;
309 XI3 = 2.0 * S2 * (Z13 - Z11);
310 XL2 = -2.0 * S3 * Z2;
311 XL3 = -2.0 * S3 * (Z3 - Z1);
312 XL4 = -2.0 * S3 * (-21.0 - 9.0 * EQSQ) * ZE;
313 XGH2 = 2.0 * S4 * Z32;
314 XGH3 = 2.0 * S4 * (Z33 - Z31);
315 XGH4 = -18.0 * S4 * ZE;
316 XH2 = -2.0 * S2 * Z22;
317 XH3 = -2.0 * S2 * (Z23 - Z21);
318
319 if(c == 0) {
320 /* DO LUNAR TERMS */
321 SSE = SE;
322 SSI = SI;
323 SSL = SL;
324 SSH = SH / SINIQ;
325 SSG = SGH - COSIQ * SSH;
326 SE2 = EE2;
327 SI2 = XI2;
328 SL2 = XL2;
329 SGH2 = XGH2;
330 SH2 = XH2;
331 SE3 = E3;
332 SI3 = XI3;
333 SL3 = XL3;
334 SGH3 = XGH3;
335 SH3 = XH3;
336 SL4 = XL4;
337 SGH4 = XGH4;
338
339 ZCOSG = ZCOSGL;
340 ZSING = ZSINGL;
341 ZCOSI = ZCOSIL;
342 ZSINI = ZSINIL;
343 ZCOSH = ZCOSHL * COSQ + ZSINHL * SINQ;
344 ZSINH = SINQ * ZCOSHL - COSQ * ZSINHL;
345 ZN = ZNL;
346 CC = C1L;
347 ZE = ZEL;
348 ZMO = ZMOL;
349 }
350 }
351
352 SSE = SSE + SE;
353 SSI = SSI + SI;
354 SSL = SSL + SL;
355 SSG = SSG + SGH - COSIQ / SINIQ * SH;
356 SSH = SSH + SH / SINIQ;
357
358 /* GEOPOTENTIAL RESONANCE INITIALIZATION FOR 12 HOUR ORBITS */
359
360 IRESFL = 0;
361 ISYNFL = 0;
362
363 if(XNQ <= .0034906585 || XNQ >= .0052359877) {
364
365 if(XNQ < (8.26E-3) || XNQ > (9.24E-3))
366 return;
367
368 if(EQ < 0.5)
369 return;
370
371 IRESFL = 1;
372 EOC = EQ * EQSQ;
373 G201 = -.306 - (EQ - .64) * .440;
374
375 if(EQ <= (.65)) {
376 G211 = 3.616 - 13.247 * EQ + 16.290 * EQSQ;
377 G310 = -19.302 + 117.390 * EQ - 228.419 * EQSQ + 156.591 * EOC;
378 G322 = -18.9068 + 109.7927 * EQ - 214.6334 * EQSQ + 146.5816 * EOC;
379 G410 = -41.122 + 242.694 * EQ - 471.094 * EQSQ + 313.953 * EOC;
380 G422 = -146.407 + 841.880 * EQ - 1629.014 * EQSQ + 1083.435 * EOC;
381 G520 = -532.114 + 3017.977 * EQ - 5740 * EQSQ + 3708.276 * EOC;
382 } else {
383 G211 = -72.099 + 331.819 * EQ - 508.738 * EQSQ + 266.724 * EOC;
384 G310 = -346.844 + 1582.851 * EQ - 2415.925 * EQSQ + 1246.113 * EOC;
385 G322 = -342.585 + 1554.908 * EQ - 2366.899 * EQSQ + 1215.972 * EOC;
386 G410 = -1052.797 + 4758.686 * EQ - 7193.992 * EQSQ +
387 3651.957 * EOC;
388 G422 = -3581.69 + 16178.11 * EQ - 24462.77 * EQSQ + 12422.52 * EOC;
389
390 if(EQ > (.715))
391 G520 = -5149.66 + 29936.92 * EQ - 54087.36 * EQSQ +
392 31324.56 * EOC;
393
394 G520 = 1464.74 - 4664.75 * EQ + 3763.64 * EQSQ;
395 }
396
397 if(EQ < (.7)) {
398 G533 = -919.2277 + 4988.61 * EQ - 9064.77 * EQSQ + 5542.21 * EOC;
399
400 G521 = -822.71072 + 4568.6173 * EQ - 8491.4146 * EQSQ +
401 5337.524 * EOC;
402
403 G532 = -853.666 + 4690.25 * EQ - 8624.77 * EQSQ + 5341.4 * EOC;
404 } else {
405 G533 = -37995.78 + 161616.52 * EQ - 229838.2 * EQSQ +
406 109377.94 * EOC;
407
408 G521 = -51752.104 + 218913.95 * EQ - 309468.16 * EQSQ +
409 146349.42 * EOC;
410
411 G532 = -40023.88 + 170470.89 * EQ - 242699.48 * EQSQ +
412 115605.82 * EOC;
413 }
414
415 SINI2 = SINIQ * SINIQ;
416 F220 = .75 * (1.0 + 2.0 * COSIQ + COSQ2);
417 F221 = 1.5 * SINI2;
418 F321 = 1.875 * SINIQ * (1.0 - 2.0 * COSIQ - 3.0 * COSQ2);
419 F322 = -1.875 * SINIQ * (1.0 + 2.0 * COSIQ - 3.0 * COSQ2);
420 F441 = 35.0 * SINI2 * F220;
421 F442 = 39.3750 * SINI2 * SINI2;
422
423 F522 = 9.84375 * SINIQ * (SINI2 * (1.0 - 2.0 * COSIQ - 5.0 * COSQ2) +
424 .33333333 * (-2.0 + 4.0 * COSIQ +
425 6.0 * COSQ2));
426
427 F523 = SINIQ * (4.92187512 * SINI2 * (-2.0 - 4.0 * COSIQ +
428 10.0 * COSQ2) +
429 6.56250012 * (1.0 +
430 2.0 * COSIQ -
431 3.0 * COSQ2));
432
433 F542 = 29.53125 * SINIQ * (2.0 - 8.0 * COSIQ +
434 COSQ2 * (-12.0 + 8.0 * COSIQ +
435 10.0 * COSQ2));
436
437 F543 = 29.53125 * SINIQ * (-2.0 - 8.0 * COSIQ +
438 COSQ2 * (12.0 + 8.0 * COSIQ -
439 10.0 * COSQ2));
440
441 XNO2 = XNQ * XNQ;
442 AINV2 = AQNV * AQNV;
443 TEMP1 = 3.0 * XNO2 * AINV2;
444 TEMP = TEMP1 * ROOT22;
445 D2201 = TEMP * F220 * G201;
446 D2211 = TEMP * F221 * G211;
447 TEMP1 = TEMP1 * AQNV;
448 TEMP = TEMP1 * ROOT32;
449 D3210 = TEMP * F321 * G310;
450 D3222 = TEMP * F322 * G322;
451 TEMP1 = TEMP1 * AQNV;
452 TEMP = 2.0 * TEMP1 * ROOT44;
453 D4410 = TEMP * F441 * G410;
454 D4422 = TEMP * F442 * G422;
455 TEMP1 = TEMP1 * AQNV;
456 TEMP = TEMP1 * ROOT52;
457 D5220 = TEMP * F522 * G520;
458 D5232 = TEMP * F523 * G532;
459 TEMP = 2.0 * TEMP1 * ROOT54;
460 D5421 = TEMP * F542* G521;
461 D5433 = TEMP * F543* G533;
462 XLAMO = XMAO + XNODEO + XNODEO - THGR - THGR;
463 BFACT = XLLDOT + XNODOT + XNODOT - THDT - THDT;
464 BFACT = BFACT + SSL + SSH + SSH;
465 } else {
466 /* SYNCHRONOUS RESONANCE TERMS INITIALIZATION */
467
468 IRESFL = 1;
469 ISYNFL = 1;
470 G200 = 1.0 + EQSQ * (-2.5 + .8125 * EQSQ);
471 G310 = 1.0 + 2.0 * EQSQ;
472 G300 = 1.0 + EQSQ * (-6.0 + 6.60937 * EQSQ);
473 F220 = .75 * (1.0 + COSIQ) * (1.0 + COSIQ);
474 F311 = .9375 * SINIQ * SINIQ * (1.0 + 3.0 * COSIQ) -
475 .75 * (1.0 + COSIQ);
476 F330 = 1.0 + COSIQ;
477 F330 = 1.875 * F330 * F330 * F330;
478 DEL1 = 3.0 * XNQ * XNQ * AQNV * AQNV;
479 DEL2 = 2.0 * DEL1 * F220 * G200 * Q22;
480 DEL3 = 3.0 * DEL1 * F330 * G300 * Q33 * AQNV;
481 DEL1 = DEL1 * F311 * G310 * Q31 * AQNV;
482 FASX2 = .13130908;
483 FASX4 = 2.8843198;
484 FASX6 = .37448087;
485 XLAMO = XMAO + XNODEO + OMEGAO - THGR;
486 BFACT = XLLDOT + XPIDOT - THDT;
487 BFACT = BFACT + SSL + SSG + SSH;
488
489 }
490
491 XFACT = BFACT - XNQ;
492
493 XLI = XLAMO;
494 XNI = XNQ;
495 ATIME = 0.0;
496 STEPP = 720.0;
497 STEPN = -720.0;
498 STEP2 = 259200.0;
499}
500
501/* ENTRANCE FOR DEEP SPACE SECULAR EFFECTS */
502
503void
504dpsec(SatData *sat, double *XLL, double *OMGASM, double *XNODES,
505 double *EM, double *XINC, double *XN, double T)
506{
507 double DELT, XL, TEMP, XOMI, X2OMI, X2LI, XLDOT;
508 double XNDOT, XNDDT, FT;
509 int state, iret, iretn, done;
510
511 DELT = XLDOT = XNDOT = XNDDT = FT = 0.0;
512 iret = iretn = 0;
513
514#if 0
515 DELT = XL = TEMP = XOMI = X2OMI = X2LI = XLDOT = signaling_nan();
516 XNDOT = XNDDT = FT = signaling_nan();
517#endif
518
519 *XLL = *XLL + SSL * T;
520 *OMGASM = *OMGASM + SSG * T;
521 *XNODES = *XNODES + SSH * T;
522 *EM = EO + SSE * T;
523 *XINC = XINCL + SSI * T;
524
525 if(*XINC < 0.0) {
526 *XINC = -*XINC;
527 *XNODES = *XNODES + PI;
528 *OMGASM = *OMGASM - PI;
529 }
530
531 if(IRESFL == 0)
532 return;
533
534 state = 1;
535 done = 0;
536 while(!done) {
537 /* printf("state = %d\n", state); */
538 switch(state) {
539 case 1:
540 /*
541 * Chunk #1
542 */
543 if(ATIME == 0.0 || (T >= 0.0 && ATIME < 0.0) ||
544 (T < 0.0 && ATIME >= 0.0)) {
545 /*
546 * Chunk #10
547 */
548 if(T >= 0.0)
549 DELT = STEPP;
550 else
551 DELT = STEPN;
552
553 ATIME = 0.0;
554 XNI = XNQ;
555 XLI = XLAMO;
556 state = 4;
557 break;
558 }
559
560 /* Fall through */
561 case 2:
562 /*
563 * Chunk #2
564 */
565 if(fabs(T) < fabs(ATIME)) {
566 /*
567 * Chunk #2
568 */
569 if(T >= 0.0)
570 DELT = STEPN;
571 else
572 DELT = STEPP;
573
574 iret = 1;
575 state = 8;
576 break;
577 }
578
579 /*
580 * Chunk #3
581 */
582 if(T > 0.0)
583 DELT = STEPP;
584 else
585 DELT = STEPN;
586
587 /* fall through */
588 case 4:
589 /*
590 * Chunk #4
591 */
592 if(fabs(T - ATIME) >= STEPP) {
593 iret = 4;
594 state = 8;
595 } else {
596 /*
597 * Chunk #5
598 */
599 FT = T - ATIME;
600 iretn = 6;
601 state = 7;
602 }
603
604 break;
605
606 case 6:
607 /*
608 * Chunk #6
609 */
610 *XN = XNI + XNDOT * FT + XNDDT * FT * FT * 0.5;
611 XL = XLI + XLDOT * FT + XNDOT * FT * FT * 0.5;
612 TEMP = -*XNODES + THGR + T * THDT;
613
614 if(ISYNFL == 0)
615 *XLL = XL + 2.0 * TEMP;
616 else
617 *XLL = XL - *OMGASM + TEMP;
618
619 done = 1;
620 break;
621
622 case 7:
623 /* DOT TERMS CALCULATED */
624
625 /*
626 * Chunk #7
627 */
628 if(ISYNFL != 0) {
629 XNDOT =
630 DEL1 * sin(XLI - FASX2) +
631 DEL2 * sin(2.0 * (XLI - FASX4)) +
632 DEL3 * sin(3.0 * (XLI - FASX6));
633
634 XNDDT =
635 DEL1 * cos(XLI - FASX2) +
636 2.0 * DEL2 * cos(2.0 * (XLI - FASX4)) +
637 3.0 * DEL3 * cos(3.0 * (XLI - FASX6));
638 } else {
639 XOMI = OMEGAQ + s_OMGDT * ATIME;
640 X2OMI = XOMI + XOMI;
641 X2LI = XLI + XLI;
642 XNDOT = D2201 * sin(X2OMI + XLI - G22) +
643 D2211 * sin(XLI - G22) +
644 D3210 * sin(XOMI + XLI - G32) +
645 D3222 * sin(- XOMI + XLI - G32) +
646 D4410 * sin(X2OMI + X2LI - G44) +
647 D4422 * sin(X2LI - G44) +
648 D5220 * sin(XOMI + XLI - G52) +
649 D5232 * sin(- XOMI + XLI - G52) +
650 D5421 * sin(XOMI + X2LI - G54) +
651 D5433 * sin(- XOMI + X2LI - G54);
652
653 XNDDT = D2201 * cos(X2OMI + XLI - G22) +
654 D2211 * cos(XLI - G22) +
655 D3210 * cos(XOMI + XLI - G32) +
656 D3222 * cos(- XOMI + XLI - G32) +
657 D5220 * cos(XOMI + XLI - G52) +
658 D5232 * cos(- XOMI + XLI - G52) +
659 2.*(D4410 * cos(X2OMI + X2LI - G44) +
660 D4422 * cos(X2LI - G44) +
661 D5421 * cos(XOMI + X2LI - G54) +
662 D5433 * cos(- XOMI + X2LI - G54));
663 }
664
665 XLDOT = XNI + XFACT;
666 XNDDT = XNDDT * XLDOT;
667
668 state = iretn;
669 break;
670
671 case 8:
672 /*
673 * Chunk #8
674 */
675
676 /* INTEGRATOR */
677 iretn = 9;
678 state = 7;
679 break;
680
681 case 9:
682 XLI = XLI + XLDOT * DELT + XNDOT * STEP2;
683 XNI = XNI + XNDOT * DELT + XNDDT * STEP2;
684 ATIME = ATIME + DELT;
685
686 state = iret;
687 break;
688 }
689 }
690}
691
692/* local */
693
694/* C */
695/* C ENTRANCES FOR LUNAR-SOLAR PERIODICS */
696/* C */
697/* C */
698/* ENTRY DPPER(EM,XINC,OMGASM,XNODES,XLL) */
699void
700dpper(SatData *sat, double *EM, double *XINC, double *OMGASM,
701 double *XNODES, double *XLL, double T)
702{
703 double SINIS, COSIS, ZM, ZF, SINZF, F2, F3, SES, SIS, SLS, SEL, SIL, SLL, PGH, PH, SINOK, COSOK, ALFDP, BETDP, DALF, DBET, XLS, DLS;
704
705#if 0
706 SINIS = COSIS = ZM = ZF = SINZF = F2 = F3 = SES = SIS = signaling_nan();
707 SLS = SEL = SIL = SLL = PGH = signaling_nan();
708 PH = SINOK = COSOK = ALFDP = BETDP = DALF = DBET = XLS = signaling_nan();
709 DLS = signaling_nan();;
710#endif
711 SINIS = sin(*XINC);
712 COSIS = cos(*XINC);
713
714
715/* IF (DABS(SAVTSN-T).LT.(30.D0)) GO TO 210 */
716 if(fabs(SAVTSN - T) >= (30.0)) {
717 SAVTSN = T;
718 ZM = ZMOS + ZNS * T;
719/* 205 ZF = ZM + 2.0 * ZES * sin(ZM) */
720 ZF = ZM + 2.0 * ZES * sin(ZM);
721 SINZF = sin(ZF);
722 F2 = .5 * SINZF * SINZF - .25;
723 F3 = -.5 * SINZF * cos(ZF);
724 SES = SE2 * F2 + SE3 * F3;
725 SIS = SI2 * F2 + SI3 * F3;
726 SLS = SL2 * F2 + SL3 * F3 + SL4 * SINZF;
727 SGHS = SGH2 * F2 + SGH3 * F3 + SGH4 * SINZF;
728 SHS = SH2 * F2 + SH3 * F3;
729 ZM = ZMOL + ZNL * T;
730 ZF = ZM + 2.0 * ZEL * sin(ZM);
731 SINZF = sin(ZF);
732 F2 = .5 * SINZF * SINZF -.25;
733 F3 = -.5 * SINZF * cos(ZF);
734 SEL = EE2 * F2 + E3 * F3;
735 SIL = XI2 * F2 + XI3 * F3;
736 SLL = XL2 * F2 + XL3 * F3 + XL4 * SINZF;
737 SGHL = XGH2 * F2 + XGH3 * F3 + XGH4 * SINZF;
738 SHL = XH2 * F2 + XH3 * F3;
739 PE = SES + SEL;
740 PINC = SIS + SIL;
741 PL = SLS + SLL;
742 }
743
744/* 210 PGH=SGHS+SGHL */
745 PGH = SGHS + SGHL;
746 PH = SHS + SHL;
747 *XINC = *XINC + PINC;
748 *EM = *EM + PE;
749
750/* IF(XQNCL.LT.(.2)) GO TO 220 */
751 if(XQNCL >= (.2)) {
752/* GO TO 218 */
753/* C */
754/* C APPLY PERIODICS DIRECTLY */
755/* C */
756/* 218 PH=PH/SINIQ */
757 PH = PH / s_SINIQ;
758 PGH = PGH - s_COSIQ * PH;
759 *OMGASM = *OMGASM + PGH;
760 *XNODES = *XNODES + PH;
761 *XLL = *XLL + PL;
762/* GO TO 230 */
763 } else {
764/* C */
765/* C APPLY PERIODICS WITH LYDDANE MODIFICATION */
766/* C */
767/* 220 SINOK=sin(XNODES) */
768 SINOK = sin(*XNODES);
769 COSOK = cos(*XNODES);
770 ALFDP = SINIS * SINOK;
771 BETDP = SINIS * COSOK;
772 DALF = PH * COSOK + PINC * COSIS * SINOK;
773 DBET = -PH * SINOK + PINC * COSIS * COSOK;
774 ALFDP = ALFDP + DALF;
775 BETDP = BETDP + DBET;
776 XLS = *XLL + *OMGASM + COSIS * *XNODES;
777 DLS = PL + PGH - PINC * *XNODES * SINIS;
778 XLS = XLS + DLS;
779 *XNODES = actan(ALFDP, BETDP);
780 *XLL = *XLL + PL;
781 *OMGASM = XLS - *XLL - cos(*XINC) * *XNODES;
782 }
783/* 230 CONTINUE */
784/* RETURN */
785
786}
787/* END */
788
789/* For RCS Only -- Do Not Edit */
[3477]790static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: deep.c,v $ $Date: 2008-03-25 17:45:13 $ $Revision: 1.7 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.