* PRCBDL94 *----------------------------------------------------------------------- * * Object : * Example program for subroutine PRCBDL94 (Fortran 77, DOS 5.0). * Read comments in subroutine PRCBDL94. * *----------------------------------------------------------------------- * implicit double precision (a-h,o-z) character*7 escscr dimension qp(12) data sdrad/0.4848136811095360d-5/ escscr=char(27)//char(91)//'2J'//char(27)//char(91)//'H' * * Fixed epoch : 4 dates since J2000 with -200000 days step. * Date epoch : 4 dates since fixed epoch with +100000 days step. * do 300 ief=1,4 ef=2451545.0d0-2.d5*(ief-1) do 200 ied=1,4 ed=ef+1.d5*(ied-1) call PRCBDL94 (ef,ed,qp) do 100 i=1,12 qp(i)=qp(i)/sdrad 100 continue C write (*,'(2x,a)') escscr write (*,1000) ed,ef,qp C pause ' Enter a blank line to continue' 200 continue 300 continue stop1248 * 1000 format (2x,'PRECESSION QUANTITIES (BDL94)'// . 2x,'Date : JD',f14.5, . 20x,'Fixed epoch : JD',f14.5// . 2x,'sin(pi.a)*sin(PI.a) :',f12.5,' "'/ . 2x,'sin(pi.a)*cos(PI.a) :',f12.5,' "'// . 2x,'pi.a : ',f12.5,' "'/ . 2x,'PI.a : ',f12.5,' "'/ . 2x,'p.a : ',f12.5,' "'// . 2x,'theta.a : ',f12.5,' "'/ . 2x,'zeta.a : ',f12.5,' "'/ . 2x,'z.a : ',f12.5,' "'// . 2x,'epsilon.a : ',f12.5,' "'// . 2x,'omega.a : ',f12.5,' "'/ . 2x,'psi.a : ',f12.5,' "'/ . 2x,'chi.a : ',f12.5,' "'/) end * * * subroutine PRCBDL94 (epocf,epocd,qp) *----------------------------------------------------------------------- * * Object : * Computation of precession quantities with formulae BDL94. * Precessional displacements from a fixed epoch "epocf" to the epoch * of the date "epocd". * * Reference : * Astron. Astrophys. 282, 663 (1994). * * Authors : * J.L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze, * G. Francou, J. Laskar (Bureau des Longitudes, Paris, France). * * Input parameters : * epocf Julian date fixed epoch (real double). * epocd Julian date epoch of the date (real double). * * Output parameters : * qp(12) precessional quantities in radians (real double). * * Notations : * For each array hereunder, the indices correspond respectively to * coefficients of power of time : * T = (fixed epoch - J2000 epoch) * and t = (epoch of the date - fixed epoch) * in thousands of Julian years of 365250 days. * These coefficients are double real expressed in 1/10000 ". * ss (0:5,1:5) : sin(pi.a)*sin(PI.a) ----> qp(1) * sc (0:5,1:5) : sin(pi.a)*cos(PI.a) ----> qp(2) * spi (0:5,1:4) : pi.a ----> qp(3) * cpi (0:5,0:5) : PI.a ----> qp(4) * p (0:5,1:6) : p.a (precession) ----> qp(5) * theta (0:5,1:6) : theta.a ----> qp(6) * zeta (0:5,1:6) : zeta.a ----> qp(7) * z (0:5,1:6) : z.a ----> qp(8) * epsil (0:5,0:5) : epsilon.a (obliquity) ----> qp(9) * omega (0:5,0:6) : omega.a ----> qp(10) * psi (0:5,1:6) : psi.a ----> qp(11) * chi (0:5,1:6) : chi.a ----> qp(12) * *----------------------------------------------------------------------- * * Declarations. * implicit double precision (a-h,o-z) dimension qp(12),qpl(12) dimension ct(0:5),st(0:6),lct(0:6,12),lst(2,12) dimension coef(0:5,0:6,12) dimension ss(0:5,5),sc(0:5,5) dimension spi(0:5,4),cpi(0:5,0:5),p(0:5,6) dimension theta(0:5,6),zeta(0:5,6),z(0:5,6) dimension epsil(0:5,0:5),omega(0:5,0:6),psi(0:5,6),chi(0:5,6) equivalence (coef(0,1,1),ss(0,1)),(coef(0,1,2),sc(0,1)) equivalence (coef(0,1,3),spi(0,1)),(coef(0,0,4),cpi(0,0)) equivalence (coef(0,1,5),p(0,1)),(coef(0,1,6),theta(0,1)) equivalence (coef(0,1,7),zeta(0,1)),(coef(0,1,8),z(0,1)) equivalence (coef(0,0,9),epsil(0,0)),(coef(0,0,10),omega(0,0)) equivalence (coef(0,1,11),psi(0,1)),(coef(0,1,12),chi(0,1)) save qpl,epocf0,epocd0 * * Initializations. * data epocf0/0.d0/,epocd0/0.d0/ data sdrad/0.4848136811095360d-5/ data t2000/2451545.d0/ data tunit/365250.d0/ data lst/1,5,1,5,1,4,0,5,1,6,1,6,1,6,1,6,0,5,0,6,1,6,1,6/ data lct/0,5,4,3,2,0,0,0,5,4,3,2,0,0, . 0,4,4,3,2,0,0,5,4,3,2,1,0,0, . 0,5,4,3,2,1,0,0,5,4,3,2,1,0, . 0,5,4,3,2,1,0,0,5,4,3,2,1,0, . 5,5,4,3,2,1,0,5,0,4,3,2,1,0, . 0,5,4,3,2,1,0,0,5,4,3,2,1,0/ * * Data for numerical expressions of precession with following constants : * - general precession : 50288.200 "/1000y * - obliquity in J2000 : 84381.412 " * - Ratio Sun mass / Mercury mass : 6023600 * - Ratio Sun mass / Venus mass : 408523.5 * - Ratio Sun mass / E-M Baryc. mass : 328900.5 * - Ratio Sun mass / Mars mass : 3098710 * - Ratio Sun mass / Jupiter mass : 1047.355 * - Ratio Sun mass / Saturn mass : 3498.5 * - Ratio Sun mass / Uranus mass : 22869 * - Ratio Sun mass / Neptune mass : 19314 * data ss/ . 419971.d0, -753286.d0, 3179.d0, 3178.d0, 7.d0, -4.d0, T0-5 t1 . 193971.d0, 5740.d0, -2541.d0, -5.d0, 6.d0, 0.d0, T0-4 t2 . -2235.d0, 859.d0, 33.d0, -3.d0, 0.d0, 0.d0, T0-3 t3 . -104.d0, -4.d0, 2.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . 2.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t5 data sc/ . -4680927.d0, -305.d0, 59967.d0, -205.d0, -125.d0, -2.d0, T0-5 t1 . 51043.d0, -31633.d0, -326.d0, 138.d0, -2.d0, 0.d0, T0-4 t2 . 5223.d0, 318.d0, -66.d0, -4.d0, 0.d0, 0.d0, T0-3 t3 . -57.d0, 19.d0, -1.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . -1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t5 data spi/ . 4699729.d0, -67011.d0, 448.d0, -19.d0, -1.d0, 0.d0, T0-4 t1 . -33505.d0, 448.d0, -28.d0, -2.d0, 1.d0, 0.d0, T0-4 t2 . -1237.d0, -4.d0, -2.d0, 1.d0, 0.d0, 0.d0, T0-3 t3 . 3.d0, -1.d0, 1.d0, 0.d0, 0.d0, 0.d0/ T0-2 t4 data cpi/ .6295434330.d0,329296590.d0, 953520.d0, -50.d0,-4590.d0, -100.d0, T0-5 t0 . -86792700.d0, -158510.d0, -1130.d0, -4480.d0, -190.d0, 0.d0, T0-4 t1 . 153420.d0, -190.d0, -4320.d0, -230.d0, 0.d0, 0.d0, T0-3 t2 . 50.d0, -2080.d0, -150.d0, 0.d0, 0.d0, 0.d0, T0-2 t3 . -370.d0, -50.d0, 0.d0, 0.d0, 0.d0, 0.d0, T0-1 t4 . -10.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t5 data p/ . 502882000.d0, 2224045.d0, 2095.d0, -9408.d0, -90.d0, 10.d0, T0-5 t1 . 1112022.d0, 2095.d0, -14111.d0, -180.d0, 26.d0, 0.d0, T0-4 t2 . 773.d0, -9410.d0, -180.d0, 35.d0, 0.d0, 0.d0, T0-3 t3 . -2353.d0, -90.d0, 26.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . -18.d0, 10.d0, 0.d0, 0.d0, 0.d0, 0.d0, T0-1 t5 . 2.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t6 data theta/ . 200420207.d0, -853131.d0, -2111.d0, 3642.d0, 8.d0, -5.d0, T0-5 t1 . -426566.d0, -2111.d0, 5463.d0, 17.d0, -12.d0, 0.d0, T0-4 t2 . -418238.d0, 359.d0, 27.d0, -1.d0, 0.d0, 0.d0, T0-3 t3 . -731.d0, 19.d0, 9.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . -127.d0, 11.d0, 0.d0, 0.d0, 0.d0, 0.d0, T0-1 t5 . 4.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t6 data zeta/ . 230609097.d0, 1397495.d0, -38.d0, -5918.d0, -37.d0, 7.d0, T0-5 t1 . 302226.d0, -2523.d0, -3840.d0, -14.d0, 7.d0, 0.d0, T0-4 t2 . 180183.d0, -1326.d0, 6.d0, 5.d0, 0.d0, 0.d0, T0-3 t3 . -583.d0, -1.d0, 7.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . -285.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, T0 t5 . -2.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t6 data z/ . 230609097.d0, 1397495.d0, -38.d0, -5918.d0, -37.d0, 7.d0, T0-5 t1 . 1095270.d0, 2446.d0, -13913.d0, -134.d0, 26.d0, 0.d0, T0-4 t2 . 182667.d0, -11400.d0, -173.d0, 44.d0, 0.d0, 0.d0, T0-3 t3 . -2821.d0, -93.d0, 32.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . -301.d0, 6.d0, 0.d0, 0.d0, 0.d0, 0.d0, T0-1 t5 . -1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t6 data epsil/ . 843814120.d0, -4680927.d0, -152.d0, 19989.d0, -51.d0, -25.d0, T0-5 t0 . -4680927.d0, -305.d0, 59967.d0, -205.d0, -125.d0, -2.d0, T0-5 t1 . -152.d0, 59967.d0, -308.d0, -250.d0, -6.d0, 0.d0, T0-4 t2 . 19989.d0, -205.d0, -250.d0, -8.d0, 0.d0, 0.d0, T0-3 t3 . -51.d0, -125.d0, -6.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . -25.d0, -2.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0-1 t5 data omega/ . 843814120.d0, -4680927.d0, -152.d0, 19989.d0, -51.d0, -25.d0, T0-5 t0 . 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, t1 . 51294.d0, -91954.d0, 298.d0, 389.d0, 2.d0, 0.d0, T0-4 t2 . -77276.d0, 235.d0, 987.d0, -1.d0, 0.d0, 0.d0, T0-3 t3 . -48.d0, 954.d0, -7.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . 333.d0, -9.d0, 0.d0, 0.d0, 0.d0, 0.d0, T0-1 t5 . -3.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t6 data psi/ . 503850672.d0, 492595.d0, -1344.d0, -2115.d0, 17.d0, 3.d0, T0-5 t1 . -1072374.d0, -10919.d0, 13673.d0, 137.d0, -28.d0, 0.d0, T0-4 t2 . -11424.d0, 26425.d0, 87.d0, -111.d0, 0.d0, 0.d0, T0-3 t3 . 13279.d0, -110.d0, -170.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . -94.d0, -123.d0, 0.d0, 0.d0, 0.d0, 0.d0, T0-1 t5 . -35.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t6 data chi/ . 1055794.d0, -1888214.d0, -1888.d0, +7950.d0, 101.d0, -9.d0, T0-5 t1 . -2381379.d0, -10910.d0, 30291.d0, +290.d0, -59.d0, 0.d0, T0-4 t2 . -12117.d0, 39055.d0, 229.d0, -159.d0, 0.d0, 0.d0, T0-3 t3 . 17024.d0, -38.d0, -214.d0, 0.d0, 0.d0, 0.d0, T0-2 t4 . -77.d0, -145.d0, 0.d0, 0.d0, 0.d0, 0.d0, T0-1 t5 . -40.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0/ T0 t6 * if (epocf.ne.epocf0.or.epocd.ne.epocd0) then * * Time parameters * ct(0)=1.d0 st(0)=1.d0 ct(1)=(epocf-t2000)/tunit st(1)=(epocd-epocf)/tunit do 100 i=2,5 ct(i)=ct(1)*ct(i-1) 100 continue do 110 i=2,6 st(i)=st(1)*st(i-1) 110 continue * * Computation of precession quantities. * do 220 k=1,12 qpl(k)=0.d0 jmin=lst(1,k) jmax=lst(2,k) if (epocd.eq.epocf) jmax=0 do 210 j=jmin,jmax q=0.d0 imax=lct(j,k) if (epocf.eq.t2000) imax=0 do 200 i=0,imax q=q+coef(i,j,k)*ct(i) 200 continue qpl(k)=qpl(k)+q*st(j) 210 continue qpl(k)=qpl(k)*1.d-4*sdrad 220 continue epocf0=epocf epocd0=epocd * endif * * Results. * do 300 i=1,12 qp(i)=qpl(i) 300 continue * return end