nova.vsop87 source code

1 /*
2  *  This library is free software; you can redistribute it and/or
3  *  modify it under the terms of the GNU Lesser General Public
4  *  License as published by the Free Software Foundation; either
5  *  version 2 of the License, or (at your option) any later version.
6  *
7  *  This library is distributed in the hope that it will be useful,
8  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
9  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
10  *  Lesser General Public License for more details.
11  *
12  *  You should have received a copy of the GNU General Public License
13  *  along with this program; if not, write to the Free Software
14  *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
15  *
16  *  Some functions in this file use the VSOP87 solution by
17  *  Messrs. Bretagnon and Francou.
18  *
19  *  Copyright (C) 2000 - 2005 Liam Girdwood <lgirdwood@gmail.com>
20  */
21 
22 module nova.vsop87;
23 
24 import std.math;
25 import nova.vsop87;
26 import nova.utility;
27 
28 public import nova.ln_types;
29 
30 struct ln_vsop
31 {
32     double A;
33     double B;
34     double C;
35 }
36 
37 extern (C) {
38 
39 @nogc double ln_calc_series(const ln_vsop[] data, double t) nothrow
40 {
41 	double value = 0.0;
42 
43 	foreach (d; data)
44 		value += d.A * cos(d.B + d.C * t);
45 
46 	return value;
47 }
48 
49 
50 /*! \fn void ln_vsop87_to_fk5(struct ln_helio_posn *position, double JD)
51 * \param position Position to transform.
52 * \param JD Julian day
53 *
54 * Transform from VSOP87 to FK5 reference frame.
55 */
56 /* Equation 31.3 Pg 207.
57 */
58 @nogc void ln_vsop87_to_fk5(ref ln_helio_posn position, double JD) nothrow
59 {
60 	double LL, cos_LL, sin_LL, T, delta_L, delta_B, B;
61 
62 	/* get julian centuries from 2000 */
63 	T =(JD - 2451545.0) / 36525.0;
64 
65 	LL = position.L + (- 1.397 - 0.00031 * T) * T;
66 	LL = ln_deg_to_rad(LL);
67 	cos_LL = cos(LL);
68 	sin_LL = sin(LL);
69 	B = ln_deg_to_rad(position.B);
70 
71 	delta_L = (-0.09033 / 3600.0) + (0.03916 / 3600.0) *
72 			(cos_LL + sin_LL) * tan(B);
73 	delta_B = (0.03916 / 3600.0) * (cos_LL - sin_LL);
74 
75 	position.L += delta_L;
76 	position.B += delta_B;
77 }
78 
79 }