libTheSky
Routines to compute sky positions of Sun, Moon, planets and more
Loading...
Searching...
No Matches
moon_position.f90
Go to the documentation of this file.
1!> \file moon_position.f90 Core procedures that calculate the position and magnitude of the Moon for libTheSky
2
3
4! Copyright (c) 2002-2023 Marc van der Sluys - marc.vandersluys.nl
5!
6! This file is part of the libTheSky package,
7! see: http://libthesky.sf.net/
8!
9! This is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License
10! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
11!
12! This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
13! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with this code. If not, see
16! <http://www.gnu.org/licenses/>.
17
18
19!***********************************************************************************************************************************
20!> \brief Procedures for the Moon
21
23 implicit none
24 save
25
26contains
27
28 !*********************************************************************************************************************************
29 !> \brief Calculate the apparent geocentric ecliptical position of the Moon, using the Lunar Solution ELP 2000-82B
30 !!
31 !! \param tjj Time for calculation in Julian millenia after 2000.0
32 !!
33 !! \param ll Apparent geocentric ecliptical longitude (output)
34 !! \param bb Apparent geocentric ecliptical latitude (output)
35 !! \param rr Apparent geocentric distance (output)
36 !!
37 !! \note
38 !! - This is supposed to be the ELP2000-85 version, with the corrections from the 1998 paper (mean arguments up to t^4,
39 !! etc,). However, the accuracy seems to be less than promised for historical calculations (0.1° rather than 0.01°
40 !! for CE 0). The subroutine elp_mpp02_lbr() below is supposed to give better results (but then again, so is this
41 !! one).
42 !! - Differences compared to ELP-MPP02 (ELP82b minus ELP-MPP02):
43 !! - longitude: +0.00001° in CE 1975, ~+0.11° in CE 0/4000 and ~+0.65° in 3000 BCE (systematicly drifting to positive values)
44 !! - replacing w(1,2) = -5.8883 with -6.8084 removes the systematic drift to +0.65° @3000 BCE; (data.f90 / readmoondata())
45 !! the drift now oscillates around +-0.012°;
46 !! - latitude: ~0.00000° in CE 1990, ~+/-0.01° in CE 0/4000 and +/- ~0.06° (now ~0.006°) in 3000 BCE
47 !! - distance: +/-0.03 km in CE 2000, +/- ~30 km in CE 0/4000 and +/- ~300 km (now ~20km) in 3000 BCE
48 !!
49 !! \see
50 !! - Chapront-Touzé & Chapront, A&A, 124, 50 (1983)
51 !! - Chapront-Touzé & Chapront, A&A, 190, 342 (1988)
52 !! - ftp://cdsarc.u-strasbg.fr/pub/cats/VI/79/
53
54 subroutine elp82b_lbr(tjj, ll,bb,rr)
55 use sufr_kinds, only: double, dbl
56 use sufr_constants, only: as2r, au, km !, r2d
57 use sufr_angles, only: rev
58
60
61 implicit none
62 real(double), intent(in) :: tjj
63 real(double), intent(out) :: ll,bb,rr
64
65 integer :: iv,itab,nt,k,j
66 real(double) :: r(3),x,y, pa,t4,t8 !, lon,lat,dist
67
68
69 r = 0.0_dbl
70 x = 0.0_dbl
71 y = 0.0_dbl
72
73 t(1) = tjj*10.0_dbl ! In Julian centuries since 2000.0
74 t(2) = t(1)**2 ! t^2
75 t(3) = t(2)*t(1) ! t^3
76 t(4) = t(2)**2 ! t^4
77
78 t4 = t(4) ! t^4
79 t8 = t4**2 ! t^8
80
81
82 do iv=1,3 ! 3 variables (lon, lat, dist)
83 r(iv) = 0.0_dbl
84
85 do itab = 1,12 ! 12 tables: ELP1,4,7,10,13,16,19,22,25,28,31,34 for lon
86 do nt = 1,nterm(iv,itab)
87
88 select case(iv) ! variable
89
90 case(1) ! lon
91 if(itab.eq.1) then
92 x = pc1(1,nt)
93 y = pc1(2,nt)
94 do k = 1,4
95 y = y + pc1(k+2,nt)*t(k)
96 end do
97 else
98 j = nrang(1,itab-1) + nt
99 x = per1(1,j)
100 y = per1(2,j) + per1(3,j)*t(1)
101 end if
102
103 case(2) ! lat
104 if(itab.eq.1) then
105 x = pc2(1,nt)
106 y = pc2(2,nt)
107 do k = 1,4
108 y = y + pc2(k+2,nt)*t(k)
109 end do
110 else
111 j = nrang(2,itab-1) + nt
112 x = per2(1,j)
113 y = per2(2,j) + per2(3,j)*t(1)
114 end if
115
116 case(3) ! dist
117 if(itab.eq.1) then
118 x = pc3(1,nt)
119 y = pc3(2,nt)
120 do k = 1,4
121 y = y + pc3(k+2,nt)*t(k)
122 end do
123 else
124 j = nrang(3,itab-1) + nt
125 x = per3(1,j)
126 y = per3(2,j) + per3(3,j)*t(1)
127 end if
128 end select
129
130 if(itab.eq.3.or.itab.eq.5 .or. itab.eq.7.or.itab.eq.9) x = x*t(1)
131 if(itab.eq.12) x = x*t(2)
132 r(iv) = r(iv) + x*sin(y)
133 end do ! nt
134 end do ! itab
135 end do ! iv
136
137
138 ! Change of units:
139 r(1) = r(1) * as2r + w(1,0) + w(1,1)*t(1) + w(1,2)*t(2) + w(1,3)*t(3) + w(1,4)*t(4) ! Add mean longitude (see ELP PS-file, p.3)
140 r(2) = r(2) * as2r
141 r(3) = r(3) * a0 / ath
142
143 !lon = r(1)
144 !lat = r(2)
145 !dist = r(3)
146 !write(*,'(A, F15.9,2F14.7,9F14.5)') 'ELP82b: ', t(1), rev(lon)*r2d, lat*r2d, dist, &
147 ! rev(w(1,0) + w(1,1)*t(1) + w(1,2)*t(2) + w(1,3)*t(3) + w(1,4)*t(4))*r2d, &
148 ! rev(w(1,0))*r2d, rev(w(1,1)*t(1))*r2d, rev(w(1,2)*t(2))*r2d, rev(w(1,3)*t(3))*r2d, rev(w(1,4)*t(4))*r2d
149
150 ! Precess from J2000 to EoD, see ELP PS-file, p12; Laskar 1986 - note: longitude only!:
151 pa = 2.438174835301452d-2 * t(1) + 5.391128133938040d-6 * t(2) + 3.733065344543427d-10 * t(3) &
152 - 1.140766591650738d-10 * t4 - 8.753311012432670d-14 * t4*t(1) + 8.460483549042511d-16 * t4*t(2) &
153 + 6.348635154129373d-18 * t4*t(3) + 1.175188363009515e-20 * t8 - 2.307228308400282d-22 * t8*t(1) &
154 - 4.198486478408581d-25 * t8*t(2)
155
156 r(1) = r(1) + pa
157
158 ll = rev(r(1))
159 bb = r(2)
160 rr = r(3)/au*km
161
162 ! The original code converts to rectangular coordinates and applies (more?) precession.
163 ! Since I want spherical coordinates, I don't do that here. This gives similar (though not identical) results.
164
165 !write(98,*) ll,bb,rr
166
167 end subroutine elp82b_lbr
168 !*********************************************************************************************************************************
169
170
171
172 !*********************************************************************************************************************************
173 !> \brief Compute the spherical lunar coordinates using the ELP2000/MPP02 lunar theory in the dynamical mean ecliptic and
174 !! equinox of J2000.
175 !!
176 !! \param jd Julian day to compute Moon position for
177 !! \param mode Index of the corrections to the constants: 0-Fit to LLR observations, 1-Fit to DE405 1950-2060 (historical)
178 !!
179 !! \param lon Ecliptic longitude (rad) (output)
180 !! \param lat Ecliptic latitude (rad) (output)
181 !! \param rad Distance (AU) (output)
182 !!
183 !! \note See Lunar Solution ELP 2000/MPP02, Chapront & Francou (2002): ftp://cyrano-se.obspm.fr/pub/2_lunar_solutions/2_elpmpp02/
184
185 subroutine elp_mpp02_lbr(jd, mode, lon,lat,rad)
186 use sufr_kinds, only: double
187 use sufr_constants, only: jd2000, au, km !, r2d
188 use sufr_angles, only: rev
190
191 implicit none
192 real(double), intent(in) :: jd
193 integer, intent(in) :: mode
194 real(double), intent(out) :: lon,lat,rad
195 integer :: ierr
196 real(double) :: xyz(3),vxyz(3)
197
198 call elp_mpp02_xyz(jd, mode, xyz,vxyz, ierr)
199
200 ! Compute ecliptic l,b,r:
201 rad = sqrt(sum(xyz**2))
202 lon = atan2(xyz(2),xyz(1))
203 lat = asin(xyz(3)/rad)
204
205 call precess_ecl(jd2000,jd, lon,lat)
206
207 rad = rad*km/au ! km -> AU
208
209 end subroutine elp_mpp02_lbr
210 !*********************************************************************************************************************************
211
212
213 !***************************************************************************************************
214 !> \brief Compute the rectangular lunar coordinates using the ELP/MPP02 lunar theory in the dynamical mean ecliptic and equinox of J2000.
215 !!
216 !! \param jd Julian day to compute Moon position for
217 !! \param mode Index of the corrections to the constants: 0-Fit to LLR observations, 1-Fit to DE405 1950-2060 (historical)
218 !!
219 !! \param xyz Geocentric rectangular coordinates: (output)
220 !! - xyz(1) : Position X (km)
221 !! - xyz(2) : Position Y (km)
222 !! - xyz(3) : Position Z (km)
223 !! \param vxyz Geocentric rectangular velocities: (output)
224 !! - vxyz(1) : Velocity X' (km/day)
225 !! - vxyz(2) : Velocity Y' (km/day)
226 !! - vxyz(3) : Velocity Z' (km/day)
227 !! \param ierr File error index - ierr=0: no error, ierr=1: file error (output)
228 !!
229 !! \note
230 !! - The subroutine elp_mpp02() uses two modules:
231 !! - elp_mpp02_constants: Constants of the solution ELP/MPP02 (input),
232 !! - elp_mpp02_series: Series of the solution ELP/MPP02 (input).
233 !!
234 !! - The nominal values of some constants have to be corrected. There are two sets of corrections, which can be selected
235 !! using the parameter 'mode' (used in elp_mpp02_initialise()).
236 !! - mode=0, the constants are fitted to LLR observations provided from 1970 to 2001; it is the default value;
237 !! - mode=1, the constants are fitted to DE405 ephemeris over one century (1950-2060); the lunar angles W1, W2, W3
238 !! receive also additive corrections to the secular coefficients ('historical mode').
239 !! When the mode is changed, the constants will be reinitialised and the data file reread.
240 !!
241 !! - Solutions (discussed) in the paper:
242 !! - ELP (original):
243 !! - ELP2000-82: using VSOP82 (1983)
244 !! - ELP2000-85: new mean lunar arguments, higher truncation level, longer time range (1988)
245 !! - ELP2000-82B, here called "ELP": ELP2000-82, using mean lunar arguments from ELP2000-85 (19??)
246 !! - ELP/MPP01: using latest planetary perturbations from MPP01 and VSOP2000, but simpler than MPP01
247 !! - ELP/MPP02: ELP/MPP01, but for some arguments back to ELP + different selection of perturbations + lower truncation. Good fit with DE 405 in [1950,2060]
248 !! - ELP/MPP02*: improved secular arguments, better long-term comparison to DE 405/406 [-3000,2500]
249 !! - ELP/MPP02(LLR): ELP/MPP02(*?), optimised for lunar ranging since 1970
250 !! - ELPa: ELP + few Poisson terms (tested in the current study only?)
251 !! - ELPa*: ELPa + better secular arguments (as in ELP/MPP02*)
252 !! - It is not entirely clear which version is given below, but we can hope it is ELP/MPP02*. However, the subroutine
253 !! elp82b_lbr() above is known to underperform (by a factor of 10) in accuracy.
254
255 subroutine elp_mpp02_xyz(jd, mode, xyz,vxyz, ierr)
256 use sufr_kinds, only: double
257 use sufr_constants, only: r2as, jd2000 !, r2d!,as2r
258 use sufr_system, only: quit_program_error
259 use sufr_angles, only: rev
261
264
265 implicit none
266 integer, intent(in) :: mode
267 real(double), intent(in) :: jd
268 integer, intent(out) :: ierr
269 real(double), intent(out) :: xyz(3),vxyz(3)
270
271
272 real(double), parameter :: a405=384747.9613701725d0, aelp=384747.980674318d0, sc=36525.d0 ! Moon mean distance for DE405 und ELP; Julian century in days
273
274 integer :: it,iLine,iVar, k
275 real(double) :: rjd, t(-1:4),v(6) !, lon,lat,dist
276 real(double) :: cbeta,clamb,cw, ppw,ppw2,ppwqpw,ppwra,pw,pw2,pwqw,pwra, qpw,qpw2,qpwra,qw,qw2,qwra
277 real(double) :: ra,rap,sbeta,slamb,sw, x,x1,x2,x3, xp,xp1,xp2,xp3, y,yp
278
279 ! Initialise data and read files if needed:
281 if(ierr.ne.0) call quit_program_error('Could not read ELP-MPP02 files',0)
282
283
284 ! Initialization of time powers:
285 rjd = jd - jd2000 ! Reduced JD - JD since 2000
286 t(0) = 1.d0
287 t(1) = rjd/sc ! t: time since 2000 in Julian centuries
288 t(2) = t(1)**2 ! t^2
289 t(3) = t(2)*t(1) ! t^3
290 t(4) = t(2)**2 ! t^4
291
292 ! Evaluation of the series: substitution of time in the series
293 v = 0.d0
294 do ivar=1,3 ! iVar=1,2,3: Longitude, Latitude, Distance
295
296 ! Main Problem series:
297 do iline=nmpb(ivar,2),nmpb(ivar,3)
298 x = cmpb(iline)
299 y = fmpb(0,iline)
300 yp = 0.d0
301
302 do k=1,4
303 y = y + fmpb(k,iline)*t(k)
304 yp = yp + k*fmpb(k,iline)*t(k-1)
305 end do ! k
306
307 v(ivar) = v(ivar) + x*sin(y)
308 v(ivar+3) = v(ivar+3) + x*yp*cos(y)
309 end do ! iLine
310
311 ! Perturbations series:
312 do it=0,3
313 do iline=nper(ivar,it,2),nper(ivar,it,3)
314 x = cper(iline)
315 y = fper(0,iline)
316 xp = 0.d0
317 yp = 0.d0
318 if(it.ne.0) xp = it * x * t(it-1)
319
320 do k=1,4
321 y = y + fper(k,iline)*t(k)
322 yp = yp + k*fper(k,iline)*t(k-1)
323 end do ! k
324
325 v(ivar) = v(ivar) + x * t(it) * sin(y)
326 v(ivar+3) = v(ivar+3) + xp*sin(y) + x*t(it)*yp*cos(y)
327 end do ! iLine
328 end do ! it
329
330 end do ! iVar
331
332
333 ! Compute the spherical coordinates for the mean inertial ecliptic and equinox of date:
334 v(1) = v(1)/r2as + w(1,0) + w(1,1)*t(1) + w(1,2)*t(2) + w(1,3)*t(3) + w(1,4)*t(4) ! Longitude + mean longitude (rad)
335 !v(1) = rev(v(1)/r2as) + w(1,0) + rev(w(1,1)*t(1)) + rev(w(1,2)*t(2)) + rev(w(1,3)*t(3)) + rev(w(1,4)*t(4)) ! Longitude + mean longitude (rad)
336 v(2) = v(2)/r2as ! Latitude (rad)
337 v(3) = v(3) * a405 / aelp ! Distance (km)
338
339 !lon = v(1)
340 !lat = v(2)
341 !dist = v(3)
342 !lon = lon + (5029.0966d0*t(1) + 1.1120d0*t(2) + 0.000077d0*t(3) - 0.00002353d0*t(4) - 0.29965d0*t(1)) * as2r ! Precession from J2000 to EoD(?), but only in longitude!
343 !write(*,'(A, F15.9,2F14.7,9F14.5)') 'ELP-MPP02: ', t(1), rev(lon)*r2d, lat*r2d, dist, &
344 ! rev(w(1,0) + w(1,1)*t(1) + w(1,2)*t(2) + w(1,3)*t(3) + w(1,4)*t(4))*r2d, &
345 ! rev(w(1,0))*r2d, rev(w(1,1)*t(1))*r2d, rev(w(1,2)*t(2))*r2d, rev(w(1,3)*t(3))*r2d, rev(w(1,4)*t(4))*r2d
346
347 v(1) = rev(v(1)) ! This adds a bit of CPU time, but also alters the outcome of the cos/sin, similarly to
348 ! taking the 5 rev()s before adding up the terms when computing v(1), which might
349 ! indicate that this gives a better result, especially for dates far from 2000
350
351
352 ! Compute the rectangular coordinates (for the EoD?):
353 clamb = cos(v(1))
354 slamb = sin(v(1))
355 cbeta = cos(v(2))
356 sbeta = sin(v(2))
357 cw = v(3)*cbeta
358 sw = v(3)*sbeta
359
360 x1 = cw*clamb
361 x2 = cw*slamb
362 x3 = sw
363
364 ! Is this simply precession in rectangular coordinates from EoD to J2000?
365 pw = (p1 + p2*t(1) + p3*t(2) + p4*t(3) + p5*t(4)) * t(1)
366 qw = (q1 + q2*t(1) + q3*t(2) + q4*t(3) + q5*t(4)) * t(1)
367
368 ra = 2*sqrt(1.d0 - pw**2 - qw**2)
369 pwqw = 2*pw*qw
370 pw2 = 1.d0 - 2*pw**2
371 qw2 = 1.d0 - 2*qw**2
372 pwra = pw*ra
373 qwra = qw*ra
374
375 xyz(1) = pw2*x1 + pwqw*x2 + pwra*x3
376 xyz(2) = pwqw*x1 + qw2*x2 - qwra*x3
377 xyz(3) = -pwra*x1 + qwra*x2 + (pw2+qw2-1.d0)*x3
378
379 !xyz(1) = x1
380 !xyz(2) = x2
381 !xyz(3) = x3
382
383
384 ! Compute the rectangular velocities for the equinox J2000:
385 v(4) = v(4)/r2as + w(1,1) + 2*w(1,2)*t(1) + 3*w(1,3)*t(2) + 4*w(1,4)*t(3)
386 v(5) = v(5)/r2as
387
388 xp1 = (v(6)*cbeta - v(5)*sw)*clamb - v(4)*x2
389 xp2 = (v(6)*cbeta - v(5)*sw)*slamb + v(4)*x1
390 xp3 = v(6)*sbeta + v(5)*cw
391
392 ppw = p1 + (2*p2 + 3*p3*t(1) + 4*p4*t(2) + 5*p5*t(3)) * t(1)
393 qpw = q1 + (2*q2 + 3*q3*t(1) + 4*q4*t(2) + 5*q5*t(3)) * t(1)
394 ppw2 = -4*pw*ppw
395 qpw2 = -4*qw*qpw
396 ppwqpw = 2*(ppw*qw + pw*qpw)
397 rap = (ppw2+qpw2)/ra
398 ppwra = ppw*ra + pw*rap
399 qpwra = qpw*ra + qw*rap
400
401 vxyz(1) = (pw2*xp1 + pwqw*xp2 + pwra*xp3 + ppw2*x1 + ppwqpw*x2 + ppwra*x3) / sc
402 vxyz(2) = (pwqw*xp1 + qw2*xp2 - qwra*xp3 + ppwqpw*x1 + qpw2*x2 - qpwra*x3) / sc
403 vxyz(3) = (-pwra*xp1 + qwra*xp2 + (pw2+qw2-1.d0)*xp3 - ppwra*x1 + qpwra*x2 + (ppw2+qpw2)*x3) / sc
404
405 end subroutine elp_mpp02_xyz
406 !***************************************************************************************************
407
408
409
410
411
412 !*********************************************************************************************************************************
413 !> \brief Quick, lower-accuracy lunar coordinates; ~600x faster than ELP
414 !!
415 !! \param jd Julian day for computation
416 !! \param calc Calculate: 1: l,b,r, 2: & ra,dec, 3: & gmst,agst, 4: & az,alt, nt: number of terms <=60
417 !! \param nt Number of terms to use
418 !!
419 !! \see
420 !! - Meeus, Astronomical Algorithms, 1998, Ch. 22 and 47
421 !! - Simon et al, A&A, 282, p.663 (1994)
422
423 subroutine moonpos_la(jd, calc,nt)
424 use sufr_kinds, only: double
425 use sufr_constants, only: au, pi,d2r, jd2000, pland,earthr
426 use sufr_angles, only: rev, rev2
427 use sufr_system, only: quit_program_warning
428
429 use thesky_sun, only: sunpos_la
433 use thesky_local, only: lon0,lat0
434
435 implicit none
436 real(double), intent(in) :: jd
437 integer, intent(in) :: calc,nt
438
439 integer :: iLine,iArg, mal(4,60),mab(4,60)
440 real(double) :: jde,deltat, tjc,tjc2,tjc3,tjc4, l,b,r, lm,d,ms,mm,f,e,esl(60),esb(60),a1,a2,a3
441 real(double) :: ls,omg,ra,dec,eps,eps0,deps,gmst,agst,lst,dpsi,az,alt,hh, args(4),argl,argb
442 real(double) :: moondat(nPlanpos), hcl0,hcb0,hcr0, gcl,gcb,delta, elon,pa,illfr
443
444 if(sum(moonla_lrb(1:2,1)) .ne. -14616581) call quit_program_warning('moonpos_la(): moon_la.dat not read', 0)
445
446 deltat = calc_deltat(jd)
447 jde = jd + deltat/86400.d0
448
449 tjc = (jde-jd2000)/36525.d0 ! Julian Centuries after 2000.0 in dynamical time, the T in Meeus, p.163
450 tjc2 = tjc**2
451 tjc3 = tjc*tjc2
452 tjc4 = tjc2**2
453
454 ! Moon's mean longitude, Meeus p.338:
455 lm = rev(3.8103408236d0 + 8399.7091116339958d0*tjc - 2.755176757d-5*tjc2 + 3.239043d-8*tjc3 - 2.6771d-10*tjc4)
456
457 ! Delaunay arguments [d, ms, mm, f] (Meeus p.144) = [D, l', l, F] in Simon et al. 1994, Sect. 3.5:
458 ! d = rev(5.19846946025d0 + 7771.37714617d0*tjc - 3.340909d-5*tjc2 + 9.2114446d-8*tjc3) ! Moon's mean elongation
459 ! ms = rev(6.24003588115d0 + 628.301956024d0*tjc - 2.79776d-6*tjc2 - 5.8177641733d-8*tjc3) ! Sun's mean anomaly
460 ! mm = rev(2.3555483693d0 + 8328.69142288d0*tjc + 1.517947757d-4*tjc2 + 3.102807559d-7*tjc3) ! Moon's mean anomaly
461 ! f = rev(1.62790192912d0 + 8433.46615832d0*tjc - 6.42717497d-5*tjc2 + 5.3329949d-8*tjc3) ! Moon's argument of latitude
462
463 ! Meeus p.338, compares better to ELP-MPP02:
464 d = rev(5.1984665298d0 + 7771.377144834d0*tjc - 3.2845d-5*tjc2 + 3.197347d-8*tjc3 - 1.5436512d-10*tjc4) ! Moon's mean elongation
465 ms = rev(6.240060127d0 + 628.301955167d0*tjc - 2.681d-6*tjc2 + 7.1267017d-10*tjc3) ! Sun's mean anomaly
466 mm = rev(2.355555637d0 + 8328.691424759d0*tjc + 1.52566d-4*tjc2 + 2.5041d-7*tjc3 - 1.18633d-9*tjc4) ! Moon's mean anomaly
467 f = rev(1.627905158d0 + 8433.466158061d0*tjc - 6.3773d-5*tjc2 - 4.94988d-9*tjc3 + 2.02167d-11*tjc4) ! Moon's argument of latitude
468 args = [d,ms,mm,f] ! Delaunay arguments
469
470 e = 1.d0 - 0.002516d0*tjc - 0.0000074d0*tjc2
471
472 l = 0.d0
473 r = 0.d0
474 b = 0.d0
475
476 mal = moonla_arg(1:4,:)
477 mab = moonla_arg(5:8,:)
478 esl = e**abs(moonla_arg(2,:))
479 esb = e**abs(moonla_arg(6,:))
480
481 do iline=1,min(nt,60)
482 argl = 0.d0
483 argb = 0.d0
484 do iarg=1,4
485 argl = argl + mal(iarg,iline)*args(iarg)
486 argb = argb + mab(iarg,iline)*args(iarg)
487 end do
488 l = l + sin(argl) * moonla_lrb(1,iline) * esl(iline)
489 r = r + cos(argl) * moonla_lrb(2,iline) * esl(iline)
490 b = b + sin(argb) * moonla_lrb(3,iline) * esb(iline)
491 end do
492
493
494 ! Perturbations by other planets, and flattening of the Earth:
495 ! Meeus, p.338:
496 a1 = rev(2.090032d0 + 2.301199d0 * tjc) ! Influence from Venus
497 a2 = rev(0.926595d0 + 8364.7398477d0 * tjc) ! Influence from Jupiter
498 a3 = rev(5.4707345d0 + 8399.6847253d0 * tjc)
499
500 ! Meeus, p.342:
501 l = l + 3958*sin(a1) + 1962*sin(lm-f) + 318*sin(a2)
502 b = b - 2235*sin(lm) + 382*sin(a3) + 175*sin(a1-f) + 175*sin(a1+f) + 127*sin(lm-mm) - 115*sin(lm+mm)
503
504
505 ! Compute nutation:
506 omg = 2.18243858558d0 - 33.7570459367d0*tjc + 3.6142278d-5*tjc2 + 3.87850944888d-8*tjc3 ! Moon's mean lon. of asc.node, Meeus p.144
507 ls = 4.89506386655d0 + 62.84528862d0*tjc ! Mean long. Sun, Meeus p.144
508 dpsi = -8.338795d-5*sin(omg) - 6.39954d-6*sin(2*ls) - 1.115d-6*sin(2*lm) + 1.018d-6*sin(2*omg)
509
510 l = rev( dble(l) * 1.d-6*d2r + lm + dpsi)
511 b = rev2(dble(b) * 1.d-6*d2r)
512 r = (dble(r*100) + 3.8500056d10)/au
513
514
515 ! Store results:
516 planpos(1) = l ! Geocentric ecliptic longitude
517 planpos(2) = b ! Geocentric ecliptic latitude
518 planpos(4) = r ! Geocentric distance
519
520 planpos(7) = 5.77551830441d-3 * r ! Light time in days
521 planpos(12) = pland(0)/(r*au) ! Apparent diameter of the Moon
522
523 planpos(40) = jde ! JDE
524 planpos(46) = tjc ! App. dyn. time in Julian Centuries since 2000.0
525 planpos(47) = dpsi ! Nutation in longitude
526
527 if(calc.eq.1) return
528
529
530 ! Obliquity of the ecliptic:
531 eps0 = 0.409092804222d0 - 2.26965525d-4*tjc - 2.86d-9*tjc2 + 8.78967d-9*tjc3 ! Mean obliquity of the ecliptic (Lieske et al, 1977)
532 deps = 4.468d-5*cos(omg) + 2.76d-6*cos(2*ls) + 4.848d-7*cos(2*lm) - 4.36d-7*cos(2*omg) ! Nutation in obliquity
533 eps = eps0 + deps ! True obliquity of the ecliptic
534
535 call ecl_2_eq(l,b,eps, ra,dec) ! Ecliptical -> equatorial coordinates
536
537 ! Store results:
538 planpos(5) = ra ! Geocentric right ascension
539 planpos(6) = dec ! Geocentric declination
540
541 planpos(48) = eps ! True obliquity of the ecliptic; corrected for nutation
542 planpos(50) = eps0 ! Mean obliquity of the ecliptic; without nutation
543
544 if(calc.eq.2) return
545
546
547 ! Sidereal time:
548 gmst = calc_gmst(jd) ! Greenwich mean sidereal time
549 agst = rev(gmst + dpsi*cos(eps)) ! Correction for equation of the equinoxes -> Gr. apparent sid. time
550 lst = rev(agst + lon0) ! Local apparent sidereal time, lon0 > 0 for E
551
552 planpos(44) = lst ! Local APPARENT sidereal time
553 planpos(45) = agst ! Greenwich APPARENT sidereal time (in radians)
554 planpos(49) = gmst ! Greenwich MEAN sidereal time (in radians)
555
556 if(calc.eq.3) return
557
558
559 call eq2horiz(ra,dec,agst, hh,az,alt) ! Equatorial -> horizontal coordinates
560
561 planpos(8) = hh ! Geocentric hour angle
562 planpos(9) = az ! Geocentric azimuth
563 planpos(10) = alt ! Geocentric altitude
564
565 if(calc.eq.4) return
566
567
568 moondat = planpos ! Store Moon data
569
570 ! Get some Sun data:
571 call sunpos_la(jd,1)
572 hcl0 = rev(planpos(1)+pi) ! Heliocentric longitude of the Earth - want geocentric lon of Sun?
573 hcb0 = -planpos(2) ! Heliocentric latitude of the Earth - want geocentric lat of Sun?
574 hcr0 = planpos(3)
575
576 planpos = moondat ! Restore Moon data:
577 gcl = planpos(1) ! l
578 gcb = planpos(2) ! b
579 delta = planpos(4) ! r
580
581
582 elon = acos(cos(gcb)*cos(hcb0)*cos(gcl-rev(hcl0+pi)))
583 pa = atan2( hcr0*sin(elon) , delta - hcr0*cos(elon) ) ! Phase angle
584 illfr = 0.5d0*(1.d0 + cos(pa)) ! Illuminated fraction of the disc
585
586 planpos(11) = rev(elon) ! Elongation
587 planpos(13) = moonmagn(pa,delta) ! Magnitude
588 planpos(14) = illfr ! Illuminated fraction
589 planpos(15) = rev(pa) ! Phase angle
590 planpos(16) = atan2(sin(planpos(8)),tan(lat0)*cos(planpos(6)) - sin(planpos(6))*cos(planpos(8))) ! Parallactic angle
591 planpos(17) = asin(earthr/(planpos(4)*au)) ! Horizontal parallax
592
593 planpos(41) = rev(hcl0+pi) ! Geocentric, true L,B,R for the Sun
594 planpos(42) = rev2(-hcb0)
595 planpos(43) = hcr0
596
597 end subroutine moonpos_la
598 !*********************************************************************************************************************************
599
600
601
602
603
604
605
606 !*********************************************************************************************************************************
607 !> \brief Calculate the magnitude of the Moon
608 !!
609 !! \param pa Phase angle (rad)
610 !! \param delta Geocentric(?) distance (AU)
611 !!
612 !! \see
613 !! - Allen, 1976, par.66
614 !! - http://cleardarksky.com/others/BenSugerman/star.htm
615
616 function moonmagn(pa, delta)
617 use sufr_kinds, only: double
618
619 implicit none
620 real(double), intent(in) :: pa,delta
621 real(double) :: moonmagn
622
623 moonmagn = -12.73d0 + 1.489690267d0*abs(pa) + 0.04310727d0*pa**4 ! Allen, 1976, par.66 (d2r)
624
625 ! Correct for variable distance and the 'opposition effect' or 'opposition surge' which occurs when pa < 7d:
626 ! Not sure what the original reference is. See e.g.:
627 ! - Whitaker, E. A. 1969, in Analysis of Apollo 8 Photography and Visual Observations, NASA SP-201, 38, p.39
628 ! - Schaefer, PASP 102 212 (1990)
629 ! - Krisciunas & Schaefer, PASP 103, 1033 (1991): ""
630 moonmagn = moonmagn - 2.5*log10( (2.5696d-3/delta)**2 * max(1.d0, 1.35d0 - 2.864789d0*abs(pa) ) ) ! (1-1.35)/2.865 ~ 7deg
631
632 end function moonmagn
633 !*********************************************************************************************************************************
634
635
636
637end module thesky_moon
638!***********************************************************************************************************************************
Procedures for coordinates.
subroutine eq2horiz(ra, dec, agst, hh, az, alt, lat, lon)
Convert spherical equatorial coordinates (RA, dec, agst) to spherical horizontal coordinates (hh,...
subroutine ecl_2_eq(l, b, eps, ra, dec)
Convert (geocentric) spherical ecliptical coordinates l,b (and eps) to spherical equatorial coordinat...
subroutine precess_ecl(jd1, jd2, l, b)
Compute the precession of the equinoxes in geocentric ecliptical coordinates.
Procedures to set constants and read data files.
Definition data.f90:22
subroutine elp_mpp02_initialise_and_read_files(mode, ierr)
Definition data.f90:730
Date and time procedures.
Definition date_time.f90:48
real(double) function calc_gmst(jd, deltat)
Calculate Greenwich Mean Sidereal Time for any instant, in radians.
real(double) function calc_deltat(jd, force_recompute)
Compute DeltaT for a given JD.
real(double), dimension(3, 0:4) w
Definition modules.f90:195
real(double), dimension(max1) cmpb
Definition modules.f90:211
real(double), dimension(0:4, max1) fmpb
Definition modules.f90:211
real(double), dimension(max2) cper
Definition modules.f90:211
integer, dimension(3, 3) nmpb
Definition modules.f90:210
real(double), dimension(0:4, max2) fper
Definition modules.f90:211
integer, dimension(3, 0:3, 3) nper
Definition modules.f90:210
Local parameters for libTheSky: location, date, time.
Definition modules.f90:44
real(double) lon0
Definition modules.f90:51
real(double) lat0
Definition modules.f90:51
Procedures for the Moon.
subroutine elp82b_lbr(tjj, ll, bb, rr)
Calculate the apparent geocentric ecliptical position of the Moon, using the Lunar Solution ELP 2000-...
real(double) function moonmagn(pa, delta)
Calculate the magnitude of the Moon.
subroutine elp_mpp02_xyz(jd, mode, xyz, vxyz, ierr)
Compute the rectangular lunar coordinates using the ELP/MPP02 lunar theory in the dynamical mean ecli...
subroutine elp_mpp02_lbr(jd, mode, lon, lat, rad)
Compute the spherical lunar coordinates using the ELP2000/MPP02 lunar theory in the dynamical mean ec...
subroutine moonpos_la(jd, calc, nt)
Quick, lower-accuracy lunar coordinates; ~600x faster than ELP.
ELP-82B Moon data, needed to compute Moon positions.
Definition modules.f90:164
real(double), parameter a0
Definition modules.f90:171
real(double), dimension(3, 0:4) w
Definition modules.f90:174
integer, dimension(3, 0:12) nrang
Definition modules.f90:177
real(double), dimension(6, 704) pc3
Definition modules.f90:175
real(double), dimension(6, 918) pc2
Definition modules.f90:175
real(double), dimension(6, 1023) pc1
Definition modules.f90:175
real(double), dimension(3, 19537) per1
Definition modules.f90:175
real(double), parameter ath
Definition modules.f90:171
real(double), dimension(0:4) t
Definition modules.f90:174
integer, dimension(3, 12) nterm
Definition modules.f90:177
real(double), dimension(3, 6766) per2
Definition modules.f90:175
real(double), dimension(3, 8924) per3
Definition modules.f90:175
Planet data, needed to compute planet positions.
Definition modules.f90:61
integer(long), dimension(3, 60) moonla_lrb
Definition modules.f90:70
integer, parameter nplanpos
Definition modules.f90:67
real(double), dimension(nplanpos) planpos
Planpos[] is an array with many different types of coordinates and related variables:
Definition modules.f90:153
integer, dimension(8, 60) moonla_arg
Definition modules.f90:72
Low-accuracy procedures for the Sun.
Definition sun.f90:22
subroutine sunpos_la(jd, calc, lat, lon)
Low-accuracy solar coordinates.
Definition sun.f90:57