kstars/kstars/skypoint.h

   1 /***************************************************************************
   2                           skypoint.h  -  K Desktop Planetarium
   3                              -------------------
   4     begin                : Sun Feb 11 2001
   5     copyright            : (C) 2001-2005 by Jason Harris
   6     email                : jharris@30doradus.org
   7     copyright            : (C) 2004-2005 by Pablo de Vicente
   8     email                : p.devicente@wanadoo.es
   9  ***************************************************************************/
  10
  11 /***************************************************************************
  12  *                                                                         *
  13  *   This program is free software; you can redistribute it and/or modify  *
  14  *   it under the terms of the GNU General Public License as published by  *
  15  *   the Free Software Foundation; either version 2 of the License, or     *
  16  *   (at your option) any later version.                                   *
  17  *                                                                         *
  18  ***************************************************************************/
  19
  20 #ifndef SKYPOINT_H
  21 #define SKYPOINT_H
  22
  23 #include <qstring.h>
  24 #include <qptrlist.h>
  25
  26 #include "dms.h"
  27
  28 /**@class SkyPoint
  29         *
  30         *The sky coordinates of a point in the sky.  The
  31         *coordinates are stored in both Equatorial (Right Ascension,
  32         *Declination) and Horizontal (Azimuth, Altitude) coordinate systems.
  33         *Provides set/get functions for each coordinate angle, and functions
  34         *to convert between the Equatorial and Horizon coordinate systems.
  35         *
  36         *Because the coordinate values change slowly over time (due to
  37         *precession, nutation), the "catalog coordinates" are stored
  38         *(RA0, Dec0), which were the true coordinates on Jan 1, 2000.
  39         *The true coordinates (RA, Dec) at any other epoch can be found
  40         *from the catalog coordinates using updateCoords().
  41         *@short Stores dms coordinates for a point in the sky.
  42         *for converting between coordinate systems.
  43         *@author Jason Harris
  44         *@version 1.0
  45         */
  46
  47 class KSNumbers;
  48 class CSegment;
  49 class SkyObject;
  50
  51 class SkyPoint {
  52 public:
  53 /**Default constructor: Sets RA, Dec and RA0, Dec0 according
  54         *to arguments.  Does not set Altitude or Azimuth.
  55         *@param r Right Ascension
  56         *@param d Declination
  57         */
  58         SkyPoint( const dms& r, const dms& d ) { set( r, d ); }
  59
  60 /**Alternate constructor using pointer arguments, for convenience.
  61         *It behaves essentially like the default constructor.
  62         *@param r Right Ascension pointer
  63         *@param d Declination pointer
  64         */
  65         SkyPoint( const dms *r, const dms *d ) { set( dms(*r), dms(*d) ); }
  66
  67 /**Alternate constructor using double arguments, for convenience.
  68         *It behaves essentially like the default constructor.
  69         *@param r Right Ascension, expressed as a double
  70         *@param d Declination, expressed as a double
  71         */
  72         SkyPoint( double r=0.0, double d=0.0 ) { set( r, d ); }
  73
  74 /**
  75         *Empty destructor.
  76         */
  77         virtual ~SkyPoint();
  78
  79 ////
  80 //// 1.  Setting Coordinates
  81 //// =======================
  82
  83 /**Sets RA, Dec and RA0, Dec0 according to arguments.
  84         *Does not set Altitude or Azimuth.
  85         *@param r Right Ascension
  86         *@param d Declination
  87         */
  88         void set( const dms& r, const dms& d );
  89
  90 /**Overloaded member function, provided for convenience.
  91         *It behaves essentially like the above function.
  92         *@param r Right Ascension
  93         *@param d Declination
  94         */
  95         void set( const dms *r, const dms *d ) { set( *r, *d ); }
  96
  97 /**Overloaded member function, provided for convenience.
  98         *It behaves essentially like the above function.
  99         *@param r Right Ascension
 100         *@param d Declination
 101         */
 102         void set( double r, double d );
 103
 104 /**Sets RA0, the catalog Right Ascension.
 105         *@param r catalog Right Ascension.
 106         */
 107         void setRA0( dms r ) { RA0.set( r ); }
 108
 109 /**Overloaded member function, provided for convenience.
 110         *It behaves essentially like the above function.
 111         *@param r Right Ascension, expressed as a double.
 112         */
 113         void setRA0( double r ) { RA0.setH( r ); }
 114
 115 /**Sets Dec0, the catalog Declination.
 116         *@param d catalog Declination.
 117         */
 118         void setDec0( dms d ) { Dec0.set( d ); }
 119
 120 /**Overloaded member function, provided for convenience.
 121         *It behaves essentially like the above function.
 122         *@param d Declination, expressed as a double.
 123         */
 124         void setDec0( double d ) { Dec0.setD( d ); }
 125
 126 /**Sets RA, the current Right Ascension.
 127         *@param r Right Ascension.
 128         */
 129         void setRA( dms r ) { RA.set( r ); }
 130
 131 /**Overloaded member function, provided for convenience.
 132         *It behaves essentially like the above function.
 133         *@param r Right Ascension, expressed as a double.
 134         */
 135         void setRA( double r ) { RA.setH( r ); }
 136
 137 /**Sets Dec, the current Declination
 138         *@param d Declination.
 139         */
 140         void setDec( dms d ) { Dec.set( d ); }
 141
 142 /**Overloaded member function, provided for convenience.
 143         *It behaves essentially like the above function.
 144         *@param d Declination, expressed as a double.
 145         */
 146         void setDec( double d ) { Dec.setD( d ); }
 147
 148 /**Sets Alt, the Altitude.
 149         *@param alt Altitude.
 150         */
 151         void setAlt( dms alt ) { Alt.set( alt ); }
 152
 153 /**Overloaded member function, provided for convenience.
 154         *It behaves essentially like the above function.
 155         *@param alt Altitude, expressed as a double.
 156         */
 157         void setAlt( double alt ) { Alt.setD( alt ); }
 158
 159 /**Sets Az, the Azimuth.
 160         *@param az Azimuth.
 161         */
 162         void setAz( dms az ) { Az.set( az ); }
 163
 164 /**Overloaded member function, provided for convenience.
 165         *It behaves essentially like the above function.
 166         *@param az Azimuth, expressed as a double.
 167         */
 168         void setAz( double az ) { Az.setD( az ); }
 169
 170 /**Sets Galactic Longitude.
 171         *@param glo Galactic Longitude.
 172         */
 173 //      void setGalLong( dms glo ) { galLong.set( glo ); }
 174
 175 /**Overloaded member function, provided for convenience.
 176         *It behaves essentially like the above function.
 177         *@param glo Galactic Longitude, expressed as a double.
 178         */
 179 //      void setGalLong( double glo ) { galLong.setD( glo ); }
 180
 181 /**Sets Galactic Longitude.
 182         *@param gla Galactic Longitude.
 183         */
 184 //      void setGalLat( dms gla ) { galLat.set( gla ); }
 185
 186 /**Overloaded member function, provided for convenience.
 187         *It behaves essentially like the above function.
 188         *@param gla Galactic Longitude, expressed as a double.
 189         */
 190 //      void setGalLat( double gla ) { galLat.setD( gla ); }
 191
 192 ////
 193 //// 2. Returning coordinates.
 194 //// =========================
 195
 196 /**@return a pointer to the catalog Right Ascension.
 197         */
 198         const dms* ra0() const { return &RA0; };
 199
 200 /**@return a pointer to the catalog Declination.
 201         */
 202         const dms* dec0() const { return &Dec0; };
 203
 204 /**@returns a pointer to the current Right Ascension.
 205         */
 206         const dms* ra() const { return &RA; }
 207
 208 /**@return a pointer to the current Declination.
 209         */
 210         const dms* dec() const { return &Dec; }
 211
 212 /**@return a pointer to the current Azimuth.
 213         */
 214         const dms* az() const { return &Az; }
 215
 216 /**@return a pointer to the current Altitude.
 217         */
 218         const dms* alt() const { return &Alt; }
 219
 220 /**@return a pointer to the current galactic latitude.
 221         */
 222 //      const dms* gLat() const { return &galLat; }
 223
 224 /**@return a pointer to the current galactic longitude.
 225         */
 226 //      const dms* gLong() const { return &galLong; }
 227
 228 ////
 229 //// 3. Coordinate conversions.
 230 //// ==========================
 231
 232 /**Determine the (Altitude, Azimuth) coordinates of the
 233         *SkyPoint from its (RA, Dec) coordinates, given the local
 234         *sidereal time and the observer's latitude.
 235         *@param LST pointer to the local sidereal time
 236         *@param lat pointer to the geographic latitude
 237         */
 238         void EquatorialToHorizontal( const dms* LST, const dms* lat );
 239
 240 /**Determine the (RA, Dec) coordinates of the
 241         *SkyPoint from its (Altitude, Azimuth) coordinates, given the local
 242         *sidereal time and the observer's latitude.
 243         *@param LST pointer to the local sidereal time
 244         *@param lat pointer to the geographic latitude
 245         */
 246         void HorizontalToEquatorial( const dms* LST, const dms* lat );
 247
 248         /**@short Convert Right Ascension/Declination to Ecliptic logitude/latitude.
 249         *@param Obliquity current Obliquity of the Ecliptic (angle from Equator)
 250         */
 251         void EquatorialToEcliptic( const KSNumbers *num );
 252
 253         /**@short Convert Ecliptic logitude/latitude to Right Ascension/Declination.
 254         *@param Obliquity current Obliquity of the Ecliptic (angle from Equator)
 255         */
 256
 257         void EclipticToEquatorial( const KSNumbers *num );
 258
 259 /**Determine the Ecliptic coordinates of the SkyPoint, given the Julian Date.
 260         *The ecliptic coordinates are returned as reference arguments (since
 261         *they are not stored internally)
 262         */
 263         void findEcliptic( const dms *Obliquity, dms &EcLong, dms &EcLat );
 264
 265 /**Set the current (RA, Dec) coordinates of the
 266         *SkyPoint, given pointers to its Ecliptic (Long, Lat) coordinates, and
 267         *to the current obliquity angle (the angle between the equator and ecliptic).
 268         */
 269         void setFromEcliptic( const dms *Obliquity, const dms *EcLong, const dms *EcLat );
 270
 271 /** Computes galactic coordinates from equatorial coordinates referred to
 272         * epoch 1950. RA and Dec are, therefore assumed to be B1950
 273         * coordinates.
 274         */
 275         void Equatorial1950ToGalactic(dms &galLong, dms &galLat);
 276
 277 /** Computes equatorial coordinates referred to 1950 from galactic ones referred to
 278         * epoch B1950. RA and Dec are, therefore assumed to be B1950
 279         * coordinates.
 280         */
 281         void GalacticToEquatorial1950(const dms* galLong, const dms* galLat);
 282
 283 ////
 284 //// 4. Coordinate update/corrections.
 285 //// =================================
 286
 287 /**Determine the current coordinates (RA, Dec) from the catalog
 288         *coordinates (RA0, Dec0), accounting for both precession and nutation.
 289         *@param num pointer to KSNumbers object containing current values of
 290         *time-dependent variables.
 291         *@param includePlanets does nothing in this implementation (see KSPlanetBase::updateCoords()).
 292         *@param lat does nothing in this implementation (see KSPlanetBase::updateCoords()).
 293         *@param LST does nothing in this implementation (see KSPlanetBase::updateCoords()).
 294         */
 295         virtual void updateCoords( KSNumbers *num, bool includePlanets=true, const dms *lat=0, const dms *LST=0 );
 296
 297 /**Computes the apparent coordinates for this SkyPoint for any epoch,
 298         *accounting for the effects of precession, nutation, and aberration.
 299         *Similar to updateCoords(), but the starting epoch need not be
 300         *J2000, and the target epoch need not be the present time.
 301         *@param jd0 Julian Day which identifies the original epoch
 302         *@param jdf Julian Day which identifies the final epoch
 303         */
 304         void apparentCoord(long double jd0, long double jdf);
 305
 306 /**Determine the effects of nutation for this SkyPoint.
 307         *@param num pointer to KSNumbers object containing current values of
 308         *time-dependent variables.
 309         */
 310         void nutate(const KSNumbers *num);
 311
 312 /**Determine the effects of aberration for this SkyPoint.
 313         *@param num pointer to KSNumbers object containing current values of
 314         *time-dependent variables.
 315         */
 316         void aberrate(const KSNumbers *num);
 317
 318 /**General case of precession. It precess from an original epoch to a
 319         *final epoch. In this case RA0, and Dec0 from SkyPoint object represent
 320         *the coordinates for the original epoch and not for J2000, as usual.
 321         *@param jd0 Julian Day which identifies the original epoch
 322         *@param jdf Julian Day which identifies the final epoch
 323         */
 324         void precessFromAnyEpoch(long double jd0, long double jdf);
 325
 326         /** Determine the E-terms of aberration
 327          *In the past, the mean places of stars published in catalogs included
 328          *the contribution to the aberration due to the ellipticity of the orbit
 329          *of the Earth. These terms, known as E-terms were almost constant, and
 330          *in the newer catalogs (FK5) are not included. Therefore to convert from
 331          *FK4 to FK5 one has to compute these E-terms.
 332          */
 333         SkyPoint Eterms(void);
 334
 335         /** Exact precession from Besselian epoch 1950 to epoch J2000. The
 336         *coordinates referred to the first epoch are in the
 337         FK4 catalog, while the latter are in the Fk5 one.
 338         *Reference: Smith, C. A.; Kaplan, G. H.; Hughes, J. A.; Seidelmann,
 339         *P. K.; Yallop, B. D.; Hohenkerk, C. Y.
 340         *Astronomical Journal, vol. 97, Jan. 1989, p. 265-279
 341         *This transformation requires 4 steps:
 342         * - Correct E-terms
 343         * - Precess from B1950 to 1984, January 1st, 0h, using Newcomb expressions
 344         * - Add zero point correction in right ascension for 1984
 345         * - Precess from 1984, January 1st, 0h to J2000
 346         */
 347         void B1950ToJ2000(void);
 348
 349         /** Exact precession from epoch J2000 Besselian epoch 1950. The coordinates
 350         *referred to the first epoch are in the FK4 catalog, while the
 351         *latter are in the Fk5 one.
 352         *Reference: Smith, C. A.; Kaplan, G. H.; Hughes, J. A.; Seidelmann,
 353         *P. K.; Yallop, B. D.; Hohenkerk, C. Y.
 354         *Astronomical Journal, vol. 97, Jan. 1989, p. 265-279
 355         *This transformation requires 4 steps:
 356         * - Precess from J2000 to 1984, January 1st, 0h
 357         * - Add zero point correction in right ascension for 1984
 358         * - Precess from 1984, January 1st, 0h, to B1950 using Newcomb expressions
 359         * - Correct E-terms
 360         */
 361         void J2000ToB1950(void);
 362
 363         /** Coordinates in the FK4 catalog include the effect of aberration due
 364          *to the ellipticity of the orbit of the Earth. Coordinates in the FK5
 365          *catalog do not include these terms. In order to convert from B1950 (FK4)
 366          *to actual mean places one has to use this function.
 367         */
 368         void addEterms(void);
 369
 370         /** Coordinates in the FK4 catalog include the effect of aberration due
 371          *to the ellipticity of the orbit of the Earth. Coordinates in the FK5
 372          *catalog do not include these terms. In order to convert from
 373          * FK5 coordinates to B1950 (FK4) one has to use this function.
 374         */
 375         void subtractEterms(void);
 376
 377         /** Computes the angular distance between two SkyObjects. The algorithm
 378          *  to compute this distance is:
 379          *  cos(distance) = sin(d1)*sin(d2) + cos(d1)*cos(d2)*cos(a1-a2)
 380          *  where a1,d1 are the coordinates of the first object and a2,d2 are
 381          *  the coordinates of the second object.
 382          *  However this algorithm is not accurate when the angular separation
 383          *  is small.
 384          *  Meeus provides a different algorithm in page 111 which we
 385          *  implement here.
 386          *  @param SkyPoint to which distance is to be calculated
 387          *  @return dms angle representing angular separation.
 388          **/
 389
 390         dms angularDistanceTo( SkyPoint *sp);
 391
 392         bool operator == ( SkyPoint &p ) { return ( ra()->Degrees() == p.ra()->Degrees() && dec()->Degrees() == p.dec()->Degrees() ); }
 393
 394         /** Computes the velocity of the Sun projected on the direction of the source.
 395          *
 396          * @param jd Epoch expressed as julian day to which the source coordinates refer to.
 397          * @return Radial velocity of the source referred to the barycenter of the solar system in km/s
 398          **/
 399         double vRSun(long double jd);
 400
 401         /** Computes the radial velocity of a source referred to the solar system barycenter
 402          * from the radial velocity referred to the
 403          * Local Standard of Rest, aka known as VLSR. To compute it we need the coordinates of the
 404          * source the VLSR and the epoch for the source coordinates.
 405          *
 406          * @param vlsr radial velocity of the source referred to the LSR in km/s
 407          * @param jd Epoch expressed as julian day to which the source coordinates refer to.
 408          * @return Radial velocity of the source referred to the barycenter of the solar system in km/s
 409          **/
 410         double vHeliocentric(double vlsr, long double jd);
 411
 412         /** Computes the radial velocity of a source referred to the Local Standard of Rest, aka known as VLSR
 413          * from the radial velocity referred to the solar system barycenter
 414          *
 415          * @param vlsr radial velocity of the source referred to the LSR in km/s
 416          * @param jd Epoch expressed as julian day to which the source coordinates refer to.
 417          * @return Radial velocity of the source referred to the barycenter of the solar system in km/s
 418          **/
 419         double vHelioToVlsr(double vhelio, long double jd);
 420
 421         /** Computes the velocity of any object projected on the direction of the source.
 422          *  @param jd0 Julian day for which we compute the direction of the source
 423          *  @return velocity of the Earth projected on the direction of the source kms-1
 424          */
 425         double vREarth(long double jd0);
 426
 427         /** Computes the radial velocity of a source referred to the center of the earth
 428          * from the radial velocity referred to the solar system barycenter
 429          *
 430          * @param vhelio radial velocity of the source referred to the barycenter of the
 431          *               solar system in km/s
 432          * @param jd     Epoch expressed as julian day to which the source coordinates refer to.
 433          * @return Radial velocity of the source referred to the center of the Earth in km/s
 434          **/
 435         double vGeocentric(double vhelio, long double jd);
 436
 437         /** Computes the radial velocity of a source referred to the solar system barycenter
 438          * from the velocity referred to the center of the earth
 439          *
 440          * @param vgeo   radial velocity of the source referred to the center of the Earth
 441          *               [km/s]
 442          * @param jd     Epoch expressed as julian day to which the source coordinates refer to.
 443          * @return Radial velocity of the source referred to the solar system barycenter in km/s
 444          **/
 445         double vGeoToVHelio(double vgeo, long double jd);
 446
 447         /** Computes the velocity of any object (oberver's site) projected on the
 448          * direction of the source.
 449          *  @param vsite velocity of that object in cartesian coordinates
 450          *  @return velocity of the object projected on the direction of the source kms-1
 451          */
 452         double vRSite(double vsite[3]);
 453
 454         /** Computes the radial velocity of a source referred to the observer site on the surface
 455          * of the earth from the geocentric velovity and the velocity of the site referred to the center
 456          * of the Earth.
 457          *
 458          * @param vgeo radial velocity of the source referred to the center of the earth in km/s
 459          * @param vsite Velocity at which the observer moves referred to the center of the earth.
 460          * @return Radial velocity of the source referred to the observer's site in km/s
 461          **/
 462         double vTopocentric(double vgeo, double vsite[3]);
 463
 464         /** Computes the radial velocity of a source referred to the center of the Earth from
 465          * the radial velocity referred to an observer site on the surface of the earth
 466          *
 467          * @param vtopo radial velocity of the source referred to the observer's site in km/s
 468          * @param vsite Velocity at which the observer moves referred to the center of the earth.
 469          * @return Radial velocity of the source referred the center of the earth in km/s
 470          **/
 471         double vTopoToVGeo(double vtopo, double vsite[3]);
 472
 473 ////
 474 //// 5. Calculating Rise/Set/Transit data.
 475 //// =====================================
 476
 477 //// To Be Moved from SkyObject....
 478
 479 ////
 480 //// 6. Constellation Identification
 481 //// =====================================
 482
 483         QString constellation( QPtrList<CSegment> &seglist, QPtrList<SkyObject> &cnames ) const;
 484
 485
 486 protected:
 487 /**Precess this SkyPoint's catalog coordinates to the epoch described by the
 488         *given KSNumbers object.
 489         *@param num pointer to a KSNumbers object describing the target epoch.
 490         */
 491         void precess(const KSNumbers *num);
 492
 493
 494 private:
 495         dms RA0, Dec0; //catalog coordinates
 496         dms RA, Dec; //current true sky coordinates
 497         dms Alt, Az;
 498 };
 499
 500 #endif