00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "lunarphase.h"
00026 #include <config-kholidays.h>
00027
00028 #include <KDebug>
00029 #include <KGlobal>
00030 #include <KLocale>
00031
00032 #include <QtCore/QDateTime>
00033
00034 using namespace KHolidays;
00035
00036 static double percentFull( uint tmpt );
00037 static double degreesToRadians( double degree );
00038 static void adj360( double *degree );
00039
00040 QString LunarPhase::phaseNameAtDate( const QDate &date )
00041 {
00042 return phaseName( phaseAtDate( date ) );
00043 }
00044
00045 QString LunarPhase::phaseName( LunarPhase::Phase phase )
00046 {
00047 switch ( phase ) {
00048 case NewMoon:
00049 return( i18n( "New Moon" ) );
00050 case FullMoon:
00051 return( i18n( "Full Moon" ) );
00052 case FirstQuarter:
00053 return( i18n( "First Quarter Moon" ) );
00054 case LastQuarter:
00055 return( i18n( "Last Quarter Moon" ) );
00056 default:
00057 case None:
00058 return QString();
00059 }
00060 }
00061
00062 LunarPhase::Phase LunarPhase::phaseAtDate( const QDate &date )
00063 {
00064 Phase retPhase = None;
00065
00066
00067 const QTime anytime( 12, 0, 0 );
00068 const QDateTime today( date, anytime, Qt::UTC );
00069 const double todayPer = percentFull( today.toTime_t() ) + 0.5;
00070
00071 const QDateTime tomorrow( date.addDays( 1 ), anytime, Qt::UTC );
00072 const double tomorrowPer = percentFull( tomorrow.toTime_t() ) + 0.5;
00073
00074 if ( static_cast<int>( todayPer ) == 100 &&
00075 static_cast<int>( tomorrowPer ) != 100 ) {
00076 retPhase = FullMoon;
00077 } else if ( static_cast<int>( todayPer ) == 0 &&
00078 static_cast<int>( tomorrowPer ) != 0 ) {
00079 retPhase = NewMoon;
00080 } else {
00081 if ( todayPer > 50 && tomorrowPer < 50 ) {
00082 retPhase = LastQuarter;
00083 }
00084 if ( todayPer < 50 && tomorrowPer > 50 ) {
00085 retPhase = FirstQuarter;
00086 }
00087
00088
00089
00090 }
00091
00092 return( retPhase );
00093 }
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 #ifdef HAVE_SYS_CDEFS_H
00131 #include <sys/cdefs.h>
00132 #endif
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 #include <ctype.h>
00147 #ifdef HAVE_ERR_H
00148 #include <err.h>
00149 #endif
00150 #include <math.h>
00151 #include <string.h>
00152 #include <stdlib.h>
00153 #include <time.h>
00154 #include <unistd.h>
00155
00156 static double PI = 3.14159265358979323846;
00157
00158
00159
00160
00161
00162
00163
00164 static int EPOCH_MINUS_1970 = ( 20 * 365 + 5 - 1 );
00165 static double EPSILONg = 279.403303;
00166 static double RHOg = 282.768422;
00167 static double ECCEN = 0.016713;
00168 static double lzero = 318.351648;
00169 static double Pzero = 36.340410;
00170 static double Nzero = 318.510107;
00171
00172
00173
00174
00175
00176 static double percentFull( uint tmpt )
00177 {
00178 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
00179 double A4, lprime, V, ldprime, D, Nm;
00180
00181 double days;
00182 days = ( tmpt - EPOCH_MINUS_1970 * 86400 ) / 86400.0;
00183
00184 N = 360 * days / 365.242191;
00185 adj360( &N );
00186 Msol = N + EPSILONg - RHOg;
00187 adj360( &Msol );
00188 Ec = 360 / PI * ECCEN * sin( degreesToRadians( Msol ) );
00189 LambdaSol = N + Ec + EPSILONg;
00190 adj360( &LambdaSol );
00191 l = 13.1763966 * days + lzero;
00192 adj360( &l );
00193 Mm = l - ( 0.1114041 * days ) - Pzero;
00194 adj360( &Mm );
00195 Nm = Nzero - ( 0.0529539 * days );
00196 adj360( &Nm );
00197 Ev = 1.2739 * sin( degreesToRadians( 2 * ( l - LambdaSol ) - Mm ) );
00198 Ac = 0.1858 * sin( degreesToRadians( Msol ) );
00199 A3 = 0.37 * sin( degreesToRadians( Msol ) );
00200 Mmprime = Mm + Ev - Ac - A3;
00201 Ec = 6.2886 * sin( degreesToRadians( Mmprime ) );
00202 A4 = 0.214 * sin( degreesToRadians( 2 * Mmprime ) );
00203 lprime = l + Ev + Ec - Ac + A4;
00204 V = 0.6583 * sin( degreesToRadians( 2 * ( lprime - LambdaSol ) ) );
00205 ldprime = lprime + V;
00206 D = ldprime - LambdaSol;
00207 D = 50.0 * ( 1 - cos( degreesToRadians( D ) ) );
00208 return D;
00209 }
00210
00211
00212
00213
00214
00215 static double degreesToRadians( double degree )
00216 {
00217 return ( degree * PI ) / 180.00;
00218 }
00219
00220
00221
00222
00223
00224 static void adj360( double *degree )
00225 {
00226 for ( ;; ) {
00227 if ( *degree < 0 ) {
00228 *degree += 360;
00229 } else if ( *degree > 360 ) {
00230 *degree -= 360;
00231 } else {
00232 break;
00233 }
00234 }
00235 }