* * ---------------- * PROGRAM SERIES96 BDL-GF9612 * ---------------- * * This program is an example of use of the subroutine SERIES * and the files of planetary series (96) in the computation of the * astrometric coordinates of a planet (equinox and equator J2000). * * If the user wants to improve the rapidity of computation, it is * recommended to transform the files in direct access and to read * the series in memory once for all in the subroutine SERIES. * * The subroutine SERIES computes the heliocentric coordinates of * planets (equinox and equator J2000). * The subroutine EARTH is necessary for geocentric coordinates. * * ------------ * Declarations * ------------ * implicit double precision (a-h,o-z) * character*3 ext(9) character*8 pla(9) * character*40 fich * dimension v1(6),v2(6),w(3) * data tlight/0.577551831d-2/ data dpi/6.283185307179586d0/ * data ext/'mer','ven','*','mar','jup','sat','ura','nep','plu'/ data pla/'Mercury','Venus','*','Mars','Jupiter','Saturn', . 'Uranus','Neptune','Pluto'/ * * ------------ * Dates limits * ------------ * data t1/2415020.5d0/ data t2/2487980.5d0/ * * ------------------- * Initial julian date * ------------------- * * From : JD2415020.5d0 (1 Jan 1900 0h) * To : JD2487980.5d0 (4 Oct 2099 0h) * t0=2415020.5d0 * * --------------- * Number of dates * --------------- ndat=5 * * ------------------ * Tab interval (day) * ------------------ * step=14500.25d0 * * ---------- * Test dates * ---------- * tlim=t0+(ndat-1)*step if (t0.lt.t1.or.t0.gt.t2) stop ' Error : Date' * * ------------ * Planet index * ------------ * * ipla=1 Mercury * ipla=2 Venus * ipla=4 Mars * ipla=5 Jupiter * ipla=6 Saturn * ipla=7 Uranus * ipla=8 Neptune * ipla=9 Pluto * ipla=1 * * ----------------- * Test planet index * ----------------- * if (ipla.lt.1.or.ipla.eq.3.or.ipla.gt.9) stop ' Error : Planet' * * ---------------- * Open planet file * ---------------- * lupla=ipla+10 fich='series96.'//ext(ipla) open (lupla,file=fich,status='old',iostat=nerr) if (nerr.ne.0) stop ' Error : planet File' * * ------------- * Open EMB file * ------------- * luemb=13 fich='series96.emb' open (luemb,file=fich,status='old',iostat=nerr) if (nerr.ne.0) stop ' Error : EMB File' * * ---------- * Dates loop * ---------- * do n=1,ndat * if (mod(n,10).eq.1) write (*,1001) * tjd=t0+(n-1)*step * * ------------------------------- * Compute astrometric coordinates * ------------------------------- * call SERIES (tjd,3,luemb,v1,ierr) if (ierr.ne.0) stop ' Error : EMB series' call EARTH (tjd,v2) do i=1,3 v1(i)=v1(i)-v2(i) enddo * call SERIES (tjd,ipla,lupla,v2,ierr) if (ierr.ne.0) stop ' Error : planetary series' * do k=1,2 do i=1,3 w(i)=v2(i)-v1(i) enddo dist=dsqrt(w(1)*w(1)+w(2)*w(2)+w(3)*w(3)) t=tjd-dist*tlight call SERIES (t,ipla,lupla,v2,ierr) if (ierr.ne.0) stop ' Error : planetary series' enddo * do i=1,3 w(i)=v2(i)-v1(i) enddo * * --------------- * Display results * --------------- * alpha=datan2(w(2),w(1)) if (alpha.lt.0.0) alpha=alpha+dpi p=dsqrt(w(1)*w(1)+w(2)*w(2)) delta=datan(w(3)/p) * write (*,1002) pla(ipla),tjd,alpha,delta * if (mod(ndat,10).eq.1) then write (*,1003) read (*,*) endif * enddo * * ----------------- * Close file SERIES * ----------------- * close (lupla) close (luemb) * * --- * End * --- * write (*,1004) read (*,*) stop * * ------- * Formats * ------- * 1001 format (//// . 2x,'--------------------------------------------'/ . 2x,'PLANEPH : ASTROMETRIC COORDINATES OF PLANETS'/ . 2x,'--------------------------------------------'/) 1002 format (2x,a8,2x,'JD',f13.5, . 3x,'alpha: ',f12.9,' rad', . 3x,'delta: ',f12.9,' rad') 1003 format (//2x,'Continue : Press Enter'//) 1004 format (//2x,'Program Ended : Press Enter'//) * end * * * subroutine SERIES (tjd,ipla,lu,r,ierr) * ====================================== * * * Ref : Bureau des Longitudes - 96.12 * J. Chapront, G. Francou (BDL) * * * Object * ------ * * Compute planetary ephemerides with series built by a method of * representation using frequency analysis. * Ref : J. Chapront, 1995, Astron. Atrophys. Sup. Ser., 109, 181. * * Planetary Ephemerides : rectangular heliocentric coordinates. * Source : DE403 (Jet Propulsion Laboratory). * Frame : Dynamical Mean Equinox and Equator J2000 (DE403). * Time : Barycentric Dynamical Time (TDB). * * Input * ----- * * tjd : Julian date TDB (double real). * * ipla : Planet index (integer). * 1 : Mercury. * 2 : Venus. * 3 : Earth-Moon Barycenter. * 4 : Mars. * 5 : Jupiter. * 6 : Saturn. * 7 : Uranus. * 8 : Neptune. * 9 : Pluto. * * lu : Logical unit of planet file (integer). * * Output * ------ * * r(6) : Table of rectangular coordinates (double real). * r(1) : X position (au). * r(2) : Y position (au). * r(3) : Z position (au). * r(4) : X' velocity (au/day). * r(5) : Y' velovity (au/day). * r(6) : Z' velocity (au/day). * * ierr : Error index (integer). * 0 : no error. * 1 : date error. * 2 : planet index error. * * Files * ----- * * According to the planet the series are in a file named SERIES.xxx * where xxx is : mer, ven, emb, mar, jup, sat, ura, nep, plu. * The file is opened in the monitor (logical unit lu). * * * Declarations * ------------ * implicit double precision (a-h,o-z) character*3 nom(9) character*12 cname * parameter (maxor=2,maxfq=250) * dimension r(6) dimension v(3),w(3),ws(3) dimension nf(0:maxor),fq(maxfq,0:maxor) dimension sec(0:3,3) dimension ct(maxfq,0:maxor,3) dimension st(maxfq,0:maxor,3) * data nom/'mer','ven','emb','mar','jup','sat','ura','nep','plu'/ * * * Clear results table * ------------------- * do i=1,6 r(i)=0.d0 enddo * * * Read file parameters * --------------------- * rewind lu read (lu,fmt='(a12,f10.1,f8.0,3i3)') cname,tzero,dt,mx,imax,iblock * * * Check up date * ------------- * ierr=1 tmax=tzero+dt*iblock if (tjd.lt.tzero-0.5d0.or.tjd.gt.tmax+0.5d0) return * nb=(tjd-tzero)/dt+1 if (tjd.le.tzero) nb=1 if (tjd.ge.tmax) nb=iblock * * * Check up planet * --------------- * ierr=2 if (ipla.lt.1.or.ipla.gt.9) return if (cname(10:12).ne.nom(ipla)) return * * * Read frequencies * ---------------- * do m=0,mx read (lu,fmt='(i4)') nf(m) do i=1,nf(m) read (lu,fmt='(f20.16)') fq(i,m) enddo enddo * * * Read separator * -------------- * read (lu,*) * * * Read series * ----------- * do k=1,nb do iv=1,3 * * * Read block and variable numbers * ------------------------------- * read (lu,*) * * * Read secular terms * ------------------ * do i=0,imax,2 read (lu,fmt='(2f14.0)') sec(i,iv),sec(i+1,iv) enddo * * * Read Poisson terms * ------------------ * do m=0,mx ip=mod(m,2) do i=1,nf(m) read (lu,fmt='(2f14.0)') a,b if (ip.eq.0) then ct(i,m,iv)=a st(i,m,iv)=b else ct(i,m,iv)=b st(i,m,iv)=a endif enddo enddo * enddo enddo * * * Change variable * --------------- * tdeb=tzero+(nb-1)*dt x=2.d0*(tjd-tdeb)/dt-1.d0 fx=x*dt/2.d0 * * Compute positions (secular terms) * --------------------------------- * do iv=1,3 v(iv)=0.d0 wx=1.d0 max=2*imax-1 do i=0,max v(iv)=v(iv)+sec(i,iv)*wx wx=wx*x enddo enddo * * Compute positions (Poisson terms) * --------------------------------- * wx=1.d0 do m=0,mx nw=nf(m) do iv=1,3 ws(iv)=0.d0 enddo do i=1,nw f=fq(i,m)*fx cf=cos(f) sf=sin(f) do iv=1,3 ws(iv)=ws(iv)+ct(i,m,iv)*cf+st(i,m,iv)*sf enddo enddo do iv=1,3 v(iv)=v(iv)+ws(iv)*wx enddo wx=wx*x enddo * * * Compute velocities (secular terms) * ---------------------------------- * wt=2.d0/dt do iv=1,3 w(iv)=0.d0 wx=1.d0 max=2*imax-1 do i=1,max w(iv)=w(iv)+i*sec(i,iv)*wx wx=wx*x enddo w(iv)=wt*w(iv) enddo * * * Compute velocities (Poisson terms) * ---------------------------------- * wx=1.d0 do m=0,mx nw=nf(m) do i=1,nw fw=fq(i,m) f=fw*fx cf=cos(f) sf=sin(f) do iv=1,3 stw=st(i,m,iv) ctw=ct(i,m,iv) w(iv)=w(iv)+fw*(stw*cf-ctw*sf)*wx if (m.gt.0) w(iv)=w(iv)+m*wt*(ctw*cf+stw*sf)*wy enddo enddo wy=wx wx=wx*x enddo * * Stock results * ------------- * do iv=1,3 r(iv)=v(iv)/1.d10 r(iv+3)=w(iv)/1.d10 enddo * * Exit * ---- * ierr=0 return end * * * subroutine EARTH (tjd,r) * ======================== * * * Ref : Bureau des Longitudes - 96.12 * J. Chapront, G. Francou (BDL) * * * Object * ------ * * Rectangular coordinates of geocentric Earth-Moon barycenter * (equinox and equateur J2000). * * Input * ----- * * tjd : Julian date TDB (double real). * * * Output * ------ * * r(3) : Table of rectangular coordinates (double real). * r(1) : X position (au). * r(2) : Y position (au). * r(3) : Z position (au). * *----------------------------------------------------------------------- * implicit double precision (a-h,o-z) * dimension r(3),v(3),n1(3),n2(3) dimension c(138),s(138),f(138) * data n1/01,44,93/ data n2/43,92,138/ * data c/ .-244075., -2965., 8528., 2345., -2486., 1426., 527., . -43., -393., 394., -218., 73., 91., -173., . 25., -20., 75., 72., 6., 72., -40., . -58., 56., -53., 46., -44., -5., 0., . -1., -12., 21., -4., -4., -2., 9., . 8., 10., 2., -10., -12., -11., 10., . -10.,-176962., -23344., -11109., -922., -4118., 714., . -1135., -601., 299., 564., -311., 261., -251., . 254., 229., 213., -179., 57., -19., -125., . -113., -87., -52., 75., 16., -5., -42., . -4., -10., 8., -19., 9., 40., 29., . -25., 11., 19., -12., -5., 18., -15., . -16., 13., 13., 9., -1., -3., -5., . -1., -76714., 25611., -10120., -400., 1387., -1785., . 310., 580., -492., -527., 44., 130., 244., . -135., 113., -38., 110., 92., -78., 25., . -54., -26., -49., -38., 27., -23., 32., . 2., -2., 1., -18., -2., -23., -4., . 3., -8., 4., -11., 17., 3., -11., . 13., 8., -11., -7., 4./ * data s/ . 192874., 25444., 1005., 4489., -778., 1238., -326., . -614., 339., -285., -276., -232., 195., -63., . 136., 124., 95., 57., -82., 5., 45., . 4., 11., -9., 20., -10., -44., -32., . 27., -20., 6., -20., 17., 17., -14., . -14., 4., -11., -10., 3., 6., -7., . 4.,-223938., -2720., 635., 7824., 2151., -2281., . 1309., -675., 483., -40., -360., 362., 327., . -200., 204., 67., 84., -159., -145., 23., . -18., 69., 66., 5., 65., 66., -36., . -53., 51., -48., 42., -40., -5., 0., . -1., -19., -11., 17., 20., -4., -4., . -2., 9., 7., -9., -13., -11., -10., . 11., -97079., -1464., -1179., 3392., 1557., 933., . -989., -754., 567., -470., 334., 210., -17., . -156., 157., -151., -87., 29., 36., -69., . 10., 43., -8., 30., -39., 29., 2., . 29., 29., -25., -16., -23., 3., 22., . -21., 18., -17., 13., -2., 14., -8., . 0., -9., 0., -8., 9./ * data f/ . 0.2299708345453799, 0.0019436907548255, 0.4579979783362081, . 0.0324605575663244,-0.1955665862245038, 0.4274811115247091, .-0.2318206046403833, 0.6555082553155372, 0.2127688645204655, . 0.2471728045705681, 0.6860251221267625, 0.2604877013568789, . 0.0496625275912389,-0.1783646161993155, 0.0172021241604381, . 0.0191456607800137,-0.2260834530360027, 0.4102791414995209, .-0.0152582792703628, 0.4407960083110198, 0.8835353991060917, . 0.1994539677341547, 0.2662248529612594, 0.4751999483613963, . 0.4427395449305955,-0.0037934608495551, 0.6383062852903491, . 0.4637351299405887, 0.1937168161297741, 0.0152585875411362, .-0.2127685562496920, 0.0000001541352498,-0.4598477484309377, . 0.9140522659175908, 0.2452292679512662, 0.6249913885040383, . 0.2300338378809035,-0.2299078312098563, 0.4446830815498973, . 0.8530185322948665, 0.4885148451477070, 0.0381977091707050, . 0.2147124011397673, 0.2299708345453799, 0.0019436907548255, . 0.2308957195928816, 0.4579979783362081, 0.0324605575663244, .-0.1955665862245038, 0.4274811115247091,-0.0028685758023272, .-0.2318206046403833, 0.6555082553155372, 0.2127688645204655, . 0.2471728045705681, 0.1946417011770021, 0.6860251221267625, . 0.4589228633837098, 0.2604877013568789, 0.0496625275912389, .-0.1783646161993155,-0.0333854426135524, 0.0172021241604381, . 0.0191456607800137,-0.2260834530360027, 0.4102791414995209, .-0.0152582792703628, 0.4284059965722108, 0.4407960083110198, . 0.8835353991060917, 0.1994539677341547, 0.2662248529612594, . 0.4751999483613963, 0.4427395449305955,-0.0037934608495551, . 0.6383062852903491, 0.4637351299405887, 0.1937168161297741, . 0.6564331403627652, 0.0152585875411362, 0.1774397311520876, .-0.2127685562496920, 0.0000001541352498,-0.4598477484309377, . 0.9140522659175908, 0.2452292679512662, 0.6249913885040383, . 0.4446830815498973, 0.6869500071742641, 0.8530185322948665, . 0.4885148451477070, 0.2251585679885010, 0.2299708345453799, . 0.2308957195928816, 0.0019436907548255, 0.4579979783362081, .-0.0028685758023272, 0.0324605575663244,-0.1955665862245038, . 0.1946417011770021, 0.4274811115247091, 0.4589228633837098, .-0.0333854426135524,-0.2318206046403833, 0.6555082553155372, . 0.2127688645204655, 0.2471728045705681, 0.4284059965722108, . 0.6860251221267625, 0.2604877013568789, 0.0496625275912389, .-0.1783646161993155, 0.0172021241604381, 0.6564331403627652, . 0.0191456607800137,-0.2260834530360027, 0.1774397311520876, . 0.4102791414995209,-0.0152582792703628, 0.6869500071742641, . 0.4407960083110198, 0.2251585679885010, 0.8835353991060917, . 0.1994539677341547, 0.4226688449678302, 0.2662248529612594, . 0.4751999483613963, 0.4427395449305955,-0.0037934608495551, . 0.2118436712021903, 0.6383062852903491,-0.0505874126387406, .-0.2614125864043805, 0.4637351299405887, 0.0200705458272416, . 0.1937168161297741, 0.0143333942228611,-0.0181270092079398/ * t=tjd-2451545.0d0 * do iv=1,3 v(iv)=0.d0 do n=n1(iv),n2(iv) x=f(n)*t cx=cos(x) sx=sin(x) v(iv)=v(iv)+c(n)*cx+s(n)*sx enddo v(iv)=v(iv)/1.d10 r(iv)=v(iv) enddo * return end