← Back to branch summary

~ubuntu-branches/ubuntu/warty/xplanet/warty

« back to all changes in this revision

Viewing changes to src/libsgp4sdp4/sgp4sdp4.c

  • Committer: Bazaar Package Importer
  • Author(s): LaMont Jones
  • Date: 2004-08-24 07:14:00 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040824071400-2dr4qnjbjmm8z3ia
Tags: 1.0.6-1ubuntu1
Build-depend: libtiff4-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/libsgp4sdp4/sgp4sdp4.c
 
1
/*
 
2
 *  Unit SGP4SDP4
 
3
 *           Author:  Dr TS Kelso 
 
4
 * Original Version:  1991 Oct 30
 
5
 * Current Revision:  1992 Sep 03
 
6
 *          Version:  1.50 
 
7
 *        Copyright:  1991-1992, All Rights Reserved 
 
8
 *
 
9
 *   Ported to C by:  Neoklis Kyriazis  April 10  2001
 
10
 */
 
11
 
 
12
#include "sgp4sdp4.h"
 
13
 
 
14
/* SGP4 */
 
15
/* This function is used to calculate the position and velocity */
 
16
/* of near-earth (period < 225 minutes) satellites. tsince is   */
 
17
/* time since epoch in minutes, tle is a pointer to a tle_t     */
 
18
/* structure with Keplerian orbital elements and pos and vel    */
 
19
/* are vector_t structures returning ECI satellite position and */
 
20
/* velocity. Use Convert_Sat_State() to convert to km and km/s.*/
 
21
void
 
22
SGP4(double tsince, tle_t *tle, vector_t *pos, vector_t *vel)
 
