kstars/kstars/ksmoon.cpp

   1 /***************************************************************************
   2                           ksmoon.cpp  -  K Desktop Planetarium
   3                              -------------------
   4     begin                : Sun Aug 26 2001
   5     copyright            : (C) 2001 by Jason Harris
   6     email                : kstars@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 <stdlib.h>
  19 #include <math.h>
  20
  21 #include <qfile.h>
  22
  23 #include "ksnumbers.h"
  24 #include "ksutils.h"
  25 #include "kssun.h"
  26 #include "ksmoon.h"
  27 #include "kstarsdata.h"
  28
  29 KSMoon::KSMoon(KStarsData *kd)
  30  : KSPlanetBase( kd, I18N_NOOP( "Moon" ), "", 3474.8 /*diameter in km*/ ) {
  31         BData.setAutoDelete(true);
  32         LRData.setAutoDelete(true);
  33 }
  34
  35 KSMoon::~KSMoon() {
  36 }
  37
  38 bool KSMoon::data_loaded = false;
  39 QPtrList<KSMoon::MoonLRData> KSMoon::LRData;
  40 QPtrList<KSMoon::MoonBData> KSMoon::BData;
  41
  42 bool KSMoon::loadData() {
  43         if (data_loaded) return true;
  44
  45         QString line;
  46         QFile f;
  47         int nd, nm, nm1, nf;
  48         double Li, Ri, Bi; //coefficients of the sums
  49
  50         if ( KSUtils::openDataFile( f, "moonLR.dat" ) ) {
  51                 QTextStream stream( &f );
  52                 while ( !stream.eof() ) {
  53                         line = stream.readLine();
  54                         QTextIStream instream( &line );
  55                         instream >> nd >> nm >> nm1 >> nf >> Li >> Ri;
  56                         LRData.append(new MoonLRData(nd, nm, nm1, nf, Li, Ri));
  57                 }
  58                 f.close();
  59         } else
  60                 return false;
  61
  62
  63         if ( KSUtils::openDataFile( f, "moonB.dat" ) ) {
  64                 QTextStream stream( &f );
  65                 while ( !stream.eof() ) {
  66                         line = stream.readLine();
  67                         QTextIStream instream( &line );
  68                         instream >> nd >> nm >> nm1 >> nf >> Bi;
  69                         BData.append(new MoonBData(nd, nm, nm1, nf, Bi));
  70                 }
  71                 f.close();
  72         }
  73
  74         data_loaded = true;
  75         return true;
  76 }
  77
  78 bool KSMoon::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase* ) {
  79         //Algorithms in this subroutine are taken from Chapter 45 of "Astronomical Algorithms"
  80   //by Jean Meeus (1991, Willmann-Bell, Inc. ISBN 0-943396-35-2.  http://www.willbell.com/math/mc1.htm)
  81         QString fname, snum, line;
  82         QFile f;
  83         double DegtoRad;
  84         double T, L, D, M, M1, F, E, A1, A2, A3;
  85         double sumL, sumR, sumB;
  86
  87         DegtoRad = acos( -1.0 )/180.0;
  88
  89         //Julian centuries since J2000
  90         T = num->julianCenturies();
  91
  92         double Et = 1.0 - 0.002516*T - 0.0000074*T*T;
  93
  94         //Moon's mean longitude
  95         L = 218.3164591 + 481267.88134236*T - 0.0013268*T*T + T*T*T/538841.0 - T*T*T*T/65194000.0;
  96         while ( L > 360.0 ) L -= 360.0;
  97         while ( L < 0.0 ) L += 360.0;
  98         L *= DegtoRad;
  99         //Moon's mean elongation
 100         D = 297.8502042 + 445267.1115168*T - 0.0016300*T*T + T*T*T/545868.0 - T*T*T*T/113065000.0;
 101         while ( D > 360.0 ) D -= 360.0;
 102         while ( D < 0.0 ) D += 360.0;
 103         D *= DegtoRad;
 104         //Sun's mean anomaly
 105   M = 357.5291092 + 35999.0502909*T - 0.0001536*T*T + T*T*T/24490000.0;
 106         while ( M > 360.0 ) M -= 360.0;
 107         while ( M < 0.0 ) M += 360.0;
 108         M *= DegtoRad;
 109         //Moon's mean anomaly
 110         M1= 134.9634114 + 477198.8676313*T + 0.0089970*T*T + T*T*T/69699.0 - T*T*T*T/14712000.0;
 111         while ( M1 > 360.0 ) M1 -= 360.0;
 112         while ( M1 < 0.0 ) M1 += 360.0;
 113         M1 *= DegtoRad;
 114   //Moon's argument of latitude (angle from ascending node)
 115         F = 93.2720993 + 483202.0175273*T - 0.0034029*T*T - T*T*T/3526000.0 + T*T*T*T/863310000.0;
 116         while ( F > 360.0 ) F -= 360.0;
 117         while ( F < 0.0 ) F += 360.0;
 118         F *= DegtoRad;
 119
 120         A1 = 119.75 +    131.849*T;
 121         A2 =  53.09 + 479264.290*T;
 122         A3 = 313.45 + 481226.484*T;
 123         while ( A1 > 360.0 ) A1 -= 360.0;
 124         while ( A1 < 0.0 ) A1 += 360.0;
 125         while ( A2 > 360.0 ) A2 -= 360.0;
 126         while ( A2 < 0.0 ) A2 += 360.0;
 127         while ( A3 > 360.0 ) A3 -= 360.0;
 128         while ( A3 < 0.0 ) A3 += 360.0;
 129         A1 *= DegtoRad;
 130         A2 *= DegtoRad;
 131         A3 *= DegtoRad;
 132
 133         //Calculate the series expansions stored in moonLR.txt and moonB.txt.
 134         //
 135         sumL = 0.0;
 136         sumR = 0.0;
 137
 138         if (!loadData()) return false;
 139
 140         for (MoonLRData *mlrd = LRData.first(); mlrd != 0; mlrd = LRData.next()) {
 141
 142                 E = 1.0;
 143                 if ( mlrd->nm ) { //if M != 0, include changing eccentricity of Earth's orbit
 144                         E = Et;
 145                         if ( abs( mlrd->nm )==2 ) E = E*E; //use E^2
 146                 }
 147                 sumL += E*mlrd->Li*sin( mlrd->nd*D + mlrd->nm*M + mlrd->nm1*M1 + mlrd->nf*F );
 148                 sumR += E*mlrd->Ri*cos( mlrd->nd*D + mlrd->nm*M + mlrd->nm1*M1 + mlrd->nf*F );
 149         }
 150
 151         sumB = 0.0;
 152         for (MoonBData *mbd = BData.first(); mbd != 0; mbd = BData.next()) {
 153
 154                 E = 1.0;
 155                 if ( mbd->nm ) { //if M != 0, include changing eccentricity of Earth's orbit
 156                         E = Et;
 157                         if ( abs( mbd->nm )==2 ) E = E*E; //use E^2
 158                 }
 159                 sumB += E*mbd->Bi*sin( mbd->nd*D + mbd->nm*M + mbd->nm1*M1 + mbd->nf*F );
 160         }
 161
 162         //Additive terms for sumL and sumB
 163         sumL += ( 3958.0*sin( A1 ) + 1962.0*sin( L-F ) + 318.0*sin( A2 ) );
 164         sumB += ( -2235.0*sin( L ) + 382.0*sin( A3 ) + 175.0*sin( A1-F ) + 175.0*sin( A1+F ) + 127.0*sin( L-M1 ) - 115.0*sin( L+M1 ) );
 165
 166         //Geocentric coordinates
 167         setEcLong( ( L + DegtoRad*sumL/1000000.0 ) * 180./dms::PI );  //convert radians to degrees
 168         setEcLat( ( DegtoRad*sumB/1000000.0 ) * 180./dms::PI );
 169         Rearth = ( 385000.56 + sumR/1000.0 )/AU_KM; //distance from Earth, in AU
 170         setAngularSize( asin(physicalSize()/Rearth/AU_KM)*60.*180./dms::PI ); //angular size in arcmin
 171
 172         EclipticToEquatorial( num->obliquity() );
 173
 174         //Determine position angle
 175         findPA( num );
 176
 177         return true;
 178 }
 179
 180 void KSMoon::findPhase( const KSSun *Sun ) {
 181         Phase.setD( ecLong()->Degrees() - Sun->ecLong()->Degrees() );
 182         Phase.setD( Phase.reduce().Degrees() );
 183         int iPhase = int( 0.1*Phase.Degrees()+0.5 );
 184         if (iPhase==36) iPhase = 0;
 185         QString sPhase;
 186         sPhase = sPhase.sprintf( "%02d", iPhase );
 187         QString imName = "moon" + sPhase + ".png";
 188
 189         QFile imFile;
 190         if ( KSUtils::openDataFile( imFile, imName ) ) {
 191                 imFile.close();
 192                 image0()->load( imFile.name() );
 193                 image()->load( imFile.name() );
 194
 195         }
 196 }
 197
 198 QString KSMoon::phaseName() const {
 199         double f = illum();
 200         double p = phase().Degrees();
 201
 202         //First, handle the major phases
 203         if ( f > 0.99 ) return i18n( "moon phase, 100 percent illuminated", "Full moon" );
 204         if ( f < 0.01 ) return i18n( "moon phase, 0 percent illuminated", "New moon" );
 205         if ( fabs( f - 0.50 ) < 0.01 ) {
 206                 if ( p < 180.0 ) return i18n( "moon phase, half-illuminated and growing", "First quarter" );
 207                 else return i18n( "moon phase, half-illuminated and shrinking", "Third quarter" );
 208         }
 209
 210         //Next, handle the more general cases
 211         if ( p < 90.0 ) return i18n( "moon phase between new moon and 1st quarter", "Waxing crescent" );
 212         else if ( p < 180.0 ) return i18n( "moon phase between 1st quarter and full moon", "Waxing gibbous" );
 213         else if ( p < 270.0 ) return i18n( "moon phase between full moon and 3rd quarter", "Waning gibbous" );
 214         else if ( p < 360.0 ) return i18n( "moon phase between 3rd quarter and new moon", "Waning crescent" );
 215
 216         else return i18n( "unknown" );
 217 }