OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
der.f
Go to the documentation of this file.
1  SUBROUTINE der(T,X,DX)
2 C VERSION OF 4/18/28
3 C PURPOSE
4 C COMPUTE DERIVATIVES OF THE EQUATIONS OF MOTION
5 C INPUT ARGUMENTS
6 C T = CURRENT JULIAN TIME (SEC)
7 C X = 6-D CURRENT CARTESIAN STATE, X,Y,Z,XD,YD,ZD, (KM,KM/SEC)
8 C INPUT COMMON BLOCKS
9 C OPTION, TIME, PLTCON, ATMCON, SPCCON, SUNCON, MUNCON, HARMON
10 C OUTPUT
11 C DX = 6-D DERIVATIVES OF POSTION AND VELOCITY
12 C CALL SUBROUTINES
13 C LEGEND, ANGLES, EPHEM, XTHIRD, SETTHD, DENS76
14 C REFERENCES
15 C JPL EM 312/87-153, 20 APRIL 1987
16 C ANALYSIS
17 C JOHNNY H. KWOK - JPL
18 C PROGRAMMER
19 C JOHNNY H. KWOK - JPL
20 C PROGRAM MODIFICATIONS
21 C NONE
22 C COMMENTS
23 C THE HARMONICS ARE DIMENSIONED FOR 40X40, WHERE COEFS. ARE
24 C STORED AS CIJ=C(I+1,J+1), ETC. FOR HIGHER FIELDS, CHANGE
25 C DIMENSION ON P, CN, SN, TN, C, S.
26 C
27 C MODIFICATION HISTORY
28 C Changed TIME common block name to TIMECMN, to eliminate linker
29 C linker warnings. B. A. Franz, GSC, November 14, 1997.
30 C
31 
32 
33  IMPLICIT DOUBLE PRECISION (a-h,o-z)
34  dimension x(6),dx(6)
35  dimension p(41,41),cn(41),sn(41),tn(41)
36  dimension xs(3),xt(3),xm(3),xn(3),vb(3)
37  common/option/l,m,ires,isun,imoon,iephem,idrag,idens,isrp,iorb
38  1 ,iprint,inode,iplot
39  common/timecmn/ti,tf,tr
40  common/pltcon/ge,re,rate,pm,aj2,ellip,ratm
41  common/atmcon/rdens,rht,sht,altmax,wt
42  common/spccon/aread,areas,scmass,cdrag,csrp
43  common/suncon/gs,es(7),et(7)
44  common/muncon/gm,em(7),en(7)
45  common/harmon/c(41,41),s(41,41)
46  DATA tpi/6.283185307179586d0/
47  DATA zero,half,one,two,trhalf/0.d0,0.5d0,1.d0,2.d0,1.5d0/
48  r2=x(1)**2+x(2)**2+x(3)**2
49  r1=dsqrt(r2)
50  r3=r1*r2
51  dx(1)=x(4)
52  dx(2)=x(5)
53  dx(3)=x(6)
54  DO 10 i=4,6
55  10 dx(i)=zero
56  IF (l.EQ.0.AND.m.EQ.0) GO TO 200
57  IF (l.EQ.2.AND.m.EQ.0) GO TO 100
58 C
59 C COMPUTE DERIVATIVE DUE TO SPHERICAL HARMONICS OF DEGREE L AND M
60 C
61  sr2=x(1)**2+x(2)**2
62  sr1=dsqrt(sr2)
63  sphi=x(3)/r1
64  phi=dasin(sphi)
65  amda=datan2(x(2),x(1))-dmod(pm+rate*(t-tr),tpi)
66  im=m
67  IF (m.LT.l) im=m+1
68  CALL legend(l,im,sphi,p)
69  CALL angles(m,amda,phi,cn,sn,tn)
70  e1=zero
71  e2=zero
72  e3=zero
73  d1=re/r1
74  d2=d1
75  DO 110 il=2,l
76  f1=zero
77  f2=zero
78  f3=zero
79  il1=il+1
80  DO 120 im=0,il,ires
81  IF (im.GT.m) GO TO 130
82  im2=im+2
83  im1=im+1
84  d3=c(il1,im1)*cn(im1)+s(il1,im1)*sn(im1)
85  f1=f1+d3*p(il1,im1)
86  IF (im2.LE.il1) f2=f2+d3*(p(il1,im2)-tn(im1)*p(il1,im1))
87  IF (im2.GT.il1) f2=f2-d3*tn(im1)*p(il1,im1)
88  IF (im.NE.0)
89  xf3=f3+im*(s(il1,im1)*cn(im1)-c(il1,im1)*sn(im1))*p(il1,im1)
90  120 CONTINUE
91  130 CONTINUE
92  d2=d2*d1
93  e1=e1+d2*il1*f1
94  e2=e2+d2*f2
95  e3=e3+d2*f3
96  110 CONTINUE
97  d3=ge/r1
98  e1=-e1*d3/r1
99  e2=e2*d3
100  e3=e3*d3
101  d1=e1/r1
102  d2=x(3)/(r2*sr1)*e2
103  d3=e3/sr2
104  dx(4)=dx(4)+(d1-d2)*x(1)-d3*x(2)
105  dx(5)=dx(5)+(d1-d2)*x(2)+d3*x(1)
106  dx(6)=dx(6)+d1*x(3)+sr1*e2/r2
107  GO TO 200
108 C
109 C COMPUTE DERIVATIVE DUE TO J2 ONLY. IF ONLY J2 IS NEEDED, THIS IS A
110 C MUCH SIMPLER SET OF EQUATIONS TO USE. NOTE THAT AJ2 IS -C20
111 C
112  100 CONTINUE
113  r5=r2*r3
114  d1=-trhalf*aj2*re*re*ge/r5
115  d2=one-5.d0*x(3)**2/r2
116  dx(4)=dx(4)+x(1)*d1*d2
117  dx(5)=dx(5)+x(2)*d1*d2
118  dx(6)=dx(6)+x(3)*d1*(d2+two)
119  200 CONTINUE
120 C
121 C *** SETUP LUNI-SOLAR POSITION IF IEPHEM.EQ.1
122 C
123  IF (iephem.EQ.0) THEN
124  trf=tr
125  ELSE
126  IF (isun.EQ.1.OR.imoon.EQ.1.OR.isrp.EQ.1)
127  1 CALL ephem(isun,imoon,isrp,t,es,em)
128  trf=t
129  ENDIF
130 C
131 C COMPUTE DERIVATIVE DUE TO THE SUN
132 C
133  IF (isun.EQ.0) GO TO 300
134  CALL xthird(t,trf,es,et,xs)
135  rs=dsqrt(xs(1)**2+xs(2)**2+xs(3)**2)
136  rs3=rs**3
137  DO 210 i=1,3
138  210 xt(i)=x(i)-xs(i)
139  rt=dsqrt(xt(1)**2+xt(2)**2+xt(3)**2)
140  rt3=rt**3
141  DO 220 i=1,3
142  220 dx(i+3)=dx(i+3)-gs*(xt(i)/rt3+xs(i)/rs3)
143  300 CONTINUE
144 C
145 C COMPUTE DERIVATIVE DUE TO SOLAR RADIATION PRESSURE
146 C
147  IF (isrp.EQ.0) GO TO 400
148  IF (isun.EQ.0) THEN
149  CALL xthird(t,trf,es,et,xs)
150  rs=dsqrt(xs(1)**2+xs(2)**2+xs(3)**2)
151  ENDIF
152  DO 310 i=1,3
153  310 xs(i)=xs(i)/rs
154 C
155 C CHECK FOR SHADOW CONDITION
156 C
157  rdrs=x(1)*xs(1)+x(2)*xs(2)+x(3)*xs(3)
158  IF (rdrs.GE.zero) GO TO 320
159  theta=dacos(rdrs/r1)
160  IF (r1*dsin(theta).LT.ratm) GO TO 400
161  320 CONTINUE
162  DO 330 i=1,3
163  330 dx(i+3)=dx(i+3)-csrp*areas*xs(i)/scmass
164 C
165 C COMPUTE DERIVATIVE DUE TO THE MOON
166 C
167  400 CONTINUE
168  IF (imoon.EQ.0) GO TO 500
169  IF (iephem.EQ.1) CALL setthd(em,en)
170  CALL xthird(t,trf,em,en,xm)
171  rm=(xm(1)**2+xm(2)**2+xm(3)**2)**trhalf
172  DO 410 i=1,3
173  410 xn(i)=x(i)-xm(i)
174  rn=(xn(1)**2+xn(2)**2+xn(3)**2)**trhalf
175  DO 420 i=1,3
176  420 dx(i+3)=dx(i+3)-gm*(xn(i)/rn+xm(i)/rm)
177 C
178 C COMPUTE DERIVATIVE DUE TO DRAG
179 C
180  500 CONTINUE
181  IF (idrag.EQ.0) GO TO 600
182  IF (ellip.EQ.zero) THEN
183  ealt=r1-re
184  ELSE
185  phi=dasin(x(3)/r1)
186  ealt=r1-dsqrt(re**2*(one-ellip**2)/(one-(ellip*dcos(phi))**2))
187  ENDIF
188  IF (ealt.GT.altmax) GO TO 600
189 C
190 C COMPUTE DENSITY AT ALTITUDE. THIS PORTION CAN BE REPLACED BY MORE
191 C SOPHISCATED MODEL
192 C
193  IF (idens.EQ.0) dens=rdens*dexp((rht-ealt)/sht)
194  IF (idens.EQ.1) CALL dens76(ealt,dens)
195  dens=dens*wt
196 C
197  vb(1)=x(4)+x(2)*rate
198  vb(2)=x(5)-x(1)*rate
199  vb(3)=x(6)
200  vbm2=vb(1)**2+vb(2)**2+vb(3)**2
201  vbm1=dsqrt(vbm2)
202  DO 510 i=1,3
203  510 vb(i)=vb(i)/vbm1
204  DO 520 i=1,3
205  520 dx(i+3)=dx(i+3)-half*cdrag*aread/scmass*dens*vbm2*vb(i)
206 C
207 C ADD THE TWO BODY DERIVATIVE TOO
208 C
209  600 CONTINUE
210  DO 610 i=1,3
211  610 dx(i+3)=dx(i+3)-ge*x(i)/r3
212  RETURN
213  END
float f1(float x)
subroutine ephem(ISUN, IMOON, ISRP, T, ES, EM)
Definition: ephem.f:2
float f3(float z)
float f2(float y)
#define re
Definition: l1_czcs_hdf.c:701
subroutine legend(N, M, X, P)
Definition: legend.f:2
subroutine xthird(T, TR, ES, ET, XS)
Definition: xthird.f:2
subroutine dens76(H, DENS)
Definition: dens76.f:2
subroutine der(T, X, DX)
Definition: der.f:2
subroutine setthd(X, Y)
Definition: setthd.f:2