subroutine gust86 (tjj,isat,pos,vel) *---- satellites d'uranus (laskar 1986, laskar and jacobson, 1987) -------- * * version preliminaire 0.0 - g. francou juin 88. * * en entree : * * tjj date julienne temps dynamique (reel dp). * * isat indice du satellite (entier). * 1 miranda. * 2 ariel. * 3 umbriel. * 4 titania. * 5 oberon. * * en sortie : * pos(3),vel(3) = uranocentric j2000 state vector c..declare implicit double precision (a-h,o-z) save logical b0 dimension pos(3),vel(3),el(6),xu(6),xe(6) dimension xp(3),xs(3),yp(3),ys(3),w(3),posvel(6),vp(3), 1 vs(3) dimension fqn(5),fqe(5),fqi(5),phn(5),phe(5),phi(5) dimension gms(5),rmu(5) dimension trans(3,3),rmat(3,3) common/asat/an(5),ae(5),ai(5) c..various data statements data fqn/4445190.550d-06,2492952.519d-06,1516148.111d-06, . 721718.509d-06,466692.120d-06/ data fqe/20.082d0,6.217d0,2.865d0,2.078d0,0.386d0/ data fqi/-20.309d0,-6.288d0,-2.836d0,-1.843d0,-0.259d0/ data phn/-238051.d-06,3098046.d-06,2285402.d-06, . 856359.d-06,-915592.d-06/ data phe/0.611392d0,2.408974d0,2.067774d0,0.735131d0,0.426767d0/ data phi/5.702313d0,0.395757d0,0.589326d0,1.746237d0,4.206896d0/ data gms/4.4d0,86.1d0,84.0d0,230.0d0,200.0d0/ data gmsu/5794554.5d0/ data alf/76.60666666666667d0/ data del/15.03222222222222d0/ data ua/149597870.0d0/ data vl/299792.458d0/ data t1950/2433282.423d0/ data t2000/2451545.0d0/ data t0/2444239.5d0/ data maxit/2/ data nul/10/ data ((rmat(i,j),j=1,3),i=1,3)/ 1 0.9999256791774783d0,-0.0111815116768724d0,-0.0048590038154553d0, 2 0.0111815116959975d0, 0.9999374845751042d0,-0.0000271625775175d0, 3 0.0048590037714450d0,-0.0000271704492210d0, 0.9999881946023742d0/ data b0/.true./ c..one time initializations if (b0) then b0 = .false. pi = 3.1415926535897932384d0 dpi = 2.d0*pi dgrad = pi/180.0d0 sej = 86400.0d0 sej2 = sej*sej anj = 365.25d0 gmu = gmsu do i=1,5 gmu=gmu-gms(i) enddo do i=1,5 rmu(i)=gmu+gms(i) enddo alf = alf*dgrad del = del*dgrad sa = dsin(alf) ca = dcos(alf) sd = dsin(del) cd = dcos(del) trans(1,1) = sa trans(2,1) = -ca trans(3,1) = 0.d0 trans(1,2) = ca*sd trans(2,2) = sa*sd trans(3,2) = -cd trans(1,3) = ca*cd trans(2,3) = sa*cd trans(3,3) = sd endif c..initialize do i=1,3 pos(i) = 0.0d0 vel(i) = 0.0d0 enddo c..calculate the arguments t = tjj - t0 do i=1,5 an(i) = fqn(i)*t+phn(i) ae(i) = fqe(i)*dgrad/anj*t+phe(i) ai(i) = fqi(i)*dgrad/anj*t+phi(i) an(i) = dmod(an(i),dpi) ae(i) = dmod(ae(i),dpi) ai(i) = dmod(ai(i),dpi) enddo c..calculate the ume50 elliptical elements if (isat .eq. 1) then call mirel (t,rn,rl,rk,rh,rq,rp) else if (isat .eq. 2) then call ariel (t,rn,rl,rk,rh,rq,rp) else if (isat .eq. 3) then call umbel (t,rn,rl,rk,rh,rq,rp) else if (isat .eq. 4) then call titel (t,rn,rl,rk,rh,rq,rp) else if (isat .eq. 5) then call obrel (t,rn,rl,rk,rh,rq,rp) endif el(1)=(rmu(isat)*sej2/rn/rn)**(1.d0/3.d0) rl=dmod(rl,dpi) if (rl.lt.0.d0) rl=rl+dpi el(2)=rl el(3)=rk el(4)=rh el(5)=rq el(6)=rp c..calculate the rectangular ume50 state vector xu call ellipx (el,rmu(isat),xu,rien,0,0) c..calculate the rectangular eme50 state vector do iv=1,3 xe(iv)=0.d0 xe(iv+3)=0.d0 do j=1,3 xe(iv)=xe(iv)+trans(iv,j)*xu(j) xe(iv+3)=xe(iv+3)+trans(iv,j)*xu(j+3) enddo enddo do iv=1,3 pos(iv) = xe(iv) vel(iv) = xe(iv+3) enddo c..fk5 solution do j=1,3 a = 0.0d0 b = 0.0d0 do k=1,3 a = a + rmat(j,k)*xe(k) b = b + rmat(j,k)*xe(k+3) enddo pos(j) = a vel(j) = b enddo c..calculate the geocentric state vector, and return c do iv=1,3 c xs(iv)=xp(iv) + xe(iv) c vs(iv)=vp(iv) + xe(iv+3) c enddo c do iv=1,3 c r(iv) = xs(iv) c r(iv+3) = vs(iv) c enddo return end subroutine mirel (t,rn,rl,rk,rh,rq,rp) * *---- calcul des elements elliptiques de miranda (gust86) -------------- * implicit double precision (a-h,o-z) common/asat/an(5),ae(5),ai(5) save * *---- rn => moyen mouvement (radian/jour) ------------------------------ * rn = 4443522.67d-06 . -34.92d-06*cos(an(1)-3.d0*an(2)+2.d0*an(3)) . +8.47d-06*cos(2.d0*an(1)-6.d0*an(2)+4.d0*an(3)) . +1.31d-06*cos(3.d0*an(1)-9.d0*an(2)+6.d0*an(3)) . -52.28d-06*cos(an(1)-an(2)) . -136.65d-06*cos(2.d0*an(1)-2.d0*an(2)) * *---- rl => longitude moyenne (radian) --------------------------------- * rl = -238051.58d-06 . +4445190.55d-06*t . +25472.17d-06*sin(an(1)-3.d0*an(2)+2.d0*an(3)) . -3088.31d-06*sin(2.d0*an(1)-6.d0*an(2)+4.d0*an(3)) . -318.10d-06*sin(3.d0*an(1)-9.d0*an(2)+6.d0*an(3)) . -37.49d-06*sin(4.d0*an(1)-12.d0*an(2)+8.d0*an(3)) . -57.85d-06*sin(an(1)-an(2)) . -62.32d-06*sin(2.d0*an(1)-2.d0*an(2)) . -27.95d-06*sin(3.d0*an(1)-3.d0*an(2)) * *---- z = k + ih ------------------------------------------------------ * rk = 1312.38d-06*cos(ae(1)) . +71.81d-06*cos(ae(2)) . +69.77d-06*cos(ae(3)) . +6.75d-06*cos(ae(4)) . +6.27d-06*cos(ae(5)) . -123.31d-06*cos(-an(1)+2.d0*an(2)) . +39.52d-06*cos(-2.d0*an(1)+3.d0*an(2)) . +194.10d-06*cos(an(1)) * rh = 1312.38d-06*sin(ae(1)) . +71.81d-06*sin(ae(2)) . +69.77d-06*sin(ae(3)) . +6.75d-06*sin(ae(4)) . +6.27d-06*sin(ae(5)) . -123.31d-06*sin(-an(1)+2.d0*an(2)) . +39.52d-06*sin(-2.d0*an(1)+3.d0*an(2)) . +194.10d-06*sin(an(1)) * *---- zeta = q + ip ---------------------------------------------------- * rq = 37871.71d-06*cos(ai(1)) . +27.01d-06*cos(ai(2)) . +30.76d-06*cos(ai(3)) . +12.18d-06*cos(ai(4)) . +5.37d-06*cos(ai(5)) * rp = 37871.71d-06*sin(ai(1)) . +27.01d-06*sin(ai(2)) . +30.76d-06*sin(ai(3)) . +12.18d-06*sin(ai(4)) . +5.37d-06*sin(ai(5)) * return end subroutine ariel (t,rn,rl,rk,rh,rq,rp) * *---- calcul des elements elliptiques d'ariel (gust86) ----------------- * implicit double precision (a-h,o-z) common/asat/an(5),ae(5),ai(5) save * *---- rn => moyen mouvement (radian/jour) ------------------------------ * rn = 2492542.57d-06 . +2.55d-06*cos(an(1)-3.d0*an(2)+2.d0*an(3)) . -42.16d-06*cos(an(2)-an(3)) . -102.56d-06*cos(2.d0*an(2)-2.d0*an(3)) * *---- rl => longitude moyenne (radian) --------------------------------- * rl = 3098046.41d-06 . +2492952.52d-06*t . -1860.50d-06*sin(an(1)-3.d0*an(2)+2.d0*an(3)) . +219.99d-06*sin(2.d0*an(1)-6.d0*an(2)+4.d0*an(3)) . +23.10d-06*sin(3.d0*an(1)-9.d0*an(2)+6.d0*an(3)) . +4.30d-06*sin(4.d0*an(1)-12.d0*an(2)+8.d0*an(3)) . -90.11d-06*sin(an(2)-an(3)) . -91.07d-06*sin(2.d0*an(2)-2.d0*an(3)) . -42.75d-06*sin(3.d0*an(2)-3.d0*an(3)) . -16.49d-06*sin(2.d0*an(2)-2.d0*an(4)) * *---- z = k + ih ------------------------------------------------------- * rk = -3.35d-06*cos(ae(1)) . +1187.63d-06*cos(ae(2)) . +861.59d-06*cos(ae(3)) . +71.50d-06*cos(ae(4)) . +55.59d-06*cos(ae(5)) . -84.60d-06*cos(-an(2)+2.d0*an(3)) . +91.81d-06*cos(-2.d0*an(2)+3.d0*an(3)) . +20.03d-06*cos(-an(2)+2.d0*an(4)) . +89.77d-06*cos(an(2)) * rh = -3.35d-06*sin(ae(1)) . +1187.63d-06*sin(ae(2)) . +861.59d-06*sin(ae(3)) . +71.50d-06*sin(ae(4)) . +55.59d-06*sin(ae(5)) . -84.60d-06*sin(-an(2)+2.d0*an(3)) . +91.81d-06*sin(-2.d0*an(2)+3.d0*an(3)) . +20.03d-06*sin(-an(2)+2.d0*an(4)) . +89.77d-06*sin(an(2)) * *---- zeta = q + ip ---------------------------------------------------- * rq = -121.75d-06*cos(ai(1)) . +358.25d-06*cos(ai(2)) . +290.08d-06*cos(ai(3)) . +97.78d-06*cos(ai(4)) . +33.97d-06*cos(ai(5)) * rp = -121.75d-06*sin(ai(1)) . +358.25d-06*sin(ai(2)) . +290.08d-06*sin(ai(3)) . +97.78d-06*sin(ai(4)) . +33.97d-06*sin(ai(5)) * return end subroutine umbel (t,rn,rl,rk,rh,rq,rp) * *---- calcul des elements elliptiques d'umbriel (gust86) --------------- * implicit double precision (a-h,o-z) common/asat/an(5),ae(5),ai(5) save * *---- rn => moyen mouvement (radian/jour) ------------------------------ * rn = 1515954.90d-06 . +9.74d-06*cos(an(3)-2.d0*an(4)+ae(3)) . -106.00d-06*cos(an(2)-an(3)) . +54.16d-06*cos(2.d0*an(2)-2.d0*an(3)) . -23.59d-06*cos(an(3)-an(4)) . -70.70d-06*cos(2.d0*an(3)-2.d0*an(4)) . -36.28d-06*cos(3.d0*an(3)-3.d0*an(4)) * *---- rl => longitude moyenne (radian) --------------------------------- * rl1 = 2285401.69d-06 . +1516148.11d-06*t . +660.57d-06*sin(an(1)-3.d0*an(2)+2.d0*an(3)) . -76.51d-06*sin(2.d0*an(1)-6.d0*an(2)+4.d0*an(3)) . -8.96d-06*sin(3.d0*an(1)-9.d0*an(2)+6.d0*an(3)) . -2.53d-06*sin(4.d0*an(1)-12.d0*an(2)+8.d0*an(3)) . -52.91d-06*sin(an(3)-4.d0*an(4)+3.d0*an(5)) . -7.34d-06*sin(an(3)-2.d0*an(4)+ae(5)) . -1.83d-06*sin(an(3)-2.d0*an(4)+ae(4)) . +147.91d-06*sin(an(3)-2.d0*an(4)+ae(3)) rl2 = -7.77d-06*sin(an(3)-2.d0*an(4)+ae(2)) . +97.76d-06*sin(an(2)-an(3)) . +73.13d-06*sin(2.d0*an(2)-2.d0*an(3)) . +34.71d-06*sin(3.d0*an(2)-3.d0*an(3)) . +18.89d-06*sin(4.d0*an(2)-4.d0*an(3)) . -67.89d-06*sin(an(3)-an(4)) . -82.86d-06*sin(2.d0*an(3)-2.d0*an(4)) rl3 = -33.81d-06*sin(3.d0*an(3)-3.d0*an(4)) . -15.79d-06*sin(4.d0*an(3)-4.d0*an(4)) . -10.21d-06*sin(an(3)-an(5)) . -17.08d-06*sin(2.d0*an(3)-2.d0*an(5)) * rl=rl1+rl2+rl3 * *---- z = k + ih ------------------------------------------------------- * rk1 = -0.21d-06*cos(ae(1)) . -227.95d-06*cos(ae(2)) . +3904.69d-06*cos(ae(3)) . +309.17d-06*cos(ae(4)) . +221.92d-06*cos(ae(5)) . +29.34d-06*cos(an(2)) . +26.20d-06*cos(an(3)) . +51.19d-06*cos(-an(2)+2.d0*an(3)) . -103.86d-06*cos(-2.d0*an(2)+3.d0*an(3)) . -27.16d-06*cos(-3.d0*an(2)+4.d0*an(3)) rk2 = -16.22d-06*cos(an(4)) . +549.23d-06*cos(-an(3)+2.d0*an(4)) . +34.70d-06*cos(-2.d0*an(3)+3.d0*an(4)) . +12.81d-06*cos(-3.d0*an(3)+4.d0*an(4)) . +21.81d-06*cos(-an(3)+2.d0*an(5)) . +46.25d-06*cos(an(3)) * rk=rk1+rk2 * rh1 = -0.21d-06*sin(ae(1)) . -227.95d-06*sin(ae(2)) . +3904.69d-06*sin(ae(3)) . +309.17d-06*sin(ae(4)) . +221.92d-06*sin(ae(5)) . +29.34d-06*sin(an(2)) . +26.20d-06*sin(an(3)) . +51.19d-06*sin(-an(2)+2.d0*an(3)) . -103.86d-06*sin(-2.d0*an(2)+3.d0*an(3)) . -27.16d-06*sin(-3.d0*an(2)+4.d0*an(3)) rh2 = -16.22d-06*sin(an(4)) . +549.23d-06*sin(-an(3)+2.d0*an(4)) . +34.70d-06*sin(-2.d0*an(3)+3.d0*an(4)) . +12.81d-06*sin(-3.d0*an(3)+4.d0*an(4)) . +21.81d-06*sin(-an(3)+2.d0*an(5)) . +46.25d-06*sin(an(3)) * rh=rh1+rh2 * *---- zeta = q + ip ---------------------------------------------------- * rq = -10.86d-06*cos(ai(1)) . -81.51d-06*cos(ai(2)) . +1113.36d-06*cos(ai(3)) . +350.14d-06*cos(ai(4)) . +106.50d-06*cos(ai(5)) * rp = -10.86d-06*sin(ai(1)) . -81.51d-06*sin(ai(2)) . +1113.36d-06*sin(ai(3)) . +350.14d-06*sin(ai(4)) . +106.50d-06*sin(ai(5)) * return end subroutine titel (t,rn,rl,rk,rh,rq,rp) * *---- calcul des elements elliptiques de titania (gust86) -------------- * implicit double precision (a-h,o-z) common/asat/an(5),ae(5),ai(5) save * *---- rn => moyen mouvement (radian/jour) ------------------------------ * rn1 = 721663.16d-06 . -2.64d-06*cos(an(3)-2.d0*an(4)+ae(3)) . -2.16d-06*cos(2.d0*an(4)-3.d0*an(5)+ae(5)) . +6.45d-06*cos(2.d0*an(4)-3.d0*an(5)+ae(4)) . -1.11d-06*cos(2.d0*an(4)-3.d0*an(5)+ae(3)) rn2 = -62.23d-06*cos(an(2)-an(4)) . -56.13d-06*cos(an(3)-an(4)) . -39.94d-06*cos(an(4)-an(5)) . -91.85d-06*cos(2.d0*an(4)-2.d0*an(5)) . -58.31d-06*cos(3.d0*an(4)-3.d0*an(5)) . -38.60d-06*cos(4.d0*an(4)-4.d0*an(5)) . -26.18d-06*cos(5.d0*an(4)-5.d0*an(5)) . -18.06d-06*cos(6.d0*an(4)-6.d0*an(5)) * rn=rn1+rn2 * *---- rl => longitude moyenne (radian) -------------------------------- * rl1 = 856358.79d-06 . +721718.51d-06*t . +20.61d-06*sin(an(3)-4.d0*an(4)+3.d0*an(5)) . -2.07d-06*sin(an(3)-2.d0*an(4)+ae(5)) . -2.88d-06*sin(an(3)-2.d0*an(4)+ae(4)) . -40.79d-06*sin(an(3)-2.d0*an(4)+ae(3)) . +2.11d-06*sin(an(3)-2.d0*an(4)+ae(2)) . -51.83d-06*sin(2.d0*an(4)-3.d0*an(5)+ae(5)) . +159.87d-06*sin(2.d0*an(4)-3.d0*an(5)+ae(4)) rl2 = -35.05d-06*sin(2.d0*an(4)-3.d0*an(5)+ae(3)) . -1.56d-06*sin(3.d0*an(4)-4.d0*an(5)+ae(5)) . +40.54d-06*sin(an(2)-an(4)) . +46.17d-06*sin(an(3)-an(4)) . -317.76d-06*sin(an(4)-an(5)) . -305.59d-06*sin(2.d0*an(4)-2.d0*an(5)) . -148.36d-06*sin(3.d0*an(4)-3.d0*an(5)) . -82.92d-06*sin(4.d0*an(4)-4.d0*an(5)) rl3 = -49.98d-06*sin(5.d0*an(4)-5.d0*an(5)) . -31.56d-06*sin(6.d0*an(4)-6.d0*an(5)) . -20.56d-06*sin(7.d0*an(4)-7.d0*an(5)) . -13.69d-06*sin(8.d0*an(4)-8.d0*an(5)) * rl=rl1+rl2+rl3 * *---- z= k + ih ------------------------------------------------------- * rk1 = -0.02d-06*cos(ae(1)) . -1.29d-06*cos(ae(2)) . -324.51d-06*cos(ae(3)) . +932.81d-06*cos(ae(4)) . +1120.89d-06*cos(ae(5)) . +33.86d-06*cos(an(2)) . +17.46d-06*cos(an(4)) . +16.58d-06*cos(-an(2)+2.d0*an(4)) . +28.89d-06*cos(an(3)) . -35.86d-06*cos(-an(3)+2.d0*an(4)) rk2 = -17.86d-06*cos(an(4)) . -32.10d-06*cos(an(5)) . -177.83d-06*cos(-an(4)+2.d0*an(5)) . +793.43d-06*cos(-2.d0*an(4)+3.d0*an(5)) . +99.48d-06*cos(-3.d0*an(4)+4.d0*an(5)) . +44.83d-06*cos(-4.d0*an(4)+5.d0*an(5)) . +25.13d-06*cos(-5.d0*an(4)+6.d0*an(5)) . +15.43d-06*cos(-6.d0*an(4)+7.d0*an(5)) * rk=rk1+rk2 * rh1 = -0.02d-06*sin(ae(1)) . -1.29d-06*sin(ae(2)) . -324.51d-06*sin(ae(3)) . +932.81d-06*sin(ae(4)) . +1120.89d-06*sin(ae(5)) . +33.86d-06*sin(an(2)) . +17.46d-06*sin(an(4)) . +16.58d-06*sin(-an(2)+2.d0*an(4)) . +28.89d-06*sin(an(3)) . -35.86d-06*sin(-an(3)+2.d0*an(4)) rh2 = -17.86d-06*sin(an(4)) . -32.10d-06*sin(an(5)) . -177.83d-06*sin(-an(4)+2.d0*an(5)) . +793.43d-06*sin(-2.d0*an(4)+3.d0*an(5)) . +99.48d-06*sin(-3.d0*an(4)+4.d0*an(5)) . +44.83d-06*sin(-4.d0*an(4)+5.d0*an(5)) . +25.13d-06*sin(-5.d0*an(4)+6.d0*an(5)) . +15.43d-06*sin(-6.d0*an(4)+7.d0*an(5)) * rh=rh1+rh2 * *---- zeta= q + ip ---------------------------------------------------- * rq = -1.43d-06*cos(ai(1)) . -1.06d-06*cos(ai(2)) . -140.13d-06*cos(ai(3)) . +685.72d-06*cos(ai(4)) . +378.32d-06*cos(ai(5)) * rp = -1.43d-06*sin(ai(1)) . -1.06d-06*sin(ai(2)) . -140.13d-06*sin(ai(3)) . +685.72d-06*sin(ai(4)) . +378.32d-06*sin(ai(5)) * return end subroutine obrel (t,rn,rl,rk,rh,rq,rp) * *---- calcul des elements elliptiques d'oberon (gust86) --------------- * implicit double precision (a-h,o-z) common/asat/an(5),ae(5),ai(5) save * *---- rn => moyen mouvement (radian/jour) ----------------------------- * rn1 = 466580.54d-06 . +2.08d-06*cos(2.d0*an(4)-3.d0*an(5)+ae(5)) . -6.22d-06*cos(2.d0*an(4)-3.d0*an(5)+ae(4)) . +1.07d-06*cos(2.d0*an(4)-3.d0*an(5)+ae(3)) . -43.10d-06*cos(an(2)-an(5)) rn2 = -38.94d-06*cos(an(3)-an(5)) . -80.11d-06*cos(an(4)-an(5)) . +59.06d-06*cos(2.d0*an(4)-2.d0*an(5)) . +37.49d-06*cos(3.d0*an(4)-3.d0*an(5)) . +24.82d-06*cos(4.d0*an(4)-4.d0*an(5)) . +16.84d-06*cos(5.d0*an(4)-5.d0*an(5)) * rn=rn1+rn2 * *---- rl => longitude moyenne (radian) -------------------------------- * rl1 = -915591.80d-06 . +466692.12d-06*t . -7.82d-06*sin(an(3)-4.d0*an(4)+3.d0*an(5)) . +51.29d-06*sin(2.d0*an(4)-3.d0*an(5)+ae(5)) . -158.24d-06*sin(2.d0*an(4)-3.d0*an(5)+ae(4)) . +34.51d-06*sin(2.d0*an(4)-3.d0*an(5)+ae(3)) . +47.51d-06*sin(an(2)-an(5)) . +38.96d-06*sin(an(3)-an(5)) . +359.73d-06*sin(an(4)-an(5)) rl2 = 282.78d-06*sin(2.d0*an(4)-2.d0*an(5)) . +138.60d-06*sin(3.d0*an(4)-3.d0*an(5)) . +78.03d-06*sin(4.d0*an(4)-4.d0*an(5)) . +47.29d-06*sin(5.d0*an(4)-5.d0*an(5)) . +30.00d-06*sin(6.d0*an(4)-6.d0*an(5)) . +19.62d-06*sin(7.d0*an(4)-7.d0*an(5)) . +13.11d-06*sin(8.d0*an(4)-8.d0*an(5)) * rl=rl1+rl2 * *---- z = k + ih ------------------------------------------------------ * rk1 = 0.00d-06*cos(ae(1)) . -0.35d-06*cos(ae(2)) . +74.53d-06*cos(ae(3)) . -758.68d-06*cos(ae(4)) . +1397.34d-06*cos(ae(5)) . +39.00d-06*cos(an(2)) . +17.66d-06*cos(-an(2)+2.d0*an(5)) rk2 = 32.42d-06*cos(an(3)) . +79.75d-06*cos(an(4)) . +75.66d-06*cos(an(5)) . +134.04d-06*cos(-an(4)+2.d0*an(5)) . -987.26d-06*cos(-2.d0*an(4)+3.d0*an(5)) . -126.09d-06*cos(-3.d0*an(4)+4.d0*an(5)) . -57.42d-06*cos(-4.d0*an(4)+5.d0*an(5)) . -32.41d-06*cos(-5.d0*an(4)+6.d0*an(5)) . -19.99d-06*cos(-6.d0*an(4)+7.d0*an(5)) . -12.94d-06*cos(-7.d0*an(4)+8.d0*an(5)) * rk=rk1+rk2 * rh1 = 0.00d-06*sin(ae(1)) . -0.35d-06*sin(ae(2)) . +74.53d-06*sin(ae(3)) . -758.68d-06*sin(ae(4)) . +1397.34d-06*sin(ae(5)) . +39.00d-06*sin(an(2)) . +17.66d-06*sin(-an(2)+2.d0*an(5)) rh2 = 32.42d-06*sin(an(3)) . +79.75d-06*sin(an(4)) . +75.66d-06*sin(an(5)) . +134.04d-06*sin(-an(4)+2.d0*an(5)) . -987.26d-06*sin(-2.d0*an(4)+3.d0*an(5)) . -126.09d-06*sin(-3.d0*an(4)+4.d0*an(5)) . -57.42d-06*sin(-4.d0*an(4)+5.d0*an(5)) . -32.41d-06*sin(-5.d0*an(4)+6.d0*an(5)) . -19.99d-06*sin(-6.d0*an(4)+7.d0*an(5)) . -12.94d-06*sin(-7.d0*an(4)+8.d0*an(5)) * rh=rh1+rh2 * *---- zeta = q + ip --------------------------------------------------- * rq = -0.44d-06*cos(ai(1)) . -0.31d-06*cos(ai(2)) . +36.89d-06*cos(ai(3)) . -596.33d-06*cos(ai(4)) . +451.69d-06*cos(ai(5)) * rp = -0.44d-06*sin(ai(1)) . -0.31d-06*sin(ai(2)) . +36.89d-06*sin(ai(3)) . -596.33d-06*sin(ai(4)) . +451.69d-06*sin(ai(5)) * return end subroutine ellipx (ell,rmu,xyz,dxyz,ider,iprt) * *---- ellipx 1.1 18 mars 1986 j. laskar ----------------------------- * * calcul des coordonnees rectangulaires (positions et vitesses) et * de leurs derivees partielles par rapport aux elements elliptiques * a partir des elements elliptiques. * * ell(6) : elements elliptiques a: demi-grand axe * l: longitude moyenne * k: exc*cos(long noeud+ arg peri) * h: exc*sin(long noeud+ arg peri) * q: sin(i/2)*cos(long noeud) * p: sin(i/2)*sin(long noeud) * rmu : constante de gravitation du probleme de deux corps * rmu = g*m1*(1+m2/m1) m1 masse centrale * m2 masse du corps considere * xyz(6) : (1:3) positions et (4:6) vitesses * dxyz(6,7) : derivees partielles * dxyz(i,j)=dron(xyz(i))/dron(ell(j)) * dxyz(i,7)=dron(xyz(i))/dron(rmu) * ider : 0 calcul des xyz seulement * 1 calcul des derivees partielles aussi * iprt : impressions si iprt.ge.1 * * subroutine utilisee : keplkh * * *---- declarations ----------------------------------------------------- * implicit double precision (a-h,o-y) save dimension ell(6),xyz(6),dxyz(6,7) dimension rot(3,2),drotp(3,2),drotq(3,2) dimension tx1(2),tx1t(2),dtx1k(2),dtx1h(2),dtx1tk(2),dtx1th(2) dimension dexy1(2,7),dexy1t(2,7) * *---- elements utiles -------------------------------------------------- * ra=ell(1) rl=ell(2) rk=ell(3) rh=ell(4) rq=ell(5) rp=ell(6) rn=sqrt(rmu/ra**3) phi=sqrt(1.d0-rk**2-rh**2) rki=sqrt(1.d0-rq**2-rp**2) psi=1.d0/(1.d0+phi) * *---- matrice de rotation ---------------------------------------------- * rot(1,1)=1.d0-2.d0*rp**2 rot(1,2)=2.d0*rp*rq rot(2,1)=2.d0*rp*rq rot(2,2)=1.d0-2.d0*rq**2 rot(3,1)=-2.d0*rp*rki rot(3,2)= 2.d0*rq*rki * *---- calcul de la longitude excentrique f ----------------------------- *---- f = anomalie excentrique e + longitude du periapse omegapi ------- * call keplkh (rl,rk,rh,f,it,0) sf =sin(f) cf =cos(f) rlmf =-rk*sf+rh*cf umrsa =rk*cf+rh*sf asr =1.d0/(1.d0-umrsa) rsa =1.d0/asr rna2sr=rn*ra*asr * *---- calcul de tx1 et tx1t -------------------------------------------- * tx1(1) =ra*(cf-psi*rh*rlmf-rk) tx1(2) =ra*(sf+psi*rk*rlmf-rh) tx1t(1)=rna2sr*(-sf+psi*rh*umrsa) tx1t(2)=rna2sr*( cf-psi*rk*umrsa) * *---- calcul de xyz ---------------------------------------------------- * do 110 i=1,3 xyz(i) =0.d0 xyz(i+3)=0.d0 do 100 j=1,2 xyz(i) =xyz(i) +rot(i,j)*tx1(j) xyz(i+3)=xyz(i+3)+rot(i,j)*tx1t(j) 100 continue 110 continue * *---- fin du calcul des coordonnees ------------------------------------ * if (ider.ne.1) then return end if * *---- calcul des derivees partielles ----------------------------------- * *---- calcul de drotp et drotq ----------------------------------------- * drotp(1,1)=-4.d0*rp drotp(1,2)= 2.d0*rq drotp(2,1)= 2.d0*rq drotp(2,2)= 0.d0 drotp(3,1)=-2.d0*rki+2.d0*rp**2/rki drotp(3,2)=-2.d0*rp*rq/rki drotq(1,1)= 0.d0 drotq(1,2)= 2.d0*rp drotq(2,1)= 2.d0*rp drotq(2,2)=-4.d0*rq drotq(3,1)= 2.d0*rp*rq/rki drotq(3,2)= 2.d0*rki-2.d0*rq**2/rki * *---- derivees partielles dans les variales f,k,h ---------------------- * psi2fi=psi**2/phi dkh=rk*rh*psi2fi dk2=rk*rk*psi2fi+psi dh2=rh*rh*psi2fi+psi dtx1k(1)=ra*(-rlmf*dkh-1.d0+psi*rh*sf) dtx1k(2)=ra*( rlmf*dk2 -psi*rk*sf) dtx1h(1)=ra*(-rlmf*dh2 -psi*rh*cf) dtx1h(2)=ra*( rlmf*dkh-1.d0+psi*rk*cf) dtx1tk(1)=rna2sr*( umrsa*dkh+psi*rh*cf)+asr*cf*tx1t(1) dtx1tk(2)=rna2sr*(-umrsa*dk2-psi*rk*cf)+asr*cf*tx1t(2) dtx1th(1)=rna2sr*( umrsa*dh2+psi*rh*sf)+asr*sf*tx1t(1) dtx1th(2)=rna2sr*(-umrsa*dkh-psi*rk*sf)+asr*sf*tx1t(2) * *---- derivees partielles dans les variables l,k,h --------------------- * do 200 i=1,2 dexy1 (i,1)= tx1 (i)/ra dexy1t(i,1)=-tx1t(i)/2.d0/ra dexy1 (i,2)= tx1t(i)/rn dexy1t(i,2)=-rn*asr**3*tx1(i) dexy1 (i,3)= sf*dexy1 (i,2)+dtx1k (i) dexy1t(i,3)= sf*dexy1t(i,2)+dtx1tk(i) dexy1 (i,4)=-cf*dexy1 (i,2)+dtx1h (i) dexy1t(i,4)=-cf*dexy1t(i,2)+dtx1th(i) 200 continue * *---- derivees partielles des variables x,y,z,xdot,ydot,zdot ----------- * do 320 iv=1,4 do 310 i=1,3 dxyz(i ,iv)=0.d0 dxyz(i+3,iv)=0.d0 do 300 j=1,2 dxyz(i ,iv)=dxyz(i ,iv)+rot(i,j)*dexy1 (j,iv) dxyz(i+3,iv)=dxyz(i+3,iv)+rot(i,j)*dexy1t(j,iv) 300 continue 310 continue 320 continue iv=5 do 340 i=1,3 dxyz(i ,iv)=0.d0 dxyz(i+3,iv)=0.d0 do 330 j=1,2 dxyz(i ,iv)=dxyz(i ,iv)+drotq(i,j)*tx1(j) dxyz(i+3,iv)=dxyz(i+3,iv)+drotq(i,j)*tx1t(j) 330 continue 340 continue iv=6 do 360 i=1,3 dxyz(i ,iv)=0.d0 dxyz(i+3,iv)=0.d0 do 350 j=1,2 dxyz(i ,iv)=dxyz(i ,iv)+drotp(i,j)*tx1(j) dxyz(i+3,iv)=dxyz(i+3,iv)+drotp(i,j)*tx1t(j) 350 continue 360 continue * *---- derivees partielles par rapport a rmu ---------------------------- * do 400 i=1,3 dxyz(i ,7)=0.d0 dxyz(i+3,7)=0.5d0*xyz(i+3)/rmu 400 continue * end subroutine keplkh (rl,rk,rh,f,it,iprt) * *---- keplkh 1.0 12 decembre 1985 j. laskar -------------------------- * * resolution de l'equation de kepler en variables longitudes, k, h * *----------------------------------------------------------------------- * implicit double precision (a-h,o-y) save if (rl.eq.0.d0) then f=0.d0 return end if itmax=20 eps=1.0d-16 f0=rl e0=abs(rl) do 110 it=1,itmax k=0 sf=sin(f0) cf=cos(f0) ff0 =f0-rk*sf+rh*cf-rl fpf0=1.d0-rk*cf-rh*sf sdir=ff0/fpf0 100 continue f=f0-sdir*(0.5d0)**k e=abs(f-f0) if (iprt.ge.1) then write (*,*) 'it=',it,'k=',k,'ecart=',e end if if (e.gt.e0) then k=k+1 goto 100 else if (k.eq.0.and.e.le.eps.and.ff0.le.eps) then return else f0=f e0=e end if end if 110 continue end