* PLANETAP *------------------------------------------------------------------------------- * * Object : * Example program for subroutine PLANETAP (Fortran 77, DOS 5.0). * Read comments in subroutine PLANETAP. * *------------------------------------------------------------------------------- * implicit double precision (a-h,o-z) character*14 corps(8) character*7 escscr dimension xyz(6) data corps/'MERCURY','VENUS','E-M BARYCENTER','MARS', . 'JUPITER','SATURN','URANUS','NEPTUNE'/ escscr=char(27)//char(91)//'2J'//char(27)//char(91)//'H' * * Compute the longitudes, latitudes, radius vectors and rectangular * coordinates with velocities of the eight major planets by step of * -150000 days for 4 dates from dj0=2451545.0d0 * dj0=2451545.0d0 pas=-150000.d0 itmax=3 * do 110 ipla=1,8 do 100 it=0,itmax dj=dj0+it*pas C write (*,'(2x,a)') escscr write (*,1000) corps(ipla) call PLANETAP (dj,ipla,dlong,dlat,rvec,xyz,ierr) if (ierr.eq.1) stop ' Wrong planet index' if (ierr.eq.2) stop ' Error in Kepler equation' write (*,1001) dj,dlong,dlat,rvec,(xyz(k),k=1,6) if (ierr.eq.3) write (*,1002) C pause ' Enter a blank line to continue' 100 continue 110 continue * stop1248 * 1000 format (2x,'PROGRAM PLANETAP PLANET : ',a/// . 2x,'Longitude (L), Latitude (B), Radius vector (R)'/ . 2x,'Rectangular heliocentric coordinates X Y Z'/ . 2x,'Rectangular heliocentric velocities X'' Y'' Z'''/ . 2x,'Equinox and ecliptic J2000'///) 1001 format (2x,'JD= ',f14.5// . 2x,'L = ',f12.7,' rad ',3x,'B = ',f12.7,' rad ', . 3x,'R = ',f12.7,' au'// . 2x,'X = ',f12.7,' au ',3x,'Y = ',f12.7,' au ', . 3x,'Z = ',f12.7,' au'// . 2x,'X''= ',f12.8,' au/d',3x,'Y''= ',f12.8,' au/d', . 3x,'Z''= ',f12.8,' au/d'//) 1002 format (2x,'You are outside of the interval 1000-3000;'/ . 2x,'the precision decreases outside of that interval.') * end * * * subroutine PLANETAP (dj,ipla,dlong,dlat,rvec,xyz,ierr) *----------------------------------------------------------------------- * * Object : * Computation of approximative ephemerides of the 8 major planets. * - longitude, latitude, radius vector. * - rectangular heliocentric coordinates and velocities. * Reference frame : Equinox and ecliptic J2000 (JD2451545.0). * * 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 : * dj Julian date tdb (double real). * ipla planet index (integer). * ipla=1 Mercury ipla=5 Jupiter * ipla=2 Venus ipla=6 Saturn * ipla=3 Earth-Moon Barycenter ipla=7 Uranus * ipla=4 Mars ipla=8 Neptune * * Output parameters : * dlong longitude in radians (double real). * dlat latitude in radians (double real). * rvec radius vector in au (double real). * xyz(6) rectangular heliocentric coordinates (double real). * xyz(1) : X in au xyz(4) : X' in au/day. * xyz(2) : Y in au xyz(5) : Y' in au/day. * xyz(3) : Z in au xyz(6) : Z' in au/day. * ierr error index (integer) * ierr=0 no error. * ierr=1 error in planet index ipla. * ierr=2 error in Kepler equation. * ierr=3 date outside of the interval (1000-3000). * * Data tables : * * The tables a, dlm, e, pi, dinc, omega give the mean elements of * classical Keplerian variables of planets limited to t**2 terms. * - a: semi-major axis, * - dlm : mean longitude (degree and arcsecond), * - e : eccentricity, * - pi : longitude of the perihelion (degree and arcsecond), * - dinc : inclination (degree and arcsecond), * - omega : longitude of the ascending node (degree and arcsecond). * * The tables kp, ca, sa give the trigonometric terms to be added * to the mean elements of the semi-major axes. * * The tables kq, cl, sl give the trigonometric terms to be added * to the mean elements of the mean longitudes. * * Precision : * Over the interval 1800-2200, the precision is : * for the longitude of Mercury 4" * for the longitude of Venus 5" * for the longitude of E-M barycenter 6" * for the longitude of Mars 17" * for the longitude of Jupiter 71" * for the longitude of Saturn 81" * for the longitude of Uranus 86" * for the longitude of Neptune 11" * Over the interval 1000-3000, the precision is better than 1.5 * times the precision over 1800-2200. * Outside of the interval 1000-3000, the precision decreases. * * Subroutine called : ELLIPX. * *----------------------------------------------------------------------- * * Declarations. * implicit double precision (a-h,o-z) dimension xyz(6),xxp(6),amas(8) dimension a(3,8),dlm(3,8),e(3,8),pi(3,8),dinc(3,8),omega(3,8) dimension kp(9,8),ca(9,8),sa(9,8),kq(10,8),cl(10,8),sl(10,8) data amas/6023600.d0,408523.5d0,328900.5d0,3098710.d0, . 1047.355d0, 3498.5d0, 22869.d0, 19314.d0/ data sdrad/0.4848136811095360d-5/ data dpi/6.283185307179586d0/ ierr=0 * * Tables for mean elements. * data a/ . 0.3870983098d0, 0.d0, 0.d0, . 0.7233298200d0, 0.d0, 0.d0, . 1.0000010178d0, 0.d0, 0.d0, . 1.5236793419d0, 3.d-10, 0.d0, . 5.2026032092d0, 19132.d-10, -39.d-10, . 9.5549091915d0, -0.0000213896d0, 444.d-10, . 19.2184460618d0 , -3716.d-10, 979.d-10, . 30.1103868694d0 , -16635.d-10, 686.d-10/ * data dlm/ . 252.25090552d0, 5381016286.88982d0, -1.92789d0, . 181.97980085d0, 2106641364.33548d0, 0.59381d0, . 100.46645683d0, 1295977422.83429d0, -2.04411d0, . 355.43299958d0, 689050774.93988d0, 0.94264d0, . 34.35151874d0, 109256603.77991d0, -30.60378d0, . 50.07744430d0, 43996098.55732d0, 75.61614d0, . 314.05500511d0, 15424811.93933d0, -1.75083d0, . 304.34866548d0, 7865503.20744d0, 0.21103d0/ * data e/ . 0.2056317526d0, 0.0002040653d0, -28349.d-10, . 0.0067719164d0, -0.0004776521d0, 98127.d-10, . 0.0167086342d0, -0.0004203654d0, -0.0000126734d0, . 0.0934006477d0, 0.0009048438d0, -80641.d-10, . 0.0484979255d0, 0.0016322542d0, -0.0000471366d0, . 0.0555481426d0, -0.0034664062d0, -0.0000643639d0, . 0.0463812221d0, -0.0002729293d0, 0.0000078913d0, . 0.0094557470d0, 0.0000603263d0, 0.d0/ * data pi/ . 77.45611904d0, 5719.11590d0, -4.83016d0, . 131.56370300d0, 175.48640d0, -498.48184d0, . 102.93734808d0, 11612.35290d0, 53.27577d0, . 336.06023395d0, 15980.45908d0, -62.32800d0, . 14.33120687d0, 7758.75163d0, 259.95938d0, . 93.05723748d0, 20395.49439d0, 190.25952d0, . 173.00529106d0, 3215.56238d0, -34.09288d0, . 48.12027554d0, 1050.71912d0, 27.39717d0/ * data dinc/ . 7.00498625d0, -214.25629d0, 0.28977d0, . 3.39466189d0, -30.84437d0, -11.67836d0, . 0.d0, 469.97289d0, -3.35053d0, . 1.84972648d0, -293.31722d0, -8.11830d0, . 1.30326698d0, -71.55890d0, 11.95297d0, . 2.48887878d0, 91.85195d0, -17.66225d0, . 0.77319689d0, -60.72723d0, 1.25759d0, . 1.76995259d0, 8.12333d0, 0.08135d0/ * data omega/ . 48.33089304d0, -4515.21727d0, -31.79892d0, . 76.67992019d0, -10008.48154d0, -51.32614d0, . 174.87317577d0, -8679.27034d0, 15.34191d0, . 49.55809321d0, -10620.90088d0, -230.57416d0, . 100.46440702d0, 6362.03561d0, 326.52178d0, . 113.66550252d0, -9240.19942d0, -66.23743d0, . 74.00595701d0, 2669.15033d0, 145.93964d0, . 131.78405702d0, -221.94322d0 , -0.78728d0/ * * Tables for trigonometric terms * data kp/ . 69613, 75645, 88306, 59899, 15746, 71087, 142173, 3086, 0, . 21863, 32794, 26934, 10931, 26250, 43725, 53867, 28939, 0, . 16002, 21863, 32004, 10931, 14529, 16368, 15318, 32794, 0, . 6345, 7818, 15636, 7077, 8184, 14163, 1107, 4872, 0, . 1760, 1454, 1167, 880, 287, 2640, 19, 2047, 1454, . 574, 0, 880, 287, 19, 1760, 1167, 306, 574, . 204, 0, 177, 1265, 4, 385, 200, 208, 204, . 0, 102, 106, 4, 98, 1367, 487, 204, 0/ * data ca/ . 4, -13, 11, -9, -9, -3, -1, 4, 0, . -156, 59, -42, 6, 19, -20, -10, -12, 0, . 64, -152, 62, -8, 32, -41, 19, -11, 0, . 124, 621, -145, 208, 54, -57, 30, 15, 0, . -23437, -2634, 6601, 6259, -1507, -1821, 2620, -2115,-1489, . 62911,-119919, 79336, 17814,-24241, 12068, 8306, -4893, 8902, . 389061,-262125,-44088, 8387,-22976, -2093, -615, -9720, 6633, . -412235,-157046,-31430, 37817, -9740, -13, -7449, 9644, 0/ * data sa/ . -29, -1, 9, 6, -6, 5, 4, 0, 0, . -48, -125, -26, -37, 18, -13, -20, -2, 0, . -150, -46, 68, 54, 14, 24, -28, 22, 0, . -621, 532, -694, -20, 192, -94, 71, -73, 0, . -14614, -19828, -5869, 1881, -4372, -2255, 782, 930, 913, . 139737, 0, 24667, 51123, -5102, 7429, -4095, -1976,-9566, . -138081, 0, 37205,-49039,-41901,-33872,-27037,-12474,18797, . 0, 28492,133236, 69654, 52322,-49577,-26430, -3593, 0/ * data kq/ . 3086, 15746, 69613, 59899, 75645, 88306, 12661, 2658, 0, 0, .21863, 32794, 10931, 73, 4387, 26934, 1473, 2157, 0, 0, . 10, 16002, 21863, 10931, 1473, 32004, 4387, 73, 0, 0, . 10, 6345, 7818, 1107, 15636, 7077, 8184, 532, 10, 0, . 19, 1760, 1454, 287, 1167, 880, 574, 2640, 19,1454, . 19, 574, 287, 306, 1760, 12, 31, 38, 19, 574, . 4, 204, 177, 8, 31, 200, 1265, 102, 4, 204, . 4, 102, 106, 8, 98, 1367, 487, 204, 4, 102/ * data cl/ . 21, -95, -157, 41, -5, 42, 23, 30, 0, 0, . -160, -313, -235, 60, -74, -76, -27, 34, 0, 0, . -325, -322, -79, 232, -52, 97, 55, -41, 0, 0, . 2268, -979, 802, 602, -668, -33, 345, 201, -55, 0, . 7610, -4997,-7689,-5841,-2617, 1115, -748, -607, 6074, 354, . -18549, 30125,20012, -730, 824, 23, 1289, -352,-14767,-2062, .-135245,-14594, 4197,-4030,-5630,-2898, 2540, -306, 2939, 1986, . 89948, 2103, 8963, 2695, 3682, 1648, 866, -154, -1963, -283/ * data sl/ . -342, 136, -23, 62, 66, -52, -33, 17, 0, 0, . 524, -149, -35, 117, 151, 122, -71, -62, 0, 0, . -105, -137, 258, 35, -116, -88, -112, -80, 0, 0, . 854, -205, -936, -240, 140, -341, -97, -232, 536, 0, . -56980, 8016, 1012, 1448,-3024,-3710, 318, 503, 3767, 577, . 138606,-13478,-4964, 1441,-1319,-1482, 427, 1236, -9167,-1918, . 71234,-41116, 5334,-4935,-1848, 66, 434,-1748, 3780, -701, . -47645, 11647, 2166, 3194, 679, 0, -244, -419, -2531, 48/ * * Wrong planet index. * if (ipla.lt.1.or.ipla.gt.8) then ierr=1 else * * Time parameters. * t=(dj-2451545.d0)/365250.d0 t2=t**2 * * Computation of the mean elements. * da=a(1,ipla)+t*a(2,ipla)+t2*a(3,ipla) dl=(3600.d0*dlm(1,ipla)+t*dlm(2,ipla)+t2*dlm(3,ipla))*sdrad de=e(1,ipla)+t*e(2,ipla)+t2*e(3,ipla) dp=(3600.d0*pi(1,ipla)+t*pi(2,ipla)+t2*pi(3,ipla))*sdrad dp=mod(dp,dpi) if (dp.lt.0.d0) dp=dp+dpi di=(3600.d0*dinc(1,ipla)+t*dinc(2,ipla)+t2*dinc(3,ipla))*sdrad do=(3600.d0*omega(1,ipla)+t*omega(2,ipla) . +t2*omega(3,ipla))*sdrad do=mod(do,dpi) if (do.lt.0.d0) do=do+dpi * * Computation of the trigonometric terms. * dmu=0.35953620*t do 100 j=1,8 arga=kp(j,ipla)*dmu argl=kq(j,ipla)*dmu da=da+(ca(j,ipla)*cos(arga)+sa(j,ipla)*sin(arga))*1.d-7 dl=dl+(cl(j,ipla)*cos(argl)+sl(j,ipla)*sin(argl))*1.d-7 100 continue arga=kp(9,ipla)*dmu da=da+t*(ca(9,ipla)*cos(arga)+sa(9,ipla)*sin(arga))*1.d-7 do 200 j=9,10 argl=kq(j,ipla)*dmu dl=dl+t*(cl(j,ipla)*cos(argl)+sl(j,ipla)*sin(argl))*1.d-7 200 continue dl=mod(dl,dpi) if (dl.lt.0.d0) dl=dl+dpi * * Computation of the rectangular coordinates and velocities. * dmas=1.d0/amas(ipla) call ELLIPX (da,dl,de,dp,di,do,xxp,dmas,ierr) if (ierr.eq.0) then do 300 i=1,6 xyz(i)=xxp(i) 300 continue * * Date outside of the interval (1000-3000). * if (dj.gt.2817151.5d0.or.dj.lt.2086307.5d0) ierr=3 * * Computation of longitude, latitude, radius vector. * dlong=datan2(xxp(2),xxp(1)) rvec=sqrt(xxp(1)**2+xxp(2)**2+xxp(3)**2) dlat=dasin(xxp(3)/rvec) * else ierr=2 dlong=0.d0 dlat=0.d0 rvec=0.d0 xyz(1)=0.d0 xyz(2)=0.d0 xyz(3)=0.d0 endif endif * return end * * * subroutine ELLIPX (a,dlm,e,p,dia,omega,xxp,dmas,ierr) *----------------------------------------------------------------------- * * Object : * This subroutine computes the variables x, y, z, x',y',z' * from the elliptic variables. * * Error index : * ierr=0 : no error. * ierr=1 : error in Kepler equation. * * Intrinsic complex functions : * dcmplex : complex number (*16) with two real numbers (*8). * dimag : imaginary part (*8) of complex number (*16). * dreal : real part (*8) of complex number (*16). * * Subroutine called : KEPLER. * *----------------------------------------------------------------------- * implicit double precision (a-h,o-y) complex*16 z,zto,dcmplx dimension xxp(6) ierr=0 * xa=a xl=dlm xk=e*cos(p) xh=e*sin(p) xq=sin(dia/2)*cos(omega) xp=sin(dia/2)*sin(omega) * asr=648000/3.141592653589793d0 gk=3548.1876069651d0 xfi=sqrt(1.d0-xk**2-xh**2) xki=sqrt(1.d0-xq**2-xp**2) u=1.d0/(1.d0+xfi) z=dcmplx(xk,xh) call KEPLER (xl,z,zto,r,u,ierr) if (ierr.ne.0) return xcw=dreal(zto) xsw=dimag(zto) xm=xp*xcw-xq*xsw xxx=xa*r xxp(1)=xxx*(xcw-2*xp*xm) xxp(2)=xxx*(xsw+2*xq*xm) xxp(3)=-2*xxx*xki*xm xms=xa*(xh+xsw)/xfi xmc=xa*(xk+xcw)/xfi xn=gk*((1.d0+dmas)/xa**3)**(0.5d0)/asr xxp(4)=xn*((-1.d0+2*xp**2)*xms+2*xp*xq*xmc) xxp(5)=xn*((1.d0-2*xq**2)*xmc-2*xp*xq*xms) xxp(6)=2*xn*xki*(xp*xms+xq*xmc) * return end * * * subroutine KEPLER (al,z,zto,r,u,ierr) *----------------------------------------------------------------------- * * Object : * Kepler equation. * * Error index : * ierr=0 : no error. * ierr=1 : error in algorithm. * * Intrinsic complex functions : * dcmplex : complex number (*16) with two real numbers (*8). * dimag : imaginary part (*8) of complex number (*16). * dreal : real part (*8) of complex number (*16). * dconjg : conjugate (*16) of complex number (*16). * cdabs : absolute value (*8) of complex number (*16). * cdexp : exponent (*16) of complex number (*16). * *----------------------------------------------------------------------- * integer ierr,k,kmax complex*16 dcmplx,dconjg,cdexp double precision dreal,dimag double precision ex,al,e,dl,u,r complex*16 z1,z2,z3,zteta,z,zto data kmax/10/ k=0 ierr=0 * al=mod(al,2*3.141592653589793d0) ex=cdabs(z) z1=dconjg(z) e=al+(ex-0.125*ex**3)*sin(al)+0.5*ex**2*sin(2*al) . +0.375*ex**3*sin(al*3) * 1 continue k=k+1 z2=dcmplx(0.d0,e) zteta=cdexp(z2) z3=z1*zteta dl=al-e+dimag(z3) r=1-dreal(z3) if (abs(dl).lt.1.d-12) then z1=u*z*dimag(z3) z2=dcmplx(dimag(z1),-dreal(z1)) zto=(-z+zteta+z2)/r return else if (k.gt.kmax) then ierr=1 return endif endif e=e+dl/r goto 1 * end