source: Sophya/trunk/SophyaExt/XephemAstroLib/aberration.c@ 3007

Last change on this file since 3007 was 2818, checked in by cmv, 20 years ago

Update de Xephem 3.7 cmv 21/08/2005
File size: 3.3 KB
Line 
1/* aberration, Jean Meeus, "Astronomical Algorithms", Willman-Bell, 1995;
2 * based on secular unperturbed Kepler orbit
3 *
4 * the corrections should be applied to ra/dec and lam/beta at the
5 * epoch of date.
6 */
7
8#include <stdio.h>
9#include <math.h>
10#include <stdlib.h>
11
12#include "astro.h"
13
14#define ABERR_CONST (20.49552/3600./180.*PI) /* aberr const in rad */
15#define AB_ECL_EOD 0
16#define AB_EQ_EOD 1
17
18static void ab_aux (double mj, double *x, double *y, double lsn, int mode);
19
20/* apply aberration correction to ecliptical coordinates *lam and *bet
21 * (in radians) for a given time m and handily supplied longitude of sun,
22 * lsn (in radians)
23 */
24void
25ab_ecl (double mj, double lsn, double *lam, double *bet)
26{
27 ab_aux(mj, lam, bet, lsn, AB_ECL_EOD);
28}
29
30/* apply aberration correction to equatoreal coordinates *ra and *dec
31 * (in radians) for a given time m and handily supplied longitude of sun,
32 * lsn (in radians)
33 */
34void
35ab_eq (double mj, double lsn, double *ra, double *dec)
36{
37 ab_aux(mj, ra, dec, lsn, AB_EQ_EOD);
38}
39
40/* because the e-terms are secular, keep the real transformation for both
41 * coordinate systems in here with the secular variables cached.
42 * mode == AB_ECL_EOD: x = lam, y = bet (ecliptical)
43 * mode == AB_EQ_EOD: x = ra, y = dec (equatoreal)
44 */
45static void
46ab_aux (double mj, double *x, double *y, double lsn, int mode)
47{
48 static double lastmj = -10000;
49 static double eexc; /* earth orbit excentricity */
50 static double leperi; /* ... and longitude of perihelion */
51 static char dirty = 1; /* flag for cached trig terms */
52
53 if (mj != lastmj) {
54 double T; /* centuries since J2000 */
55
56 T = (mj - J2000)/36525.;
57 eexc = 0.016708617 - (42.037e-6 + 0.1236e-6 * T) * T;
58 leperi = degrad(102.93735 + (0.71953 + 0.00046 * T) * T);
59 lastmj = mj;
60 dirty = 1;
61 }
62
63 switch (mode) {
64 case AB_ECL_EOD: /* ecliptical coords */
65 {
66 double *lam = x, *bet = y;
67 double dlsun, dlperi;
68
69 dlsun = lsn - *lam;
70 dlperi = leperi - *lam;
71
72 /* valid only for *bet != +-PI/2 */
73 *lam -= ABERR_CONST/cos(*bet) * (cos(dlsun) -
74 eexc*cos(dlperi));
75 *bet -= ABERR_CONST*sin(*bet) * (sin(dlsun) -
76 eexc*sin(dlperi));
77 }
78 break;
79
80 case AB_EQ_EOD: /* equatoreal coords */
81 {
82 double *ra = x, *dec = y;
83 double sr, cr, sd, cd, sls, cls;/* trig values coords */
84 static double cp, sp, ce, se; /* .. and perihel/eclipic */
85 double dra, ddec; /* changes in ra and dec */
86
87 if (dirty) {
88 double eps;
89
90 cp = cos(leperi);
91 sp = sin(leperi);
92 obliquity(mj, &eps);
93 se = sin(eps);
94 ce = cos(eps);
95 dirty = 0;
96 }
97
98 sr = sin(*ra);
99 cr = cos(*ra);
100 sd = sin(*dec);
101 cd = cos(*dec);
102 sls = sin(lsn);
103 cls = cos(lsn);
104
105 dra = ABERR_CONST/cd * ( -(cr * cls * ce + sr * sls) +
106 eexc * (cr * cp * ce + sr * sp));
107
108 ddec = se/ce * cd - sr * sd; /* tmp use */
109 ddec = ABERR_CONST * ( -(cls * ce * ddec + cr * sd * sls) +
110 eexc * (cp * ce * ddec + cr * sd * sp) );
111
112 *ra += dra;
113 range (ra, 2*PI);
114 *dec += ddec;
115 }
116 break;
117
118 default:
119 printf ("ab_aux: bad mode: %d\n", mode);
120 abort();
121 break;
122
123 } /* switch (mode) */
124}
125
126/* For RCS Only -- Do Not Edit */
127static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: aberration.c,v $ $Date: 2005-08-21 10:02:36 $ $Revision: 1.5 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.