60 subroutine planet_position(jd,pl, lat,lon,hgt, LBaccur,Raccur, ltime,aber,to_fk5,assume_jde, &
61 lunar_theory,nutat,magmdl, verbosity)
63 use sufr_kinds,
only: double
64 use sufr_constants,
only: pi,pi2, r2d,r2h, au,earthr,pland, enpname, jd2000
65 use sufr_system,
only: warn, quit_program_error
66 use sufr_angles,
only: rev, rev2
67 use sufr_dummy,
only: dumdbl1,dumdbl2
68 use sufr_text,
only: d2s
69 use sufr_numerics,
only: deq0
85 real(double),
intent(in) :: jd
86 integer,
intent(in) :: pl
87 real(double),
intent(in),
optional :: lat,lon,hgt, LBaccur,Raccur
88 logical,
intent(in),
optional :: ltime,aber,to_fk5, assume_jde
89 integer,
intent(in),
optional :: lunar_theory, nutat,magmdl, verbosity
91 integer :: j, llunar_theory, lnutat,lmagmdl, lverb
92 real(double) :: jdl,jde,jde_lt, tjm,tjm0, llat,llon,lhgt, lLBaccur,lRaccur, dpsi,eps0,deps,eps,tau,tau1
93 real(double) :: hcl0,hcb0,hcr0, hcl,hcb,hcr, hcl00,hcb00,hcr00, sun_gcl,sun_gcb, gcx,gcy,gcz, gcx0,gcy0,gcz0, dhcr
94 real(double) :: gcl,gcb,delta,gcl0,gcb0,delta0
95 real(double) :: ra,dec,gmst,agst,lst,hh,az,alt,elon, topra,topdec,topl,topb,topdiam,topdelta,tophh,topaz,topalt
96 real(double) :: diam,illfr,pa,magn, parang,parang_ecl,hp, rES1,rES2
97 logical :: lltime,laber,lto_fk5
100 gcl0 = 0.d0; gcb0 = 0.d0; hcr00 = 0.d0
101 gcx = 0.d0; gcy = 0.d0; gcz = 0.d0
102 gcx0 = 0.d0; gcy0 = 0.d0; gcz0 = 0.d0
103 topdelta = 0.d0; delta0 = 0.d0
115 if(
present(lat)) llat = lat
116 if(
present(lon)) llon = lon
117 if(
present(hgt)) lhgt = hgt
122 if(
present(lbaccur)) llbaccur = lbaccur
123 if(
present(raccur)) lraccur = raccur
127 if(
present(ltime)) lltime = ltime
130 if(
present(aber)) laber = aber
133 if(
present(to_fk5)) lto_fk5 = to_fk5
138 if(
present(lunar_theory)) llunar_theory = lunar_theory
139 if(llunar_theory.lt.1 .or. llunar_theory.gt.3) &
140 call quit_program_error(
'planet_position(): lunar_theory must be 1, 2 or 3', 1)
144 if(
present(nutat))
then
146 if(lnutat.ne.1980 .and. lnutat.ne.2000 .and. lnutat.ne.0) &
147 call quit_program_error(
'planet_position(): nutat must be 2000, 1980 or 0', 1)
152 if(
present(magmdl))
then
154 if(lmagmdl.lt.1 .or. lmagmdl.gt.4) &
155 call quit_program_error(
'planet_position(): magmdl must be 1-4', 1)
160 if(
present(verbosity)) lverb = verbosity
169 if(
present(assume_jde) .and. assume_jde)
then
171 jdl = jde -
deltat/86400.d0
173 jde = jdl +
deltat/86400.d0
176 tjm = (jde-jd2000)/365250.d0
179 print*,
'JD: ', d2s(jdl,8)
180 print*,
'DeltaT: ', d2s(
deltat,3)
181 print*,
'JDE: ', d2s(jde,8)
182 print*,
't_jm: ', d2s(tjm,10)
185 call vsop87d_lbr(tjm,3, hcl0,hcb0,hcr0, llbaccur,lraccur)
194 print*,
't_jm: ', d2s(tjm, 10)
195 print*,
'Heliocentric longitude Earth: ', d2s(modulo(hcl0,pi2)*r2d, 9),
' deg'
196 print*,
'Heliocentric latitude Earth: ', d2s(hcb0*r2d, 9),
' deg'
197 print*,
'Heliocentric distance Earth: ', d2s(hcr0, 9),
' au'
202 tau = 0.d0; tau1 = 0.d0
204 do while(deq0(tau1) .or. abs((tau1-tau)/tau1).gt.1.d-10)
207 tjm = tjm0 - tau/365250.d0
208 jde_lt = jd2000 + tjm*365250
209 hcl = 0.d0; hcb = 0.d0; hcr = 0.d0
213 if(pl.gt.0.and.pl.lt.9)
call vsop87d_lbr(tjm,pl, hcl,hcb,hcr, llbaccur,lraccur)
214 if(pl.eq.9)
call plutolbr(tjm*10.d0, hcl,hcb,hcr)
215 if(pl.gt.10000)
call asteroid_lbr(tjm,pl-10000, hcl,hcb,hcr)
217 call vsop87d_lbr(tjm,3, hcl,hcb,hcr, llbaccur,lraccur)
219 select case(llunar_theory)
234 select case(llunar_theory)
244 if(pl.ne.0.and.pl.lt.10 .or. pl.gt.10000)
then
251 if(pl.eq.3) delta = hcr
255 if(pl.gt.10.and.pl.lt.10000) &
256 call cometgc(tjm,tjm0, pl, hcr,gcl,gcb,delta)
261 if(abs(tau).gt.1.d-10)
call warn(
'planet_position(): tau != 0 on first light-time iteration', 0)
276 tau1 = 5.77551830441d-3 * delta
281 print*,
'Iteration: ', j
289 print*,
't_jm: ', d2s(tjm, 10)
290 print*,
'Heliocentric longitude: ', d2s(modulo(hcl,pi2)*r2d, 9),
' deg'
291 print*,
'Heliocentric latitude: ', d2s(hcb*r2d, 9),
' deg'
292 print*,
'Heliocentric distance: ', d2s(hcr, 9),
' au'
294 print*,
'Geocentric x: ', d2s(gcx, 9),
' au'
295 print*,
'Geocentric y: ', d2s(gcy, 9),
' au'
296 print*,
'Geocentric z: ', d2s(gcz, 9),
' au'
298 print*,
'Geocentric longitude: ', d2s(modulo(gcl,pi2)*r2d, 9),
' deg'
299 print*,
'Geocentric latitude: ', d2s(gcb*r2d, 9),
' deg'
300 print*,
'True geocentric distance: ', d2s(delta0, 9),
' au'
301 print*,
'Apparent geocentric distance: ', d2s(delta, 9),
' au'
303 print*,
'JDE_i: ', d2s(jde_lt, 9)
304 print*,
't_jm: ', d2s(tjm, 10)
305 print*,
'Light time: ', d2s(tau,9),
' day'
306 print*,
'Light time1: ', d2s(tau1,9),
' day'
315 call warn(
'planet_position(): Light time failed to converge for '//trim(enpname(pl)), 0)
326 print*,
'Number of iterations: ', j
327 print*,
'Final light time: ', d2s(tau, 9)
328 print*,
'Final t_Jm: ', d2s(tjm, 10)
333 select case(llunar_theory)
364 print*,
'Planet after light-time iterations, before correction for nutation, aberration, FK5:'
365 print*,
'Final t_jm: ', d2s(tjm, 10)
366 print*,
'Final heliocentric longitude: ', d2s(modulo(hcl,pi2)*r2d, 9),
' deg'
367 print*,
'Final heliocentric latitude: ', d2s(hcb*r2d, 9),
' deg'
368 print*,
'Final heliocentric distance: ', d2s(hcr, 9),
' au'
370 print*,
'Final geocentric x: ', d2s(gcx, 9),
' au'
371 print*,
'Final geocentric y: ', d2s(gcy, 9),
' au'
372 print*,
'Final geocentric z: ', d2s(gcz, 9),
' au'
374 print*,
'Final geocentric longitude: ', d2s(modulo(gcl,pi2)*r2d, 9),
' deg'
375 print*,
'Final geocentric latitude: ', d2s(gcb*r2d, 9),
' deg'
376 print*,
'Final (true) geocentric distance: ', d2s(delta, 9),
' au'
378 print*,
'Final jde_i: ', d2s(jde_lt, 9)
379 print*,
'Final t_jm: ', d2s(tjm, 10)
380 print*,
'Final light time: ', d2s(tau,9),
' day'
381 print*,
'Final light time1: ', d2s(tau1,9),
' day'
386 if(lnutat.eq.1980)
then
402 if(lto_fk5)
call fk5(tjm, gcl,gcb)
407 sun_gcl = rev(hcl0+pi)
412 if(lto_fk5)
call fk5(tjm, sun_gcl,sun_gcb)
413 sun_gcl = sun_gcl + dpsi
417 call fk5(tjm, hcl,hcb)
418 call fk5(tjm, hcl0,hcb0)
419 call fk5(tjm, hcl00,hcb00)
424 print*,
'Convert to the FK5 system: ', lto_fk5
425 print*,
'True heliocentric longitude, FK5: ', d2s(modulo(hcl00,pi2)*r2d, 9),
' deg'
426 print*,
'True heliocentric latitude, FK5: ', d2s(hcb00*r2d, 9),
' deg'
427 print*,
'True heliocentric distance: ', d2s(hcr00, 9),
' au'
438 agst = rev(gmst + dpsi*cos(eps))
439 lst = rev(agst + llon)
443 if(pl.ge.0.and.pl.lt.10) diam = atan(pland(pl)/(delta*au))
451 call geoc2topoc_ecl(gcl,gcb, delta,diam/2.d0, eps,lst, topl,topb,topdiam, lat=llat,hgt=lhgt)
452 if(pl.eq.-1) topdiam = res2
454 if(pl.ge.0.and.pl.lt.10) topdelta = pland(pl)/tan(topdiam)/au
458 call eq2horiz(ra,dec,agst, hh,az,alt, lat=llat,lon=llon)
462 elon = acos(cos(gcb)*cos(hcb0)*cos(gcl-rev(hcl0+pi)))
463 if(pl.gt.10) elon = acos((hcr0**2 + delta**2 - hcr**2)/(2*hcr0*delta))
467 call ecl_2_eq(topl,topb,eps, topra,topdec)
468 call eq2horiz(topra,topdec,agst, tophh,topaz,topalt, lat=llat,lon=llon)
473 pa = atan2( hcr0*sin(elon) , delta - hcr0*cos(elon) )
474 else if(pl.gt.0)
then
475 pa = acos( (hcr**2 + delta**2 - hcr0**2) / (2*hcr*delta) )
477 illfr = 0.5d0*(1.d0 + cos(pa))
481 if(pl.gt.0.and.pl.lt.10)
then
483 if(pl.eq.6) magn = magn +
dsatmagn(tjm*10., gcl,gcb)
486 if(pl.eq.0) magn =
moonmagn(pa,delta)
490 if(pl.gt.10 .and. pl.lt.10000)
then
504 parang = atan2(sin(tophh),tan(llat)*cos(topdec) - sin(topdec)*cos(tophh))
505 parang_ecl = atan( (cos(topl)*tan(eps)) / (sin(topl)*sin(topb)*tan(eps) - cos(topb)) )
508 hp = asin(earthr/(delta*au))
512 if(pl.eq.0) hcr = 0.d0
546 print*,
'Final apparent geocentric longitude: ', d2s(rev(gcl)*r2d, 9),
' deg'
547 print*,
'Final apparent geocentric latitude: ', d2s(rev2(gcb)*r2d, 9),
' deg'
548 print*,
'Final apparent geocentric distance: ', d2s(delta, 9),
' au'
549 print*,
'Final apparent heliocentric distance: ', d2s(hcr, 9),
' au'
551 print*,
'Final apparent geocentric right ascension: ', d2s(rev(ra)*r2h, 9),
' hr'
552 print*,
'Final apparent geocentric declination: ', d2s(dec*r2d, 9),
' deg'
594 planpos(41) = rev(hcl0+pi+dpsi)
771 use sufr_kinds,
only: double
772 use sufr_constants,
only: r2d
773 use sufr_system,
only: quit_program_error
776 integer,
intent(in) :: pl
777 integer,
intent(in),
optional :: model
778 integer :: lpl, lmodel
779 real(double),
intent(in) :: hc_dist,gc_dist,phang
784 if(
present(model)) lmodel = model
788 if(lmodel.eq.1 .and. pl.eq.1) pa = pa - 50.d0
796 a0 = [-1.16d0, 4.d0, 0.d0, 1.3d0, 8.93d0, 8.68d0, 6.85d0, 7.05d0, 1.d0] * (-1)
797 a1 = [ 2.838d0, 1.322d0, 0.d0, 1.486d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-2
798 a2 = [ 1.023d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-4
799 a3 = [ 0.d0, 0.4247d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-6
801 else if(lmodel.eq.2)
then
803 a0 = [ 0.42d0, 4.40d0, 0.d0, 1.52d0, 9.40d0, 8.88d0, 7.19d0, 6.87d0, 1.00d0] * (-1)
804 a1 = [ 3.80d0, 0.09d0, 0.d0, 1.6d0, 0.5d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-2
805 a2 = [-2.73d0, 2.39d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-4
806 a3 = [ 2.d0, -0.65d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-6
808 else if(lmodel.eq.3)
then
810 a0 = [ 0.36d0, 4.29d0, 0.d0, 1.52d0, 9.25d0, 8.88d0, 7.19d0, 6.87d0, 1.00d0] * (-1)
811 a1 = [ 3.80d0, 0.09d0, 0.d0, 1.6d0, 0.5d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-2
812 a2 = [-2.73d0, 2.39d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-4
813 a3 = [ 2.d0, -0.65d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-6
815 else if(lmodel.eq.4)
then
817 a0 = [ 0.60d0, 4.47d0, -0.98d0, 1.52d0, 9.40d0, 8.88d0, 7.19d0, 6.87d0, 1.01d0] * (-1)
818 a1 = [ 4.98d0, 1.03d0, -1.02d0, 1.6d0, 0.5d0, 4.4d0, 0.2d0, 0.d0, 0.d0] * 1.d-2
819 a2 = [-4.88d0, 0.57d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-4
820 a3 = [ 3.02d0, 0.13d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0] * 1.d-6
822 if(pl.eq.2 .and. pa.gt.163.6d0) lpl = 3
824 call quit_program_error(
'planet_magnitude(): model must be 1-4', 0)
827 planet_magnitude = 5*log10(hc_dist*gc_dist) + a0(lpl) + a1(lpl)*pa + a2(lpl)*pa2 + a3(lpl)*pa3
976 use sufr_kinds,
only: double
977 use sufr_constants,
only: d2r, pi, jd1900
978 use sufr_angles,
only: rev
986 real(double),
intent(in) :: jd
987 real(double),
intent(out) :: de,ds, omg1,omg2, dphase, pa, in,om
988 real(double) :: jde,d,t1,t2,t3
989 real(double) :: alp0,del0,w1,w2,x,y,z,eps0,eps,dpsi,l,b,alps,dels, r,r0,l0
990 real(double) :: u,v,alp,del,dze,delta,l1,b1
1001 jde = jd +
deltat/86400.d0
1004 d = jde - 2433282.5d0
1007 t2 = (jde - jd1900)/36525.d0
1008 in = (3.120262d0 + 0.0006d0 * t2) * d2r
1010 t3 = jde - 2443000.5d0 -
planpos(7)
1011 om = (316.518203d0 - 2.08362d-6*t3) * d2r
1012 om = om + (1.3966626d0*t1 + 3.088d-4*t1**2)*d2r
1014 alp0 = rev((268.d0 + 0.1061d0*t1)*d2r)
1015 del0 = rev((64.5d0 - 0.0164d0*t1)*d2r)
1018 w1 = rev((17.710d0 + 877.90003539d0*d)*d2r)
1019 w2 = rev((16.838d0 + 870.27003539d0*d)*d2r)
1042 ds = -sin(del0)*sin(dels) - cos(del0)*cos(dels)*cos(alp0-alps)
1045 u = y*cos(eps0) - z*sin(eps0)
1046 v = y*sin(eps0) + z*cos(eps0)
1048 alp = rev(atan2(u,x))
1049 del = atan2(v,sqrt(x**2+u**2))
1050 dze = rev( atan2( sin(del0)*cos(del)*cos(alp0-alp) - sin(del)*cos(del0) , cos(del)*sin(alp0-alp) ) )
1053 de = -sin(del0)*sin(del) - cos(del0)*cos(del)*cos(alp0-alp)
1056 dphase = abs( (2*r*delta + r0**2 - r**2 - delta**2)/(4*r*delta) ) * sin(l-l0)/abs(sin(l-l0))
1059 omg1 = rev(w1 - dze - 5.07033d0*d2r * delta + dphase)
1060 omg2 = rev(w2 - dze - 5.02626d0*d2r * delta + dphase)
1063 call eq_2_ecl(alp0,del0,eps0, l1,b1)
1064 call ecl_2_eq(l1+dpsi,b1,eps, alp0,del0)
1067 pa = atan2( cos(del0)*sin(alp0-alp) , sin(del0)*cos(del) - cos(del0)*sin(del)*cos(alp0-alp) )
1099 subroutine saturnphys(jd, be,bs, pa, in,om, ar,br, du, pa_s, ar_s,br_s)
1100 use sufr_kinds,
only: double
1101 use sufr_constants,
only: pi, as2r,d2r
1102 use sufr_angles,
only: rev
1108 real(double),
intent(in) :: jd
1109 real(double),
intent(out) :: be,bs, pa, in,om, ar,br, du, pa_s, ar_s,br_s
1111 real(double) :: t, lam,bet, d,l,b,r, l0,b0,r0, x,y,z, dpsi,eps, ascn
1112 real(double) :: l2,b2, u1,u2, lam0,bet0, dl,db, ra,dec, ra0,dec0, sinbe
1145 in = (28.075216d0 - 0.012998d0*t + 4.d-6 * t**2) * d2r
1146 om = (169.508470d0 + 1.394681d0*t + 4.12d-4 * t**2) * d2r
1149 sinbe = sin(in)*cos(bet)*sin(lam-om) - cos(in)*sin(bet)
1150 ar = 375.35d0*as2r/d
1158 l2 = l - 0.01759d0*d2r / r
1159 b2 = b - 7.64d-4*d2r * cos(l-ascn) / r
1162 bs = asin(sin(in)*cos(b2)*sin(l2-om) - cos(in)*sin(b2))
1165 u1 = atan2(sin(in)*sin(b2) + cos(in)*cos(b2)*sin(l2-om), cos(b2)*cos(l2-om))
1166 u2 = atan2(sin(in)*sin(bet) + cos(in)*cos(bet)*sin(lam-om), cos(bet)*cos(lam-om))
1176 dl = 0.005693d0*d2r * cos(l0-lam)/cos(bet)
1177 db = 0.005693d0*d2r * sin(l0-lam)*sin(bet)
1188 call ecl_2_eq(lam0,bet0,eps, ra0,dec0)
1191 pa = atan2(cos(dec0)*sin(ra0-ra),sin(dec0)*cos(dec) - cos(dec0)*sin(dec)*cos(ra0-ra))