libTheSky
Routines to compute sky positions of Sun, Moon, planets and more
Loading...
Searching...
No Matches
coordinates.f90
Go to the documentation of this file.
1!> \file coordinates.f90 Procedures to perform coordinate transformations, apply precession, and more 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 coordinates
21
23 implicit none
24 save
25
26
27contains
28
29
30 !*********************************************************************************************************************************
31 !> \brief Compute the geocentric rectangular coordinates of a planet, from its and the Earth's heliocentric spherical position
32 !!
33 !! \param l Heliocentric longitude of planet
34 !! \param b Heliocentric latitude of planet
35 !! \param r Heliocentric distance of planet
36 !!
37 !! \param l0 Heliocentric longitude of Earth
38 !! \param b0 Heliocentric latitude of Earth
39 !! \param r0 Heliocentric distance of Earth
40 !!
41 !! \param x Geocentric rectangular X of planet (output)
42 !! \param y Geocentric rectangular Y of planet (output)
43 !! \param z Geocentric rectangular Z of planet (output)
44
45 subroutine hc_spher_2_gc_rect(l,b,r, l0,b0,r0, x,y,z)
46 use sufr_kinds, only: double
47
48 implicit none
49 real(double), intent(in) :: l,b,r, l0,b0,r0
50 real(double), intent(out) :: x,y,z
51
52 x = r * cos(b) * cos(l) - r0 * cos(b0) * cos(l0)
53 y = r * cos(b) * sin(l) - r0 * cos(b0) * sin(l0)
54 z = r * sin(b) - r0 * sin(b0)
55
56 end subroutine hc_spher_2_gc_rect
57 !*********************************************************************************************************************************
58
59
60 !*********************************************************************************************************************************
61 !> \brief Convert spherical, ecliptical coordinates to rectangular, equatorial coordinates of an object - both geocentric
62 !!
63 !! \param l Heliocentric longitude of planet
64 !! \param b Heliocentric latitude of planet
65 !! \param r Heliocentric distance of planet
66 !! \param eps Obliquity of the ecliptic
67 !!
68 !! \param x Geocentric rectangular X of planet (output)
69 !! \param y Geocentric rectangular Y of planet (output)
70 !! \param z Geocentric rectangular Z of planet (output)
71
72 subroutine ecl_spher_2_eq_rect(l,b,r, eps, x,y,z)
73 use sufr_kinds, only: double
74
75 implicit none
76 real(double), intent(in) :: l,b,r,eps
77 real(double), intent(out) :: x,y,z
78
79 x = r * cos(b) * cos(l)
80 y = r * (cos(b) * sin(l) * cos(eps) - sin(b) * sin(eps))
81 z = r * (cos(b) * sin(l) * sin(eps) + sin(b) * cos(eps))
82
83 end subroutine ecl_spher_2_eq_rect
84 !*********************************************************************************************************************************
85
86
87
88 !*********************************************************************************************************************************
89 !> \brief Compute the geocentric equatorial rectangular coordinates of the Sun, from Earth's heliocentric, spherical position
90 !!
91 !! \param t1 Dynamical time in Julian Millennia after 2000.0
92 !! \param l0 Heliocentric longitude of Earth
93 !! \param b0 Heliocentric latitude of Earth
94 !! \param r0 Heliocentric distance of Earth
95 !!
96 !! \param x Geocentric rectangular X of planet (output)
97 !! \param y Geocentric rectangular Y of planet (output)
98 !! \param z Geocentric rectangular Z of planet (output)
99 !!
100 !! \todo To be disposed, has been replcaed by ecl_spher_2_eq_rect above
101
102 subroutine calcsunxyz(t1, l0,b0,r0, x,y,z)
103 use sufr_kinds, only: double
104 use sufr_constants, only: pi
105 use sufr_angles, only: rev
106
107 use thesky_nutation, only: nutation
108
109 implicit none
110 real(double), intent(in) :: t1, l0,b0,r0
111 real(double), intent(out) :: x,y,z
112 real(double) :: t, l,b,r, dpsi,eps,eps0,deps
113
114 t = t1
115 l = rev(l0+pi) ! Heliocentric longitude of the Earth -> geocentric longitude of the Sun
116 b = -b0 ! Heliocentric latitude of the Earth -> geocentric latitude of the Sun
117 r = r0 ! Heliocentric distance of the Earth = geocentric distance of the Sun
118
119 call fk5(t, l,b)
120 call nutation(t,dpsi,eps0,deps)
121 eps = eps0
122
123 x = r * cos(b) * cos(l)
124 y = r * (cos(b) * sin(l) * cos(eps) - sin(b) * sin(eps))
125 z = r * (cos(b) * sin(l) * sin(eps) + sin(b) * cos(eps))
126
127 !jd1 = t1*365250.d0 + jd2000 !From equinox of date...
128 !jd2 = comepoche ! ... to equinox of elements
129 !call precess_xyz(jd1,jd2,x,y,z)
130
131 end subroutine calcsunxyz
132 !*********************************************************************************************************************************
133
134
135 !*********************************************************************************************************************************
136 !> \brief Compute the precession of the equinoxes in rectangular coordinates, from jd1 to jd2
137 !!
138 !! \param jd1 Original Julian day
139 !! \param jd2 Target Julian day
140 !!
141 !! \param x Geocentric rectangular X (output)
142 !! \param y Geocentric rectangular Y (output)
143 !! \param z Geocentric rectangular Z (output)
144 !!
145 !! \see Meeus, Astronomical Algorithms, 1998, Ch. 21
146
147 subroutine precess_xyz(jd1,jd2, x,y,z)
148 use sufr_kinds, only: double
149 use sufr_constants, only: jd2000
150
151 implicit none
152 real(double), intent(in) :: jd1,jd2
153 real(double), intent(inout) :: x,y,z
154 real(double) :: t1,t2,t22,t23, x1,y1,z1
155 real(double) :: dz,ze,th,sd,cd,sz,cz,st,ct
156 real(double) :: xx,xy,xz,yx,yy,yz,zx,zy,zz
157
158 t1 = (jd1 - jd2000)/36525.d0 ! t since 2000.0 in Julian centuries
159 t2 = (jd2 - jd1)/36525.d0 ! dt in Julian centuries
160 t22 = t2*t2 ! t2^2
161 t23 = t22*t2 ! t2^3
162
163 ! Meeus, Eq. 21.2:
164 dz = (1.11808609d-2 + 6.770714d-6*t1 - 6.739d-10*t1*t1)*t2 + (1.463556d-6 - 1.668d-9*t1)*t22 + 8.725677d-8*t23
165 ze = (1.11808609d-2 + 6.770714d-6*t1 - 6.739d-10*t1*t1)*t2 + (5.307158d-6 + 3.20d-10*t1)*t22 + 8.825063d-8*t23
166 th = (9.71717346d-3 - 4.136915d-6*t1 - 1.0520d-9*t1*t1)*t2 - (2.068458d-6 + 1.052d-9*t1)*t22 - 2.028121d-7*t23
167
168 sd = sin(dz)
169 cd = cos(dz)
170 sz = sin(ze)
171 cz = cos(ze)
172 st = sin(th)
173 ct = cos(th)
174
175 ! Ref?
176 xx = cd*cz*ct - sd*sz
177 xy = sd*cz + cd*sz*ct
178 xz = cd*st
179 yx = -cd*sz - sd*cz*ct
180 yy = cd*cz - sd*sz*ct
181 yz = -sd*st
182 zx = -cz*st
183 zy = -sz*st
184 zz = ct
185
186 x1 = xx*x + yx*y + zx*z
187 y1 = xy*x + yy*y + zy*z
188 z1 = xz*x + yz*y + zz*z
189
190 x = x1
191 y = y1
192 z = z1
193
194 end subroutine precess_xyz
195 !*********************************************************************************************************************************
196
197
198
199 !*********************************************************************************************************************************
200 !> \brief Compute the precession of the equinoxes in equatorial coordinates, from jd1 to jd2
201 !!
202 !! \param jd1 Original Julian day
203 !! \param jd2 Target Julian day
204 !! \param a1 Right ascension (IO, rad) (output)
205 !! \param d1 Declination (IO, rad) (output)
206 !!
207 !! \see Meeus, Astronomical Algorithms, 1998, Ch.21, p.134
208
209 subroutine precess_eq(jd1,jd2, a1,d1)
210 use sufr_kinds, only: double
211 use sufr_constants, only: jd2000
212
213 implicit none
214 real(double), intent(in) :: jd1,jd2
215 real(double), intent(inout) :: a1,d1
216 real(double) :: t1,t12, t2,t22,t23, dz,ze,th,a,b,c
217
218 t1 = (jd1 - jd2000)/36525.d0 ! t since 2000.0 in Julian centuries
219 t12 = t1**2 ! t1^2
220 t2 = (jd2 - jd1)/36525.d0 ! dt in Julian centuries
221 t22 = t2**2 ! t2^2
222 t23 = t22*t2 ! t2^3
223
224 ! Meeus, Eq. 21.2:
225 dz = (1.11808609d-2 + 6.770714d-6*t1 - 6.739d-10*t12)*t2 + (1.463556d-6 - 1.668d-9*t1)*t22 + 8.725677d-8*t23
226 ze = (1.11808609d-2 + 6.770714d-6*t1 - 6.739d-10*t12)*t2 + (5.307158d-6 + 3.20d-10*t1)*t22 + 8.825063d-8*t23
227 th = (9.71717346d-3 - 4.136915d-6*t1 - 1.0520d-9*t12)*t2 - (2.068458d-6 + 1.052d-9*t1)*t22 - 2.028121d-7*t23
228
229 ! Meeus, Eq. 21.4:
230 a = cos(d1) * sin(a1+dz)
231 b = cos(th) * cos(d1) * cos(a1+dz) - sin(th) * sin(d1)
232 c = sin(th) * cos(d1) * cos(a1+dz) + cos(th) * sin(d1)
233
234 a1 = atan2(a,b) + ze
235 d1 = asin(c)
236
237 end subroutine precess_eq
238 !*********************************************************************************************************************************
239
240
241
242 !*********************************************************************************************************************************
243 !> \brief Compute the precession of the equinoxes in equatorial coordinates, from yr1 to yr2
244 !!
245 !! \param yr1 Original year
246 !! \param yr2 Target year
247 !! \param a1 Right ascension (IO, rad) (output)
248 !! \param d1 Declination (IO, rad) (output)
249 !!
250 !! \see Meeus, Astronomical Algorithms, 1998, Ch.21, p.134
251
252 subroutine precess_eq_yr(yr1,yr2,a1,d1)
253 use sufr_kinds, only: double
254
255 implicit none
256 real(double), intent(in) :: yr1,yr2
257 real(double), intent(inout) :: a1,d1
258 real(double) :: t1,t12, t2,t22,t23, dz,ze,th,a,b,c
259
260 t1 = (yr1 - 2000.d0)/100.d0 ! t since 2000.0 in Julian centuries
261 t12 = t1**2 ! t1^2
262 t2 = (yr2 - yr1)/100.d0 ! dt in Julian centuries
263 t22 = t2*t2 ! t2^2
264 t23 = t22*t2 ! t2^3
265
266 ! Meeus, Eq. 21.2:
267 dz = (1.11808609d-2 + 6.770714d-6*t1 - 6.739d-10*t12)*t2 + (1.463556d-6 - 1.668d-9*t1)*t22 + 8.725677d-8*t23
268 ze = (1.11808609d-2 + 6.770714d-6*t1 - 6.739d-10*t12)*t2 + (5.307158d-6 + 3.20d-10*t1)*t22 + 8.825063d-8*t23
269 th = (9.71717346d-3 - 4.136915d-6*t1 - 1.0520d-9*t12)*t2 - (2.068458d-6 + 1.052d-9*t1)*t22 - 2.028121d-7*t23
270
271 ! Meeus, Eq. 21.4:
272 a = cos(d1) * sin(a1+dz)
273 b = cos(th) * cos(d1) * cos(a1+dz) - sin(th) * sin(d1)
274 c = sin(th) * cos(d1) * cos(a1+dz) + cos(th) * sin(d1)
275
276 a1 = atan2(a,b) + ze
277 d1 = asin(c)
278
279 end subroutine precess_eq_yr
280 !*********************************************************************************************************************************
281
282
283
284 !*********************************************************************************************************************************
285 !> \brief Compute the precession of the equinoxes in geocentric ecliptical coordinates
286 !!
287 !! \param jd1 Original Julian day
288 !! \param jd2 Target Julian day
289 !! \param l Ecliptic longitude (I/O, rad) (output)
290 !! \param b Ecliptic latitude (I/O, rad) (output)
291 !!
292 !! \see Meeus, Astronomical Algorithms, 1998, Ch.21, p.136
293
294 subroutine precess_ecl(jd1,jd2, l,b)
295 use sufr_kinds, only: double
296 use sufr_constants, only: jd2000
297
298 implicit none
299 real(double), intent(in) :: jd1,jd2
300 real(double), intent(inout) :: l,b
301 real(double) :: t1,t12, t2,t22,t23, eta,pii,p,aa,bb,cc
302
303 t1 = (jd1 - jd2000)/36525.d0 ! t since 2000.0 in Julian centuries
304 t12 = t1**2 ! t1^2
305 t2 = (jd2 - jd1)/36525.d0 ! dt in Julian centuries
306 t22 = t2*t2 ! t2^2
307 t23 = t22*t2 ! t2^3
308
309 ! Meeus, Eq. 21.5 - the different powers of t1/t2 have been checked and are ok:
310 eta = (2.278765d-4 - 3.2012d-7 *t1 + 2.899d-9*t12) *t2 + (-1.60085d-7 + 2.899d-9*t1) *t22 + 2.909d-10*t23
311 pii = 3.0521686858d0 + 1.59478437d-2 *t1 + 2.939037d-6*t12 - (4.2169525d-3 + 2.44787d-6*t1) *t2 + 1.7143d-7*t22
312 p = (2.438174835d-2 + 1.077382d-5 *t1 - 2.036d-10*t12) *t2 + (5.38691d-6 - 2.036d-10*t1) *t22 - 2.91d-11*t23
313
314 ! Meeus, Eq. 21.7:
315 aa = cos(eta) * cos(b) * sin(pii-l) - sin(eta) * sin(b)
316 bb = cos(b) * cos(pii-l)
317 cc = cos(eta) * sin(b) + sin(eta) * cos(b) * sin(pii-l)
318
319 l = p + pii - atan2(aa,bb)
320 b = asin(cc)
321
322 end subroutine precess_ecl
323 !*********************************************************************************************************************************
324
325
326
327 !*********************************************************************************************************************************
328 !> \brief Compute the precession of the equinoxes in orbital elements
329 !!
330 !! \param jd1 Original Julian day
331 !! \param jd2 Target Julian day
332 !!
333 !! \param i Inclination (output)
334 !! \param o1 Argument of perihelion (output)
335 !! \param o2 Longitude of ascending node (output)
336 !!
337 !! \see Meeus, Astronomical Algorithms, 1998, Ch.24
338
339
340 subroutine precess_orb(jd1,jd2,i,o1,o2)
341 use sufr_kinds, only: double
342 use sufr_constants, only: jd2000
343
344 implicit none
345 real(double), intent(in) :: jd1,jd2
346 real(double), intent(inout) :: i,o1,o2
347 real(double) :: t1,t12,t2,t22,t23, eta,pii,p,psi,aa,bb,cc,dd
348
349 t1 = (jd1 - jd2000)/36525.d0 ! t since 2000.0 in Julian centuries
350 t12 = t1**2 ! t1^2
351 t2 = (jd2 - jd1)/36525.d0 ! dt in Julian centuries
352 t22 = t2*t2 ! t2^2
353 t23 = t22*t2 ! t2^3
354
355 ! Meeus, Eq. 21.5:
356 eta = (2.278765d-4 - 3.2012d-7 * t1 + 2.899d-9*t12) *t2 + (-1.60085d-7 + 2.899d-9*t1) * t22 + 2.909d-10 * t23
357 pii = 3.0521686858d0 + 1.59478437d-2 * t1 + 2.939037d-6*t12 - (4.2169525d-3 + 2.44787d-6*t1) * t2 + 1.7143d-7 * t22
358 p = (2.438174835d-2 + 1.077382d-5 * t1 - 2.036d-10*t12) *t2 + (5.38691d-6 - 2.036d-10*t1) * t22 - 2.91d-11 * t23
359 psi = pii + p
360
361 ! Meeus, Eq. 24.2:
362 aa = sin(i) * sin(o2-pii)
363 bb = -sin(eta) * cos(i) + cos(eta) * sin(i) * cos(o2-pii)
364 cc = -sin(eta) * sin(o2-pii)
365 dd = sin(i) * cos(eta) - cos(i) * sin(eta) * cos(o2-pii)
366
367 i = asin(sqrt(aa*aa+bb*bb)) ! Inclination
368 o2 = atan2(aa,bb) + psi ! Longitude of ascending node (Omega)
369 o1 = o1 + atan2(cc,dd) ! Argument of perihelion (omega)
370
371 end subroutine precess_orb
372 !*********************************************************************************************************************************
373
374
375
376 !*********************************************************************************************************************************
377 !> \brief Convert rectangular coordinates x,y,z to spherical coordinates l,b,r
378 !!
379 !! \param x Rectangular x coordinate (same unit as r)
380 !! \param y Rectangular y coordinate (same unit as r)
381 !! \param z Rectangular z coordinate (same unit as r)
382 !!
383 !! \param l Longitude (rad) (output)
384 !! \param b Latitude (rad) (output)
385 !! \param r Distance (same unit as x,y,z) (output)
386
387 subroutine rect_2_spher(x,y,z, l,b,r)
388 use sufr_kinds, only: double
389 use sufr_angles, only: rev
390 use sufr_numerics, only: deq0
391
392 implicit none
393 real(double), intent(in) :: x,y,z
394 real(double), intent(out) :: l,b,r
395 real(double) :: x2,y2
396
397 if(deq0(x) .and. deq0(y) .and. deq0(z)) then
398 l = 0.d0
399 b = 0.d0
400 r = 0.d0
401 else
402 x2 = x**2
403 y2 = y**2
404
405 l = rev( atan2(y,x) ) ! Longitude
406 b = atan2(z, sqrt(x2 + y2)) ! Latitude
407 r = sqrt(x2 + y2 + z**2) ! Distance
408 end if
409
410 end subroutine rect_2_spher
411 !*********************************************************************************************************************************
412
413
414 !*********************************************************************************************************************************
415 !> \brief Correct ecliptic longitude and latitiude for annual aberration
416 !!
417 !! \param t Dynamical time in Julian Millennia since 2000.0
418 !! \param l0 Earth longitude (rad)
419 !!
420 !! \param l Longitude of the object (rad, I/O) (output)
421 !! \param b Latitude of the object (rad, I/O) (output)
422 !!
423 !! \see Meeus, Astronomical Algorithms, 1998, Ch.23
424
425 subroutine aberration_ecl(t,l0, l,b)
426 use sufr_kinds, only: double
427 use sufr_constants, only: pi
428 use sufr_angles, only: rev
429
430 implicit none
431 real(double), intent(in) :: t,l0
432 real(double), intent(inout) :: l,b
433 real(double) :: tt,tt2, k,odot,e,pii,dl,db
434
435 tt = t*10 ! Time in Julian Centuries
436 tt2 = tt*tt
437
438 k = 9.93650849745d-5 ! Constant of aberration, radians
439 odot = rev(l0+pi) ! Longitude of the Sun
440
441 e = 0.016708634d0 - 4.2037d-5 * tt - 1.267d-7 * tt2 ! Eccentricity of the Earth's orbit
442 pii = 1.79659568d0 + 3.001024d-2 * tt + 8.0285d-6 * tt2 ! Longitude of perihelion of "
443
444 dl = k * ( -cos(odot-l) + e * cos(pii-l)) / cos(b) ! Meeus, Eq. 23.2
445 db = -k * sin(b) * (sin(odot-l) - e * sin(pii-l))
446
447 l = l + dl
448 b = b + db
449
450 end subroutine aberration_ecl
451 !*********************************************************************************************************************************
452
453 !*********************************************************************************************************************************
454 !> \brief Correct equatorial coordinates for annual aberration - moderate accuracy, use for stars
455 !!
456 !! \param jd Julian day of epoch
457 !!
458 !! \param ra Right ascension of the object (rad)
459 !! \param dec Declination of the object (rad)
460 !!
461 !! \param dra Correction in right ascension due to aberration (rad) (output)
462 !! \param ddec Correction in declination due to aberration(rad) (output)
463 !!
464 !! \param eps0 The mean obliquity of the ecliptic
465 !!
466 !!
467 !! \see Meeus, Astronomical Algorithms, 1998, Ch.23
468
469 subroutine aberration_eq(jd, ra,dec, dra, ddec, eps0)
470 use sufr_kinds, only: double
471 use sufr_constants, only: jd2000
472 use sufr_angles, only: rev
473 use sufr_dummy, only: dumdbl1,dumdbl2
474 use sufr_numerics, only: dne
475 use thesky_nutation, only: nutation
476
477 implicit none
478 real(double), intent(in) :: jd, ra,dec
479 real(double), intent(out) :: dra,ddec
480 real(double), intent(in), optional :: eps0
481 real(double) :: jdold, tjc,tjc2, l0,mas,sec,odot, k,ee,pii, eps
482 save jdold, tjc,odot,k,ee,pii, eps
483
484
485 ! Compute these constants (independent of ra,dec) only if jd has changed since the last call
486 if(dne(jd,jdold)) then
487
488 ! Meeus, Ch. 25:
489 tjc = (jd-jd2000)/36525.0d0 ! Julian centuries since 2000.0, Meeus Eq. 25.1
490 tjc2 = tjc**2
491 l0 = 4.8950631684d0 + 628.331966786d0 * tjc + 5.291838d-6 *tjc2 ! Mean longitude of the Sun, Meeus Eq. 25.2
492
493 ! Mean anomaly of the Sun, Meeus Eq. 25.3 -> p.144 (nutation):
494 mas = 6.240035881d0 + 628.301956024d0 * tjc - 2.79776d-6 * tjc2 - 5.817764d-8 * tjc**3
495
496 ! Sun's equation of the centre, Meeus p.164:
497 sec = (3.3416109d-2 - 8.40725d-5*tjc - 2.443d-7*tjc2)*sin(mas) + (3.489437d-4-1.763d-6*tjc)*sin(2*mas) + 5.044d-6*sin(3*mas)
498 odot = rev(l0 + sec) ! Sun's true longitude
499
500 ! Meeus, p.151:
501 k = 9.93650849745d-5 ! Constant of aberration, in radians (kappa in Meeus)
502 ee = 0.016708634d0 - 4.2037d-5 * tjc - 1.267d-7 * tjc2 ! Eccentricity of the Earth's orbit
503 pii = 1.79659568d0 + 3.001024d-2 * tjc + 8.0285d-6 * tjc2 ! Longitude of perihelion of "
504
505
506 ! If eps0 is provided, use it - otherwise compute it:
507 if(present(eps0)) then
508 eps = eps0
509 else
510 call nutation(tjc/10.d0, dumdbl1, eps, dumdbl2) ! Note: eps is actually eps0
511 end if
512
513 end if
514 jdold = jd
515
516
517 ! Note that eps in Meeus is actually eps0, the *mean* obliquity of the ecliptic
518 ! Meeus, Eq. 23.3:
519 dra = k * ( -(cos(ra) * cos(odot) * cos(eps) + sin(ra) * sin(odot)) / cos(dec) &
520 + ee * (cos(ra) * cos(pii) * cos(eps) + sin(ra) * sin(pii)) / cos(dec) )
521
522 ddec = k * ( -(cos(odot) * cos(eps) * (tan(eps) * cos(dec) - sin(ra) * sin(dec)) + cos(ra) * sin(dec) * sin(odot)) &
523 + ee * (cos(pii) * cos(eps) * (tan(eps) * cos(dec) - sin(ra) * sin(dec)) + cos(ra) * sin(dec) * sin(pii)) )
524
525 end subroutine aberration_eq
526 !*********************************************************************************************************************************
527
528
529
530 !*********************************************************************************************************************************
531 !> \brief Convert coordinates to the FK5 system
532 !!
533 !! \param t Dynamical time in Julian Millennia since 2000.0
534 !! \param l Longitude (rad, I/O) (output)
535 !! \param b Latitude (rad, I/O) (output)
536
537 subroutine fk5(t, l,b)
538 use sufr_kinds, only: double
539
540 implicit none
541 real(double), intent(in) :: t
542 real(double), intent(inout) :: l,b
543 real(double) :: tt,l2,dl,db
544
545 tt = t*10 ! Julian Centuries
546
547 l2 = l - 0.02438225d0*tt - 5.41052d-6*tt*tt
548 dl = -4.379322d-7 + 1.89853d-7*(cos(l2) + sin(l2))*tan(b)
549 db = 1.89853d-7*(cos(l2) - sin(l2))
550
551 l = l + dl
552 b = b + db
553
554 end subroutine fk5
555 !*********************************************************************************************************************************
556
557
558 !*********************************************************************************************************************************
559 !> \brief Convert (geocentric) spherical ecliptical coordinates l,b (and eps) to spherical equatorial coordinates RA, Dec
560 !!
561 !! \param l Longitude (rad)
562 !! \param b Latitude (rad)
563 !! \param eps Obliquity of the ecliptic
564 !!
565 !! \param ra Right ascension (rad) (output)
566 !! \param dec Declination (rad) (output)
567 !!
568 !! \see Expl. Supl. t.t. Astronimical Almanac 3rd Ed, Eq.14.43
569
570 subroutine ecl_2_eq(l,b,eps, ra,dec)
571 use sufr_kinds, only: double
572 use sufr_angles, only: rev
573
574 implicit none
575 real(double), intent(in) :: l,b,eps
576 real(double), intent(out) :: ra,dec
577
578 ra = rev(atan2( sin(l) * cos(eps) - tan(b) * sin(eps), cos(l) ))
579 dec = asin( sin(b) * cos(eps) + cos(b) * sin(eps) * sin(l) )
580
581 end subroutine ecl_2_eq
582 !*********************************************************************************************************************************
583
584
585 !*********************************************************************************************************************************
586 !> \brief Convert (geocentric) spherical equatorial coordinates RA, Dec (and eps) to spherical ecliptical coordinates l,b
587 !!
588 !! \param ra Right ascension (rad)
589 !! \param dec Declination (rad)
590 !! \param eps Obliquity of the ecliptic
591 !!
592 !! \param l Longitude (rad) (output)
593 !! \param b Latitude (rad) (output)
594 !!
595 !! \see Expl. Supl. t.t. Astronimical Almanac 3rd Ed, Eq.14.42
596
597 subroutine eq_2_ecl(ra,dec,eps, l,b)
598 use sufr_kinds, only: double
599 use sufr_angles, only: rev
600
601 implicit none
602 real(double), intent(in) :: ra,dec,eps
603 real(double), intent(out) :: l,b
604
605 l = rev(atan2( sin(ra) * cos(eps) + tan(dec) * sin(eps), cos(ra) ))
606 b = asin( sin(dec) * cos(eps) - cos(dec) * sin(eps) * sin(ra) )
607
608 end subroutine eq_2_ecl
609 !*********************************************************************************************************************************
610
611
612 !*********************************************************************************************************************************
613 !> \brief Convert spherical equatorial coordinates (RA, dec, agst) to spherical horizontal coordinates (hh, az, alt)
614 !!
615 !! \param ra Right ascension
616 !! \param dec Declination
617 !! \param agst Greenwich sidereal time
618 !!
619 !! \param hh Local hour angle (output)
620 !! \param az Azimuth (output)
621 !! \param alt Altitude (output)
622 !!
623 !! \param lat Geographical latitude (rad), optional
624 !! \param lon Geographical longitude (rad), optional
625 !!
626 !! \see Expl. Supl. t.t. Astronimical Almanac 3rd Ed, Eq.14.44, 14.45
627
628 subroutine eq2horiz(ra,dec,agst, hh,az,alt, lat,lon)
629 use sufr_kinds, only: double
630 use sufr_angles, only: rev,rev2
631 use thesky_local, only: lat0,lon0
632
633 implicit none
634 real(double), intent(in) :: ra,dec,agst
635 real(double), intent(in), optional :: lat,lon
636 real(double), intent(out) :: hh,az,alt
637 real(double) :: llat,llon, sinhh,coshh, sinlat,coslat, sindec,cosdec,tandec
638
639 ! Lat, lon added later as optional variables:
640 llat = lat0
641 llon = lon0
642 if(present(lat)) llat = lat
643 if(present(lon)) llon = lon
644
645 hh = rev( agst + llon - ra ) ! Local Hour Angle (agst since ra is also corrected for nutation?)
646
647 ! Some preparation, saves ~29%:
648 sinhh = sin(hh)
649 coshh = cos(hh)
650 sinlat = sin(llat)
651 coslat = sqrt(1.d0-sinlat**2) ! Cosine of a latitude is always positive
652 sindec = sin(dec)
653 cosdec = sqrt(1.d0-sindec**2) ! Cosine of a declination is always positive
654 tandec = sindec/cosdec
655
656 az = rev( atan2( sinhh, coshh * sinlat - tandec * coslat )) ! Azimuth
657 alt = rev2( asin( sinlat * sindec + coslat * cosdec * coshh )) ! Altitude
658
659 end subroutine eq2horiz
660 !*********************************************************************************************************************************
661
662
663 !*********************************************************************************************************************************
664 !> \brief Convert spherical horizontal coordinates (az, alt, agst) to spherical equatorial coordinates (hh, RA, dec)
665 !!
666 !! \param az Azimuth
667 !! \param alt Altitude
668 !! \param agst Greenwich sidereal time
669 !!
670 !! \param hh Local hour angle (output)
671 !! \param ra Right ascension (output)
672 !! \param dec Declination (output)
673 !!
674 !! \param lat Geographical latitude (rad), optional
675 !! \param lon Geographical longitude (rad), optional
676 !!
677 !! \see Expl. Supl. t.t. Astronimical Almanac 3rd Ed, Eq.14.46, 14.47
678
679 subroutine horiz2eq(az,alt, agst, hh,ra,dec, lat,lon)
680 use sufr_kinds, only: double
681 use sufr_angles, only: rev,rev2
682 use thesky_local, only: lat0,lon0
683
684 implicit none
685 real(double), intent(in) :: az,alt,agst
686 real(double), intent(out) :: ra,dec,hh
687 real(double), intent(in), optional :: lat,lon
688 real(double) :: sinlat,coslat, cosaz,sinalt,cosalt,tanalt, llat,llon
689
690 ! Lat, lon added later as optional variables:
691 llat = lat0
692 llon = lon0
693 if(present(lat)) llat = lat
694 if(present(lon)) llon = lon
695
696 ! Some preparation to save time:
697 sinlat = sin(llat)
698 coslat = sqrt(1.d0-sinlat**2) ! Cosine of a latitude is always positive
699 cosaz = cos(az)
700 sinalt = sin(alt)
701 cosalt = sqrt(1.d0-sinalt**2) ! Cosine of a altitude is always positive
702 tanalt = sinalt/cosalt
703
704 hh = rev( atan2( sin(az), cosaz * sinlat + tanalt * coslat )) ! Local Hour Angle
705 dec = rev2( asin( sinlat * sinalt - coslat * cosalt * cosaz )) ! Declination
706 ra = rev( agst + llon - hh ) ! Right ascension
707
708 end subroutine horiz2eq
709 !*********************************************************************************************************************************
710
711
712 !*********************************************************************************************************************************
713 !> \brief Convert spherical equatorial coordinates (RA, dec) to spherical galactic coordinates (l,b), for J2000.0!!!
714 !!
715 !! \param ra Right ascension
716 !! \param dec Declination
717 !!
718 !! \param l Longitude (output)
719 !! \param b Latitude (output)
720 !!
721 !! \see Meeus, Astronomical Algorithms, 1998, Ch.13
722
723 subroutine eq2gal(ra,dec, l,b)
724 use sufr_kinds, only: double
725 use sufr_angles, only: rev
726
727 implicit none
728 real(double), intent(in) :: ra,dec
729 real(double), intent(out) :: l,b
730 real(double) :: ra0,dec0,l0
731
732 ra0 = 3.36603472d0 ! RA of GNP, in rad, J2000.0 (12h51m26.3s ~ 192.8596deg)
733 dec0 = 0.473478737d0 ! Decl. of gal. NP in rad, J2000.0 (27d07'42"=27.1283deg)
734 l0 = 5.287194d0 ! J2000.0?, = 302.9339deg (the 303deg in Meeus, p.94) = l0+180deg
735
736 l = rev( l0 - atan2( sin(ra0-ra), cos(ra0-ra) * sin(dec0) - tan(dec) * cos(dec0) ))
737 b = asin( sin(dec) * sin(dec0) + cos(dec) * cos(dec0) * cos(ra0-ra) )
738
739 end subroutine eq2gal
740 !*********************************************************************************************************************************
741
742
743 !*********************************************************************************************************************************
744 !> \brief Convert spherical galactic coordinates (l,b) to spherical equatorial coordinates (RA, dec), for J2000.0!!!
745 !!
746 !! \param l Longitude
747 !! \param b Latitude
748 !!
749 !! \param ra Right ascension (output)
750 !! \param dec Declination (output)
751
752 subroutine gal2eq(l,b, ra,dec)
753 use sufr_kinds, only: double
754 use sufr_angles, only: rev
755
756 real(double), intent(in) :: l,b
757 real(double), intent(out) :: ra,dec
758 real(double) :: ra0,dec0,l0
759
760 ra0 = 0.22444207d0 ! RA of GNP - 12h, in rad, J2000.0 (12h51m26.3s - 12h ~ 12.8596deg)
761 dec0 = 0.473478737d0 ! Decl. of gal. NP in rad, J2000.0 (from Sterrengids?: 27d07'42"=27.1283deg)
762 l0 = 2.145601d0 ! J2000.0?, = 123.9339deg (the 123deg in Meeus, p.94)
763
764 ra = rev(atan2( sin(l-l0), cos(l-l0) * sin(dec0) - tan(b) * cos(dec0)) + ra0)
765 dec = asin( sin(b) * sin(dec0) + cos(b) * cos(dec0) * cos(l-l0))
766
767 end subroutine gal2eq
768 !*********************************************************************************************************************************
769
770
771
772 !*********************************************************************************************************************************
773 !> \brief Convert spherical ecliptical coordinates from the geocentric to the topocentric system
774 !!
775 !! \param gcl Geocentric ecliptic longitude of the object (rad)
776 !! \param gcb Geocentric ecliptic latitude of the object (rad)
777 !! \param gcr Geocentric distance of the object
778 !! \param gcs Geocentric semi-diameter of the object (rad)
779 !!
780 !! \param eps Obliquity of the ecliptic (rad)
781 !! \param lst Local sidereal time (rad)
782 !!
783 !! \param tcl Topocentric ecliptic longitude of the object (rad) (output)
784 !! \param tcb Topocentric ecliptic latitude of the object (rad) (output)
785 !! \param tcs Topocentric semi-diameter of the object (rad) (output)
786 !!
787 !!
788 !! \param lat Latitude of the observer (rad, optional)
789 !! \param hgt Altitude/elevation of the observer above sea level (metres, optional)
790 !!
791 !! \note lat/lat0 and hgt/height can be provided through the module TheSky_local (rad, and m), or through the optional arguments.
792 !! Note that using the latter will update the former!
793 !!
794 !! \see Meeus, Astronomical Algorithms, 1998, Ch. 11 and 40
795
796 subroutine geoc2topoc_ecl(gcl,gcb, gcr,gcs, eps,lst, tcl,tcb, tcs, lat,hgt)
797 use sufr_kinds, only: double
798 use sufr_constants, only: earthr,au
799 use sufr_angles, only: rev
800 use thesky_local, only: lat0, height
801
802 implicit none
803 real(double), intent(in) :: gcl,gcb,gcr,gcs,eps,lst
804 real(double), intent(in), optional :: lat,hgt
805 real(double), intent(out) :: tcl,tcb,tcs
806 real(double) :: llat,lhgt, ba,re,u,rs,rc,sinHp,n
807
808 ! Handle optional variables:
809 llat = lat0
810 lhgt = height
811 if(present(lat)) llat = lat
812 if(present(hgt)) lhgt = hgt
813
814 ! Meeus, Ch.11, p.82:
815 ba = 0.996647189335d0 ! b/a = 1-f: flattening of the Earth - WGS84 ellipsoid
816 re = earthr*1.d-2 ! Equatorial radius of the Earth in metres
817 ! (http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf)
818
819 u = atan(ba*tan(llat))
820 rs = ba*sin(u) + lhgt/re*sin(llat)
821 rc = cos(u) + lhgt/re*cos(llat)
822
823
824 sinhp = sin(earthr/au)/gcr ! Sine of the horizontal parallax, Meeus, Eq. 40.1
825
826 ! Meeus, Ch.40, p.282:
827 n = cos(gcl)*cos(gcb) - rc*sinhp*cos(lst)
828
829 tcl = rev( atan2( sin(gcl)*cos(gcb) - sinhp*(rs*sin(eps) + rc*cos(eps)*sin(lst)) , n ) ) ! Topocentric longitude
830 tcb = atan((cos(tcl)*(sin(gcb) - sinhp*(rs*cos(eps) - rc*sin(eps)*sin(lst))))/n) ! Topocentric latitude
831 tcs = asin(cos(tcl)*cos(tcb)*sin(gcs)/n) ! Topocentric semi-diameter
832
833 end subroutine geoc2topoc_ecl
834 !*********************************************************************************************************************************
835
836
837 !*********************************************************************************************************************************
838 !> \brief Convert geocentric equatorial coordinates to topocentric
839 !!
840 !! \param gcra Geocentric right ascension
841 !! \param gcd Geocentric declination
842 !! \param gcr Geocentric distance
843 !! \param gch Geocentric hour angle
844 !!
845 !! \param tcra Topocentric right ascension (output)
846 !! \param tcd Topocentric declination (output)
847 !!
848 !! \param lat Latitude of the observer (rad, optional)
849 !! \param hgt Altitude/elevation of the observer above sea level (metres, optional)
850 !!
851 !! \note lat/lat0 and hgt/height can be provided through the module TheSky_local (rad, and m), or through the optional arguments.
852 !! Note that using the latter will update the former!
853 !!
854 !!
855 !! \see Meeus, Astronomical Algorithms, 1998, Ch. 11 and 40
856
857 subroutine geoc2topoc_eq(gcra,gcd,gcr,gch, tcra,tcd, lat,hgt)
858 use sufr_kinds, only: double
859 use sufr_constants, only: earthr,au
860 use thesky_local, only: lat0, height
861
862 implicit none
863 real(double), intent(in) :: gcra,gcd,gcr,gch
864 real(double), intent(in), optional :: lat,hgt
865 real(double), intent(out) :: tcra,tcd
866 real(double) :: llat,lhgt, ba,re,u,rs,rc,sinHp,dra
867
868 ! Handle optional variables:
869 llat = lat0
870 lhgt = height
871 if(present(lat)) llat = lat
872 if(present(hgt)) lhgt = hgt
873
874 ! Meeus, Ch.11, p.82:
875 ba = 0.996647189335d0 ! b/a = 1-f: flattening of the Earth - WGS84 ellipsoid
876 re = earthr*1.d-2 ! Equatorial radius of the Earth in metres
877 ! (http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf)
878
879 u = atan(ba*tan(llat))
880 rs = ba*sin(u) + lhgt/re*sin(llat)
881 rc = cos(u) + lhgt/re*cos(llat)
882
883 ! Meeus Ch.40:
884 sinhp = sin(earthr/au)/gcr ! Sine of the horizontal parallax, Meeus, Eq. 40.1
885
886 dra = atan2( -rc*sinhp*sin(gch) , cos(gcd)-rc*sinhp*cos(gch) ) ! Meeus, Eq. 40.2
887 tcra = gcra + dra ! Topocentric right ascension
888 tcd = atan2( (sin(gcd)-rs*sinhp)*cos(dra) , cos(gcd)-rc*sinhp*cos(gch) ) ! Topocentric declination - Meeus, Eq. 40.3
889
890 end subroutine geoc2topoc_eq
891 !*********************************************************************************************************************************
892
893
894 !*********************************************************************************************************************************
895 !> \brief Compute the atmospheric refraction for a given true altitude. You should add the result to the uncorrected altitude
896 !! in order to obtain the observed altitude. Return 0 if alt + refract < -0.3°.
897 !!
898 !! \param alt True (computed) altitude (rad)
899 !! \param press Air pressure (hPa; optional)
900 !! \param temp Air temperature (degrees Celsius; optional)
901 !!
902 !! \retval refract Refraction in altitude (rad). You should *add* the result to the uncorrected altitude.
903 !!
904 !! \see
905 !! - T. Saemundsson, "Atmospheric refraction", Sky & Telescope vol.72, p.70 (1986), converted to radians;
906 !! - Meeus (1998), Eq. 16.4 ff.
907 !!
908
909 function refract(alt, press,temp)
910 use sufr_kinds, only: double
911 use sufr_constants, only: pio2, d2r
912
913 implicit none
914 real(double), intent(in) :: alt
915 real(double), intent(in), optional :: press,temp
916 real(double) :: refract
917
918 if(abs(alt).gt.pio2) then ! for |alt| > 90° refraction is meaningless
919 refract = 0.d0
920 else
921 if(alt.lt.-1.102d0*d2r) then ! Maximum possible refraction on Earth, for T=-90°C and P=1085 mbar
922 ! Without this, refract reaches a maximum around alt=-1.9° and becomes similar to the value for |alt| for lower alt (is this bad?).
923 ! See https://en.wikipedia.org/wiki/Lowest_temperature_recorded_on_Earth, https://en.wikipedia.org/wiki/Atmospheric_pressure
924 refract = 0.663455d0*d2r ! Refraction for alt=-1.102°, but P=1010 mbar and T=10°C.
925 else
926 refract = 2.9670597d-4/tan(alt + 3.137559d-3/(alt + 8.91863d-2)) ! Overstated accuracy in translation from degrees
927 end if
928
929 if(present(press)) refract = refract * press/1010.d0 ! Correct for pressure
930 if(present(temp)) refract = refract * 283.d0/(273.d0 + temp) ! Correct for temperature
931 end if
932
933 end function refract
934 !*********************************************************************************************************************************
935
936
937
938 !*********************************************************************************************************************************
939 !> \brief Compute the atmospheric refraction for a given true altitude, using single-precision values. This
940 !! is a wrapper for refract().
941 !!
942 !! \param alt True (computed) altitude (rad)
943 !! \param press Air pressure (hPa; optional)
944 !! \param temp Air temperature (degrees Celsius; optional)
945 !!
946 !! \retval refract Refraction in altitude (rad). You should add the result to the uncorrected altitude.
947 !!
948 !! \see refract() for details.
949
950 function refract_sp(alt, press,temp)
951 use sufr_kinds, only: double
952 use sufr_constants, only: rpio2
953
954 implicit none
955 real, intent(in) :: alt
956 real, intent(in), optional :: press,temp
957 real :: refract_sp
958 real(double) :: lalt, lpress,ltemp
959
960 if(abs(alt).ge.rpio2) then ! |alt| >= 90° refraction is meaningless
961 refract_sp = 0.0
962 else
963 lalt = dble(alt)
964 lpress = 1010.d0 ! Default value for refract()
965 if(present(press)) lpress = dble(press)
966 ltemp = 10.d0 ! Default value for refract()
967 if(present(temp)) ltemp = dble(temp)
968
969 refract_sp = real(refract(lalt, lpress,ltemp))
970 end if
971
972 end function refract_sp
973 !*********************************************************************************************************************************
974
975
976
977 !*********************************************************************************************************************************
978 !> \brief Compute the atmospheric refraction of light for a given *true* altitude. Return 0 for alt<-0.9°.
979 !! This is a wrapper for aref(), which does the opposite (compute refraction for an *apparent* zenith angle).
980 !! This is an expensive way to go about(!)
981 !!
982 !! \param alt0 The true (theoretical, computed) altitude of the object in radians
983 !!
984 !! \param h0 The height of the observer above sea level in metres
985 !! \param lat0 The latitude of the observer in radians
986 !!
987 !! \param t0 The temperature at the observer in degrees Celsius (e.g. 10)
988 !! \param p0 The pressure at the observer in millibars (e.g. 1010)
989 !! \param rh The relative humidity at the observer in percent (e.g. 50)
990 !!
991 !! \param lam The wavelength of the light at the observer in nanometres (e.g. 550)
992 !! \param dTdh The temperature lapse rate dT/dh in Kelvin/metre in the troposphere (only the absolute value is used; e.g. 0.0065)
993 !!
994 !! \param eps The desired precision in *arcseconds* (e.g. 1.d-3)
995 !!
996 !! \retval atmospheric_refraction The refraction at the observer in radians
997 !!
998 !!
999 !! \todo Adapt aref() to compute the integral in the other direction for a direct method(?)
1000 !!
1001
1002 function atmospheric_refraction(alt0, h0,lat0, t0,p0,rh, lam,dTdh, eps)
1003 use sufr_kinds, only: double
1004 use sufr_constants, only: d2r,r2d, pio2
1005
1006 implicit none
1007 real(double), intent(in) :: alt0, h0,lat0, t0,p0, rh,lam, dtdh, eps
1008 real(double) :: atmospheric_refraction, z0,zi,zio,dz, ph,lamm,t0k,rhfr
1009
1011 if(abs(alt0).ge.pio2) return ! No refraction if |alt| >= 90°
1012 if(alt0.lt.-0.9d0*d2r) return ! No refraction if alt < -0.9° - allowing this can cause aref() to fail
1013
1014 ! Convert some variables/units:
1015 z0 = 90.d0 - alt0*r2d ! Altitude (rad) -> zenith angle (in *degrees!*)
1016 ph = lat0*r2d ! Latitude: rad -> *degrees!*
1017 lamm = lam / 1000.d0 ! Wavelength: nanometre -> micrometre
1018 t0k = t0 + 273.15d0 ! Temperature: degC -> K
1019 rhfr = rh/100.d0 ! Relative humidity: % -> fraction
1020
1021 ! Set initial values for the iteration:
1022 zi = z0
1023 zio = huge(zio)
1024 dz = huge(dz)
1025
1026 do while(dz.gt.eps/3600.d0) ! Iterate until the desired accuracy is achieved
1027 atmospheric_refraction = aref(zi, h0,ph, t0k,p0,rhfr, lamm,dtdh, eps) ! Do the *inverse* calculation: refraction for given *apparent* alt (in *degrees!*)
1028 zi = z0 - atmospheric_refraction ! Compute the new zenith angle
1029 dz = abs(zi-zio) ! Compute the absolute difference with the previous iteration
1030 zio = zi ! Store the current value as the 'old value' for the next iteration
1031 end do
1032
1033 atmospheric_refraction = atmospheric_refraction * d2r ! Convert result back to radians
1034
1035 end function atmospheric_refraction
1036 !*********************************************************************************************************************************
1037
1038
1039 !*********************************************************************************************************************************
1040 !> \brief Compute the atmospheric refraction of light for a given *apparent* (observed) altitude of a celestial object
1041 !! and a given observer. Note that you will usually want to use the *true* altitude as input instead.
1042 !! The method is based on N.A.O Technical Notes 59 and 63 and a paper by Standish and Auer 'Astronomical Refraction:
1043 !! Computational Method for all Zenith Angles'. Return 0 if alt<-1.102°.
1044 !!
1045 !! \param alt0 The observed (apparent) altitude of the object (rad).
1046 !!
1047 !! \param h0 The height of the observer above sea level (metres).
1048 !! \param ph The latitude of the observer (rad).
1049 !!
1050 !! \param t0 The temperature at the observer's site; optional, default 10 (°C).
1051 !! \param p0 The air pressure at the observer's site; optional, default: 1010 (mbar).
1052 !! \param rh The relative humidity at the observer's site; optional, default: 0.5 (fraction: 0-1).
1053 !!
1054 !! \param lam The wavelength of the light; optional, default: 550 (nm).
1055 !! \param dTdh The temperature gradient dT/dh in the troposphere; the absolute value is used; optional, default: 6.5e-3 (Kelvin/metre).
1056 !!
1057 !! \param eps The desired precision; optional, default: 1e-4; don't make this smaller than ~1e-8 (rad).
1058 !!
1059 !! \retval atmospheric_refraction_apparent The refraction at the observer (rad).
1060 !!
1061 !! \see Hohenkerk & Sinclair, HMNAO technical note 63 (1985).
1062
1063 function atmospheric_refraction_apparent(alt0, h0,ph, t0,p0,rh, lam,dTdh,eps)
1064 use sufr_kinds, only: double
1065 use sufr_constants, only: pio2, d2r,r2d, earthr
1066
1067 implicit none
1068 real(double), intent(in) :: alt0, h0,ph
1069 real(double), intent(in), optional :: t0,p0, rh,lam, dtdh, eps
1070
1071 real(double), parameter :: gcr=8314.36d0, md=28.966d0, mw=18.016d0, gamma=18.36d0
1072 real(double), parameter :: ht=11.d3, hs=80.d3, z2=11.2684d-06
1073
1074 integer :: i,in,is,istart, j,k
1075 real(double) :: atmospheric_refraction_apparent, z0, t0l,p0l,rhl, laml, dtdhl, epsl
1076 real(double) :: n,n0,nt,nts,ns,a(10), dndr,dndr0,dndrs,dndrt,dndrts, f,f0, fb,fe,ff,fo,fs,ft,fts,gb, h
1077 real(double) :: pw0,r,r0,refo,refp,reft,rg,rs,rt, sk0,step,t,t0o,tg,tt,z,z1,zs,zt,zts, earthrm
1078
1080 if(alt0.lt.-1.102d0*d2r) return ! Cannot *observe* object below horizon. However, for the *inverse* case, using atmospheric_refraction() to iterate on this function, we may start slightly below h=0.
1081 ! 1.102° is the maximum possible refraction on Earth, for T=-90°C and P=1085 mbar, according to refract().
1082
1083 ! Original algorithm was in degrees and kelvin and used the zenith angle as input iso altitude:
1084 z0 = (pio2 - alt0)*r2d ! alt -> za; rad -> deg
1085
1086 ! Deal with optional parameters:
1087 if(present(t0)) then ! Temperature
1088 t0l = t0 + 273.15d0 ! Use user-specified value; °C -> K
1089 else
1090 t0l = 283.15d0 ! 10°C -> K by default
1091 end if
1092 if(present(p0)) then ! Air pressure
1093 p0l = p0 ! Use user-specified value
1094 else
1095 p0l = 1010 ! 1010 mbar by default
1096 end if
1097 if(present(rh)) then ! Relative humidity
1098 rhl = rh ! Use user-specified value
1099 else
1100 rhl = 0.5 ! 50% by default
1101 end if
1102 if(present(lam)) then ! Wavelength
1103 laml = lam*1.d-3 ! Use user-specified value; nm -> μm
1104 else
1105 laml = 550*1.d-3 ! 550 nm -> μm by default
1106 end if
1107 if(present(dtdh)) then ! Temperature gradient
1108 dtdhl = dtdh ! Use user-specified value
1109 else
1110 dtdhl = 6.5d-3 ! Use default value: 6.5e-3 K/m
1111 end if
1112 if(present(eps)) then ! Desired accuracy
1113 epsl = eps ! Use user-specified value
1114 else
1115 epsl = 1d-4 ! Use default value: 1e-4 rad ~ 0.006°
1116 end if
1117
1118
1119 ! Always defined:
1120 z = 0.d0; reft=0.d0
1121
1122 ! Set up parameters defined at the observer for the atmosphere:
1123 gb = 9.784d0*(1.d0 - 0.0026d0*cos(2.d0*ph) - 0.00000028d0*h0)
1124 z1 = (287.604d0 + 1.6288d0/(laml**2) + 0.0136d0/(laml**4)) * (273.15d0/1013.25d0)*1.d-6
1125
1126 a(1) = abs(dtdhl)
1127 a(2) = (gb*md)/gcr
1128 a(3) = a(2)/a(1)
1129 a(4) = gamma
1130 pw0 = rhl*(t0l/247.1d0)**a(4)
1131 a(5) = pw0*(1.d0 - mw/md)*a(3)/(a(4)-a(3))
1132 a(6) = p0l + a(5)
1133 a(7) = z1*a(6)/t0l
1134 a(8) = ( z1*a(5) + z2*pw0)/t0l
1135 a(9) = (a(3) - 1.d0)*a(1)*a(7)/t0l
1136 a(10) = (a(4) - 1.d0)*a(1)*a(8)/t0l
1137
1138 ! At the Observer:
1139 earthrm = earthr*1.d-2 ! cm -> m
1140 r0 = earthrm + h0
1141 call troposphere_model(r0,t0l,a,r0,t0o,n0,dndr0)
1142 sk0 = n0 * r0 * sin(z0*d2r)
1143
1144
1145 f0 = refi(r0,n0,dndr0)
1146
1147 ! At the Tropopause in the Troposphere:
1148 rt = earthrm + ht
1149 call troposphere_model(r0,t0l,a,rt,tt,nt,dndrt)
1150 zt = asin(sk0/(rt*nt))*r2d
1151 ft = refi(rt,nt,dndrt)
1152
1153 ! At the Tropopause in the Stratosphere:
1154 call stratosphere_model(rt,tt,nt,a(2),rt,nts,dndrts)
1155 zts = asin(sk0/(rt*nts))*r2d
1156 fts = refi(rt,nts,dndrts)
1157
1158 ! At the stratosphere limit:
1159 rs = earthrm + hs
1160 call stratosphere_model(rt,tt,nt,a(2),rs,ns,dndrs)
1161 zs = asin(sk0/(rs*ns))*r2d
1162 fs = refi(rs,ns,dndrs)
1163
1164 ! Integrate the refraction integral in the troposhere and stratosphere,
1165 ! i.e. total refraction = refraction troposhere + refraction stratopshere
1166
1167 ! Initial step lengths etc:
1168 refo = -huge(refo)
1169 is = 16
1170 do k = 1,2
1171 istart = 0
1172 fe = 0.d0
1173 fo = 0.d0
1174
1175 if(k.eq.1) then
1176 h = (zt - z0)/dble(is)
1177 fb = f0
1178 ff = ft
1179 else if(k.eq.2) then
1180 h = (zs - zts )/dble(is)
1181 fb = fts
1182 ff = fs
1183 end if
1184
1185 in = is - 1
1186 is = is/2
1187 step = h
1188
1189200 continue
1190
1191 do i = 1,in
1192 if(i.eq.1.and.k.eq.1) then
1193 z = z0 + h
1194 r = r0
1195 else if(i.eq.1.and.k.eq.2) then
1196 z = zts + h
1197 r = rt
1198 else
1199 z = z + step
1200 end if
1201
1202 ! Given the Zenith distance (Z) find R:
1203 rg = r
1204 do j = 1,4
1205 if(k.eq.1) then
1206 call troposphere_model(r0,t0l,a,rg,tg,n,dndr)
1207 else if(k.eq.2) then
1208 call stratosphere_model(rt,tt,nt,a(2),rg,n,dndr)
1209 end if
1210 rg = rg - ( (rg*n - sk0/sin(z*d2r))/(n + rg*dndr) )
1211 end do
1212 r = rg
1213
1214 ! Find Refractive index and Integrand at R:
1215 if(k.eq.1) then
1216 call troposphere_model(r0,t0l,a,r,t,n,dndr)
1217 else if(k.eq.2) then
1218 call stratosphere_model(rt,tt,nt,a(2),r,n,dndr)
1219 end if
1220
1221 f = refi(r,n,dndr)
1222
1223 if(istart.eq.0.and.mod(i,2).eq.0) then
1224 fe = fe + f
1225 else
1226 fo = fo + f
1227 end if
1228 end do
1229
1230 ! Evaluate the integrand using Simpson's Rule:
1231 refp = h*(fb + 4.d0*fo + 2.d0*fe + ff)/3.d0
1232
1233 if(abs(refp-refo).gt.0.5d0*epsl*r2d) then
1234 is = 2*is
1235 in = is
1236 step = h
1237 h = h/2.d0
1238 fe = fe + fo
1239 fo = 0.d0
1240 refo = refp
1241 if(istart.eq.0) istart = 1
1242 goto 200
1243 end if
1244 if(k.eq.1) reft = refp
1245 end do
1246
1247 atmospheric_refraction_apparent = (reft + refp)*d2r ! Convert result from degrees to radians
1248
1250 !*********************************************************************************************************************************
1251
1252
1253 !*********************************************************************************************************************************
1254 !> \brief Compute the atmospheric refraction of light for a given *apparent* (observed) zenith angle.
1255 !! The method is based on N.A.O Technical Notes 59 and 63 and a paper by Standish and Auer 'Astronomical Refraction:
1256 !! Computational Method for all Zenith Angles'. Return 0 if z0>91.102°.
1257 !!
1258 !! \param z0 The observed zenith distance of the object in *degrees*(!)
1259 !!
1260 !! \param h0 The height of the observer above sea level in metres
1261 !! \param ph The latitude of the observer in *degrees*(!)
1262 !!
1263 !! \param t0 The temperature at the observer in Kelvin
1264 !! \param p0 The pressure at the observer in millibars
1265 !! \param rh The relative humidity at the observer as a fraction (0-1)
1266 !!
1267 !! \param lam The wavelength of the light at the observer in micrometres
1268 !! \param dTdh The temperature lapse rate dT/dh in Kelvin/metre in the troposphere; the absolute value is used
1269 !!
1270 !! \param eps The desired precision in *arcseconds*(!)
1271 !!
1272 !! \retval aref The refraction at the observer in *degrees*(!)
1273 !!
1274 !! \see Hohenkerk & Sinclair, HMNAO technical note 63 (1985)
1275
1276 function aref(z0, h0,ph, t0,p0,rh, lam,dTdh,eps)
1277 use sufr_kinds, only: double
1278 use sufr_constants, only: d2r,r2d
1279
1280 implicit none
1281 real(double), intent(in) :: z0, h0,ph, t0,p0, rh,lam, dtdh, eps
1282 real(double) :: aref, alt0
1283
1284 alt0 = (90.d0-z0)*d2r ! za -> alt; deg -> rad
1285 ! input: deg->rad, K->°C, μm->nm; output: rad -> deg:
1286 aref = atmospheric_refraction_apparent(alt0, h0,ph*d2r, t0-273.15d0,p0,rh, lam*1.d3, dtdh, eps*d2r)*r2d
1287
1288 end function aref
1289 !*********************************************************************************************************************************
1290
1291
1292 !*********************************************************************************************************************************
1293 !> \brief The refraction integrand
1294 !!
1295 !! \param r The current distance from the centre of the Earth in metres
1296 !! \param n The refractive index at R
1297 !! \param dndr The rate the refractive index is changing at R
1298 !!
1299 !! \retval refi The integrand of the refraction function
1300
1301 function refi(r,n,dndr)
1302 use sufr_kinds, only: double
1303 implicit none
1304 real(double), intent(in) :: r,n,dndr
1305 real(double) :: refi
1306
1307 refi = r*dndr/(n + r*dndr)
1308 end function refi
1309 !*********************************************************************************************************************************
1310
1311
1312 !*********************************************************************************************************************************
1313 !> \brief Atmospheric model for the troposphere
1314 !!
1315 !! \param r0 The height of the observer from the centre of the Earth
1316 !! \param t0 The temperature at the observer in Kelvin
1317 !! \param a Constants defined at the observer
1318 !! \param r The current distance from the centre of the Earth in metres
1319 !!
1320 !! \param t The temperature at R in Kelvin (output)
1321 !! \param n The refractive index at R (output)
1322 !! \param dndr The rate the refractive index is changing at R (output)
1323
1324 subroutine troposphere_model(r0,t0, a,r, t,n,dndr)
1325 use sufr_kinds, only: double
1326 implicit none
1327 real(double), intent(in) :: r0,t0, a(10),r
1328 real(double), intent(out) :: t,n,dndr
1329 real(double) :: tt0, tt01,tt02
1330
1331 t = t0 - a(1)*(r-r0)
1332 tt0 = t/t0
1333 tt01 = tt0**(a(3)-2.d0)
1334 tt02 = tt0**(a(4)-2.d0)
1335
1336 n = 1.d0 + ( a(7)*tt01 - a(8)*tt02 )*tt0
1337 dndr = -a(9)*tt01 + a(10)*tt02
1338
1339 end subroutine troposphere_model
1340 !*********************************************************************************************************************************
1341
1342
1343 !*********************************************************************************************************************************
1344 !> \brief Atmospheric model for the stratosphere
1345 !!
1346 !! \param rt The height of the tropopause from the centre of the Earth in metres
1347 !! \param tt The temperature at the tropopause in Kelvin
1348 !! \param nt The refractive index at the tropopause
1349 !! \param a Constant of the atmospheric model = G*MD/R
1350 !! \param r The current distance from the centre of the Earth in metres
1351 !!
1352 !! \param n The refractive index at R (output)
1353 !! \param dndr The rate the refractive index is changing at R (output)
1354
1355 subroutine stratosphere_model(rt,tt,nt, a,r, n,dndr)
1356 use sufr_kinds, only: double
1357 implicit none
1358 real(double), intent(in) :: rt,tt,nt, a,r
1359 real(double), intent(out) :: n,dndr
1360 real(double) :: b
1361
1362 b = a/tt
1363 n = 1.d0 + (nt - 1.d0)*exp(-b*(r-rt))
1364 dndr = -b*(nt-1.d0)*exp(-b*(r-rt))
1365
1366 end subroutine stratosphere_model
1367 !*********************************************************************************************************************************
1368
1369end module thesky_coordinates
1370!***********************************************************************************************************************************
1371
Procedures for coordinates.
subroutine precess_xyz(jd1, jd2, x, y, z)
Compute the precession of the equinoxes in rectangular coordinates, from jd1 to jd2.
real(double) function refract(alt, press, temp)
Compute the atmospheric refraction for a given true altitude. You should add the result to the uncorr...
subroutine eq2horiz(ra, dec, agst, hh, az, alt, lat, lon)
Convert spherical equatorial coordinates (RA, dec, agst) to spherical horizontal coordinates (hh,...
real(double) function refi(r, n, dndr)
The refraction integrand.
subroutine horiz2eq(az, alt, agst, hh, ra, dec, lat, lon)
Convert spherical horizontal coordinates (az, alt, agst) to spherical equatorial coordinates (hh,...
subroutine eq_2_ecl(ra, dec, eps, l, b)
Convert (geocentric) spherical equatorial coordinates RA, Dec (and eps) to spherical ecliptical coord...
real(double) function atmospheric_refraction(alt0, h0, lat0, t0, p0, rh, lam, dtdh, eps)
Compute the atmospheric refraction of light for a given true altitude. Return 0 for alt<-0....
subroutine stratosphere_model(rt, tt, nt, a, r, n, dndr)
Atmospheric model for the stratosphere.
subroutine gal2eq(l, b, ra, dec)
Convert spherical galactic coordinates (l,b) to spherical equatorial coordinates (RA,...
real(double) function atmospheric_refraction_apparent(alt0, h0, ph, t0, p0, rh, lam, dtdh, eps)
Compute the atmospheric refraction of light for a given apparent (observed) altitude of a celestial o...
subroutine aberration_eq(jd, ra, dec, dra, ddec, eps0)
Correct equatorial coordinates for annual aberration - moderate accuracy, use for stars.
subroutine precess_orb(jd1, jd2, i, o1, o2)
Compute the precession of the equinoxes in orbital elements.
real function refract_sp(alt, press, temp)
Compute the atmospheric refraction for a given true altitude, using single-precision values....
subroutine ecl_spher_2_eq_rect(l, b, r, eps, x, y, z)
Convert spherical, ecliptical coordinates to rectangular, equatorial coordinates of an object - both ...
subroutine eq2gal(ra, dec, l, b)
Convert spherical equatorial coordinates (RA, dec) to spherical galactic coordinates (l,...
subroutine rect_2_spher(x, y, z, l, b, r)
Convert rectangular coordinates x,y,z to spherical coordinates l,b,r.
subroutine precess_eq(jd1, jd2, a1, d1)
Compute the precession of the equinoxes in equatorial coordinates, from jd1 to jd2.
subroutine ecl_2_eq(l, b, eps, ra, dec)
Convert (geocentric) spherical ecliptical coordinates l,b (and eps) to spherical equatorial coordinat...
real(double) function aref(z0, h0, ph, t0, p0, rh, lam, dtdh, eps)
Compute the atmospheric refraction of light for a given apparent (observed) zenith angle....
subroutine fk5(t, l, b)
Convert coordinates to the FK5 system.
subroutine hc_spher_2_gc_rect(l, b, r, l0, b0, r0, x, y, z)
Compute the geocentric rectangular coordinates of a planet, from its and the Earth's heliocentric sph...
subroutine calcsunxyz(t1, l0, b0, r0, x, y, z)
Compute the geocentric equatorial rectangular coordinates of the Sun, from Earth's heliocentric,...
subroutine aberration_ecl(t, l0, l, b)
Correct ecliptic longitude and latitiude for annual aberration.
subroutine troposphere_model(r0, t0, a, r, t, n, dndr)
Atmospheric model for the troposphere.
subroutine precess_eq_yr(yr1, yr2, a1, d1)
Compute the precession of the equinoxes in equatorial coordinates, from yr1 to yr2.
subroutine precess_ecl(jd1, jd2, l, b)
Compute the precession of the equinoxes in geocentric ecliptical coordinates.
subroutine geoc2topoc_eq(gcra, gcd, gcr, gch, tcra, tcd, lat, hgt)
Convert geocentric equatorial coordinates to topocentric.
subroutine geoc2topoc_ecl(gcl, gcb, gcr, gcs, eps, lst, tcl, tcb, tcs, lat, hgt)
Convert spherical ecliptical coordinates from the geocentric to the topocentric system.
Local parameters for libTheSky: location, date, time.
Definition modules.f90:44
real(double) height
Definition modules.f90:51
real(double) lon0
Definition modules.f90:51
real(double) lat0
Definition modules.f90:51
Procedures for nutation.
Definition nutation.f90:24
subroutine nutation(t, dpsi, eps0, deps)
Calculate nutation - cheap routine from Meeus - as well as the mean obliquity of the ecliptic.
Definition nutation.f90:45