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 * Copyright (C) 2000 - 2005 Liam Girdwood 17 */ 18 19 module nova.precession; 20 21 import std.math; 22 23 import nova.utility; 24 import nova.ln_types; 25 26 // fixed missing *l functions 27 28 extern (C) { 29 30 /* 31 ** Precession 32 */ 33 34 /*! \fn void ln_get_equ_prec(struct ln_equ_posn *mean_position, double JD, struct ln_equ_posn *position) 35 * \param mean_position Mean object position as of JD2000 36 * \param JD Julian day 37 * \param position Pointer to store new object position as of the given Julian day. 38 * 39 * Calculate equatorial coordinates with the effects of precession for a given Julian Day. 40 * Uses mean equatorial coordinates and is 41 * only for initial epoch J2000.0 42 */ 43 /* Equ 20.3, 20.4 pg 126 44 */ 45 @nogc void ln_get_equ_prec(const ref ln_equ_posn mean_position, double JD, 46 ref ln_equ_posn position) nothrow 47 { 48 real t, t2, t3, A, B, C, zeta, eta, theta, 49 ra, dec, mean_ra, mean_dec; 50 51 /* change original ra and dec to radians */ 52 mean_ra = ln_deg_to_rad(mean_position.ra); 53 mean_dec = ln_deg_to_rad(mean_position.dec); 54 55 /* calc t, zeta, eta and theta for J2000.0 Equ 20.3 */ 56 t =(JD - JD2000) / 36525.0; 57 t *= 1.0 / 3600.0; 58 t2 = t * t; 59 t3 = t2 *t; 60 zeta = 2306.2181 * t + 0.30188 * t2 + 0.017998 * t3; 61 eta = 2306.2181 * t + 1.09468 * t2 + 0.041833 * t3; 62 theta = 2004.3109 * t - 0.42665 * t2 - 0.041833 * t3; 63 zeta = ln_deg_to_rad(zeta); 64 eta = ln_deg_to_rad(eta); 65 theta = ln_deg_to_rad(theta); 66 67 /* calc A,B,C equ 20.4 */ 68 A = cos(mean_dec) * sin(mean_ra + zeta); 69 B = cos(theta) * cos(mean_dec) * 70 cos(mean_ra + zeta) - sin(theta) * sin(mean_dec); 71 C = sin(theta) * cos(mean_dec) * 72 cos(mean_ra + zeta) + cos(theta) * sin(mean_dec); 73 74 ra = atan2(A, B) + eta; 75 76 /* check for object near celestial pole */ 77 if (mean_dec > (0.4 * PI) || mean_dec < (-0.4 * PI)) { 78 /* close to pole */ 79 dec = acos(sqrt(A * A + B * B)); 80 if (mean_dec < 0.) 81 dec *= -1; /* 0 <= acos() <= PI */ 82 } else { 83 /* not close to pole */ 84 dec = asin(C); 85 } 86 87 /* change to degrees */ 88 position.ra = ln_range_degrees(ln_rad_to_deg(ra)); 89 position.dec = ln_rad_to_deg(dec); 90 } 91 92 /*! \fn void ln_get_equ_prec2(struct ln_equ_posn *mean_position, double fromJD, double toJD, struct ln_equ_posn *position); 93 * 94 * \param mean_position Mean object position as of fromJD 95 * \param fromJD Julian day (start) 96 * \param toJD Julian day (end) 97 * \param position Pointer to store new object position as of toJD. 98 * 99 * Calculate the effects of precession on equatorial coordinates, between arbitary Jxxxx epochs. 100 * Use fromJD and toJD parameters to specify required Jxxxx epochs. 101 */ 102 103 /* Equ 20.2, 20.4 pg 126 */ 104 @nogc void ln_get_equ_prec2(const ref ln_equ_posn mean_position, double fromJD, 105 double toJD, ref ln_equ_posn position) nothrow 106 { 107 real t, t2, t3, A, B, C, zeta, eta, theta, ra, dec, mean_ra, 108 mean_dec, T, T2; 109 110 /* change original ra and dec to radians */ 111 mean_ra = ln_deg_to_rad(mean_position.ra); 112 mean_dec = ln_deg_to_rad(mean_position.dec); 113 114 /* calc t, T, zeta, eta and theta Equ 20.2 */ 115 T = (cast(real) (fromJD - JD2000)) / 36525.0; 116 T *= 1.0 / 3600.0; 117 t = (cast(real) (toJD - fromJD)) / 36525.0; 118 t *= 1.0 / 3600.0; 119 T2 = T * T; 120 t2 = t * t; 121 t3 = t2 *t; 122 zeta = (2306.2181 + 1.39656 * T - 0.000139 * T2) * 123 t + (0.30188 - 0.000344 * T) * t2 + 0.017998 * t3; 124 eta = (2306.2181 + 1.39656 * T - 0.000139 * T2) * 125 t + (1.09468 + 0.000066 * T) * t2 + 0.018203 * t3; 126 theta = (2004.3109 - 0.85330 * T - 0.000217 * T2) * 127 t - (0.42665 + 0.000217 * T) * t2 - 0.041833 * t3; 128 zeta = ln_deg_to_rad(zeta); 129 eta = ln_deg_to_rad(eta); 130 theta = ln_deg_to_rad(theta); 131 132 /* calc A,B,C equ 20.4 */ 133 A = cos(mean_dec) * sin(mean_ra + zeta); 134 B = cos(theta) * cos(mean_dec) * cos(mean_ra + zeta) - 135 sin(theta) * sin(mean_dec); 136 C = sin(theta) * cos(mean_dec) * cos(mean_ra + zeta) + 137 cos(theta) * sin(mean_dec); 138 139 ra = atan2(A, B) + eta; 140 141 /* check for object near celestial pole */ 142 if (mean_dec > (0.4 * PI) || mean_dec < (-0.4 * PI)) { 143 /* close to pole */ 144 dec = acos(sqrt(A * A + B * B)); 145 if (mean_dec < 0.) 146 dec *= -1; /* 0 <= acos() <= PI */ 147 } else { 148 /* not close to pole */ 149 dec = asin(C); 150 } 151 152 /* change to degrees */ 153 position.ra = ln_range_degrees(ln_rad_to_deg(ra)); 154 position.dec = ln_rad_to_deg(dec); 155 } 156 157 /*! \fn void ln_get_ecl_prec(struct ln_lnlat_posn *mean_position, double JD, struct ln_lnlat_posn *position) 158 * \param mean_position Mean object position 159 * \param JD Julian day 160 * \param position Pointer to store new object position. 161 * 162 * Calculate ecliptical coordinates with the effects of precession for a given Julian Day. 163 * Uses mean ecliptical coordinates and is 164 * only for initial epoch J2000.0 165 * \todo To be implemented. 166 */ 167 /* Equ 20.5, 20.6 pg 128 168 */ 169 @nogc void ln_get_ecl_prec(const ln_lnlat_posn *mean_position, double JD, 170 ln_lnlat_posn *position) nothrow 171 { 172 } 173 174 }