NASA Ocean Color

ocssw V2022
       SUBROUTINE der(T,X,DX)
 C  VERSION OF 4/18/28
 C  PURPOSE   
 C    COMPUTE DERIVATIVES OF THE EQUATIONS OF MOTION
 C  INPUT ARGUMENTS
 C    T      = CURRENT JULIAN TIME (SEC)
 C    X      = 6-D CURRENT CARTESIAN STATE, X,Y,Z,XD,YD,ZD, (KM,KM/SEC)
 C  INPUT COMMON BLOCKS
 C    OPTION, TIME, PLTCON, ATMCON, SPCCON, SUNCON, MUNCON, HARMON
 C  OUTPUT
 C    DX     = 6-D DERIVATIVES OF POSTION AND VELOCITY  
 C  CALL SUBROUTINES  
 C    LEGEND, ANGLES, EPHEM, XTHIRD, SETTHD, DENS76
 C  REFERENCES
 C    JPL EM 312/87-153, 20 APRIL 1987
 C  ANALYSIS  
 C    JOHNNY H. KWOK - JPL  
 C  PROGRAMMER
 C    JOHNNY H. KWOK - JPL  
 C  PROGRAM MODIFICATIONS 
 C    NONE  
 C  COMMENTS  
 C    THE HARMONICS ARE DIMENSIONED FOR 40X40, WHERE COEFS. ARE 
 C    STORED AS CIJ=C(I+1,J+1), ETC. FOR HIGHER FIELDS, CHANGE  
 C    DIMENSION ON P, CN, SN, TN, C, S. 
 C   
 C  MODIFICATION HISTORY
 C   Changed TIME common block name to TIMECMN, to eliminate linker 
 C   linker warnings.  B. A. Franz, GSC, November 14, 1997.
 C   
  
  
       IMPLICIT DOUBLE PRECISION (a-h,o-z)   
       dimension x(6),dx(6)  
       dimension p(41,41),cn(41),sn(41),tn(41)   
       dimension xs(3),xt(3),xm(3),xn(3),vb(3)
       common/option/l,m,ires,isun,imoon,iephem,idrag,idens,isrp,iorb
      1  ,iprint,inode,iplot
       common/timecmn/ti,tf,tr
       common/pltcon/ge,re,rate,pm,aj2,ellip,ratm
       common/atmcon/rdens,rht,sht,altmax,wt
       common/spccon/aread,areas,scmass,cdrag,csrp
       common/suncon/gs,es(7),et(7)
       common/muncon/gm,em(7),en(7)
       common/harmon/c(41,41),s(41,41)   
       DATA tpi/6.283185307179586d0/
       DATA zero,half,one,two,trhalf/0.d0,0.5d0,1.d0,2.d0,1.5d0/
       r2=x(1)**2+x(2)**2+x(3)**2
       r1=dsqrt(r2)  
       r3=r1*r2
       dx(1)=x(4)
       dx(2)=x(5)
       dx(3)=x(6)
       DO 10 i=4,6   
    10 dx(i)=zero
       IF (l.EQ.0.AND.m.EQ.0) GO TO 200
       IF (l.EQ.2.AND.m.EQ.0) GO TO 100
 C
 C  COMPUTE DERIVATIVE DUE TO SPHERICAL HARMONICS OF DEGREE L AND M
 C
       sr2=x(1)**2+x(2)**2   
       sr1=dsqrt(sr2)
       sphi=x(3)/r1  
       phi=dasin(sphi)   
       amda=datan2(x(2),x(1))-dmod(pm+rate*(t-tr),tpi)   
       im=m  
       IF (m.LT.l) im=m+1
       CALL legend(l,im,sphi,p)  
       CALL angles(m,amda,phi,cn,sn,tn)
       e1=zero
       e2=zero
       e3=zero
       d1=re/r1  
       d2=d1
       DO 110 il=2,l 
       f1=zero
       f2=zero
       f3=zero
       il1=il+1  
       DO 120 im=0,il,ires
       IF (im.GT.m) GO TO 130
       im2=im+2
       im1=im+1
       d3=c(il1,im1)*cn(im1)+s(il1,im1)*sn(im1)  
       f1=f1+d3*p(il1,im1)   
       IF (im2.LE.il1) f2=f2+d3*(p(il1,im2)-tn(im1)*p(il1,im1))  
       IF (im2.GT.il1) f2=f2-d3*tn(im1)*p(il1,im1)   
       IF (im.NE.0)  
      xf3=f3+im*(s(il1,im1)*cn(im1)-c(il1,im1)*sn(im1))*p(il1,im1)
   120 CONTINUE  
   130 CONTINUE  
       d2=d2*d1  
       e1=e1+d2*il1*f1   
       e2=e2+d2*f2
       e3=e3+d2*f3
   110 CONTINUE  
       d3=ge/r1
       e1=-e1*d3/r1  
       e2=e2*d3  
       e3=e3*d3  
       d1=e1/r1  
       d2=x(3)/(r2*sr1)*e2   
       d3=e3/sr2 
       dx(4)=dx(4)+(d1-d2)*x(1)-d3*x(2)  
       dx(5)=dx(5)+(d1-d2)*x(2)+d3*x(1)  
       dx(6)=dx(6)+d1*x(3)+sr1*e2/r2
       GO TO 200
 C
 C  COMPUTE DERIVATIVE DUE TO J2 ONLY.  IF ONLY J2 IS NEEDED, THIS IS A
 C  MUCH SIMPLER SET OF EQUATIONS TO USE.  NOTE THAT AJ2 IS -C20
 C
   100 CONTINUE
       r5=r2*r3
       d1=-trhalf*aj2*re*re*ge/r5
       d2=one-5.d0*x(3)**2/r2
       dx(4)=dx(4)+x(1)*d1*d2
       dx(5)=dx(5)+x(2)*d1*d2
       dx(6)=dx(6)+x(3)*d1*(d2+two)
   200 CONTINUE
 C
 C *** SETUP LUNI-SOLAR POSITION IF IEPHEM.EQ.1
 C
       IF (iephem.EQ.0) THEN
         trf=tr
       ELSE
         IF (isun.EQ.1.OR.imoon.EQ.1.OR.isrp.EQ.1)
      1  CALL ephem(isun,imoon,isrp,t,es,em)
         trf=t
       ENDIF
 C
 C  COMPUTE DERIVATIVE DUE TO THE SUN
 C
       IF (isun.EQ.0) GO TO 300  
       CALL xthird(t,trf,es,et,xs) 
       rs=dsqrt(xs(1)**2+xs(2)**2+xs(3)**2)
       rs3=rs**3
       DO 210 i=1,3
   210 xt(i)=x(i)-xs(i)
       rt=dsqrt(xt(1)**2+xt(2)**2+xt(3)**2)
       rt3=rt**3
       DO 220 i=1,3
   220 dx(i+3)=dx(i+3)-gs*(xt(i)/rt3+xs(i)/rs3)
   300 CONTINUE
 C
 C  COMPUTE DERIVATIVE DUE TO SOLAR RADIATION PRESSURE
 C
       IF (isrp.EQ.0) GO TO 400
       IF (isun.EQ.0) THEN
         CALL xthird(t,trf,es,et,xs)
         rs=dsqrt(xs(1)**2+xs(2)**2+xs(3)**2)
       ENDIF
       DO 310 i=1,3
   310 xs(i)=xs(i)/rs
 C
 C  CHECK FOR SHADOW CONDITION
 C
       rdrs=x(1)*xs(1)+x(2)*xs(2)+x(3)*xs(3)
       IF (rdrs.GE.zero) GO TO 320
       theta=dacos(rdrs/r1)
       IF (r1*dsin(theta).LT.ratm) GO TO 400
   320 CONTINUE
       DO 330 i=1,3
   330 dx(i+3)=dx(i+3)-csrp*areas*xs(i)/scmass
 C
 C  COMPUTE DERIVATIVE DUE TO THE MOON
 C
   400 CONTINUE
       IF (imoon.EQ.0) GO TO 500
       IF (iephem.EQ.1) CALL setthd(em,en)
       CALL xthird(t,trf,em,en,xm) 
       rm=(xm(1)**2+xm(2)**2+xm(3)**2)**trhalf
       DO 410 i=1,3
   410 xn(i)=x(i)-xm(i)  
       rn=(xn(1)**2+xn(2)**2+xn(3)**2)**trhalf
       DO 420 i=1,3  
   420 dx(i+3)=dx(i+3)-gm*(xn(i)/rn+xm(i)/rm)
 C
 C  COMPUTE DERIVATIVE DUE TO DRAG
 C 
   500 CONTINUE
       IF (idrag.EQ.0) GO TO 600
       IF (ellip.EQ.zero) THEN
         ealt=r1-re
       ELSE
         phi=dasin(x(3)/r1)
         ealt=r1-dsqrt(re**2*(one-ellip**2)/(one-(ellip*dcos(phi))**2))
       ENDIF
       IF (ealt.GT.altmax) GO TO 600
 C
 C  COMPUTE DENSITY AT ALTITUDE.  THIS PORTION CAN BE REPLACED BY MORE
 C  SOPHISCATED MODEL
 C
       IF (idens.EQ.0) dens=rdens*dexp((rht-ealt)/sht)
       IF (idens.EQ.1) CALL dens76(ealt,dens)
       dens=dens*wt
 C
       vb(1)=x(4)+x(2)*rate
       vb(2)=x(5)-x(1)*rate
       vb(3)=x(6)
       vbm2=vb(1)**2+vb(2)**2+vb(3)**2
       vbm1=dsqrt(vbm2)
       DO 510 i=1,3
   510 vb(i)=vb(i)/vbm1
       DO 520 i=1,3
   520 dx(i+3)=dx(i+3)-half*cdrag*aread/scmass*dens*vbm2*vb(i)
 C
 C  ADD THE TWO BODY DERIVATIVE TOO
 C
   600 CONTINUE
       DO 610 i=1,3  
   610 dx(i+3)=dx(i+3)-ge*x(i)/r3
       RETURN
       END   