kstars/kstars/indi/indicom.c

   1 /*
   2     INDI LIB
   3     Common routines used by all drivers
   4     Copyright (C) 2003 by Jason Harris (jharris@30doradus.org)
   5                           Elwood C. Downey
   6
   7     This is the C version of the astronomical library in KStars
   8     modified by Jasem Mutlaq (mutlaqja@ikarustech.com)
   9
  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.
  14
  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.
  19
  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
  23
  24 */
  25
  26 /* needed for sincos() in math.h */
  27 #define _GNU_SOURCE
  28
  29 #include <stdlib.h>
  30 #include <math.h>
  31 #include <stdio.h>
  32 #include <string.h>
  33
  34 #include "indicom.h"
  35
  36 const char * Direction[] = { "North", "West", "East", "South", "All"};
  37 const char * SolarSystem[] = { "Mercury", "Venus", "Moon", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"};
  38
  39 /* make it compile on solaris */
  40 #ifndef M_PI
  41 #define M_PI 3.14159265358979323846264338327950288419716939937510582097494459
  42 #endif
  43
  44 /******** Prototypes ***********/
  45 double DegToRad( double num );
  46 double RadToDeg( double num );
  47 void SinCos( double Degrees, double *sina, double *cosa );
  48
  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);
  57 double dObliq(void);
  58 double dEcLong(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);
  68
  69 double K = ( 20.49552 / 3600. );
  70 double Obliquity, K, L, L0, LM, M, M0, O, P;
  71 double XP, YP, ZP;
  72 double CX, SX, CY, SY, CZ, SZ;
  73 double P1[3][3], P2[3][3];
  74 double deltaObliquity, deltaEcLong;
  75 double e, T;
  76
  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]; }
  90
  91 void updateAstroValues( double jd )
  92 {
  93         static double days = 0;
  94         double jm;
  95         double C, U, dObliq2;
  96         /* Nutation parameters */
  97         double  L2, M2, O2;
  98         double sin2L, cos2L, sin2M, cos2M;
  99         double sinO, cosO, sin2O, cos2O;
 100
 101         /* no need to recalculate if same as previous JD */
 102         if (days == jd)
 103          return;
 104
 105         days = jd;
 106
 107         /* Julian Centuries since J2000.0 */
 108         T = ( jd - J2000 ) / 36525.;
 109
 110         /* Julian Millenia since J2000.0 */
 111         jm = T / 10.0;
 112
 113         /* Sun's Mean Longitude */
 114         L = ( 280.46645 + 36000.76983*T + 0.0003032*T*T );
 115
 116         /* Sun's Mean Anomaly */
 117         M = ( 357.52910 + 35999.05030*T - 0.0001559*T*T - 0.00000048*T*T*T );
 118
 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.);
 121
 122         /* Longitude of Moon's Ascending Node */
 123         O = ( 125.04452 - 1934.136261*T + 0.0020708*T*T + T*T*T/450000.0 );
 124
 125         /* Earth's orbital eccentricity */
 126         e = 0.016708617 - 0.000042037*T - 0.0000001236*T*T;
 127
 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));
 131
 132         /* Sun's True Longitude */
 133         L0 = ( L + C );
 134
 135         /* Sun's True Anomaly */
 136         M0 = ( M + C );
 137
 138         /* Obliquity of the Ecliptic */
 139         U = T/100.0;
 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);
 146
 147         O2 = ( 2.0 * O );
 148         L2 = ( 2.0 * L );  /* twice mean ecl. long. of Sun */
 149         M2 = ( 2.0 * LM);  /* twice mean ecl. long. of Moon */
 150
 151         SinCos( O, &sinO, &cosO );
 152         SinCos( O2, &sin2O, &cos2O );
 153         SinCos( L2, &sin2L, &cos2L );
 154         SinCos( M2, &sin2M, &cos2M );
 155
 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 */
 158
 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 );
 163
 164         SinCos(XP, &SX, &CX );
 165         SinCos(YP, &SY, &CY );
 166         SinCos(ZP, &SZ, &CZ );
 167
 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;
 171         P1[2][0] = CX*SY;
 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;
 177         P1[2][2] = CY;
 178
 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;
 186         P2[0][2] = CX*SY;
 187         P2[1][2] = -1.0*SX*SY;
 188         P2[2][2] = CY;
 189
 190 }
 191
 192 double DegToRad( double num )
 193 {
 194   /*return (num * deg_rad);*/
 195   return (num * (M_PI / 180.0));
 196 }
 197
 198 double RadToDeg( double num )
 199 {
 200   return (num / (M_PI / 180.0));
 201 }
 202
 203 void SinCos( double Degrees, double *sina, double *cosa )
 204 {
 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.
 209                 */
 210                 static int rDirty = 1;
 211                 double Sin, Cos;
 212
 213
 214                 if (rDirty)
 215                 {
 216                 #ifdef __GLIBC__
 217                 #if ( __GLIBC__ >= 2 && __GLIBC_MINOR__ >=1 )
 218                 /* GNU version */
 219                 sincos( DegToRad(Degrees), &Sin, &Cos );
 220                 #else
 221                 /* For older GLIBC versions */
 222                 Sin = ::sin( DegToRad(Degrees) );
 223                 Cos = ::cos( DegToRad(Degrees) );
 224                 #endif
 225                 #else
 226                 /* ANSI-compliant version */
 227                 Sin = sin( DegToRad(Degrees) );
 228                 Cos = cos( DegToRad(Degrees) );
 229                 rDirty = 0;
 230                 #endif
 231                 }
 232                 else
 233                 {
 234                  Sin = sin( DegToRad(Degrees) );
 235                  Cos = cos( DegToRad(Degrees) );
 236                 }
 237
 238                 *sina = Sin;
 239                 *cosa = Cos;
 240 }
 241
 242 double calculateDec(double latitude, double SDTime)
 243 {
 244   if (SDTime > 12) SDTime -= 12;
 245
 246   return RadToDeg ( atan ( (cos (DegToRad (latitude)) / sin (DegToRad (latitude))) * cos (DegToRad (SDTime))) );
 247
 248 }
 249
 250 double calculateRA(double SDTime)
 251 {
 252   double ra;
 253
 254   ra = SDTime + 12;
 255
 256   if (ra < 24)
 257     return ra;
 258   else
 259     return (ra - 24);
 260 }
 261
 262
 263 void nutate(double *RA, double *Dec)
 264 {
 265         double cosRA, sinRA, cosDec, sinDec, tanDec;
 266         double dRA1, dDec1;
 267         double cosOb, sinOb;
 268
 269         SinCos( *RA, &sinRA, &cosRA );
 270         SinCos( *Dec, &sinDec, &cosDec );
 271
 272         SinCos( obliquity(), &sinOb, &cosOb );
 273
 274         /* Step 2: Nutation */
 275         /* approximate method */
 276         tanDec = sinDec/cosDec;
 277
 278         dRA1  = ( dEcLong()*( cosOb + sinOb*sinRA*tanDec ) - dObliq()*cosRA*tanDec );
 279         dDec1 = ( dEcLong()*( sinOb*cosRA ) + dObliq()*sinRA );
 280
 281         *RA  = ( *RA + dRA1 );
 282         *Dec = ( *Dec+ dDec1);
 283 }
 284
 285 void aberrate(double *RA, double *Dec)
 286 {
 287         double cosRA, sinRA, cosDec, sinDec;
 288         double dRA2, dDec2;
 289         double cosOb, sinOb, cosL, sinL, cosP, sinP;
 290         double K2, e2, tanOb;
 291
 292         K2 = constAberr();  /* constant of aberration */
 293         e2 = earthEccentricity();
 294
 295         SinCos( *RA, &sinRA, &cosRA );
 296         SinCos( *Dec, &sinDec, &cosDec );
 297
 298         SinCos( obliquity(), &sinOb, &cosOb );
 299         tanOb = sinOb/cosOb;
 300
 301         SinCos( sunTrueLongitude(), &sinL, &cosL );
 302         SinCos( earthPerihelionLongitude(), &sinP, &cosP );
 303
 304         /* Step 3: Aberration */
 305         dRA2 = ( -1.0 * K2 * ( cosRA * cosL * cosOb + sinRA * sinL )/cosDec
 306                 + e2 * K2 * ( cosRA * cosP * cosOb + sinRA * sinP )/cosDec );
 307
 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 ) );
 310
 311         *RA  = ( *RA  + dRA2 );
 312         *Dec = ( *Dec + dDec2);
 313 }
 314
 315 void precessFromAnyEpoch(double jd0, double jdf, double *RA, double *Dec)
 316 {
 317         double cosRA0, sinRA0, cosDec0, sinDec0;
 318         double v[3], s[3];
 319         unsigned int i=0;
 320
 321         SinCos( *RA, &sinRA0, &cosRA0 );
 322         SinCos( *Dec, &sinDec0, &cosDec0 );
 323
 324         /* Need first to precess to J2000.0 coords */
 325
 326         if ( jd0 != J2000 )
 327         {
 328
 329         /* v is a column vector representing input coordinates. */
 330
 331                 updateAstroValues(jd0);
 332
 333                 v[0] = cosRA0*cosDec0;
 334                 v[1] = sinRA0*cosDec0;
 335                 v[2] = sinDec0;
 336
 337
 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] +
 342                                 p1( 2, i )*v[2];
 343                 }
 344
 345         } else
 346         {
 347
 348         /* Input coords already in J2000, set s accordingly. */
 349                 s[0] = cosRA0*cosDec0;
 350                 s[1] = sinRA0*cosDec0;
 351                 s[2] = sinDec0;
 352        }
 353
 354        updateAstroValues(jdf);
 355
 356         /* Multiply P2 and s to get v, the vector representing the new coords. */
 357
 358         for ( i=0; i<3; ++i ) {
 359                 v[i] = p2( 0, i )*s[0] + p2( 1, i )*s[1] +
 360                 p2( 2, i )*s[2];
 361         }
 362
 363         /*Extract RA, Dec from the vector:
 364         RA.setRadians( atan( v[1]/v[0] ) ); */
 365
 366         *RA  = RadToDeg( atan2( v[1],v[0] ) );
 367         *Dec = RadToDeg( asin( v[2] ) );
 368
 369         if (*RA < 0.0 )
 370             *RA = ( *RA + 360.0 );
 371 }
 372
 373 void apparentCoord(double jd0, double jdf, double *RA, double *Dec)
 374 {
 375         *RA = *RA * 15.0;
 376
 377         precessFromAnyEpoch(jd0,jdf, RA, Dec);
 378
 379         updateAstroValues(jdf);
 380
 381         nutate(RA, Dec);
 382         aberrate(RA, Dec);
 383
 384         *RA = *RA / 15.0;
 385 }
 386
 387 double UTtoJD(struct tm *utm)
 388 {
 389   int year, month, day, hour, minute, second;
 390   int m, y, A, B, C, D;
 391   double d, jd;
 392
 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 */
 396   year   = utm->tm_year;
 397   month  = utm->tm_mon;
 398   day    = utm->tm_mday;
 399   hour   = utm->tm_hour;
 400   minute = utm->tm_min;
 401   second = utm->tm_sec;
 402
 403
 404   if (month > 2)
 405   {
 406           m = month;
 407           y = year;
 408   } else
 409   {
 410           y = year - 1;
 411           m = month + 12;
 412   }
 413
 414  /*  If the date is after 10/15/1582, then take Pope Gregory's modification
 415          to the Julian calendar into account */
 416
 417          if (( year  >1582 ) ||
 418                  ( year ==1582 && month  >9 ) ||
 419                  ( year ==1582 && month ==9 && day >15 ))
 420          {
 421                  A = (int) y/100;
 422                  B = 2 - A + (int) A/4;
 423          } else {
 424                  B = 0;
 425          }
 426
 427   if (y < 0) {
 428           C = (int) ((365.25*y) - 0.75);
 429   } else {
 430           C = (int) (365.25*y);
 431   }
 432
 433   D = (int) (30.6001*(m+1));
 434
 435   d = (double) day + ( (double) hour + ( (double) minute + (double) second/60.0)/60.0)/24.0;
 436   jd = B + C + D + d + 1720994.5;
 437
 438   return jd;
 439 }
 440
 441 double JDtoGMST( double jd )
 442 {
 443         double Gmst;
 444
 445         /* Julian Centuries since J2000.0 */
 446         T = ( jd - J2000 ) / 36525.;
 447
 448         /* Greewich Mean Sidereal Time */
 449
 450         Gmst = 280.46061837
 451                 + 360.98564736629*(jd-J2000)
 452                 + 0.000387933*T*T
 453                 - T*T*T/38710000.0;
 454
 455         return Gmst;
 456 }
 457
 458 int extractISOTime(char *timestr, struct tm *utm)
 459 {
 460
 461   if (strptime(timestr, "%Y-%m-%dT%H:%M:%S", utm))
 462    return (0);
 463   if (strptime(timestr, "%Y/%m/%dT%H:%M:%S", utm))
 464    return (0);
 465   if (strptime(timestr, "%Y:%m:%dT%H:%M:%S", utm))
 466    return (0);
 467   if (strptime(timestr, "%Y-%m-%dT%H-%M-%S", utm))
 468    return (0);
 469
 470    return (-1);
 471
 472 }
 473
 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
 478  *      36000:  <w>:mm:ss.s
 479  *      3600:   <w>:mm:ss
 480  *      600:    <w>:mm.m
 481  *      60:     <w>:mm
 482  * return number of characters written to out, not counting final '\0'.
 483  */
 484 int
 485 fs_sexa (char *out, double a, int w, int fracbase)
 486 {
 487         char *out0 = out;
 488         unsigned long n;
 489         int d;
 490         int f;
 491         int m;
 492         int s;
 493         int isneg;
 494
 495         /* save whether it's negative but do all the rest with a positive */
 496         isneg = (a < 0);
 497         if (isneg)
 498             a = -a;
 499
 500         /* convert to an integral number of whole portions */
 501         n = (unsigned long)(a * fracbase + 0.5);
 502         d = n/fracbase;
 503         f = n%fracbase;
 504
 505         /* form the whole part; "negative 0" is a special case */
 506         if (isneg && d == 0)
 507             out += sprintf (out, "%*s-0", w-2, "");
 508         else
 509             out += sprintf (out, "%*d", w, isneg ? -d : d);
 510
 511         /* do the rest */
 512         switch (fracbase) {
 513         case 60:        /* dd:mm */
 514             m = f/(fracbase/60);
 515             out += sprintf (out, ":%02d", m);
 516             break;
 517         case 600:       /* dd:mm.m */
 518             out += sprintf (out, ":%02d.%1d", f/10, f%10);
 519             break;
 520         case 3600:      /* dd:mm:ss */
 521             m = f/(fracbase/60);
 522             s = f%(fracbase/60);
 523             out += sprintf (out, ":%02d:%02d", m, s);
 524             break;
 525         case 36000:     /* dd:mm:ss.s*/
 526             m = f/(fracbase/60);
 527             s = f%(fracbase/60);
 528             out += sprintf (out, ":%02d:%02d.%1d", m, s/10, s%10);
 529             break;
 530         case 360000:    /* dd:mm:ss.ss */
 531             m = f/(fracbase/60);
 532             s = f%(fracbase/60);
 533             out += sprintf (out, ":%02d:%02d.%02d", m, s/100, s%100);
 534             break;
 535         default:
 536             printf ("fs_sexa: unknown fracbase: %d\n", fracbase);
 537             exit(1);
 538         }
 539
 540         return (out - out0);
 541 }
 542
 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.
 547  */
 548 int
 549 f_scansexa (
 550 const char *str0,       /* input string */
 551 double *dp)             /* cracked value, if return 0 */
 552 {
 553         double a = 0, b = 0, c = 0;
 554         char str[128];
 555         char *neg;
 556         int r;
 557
 558         /* copy str0 so we can play with it */
 559         strncpy (str, str0, sizeof(str)-1);
 560         str[sizeof(str)-1] = '\0';
 561
 562         neg = strchr(str, '-');
 563         if (neg)
 564             *neg = ' ';
 565         r = sscanf (str, "%lf%*[^0-9]%lf%*[^0-9]%lf", &a, &b, &c);
 566         if (r < 1)
 567             return (-1);
 568         *dp = a + b/60 + c/3600;
 569         if (neg)
 570             *dp *= -1;
 571         return (0);
 572 }
 573
 574 void getSexComponents(double value, int *d, int *m, int *s)
 575 {
 576
 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);
 580
 581   if (value < 0)
 582    *d *= -1;
 583 }
 584
 585 /* fill buf with properly formatted INumber string. return length */
 586 int
 587 numberFormat (char *buf, const char *format, double value)
 588 {
 589         int w, f, s;
 590         char m;
 591
 592         if (sscanf (format, "%%%d.%d%c", &w, &f, &m) == 3 && m == 'm') {
 593             /* INDI sexi format */
 594             switch (f) {
 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;
 600             }
 601             return (fs_sexa (buf, value, w-f, s));
 602         } else {
 603             /* normal printf format */
 604             return (sprintf (buf, format, value));
 605         }
 606 }
 607
 608 double angularDistance(double fromRA, double fromDEC, double toRA, double toDEC)
 609 {
 610
 611         double dalpha = DegToRad(fromRA) - DegToRad(toRA);
 612         double ddelta = DegToRad(fromDEC) - DegToRad(toDEC);
 613
 614         double sa = sin(dalpha/2.);
 615         double sd = sin(ddelta/2.);
 616
 617         double hava = sa*sa;
 618         double havd = sd*sd;
 619
 620         double aux = havd + cos (DegToRad(fromDEC)) * cos(DegToRad(toDEC)) * hava;
 621
 622         return (RadToDeg ( 2 * fabs(asin( sqrt(aux) ))));
 623 }
 624
 625 /* return current system time in message format */
 626 const char *
 627 timestamp()
 628 {
 629         static char ts[32];
 630         struct tm *tp;
 631         time_t t;
 632
 633         time (&t);
 634         tp = gmtime (&t);
 635         strftime (ts, sizeof(ts), "%Y-%m-%dT%H:%M:%S", tp);
 636         return (ts);
 637 }
 638