23
{
 
24
  static double
 
25
    aodp,aycof,c1,c4,c5,cosio,d2,d3,d4,delmo,omgcof,
 
26
    eta,omgdot,sinio,xnodp,sinmo,t2cof,t3cof,t4cof,t5cof,
 
27
    x1mth2,x3thm1,x7thm1,xmcof,xmdot,xnodcf,xnodot,xlcof;
 
28
 
 
29
  double
 
30
    cosuk,sinuk,rfdotk,vx,vy,vz,ux,uy,uz,xmy,xmx,
 
31
    cosnok,sinnok,cosik,sinik,rdotk,xinck,xnodek,uk,
 
32
    rk,cos2u,sin2u,u,sinu,cosu,betal,rfdot,rdot,r,pl,
 
33
    elsq,esine,ecose,epw,cosepw,x1m5th,xhdot1,tfour,
 
34
    sinepw,capu,ayn,xlt,aynl,xll,axn,xn,beta,xl,e,a,
 
35
    tcube,delm,delomg,templ,tempe,tempa,xnode,tsq,xmp,
 
36
    omega,xnoddf,omgadf,xmdf,a1,a3ovk2,ao,betao,betao2,
 
37
    c1sq,c2,c3,coef,coef1,del1,delo,eeta,eosq,etasq,
 
38
    perige,pinvsq,psisq,qoms24,s4,temp,temp1,temp2,
 
39
    temp3,temp4,temp5,temp6,theta2,theta4,tsi;
 
40
 
 
41
  int i;  
 
42
 
 
43
  /* Initialization */
 
44
  if (isFlagClear(SGP4_INITIALIZED_FLAG))
 
45
    {
 
46
      SetFlag(SGP4_INITIALIZED_FLAG);
 
47
 
 
48
      /* Recover original mean motion (xnodp) and   */
 
49
      /* semimajor axis (aodp) from input elements. */
 
50
      a1 = pow(xke/tle->xno,tothrd);
 
51
      cosio = cos(tle->xincl);
 
52
      theta2 = cosio*cosio;
 
53
      x3thm1 = 3*theta2-1.0;
 
54
      eosq = tle->eo*tle->eo;
 
55
      betao2 = 1-eosq;
 
56
      betao = sqrt(betao2);
 
57
      del1 = 1.5*ck2*x3thm1/(a1*a1*betao*betao2);
 
58
      ao = a1*(1-del1*(0.5*tothrd+del1*(1+134/81*del1)));
 
59
      delo = 1.5*ck2*x3thm1/(ao*ao*betao*betao2);
 
60
      xnodp = tle->xno/(1+delo);
 
61
      aodp = ao/(1-delo);
 
62
 
 
63
      /* For perigee less than 220 kilometers, the "simple" flag is set */
 
64
      /* and the equations are truncated to linear variation in sqrt a  */
 
65
      /* and quadratic variation in mean anomaly.  Also, the c3 term,   */
 
66
      /* the delta omega term, and the delta m term are dropped.        */
 
67
      if((aodp*(1-tle->eo)/ae) < (220/xkmper+ae))
 
68
        SetFlag(SIMPLE_FLAG);
 
69
      else
 
70
        ClearFlag(SIMPLE_FLAG);
 
71
 
 
72
      /* For perigee below 156 km, the       */ 
 
73
      /* values of s and qoms2t are altered. */
 
74
      s4 = s;
 
75
      qoms24 = qoms2t;
 
76
      perige = (aodp*(1-tle->eo)-ae)*xkmper;
 
77
      if(perige < 156)
 
78
        {
 
79
          if(perige <= 98)
 
80
            s4 = 20;
 
81
          else
 
82
            s4 = perige-78;
 
83
          qoms24 = pow((120-s4)*ae/xkmper,4);
 
84
          s4 = s4/xkmper+ae;
 
85
        }; /* End of if(perige <= 98) */
 
86
 
 
87
      pinvsq = 1/(aodp*aodp*betao2*betao2);
 
88
      tsi = 1/(aodp-s4);
 
89
      eta = aodp*tle->eo*tsi;
 
90
      etasq = eta*eta;
 
91
      eeta = tle->eo*eta;
 
92
      psisq = fabs(1-etasq);
 
93
      coef = qoms24*pow(tsi,4);
 
94
      coef1 = coef/pow(psisq,3.5);
 
95
      c2 = coef1*xnodp*(aodp*(1+1.5*etasq+eeta*(4+etasq))+
 
96
           0.75*ck2*tsi/psisq*x3thm1*(8+3*etasq*(8+etasq)));
 
97
      c1 = tle->bstar*c2;
 
98
      sinio = sin(tle->xincl);
 
99
      a3ovk2 = -xj3/ck2*pow(ae,3);
 
100
      c3 = coef*tsi*a3ovk2*xnodp*ae*sinio/tle->eo;
 
101
      x1mth2 = 1-theta2;
 
102
      c4 = 2*xnodp*coef1*aodp*betao2*(eta*(2+0.5*etasq)+
 
103
           tle->eo*(0.5+2*etasq)-2*ck2*tsi/(aodp*psisq)*
 
104
           (-3*x3thm1*(1-2*eeta+etasq*(1.5-0.5*eeta))+0.75*
 
105
           x1mth2*(2*etasq-eeta*(1+etasq))*cos(2*tle->omegao)));
 
106
      c5 = 2*coef1*aodp*betao2*(1+2.75*(etasq+eeta)+eeta*etasq);
 
107
      theta4 = theta2*theta2;
 
108
      temp1 = 3*ck2*pinvsq*xnodp;
 
109
      temp2 = temp1*ck2*pinvsq;
 
110
      temp3 = 1.25*ck4*pinvsq*pinvsq*xnodp;
 
111
      xmdot = xnodp+0.5*temp1*betao*x3thm1+
 
112
              0.0625*temp2*betao*(13-78*theta2+137*theta4);
 
113
      x1m5th = 1-5*theta2;
 
114
      omgdot = -0.5*temp1*x1m5th+0.0625*temp2*(7-114*theta2+
 
115
               395*theta4)+temp3*(3-36*theta2+49*theta4);
 
116
      xhdot1 = -temp1*cosio;
 
117
      xnodot = xhdot1+(0.5*temp2*(4-19*theta2)+
 
118
               2*temp3*(3-7*theta2))*cosio;
 
119
      omgcof = tle->bstar*c3*cos(tle->omegao);
 
120
      xmcof = -tothrd*coef*tle->bstar*ae/eeta;
 
121
      xnodcf = 3.5*betao2*xhdot1*c1;
 
122
      t2cof = 1.5*c1;
 
123
      xlcof = 0.125*a3ovk2*sinio*(3+5*cosio)/(1+cosio);
 
124
      aycof = 0.25*a3ovk2*sinio;
 
125
      delmo = pow(1+eta*cos(tle->xmo),3);
 
126
      sinmo = sin(tle->xmo);
 
127
      x7thm1 = 7*theta2-1;
 
128
      if (isFlagClear(SIMPLE_FLAG))
 
129
        {
 
130
          c1sq = c1*c1;
 
131
          d2 = 4*aodp*tsi*c1sq;
 
132
          temp = d2*tsi*c1/3;
 
133
          d3 = (17*aodp+s4)*temp;
 
134
          d4 = 0.5*temp*aodp*tsi*(221*aodp+31*s4)*c1;
 
135
          t3cof = d2+2*c1sq;
 
136
          t4cof = 0.25*(3*d3+c1*(12*d2+10*c1sq));
 
137
          t5cof = 0.2*(3*d4+12*c1*d3+6*d2*d2+15*c1sq*(2*d2+c1sq));
 
138
        }; /* End of if (isFlagClear(SIMPLE_FLAG)) */
 
139
    }; /* End of SGP4() initialization */
 
140
 
 
141
  /* Update for secular gravity and atmospheric drag. */
 
142
  xmdf = tle->xmo+xmdot*tsince;
 
143
  omgadf = tle->omegao+omgdot*tsince;
 
144
  xnoddf = tle->xnodeo+xnodot*tsince;
 
145
  omega = omgadf;
 
146
  xmp = xmdf;
 
147
  tsq = tsince*tsince;
 
148
  xnode = xnoddf+xnodcf*tsq;
 
149
  tempa = 1-c1*tsince;
 
150
  tempe = tle->bstar*c4*tsince;
 
151
  templ = t2cof*tsq;
 
152
  if (isFlagClear(SIMPLE_FLAG))
 
153
    {
 
154
      delomg = omgcof*tsince;
 
155
      delm = xmcof*(pow(1+eta*cos(xmdf),3)-delmo);
 
156
      temp = delomg+delm;
 
157
      xmp = xmdf+temp;
 
158
      omega = omgadf-temp;
 
159
      tcube = tsq*tsince;
 
160
      tfour = tsince*tcube;
 
161
      tempa = tempa-d2*tsq-d3*tcube-d4*tfour;
 
162
      tempe = tempe+tle->bstar*c5*(sin(xmp)-sinmo);
 
163
      templ = templ+t3cof*tcube+tfour*(t4cof+tsince*t5cof);
 
164
    }; /* End of if (isFlagClear(SIMPLE_FLAG)) */
 
165
 
 
166
  a = aodp*pow(tempa,2);
 
167
  e = tle->eo-tempe;
 
168
  xl = xmp+omega+xnode+xnodp*templ;
 
169
  beta = sqrt(1-e*e);
 
170
  xn = xke/pow(a,1.5);
 
171
 
 
172
  /* Long period periodics */
 
173
  axn = e*cos(omega);
 
174
  temp = 1/(a*beta*beta);
 
175
  xll = temp*xlcof*axn;
 
176
  aynl = temp*aycof;
 
177
  xlt = xl+xll;
 
178
  ayn = e*sin(omega)+aynl;
 
179
 
 
180
  /* Solve Kepler's' Equation */
 
181
  capu = FMod2p(xlt-xnode);
 
182
  temp2 = capu;
 
183
 
 
184
  i = 0;
 
185
  do
 
186
    {
 
187
      sinepw = sin(temp2);
 
188
      cosepw = cos(temp2);
 
189
      temp3 = axn*sinepw;
 
190
      temp4 = ayn*cosepw;
 
191
      temp5 = axn*cosepw;
 
192
      temp6 = ayn*sinepw;
 
193
      epw = (capu-temp4+temp3-temp2)/(1-temp5-temp6)+temp2;
 
194
      if(fabs(epw-temp2) <= e6a) break;
 
195
      temp2 = epw;
 
196
    }
 
197
  while( i++ < 10 );
 
198
 
 
199
  /* Short period preliminary quantities */
 
200
  ecose = temp5+temp6;
 
201
  esine = temp3-temp4;
 
202
  elsq = axn*axn+ayn*ayn;
 
203
  temp = 1-elsq;
 
204
  pl = a*temp;
 
205
  r = a*(1-ecose);
 
206
  temp1 = 1/r;
 
207
  rdot = xke*sqrt(a)*esine*temp1;
 
208
  rfdot = xke*sqrt(pl)*temp1;
 
209
  temp2 = a*temp1;
 
210
  betal = sqrt(temp);
 
211
  temp3 = 1/(1+betal);
 
212
  cosu = temp2*(cosepw-axn+ayn*esine*temp3);
 
213
  sinu = temp2*(sinepw-ayn-axn*esine*temp3);
 
214
  u = AcTan(sinu, cosu);
 
215
  sin2u = 2*sinu*cosu;
 
216
  cos2u = 2*cosu*cosu-1;
 
217
  temp = 1/pl;
 
218
  temp1 = ck2*temp;
 
219
  temp2 = temp1*temp;
 
220
 
 
221
  /* Update for short periodics */
 
222
  rk = r*(1-1.5*temp2*betal*x3thm1)+0.5*temp1*x1mth2*cos2u;
 
223
  uk = u-0.25*temp2*x7thm1*sin2u;
 
224
  xnodek = xnode+1.5*temp2*cosio*sin2u;
 
225
  xinck = tle->xincl+1.5*temp2*cosio*sinio*cos2u;
 
226
  rdotk = rdot-xn*temp1*x1mth2*sin2u;
 
227
  rfdotk = rfdot+xn*temp1*(x1mth2*cos2u+1.5*x3thm1);
 
228
 
 
229
  /* Orientation vectors */
 
230
  sinuk = sin(uk);
 
231
  cosuk = cos(uk);
 
232
  sinik = sin(xinck);
 
233
  cosik = cos(xinck);
 
234
  sinnok = sin(xnodek);
 
235
  cosnok = cos(xnodek);
 
236
  xmx = -sinnok*cosik;
 
237
  xmy = cosnok*cosik;
 
238
  ux = xmx*sinuk+cosnok*cosuk;
 
239
  uy = xmy*sinuk+sinnok*cosuk;
 
240
  uz = sinik*sinuk;
 
241
  vx = xmx*cosuk-cosnok*sinuk;
 
242
  vy = xmy*cosuk-sinnok*sinuk;
 
243
  vz = sinik*cosuk;
 
244
 
 
245
  /* Position and velocity */
 
246
  pos->x = rk*ux;
 
247
  pos->y = rk*uy;
 
248
  pos->z = rk*uz;
 
249
  vel->x = rdotk*ux+rfdotk*vx;
 
250
  vel->y = rdotk*uy+rfdotk*vy;
 
251
  vel->z = rdotk*uz+rfdotk*vz;
 
252
 
 
253
} /*SGP4*/
 
