kstars/kstars/kssun.cpp

   1 /***************************************************************************
   2                           kssun.cpp  -  K Desktop Planetarium
   3                              -------------------
   4     begin                : Sun Jul 22 2001
   5     copyright            : (C) 2001 by Jason Harris
   6     email                : jharris@30doradus.org
   7  ***************************************************************************/
   8
   9 /***************************************************************************
  10  *                                                                         *
  11  *   This program is free software; you can redistribute it and/or modify  *
  12  *   it under the terms of the GNU General Public License as published by  *
  13  *   the Free Software Foundation; either version 2 of the License, or     *
  14  *   (at your option) any later version.                                   *
  15  *                                                                         *
  16  ***************************************************************************/
  17
  18 #include <math.h>
  19 #include <qdatetime.h>
  20
  21 #include "kssun.h"
  22 #include "ksutils.h"
  23 #include "ksnumbers.h"
  24 #include "kstarsdatetime.h"
  25
  26 KSSun::KSSun( KStarsData *kd, QString fn, double pSize ) : KSPlanet( kd, I18N_NOOP( "Sun" ), fn, pSize ) {
  27         /*
  28         JD0 = 2447892.5; //Jan 1, 1990
  29         eclong0 = 279.403303; //mean ecliptic longitude at JD0
  30         plong0 = 282.768422; //longitude of sun at perigee for JD0
  31         e0 = 0.016713; //eccentricity of Earth's orbit at JD0
  32         */
  33 }
  34
  35 bool KSSun::loadData() {
  36 //      kdDebug() << k_funcinfo << endl;
  37         return (odm.loadData("earth") != 0);
  38 }
  39
  40 bool KSSun::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase *Earth ) {
  41         if (Earth) {
  42                 //
  43                 // For the precision we need, the earth's orbit is circular.
  44                 // So don't bother to iterate like KSPlanet does. Just subtract
  45                 // The current delay and recompute (once).
  46                 //
  47                 double delay = (.0057755183 * Earth->rsun()) / 365250.0;
  48
  49                 //
  50                 // MHH 2002-02-04 I don't like this. But it avoids code duplication.
  51                 // Maybe we can find a better way.
  52                 //
  53                 const KSPlanet *pEarth = dynamic_cast<const KSPlanet *>(Earth);
  54                 /* FIXME: if you use pEarth at some point again, make sure you
  55                         check for 0L after the dynamic_cast! */
  56                 EclipticPosition trialpos;
  57                 pEarth->calcEcliptic(num->julianMillenia() - delay, trialpos);
  58
  59                 setEcLong( trialpos.longitude.Degrees() + 180.0 );
  60                 setEcLong( ecLong()->reduce().Degrees() );
  61                 setEcLat( -1.0*trialpos.latitude.Degrees() );
  62
  63         } else {
  64                 double sum[6];
  65                 dms EarthLong, EarthLat; //heliocentric coords of Earth
  66                 OrbitDataColl * odc;
  67                 double T = num->julianMillenia(); //Julian millenia since J2000
  68                 double Tpow[6];
  69
  70                 Tpow[0] = 1.0;
  71                 for (int i=1; i<6; ++i) {
  72                         Tpow[i] = Tpow[i-1] * T;
  73                 }
  74                         //First, find heliocentric coordinates
  75
  76                 if (!(odc =  odm.loadData("earth"))) return false;
  77
  78                 //Ecliptic Longitude
  79                 for (int i=0; i<6; ++i) {
  80                         sum[i] = 0.0;
  81                         for (uint j = 0; j < odc->Lon[i].size(); ++j) {
  82                                 sum[i] += odc->Lon[i][j]->A * cos( odc->Lon[i][j]->B + odc->Lon[i][j]->C*T );
  83                         }
  84                         sum[i] *= Tpow[i];
  85                         //kdDebug() << name() << " : sum[" << i << "] = " << sum[i] <<endl;
  86                 }
  87
  88                 EarthLong.setRadians( sum[0] + sum[1] + sum[2] +
  89                                 sum[3] + sum[4] + sum[5] );
  90                 EarthLong.setD( EarthLong.reduce().Degrees() );
  91
  92                 //Compute Ecliptic Latitude
  93                 for (int i=0; i<6; ++i) {
  94                         sum[i] = 0.0;
  95                         for (uint j = 0; j < odc->Lat[i].size(); ++j) {
  96                                 sum[i] += odc->Lat[i][j]->A * cos( odc->Lat[i][j]->B + odc->Lat[i][j]->C*T );
  97                         }
  98                         sum[i] *= Tpow[i];
  99                 }
 100
 101
 102                 EarthLat.setRadians( sum[0] + sum[1] + sum[2] + sum[3] +
 103                                 sum[4] + sum[5] );
 104
 105                 //Compute Heliocentric Distance
 106                 for (int i=0; i<6; ++i) {
 107                         sum[i] = 0.0;
 108                         for (uint j = 0; j < odc->Dst[i].size(); ++j) {
 109                                 sum[i] += odc->Dst[i][j]->A * cos( odc->Dst[i][j]->B + odc->Dst[i][j]->C*T );
 110                         }
 111                         sum[i] *= Tpow[i];
 112                 }
 113
 114                 ep.radius = sum[0] + sum[1] + sum[2] + sum[3] + sum[4] + sum[5];
 115
 116                 setEcLong( EarthLong.Degrees() + 180.0 );
 117                 setEcLong( ecLong()->reduce().Degrees() );
 118                 setEcLat( -1.0*EarthLat.Degrees() );
 119         }
 120
 121         //Finally, convert Ecliptic coords to Ra, Dec.  Ecliptic latitude is zero, by definition
 122         EclipticToEquatorial( num->obliquity() );
 123
 124         nutate(num);
 125         aberrate(num);
 126
 127         // We obtain the apparent geocentric ecliptic coordinates. That is, after
 128         // nutation and aberration have been applied.
 129         EquatorialToEcliptic( num->obliquity() );
 130
 131         //Determine the position angle
 132         findPA( num );
 133
 134         return true;
 135 }
 136
 137 long double KSSun::springEquinox(int year) {
 138         return equinox(year, 18, 3, 0.);
 139 }
 140
 141 long double KSSun::summerSolstice(int year) {
 142         return equinox(year, 18, 6, 90.);
 143 }
 144
 145 long double KSSun::autumnEquinox(int year) {
 146         return equinox(year, 19, 9, 180.);
 147 }
 148
 149 long double KSSun::winterSolstice(int year) {
 150         return equinox(year, 18, 12, 270.);
 151 }
 152
 153 long double KSSun::equinox(int year, int d, int m, double angle) {
 154         long double jd0[5];
 155         long double eclipticLongitude[5];
 156
 157         for(int i = 0; i<5; ++i) {
 158                 jd0[i] = KStarsDateTime( ExtDate(year,m,d+i), QTime(0,0,0) ).djd();
 159                 KSNumbers *ksn = new KSNumbers(jd0[i]);
 160                 //FIXME this is the Earth position
 161                 findGeocentricPosition( ksn );
 162                 delete ksn;
 163                 eclipticLongitude[i] = (long double)ecLong()->Degrees();
 164         }
 165
 166         return KSUtils::lagrangeInterpolation( eclipticLongitude, jd0, 5, angle );
 167 }