3 Common routines used by all drivers
4 Copyright (C) 2003 by Jason Harris (jharris@30doradus.org)
7 This is the C version of the astronomical library in KStars
8 modified by Jasem Mutlaq (mutlaqja@ikarustech.com)
10 This library is free software; you can redistribute it and/or
11 modify it under the terms of the GNU Lesser General Public
12 License as published by the Free Software Foundation; either
13 version 2.1 of the License, or (at your option) any later version.
15 This library is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 Lesser General Public License for more details.
20 You should have received a copy of the GNU Lesser General Public
21 License along with this library; if not, write to the Free Software
22 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 /* needed for sincos() in math.h */
36 const char * Direction
[] = { "North", "West", "East", "South", "All"};
37 const char * SolarSystem
[] = { "Mercury", "Venus", "Moon", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"};
39 /* make it compile on solaris */
41 #define M_PI 3.14159265358979323846264338327950288419716939937510582097494459
44 /******** Prototypes ***********/
45 double DegToRad( double num
);
46 double RadToDeg( double num
);
47 void SinCos( double Degrees
, double *sina
, double *cosa
);
49 double obliquity(void);
50 double constAberr(void);
51 double sunMeanAnomaly(void);
52 double sunMeanLongitude(void);
53 double sunTrueAnomaly(void);
54 double sunTrueLongitude(void);
55 double earthPerihelionLongitude(void);
56 double earthEccentricity(void);
59 double julianCenturies(void);
60 double p1( int i1
, int i2
);
61 double p2( int i1
, int i2
);
62 void updateAstroValues( double jd
);
63 void nutate(double *RA
, double *Dec
);
64 void aberrate(double *RA
, double *Dec
);
65 void precessFromAnyEpoch(double jd0
, double jdf
, double *RA
, double *Dec
);
66 void apparentCoord(double jd0
, double jdf
, double *RA
, double *Dec
);
67 void getSexComponents(double value
, int *d
, int *m
, int *s
);
69 double K
= ( 20.49552 / 3600. );
70 double Obliquity
, K
, L
, L0
, LM
, M
, M0
, O
, P
;
72 double CX
, SX
, CY
, SY
, CZ
, SZ
;
73 double P1
[3][3], P2
[3][3];
74 double deltaObliquity
, deltaEcLong
;
77 double obliquity() { return Obliquity
; }
78 double constAberr() { return K
; }
79 double sunMeanAnomaly() { return M
; }
80 double sunMeanLongitude() { return L
; }
81 double sunTrueAnomaly() { return M0
; }
82 double sunTrueLongitude() { return L0
; }
83 double earthPerihelionLongitude() { return P
; }
84 double earthEccentricity() { return e
; }
85 double dObliq() { return deltaObliquity
; }
86 double dEcLong() { return deltaEcLong
; }
87 double julianCenturies() { return T
; }
88 double p1( int i1
, int i2
) { return P1
[i1
][i2
]; }
89 double p2( int i1
, int i2
) { return P2
[i1
][i2
]; }
91 void updateAstroValues( double jd
)
93 static double days
= 0;
96 /* Nutation parameters */
98 double sin2L
, cos2L
, sin2M
, cos2M
;
99 double sinO
, cosO
, sin2O
, cos2O
;
101 /* no need to recalculate if same as previous JD */
107 /* Julian Centuries since J2000.0 */
108 T
= ( jd
- J2000
) / 36525.;
110 /* Julian Millenia since J2000.0 */
113 /* Sun's Mean Longitude */
114 L
= ( 280.46645 + 36000.76983*T
+ 0.0003032*T
*T
);
116 /* Sun's Mean Anomaly */
117 M
= ( 357.52910 + 35999.05030*T
- 0.0001559*T
*T
- 0.00000048*T
*T
*T
);
119 /* Moon's Mean Longitude */
120 LM
= ( 218.3164591 + 481267.88134236*T
- 0.0013268*T
*T
+ T
*T
*T
/538841. - T
*T
*T
*T
/6519400.);
122 /* Longitude of Moon's Ascending Node */
123 O
= ( 125.04452 - 1934.136261*T
+ 0.0020708*T
*T
+ T
*T
*T
/450000.0 );
125 /* Earth's orbital eccentricity */
126 e
= 0.016708617 - 0.000042037*T
- 0.0000001236*T
*T
;
128 C
= ( 1.914600 - 0.004817*T
- 0.000014*T
*T
) * sin( DegToRad(M
) )
129 + ( 0.019993 - 0.000101*T
) * sin( 2.0* DegToRad(M
) )
130 + 0.000290 * sin( 3.0* DegToRad(M
));
132 /* Sun's True Longitude */
135 /* Sun's True Anomaly */
138 /* Obliquity of the Ecliptic */
140 dObliq2
= -4680.93*U
- 1.55*U
*U
+ 1999.25*U
*U
*U
141 - 51.38*U
*U
*U
*U
- 249.67*U
*U
*U
*U
*U
142 - 39.05*U
*U
*U
*U
*U
*U
+ 7.12*U
*U
*U
*U
*U
*U
*U
143 + 27.87*U
*U
*U
*U
*U
*U
*U
*U
+ 5.79*U
*U
*U
*U
*U
*U
*U
*U
*U
144 + 2.45*U
*U
*U
*U
*U
*U
*U
*U
*U
*U
;
145 Obliquity
= ( 23.43929111 + dObliq2
/3600.0);
148 L2
= ( 2.0 * L
); /* twice mean ecl. long. of Sun */
149 M2
= ( 2.0 * LM
); /* twice mean ecl. long. of Moon */
151 SinCos( O
, &sinO
, &cosO
);
152 SinCos( O2
, &sin2O
, &cos2O
);
153 SinCos( L2
, &sin2L
, &cos2L
);
154 SinCos( M2
, &sin2M
, &cos2M
);
156 deltaEcLong
= ( -17.2*sinO
- 1.32*sin2L
- 0.23*sin2M
+ 0.21*sin2O
)/3600.0; /* Ecl. long. correction */
157 deltaObliquity
= ( 9.2*cosO
+ 0.57*cos2L
+ 0.10*cos2M
- 0.09*cos2O
)/3600.0; /* Obliq. correction */
159 /* Compute Precession Matrices: */
160 XP
= ( 0.6406161*T
+ 0.0000839*T
*T
+ 0.0000050*T
*T
*T
);
161 YP
= ( 0.5567530*T
- 0.0001185*T
*T
- 0.0000116*T
*T
*T
);
162 ZP
= ( 0.6406161*T
+ 0.0003041*T
*T
+ 0.0000051*T
*T
*T
);
164 SinCos(XP
, &SX
, &CX
);
165 SinCos(YP
, &SY
, &CY
);
166 SinCos(ZP
, &SZ
, &CZ
);
168 /* P1 is used to precess from any epoch to J2000 */
169 P1
[0][0] = CX
*CY
*CZ
- SX
*SZ
;
170 P1
[1][0] = CX
*CY
*SZ
+ SX
*CZ
;
172 P1
[0][1] = -1.0*SX
*CY
*CZ
- CX
*SZ
;
173 P1
[1][1] = -1.0*SX
*CY
*SZ
+ CX
*CZ
;
174 P1
[2][1] = -1.0*SX
*SY
;
175 P1
[0][2] = -1.0*SY
*CZ
;
176 P1
[1][2] = -1.0*SY
*SZ
;
179 /* P2 is used to precess from J2000 to any other epoch (it is the transpose of P1) */
180 P2
[0][0] = CX
*CY
*CZ
- SX
*SZ
;
181 P2
[1][0] = -1.0*SX
*CY
*CZ
- CX
*SZ
;
182 P2
[2][0] = -1.0*SY
*CZ
;
183 P2
[0][1] = CX
*CY
*SZ
+ SX
*CZ
;
184 P2
[1][1] = -1.0*SX
*CY
*SZ
+ CX
*CZ
;
185 P2
[2][1] = -1.0*SY
*SZ
;
187 P2
[1][2] = -1.0*SX
*SY
;
192 double DegToRad( double num
)
194 /*return (num * deg_rad);*/
195 return (num
* (M_PI
/ 180.0));
198 double RadToDeg( double num
)
200 return (num
/ (M_PI
/ 180.0));
203 void SinCos( double Degrees
, double *sina
, double *cosa
)
205 /**We have two versions of this function. One is ANSI standard, but
206 *slower. The other is faster, but is GNU only.
207 *Using the __GLIBC_ and __GLIBC_MINOR_ constants to determine if the
208 * GNU extension sincos() is defined.
210 static int rDirty
= 1;
217 #if ( __GLIBC__ >= 2 && __GLIBC_MINOR__ >=1 )
219 sincos( DegToRad(Degrees
), &Sin
, &Cos
);
221 /* For older GLIBC versions */
222 Sin
= ::sin( DegToRad(Degrees
) );
223 Cos
= ::cos( DegToRad(Degrees
) );
226 /* ANSI-compliant version */
227 Sin
= sin( DegToRad(Degrees
) );
228 Cos
= cos( DegToRad(Degrees
) );
234 Sin
= sin( DegToRad(Degrees
) );
235 Cos
= cos( DegToRad(Degrees
) );
242 double calculateDec(double latitude
, double SDTime
)
244 if (SDTime
> 12) SDTime
-= 12;
246 return RadToDeg ( atan ( (cos (DegToRad (latitude
)) / sin (DegToRad (latitude
))) * cos (DegToRad (SDTime
))) );
250 double calculateRA(double SDTime
)
263 void nutate(double *RA
, double *Dec
)
265 double cosRA
, sinRA
, cosDec
, sinDec
, tanDec
;
269 SinCos( *RA
, &sinRA
, &cosRA
);
270 SinCos( *Dec
, &sinDec
, &cosDec
);
272 SinCos( obliquity(), &sinOb
, &cosOb
);
274 /* Step 2: Nutation */
275 /* approximate method */
276 tanDec
= sinDec
/cosDec
;
278 dRA1
= ( dEcLong()*( cosOb
+ sinOb
*sinRA
*tanDec
) - dObliq()*cosRA
*tanDec
);
279 dDec1
= ( dEcLong()*( sinOb
*cosRA
) + dObliq()*sinRA
);
281 *RA
= ( *RA
+ dRA1
);
282 *Dec
= ( *Dec
+ dDec1
);
285 void aberrate(double *RA
, double *Dec
)
287 double cosRA
, sinRA
, cosDec
, sinDec
;
289 double cosOb
, sinOb
, cosL
, sinL
, cosP
, sinP
;
290 double K2
, e2
, tanOb
;
292 K2
= constAberr(); /* constant of aberration */
293 e2
= earthEccentricity();
295 SinCos( *RA
, &sinRA
, &cosRA
);
296 SinCos( *Dec
, &sinDec
, &cosDec
);
298 SinCos( obliquity(), &sinOb
, &cosOb
);
301 SinCos( sunTrueLongitude(), &sinL
, &cosL
);
302 SinCos( earthPerihelionLongitude(), &sinP
, &cosP
);
304 /* Step 3: Aberration */
305 dRA2
= ( -1.0 * K2
* ( cosRA
* cosL
* cosOb
+ sinRA
* sinL
)/cosDec
306 + e2
* K2
* ( cosRA
* cosP
* cosOb
+ sinRA
* sinP
)/cosDec
);
308 dDec2
= ( -1.0 * K2
* ( cosL
* cosOb
* ( tanOb
* cosDec
- sinRA
* sinDec
) + cosRA
* sinDec
* sinL
)
309 + e2
* K2
* ( cosP
* cosOb
* ( tanOb
* cosDec
- sinRA
* sinDec
) + cosRA
* sinDec
* sinP
) );
311 *RA
= ( *RA
+ dRA2
);
312 *Dec
= ( *Dec
+ dDec2
);
315 void precessFromAnyEpoch(double jd0
, double jdf
, double *RA
, double *Dec
)
317 double cosRA0
, sinRA0
, cosDec0
, sinDec0
;
321 SinCos( *RA
, &sinRA0
, &cosRA0
);
322 SinCos( *Dec
, &sinDec0
, &cosDec0
);
324 /* Need first to precess to J2000.0 coords */
329 /* v is a column vector representing input coordinates. */
331 updateAstroValues(jd0
);
333 v
[0] = cosRA0
*cosDec0
;
334 v
[1] = sinRA0
*cosDec0
;
338 /*s is the product of P1 and v; s represents the
339 coordinates precessed to J2000 */
340 for ( i
=0; i
<3; ++i
) {
341 s
[i
] = p1( 0, i
)*v
[0] + p1( 1, i
)*v
[1] +
348 /* Input coords already in J2000, set s accordingly. */
349 s
[0] = cosRA0
*cosDec0
;
350 s
[1] = sinRA0
*cosDec0
;
354 updateAstroValues(jdf
);
356 /* Multiply P2 and s to get v, the vector representing the new coords. */
358 for ( i
=0; i
<3; ++i
) {
359 v
[i
] = p2( 0, i
)*s
[0] + p2( 1, i
)*s
[1] +
363 /*Extract RA, Dec from the vector:
364 RA.setRadians( atan( v[1]/v[0] ) ); */
366 *RA
= RadToDeg( atan2( v
[1],v
[0] ) );
367 *Dec
= RadToDeg( asin( v
[2] ) );
370 *RA
= ( *RA
+ 360.0 );
373 void apparentCoord(double jd0
, double jdf
, double *RA
, double *Dec
)
377 precessFromAnyEpoch(jd0
,jdf
, RA
, Dec
);
379 updateAstroValues(jdf
);
387 double UTtoJD(struct tm
*utm
)
389 int year
, month
, day
, hour
, minute
, second
;
390 int m
, y
, A
, B
, C
, D
;
393 /* Note: The tm_year was modified by adding +1900 to it since tm_year refers
394 to the number of years after 1900. The month field was also modified by adding 1 to it
395 since the tm_mon range is from 0 to 11 */
400 minute
= utm
->tm_min
;
401 second
= utm
->tm_sec
;
414 /* If the date is after 10/15/1582, then take Pope Gregory's modification
415 to the Julian calendar into account */
417 if (( year
>1582 ) ||
418 ( year
==1582 && month
>9 ) ||
419 ( year
==1582 && month
==9 && day
>15 ))
422 B
= 2 - A
+ (int) A
/4;
428 C
= (int) ((365.25*y
) - 0.75);
430 C
= (int) (365.25*y
);
433 D
= (int) (30.6001*(m
+1));
435 d
= (double) day
+ ( (double) hour
+ ( (double) minute
+ (double) second
/60.0)/60.0)/24.0;
436 jd
= B
+ C
+ D
+ d
+ 1720994.5;
441 double JDtoGMST( double jd
)
445 /* Julian Centuries since J2000.0 */
446 T
= ( jd
- J2000
) / 36525.;
448 /* Greewich Mean Sidereal Time */
451 + 360.98564736629*(jd
-J2000
)
458 int extractISOTime(char *timestr
, struct tm
*utm
)
461 if (strptime(timestr
, "%Y-%m-%dT%H:%M:%S", utm
))
463 if (strptime(timestr
, "%Y/%m/%dT%H:%M:%S", utm
))
465 if (strptime(timestr
, "%Y:%m:%dT%H:%M:%S", utm
))
467 if (strptime(timestr
, "%Y-%m-%dT%H-%M-%S", utm
))
474 /* sprint the variable a in sexagesimal format into out[].
475 * w is the number of spaces for the whole part.
476 * fracbase is the number of pieces a whole is to broken into; valid options:
477 * 360000: <w>:mm:ss.ss
482 * return number of characters written to out, not counting final '\0'.
485 fs_sexa (char *out
, double a
, int w
, int fracbase
)
495 /* save whether it's negative but do all the rest with a positive */
500 /* convert to an integral number of whole portions */
501 n
= (unsigned long)(a
* fracbase
+ 0.5);
505 /* form the whole part; "negative 0" is a special case */
507 out
+= sprintf (out
, "%*s-0", w
-2, "");
509 out
+= sprintf (out
, "%*d", w
, isneg
? -d
: d
);
515 out
+= sprintf (out
, ":%02d", m
);
517 case 600: /* dd:mm.m */
518 out
+= sprintf (out
, ":%02d.%1d", f
/10, f
%10);
520 case 3600: /* dd:mm:ss */
523 out
+= sprintf (out
, ":%02d:%02d", m
, s
);
525 case 36000: /* dd:mm:ss.s*/
528 out
+= sprintf (out
, ":%02d:%02d.%1d", m
, s
/10, s
%10);
530 case 360000: /* dd:mm:ss.ss */
533 out
+= sprintf (out
, ":%02d:%02d.%02d", m
, s
/100, s
%100);
536 printf ("fs_sexa: unknown fracbase: %d\n", fracbase
);
543 /* convert sexagesimal string str AxBxC to double.
544 * x can be anything non-numeric. Any missing A, B or C will be assumed 0.
545 * optional - and + can be anywhere.
546 * return 0 if ok, -1 if can't find a thing.
550 const char *str0
, /* input string */
551 double *dp
) /* cracked value, if return 0 */
553 double a
= 0, b
= 0, c
= 0;
558 /* copy str0 so we can play with it */
559 strncpy (str
, str0
, sizeof(str
)-1);
560 str
[sizeof(str
)-1] = '\0';
562 neg
= strchr(str
, '-');
565 r
= sscanf (str
, "%lf%*[^0-9]%lf%*[^0-9]%lf", &a
, &b
, &c
);
568 *dp
= a
+ b
/60 + c
/3600;
574 void getSexComponents(double value
, int *d
, int *m
, int *s
)
577 *d
= (int) fabs(value
);
578 *m
= (int) ((fabs(value
) - *d
) * 60.0);
579 *s
= (int) rint(((fabs(value
) - *d
) * 60.0 - *m
) *60.0);
585 /* fill buf with properly formatted INumber string. return length */
587 numberFormat (char *buf
, const char *format
, double value
)
592 if (sscanf (format
, "%%%d.%d%c", &w
, &f
, &m
) == 3 && m
== 'm') {
593 /* INDI sexi format */
595 case 9: s
= 360000; break;
596 case 8: s
= 36000; break;
597 case 6: s
= 3600; break;
598 case 5: s
= 600; break;
599 default: s
= 60; break;
601 return (fs_sexa (buf
, value
, w
-f
, s
));
603 /* normal printf format */
604 return (sprintf (buf
, format
, value
));
608 double angularDistance(double fromRA
, double fromDEC
, double toRA
, double toDEC
)
611 double dalpha
= DegToRad(fromRA
) - DegToRad(toRA
);
612 double ddelta
= DegToRad(fromDEC
) - DegToRad(toDEC
);
614 double sa
= sin(dalpha
/2.);
615 double sd
= sin(ddelta
/2.);
620 double aux
= havd
+ cos (DegToRad(fromDEC
)) * cos(DegToRad(toDEC
)) * hava
;
622 return (RadToDeg ( 2 * fabs(asin( sqrt(aux
) ))));
625 /* return current system time in message format */
635 strftime (ts
, sizeof(ts
), "%Y-%m-%dT%H:%M:%S", tp
);