254
 
 
255
/*------------------------------------------------------------------*/
 
256
 
 
257
/* SDP4 */
 
258
/* This function is used to calculate the position and velocity */
 
259
/* of deep-space (period > 225 minutes) satellites. tsince is   */
 
260
/* time since epoch in minutes, tle is a pointer to a tle_t     */
 
261
/* structure with Keplerian orbital elements and pos and vel    */
 
262
/* are vector_t structures returning ECI satellite position and */
 
263
/* velocity. Use Convert_Sat_State() to convert to km and km/s. */
 
264
void 
 
265
SDP4(double tsince, tle_t *tle, vector_t *pos, vector_t *vel)
 
266
{
 
267
  int i;
 
268
 
 
269
  static double
 
270
    x3thm1,c1,x1mth2,c4,xnodcf,t2cof,xlcof,aycof,x7thm1;
 
271
 
 
272
  double
 
273
    a,axn,ayn,aynl,beta,betal,capu,cos2u,cosepw,cosik,
 
274
    cosnok,cosu,cosuk,ecose,elsq,epw,esine,pl,theta4,
 
275
    rdot,rdotk,rfdot,rfdotk,rk,sin2u,sinepw,sinik,
 
276
    sinnok,sinu,sinuk,tempe,templ,tsq,u,uk,ux,uy,uz,
 
277
    vx,vy,vz,xinck,xl,xlt,xmam,xmdf,xmx,xmy,xnoddf,
 
278
    xnodek,xll,a1,a3ovk2,ao,c2,coef,coef1,x1m5th,
 
279
    xhdot1,del1,r,delo,eeta,eta,etasq,perige,
 
280
    psisq,tsi,qoms24,s4,pinvsq,temp,tempa,temp1,
 
281
    temp2,temp3,temp4,temp5,temp6;
 
282
 
 
283
  static deep_arg_t deep_arg;
 
284
 
 
285
  /* Initialization */
 
286
  if (isFlagClear(SDP4_INITIALIZED_FLAG))
 
287
    {
 
288
      SetFlag(SDP4_INITIALIZED_FLAG);
 
289
 
 
290
      /* Recover original mean motion (xnodp) and   */
 
291
      /* semimajor axis (aodp) from input elements. */
 
292
      a1 = pow(xke/tle->xno,tothrd);
 
293
      deep_arg.cosio = cos(tle->xincl);
 
294
      deep_arg.theta2 = deep_arg.cosio*deep_arg.cosio;
 
295
      x3thm1 = 3*deep_arg.theta2-1;
 
296
      deep_arg.eosq = tle->eo*tle->eo;
 
297
      deep_arg.betao2 = 1-deep_arg.eosq;
 
298
      deep_arg.betao = sqrt(deep_arg.betao2);
 
299
      del1 = 1.5*ck2*x3thm1/(a1*a1*deep_arg.betao*deep_arg.betao2);
 
300
      ao = a1*(1-del1*(0.5*tothrd+del1*(1+134/81*del1)));
 
301
      delo = 1.5*ck2*x3thm1/(ao*ao*deep_arg.betao*deep_arg.betao2);
 
302
      deep_arg.xnodp = tle->xno/(1+delo);
 
303
      deep_arg.aodp = ao/(1-delo);
 
304
 
 
305
      /* For perigee below 156 km, the values */
 
306
      /* of s and qoms2t are altered.         */
 
307
      s4 = s;
 
308
      qoms24 = qoms2t;
 
309
      perige = (deep_arg.aodp*(1-tle->eo)-ae)*xkmper;
 
310
      if(perige < 156)
 
311
        {
 
312
          if(perige <= 98)
 
313
            s4 = 20;
 
314
          else
 
315
            s4 = perige-78;
 
316
          qoms24 = pow((120-s4)*ae/xkmper,4);
 
317
          s4 = s4/xkmper+ae;
 
318
        }
 
319
      pinvsq = 1/(deep_arg.aodp*deep_arg.aodp*
 
320
               deep_arg.betao2*deep_arg.betao2);
 
321
      deep_arg.sing = sin(tle->omegao);
 
322
      deep_arg.cosg = cos(tle->omegao);
 
323
      tsi = 1/(deep_arg.aodp-s4);
 
324
      eta = deep_arg.aodp*tle->eo*tsi;
 
325
      etasq = eta*eta;
 
326
      eeta = tle->eo*eta;
 
327
      psisq = fabs(1-etasq);
 
328
      coef = qoms24*pow(tsi,4);
 
329
      coef1 = coef/pow(psisq,3.5);
 
330
      c2 = coef1*deep_arg.xnodp*(deep_arg.aodp*(1+1.5*etasq+eeta*
 
331
           (4+etasq))+0.75*ck2*tsi/psisq*x3thm1*(8+3*etasq*(8+etasq)));
 
332
      c1 = tle->bstar*c2;
 
333
      deep_arg.sinio = sin(tle->xincl);
 
334
      a3ovk2 = -xj3/ck2*pow(ae,3);
 
335
      x1mth2 = 1-deep_arg.theta2;
 
336
      c4 = 2*deep_arg.xnodp*coef1*deep_arg.aodp*deep_arg.betao2*
 
337
           (eta*(2+0.5*etasq)+tle->eo*(0.5+2*etasq)-2*ck2*tsi/
 
338
           (deep_arg.aodp*psisq)*(-3*x3thm1*(1-2*eeta+etasq*
 
339
           (1.5-0.5*eeta))+0.75*x1mth2*(2*etasq-eeta*(1+etasq))*
 
340
           cos(2*tle->omegao)));
 
341
      theta4 = deep_arg.theta2*deep_arg.theta2;
 
342
      temp1 = 3*ck2*pinvsq*deep_arg.xnodp;
 
343
      temp2 = temp1*ck2*pinvsq;
 
344
      temp3 = 1.25*ck4*pinvsq*pinvsq*deep_arg.xnodp;
 
345
      deep_arg.xmdot = deep_arg.xnodp+0.5*temp1*deep_arg.betao*
 
346
                       x3thm1+0.0625*temp2*deep_arg.betao*
 
347
                       (13-78*deep_arg.theta2+137*theta4);
 
348
      x1m5th = 1-5*deep_arg.theta2;
 
349
      deep_arg.omgdot = -0.5*temp1*x1m5th+0.0625*temp2*
 
350
                        (7-114*deep_arg.theta2+395*theta4)+
 
351
                        temp3*(3-36*deep_arg.theta2+49*theta4);
 
352
      xhdot1 = -temp1*deep_arg.cosio;
 
353
      deep_arg.xnodot = xhdot1+(0.5*temp2*(4-19*deep_arg.theta2)+
 
354
                        2*temp3*(3-7*deep_arg.theta2))*deep_arg.cosio;
 
355
      xnodcf = 3.5*deep_arg.betao2*xhdot1*c1;
 
356
      t2cof = 1.5*c1;
 
357
      xlcof = 0.125*a3ovk2*deep_arg.sinio*(3+5*deep_arg.cosio)/
 
358
              (1+deep_arg.cosio);
 
359
      aycof = 0.25*a3ovk2*deep_arg.sinio;
 
360
      x7thm1 = 7*deep_arg.theta2-1;
 
361
 
 
362
      /* initialize Deep() */
 
363
      Deep(dpinit, tle, &deep_arg);
 
364
    }; /*End of SDP4() initialization */
 
365
 
 
366
  /* Update for secular gravity and atmospheric drag */
 
367
  xmdf = tle->xmo+deep_arg.xmdot*tsince;
 
368
  deep_arg.omgadf = tle->omegao+deep_arg.omgdot*tsince;
 
369
  xnoddf = tle->xnodeo+deep_arg.xnodot*tsince;
 
370
  tsq = tsince*tsince;
 
371
  deep_arg.xnode = xnoddf+xnodcf*tsq;
 
372
  tempa = 1-c1*tsince;
 
373
  tempe = tle->bstar*c4*tsince;
 
374
  templ = t2cof*tsq;
 
375
  deep_arg.xn = deep_arg.xnodp;
 
376
 
 
377
  /* Update for deep-space secular effects */
 
378
  deep_arg.xll = xmdf;
 
379
  deep_arg.t = tsince;
 
380
 
 
381
  Deep(dpsec, tle, &deep_arg);
 
382
 
 
383
  xmdf = deep_arg.xll;
 
384
  a = pow(xke/deep_arg.xn,tothrd)*tempa*tempa;
 
385
  deep_arg.em = deep_arg.em-tempe;
 
386
  xmam = xmdf+deep_arg.xnodp*templ;
 
387
 
 
388
  /* Update for deep-space periodic effects */
 
389
  deep_arg.xll = xmam;
 
390
 
 
391
  Deep(dpper, tle, &deep_arg);
 
392
 
 
393
  xmam = deep_arg.xll;
 
394
  xl = xmam+deep_arg.omgadf+deep_arg.xnode;
 
395
  beta = sqrt(1-deep_arg.em*deep_arg.em);
 
396
  deep_arg.xn = xke/pow(a,1.5);
 
397
 
 
398
  /* Long period periodics */
 
399
  axn = deep_arg.em*cos(deep_arg.omgadf);
 
400
  temp = 1/(a*beta*beta);
 
401
  xll = temp*xlcof*axn;
 
402
  aynl = temp*aycof;
 
403
  xlt = xl+xll;
 
404
  ayn = deep_arg.em*sin(deep_arg.omgadf)+aynl;
 
405
 
 
406
  /* Solve Kepler's Equation */
 
407
  capu = FMod2p(xlt-deep_arg.xnode);
 
408
  temp2 = capu;
 
409
 
 
410
  i = 0;
 
411
  do
 
412
    {
 
413
      sinepw = sin(temp2);
 
414
      cosepw = cos(temp2);
 
415
      temp3 = axn*sinepw;
 
416
      temp4 = ayn*cosepw;
 
417
      temp5 = axn*cosepw;
 
418
      temp6 = ayn*sinepw;
 
419
      epw = (capu-temp4+temp3-temp2)/(1-temp5-temp6)+temp2;
 
420
      if(fabs(epw-temp2) <= e6a) break;
 
421
      temp2 = epw;
 
422
    }
 
423
  while( i++ < 10 );
 
424
 
 
425
  /* Short period preliminary quantities */
 
426
  ecose = temp5+temp6;
 
427
  esine = temp3-temp4;
 
428
  elsq = axn*axn+ayn*ayn;
 
429
  temp = 1-elsq;
 
430
  pl = a*temp;
 
431
  r = a*(1-ecose);
 
432
  temp1 = 1/r;
 
433
  rdot = xke*sqrt(a)*esine*temp1;
 
434
  rfdot = xke*sqrt(pl)*temp1;
 
435
  temp2 = a*temp1;
 
436
  betal = sqrt(temp);
 
437
  temp3 = 1/(1+betal);
 
438
  cosu = temp2*(cosepw-axn+ayn*esine*temp3);
 
439
  sinu = temp2*(sinepw-ayn-axn*esine*temp3);
 
440
  u = AcTan(sinu,cosu);
 
441
  sin2u = 2*sinu*cosu;
 
442
  cos2u = 2*cosu*cosu-1;
 
443
  temp = 1/pl;
 
444
  temp1 = ck2*temp;
 
445
  temp2 = temp1*temp;
 
446
 
 
447
  /* Update for short periodics */
 
448
  rk = r*(1-1.5*temp2*betal*x3thm1)+0.5*temp1*x1mth2*cos2u;
 
449
  uk = u-0.25*temp2*x7thm1*sin2u;
 
450
  xnodek = deep_arg.xnode+1.5*temp2*deep_arg.cosio*sin2u;
 
451
  xinck = deep_arg.xinc+1.5*temp2*deep_arg.cosio*deep_arg.sinio*cos2u;
 
452
  rdotk = rdot-deep_arg.xn*temp1*x1mth2*sin2u;
 
453
  rfdotk = rfdot+deep_arg.xn*temp1*(x1mth2*cos2u+1.5*x3thm1);
 
454
 
 
455
  /* Orientation vectors */
 
456
  sinuk = sin(uk);
 
457
  cosuk = cos(uk);
 
458
  sinik = sin(xinck);
 
459
  cosik = cos(xinck);
 
460
  sinnok = sin(xnodek);
 
461
  cosnok = cos(xnodek);
 
462
  xmx = -sinnok*cosik;
 
463
  xmy = cosnok*cosik;
 
464
  ux = xmx*sinuk+cosnok*cosuk;
 
465
  uy = xmy*sinuk+sinnok*cosuk;
 
466
  uz = sinik*sinuk;
 
467
  vx = xmx*cosuk-cosnok*sinuk;
 
468
  vy = xmy*cosuk-sinnok*sinuk;
 
469
  vz = sinik*cosuk;
 
470
 
 
471
  /* Position and velocity */
 
472
  pos->x = rk*ux;
 
473
  pos->y = rk*uy;
 
474
  pos->z = rk*uz;
 
475
  vel->x = rdotk*ux+rfdotk*vx;
 
476
  vel->y = rdotk*uy+rfdotk*vy;
 
477
  vel->z = rdotk*uz+rfdotk*vz;
 
478
 
 
479
} /* SDP4 */
 
