00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <kdebug.h>
00019
00020 #include "kstarsdata.h"
00021 #include "kstarsdatetime.h"
00022 #include "ksnumbers.h"
00023 #include "dms.h"
00024 #include "kscomet.h"
00025
00026
00027 KSComet::KSComet( KStarsData *_kd, QString _s, QString imfile,
00028 long double _JD, double _q, double _e, dms _i, dms _w, dms _Node, double Tp )
00029 : KSPlanetBase(_kd, _s, imfile), kd(_kd), JD(_JD), q(_q), e(_e), i(_i), w(_w), N(_Node) {
00030
00031 setType( 9 );
00032
00033
00034
00035 int year = int( Tp/10000.0 );
00036 int month = int( (int(Tp) % 10000)/100.0 );
00037 int day = int( int(Tp) % 100 );
00038 double Hour = 24.0 * ( Tp - int(Tp) );
00039 int h = int( Hour );
00040 int m = int( 60.0 * ( Hour - h ) );
00041 int s = int( 60.0 * ( 60.0 * ( Hour - h) - m ) );
00042
00043 JDp = KStarsDateTime( ExtDate( year, month, day ), QTime( h, m, s ) ).djd();
00044
00045
00046 a = q/(1.0-e);
00047
00048
00049 P = 365.2568984 * pow(a, 1.5);
00050
00051
00052 if ( name().contains( "/" ) ) {
00053 setLongName( name() );
00054 setName( name().mid( name().find("/") + 1 ) );
00055 }
00056 }
00057
00058 bool KSComet::findGeocentricPosition( const KSNumbers *num, const KSPlanetBase *Earth ) {
00059 double v(0.0), r(0.0);
00060
00061
00062 dms n = dms( double(N.Degrees() - 3.82394E-5 * ( num->julianDay() - J2000 )) ).reduce();
00063
00064 if ( e > 0.98 ) {
00065
00066 double k = 0.01720209895;
00067 double a = 0.75 * ( num->julianDay() - JDp ) * k * sqrt( (1+e)/(q*q*q) );
00068 double b = sqrt( 1.0 + a*a );
00069 double W = pow((b+a),1.0/3.0) - pow((b-a),1.0/3.0);
00070 double c = 1.0 + 1.0/(W*W);
00071 double f = (1.0-e)/(1.0+e);
00072 double g = f/(c*c);
00073
00074 double a1 = (2.0/3.0) + (2.0*W*W/5.0);
00075 double a2 = (7.0/5.0) + (33.0*W*W/35.0) + (37.0*W*W*W*W/175.0);
00076 double a3 = W*W*( (432.0/175.0) + (956.0*W*W/1125.0) + (84.0*W*W*W*W/1575.0) );
00077 double w = W*(1.0 + g*c*( a1 + a2*g + a3*g*g ));
00078
00079 v = 2.0*atan(w) / dms::DegToRad;
00080 r = q*( 1.0 + w*w )/( 1.0 + w*w*f );
00081 } else {
00082
00083
00084 dms m = dms( double(360.0*( num->julianDay() - JDp )/P) ).reduce();
00085 double sinm, cosm;
00086 m.SinCos( sinm, cosm );
00087
00088
00089 double E = m.Degrees() + e*180.0/dms::PI * sinm * ( 1.0 + e*cosm );
00090
00091 if ( e > 0.05 ) {
00092 double E0;
00093 int iter(0);
00094 do {
00095 E0 = E;
00096 iter++;
00097 E = E0 - ( E0 - e*180.0/dms::PI *sin( E0*dms::DegToRad ) - m.Degrees() )/(1 - e*cos( E0*dms::DegToRad ) );
00098 } while ( fabs( E - E0 ) > 0.001 && iter < 1000 );
00099 }
00100
00101 double sinE, cosE;
00102 dms E1( E );
00103 E1.SinCos( sinE, cosE );
00104
00105 double xv = a * ( cosE - e );
00106 double yv = a * sqrt( 1.0 - e*e ) * sinE;
00107
00108
00109
00110 v = atan( yv/xv ) / dms::DegToRad;
00111
00112 if ( xv < 0.0 ) v += 180.0;
00113
00114 r = sqrt( xv*xv + yv*yv );
00115 }
00116
00117
00118 dms vw( v + w.Degrees() );
00119 double sinN, cosN, sinvw, cosvw, sini, cosi;
00120
00121 n.SinCos( sinN, cosN );
00122 vw.SinCos( sinvw, cosvw );
00123 i.SinCos( sini, cosi );
00124
00125
00126 double xh = r * ( cosN * cosvw - sinN * sinvw * cosi );
00127 double yh = r * ( sinN * cosvw + cosN * sinvw * cosi );
00128 double zh = r * ( sinvw * sini );
00129
00130
00131 double cosBe, sinBe, cosLe, sinLe;
00132 Earth->ecLong()->SinCos( sinLe, cosLe );
00133 Earth->ecLat()->SinCos( sinBe, cosBe );
00134
00135 double xe = Earth->rsun() * cosBe * cosLe;
00136 double ye = Earth->rsun() * cosBe * sinLe;
00137 double ze = Earth->rsun() * sinBe;
00138
00139
00140 xh -= xe;
00141 yh -= ye;
00142 zh -= ze;
00143
00144
00145 double ELongRad = atan( yh/xh );
00146
00147 if ( xh < 0.0 ) ELongRad += dms::PI;
00148
00149 double rr = sqrt( xh*xh + yh*yh );
00150 double ELatRad = atan( zh/rr );
00151
00152 ep.longitude.setRadians( ELongRad );
00153 ep.latitude.setRadians( ELatRad );
00154 setRsun( r );
00155 setRearth( Earth );
00156
00157 EclipticToEquatorial( num->obliquity() );
00158 nutate( num );
00159 aberrate( num );
00160
00161 return true;
00162 }
00163
00164
00165 bool KSComet::loadData() { return false; }