subroutine galsat(r,rorb,tjd,ksat,kflag) c**************************************************************** c this version calculates jupiter-centered (ksat=1,4) coordinates c or barycentric (ksat=-1,-4) satellite coordinates, c or coordinates of jupiter relative to barycenter (ksat=0) c**************************************************************** c documented in "theory of motion of jupiter's galilean satellites," c astronomy & astrophysics 56, 333-352 (1977). c see also jpl engineering memorandum 314-112 (1 mar 1977). c c c input c tjd (double precision) = tdb julian date requested c nsat (integer) = satellite i.d. c nsat = 1 is io c nsat = 2 is europa c nsat = 3 is ganymede c nsat = 4 is callisto c c kflag=1 means position only c kflag=2 means position and velocity c c output c r (dimension 6) = position and velocity in earth j2000 equatorial c frame in au, au/day implicit none integer ksat, kflag double precision r(6),tjd double precision rorb(6) logical first save first character*80 verson integer ioin, ioout, iomsg, iopun common /ioblok/ ioin, ioout, iomsg, iopun double precision earay,baray,paray common /ebblok/ earay(28),baray(22),paray(28) save /ebblok/ double precision trmcod common /trmblk/ trmcod(53) save /trmblk/ double precision angcod,ratcod,ang,rat common /angblk/ angcod(99),ratcod(99),ang(22),rat(23) save /angblk/ double precision axis(4),cxi1(10),argx1(10),ratx1(10), 1 cz1(7),argz1(7),ratz1(7),cv1(41),argv1(41),ratv1(41),cxi2(24), 2 argx2(24),ratx2(24),cz2(11),argz2(11),ratz2(11),cv2(66), 3 argv2(66),ratv2(66),cxi3(31),argx3(31),ratx3(31),cz3(13), 4 argz3(13),ratz3(13),cv3(75),argv3(75),ratv3(75),cxi4(49), 5 argx4(49),ratx4(49),cz4(18),argz4(18),ratz4(18), 6 cv4(89),argv4(89),ratv4(89),epsln integer nxi1t,nz1t,nv1t,nxi2t,nz2t,nv2t,nxi3t,nz3t, 2 nv3t,nxi4t,nz4t,nv4t integer kodx1,kodz1,kodv1, 2 kodx2,kodz2,kodv2, 3 kodx3,kodz3,kodv3, 4 kodx4,kodz4,kodv4 common /theory/ axis,cxi1,argx1,ratx1,cz1,argz1,ratz1,cv1, 1 argv1,ratv1,cxi2,argx2,ratx2,cz2,argz2,ratz2,cv2,argv2,ratv2, 2 cxi3,argx3,ratx3,cz3,argz3,ratz3,cv3,argv3,ratv3,cxi4,argx4, 3 ratx4,cz4,argz4,ratz4,cv4,argv4,ratv4,epsln,nxi1t,nz1t,nv1t, 4 nxi2t,nz2t,nv2t,nxi3t,nz3t,nv3t,nxi4t,nz4t,nv4t, 5 kodx1(2,10),kodz1(2,7),kodv1(2,41),kodx2(2,24), 6 kodz2(2,11),kodv2(2,66),kodx3(2,31),kodz3(2,13), 7 kodv3(2,75),kodx4(2,49),kodz4(2,18),kodv4(2,89) save /theory/ c--set up equivalence arrays for ease in calling double precision satx1(10,3),satz1(7,3),satv1(41,3), 1 satx2(24,3),satz2(11,3),satv2(66,3),satx3(31,3),satz3(13,3), 2 satv3(75,3),satx4(49,3),satz4(18,3),satv4(89,3) equivalence (satx1(1,1),cxi1), (satz1(1,1),cz1), 1 (satv1(1,1),cv1), (satx2(1,1),cxi2), (satz2(1,1),cz2), 2 (satv2(1,1),cv2), (satx3(1,1),cxi3), (satz3(1,1),cz3), 3 (satv3(1,1),cv3), (satx4(1,1),cxi4), (satz4(1,1),cz4), 4 (satv4(1,1),cv4) double precision tref double precision qpsi11,qpsi21 c for galsat/galsap local communication double precision twopi,degrad,t,rb(6) c these variables are for calculation of bary to jup vector c they make use of fact that cofby=cofbx with cos replaced by sin double precision cofbx(7),cofbz(5) c note that cofbx and cofbz were single precision in the @for version! double precision angbx(7,2),angbz(5,2) c--local common block for subroutines common /local/ twopi,degrad,t,rb,angbx,angbz,cofbx,cofbz save /local/ c c note: the above /local/ is the definition of the common block. c to strictly follow f77 standards, all occurrences of /local/ c must be declared the same length. i let the compiler do it. c affected routines are samjay, barcor in galsat c and samjap, dekod, barcop in galsap c and updat, chkgal, revizg in satsap c local variables to be saved for galsat/galsap and chkgal, qqdot double precision tlast,tlastg,cj,sj,ci,si,cn,sn,cp,sp,phi,phidot, $ p(3,3),q(3,3),qdot(3,3),qmat(3,3), tlastp common /svtloc/ tlast,tlastg,cj,sj,ci,si,cn,sn,cp,sp,phi,phidot, $ p,q,qdot,qmat, tlastp save /svtloc/ integer i, j, kflg, nflag, nsat c tolt is tolerance for tlast updates, tolg is tolerance for tlastg updates double precision tolt double precision tolg data tolt /1.d-4/ data tolg /50.d0/ c data first/.true./ data tref/2443000.5d0/ data verson/' galsat version 5.0 jay lieske, jpl'/ c..start computation nflag=kflag kflg=1 if(ksat.le.0) kflg=-1 nsat=iabs(ksat) if((nsat.lt.0).or.(nsat.gt.4)) then write(iomsg,100) 100 format(1h1,'error exit in galsat. run terminated.') if (iomsg .ne. ioout) write(ioout,100) stop endif t=tjd-tref c note that chkgal will re-set values for tlast, tlastg, tlastp. c galsap will check for previous use by galsat, but galsat won't check c for galsap. c hence, for rigorous check-out one should first call galsat before galsap or c else always update the ancillary angles by setting tolt and tolg equal to c zero. if (first) then first=.false. call chkgal(verson) endif if(dabs(t-tlast) .gt. tolt ) then tlast=t c mark tlastp as being 'dirty' so can intermix galsat/galsap calls c because galsap uses 'q' in a different sense for partials c first call to galsap will initialize tlastp to -4.d20, so this c is our way of telling whether or not galsap has been initialized, c as well as being able to force a tlastg update tlastp = -6.d20 call qqdot(t,qpsi11,qpsi21) call mult3g(p,qdot,qmat) call mult3g(p,q,qdot) endif c now pq is in qdot (for position) and pqdot is in qmat (for vel) c positions are r = qdot * rb, c velocities are rdot = qdot * rbdot + qmat * rb c--revise g by adding dg now c--these update arguments for dg every 50 days if(dabs(t-tlastg) .gt. tolg ) then tlastg=t call revizg endif c now compute satellite coordinates c since do not have args subscripted, set up subr do 52 i=1,6 rb(i)=0.d0 52 r(i)=0.d0 if(nsat.eq.1) then call samjay(nxi1t,nv1t,nz1t,satx1,satv1,satz1,10,41 $ ,7,nsat,ang,rat,nflag) else if(nsat.eq.2) then call samjay(nxi2t,nv2t,nz2t,satx2,satv2,satz2,24,66 $ ,11,nsat,ang,rat,nflag) else if(nsat.eq.3) then call samjay(nxi3t,nv3t,nz3t,satx3,satv3,satz3,31,75 $ ,13,nsat,ang,rat,nflag) else if(nsat.eq.4) then call samjay(nxi4t,nv4t,nz4t,satx4,satv4,satz4,49,89 $ ,18,nsat,ang,rat,nflag) endif if(kflg.lt.0) call barcor c..return the orbital plane state vector do i=1,6 rorb(i) = rb(i) enddo c--now go to earth's mean equator 53 do 5 i=1,3 do 5 j=1,3 r(i)=r(i)+qdot(i,j)*rb(j) if(kflag.ne.1) $ r(i+3)=r(i+3)+qdot(i,j)*rb(j+3)+qmat(i,j)*rb(j) 5 continue return end C = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine samjay(nx,nv,nz,cx,cv,cz,nmx,nmv,nmz,nsat,ang,rat, $ nflag) c**************************************************** implicit none save integer nx, nv, nz, nmx, nmv, nmz, nsat, nflag double precision ang(*),rat(*) double precision xi,v,s,xidot,vdot,sdot,angl,dt,ca,sa double precision cx(nmx,3),cv(nmv,3),cz(nmz,3) double precision q1,q2,q3,q4,sdfac c>> note: the q1..q4 variables are employed for consistency with samjap in galsap c>> if they are not employed, the code would be somewhat shorter. c>> they are used here to maintain compatibility and ease of program revisions. c>> for galsat/galsap local communication *j* double precision twopi,degrad,t,rb(6) *j*c--local common block for subroutines *j* common /local/ twopi,degrad,t,rb c>> c> only twopi,degrad,t,rb(6) are needed here, but we include all to meet f77 standard c****************************************** c--local common block for subroutines double precision twopi,degrad,t,rb(6) double precision cofbx(7),cofbz(5) double precision angbx(7,2),angbz(5,2) common /local/ twopi,degrad,t,rb,angbx,angbz,cofbx,cofbz c>> c****************************************** c>> c****************************************** *j* double precision axis(4) *j* common /theory/ axis c> only axis is needed here, but we include all to meet f77 standard. c****************************************** double precision axis(4),cxi1(10),argx1(10),ratx1(10), 1 cz1(7),argz1(7),ratz1(7),cv1(41),argv1(41),ratv1(41),cxi2(24), 2 argx2(24),ratx2(24),cz2(11),argz2(11),ratz2(11),cv2(66), 3 argv2(66),ratv2(66),cxi3(31),argx3(31),ratx3(31),cz3(13), 4 argz3(13),ratz3(13),cv3(75),argv3(75),ratv3(75),cxi4(49), 5 argx4(49),ratx4(49),cz4(18),argz4(18),ratz4(18), 6 cv4(89),argv4(89),ratv4(89),epsln integer nxi1t,nz1t,nv1t,nxi2t,nz2t,nv2t,nxi3t,nz3t, 2 nv3t,nxi4t,nz4t,nv4t integer kodx1,kodz1,kodv1, 2 kodx2,kodz2,kodv2, 3 kodx3,kodz3,kodv3, 4 kodx4,kodz4,kodv4 common /theory/ axis,cxi1,argx1,ratx1,cz1,argz1,ratz1,cv1, 1 argv1,ratv1,cxi2,argx2,ratx2,cz2,argz2,ratz2,cv2,argv2,ratv2, 2 cxi3,argx3,ratx3,cz3,argz3,ratz3,cv3,argv3,ratv3,cxi4,argx4, 3 ratx4,cz4,argz4,ratz4,cv4,argv4,ratv4,epsln,nxi1t,nz1t,nv1t, 4 nxi2t,nz2t,nv2t,nxi3t,nz3t,nv3t,nxi4t,nz4t,nv4t, 5 kodx1(2,10),kodz1(2,7),kodv1(2,41),kodx2(2,24), 6 kodz2(2,11),kodv2(2,66),kodx3(2,31),kodz3(2,13), 7 kodv3(2,75),kodx4(2,49),kodz4(2,18),kodv4(2,89) c****************************************** c> integer k c-- xi=0.d0 v=0.d0 s=0.d0 xidot=0.d0 vdot=0.d0 sdot=0.d0 do 1 k=1,nx angl=dmod(cx(k,2)+cx(k,3)*t,twopi) ca=dcos(angl) q1=cx(k,1)*ca xi=xi+q1 if(nflag.ne.1) then sa=dsin(angl) q2=cx(k,1)*sa xidot=xidot-q2*cx(k,3) endif 1 continue do 2 k=1,nv angl=dmod(cv(k,2)+cv(k,3)*t,twopi) sa=dsin(angl) q2=cv(k,1)*sa v=v+q2 if(nflag.ne.1) then ca=dcos(angl) q1=cv(k,1)*ca vdot=vdot+q1*cv(k,3) endif 2 continue dt=v/rat(nsat) sdfac=1.d0+vdot/rat(nsat) c>> sdfac is irrelevant (its value will be 1) when nflag=1 (position-only) do 3 k=1,nz angl=dmod(cz(k,2)+cz(k,3)*(t+dt),twopi) sa=dsin(angl) q2=cz(k,1)*sa s=s+q2 if(nflag.ne.1) then ca=dcos(angl) q1=cz(k,1)*ca *jay put factor later !sdot=sdot+q1*cz(k,3)*(1.d0+vdot/rat(nsat)) sdot=sdot+q1*cz(k,3) endif 3 continue c--this is l-psi+v angl=dmod(ang(nsat)-ang(15)+(rat(nsat)-rat(15))*t,twopi)+v q1=axis(nsat)*dcos(angl) q2=axis(nsat)*dsin(angl) q3=axis(nsat)*s q4=1.d0+xi rb(1)=q1*q4 rb(2)=q2*q4 rb(3)=q3*q4 if(nflag.eq.1) return c>> now correct for the sdot factor in time-completed: sdot=sdot*sdfac ca=rat(nsat)-rat(15)+vdot rb(4)=q1*xidot-rb(2)*ca rb(5)=q2*xidot+rb(1)*ca rb(6)=q3*xidot+axis(nsat)*q4*sdot return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine barcor c************************************************** c--calculate barycenter-to-jupiter vector implicit none save double precision angl,ca,sa,t1,t2 c>> for galsat/galsap local communication double precision twopi,degrad,t,rb(6) c--these variables are for calculation of bary to jup vector c--they make use of fact that cofby=cofbx with cos replaced by sin double precision cofbx(7),cofbz(5) c--note that cofbx and cofbz were single precision in the @for version! double precision angbx(7,2),angbz(5,2) c--local common block for subroutines common /local/ twopi,degrad,t,rb,angbx,angbz,cofbx,cofbz c>> integer i1 c>> c this routine calculates the barycenter-to-jupiter shift c vector for cases when galsap is called with negative satellite c number or zero. see lieske, jpl engineering memorandum 314-112 c "additional fortran subroutines for obtaining positions and c partial derivatives of jupiter's galilean satellites", 1 mar 1977. c the short series employed are stored in cofbx and cofbz for the coefficients c and in angbx and angbz for the angles. the series are as follows: c 10**10 xbj = (-1262 - 1267 e1) cos (l1 - psi) -1133 (1 + e2) cos (l2 - psi) c -5715 (1 + e3) cos (l3 - psi) -5668 (1 + e4) cos (l4 - psi) c +12 cos (pi3 - psi) + (68 + 68 e19 + 67 e4) cos (pi4 - psi) c -21 (1 + e4) cos (2l4 - pi4 - psi) c c 10**10 ybj = same as xbj with cos replaced by sin c c 10**10 zbj = -9 sin (l2 - omega2) -18 (1 + e3) sin (l3 - omega3) c -27 (1 + e4) sin (l4 - omega4) +9 sin (l3 - psi) c +(42 + 44 e4) sin (l4 - psi) c c>> c-- do 2 i1=1,7 angl=dmod(angbx(i1,1)+angbx(i1,2)*t,twopi) t1=cofbx(i1) t2=angbx(i1,2) ca=t1*dcos(angl)*1.d-10 sa=t1*dsin(angl)*1.d-10 rb(1)=rb(1)+ca rb(2)=rb(2)+sa rb(4)=rb(4)-sa*t2 2 rb(5)=rb(5)+ca*t2 do 3 i1=1,5 angl=dmod(angbz(i1,1)+angbz(i1,2)*t,twopi) t1=cofbz(i1) t2=angbz(i1,2) ca=t1*dcos(angl)*1.d-10 sa=t1*dsin(angl)*1.d-10 rb(3)=rb(3)+sa 3 rb(6)=rb(6)+ca*t2 return end c************************************************************* **!!i satsap40.f c************************************************************* subroutine satsap c********************************************** c subroutine satsap is not designed to be called, but in case c it is, we'll give the names of the routines and their calling sequences. write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'the following routines are contained in satsap:' write(*,*) 'call chkgal(verson)' write(*,*) ' char' write(*,*) ' in' write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'call mult3g(a, b, c)' write(*,*) ' dp(3,3), dp(3,3), dp(3,3)' write(*,*) ' in in out' write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'call qqdot(t, qpsi11, qpsi21)' write(*,*) ' dp, dp, dp' write(*,*) ' in out out' write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'call revizg()' write(*,*) ' dp argument t passed via common/local/' write(*,*) ' in' write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'call rotg(id, arg, a)' write(*,*) ' int dp dp(3,3)' write(*,*) ' in in out' write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'call updat(argx, kodx, argv, kodv,' write(*,*) ' dp(nmx), int(2,nmx), dp(nmv), int(2,nmv),' write(*,*) ' out, in, out, in,' write(*,*) write(*,*) '.. argz, kodz, nmx, nmv, nmz)' write(*,*) '.. out, in, in, in, in' write(*,*) '.. dp(nmz), int(2,nmz), int, int, int' write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'call unkod(kode, kod, kmin)' write(*,*) ' int(2), int(8), int' write(*,*) ' in, out, out' write(*,*) '- - - - - - - - - - - - - - - - - - - -' return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine chkgal(verson) c********************************************** implicit none save c>> check the common blocks to see if everything's loaded & print version character*(*) verson c****************************************** integer ioin, ioout, iomsg, iopun common /ioblok/ ioin, ioout, iomsg, iopun c****************************************** double precision earay,baray,paray common /ebblok/ earay(28),baray(22),paray(28) c****************************************** double precision trmcod common /trmblk/ trmcod(53) c****************************************** double precision angcod,ratcod,ang,rat common /angblk/ angcod(99),ratcod(99),ang(22),rat(23) c****************************************** c>> for galsat/galsap local communication double precision twopi,degrad,t,rb(6) c--these variables are for calculation of bary to jup vector c--they make use of fact that cofby=cofbx with cos replaced by sin double precision cofbx(7),cofbz(5) c--note that cofbx and cofbz were single precision in the @for version! double precision angbx(7,2),angbz(5,2) c--local common block for subroutines common /local/ twopi,degrad,t,rb,angbx,angbz,cofbx,cofbz c>> c****************************************** c>> local variables to be saved for galsat/galsap and chkgal, qqdot double precision tlast,tlastg,cj,sj,ci,si,cn,sn,cp,sp,phi,phidot, $ p(3,3),q(3,3),qdot(3,3),qmat(3,3), tlastp common /svtloc/ tlast,tlastg,cj,sj,ci,si,cn,sn,cp,sp,phi,phidot, $ p,q,qdot,qmat, tlastp c>> c****************************************** double precision orbecl,orbequ,obl c****************************************** integer k c-- twopi=6.28318530717958648d0 degrad=57.2957795130823208d0 c..kill this write for now c write(iomsg,*) verson c write(iomsg,1001) ioin,ioout,iomsg,iopun c 1001 format(1h ,'common block /ioblok/: [ioin = ',i3,'] ioout = ',i3, c $ ' iomsg = ',i3,' [iopun = ',i3,']') c if (iomsg .ne. ioout) then c write(ioout,*) verson c write(ioout,1001) ioin,ioout,iomsg,iopun c endif c--check to see if common arrays loaded if( (trmcod(52).ne.1.d4) .or. (angcod(6).ne.9.d0) ) then write(iomsg,201) 201 format(1h1,'error--loading routine was not called first') c if (iomsg .ne. ioout) write(ioout,201) stop endif c--set up bary to jupiter shift coefficients first call cofbx(1)=-1262.d0-1267.d0*earay(1) cofbx(2)=-1133.d0*(1.d0+earay(2)) cj=1.d0+earay(3) cofbx(3)=-5715.d0*cj sj=1.d0+earay(4) cofbx(4)=-5668.d0*sj cofbx(5)=12.d0 cofbx(6)=68.d0*(1.d0+earay(19))+67.d0*earay(4) cofbx(7)=-21.d0*sj cofbz(1)=-9.d0 cofbz(2)=-18.d0*cj cofbz(3)=-27.d0*sj cofbz(4)=9.d0 cofbz(5)=42.d0+44.d0*earay(4) do 50 k=1,4 angbx(k,1)=ang(k)-ang(15) angbx(k,2)=rat(k)-rat(15) if(k.ne.1) then angbz(k-1,1)=ang(k)-ang(k+10) angbz(k-1,2)=rat(k)-rat(k+10) endif 50 continue do 51 k=1,2 angbx(k+4,1)=ang(k+7)-ang(15) angbx(k+4,2)=rat(k+7)-rat(15) angbz(k+3,1)=angbx(k+2,1) 51 angbz(k+3,2)=angbx(k+2,2) angbx(7,1)=2.d0*ang(4)-ang(9)-ang(15) angbx(7,2)=2.d0*rat(4)-rat(9)-rat(15) c--this completes setup of bary-to-jupiter arrays c--set up rotation matrices orbecl=paray(26)*(1.d0+earay(26))/degrad cj=dcos(orbecl) sj=dsin(orbecl) orbequ=paray(25)*(1.d0+earay(25))/degrad ci=dcos(orbequ) si=dsin(orbequ) c--note if node has rate, it should be calculated after stat 1 cn=dcos(ang(22)) sn=dsin(ang(22)) c-- set up rotation for ecliptic of 1950 to equator of 1950 obl=paray(27)*(1.d0+earay(27))/degrad call rotg(1,obl,p ) tlast=-2.d20 tlastp=-4.d20 tlastg=-6.d20 return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine mult3g(a,b,c) c**************************************************** implicit none save double precision a(3,3),b(3,3),c(3,3),sum integer i, j, k c-- do 2 i=1,3 do 2 j=1,3 sum=0.d0 do 1 k=1,3 1 sum=sum+a(i,k)*b(k,j) 2 c(i,j)=sum return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine qqdot(t,qpsi11,qpsi21) c**************************************************** c>> calculate q and qdot matrices where c q = r(-node) p(-j) r(-phi) p(-i) c qdot = r(-node) p(-j) rdot(-phi) p(-i) c the variables qpsi11 and qpsi21 are not used by galsat, c but are calculated here so that galsap can easily use them. c 21 feb 91 jay lieske *j* implicit none double precision t,qpsi11,qpsi21 c****************************************** double precision angcod,ratcod,ang,rat common /angblk/ angcod(99),ratcod(99),ang(22),rat(23) c****************************************** c>> local variables to be saved for galsat/galsap and chkgal, qqdot double precision tlast,tlastg,cj,sj,ci,si,cn,sn,cp,sp,phi,phidot, $ p(3,3),q(3,3),qdot(3,3),qmat(3,3), tlastp common /svtloc/ tlast,tlastg,cj,sj,ci,si,cn,sn,cp,sp,phi,phidot, $ p,q,qdot,qmat, tlastp c>> c****************************************** integer l c-- phidot=rat(15) phi=phidot*t+ang(15)-ang(22) cp=dcos(phi) sp=dsin(phi) c--set up matrix to go from jup equ to 1950 ecl q(1,1)=cn*cp-(sn*cj)*sp qpsi11=-cn*sp-(sn*cj)*cp q(1,2)=qpsi11*ci+(sn*sj)*si q(1,3)=-qpsi11*si+(sn*sj)*ci q(2,1)=sn*cp+(cn*cj)*sp qdot(1,1)=qpsi11*phidot qpsi21=-sn*sp+(cn*cj)*cp qdot(2,1)=qpsi21*phidot q(2,2)=qpsi21*ci-(cn*sj)*si q(2,3)=-qpsi21*si-(cn*sj)*ci q(3,1)=sp*sj q(3,2)=(cp*sj)*ci+cj*si q(3,3)=-(cp*sj)*si+cj*ci do 3 l=1,3 qdot(l,2)=-(q(l,1)*phidot)*ci 3 qdot(l,3)=(q(l,1)*phidot)*si qdot(3,1)=cp*sj*phidot c--note if node rate .ne. 0, then place cn and sn after stat 1 c-- and define phidot=rat(15)-rat(22) (rad/day), and add c--increments qdot(1,1)=qdot(1,1)-q(2,1)*rat(22) c-- qdot(2,1)=qdot(2,1)+q(1,1)*rat(22) c-- qdot(1,2)=qdot(1,2)-q(2,2)*rat(22) c-- qdot(2,2)=qdot(2,2)+q(1,2)*rat(22) c-- qdot(1,3)=qdot(1,3)-q(2,3)*rat(22) c-- qdot(2,3)=qdot(2,3)+q(1,3)*rat(22) c-- return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine revizg c********************************************** c>> revise angles that depend on jupiter's g for jupiter/saturn inequality c see lieske, astronomy & astrophysics 56,333-352 (1977) table 3 footnote. c****************************************** implicit none save double precision angcod,ratcod,ang,rat common /angblk/ angcod(99),ratcod(99),ang(22),rat(23) c****************************************** double precision axis(4),cxi1(10),argx1(10),ratx1(10), 1 cz1(7),argz1(7),ratz1(7),cv1(41),argv1(41),ratv1(41),cxi2(24), 2 argx2(24),ratx2(24),cz2(11),argz2(11),ratz2(11),cv2(66), 3 argv2(66),ratv2(66),cxi3(31),argx3(31),ratx3(31),cz3(13), 4 argz3(13),ratz3(13),cv3(75),argv3(75),ratv3(75),cxi4(49), 5 argx4(49),ratx4(49),cz4(18),argz4(18),ratz4(18), 6 cv4(89),argv4(89),ratv4(89),epsln integer nxi1t,nz1t,nv1t,nxi2t,nz2t,nv2t,nxi3t,nz3t, 2 nv3t,nxi4t,nz4t,nv4t integer kodx1,kodz1,kodv1, 2 kodx2,kodz2,kodv2, 3 kodx3,kodz3,kodv3, 4 kodx4,kodz4,kodv4 common /theory/ axis,cxi1,argx1,ratx1,cz1,argz1,ratz1,cv1, 1 argv1,ratv1,cxi2,argx2,ratx2,cz2,argz2,ratz2,cv2,argv2,ratv2, 2 cxi3,argx3,ratx3,cz3,argz3,ratz3,cv3,argv3,ratv3,cxi4,argx4, 3 ratx4,cz4,argz4,ratz4,cv4,argv4,ratv4,epsln,nxi1t,nz1t,nv1t, 4 nxi2t,nz2t,nv2t,nxi3t,nz3t,nv3t,nxi4t,nz4t,nv4t, 5 kodx1(2,10),kodz1(2,7),kodv1(2,41),kodx2(2,24), 6 kodz2(2,11),kodv2(2,66),kodx3(2,31),kodz3(2,13), 7 kodv3(2,75),kodx4(2,49),kodz4(2,18),kodv4(2,89) c****************************************** c>> for galsat/galsap local communication *j* double precision twopi,degrad,t *j*c--local common block for subroutines *j* common /local/ twopi,degrad,t c> only twopi, degrad, t are needed, but we include all to meet f77 standard c****************************************** c--local common block for subroutines double precision twopi,degrad,t,rb(6) double precision cofbx(7),cofbz(5) double precision angbx(7,2),angbz(5,2) common /local/ twopi,degrad,t,rb,angbx,angbz,cofbx,cofbz c>> c****************************************** c>> double precision qx,dg integer k c-- *j* tlastg=t qx=2.d0*ang(16)-ang(17)+.76699/degrad+(2.d0*rat(16)-rat(17))*t qx=dmod(qx,twopi) dg=.03439d0*dsin(qx) qx=5.d0*ang(16)-2.d0*ang(17)+64.26288/degrad 1 +(5.d0*rat(16)-2.d0*rat(17))*t-.02276946941d0/degrad* 2 t/365.25d0 qx=dmod(qx,twopi) dg=(dg+.33033d0*dsin(qx))/degrad do 6 k=86,92 qx=-90+k if(qx.ge.0.d0) qx=qx+1.d0 6 angcod(k)=qx*(ang(17)+dg) call updat(argx1,kodx1,argv1,kodv1,argz1,kodz1,10,41,7) call updat(argx2,kodx2,argv2,kodv2,argz2,kodz2,24,66,11) call updat(argx3,kodx3,argv3,kodv3,argz3,kodz3,31,75,13) call updat(argx4,kodx4,argv4,kodv4,argz4,kodz4,49,89,18) return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine rotg(id,arg,a) c**************************************************** c--these are set up for negative rotations *j* implicit none integer id double precision a(3,3),arg c set up rotation matrix for x,y,z rotations p=x=1, q=y=2, r=z=3 c id 11 12 13 21 22 23 31 32 33 c px 1 1 0 0 0 c -s 0 s c c qy 2 c 0 s 0 1 0 -s 0 c c rz 3 c -s 0 s c 0 0 0 1 integer j, k c-- a(id,id)=1.d0 j=id-1 k=id+1 if (id.eq.1) then j=3 else if (id.eq.3) then k=1 endif a(j,j)=dcos(arg) a(k,k)=a(j,j) a(j,k)=dsin(arg) a(k,j)=-a(j,k) a(j,id)=0.d0 a(id,j)=0.d0 a(k,id)=0.d0 a(id,k)=0.d0 return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine updat(argx,kodx,argv,kodv,argz,kodz,nmx,nmv,nmz) c**************************************************** c called by revizg to update jupiter's mean anomaly for inequalities. implicit none save integer nmx, nmv, nmz integer kodx(2,nmx),kodv(2,nmv),kodz(2,nmz),kod(8) double precision argx(nmx),argv(nmv),argz(nmz) c--updates angles for dg changes c**************************************************** c>> for galsat/galsap local communication *j* double precision twopi *j* common /local/ twopi c> only twopi is needed, but we include all to meet f77 standard c****************************************** c--local common block for subroutines double precision twopi,degrad,t,rb(6) double precision cofbx(7),cofbz(5) double precision angbx(7,2),angbz(5,2) common /local/ twopi,degrad,t,rb,angbx,angbz,cofbx,cofbz c>> c****************************************** c>> c**************************************************** *j* double precision angcod *j* common /angblk/ angcod(99) c> only angcod is needed, but we include all for f77 standard. c> c****************************************** double precision angcod,ratcod,ang,rat common /angblk/ angcod(99),ratcod(99),ang(22),rat(23) c****************************************** c> c**************************************************** integer kl, km, km1, kmin, kmz c-- do 1 kl=1,nmx if ((kodx(1,kl).ne.0).or.(kodx(2,kl).ne.0)) then call unkod(kodx(1,kl),kod,kmin) do 2 km=kmin,8 if((kod(km).ge.86).and.(kod(km).le.92)) then argx(kl)=0.d0 do 4 km1=kmin,8 kmz=kod(km1) 4 argx(kl)=argx(kl)+angcod(kmz) argx(kl)=dmod(argx(kl),twopi) endif 2 continue endif 1 continue c do 8 kl=1,nmv if ((kodv(1,kl).ne.0).or.(kodv(2,kl).ne.0)) then call unkod(kodv(1,kl),kod,kmin) do 5 km=kmin,8 if((kod(km).ge.86).and.(kod(km).le.92)) then argv(kl)=0.d0 do 7 km1=kmin,8 kmz=kod(km1) 7 argv(kl)=argv(kl)+angcod(kmz) argv(kl)=dmod(argv(kl),twopi) endif 5 continue endif 8 continue c do 9 kl=1,nmz if ((kodz(1,kl).ne.0).or.(kodz(2,kl).ne.0)) then call unkod(kodz(1,kl),kod,kmin) do 12 km=kmin,8 if((kod(km).ge.86).and.(kod(km).le.92)) then argz(kl)=0.d0 do 10 km1=kmin,8 kmz=kod(km1) 10 argz(kl)=argz(kl)+angcod(kmz) argz(kl)=dmod(argz(kl),twopi) endif 12 continue endif 9 continue return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine unkod(kode,kod,kmin) c********************************************** implicit none save integer kod(8),kode(2) integer kmin integer kx, kb c-- kod(1)=-kode(1)/1000000 kx=-kode(1)-kod(1)*1000000 kod(2)=kx/10000 kx=kx-kod(2)*10000 kod(3)=kx/100 kod(4)=kx-kod(3)*100 kod(5)=-kode(2)/1000000 kx=-kode(2)-kod(5)*1000000 kod(6)=kx/10000 kx=kx-kod(6)*10000 kod(7)=kx/100 kod(8)=kx-kod(7)*100 do 2 kb=1,8 if(kod(kb).gt.0) go to 3 2 continue 3 kmin=kb+2 return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine qikqip c********************************************** c subroutine qikqip is not designed to be called, but in case c it is, we'll give the names of the routines and their calling sequences. write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'the following routines are contained in qikqip:' write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'call cd2com(verson, lcom)' write(*,*) ' char(*), char(*)(*)' write(*,*) ' in in' write(*,*) '- - - - - - - - - - - - - - - - - - - -' write(*,*) 'call ercard(lab1, lab2)' write(*,*) ' char(*), char(*)' write(*,*) ' in in' write(*,*) '- - - - - - - - - - - - - - - - - - - -' return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = c = = = = = = = = = = = = = = = = = = = = = = = = = = = c subroutine cd2com(verson, lcom) subroutine cd2com c********************************************** c>> read the 'card' file (output as punch file from kodlod or the c>> first 719 cards of kodlop) and fill the necessary theory c>> common blocks /ebblok/ /trmblk/ /angblk/ /theory/ c..kill the pass implicit none save c character*(*) verson c character*(*) lcom(*) c****************************************** integer ioin, ioout, iomsg, iopun common /ioblok/ ioin, ioout, iomsg, iopun c****************************************** c>> common block used by kodqik, kodqip, kodlod, kodlop double precision picon, two31 common /xconst/ picon, two31 **j save /xconst/ c****************************************** character*80 title common /titlod/ title **j save /titlod/ c****************************************** double precision earay,baray,paray common /ebblok/ earay(28),baray(22),paray(28) c****************************************** double precision trmcod common /trmblk/ trmcod(53) c****************************************** double precision angcod,ratcod,ang,rat common /angblk/ angcod(99),ratcod(99),ang(22),rat(23) c****************************************** double precision axis(4),cxi1(10),argx1(10),ratx1(10), 1 cz1(7),argz1(7),ratz1(7),cv1(41),argv1(41),ratv1(41),cxi2(24), 2 argx2(24),ratx2(24),cz2(11),argz2(11),ratz2(11),cv2(66), 3 argv2(66),ratv2(66),cxi3(31),argx3(31),ratx3(31),cz3(13), 4 argz3(13),ratz3(13),cv3(75),argv3(75),ratv3(75),cxi4(49), 5 argx4(49),ratx4(49),cz4(18),argz4(18),ratz4(18), 6 cv4(89),argv4(89),ratv4(89),epsln integer nxi1t,nz1t,nv1t,nxi2t,nz2t,nv2t,nxi3t,nz3t, 2 nv3t,nxi4t,nz4t,nv4t integer kodx1,kodz1,kodv1, 2 kodx2,kodz2,kodv2, 3 kodx3,kodz3,kodv3, 4 kodx4,kodz4,kodv4 common /theory/ axis,cxi1,argx1,ratx1,cz1,argz1,ratz1,cv1, 1 argv1,ratv1,cxi2,argx2,ratx2,cz2,argz2,ratz2,cv2,argv2,ratv2, 2 cxi3,argx3,ratx3,cz3,argz3,ratz3,cv3,argv3,ratv3,cxi4,argx4, 3 ratx4,cz4,argz4,ratz4,cv4,argv4,ratv4,epsln,nxi1t,nz1t,nv1t, 4 nxi2t,nz2t,nv2t,nxi3t,nz3t,nv3t,nxi4t,nz4t,nv4t, 5 kodx1(2,10),kodz1(2,7),kodv1(2,41),kodx2(2,24), 6 kodz2(2,11),kodv2(2,66),kodx3(2,31),kodz3(2,13), 7 kodv3(2,75),kodx4(2,49),kodz4(2,18),kodv4(2,89) c***************************************************************** character*6 lab1,lab2 c****************************************** c..add this instead of the pass character*80 verson character*6 lcom(7) c>> the lcom variables contain the following: data lcom/'ebblok', 'trmblk', 'angblk', 'theory', $ 'start ','end ', 'title ' / data verson/' kodqik version 5.0 jay lieske, jpl'/ c>> c..try a rewind c rewind(unit=ioin) picon=3.14159265358979324d0 two31=2147483647.d0 c..kill this header for now c write(iomsg,*) verson c write(iomsg,1000) ioin,ioout,iomsg,iopun c 1000 format(1h ,'common block /ioblok/: ioin = ',i3,' ioout = ',i3, c 1 ' iomsg = ',i3,' [iopun = ',i3,']') c if (iomsg .ne. ioout) then c write(ioout,*) verson c write(ioout,1000) ioin,ioout,iomsg,iopun c endif c>> add title stuff here c>> version 4 and later supports the /title / field. c>> we also are backwards-compatible with input files that c>> do not support the /title / field. c>> c attempt to read the first cards of the input file. c if they're the /title / field, then proceed. otherwise see c if they represent the pre-4.0 input files that start with /ebblok/. c-- read /title / field read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) if( (lab1.eq.lcom(7)) .and. (lab2.eq.lcom(5)) ) then c>> do title stuff and then read header of /ebblok/ field read(ioin,'(a80)') title read(ioin,100) lab1,lab2 call casedn(title) call casedn(lab1) call casedn(lab2) if((lab1.ne.lcom(7)) .or. (lab2.ne.lcom(6))) $ call ercard(lab1,lab2) c> so far so good. read /ebblok/ field. c--read /ebblok/ read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) 100 format(1x,a6,2x,a6) endif c..kill this write for now c write(ioout,1001) title c1001 format(1x,'/title /: ', a) if((lab1.ne.lcom(1)) .or. (lab2.ne.lcom(5))) $ call ercard(lab1,lab2) read(ioin,101) earay,baray,paray 101 format(3d24.17) read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) if((lab1.ne.lcom(1)) .or. (lab2.ne.lcom(6))) $ call ercard(lab1,lab2) c--read /trmblk/ read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) if((lab1.ne.lcom(2)) .or. (lab2.ne.lcom(5))) $ call ercard(lab1,lab2) read(ioin,101) trmcod read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) if((lab1.ne.lcom(2)) .or. (lab2.ne.lcom(6))) $ call ercard(lab1,lab2) c--read /angblk/ read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) if((lab1.ne.lcom(3)) .or. (lab2.ne.lcom(5))) $ call ercard(lab1,lab2) read(ioin,101) angcod,ratcod,ang,rat read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) if((lab1.ne.lcom(3)) .or. (lab2.ne.lcom(6))) $ call ercard(lab1,lab2) c--read /theory/ read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) if((lab1.ne.lcom(4)) .or. (lab2.ne.lcom(5))) $ call ercard(lab1,lab2) read(ioin,101) axis,cxi1,argx1,ratx1,cz1,argz1,ratz1,cv1,argv1, $ ratv1, 1 cxi2,argx2,ratx2,cz2,argz2,ratz2,cv2,argv2,ratv2, 2 cxi3,argx3,ratx3,cz3,argz3,ratz3,cv3,argv3,ratv3, 3 cxi4,argx4,ratx4,cz4,argz4,ratz4,cv4,argv4,ratv4,epsln read(ioin,102) nxi1t,nz1t,nv1t,nxi2t,nz2t,nv2t,nxi3t,nz3t,nv3t, 1 nxi4t,nz4t,nv4t,kodx1,kodz1,kodv1,kodx2,kodz2,kodv2,kodx3, 2 kodz3,kodv3,kodx4,kodz4,kodv4 102 format(6i11) read(ioin,100) lab1,lab2 call casedn(lab1) call casedn(lab2) if((lab1.ne.lcom(4)) .or. (lab2.ne.lcom(6))) $ call ercard(lab1,lab2) c--all common areas are now loaded for state return end c = = = = = = = = = = = = = = = = = = = = = = = = = = = subroutine ercard(lab1,lab2) c********************************************** c>> wrong cards. print error message and stop. implicit none save character *(*) lab1,lab2 c****************************************** integer ioin, ioout, iomsg, iopun common /ioblok/ ioin, ioout, iomsg, iopun c****************************************** write(iomsg,1001) lab1,lab2 1001 format(1h1 ,'********wrong cards read *****', 2a6) c if (iomsg .ne. ioout) write(ioout,1001) lab1,lab2 stop end