480
 
 
481
/*------------------------------------------------------------------*/
 
482
 
 
483
/* DEEP */
 
484
/* This function is used by SDP4 to add lunar and solar */
 
485
/* perturbation effects to deep-space orbit objects.    */
 
486
void
 
487
Deep(int ientry, tle_t *tle, deep_arg_t *deep_arg)
 
488
{
 
489
  static double
 
490
    thgr,xnq,xqncl,omegaq,zmol,zmos,savtsn,ee2,e3,xi2,
 
491
    xl2,xl3,xl4,xgh2,xgh3,xgh4,xh2,xh3,sse,ssi,ssg,xi3,
 
492
    se2,si2,sl2,sgh2,sh2,se3,si3,sl3,sgh3,sh3,sl4,sgh4,
 
493
    ssl,ssh,d3210,d3222,d4410,d4422,d5220,d5232,d5421,
 
494
    d5433,del1,del2,del3,fasx2,fasx4,fasx6,xlamo,xfact,
 
495
    xni,atime,stepp,stepn,step2,preep,pl,sghs,xli,
 
496
    d2201,d2211,sghl,sh1,pinc,pe,shs,zsingl,zcosgl,
 
497
    zsinhl,zcoshl,zsinil,zcosil;
 
498
 
 
499
  double
 
500
    a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,ainv2,alfdp,aqnv,
 
501
    sgh,sini2,sinis,sinok,sh,si,sil,day,betdp,dalf,
 
502
    bfact,c,cc,cosis,cosok,cosq,ctem,f322,zx,zy,
 
503
    dbet,dls,eoc,eq,f2,f220,f221,f3,f311,f321,xnoh,
 
504
    f330,f441,f442,f522,f523,f542,f543,g200,g201,
 
505
    g211,pgh,ph,s1,s2,s3,s4,s5,s6,s7,se,sel,ses,xls,
 
506
    g300,g310,g322,g410,g422,g520,g521,g532,g533,gam,
 
507
    sinq,sinzf,sis,sl,sll,sls,stem,temp,temp1,x1,x2,
 
508
    x2li,x2omi,x3,x4,x5,x6,x7,x8,xl,xldot,xmao,xnddt,
 
509
    xndot,xno2,xnodce,xnoi,xomi,xpidot,z1,z11,z12,z13,
 
510
    z2,z21,z22,z23,z3,z31,z32,z33,ze,zf,zm,zmo,zn,
 
511
    zsing,zsinh,zsini,zcosg,zcosh,zcosi,delt=0,ft=0;
 
512
 
 
513
  switch(ientry)
 
514
    {
 
515
    case dpinit : /* Entrance for deep space initialization */
 
516
      thgr = ThetaG(tle->epoch, deep_arg);
 
517
      eq = tle->eo;
 
518
      xnq = deep_arg->xnodp;
 
519
      aqnv = 1/deep_arg->aodp;
 
520
      xqncl = tle->xincl;
 
521
      xmao = tle->xmo;
 
522
      xpidot = deep_arg->omgdot+deep_arg->xnodot;
 
523
      sinq = sin(tle->xnodeo);
 
524
      cosq = cos(tle->xnodeo);
 
525
      omegaq = tle->omegao;
 
526
 
 
527
      /* Initialize lunar solar terms */
 
528
      day = deep_arg->ds50+18261.5;  /*Days since 1900 Jan 0.5*/
 
529
      if (day != preep)
 
530
        {
 
531
          preep = day;
 
532
          xnodce = 4.5236020-9.2422029E-4*day;
 
533
          stem = sin(xnodce);
 
534
          ctem = cos(xnodce);
 
535
          zcosil = 0.91375164-0.03568096*ctem;
 
536
          zsinil = sqrt(1-zcosil*zcosil);
 
537
          zsinhl = 0.089683511*stem/zsinil;
 
538
          zcoshl = sqrt(1-zsinhl*zsinhl);
 
539
          c = 4.7199672+0.22997150*day;
 
540
          gam = 5.8351514+0.0019443680*day;
 
541
          zmol = FMod2p(c-gam);
 
542
          zx = 0.39785416*stem/zsinil;
 
543
          zy = zcoshl*ctem+0.91744867*zsinhl*stem;
 
544
          zx = AcTan(zx,zy);
 
545
          zx = gam+zx-xnodce;
 
546
          zcosgl = cos(zx);
 
547
          zsingl = sin(zx);
 
548
          zmos = 6.2565837+0.017201977*day;
 
549
          zmos = FMod2p(zmos);
 
550
        } /* End if(day != preep) */
 
551
 
 
552
      /* Do solar terms */
 
553
      savtsn = 1E20;
 
554
      zcosg = zcosgs;
 
555
      zsing = zsings;
 
556
      zcosi = zcosis;
 
557
      zsini = zsinis;
 
558
      zcosh = cosq;
 
559
      zsinh = sinq;
 
560
      cc = c1ss;
 
561
      zn = zns;
 
562
      ze = zes;
 
563
      zmo = zmos;
 
564
      xnoi = 1/xnq;
 
565
 
 
566
      /* Loop breaks when Solar terms are done a second */
 
567
      /* time, after Lunar terms are initialized        */
 
568
      for(;;) 
 
569
        {
 
570
          /* Solar terms done again after Lunar terms are done */
 
571
          a1 = zcosg*zcosh+zsing*zcosi*zsinh;
 
572
          a3 = -zsing*zcosh+zcosg*zcosi*zsinh;
 
573
          a7 = -zcosg*zsinh+zsing*zcosi*zcosh;
 
574
          a8 = zsing*zsini;
 
575
          a9 = zsing*zsinh+zcosg*zcosi*zcosh;
 
576
          a10 = zcosg*zsini;
 
577
          a2 = deep_arg->cosio*a7+ deep_arg->sinio*a8;
 
578
          a4 = deep_arg->cosio*a9+ deep_arg->sinio*a10;
 
579
          a5 = -deep_arg->sinio*a7+ deep_arg->cosio*a8;
 
580
          a6 = -deep_arg->sinio*a9+ deep_arg->cosio*a10;
 
581
          x1 = a1*deep_arg->cosg+a2*deep_arg->sing;
 
582
          x2 = a3*deep_arg->cosg+a4*deep_arg->sing;
 
583
          x3 = -a1*deep_arg->sing+a2*deep_arg->cosg;
 
584
          x4 = -a3*deep_arg->sing+a4*deep_arg->cosg;
 
585
          x5 = a5*deep_arg->sing;
 
586
          x6 = a6*deep_arg->sing;
 
587
          x7 = a5*deep_arg->cosg;
 
588
          x8 = a6*deep_arg->cosg;
 
589
          z31 = 12*x1*x1-3*x3*x3;
 
590
          z32 = 24*x1*x2-6*x3*x4;
 
591
          z33 = 12*x2*x2-3*x4*x4;
 
592
          z1 = 3*(a1*a1+a2*a2)+z31*deep_arg->eosq;
 
593
          z2 = 6*(a1*a3+a2*a4)+z32*deep_arg->eosq;
 
594
          z3 = 3*(a3*a3+a4*a4)+z33*deep_arg->eosq;
 
595
          z11 = -6*a1*a5+deep_arg->eosq*(-24*x1*x7-6*x3*x5);
 
596
          z12 = -6*(a1*a6+a3*a5)+ deep_arg->eosq*
 
597
                (-24*(x2*x7+x1*x8)-6*(x3*x6+x4*x5));
 
598
          z13 = -6*a3*a6+deep_arg->eosq*(-24*x2*x8-6*x4*x6);
 
599
          z21 = 6*a2*a5+deep_arg->eosq*(24*x1*x5-6*x3*x7);
 
600
          z22 = 6*(a4*a5+a2*a6)+ deep_arg->eosq*
 
601
                (24*(x2*x5+x1*x6)-6*(x4*x7+x3*x8));
 
602
          z23 = 6*a4*a6+deep_arg->eosq*(24*x2*x6-6*x4*x8);
 
603
          z1 = z1+z1+deep_arg->betao2*z31;
 
604
          z2 = z2+z2+deep_arg->betao2*z32;
 
605
          z3 = z3+z3+deep_arg->betao2*z33;
 
606
          s3 = cc*xnoi;
 
607
          s2 = -0.5*s3/deep_arg->betao;
 
608
          s4 = s3*deep_arg->betao;
 
609
          s1 = -15*eq*s4;
 
610
          s5 = x1*x3+x2*x4;
 
611
          s6 = x2*x3+x1*x4;
 
612
          s7 = x2*x4-x1*x3;
 
613
          se = s1*zn*s5;
 
614
          si = s2*zn*(z11+z13);
 
615
          sl = -zn*s3*(z1+z3-14-6*deep_arg->eosq);
 
616
          sgh = s4*zn*(z31+z33-6);
 
617
          sh = -zn*s2*(z21+z23);
 
618
          if (xqncl < 5.2359877E-2) sh = 0;
 
619
          ee2 = 2*s1*s6;
 
620
          e3 = 2*s1*s7;
 
621
          xi2 = 2*s2*z12;
 
622
          xi3 = 2*s2*(z13-z11);
 
623
          xl2 = -2*s3*z2;
 
624
          xl3 = -2*s3*(z3-z1);
 
625
          xl4 = -2*s3*(-21-9*deep_arg->eosq)*ze;
 
626
          xgh2 = 2*s4*z32;
 
627
          xgh3 = 2*s4*(z33-z31);
 
628
          xgh4 = -18*s4*ze;
 
629
          xh2 = -2*s2*z22;
 
630
          xh3 = -2*s2*(z23-z21);
 
631
 
 
632
          if(isFlagSet(LUNAR_TERMS_DONE_FLAG)) break;
 
633
 
 
634
          /* Do lunar terms */
 
635
          sse = se;
 
636
          ssi = si;
 
637
          ssl = sl;
 
638
          ssh = sh/deep_arg->sinio;
 
639
          ssg = sgh-deep_arg->cosio*ssh;
 
640
          se2 = ee2;
 
641
          si2 = xi2;
 
642
          sl2 = xl2;
 
643
          sgh2 = xgh2;
 
644
          sh2 = xh2;
 
645
          se3 = e3;
 
646
          si3 = xi3;
 
647
          sl3 = xl3;
 
648
          sgh3 = xgh3;
 
649
          sh3 = xh3;
 
650
          sl4 = xl4;
 
651
          sgh4 = xgh4;
 
652
          zcosg = zcosgl;
 
653
          zsing = zsingl;
 
654
          zcosi = zcosil;
 
655
          zsini = zsinil;
 
656
          zcosh = zcoshl*cosq+zsinhl*sinq;
 
657
          zsinh = sinq*zcoshl-cosq*zsinhl;
 
658
          zn = znl;
 
659
          cc = c1l;
 
660
          ze = zel;
 
661
          zmo = zmol;
 
662
          SetFlag(LUNAR_TERMS_DONE_FLAG);
 
663
        } /* End of for(;;) */
 
664
 
 
665
      sse = sse+se;
 
666
      ssi = ssi+si;
 
667
      ssl = ssl+sl;
 
668
      ssg = ssg+sgh-deep_arg->cosio/deep_arg->sinio*sh;
 
669
      ssh = ssh+sh/deep_arg->sinio;
 
670
 
 
671
      /* Geopotential resonance initialization for 12 hour orbits */
 
672
      ClearFlag(RESONANCE_FLAG);
 
673
      ClearFlag(SYNCHRONOUS_FLAG);
 
674
 
 
675
      if( !((xnq < 0.0052359877) && (xnq > 0.0034906585)) )
 
676
        {
 
677
          if( (xnq < 0.00826) || (xnq > 0.00924) ) return;
 
678
          if (eq < 0.5) return;
 
679
          SetFlag(RESONANCE_FLAG);
 
680
          eoc = eq*deep_arg->eosq;
 
681
          g201 = -0.306-(eq-0.64)*0.440;
 
682
          if (eq <= 0.65)
 
683
            {
 
684
              g211 = 3.616-13.247*eq+16.290*deep_arg->eosq;
 
685
              g310 = -19.302+117.390*eq-228.419*
 
686
                      deep_arg->eosq+156.591*eoc;
 
687
              g322 = -18.9068+109.7927*eq-214.6334*
 
688
                     deep_arg->eosq+146.5816*eoc;
 
689
              g410 = -41.122+242.694*eq-471.094*
 
690
                     deep_arg->eosq+313.953*eoc;
 
691
              g422 = -146.407+841.880*eq-1629.014*
 
692
                     deep_arg->eosq+1083.435*eoc;
 
693
              g520 = -532.114+3017.977*eq-5740*
 
694
                     deep_arg->eosq+3708.276*eoc;
 
695
            }
 
696
          else
 
697
            {
 
698
              g211 = -72.099+331.819*eq-508.738*
 
699
                     deep_arg->eosq+266.724*eoc;
 
700
              g310 = -346.844+1582.851*eq-2415.925*
 
701
                     deep_arg->eosq+1246.113*eoc;
 
702
              g322 = -342.585+1554.908*eq-2366.899*
 
703
                     deep_arg->eosq+1215.972*eoc;
 
704
              g410 = -1052.797+4758.686*eq-7193.992*
 
705
                     deep_arg->eosq+3651.957*eoc;
 
706
              g422 = -3581.69+16178.11*eq-24462.77*
 
707
                     deep_arg->eosq+ 12422.52*eoc;
 
708
              if (eq <= 0.715)
 
709
                g520 = 1464.74-4664.75*eq+3763.64*deep_arg->eosq;
 
710
              else
 
711
                g520 = -5149.66+29936.92*eq-54087.36*
 
712
                       deep_arg->eosq+31324.56*eoc;
 
713
            } /* End if (eq <= 0.65) */
 
714
 
 
715
          if (eq < 0.7)
 
716
            {
 
717
              g533 = -919.2277+4988.61*eq-9064.77*
 
718
                     deep_arg->eosq+5542.21*eoc;
 
719
              g521 = -822.71072+4568.6173*eq-8491.4146*
 
720
                     deep_arg->eosq+5337.524*eoc;
 
721
              g532 = -853.666+4690.25*eq-8624.77*
 
722
                     deep_arg->eosq+ 5341.4*eoc;
 
723
            }
 
724
          else
 
725
            {
 
726
              g533 = -37995.78+161616.52*eq-229838.2*
 
727
                     deep_arg->eosq+109377.94*eoc;
 
728
              g521 = -51752.104+218913.95*eq-309468.16*
 
729
                     deep_arg->eosq+146349.42*eoc;
 
730
             g532 = -40023.88+170470.89*eq-242699.48*
 
731
                    deep_arg->eosq+115605.82*eoc;
 
732
            } /* End if (eq <= 0.7) */
 
733
 
 
734
          sini2 = deep_arg->sinio*deep_arg->sinio;
 
735
          f220 = 0.75*(1+2*deep_arg->cosio+deep_arg->theta2);
 
736
          f221 = 1.5*sini2;
 
737
          f321 = 1.875*deep_arg->sinio*(1-2*\
 
738
                 deep_arg->cosio-3*deep_arg->theta2);
 
739
          f322 = -1.875*deep_arg->sinio*(1+2*
 
740
                 deep_arg->cosio-3*deep_arg->theta2);
 
741
          f441 = 35*sini2*f220;
 
742
          f442 = 39.3750*sini2*sini2;
 
743
          f522 = 9.84375*deep_arg->sinio*(sini2*(1-2*deep_arg->cosio-5*
 
744
                 deep_arg->theta2)+0.33333333*(-2+4*deep_arg->cosio+
 
745
                 6*deep_arg->theta2));
 
746
          f523 = deep_arg->sinio*(4.92187512*sini2*(-2-4*
 
747
                 deep_arg->cosio+10*deep_arg->theta2)+6.56250012
 
748
                 *(1+2*deep_arg->cosio-3*deep_arg->theta2));
 
749
          f542 = 29.53125*deep_arg->sinio*(2-8*
 
750
                 deep_arg->cosio+deep_arg->theta2*
 
751
                 (-12+8*deep_arg->cosio+10*deep_arg->theta2));
 
752
          f543 = 29.53125*deep_arg->sinio*(-2-8*deep_arg->cosio+
 
753
                 deep_arg->theta2*(12+8*deep_arg->cosio-10*
 
754
                 deep_arg->theta2));
 
755
          xno2 = xnq*xnq;
 
756
          ainv2 = aqnv*aqnv;
 
757
          temp1 = 3*xno2*ainv2;
 
758
          temp = temp1*root22;
 
759
          d2201 = temp*f220*g201;
 
760
          d2211 = temp*f221*g211;
 
761
          temp1 = temp1*aqnv;
 
762
          temp = temp1*root32;
 
763
          d3210 = temp*f321*g310;
 
764
          d3222 = temp*f322*g322;
 
765
          temp1 = temp1*aqnv;
 
766
          temp = 2*temp1*root44;
 
767
          d4410 = temp*f441*g410;
 
768
          d4422 = temp*f442*g422;
 
769
          temp1 = temp1*aqnv;
 
770
          temp = temp1*root52;
 
771
          d5220 = temp*f522*g520;
 
772
          d5232 = temp*f523*g532;
 
773
          temp = 2*temp1*root54;
 
774
          d5421 = temp*f542*g521;
 
775
          d5433 = temp*f543*g533;
 
776
          xlamo = xmao+tle->xnodeo+tle->xnodeo-thgr-thgr;
 
777
          bfact = deep_arg->xmdot+deep_arg->xnodot+
 
778
                  deep_arg->xnodot-thdt-thdt;
 
779
          bfact = bfact+ssl+ssh+ssh;
 
780
        } /* if( !(xnq < 0.0052359877) && (xnq > 0.0034906585) ) */
 
781
      else
 
782
        {
 
783
          SetFlag(RESONANCE_FLAG);
 
784
          SetFlag(SYNCHRONOUS_FLAG);
 
785
          /* Synchronous resonance terms initialization */
 
786
          g200 = 1+deep_arg->eosq*(-2.5+0.8125*deep_arg->eosq);
 
787
          g310 = 1+2*deep_arg->eosq;
 
788
          g300 = 1+deep_arg->eosq*(-6+6.60937*deep_arg->eosq);
 
789
          f220 = 0.75*(1+deep_arg->cosio)*(1+deep_arg->cosio);
 
790
          f311 = 0.9375*deep_arg->sinio*deep_arg->sinio*
 
791
                 (1+3*deep_arg->cosio)-0.75*(1+deep_arg->cosio);
 
792
          f330 = 1+deep_arg->cosio;
 
793
          f330 = 1.875*f330*f330*f330;
 
794
          del1 = 3*xnq*xnq*aqnv*aqnv;
 
795
          del2 = 2*del1*f220*g200*q22;
 
796
          del3 = 3*del1*f330*g300*q33*aqnv;
 
797
          del1 = del1*f311*g310*q31*aqnv;
 
798
          fasx2 = 0.13130908;
 
799
          fasx4 = 2.8843198;
 
800
          fasx6 = 0.37448087;
 
801
          xlamo = xmao+tle->xnodeo+tle->omegao-thgr;
 
802
          bfact = deep_arg->xmdot+xpidot-thdt;
 
803
          bfact = bfact+ssl+ssg+ssh;
 
804
        } /* End if( !(xnq < 0.0052359877) && (xnq > 0.0034906585) ) */
 
805
 
 
806
      xfact = bfact-xnq;
 
807
 
 
808
      /* Initialize integrator */
 
809
      xli = xlamo;
 
810
      xni = xnq;
 
811
      atime = 0;
 
812
      stepp = 720;
 
813
      stepn = -720;
 
814
      step2 = 259200;
 
815
      /* End case dpinit: */
 
816
      return;
 
817
 
 
818
    case dpsec: /* Entrance for deep space secular effects */
 
819
      deep_arg->xll = deep_arg->xll+ssl*deep_arg->t;
 
820
      deep_arg->omgadf = deep_arg->omgadf+ssg*deep_arg->t;
 
821
      deep_arg->xnode = deep_arg->xnode+ssh*deep_arg->t;
 
822
      deep_arg->em = tle->eo+sse*deep_arg->t;
 
823
      deep_arg->xinc = tle->xincl+ssi*deep_arg->t;
 
824
      if (deep_arg->xinc < 0)
 
825
        {
 
826
          deep_arg->xinc = -deep_arg->xinc;
 
827
          deep_arg->xnode = deep_arg->xnode + pi;
 
828
          deep_arg->omgadf = deep_arg->omgadf-pi;
 
829
        }
 
830
      if( isFlagClear(RESONANCE_FLAG) ) return;
 
831
 
 
832
      do
 
833
        {
 
834
          if( (atime == 0) ||
 
835
             ((deep_arg->t >= 0) && (atime < 0)) || 
 
836
             ((deep_arg->t < 0) && (atime >= 0)) )
 
837
            {
 
838
              /* Epoch restart */
 
839
              if( deep_arg->t >= 0 )
 
840
                delt = stepp;
 
841
              else
 
842
                delt = stepn;
 
843
 
 
844
              atime = 0;
 
845
              xni = xnq;
 
846
              xli = xlamo;
 
847
            }
 
848
          else
 
849
            {     
 
850
              if( fabs(deep_arg->t) >= fabs(atime) )
 
851
                {
 
852
                  if ( deep_arg->t > 0 )
 
853
                    delt = stepp;
 
854
                  else
 
855
                    delt = stepn;
 
856
                }
 
857
            }
 
858
 
 
859
          do 
 
860
            {
 
861
              if ( fabs(deep_arg->t-atime) >= stepp )
 
862
                {
 
863
                  SetFlag(DO_LOOP_FLAG);
 
864
                  ClearFlag(EPOCH_RESTART_FLAG);
 
865
                }
 
866
              else
 
867
                {
 
868
                  ft = deep_arg->t-atime;
 
869
                  ClearFlag(DO_LOOP_FLAG);
 
870
                }
 
871
 
 
872
              if( fabs(deep_arg->t) < fabs(atime) )
 
873
                {
 
874
                  if (deep_arg->t >= 0)
 
875
                    delt = stepn;
 
876
                  else
 
877
                    delt = stepp;
 
878
                  SetFlag(DO_LOOP_FLAG | EPOCH_RESTART_FLAG);
 
879
                }
 
880
 
 
881
              /* Dot terms calculated */
 
882
              if( isFlagSet(SYNCHRONOUS_FLAG) )
 
883
                {
 
884
                  xndot = del1*sin(xli-fasx2)+del2*sin(2*(xli-fasx4))
 
885
                          +del3*sin(3*(xli-fasx6));
 
886
                  xnddt = del1*cos(xli-fasx2)+2*del2*cos(2*(xli-fasx4))
 
887
                          +3*del3*cos(3*(xli-fasx6));
 
888
                }
 
889
              else
 
890
                {
 
891
                  xomi = omegaq+deep_arg->omgdot*atime;
 
892
                  x2omi = xomi+xomi;
 
893
                  x2li = xli+xli;
 
894
                  xndot = d2201*sin(x2omi+xli-g22)
 
895
                          +d2211*sin(xli-g22)
 
896
                          +d3210*sin(xomi+xli-g32)
 
897
                          +d3222*sin(-xomi+xli-g32)
 
898
                          +d4410*sin(x2omi+x2li-g44)
 
899
                          +d4422*sin(x2li-g44)
 
900
                          +d5220*sin(xomi+xli-g52)
 
901
                          +d5232*sin(-xomi+xli-g52)
 
902
                          +d5421*sin(xomi+x2li-g54)
 
903
                          +d5433*sin(-xomi+x2li-g54);
 
904
                  xnddt = d2201*cos(x2omi+xli-g22)
 
905
                          +d2211*cos(xli-g22)
 
906
                          +d3210*cos(xomi+xli-g32)
 
907
                          +d3222*cos(-xomi+xli-g32)
 
908
                          +d5220*cos(xomi+xli-g52)
 
909
                          +d5232*cos(-xomi+xli-g52)
 
910
                          +2*(d4410*cos(x2omi+x2li-g44)
 
911
                          +d4422*cos(x2li-g44)
 
912
                          +d5421*cos(xomi+x2li-g54)
 
913
                          +d5433*cos(-xomi+x2li-g54));
 
914
                } /* End of if (isFlagSet(SYNCHRONOUS_FLAG)) */
 
915
 
 
916
              xldot = xni+xfact;
 
917
              xnddt = xnddt*xldot;
 
918
 
 
919
              if(isFlagSet(DO_LOOP_FLAG))
 
920
                {
 
921
                  xli = xli+xldot*delt+xndot*step2;
 
922
                  xni = xni+xndot*delt+xnddt*step2;
 
923
                  atime = atime+delt;
 
924
                }
 
925
            }
 
926
          while(isFlagSet(DO_LOOP_FLAG) && isFlagClear(EPOCH_RESTART_FLAG));
 
927
        }
 
928
      while(isFlagSet(DO_LOOP_FLAG) && isFlagSet(EPOCH_RESTART_FLAG));
 
929
 
 
930
      deep_arg->xn = xni+xndot*ft+xnddt*ft*ft*0.5;
 
931
      xl = xli+xldot*ft+xndot*ft*ft*0.5;
 
932
      temp = -deep_arg->xnode+thgr+deep_arg->t*thdt;
 
933
 
 
934
      if (isFlagClear(SYNCHRONOUS_FLAG))
 
935
        deep_arg->xll = xl+temp+temp;
 
936
      else
 
937
        deep_arg->xll = xl-deep_arg->omgadf+temp;
 
938
 
 
939
      return;
 
940
      /*End case dpsec: */
 
941
 
 
942
    case dpper: /* Entrance for lunar-solar periodics */
 
943
      sinis = sin(deep_arg->xinc);
 
944
      cosis = cos(deep_arg->xinc);
 
945
      if (fabs(savtsn-deep_arg->t) >= 30)
 
946
        {
 
947
          savtsn = deep_arg->t;
 
948
          zm = zmos+zns*deep_arg->t;
 
949
          zf = zm+2*zes*sin(zm);
 
950
          sinzf = sin(zf);
 
951
          f2 = 0.5*sinzf*sinzf-0.25;
 
952
          f3 = -0.5*sinzf*cos(zf);
 
953
          ses = se2*f2+se3*f3;
 
954
          sis = si2*f2+si3*f3;
 
955
          sls = sl2*f2+sl3*f3+sl4*sinzf;
 
956
          sghs = sgh2*f2+sgh3*f3+sgh4*sinzf;
 
957
          shs = sh2*f2+sh3*f3;
 
958
          zm = zmol+znl*deep_arg->t;
 
959
          zf = zm+2*zel*sin(zm);
 
960
          sinzf = sin(zf);
 
961
          f2 = 0.5*sinzf*sinzf-0.25;
 
962
          f3 = -0.5*sinzf*cos(zf);
 
963
          sel = ee2*f2+e3*f3;
 
964
          sil = xi2*f2+xi3*f3;
 
965
          sll = xl2*f2+xl3*f3+xl4*sinzf;
 
966
          sghl = xgh2*f2+xgh3*f3+xgh4*sinzf;
 
967
          sh1 = xh2*f2+xh3*f3;
 
968
          pe = ses+sel;
 
969
          pinc = sis+sil;
 
970
          pl = sls+sll;
 
971
        }
 
972
 
 
973
      pgh = sghs+sghl;
 
974
      ph = shs+sh1;
 
975
      deep_arg->xinc = deep_arg->xinc+pinc;
 
976
      deep_arg->em = deep_arg->em+pe;
 
977
 
 
978
      if (xqncl >= 0.2)
 
979
        {
 
980
          /* Apply periodics directly */
 
981
          ph = ph/deep_arg->sinio;
 
982
          pgh = pgh-deep_arg->cosio*ph;
 
983
          deep_arg->omgadf = deep_arg->omgadf+pgh;
 
984
          deep_arg->xnode = deep_arg->xnode+ph;
 
985
          deep_arg->xll = deep_arg->xll+pl;
 
986
        }
 
987
      else
 
988
        {
 
989
          /* Apply periodics with Lyddane modification */
 
990
          sinok = sin(deep_arg->xnode);
 
991
          cosok = cos(deep_arg->xnode);
 
992
          alfdp = sinis*sinok;
 
993
          betdp = sinis*cosok;
 
994
          dalf = ph*cosok+pinc*cosis*sinok;
 
995
          dbet = -ph*sinok+pinc*cosis*cosok;
 
996
          alfdp = alfdp+dalf;
 
997
          betdp = betdp+dbet;
 
998
          deep_arg->xnode = FMod2p(deep_arg->xnode);
 
999
          xls = deep_arg->xll+deep_arg->omgadf+cosis*deep_arg->xnode;
 
1000
          dls = pl+pgh-pinc*deep_arg->xnode*sinis;
 
1001
          xls = xls+dls;
 
1002
          xnoh = deep_arg->xnode;
 
1003
          deep_arg->xnode = AcTan(alfdp,betdp);
 
1004
 
 
1005
          /* This is a patch to Lyddane modification */
 
1006
          /* suggested by Rob Matson. */
 
1007
          if(fabs(xnoh-deep_arg->xnode) > pi)
 
1008
            {
 
1009
              if(deep_arg->xnode < xnoh)
 
1010
                deep_arg->xnode +=twopi;
 
1011
              else
 
1012
                deep_arg->xnode -=twopi;
 
1013
            }
 
1014
 
 
1015
          deep_arg->xll = deep_arg->xll+pl;
 
1016
          deep_arg->omgadf = xls-deep_arg->xll-cos(deep_arg->xinc)*
 
1017
                             deep_arg->xnode;
 
1018
        } /* End case dpper: */
 
1019
      return;
 
1020
 
 
1021
    } /* End switch(ientry) */
 
1022
 
 
1023
} /* End of Deep() */
 
1024
 
 
1025
/*------------------------------------------------------------------*/
 
1026
 
 
1027
/* Functions for testing and setting/clearing flags */
 
1028
 
 
1029
/* An int variable holding the single-bit flags */
 
1030
static int Flags = 0;
 
1031
 
 
1032
int
 
1033
isFlagSet(int flag)
 
1034
{
 
1035
  return (Flags & flag);
 
1036
}
 
1037
 
 
1038
int
 
1039
isFlagClear(int flag)
 
1040
{
 
1041
  return (~Flags & flag);
 
1042
}
 
1043
 
 
1044
void
 
1045
SetFlag(int flag)
 
1046
{
 
1047
  Flags |= flag;
 
1048
}
 
1049
 
 
1050
void
 
1051
ClearFlag(int flag)
 
1052
{
 
1053
  Flags &= ~flag;
 
1054
}
 
1055
 
 
1056
/*------------------------------------------------------------------*/