der.f

Go to the documentation of this file.

      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