kstars/kstars/ksplanetbase.cpp

   1 /***************************************************************************
   2                           ksplanetbase.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
  20 #include <qfile.h>
  21 #include <qpoint.h>
  22 #include <qwmatrix.h>
  23
  24 #include "ksplanetbase.h"
  25 #include "ksplanet.h"
  26 #include "kstarsdata.h"
  27 #include "ksutils.h"
  28 #include "ksnumbers.h"
  29 #include "kspopupmenu.h"
  30
  31
  32 KSPlanetBase::KSPlanetBase( KStarsData *kd, QString s, QString image_file, double pSize )
  33  : SkyObject( 2, 0.0, 0.0, 0.0, s, "" ), Rearth(0.0), Image(0), data(kd) {
  34
  35          if (! image_file.isEmpty()) {
  36                 QFile imFile;
  37
  38                 if ( KSUtils::openDataFile( imFile, image_file ) ) {
  39                         imFile.close();
  40                         Image0.load( imFile.name() );
  41                         Image = Image0.convertDepth( 32 );
  42                         Image0 = Image;
  43                 }
  44         }
  45         PositionAngle = 0.0;
  46         ImageAngle = 0.0;
  47         PhysicalSize = pSize;
  48         Trail.setAutoDelete( TRUE );
  49 }
  50
  51 void KSPlanetBase::EquatorialToEcliptic( const dms *Obliquity ) {
  52         findEcliptic( Obliquity, ep.longitude, ep.latitude );
  53 }
  54
  55 void KSPlanetBase::EclipticToEquatorial( const dms *Obliquity ) {
  56         setFromEcliptic( Obliquity, &ep.longitude, &ep.latitude );
  57 }
  58
  59 void KSPlanetBase::updateCoords( KSNumbers *num, bool includePlanets, const dms *lat, const dms *LST ){
  60         if ( includePlanets ) {
  61                 data->earth()->findPosition( num ); //since we don't pass lat & LST, localizeCoords will be skipped
  62
  63                 if ( lat && LST ) {
  64                         findPosition( num, lat, LST, data->earth() );
  65                         if ( hasTrail() ) Trail.removeLast();
  66                 } else {
  67                         findGeocentricPosition( num, data->earth() );
  68                 }
  69         }
  70 }
  71
  72 void KSPlanetBase::findPosition( const KSNumbers *num, const dms *lat, const dms *LST, const KSPlanetBase *Earth ) {
  73         findGeocentricPosition( num, Earth );  //private function, reimplemented in each subclass
  74
  75         if ( Earth ) setRearth( Earth );
  76
  77         if ( lat && LST )
  78                 localizeCoords( num, lat, LST ); //correct for figure-of-the-Earth
  79
  80         if ( hasTrail() ) {
  81                 Trail.append( new SkyPoint( ra(), dec() ) );
  82                 if ( Trail.count() > MAXTRAIL ) Trail.removeFirst();
  83         }
  84
  85         if ( isMajorPlanet() )
  86                 findMagnitude(num);
  87 }
  88
  89 bool KSPlanetBase::isMajorPlanet() const {
  90         if ( name() == "Mercury" || name() == "Venus" || name() == "Mars" ||
  91                                 name() == "Jupiter" || name() == "Saturn" || name() == "Uranus" ||
  92                                 name() == "Neptune" || name() == "Pluto" )
  93                 return true;
  94
  95         return false;
  96 }
  97
  98 void KSPlanetBase::localizeCoords( const KSNumbers *num, const dms *lat, const dms *LST ) {
  99         //convert geocentric coordinates to local apparent coordinates (topocentric coordinates)
 100         dms HA, HA2; //Hour Angle, before and after correction
 101         double rsinp, rcosp, u, sinHA, cosHA, sinDec, cosDec, D;
 102         double cosHA2;
 103         double r = Rearth * AU_KM; //distance from Earth, in km
 104         u = atan( 0.996647*tan( lat->radians() ) );
 105         rsinp = 0.996647*sin( u );
 106         rcosp = cos( u );
 107         HA.setD( LST->Degrees() - ra()->Degrees() );
 108         HA.SinCos( sinHA, cosHA );
 109         dec()->SinCos( sinDec, cosDec );
 110
 111         D = atan( ( rcosp*sinHA )/( r*cosDec/6378.14 - rcosp*cosHA ) );
 112         dms temp;
 113         temp.setRadians( ra()->radians() - D );
 114         setRA( temp );
 115
 116         HA2.setD( LST->Degrees() - ra()->Degrees() );
 117         cosHA2 = cos( HA2.radians() );
 118         temp.setRadians( atan( cosHA2*( r*sinDec/6378.14 - rsinp )/( r*cosDec*cosHA/6378.14 - rcosp ) ) );
 119         setDec( temp );
 120
 121         EquatorialToEcliptic( num->obliquity() );
 122 }
 123
 124 void KSPlanetBase::setRearth( const KSPlanetBase *Earth ) {
 125         double sinL, sinB, sinL0, sinB0;
 126         double cosL, cosB, cosL0, cosB0;
 127         double x,y,z;
 128
 129         //The Moon's Rearth is set in its findGeocentricPosition()...
 130         if ( name() == "Moon" ) {
 131                 return;
 132         }
 133
 134         if ( name() == "Earth" ) {
 135                 Rearth = 0.0;
 136                 return;
 137         }
 138
 139         if ( ! Earth && name() != "Moon" ) {
 140                 kdDebug() << i18n( "KSPlanetBase::setRearth():  Error: Need an Earth pointer.  (" ) << name() << ")" << endl;
 141                 Rearth = 1.0;
 142                 return;
 143         }
 144
 145         Earth->ecLong()->SinCos( sinL0, cosL0 );
 146         Earth->ecLat()->SinCos( sinB0, cosB0 );
 147         double eX = Earth->rsun()*cosB0*cosL0;
 148         double eY = Earth->rsun()*cosB0*sinL0;
 149         double eZ = Earth->rsun()*sinB0;
 150
 151         ecLong()->SinCos( sinL, cosL );
 152         ecLat()->SinCos( sinB, cosB );
 153         x = rsun()*cosB*cosL - eX;
 154         y = rsun()*cosB*sinL - eY;
 155         z = rsun()*sinB - eZ;
 156
 157         Rearth = sqrt(x*x + y*y + z*z);
 158
 159         //Set angular size, in arcmin
 160         AngularSize = asin(PhysicalSize/Rearth/AU_KM)*60.*180./dms::PI;
 161 }
 162
 163 void KSPlanetBase::updateTrail( dms *LST, const dms *lat ) {
 164         for ( SkyPoint *sp = Trail.first(); sp; sp = Trail.next() )
 165                 sp->EquatorialToHorizontal( LST, lat );
 166 }
 167
 168 void KSPlanetBase::findPA( const KSNumbers *num ) {
 169         //Determine position angle of planet (assuming that it is aligned with
 170         //the Ecliptic, which is only roughly correct).
 171         //Displace a point along +Ecliptic Latitude by 1 degree
 172         SkyPoint test;
 173         dms newELat( ecLat()->Degrees() + 1.0 );
 174         test.setFromEcliptic( num->obliquity(), ecLong(), &newELat );
 175         double dx = test.ra()->Degrees() - ra()->Degrees();
 176         double dy = dec()->Degrees() - test.dec()->Degrees();
 177         double pa;
 178         if ( dy ) {
 179                 pa = atan( dx/dy )*180.0/dms::PI;
 180         } else {
 181                 pa = 90.0;
 182                 if ( dx > 0 ) pa = -90.0;
 183         }
 184         setPA( pa );
 185 }
 186
 187 void KSPlanetBase::rotateImage( double imAngle ) {
 188         ImageAngle = imAngle;
 189         QWMatrix m;
 190         m.rotate( ImageAngle );
 191         Image = Image0.xForm( m );
 192 }
 193
 194 void KSPlanetBase::scaleRotateImage( int scale, double imAngle ) {
 195         ImageAngle = imAngle;
 196         QWMatrix m;
 197         m.rotate( ImageAngle );
 198         Image = Image0.xForm( m ).smoothScale( scale, scale );
 199 }
 200
 201 void KSPlanetBase::findMagnitude(const KSNumbers *num) {
 202
 203         double cosDec, sinDec;
 204         dec()->SinCos(cosDec, sinDec);
 205
 206         // I would rather like to obtain the Earth Sun distance rather than assigning it
 207         // by hand..... Still TODO
 208         // double earthSun = Earth->rsun();
 209         double earthSun = 1.;
 210
 211         /* Phase of the planet in degrees */
 212
 213         double cosPhase = (rsun()*rsun() + rearth()*rearth() - earthSun*earthSun) / (2 * rsun() * rearth() );
 214         double phase = acos ( cosPhase ) * dms::PI / 180.;
 215
 216        /* Computation of the visual magnitude (V band) of the planet.
 217         * Algorithm provided by Pere Planesas (Observatorio Astronomico Nacional)
 218         * It has some simmilarity to J. Meeus algorithm in Astronomical Algorithms, Chapter 40.
 219         * */
 220
 221         // Initialized to the faintest magnitude observable with the HST
 222         float magnitude = 30;
 223
 224         double param = 5 * log10(rsun() * rearth() );
 225         double f1 = phase/100.;
 226         if ( name() == "Mercury" ) {
 227                 if (phase > 150)
 228                         f1 = 1.5;
 229                 magnitude = -0.36 + param + 3.8*f1 - 2.73*f1*f1 + 2*f1*f1*f1;
 230         }
 231         if ( name() =="Venus")
 232                 magnitude = -4.29 + param + 0.09*f1 + 2.39*f1*f1 - 0.65*f1*f1*f1;
 233         if( name() == "Mars")
 234                 magnitude = -1.52 + param + 0.016*phase;
 235         if( name() == "Jupiter")
 236                 magnitude = -9.25 + param + 0.005*phase;
 237
 238         if( name() == "Saturn") {
 239                 double T = num->julianCenturies();
 240                 double a0 = (40.66-4.695*T)* dms::PI / 180.;
 241                 double d0 = (83.52+0.403*T)* dms::PI / 180.;
 242                 double sinx = -cos(d0)*cosDec*cos(a0 - ra()->radians());
 243                 sinx = fabs(sinx-sin(d0)*sinDec);
 244                 double rings = -2.6*sinx + 1.25*sinx*sinx;
 245                 magnitude = -8.88 + param + 0.044*phase + rings;
 246         }
 247
 248         if( name() == "Uranus")
 249                 magnitude = -7.19 + param + 0.0028*phase;
 250         if( name() == "Neptune")
 251                 magnitude = -6.87 + param;
 252         if( name() == "Pluto")
 253                 magnitude = -1.01 + param + 0.041*phase;
 254
 255         setMag(magnitude);
 256 }