libTheSky
Routines to compute sky positions of Sun, Moon, planets and more
All Namespaces Files Functions Variables Pages
moon_routines.f90
Go to the documentation of this file.
1!> \file moon_routines.f90 Procedures that calculate the physical data, phases and age of the Moon for libTheSky
2
3
4! Copyright (c) 2002-2025 Marc van der Sluys - Nikhef/Utrecht University - marc.vandersluys.nl
5!
6! This file is part of the libTheSky package,
7! see: https://libthesky.sf.net/
8!
9! This is free software: you can redistribute it and/or modify it under the terms of the
10! European Union Public Licence 1.2 (EUPL 1.2).
11!
12! This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
13! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
14! See the EU Public Licence for more details.
15!
16! You should have received a copy of the European Union Public Licence along with this code.
17! If not, see <https://www.eupl.eu/1.2/en/>.
18
19
20!***********************************************************************************************************************************
21!> \brief Procedures for the Moon
22
24 implicit none
25 save
26
27contains
28
29
30 !*********************************************************************************************************************************
31 !> \brief Get physical data for the Moon: librations, position angles, selenographic position of the Sun
32 !!
33 !! \param jd Julian day for computation (epoch)
34 !!
35 !! \param libl Libration (physical+optical) in longitude (output)
36 !! \param libb Libration (physical+optical) in latitude (output)
37 !!
38 !! \param pa Position angle of the Moon's axis/north pole w.r.t. an equatorial grid (output)
39 !! \param blpa Position angle of the Moon's bright limb w.r.t. an equatorial grid (output)
40 !!
41 !! \param sunl Selenographic longitude of the Sun (output)
42 !! \param sunb Selenographic latitude of the Sun (output)
43 !!
44 !! \param jdEqnx Julian day for equinox (optional; default: jd = JD of epoch)
45 !!
46 !!
47 !! \see Meeus, Astronomical Algorithms, 1998, Ch. 53
48 !!
49 !! \note This routine does NOT save Moon data in planpos !!!
50
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
55
59
60 implicit none
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
66
67
68 storepos = planpos ! Store current data
69
70 call planet_position(jd,0) ! Moon
71 moonpos = planpos
72
73 lm = moonpos(1) ! Geocentric longitude of the Moon
74 bm = moonpos(2) ! Geocentric latitude of the Moon
75 ram = moonpos(5) ! Geocentric right ascension of the Moon
76 decm = moonpos(6) ! Geocentric declination of the Moon
77
78 if(present(jdeqnx)) then ! Precess from EoD -> selected equinox
79 call precess_ecl(jd,jdeqnx, lm,bm)
80 call precess_eq(jd, jdeqnx, ram,decm)
81 end if
82
83 tjc = moonpos(46) ! Apparent dynamical time in Julian Centuries since 2000.0
84 tjc2 = tjc**2 ! t^2
85 tjc3 = tjc2*tjc ! t^3
86 tjc4 = tjc2**2 ! t^4
87
88 dpsi = moonpos(47) ! Nutation in longitude
89 eps = moonpos(48) ! True obliquity of the ecliptic, corrected for nutation
90
91
92 omg = 2.1824390725d0 - 33.7570464271d0 *tjc + 3.622256d-5 *tjc2 + 3.7337958d-8 *tjc3 - 2.879321d-10 *tjc4 ! Moon's longitude of mean ascending node
93 mal = 1.62790515798d0 + 8433.46615806092d0 *tjc - 6.37725855d-5 *tjc2 - 4.9498844d-9 *tjc3 + 2.0216715d-11 *tjc4 ! Moon's argument of latitude (F)
94
95 ! Lower accuracy(?) from Meeus Ch. 22:
96 mem = 5.19846946025d0 + 7771.37714617d0 *tjc - 3.340909d-5 *tjc2 + 9.2114446d-8 *tjc3 ! Mean elongation of the Moon
97 mas = 6.24003588115d0 + 628.301956024d0 *tjc - 2.79776d-6 *tjc2 - 5.8177641733d-8 *tjc3 ! Mean anomaly of the Sun (or Earth)
98 mam = 2.3555483693d0 + 8328.69142288d0 *tjc + 1.517947757d-4 *tjc2 + 3.102807559d-7 *tjc3 ! Mean anomaly of the Moon
99 ee = 1.d0 - 0.002516d0 *tjc - 0.0000074 *tjc2
100 k1 = 2.09003d0 + 2.301199d0 *tjc
101 k2 = 1.2664d0 + 0.352312d0 *tjc
102
103
104 ! Optical librations:
105 in = 2.69203d-2 ! Moon equator inclination to ecliptic, Meeus p.372
106
107 ! Meeus, Eq. 53.1:
108 ww = lm - dpsi - omg
109 aa = atan2( sin(ww)*cos(bm)*cos(in) - sin(bm)*sin(in), cos(ww)*cos(bm) )
110
111 libl = aa - mal ! Along with AA: - = W on celestial sphere, E in selenographic coordinates - opposed to Sterrengids
112 libb = asin( -sin(ww)*cos(bm)*sin(in) - sin(bm)*cos(in) )
113
114
115 ! Physical librations:
116 ! Meeus, p.373, in degrees:
117 rho = &
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))
122
123 sig = &
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)
129
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))
137
138 rho = rho*d2r
139 sig = sig*d2r
140 tau = tau*d2r
141
142 ! Meeus, Eq. 53.2:
143 libl2 = -tau + (rho*cos(aa) + sig*sin(aa)) * tan(libb)
144 libb2 = sig*cos(aa) - rho*sin(aa)
145
146 ! Total librations:
147 libl = libl + libl2 ! Along with AA ( - = W on celestial sphere, E in selenographic coordinates)
148 libb = libb + libb2
149
150
151 ! Position Angle of the axis w.r.t. equatorial grid:
152 ! Meeus, p.374:
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)
156 om = atan2(xx,yy)
157 pa = asin( sqrt(xx*xx+yy*yy) * cos(ram-om) / cos(libb) )
158
159
160 call planet_position(jd,3) ! Sun - CHECK - need full accuracy?
161 ls = planpos(1) ! Geocentric ecliptic longitude of the Sun
162 bs = planpos(2) ! Geocentric ecliptic latitude of the Sun
163 rs = planpos(3) ! Geocentric distance of the Sun
164 ras = planpos(5) ! Geocentric right ascension of the Sun
165 decs = planpos(6) ! Geocentric declination of the Sun
166
167 if(present(jdeqnx)) then ! Precess from EoD -> selected equinox
168 call precess_ecl(jd,jdeqnx, ls,bs)
169 call precess_eq(jd, jdeqnx, ras,decs)
170 end if
171
172 ! Position angle of the bright limb:
173 blpa = atan2( cos(decs)*sin(ras-ram), sin(decs)*cos(decm) - cos(decs)*sin(decm)*cos(ras-ram) )
174
175
176 ! Selenographic position of the Sun (Meeus, p.376):
177
178 ! Heliocentric l,b:
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) )
183
184 ! Selenographic coordinates of the Sun:
185 sunl = aa - mal
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)
191
192
193 ! Restore current data:
194 planpos = storepos
195
196 end subroutine moonphys
197 !*********************************************************************************************************************************
198
199
200 !*********************************************************************************************************************************
201 !> \brief Calculates Julian Day of phase of the Moon for the desired phase k0
202 !!
203 !! \param k0 Desired phase: x.00: New Moon - k = x.75: Last Quarter. k=0 ~ 2000.0
204 !! \retval moonphase JD of the desired lunar phase
205 !!
206 !! \see Meeus, Astronomical Algorithms, 1998, Ch. 49
207 !!
208 !! \note
209 !! - First version: June 2003
210
211 function moonphase(k0)
212 use sufr_kinds, only: double, dbl
213 use sufr_constants, only: d2r
214 use sufr_angles, only: rev
215
216 use thesky_local, only: deltat
217
218 implicit none
219 real(double), intent(in) :: k0
220 real(double) :: moonphase,k
221 real(double) :: t,t2,jde, ee,mas,mam,ff,omg, aa(14),ww
222 integer :: phase,i
223
224 k = k0
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
228
229 t = k/1236.85d0 ! Time in Julian Centuries
230 t2 = t*t
231
232 ! First approximation:
233 jde = 2451550.09766d0 + 29.530588861d0*k + t2*(0.00015437 + t*(-0.000000150 + 0.00000000073*t)) ! Meeus, Eq. 49.1
234
235 ! Return very rough estimate (~day):
236 !moonphase = jde - deltat/86400.d0
237 !return
238
239 ! Correction parameters:
240 ee = 1.0d0 + t * (-0.002516d0 - 0.0000074d0 * t) ! Meeus, Eq. 47.6
241
242 ! Mean anomaly of the Sun, Meeus, Eq. 49.4:
243 mas = 2.5534d0 + 29.10535669d0 * k + t2 * (-0.0000218d0 - 0.00000011d0 * t)
244
245 ! Mean anomaly of the Moon, Meeus, Eq. 49.5:
246 mam = 201.5643d0 + 385.81693528d0 * k + t2 * (0.0107438d0 + t*(0.00001239d0 - 0.000000058d0 * t))
247
248 ! Moon's argument of latitude, Meeus, Eq. 49.6:
249 ff = 160.7108d0 + 390.67050274d0 * k + t2 * (-0.0016341d0*t * (-0.00000227d0 + 0.000000011d0 * t))
250
251 ! Longitude of ascending node of lunar orbit, Meeus, Eq. 49.7:
252 omg = 124.7746d0 - 1.56375580d0 * k + t2 * (0.0020691d0 + 0.00000215d0 * t)
253
254
255 ! Planetary arguments:
256 ! Meeus, p.351:
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
271
272 ! Convert degrees -> radians:
273 mas = rev(mas*d2r)
274 mam = rev(mam*d2r)
275 ff = rev(ff*d2r)
276 omg = rev(omg*d2r)
277
278 do i=1,14
279 aa(i) = rev(aa(i)*d2r)
280 end do
281
282
283
284 !*******************************************************************************************************************************
285 ! New Moon (phase 0), Meeus p.351:
286 if(phase.eq.0) then
287 jde = jde &
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)
301 end if
302
303
304 !*******************************************************************************************************************************
305 ! Full Moon (phase 2), Meeus p.351:
306 if(phase.eq.2) then
307 jde = jde &
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)
321 end if
322
323
324 !*******************************************************************************************************************************
325 ! First, last quarter (phase 1,3), Meeus p.352:
326 if(phase.eq.1.or.phase.eq.3) then
327 jde = jde &
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)
341
342 ww = 0.00306 &
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)
346
347 if(phase.eq.3) ww = -ww
348 jde = jde + ww
349 end if
350
351
352 ! Final corrections, Meeus p.352:
353 jde = jde &
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))
361
362
363 moonphase = jde - deltat/86400.d0 ! JDE - DeltaT = JD
364
365 end function moonphase
366 !*********************************************************************************************************************************
367
368
369
370 !*********************************************************************************************************************************
371 !> \brief Compute the age of the Moon for a given JD
372 !!
373 !! \param jd Julian day for computation
374 !! \retval moon_age Age of the Moon (since last New Moon) in days
375
376 function moon_age(jd)
377 use sufr_kinds, only: double
378 use sufr_constants, only: jd2000
379
380 implicit none
381 real(double), intent(in) :: jd
382 real(double) :: moon_age, djd, jd0, k0
383
384
385 ! Moon age:
386 djd = jd - jd2000 - 0.5d0
387
388 ! Overshoot by a month, to make sure you don't miss one. Age.gt.29.530589d0 doesn't work, since this is the MEAN month:
389 k0 = floor(djd/365.25*12.3685d0) + 1.d0
390 moon_age = -1.d0
391
392 do while(moon_age.lt.0.d0)
393 jd0 = moonphase(k0)
394 moon_age = jd-jd0
395 k0 = k0 - 1.d0 ! Previous New Moon
396 end do
397
398 end function moon_age
399 !*********************************************************************************************************************************
400
401
402 !*********************************************************************************************************************************
403 !> \brief Compute the next lunar phase after a given JD
404 !!
405 !! \param JD Julian day for computation
406 !! \param phase Next lunar phase: 0-NM ... 3-LQ (output)
407 !! \param JDnext Julian day of next lunar phase (output)
408
409 subroutine moon_phase_next(JD, phase, JDnext)
410 use sufr_kinds, only: double
411 use sufr_constants, only: jd2000
412 use thesky_datetime, only: printdate1
413
414 implicit none
415 real(double), intent(in) :: JD
416 integer, intent(out) :: phase
417 real(double), intent(out) :: JDnext
418 real(double) :: dJD, k0
419
420
421 ! Undershoot by half a month, to make sure you don't miss anything:
422 djd = jd - jd2000 - 0.5d0
423 k0 = floor(djd/365.25*12.3685d0) - 0.5d0
424 jdnext = jd - 30.d0
425
426 do while(jdnext.lt.jd)
427 k0 = k0 + 0.25d0 ! Next phase
428 jdnext = moonphase(k0)
429 end do
430
431 phase = nint( (k0 - floor(k0))*4.d0 ) ! 0-3 for NM - LQ
432
433 end subroutine moon_phase_next
434 !*********************************************************************************************************************************
435
436
437end module thesky_moonroutines
438!***********************************************************************************************************************************
Procedures for coordinates.
subroutine precess_eq(jd1, jd2, a1, d1)
Compute the precession of the equinoxes in equatorial coordinates, from jd1 to jd2.
subroutine precess_ecl(jd1, jd2, l, b)
Compute the precession of the equinoxes in geocentric ecliptical coordinates.
Date and time procedures.
Definition date_time.f90:49
subroutine printdate1(jd, jde)
Prints date/time of a given Julian day (UT) to standard output, but without a newline.
Local parameters for libTheSky: location, date, time.
Definition modules.f90:56
real(double) deltat
Current value of DeltaT (s)
Definition modules.f90:66
Procedures for the Moon.
real(double) function moonphase(k0)
Calculates Julian Day of phase of the Moon for the desired phase k0.
subroutine moon_phase_next(jd, phase, jdnext)
Compute the next lunar phase after a given JD.
real(double) function moon_age(jd)
Compute the age of the Moon for a given JD.
subroutine moonphys(jd, libl, libb, pa, blpa, sunl, sunb, jdeqnx)
Get physical data for the Moon: librations, position angles, selenographic position of the Sun.
Planet data, needed to compute planet positions.
Definition modules.f90:85
integer, parameter nplanpos
Number of entries in the planpos array.
Definition modules.f90:116
real(double), dimension(nplanpos) planpos
Planpos[] is an array with many different types of coordinates and related variables:
Definition modules.f90:192
Procedures for planets.
Definition planets.f90:23
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.
Definition planets.f90:63