51 subroutine moonphys(jd, libl,libb, pa,blpa, sunl,sunb, jdEqnx)
52 use sufr_kinds,
only: double
53 use sufr_constants,
only: pi, d2r
54 use sufr_angles,
only: rev, rev2
61 real(double),
intent(in) :: jd
62 real(double),
intent(in),
optional :: jdEqnx
63 real(double),
intent(out) :: libl,libb, pa,blpa, sunl,sunb
64 real(double) :: moonpos(nplanpos),storepos(nplanpos), tjc,tjc2,tjc3,tjc4, dpsi,eps, lm,bm, ram,decm, omg,mal,mem,mas,mam,ee,k1,k2,in,ww,aa
65 real(double) :: rho,sig,tau,libl2,libb2,vv,xx,yy,om, ras,decs, ls,bs,rs, lhc,bhc, sunl2,sunb2
78 if(
present(jdeqnx))
then
92 omg = 2.1824390725d0 - 33.7570464271d0 *tjc + 3.622256d-5 *tjc2 + 3.7337958d-8 *tjc3 - 2.879321d-10 *tjc4
93 mal = 1.62790515798d0 + 8433.46615806092d0 *tjc - 6.37725855d-5 *tjc2 - 4.9498844d-9 *tjc3 + 2.0216715d-11 *tjc4
96 mem = 5.19846946025d0 + 7771.37714617d0 *tjc - 3.340909d-5 *tjc2 + 9.2114446d-8 *tjc3
97 mas = 6.24003588115d0 + 628.301956024d0 *tjc - 2.79776d-6 *tjc2 - 5.8177641733d-8 *tjc3
98 mam = 2.3555483693d0 + 8328.69142288d0 *tjc + 1.517947757d-4 *tjc2 + 3.102807559d-7 *tjc3
99 ee = 1.d0 - 0.002516d0 *tjc - 0.0000074 *tjc2
100 k1 = 2.09003d0 + 2.301199d0 *tjc
101 k2 = 1.2664d0 + 0.352312d0 *tjc
109 aa = atan2( sin(ww)*cos(bm)*cos(in) - sin(bm)*sin(in), cos(ww)*cos(bm) )
112 libb = asin( -sin(ww)*cos(bm)*sin(in) - sin(bm)*cos(in) )
118 - 2.752d-2*cos(mam) - 2.245d-2*sin(mal) + 6.84d-3*cos(mam-2*mal) &
119 - 2.93d-3*cos(2*mal) - 8.5d-4*cos(2*(mal-mem)) - 5.4d-4*cos(mam-2*mem) &
120 - 2.d-4*sin(mam+mal) - 2.d-4*cos(mam+2*mal) - 2.d-4*cos(mam-mal) &
121 + 1.4d-4*cos(mam+2*(mal-mem))
124 - 2.816d-2*sin(mam) + 2.244d-2*cos(mal) - 6.82d-3*sin(mam-2*mal) &
125 - 2.79d-3*sin(2*mal) - 8.3d-4*sin(2*(mal-mem)) + 6.9d-4*sin(mam-2*mem) &
126 + 4.d-4*cos(mam+mal) - 2.5d-4*sin(2*mam) - 2.3d-4*sin(mam+2*mal) &
127 + 2.d-4*cos(mam-mal) + 1.9d-4*sin(mam-mal) &
128 + 1.3d-4*sin(mam+2*(mal-mem)) - 1.d-4*cos(mam-3*mal)
130 tau = 2.520d-2*ee*sin(mas) + 4.74d-3*sin(2*(mam-mal)) - 4.67d-3*sin(mam) &
131 + 3.96d-3*sin(k1) + 2.76d-3*sin(2*(mam-mem)) + 1.96d-3*sin(omg) &
132 - 1.83d-3*cos(mam-mal) + 1.15d-3*sin(mam-2*mem) - 9.6d-4*sin(mam-mem) &
133 + 4.6d-4*sin(2*(mal-mem)) - 3.9d-4*sin(mam-mal) - 3.2d-4*sin(mam-mas-mem) &
134 + 2.7d-4*sin(2*(mam-mem)-mas) + 2.3d-4*sin(k2) - 1.4d-4*sin(2*mem) &
135 + 1.4d-4*cos(2*(mam-mal)) - 1.2d-4*sin(mam-2*mal) - 1.2d-4*sin(2*mam) &
136 + 1.1d-4*sin(2*(mam-mas-mem))
143 libl2 = -tau + (rho*cos(aa) + sig*sin(aa)) * tan(libb)
144 libb2 = sig*cos(aa) - rho*sin(aa)
153 vv = omg + dpsi + sig/sin(in)
154 xx = sin(in+rho) * sin(vv)
155 yy = sin(in+rho) * cos(vv) * cos(eps) - cos(in+rho) * sin(eps)
157 pa = asin( sqrt(xx*xx+yy*yy) * cos(ram-om) / cos(libb) )
167 if(
present(jdeqnx))
then
173 blpa = atan2( cos(decs)*sin(ras-ram), sin(decs)*cos(decm) - cos(decs)*sin(decm)*cos(ras-ram) )
179 lhc = ls + pi + moonpos(4)/rs*cos(bm)*sin(ls-lm)
180 bhc = moonpos(4)/rs*bm
181 ww = lhc - dpsi - omg
182 aa = atan2( sin(ww)*cos(bhc)*cos(in) - sin(bhc)*sin(in), cos(ww)*cos(bhc) )
186 sunb = asin( -sin(ww) * cos(bhc) * sin(in) - sin(bhc) * cos(in))
187 sunl2 = -tau + (rho*cos(aa) + sig*sin(aa)) * tan(sunb)
188 sunb2 = sig*cos(aa) - rho*sin(aa)
189 sunl = rev(sunl + sunl2)
190 sunb = rev2(sunb + sunb2)
212 use sufr_kinds,
only: double, dbl
213 use sufr_constants,
only: d2r
214 use sufr_angles,
only: rev
219 real(double),
intent(in) :: k0
221 real(double) :: t,t2,jde, ee,mas,mam,ff,omg, aa(14),ww
225 phase = nint((k - int(k))*4.0_dbl)
226 if(phase.lt.0) phase = phase + 4
227 if(phase.gt.3) phase = phase - 4
233 jde = 2451550.09766d0 + 29.530588861d0*k + t2*(0.00015437 + t*(-0.000000150 + 0.00000000073*t))
240 ee = 1.0d0 + t * (-0.002516d0 - 0.0000074d0 * t)
243 mas = 2.5534d0 + 29.10535669d0 * k + t2 * (-0.0000218d0 - 0.00000011d0 * t)
246 mam = 201.5643d0 + 385.81693528d0 * k + t2 * (0.0107438d0 + t*(0.00001239d0 - 0.000000058d0 * t))
249 ff = 160.7108d0 + 390.67050274d0 * k + t2 * (-0.0016341d0*t * (-0.00000227d0 + 0.000000011d0 * t))
252 omg = 124.7746d0 - 1.56375580d0 * k + t2 * (0.0020691d0 + 0.00000215d0 * t)
257 aa(1) = 299.77d0 + 0.107408d0 * k - 0.009173d0 * t2
258 aa(2) = 251.88d0 + 0.016321d0 * k
259 aa(3) = 251.83d0 + 26.651886d0 * k
260 aa(4) = 349.42d0 + 36.412478d0 * k
261 aa(5) = 84.66d0 + 18.206239d0 * k
262 aa(6) = 141.74d0 + 53.303771d0 * k
263 aa(7) = 207.14d0 + 2.453732d0 * k
264 aa(8) = 154.84d0 + 7.306860d0 * k
265 aa(9) = 34.52d0 + 27.261239d0 * k
266 aa(10) = 207.19d0 + 0.121824d0 * k
267 aa(11) = 291.34d0 + 1.844379d0 * k
268 aa(12) = 161.72d0 + 24.198154d0 * k
269 aa(13) = 239.56d0 + 25.513099d0 * k
270 aa(14) = 331.55d0 + 3.592518d0 * k
279 aa(i) = rev(aa(i)*d2r)
288 - 0.40720d0 * sin(mam) + 0.17241d0 * ee * sin(mas) &
289 + 0.01608d0 * sin(2*mam) + 0.01039d0 * sin(2*ff) &
290 + 0.00739d0 * ee * sin(mam-mas) - 0.00514d0 * ee * sin(mam+mas) &
291 + 0.00208d0 * ee**2 * sin(2*mas) - 0.00111d0 * sin(mam-2*ff) &
292 - 0.00057d0 * sin(mam+2*ff) + 0.00056d0 * ee * sin(2*mam+mas) &
293 - 0.00042d0 * sin(3*mam) + 0.00042d0 * ee * sin(mas+2*ff) &
294 + 0.00038d0 * ee * sin(mas-2*ff) - 0.00024d0 * ee * sin(2*mam-mas) &
295 - 0.00017d0 * sin(omg) - 0.00007d0 * sin(mam+2*mas) &
296 + 0.00004d0 * sin(2*mam-2*ff) + 0.00004d0 * sin(3*mas) &
297 + 0.00003d0 * sin(mam+mas-2*ff) + 0.00003d0 * sin(2*mam+2*ff) &
298 - 0.00003d0 * sin(mam+mas+2*ff) + 0.00003d0 * sin(mam-mas+2*ff) &
299 - 0.00002d0 * sin(mam-mas-2*ff) - 0.00002d0 * sin(3*mam+mas) &
300 + 0.00002d0 * sin(4*mam)
308 - 0.40614d0 * sin(mam) + 0.17302d0 * ee * sin(mas) &
309 + 0.01614d0 * sin(2*mam) + 0.01043d0 * sin(2*ff) &
310 + 0.00734d0 * ee * sin(mam-mas) - 0.00515d0 * ee * sin(mam+mas) &
311 + 0.00209d0 * ee**2 * sin(2*mas) - 0.00111d0 * sin(mam-2*ff) &
312 - 0.00057d0 * sin(mam+2*ff) + 0.00056d0 * ee * sin(2*mam+mas) &
313 - 0.00042d0 * sin(3*mam) + 0.00042d0 * ee * sin(mas+2*ff) &
314 + 0.00038d0 * ee * sin(mas-2*ff) - 0.00024d0 * ee * sin(2*mam-mas) &
315 - 0.00017d0 * sin(omg) - 0.00007d0 * sin(mam+2*mas) &
316 + 0.00004d0 * sin(2*mam-2*ff) + 0.00004d0 * sin(3*mas) &
317 + 0.00003d0 * sin(mam+mas-2*ff) + 0.00003d0 * sin(2*mam+2*ff) &
318 - 0.00003d0 * sin(mam+mas+2*ff) + 0.00003d0 * sin(mam-mas+2*ff) &
319 - 0.00002d0 * sin(mam-mas-2*ff) - 0.00002d0 * sin(3*mam+mas) &
320 + 0.00002d0 * sin(4*mam)
326 if(phase.eq.1.or.phase.eq.3)
then
328 - 0.62801d0 * sin(mam) + 0.17172d0 * ee * sin(mas) &
329 - 0.01183d0 * ee * sin(mam+mas) + 0.00862d0 * sin(2*mam) &
330 + 0.00804d0 * sin(2*ff) + 0.00454d0 * ee * sin(mam-mas) &
331 + 0.00204d0 * ee**2 * sin(2*mas) - 0.00180d0 * sin(mam-2*ff) &
332 - 0.00070d0 * sin(mam+2*ff) - 0.00040d0 * sin(3*mam) &
333 - 0.00034d0 * ee * sin(2*mam-mas) + 0.00032d0 * ee * sin(mas+2*ff) &
334 + 0.00032d0 * ee * sin(mas-2*ff) - 0.00028d0 * ee**2 * sin(mam+2*mas) &
335 + 0.00027d0 * ee * sin(2*mam+mas) - 0.00017d0 * sin(omg) &
336 - 0.00005d0 * sin(mam-mas-2*ff) &
337 + 0.00004d0 * sin(2*mam+2*ff) - 0.00004d0 * sin(mam+mas+2*ff) &
338 + 0.00004d0 * sin(mam-2*mas) + 0.00003d0 * sin(mam+mas-2*ff) &
339 + 0.00003d0 * sin(3*mas) + 0.00002d0 * sin(2*mam-2*ff) &
340 + 0.00002d0 * sin(mam-mas+2*ff) - 0.00002d0 * sin(3*mam+mas)
343 - 0.00038d0 * ee * cos(mas) + 0.00026d0 * cos(mam) &
344 - 0.00002d0 * cos(mam-mas) + 0.00002d0 * cos(mam+mas) &
345 + 0.00002d0 * cos(2*ff)
347 if(phase.eq.3) ww = -ww
354 + 0.000325d0 * sin(aa(1)) + 0.000165d0 * sin(aa(2)) &
355 + 0.000164d0 * sin(aa(3)) + 0.000126d0 * sin(aa(4)) &
356 + 0.000110d0 * sin(aa(5)) + 0.000062d0 * sin(aa(6)) &
357 + 0.000060d0 * sin(aa(7)) + 0.000056d0 * sin(aa(8)) &
358 + 0.000047d0 * sin(aa(9)) + 0.000042d0 * sin(aa(10)) &
359 + 0.000040d0 * sin(aa(11)) + 0.000037d0 * sin(aa(12)) &
360 + 0.000035d0 * sin(aa(13)) + 0.000023d0 * sin(aa(14))
subroutine planet_position(jd, pl, lat, lon, hgt, lbaccur, raccur, ltime, aber, to_fk5, assume_jde, lunar_theory, nutat, magmdl, verbosity)
Compute the position, distance, etc of a planet.