KHolidays Library
lunarphase.cpp
00001 /* 00002 This file is part of the kholidays library. 00003 00004 Copyright (c) 2004,2007 Allen Winter <winter@kde.org> 00005 00006 Copyright (c) 1989, 1993 //krazy:exclude=copyright 00007 The Regents of the University of California. All rights reserved. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Library General Public 00011 License as published by the Free Software Foundation; either 00012 version 2 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00017 GNU Library General Public License for more details. 00018 00019 You should have received a copy of the GNU Library General Public License 00020 along with this library; see the file COPYING.LIB. If not, write to the 00021 Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 00022 Boston, MA 02110-1301, USA. 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 // compute percent-full for the middle of today and yesterday. 00067 const QTime anytime( 0, 0, 0 ); 00068 const QDateTime today( date, anytime, Qt::UTC ); 00069 const double todayPer = percentFull( today.toTime_t() ) + 0.5; 00070 00071 if ( static_cast<int>( todayPer ) == 100 ) { 00072 retPhase = FullMoon; 00073 //kDebug() << date << " todayPer: " << todayPer << " " << phaseName( retPhase ); 00074 } else if ( static_cast<int>( todayPer ) == 0 ) { 00075 retPhase = NewMoon; 00076 //kDebug() << date << " todayPer: " << todayPer << " " << phaseName( retPhase ); 00077 } else { 00078 const QDateTime tomorrow( date.addDays( 1 ), anytime, Qt::UTC ); 00079 const double tomorrowPer = percentFull( tomorrow.toTime_t() ); 00080 00081 if ( static_cast<int>( todayPer ) == 50 ) { 00082 retPhase = ( tomorrowPer > todayPer ) ? FirstQuarter : LastQuarter; 00083 //kDebug() << date 00084 // << " tomorrowPer: " << tomorrowPer 00085 // << " " << phaseName( retPhase ); 00086 } else { 00087 //kDebug() << date 00088 // << " todayPer: " << todayPer 00089 // << " tomorrowPer: " << tomorrowPer 00090 // << " " << phaseName( None ); 00091 } 00092 00093 // Note: if you want to support crescent and gibbous phases then please 00094 // read the source for the original BSD 'pom' program. 00095 } 00096 00097 return( retPhase ); 00098 } 00099 00100 /* 00101 * Copyright (c) 1989, 1993 00102 * The Regents of the University of California. All rights reserved. 00103 * 00104 * This code is derived from software posted to USENET. 00105 * 00106 * Redistribution and use in source and binary forms, with or without 00107 * modification, are permitted provided that the following conditions 00108 * are met: 00109 * 1. Redistributions of source code must retain the above copyright 00110 * notice, this list of conditions and the following disclaimer. 00111 * 2. Redistributions in binary form must reproduce the above copyright 00112 * notice, this list of conditions and the following disclaimer in the 00113 * documentation and/or other materials provided with the distribution. 00114 * 3. All advertising materials mentioning features or use of this software 00115 * must display the following acknowledgement: 00116 * This product includes software developed by the University of 00117 * California, Berkeley and its contributors. 00118 * 4. Neither the name of the University nor the names of its contributors 00119 * may be used to endorse or promote products derived from this software 00120 * without specific prior written permission. 00121 * 00122 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 00123 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00124 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 00125 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 00126 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 00127 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 00128 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 00129 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 00130 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 00131 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 00132 * SUCH DAMAGE. 00133 */ 00134 00135 #ifdef HAVE_SYS_CDEFS_H 00136 #include <sys/cdefs.h> 00137 #endif 00138 00139 /* 00140 * Phase of the Moon. Calculates the current phase of the moon. 00141 * Based on routines from `Practical Astronomy with Your Calculator', 00142 * by Duffett-Smith. Comments give the section from the book that 00143 * particular piece of code was adapted from. 00144 * 00145 * -- Keith E. Brandt VIII 1984 00146 * 00147 * Updated to the Third Edition of Duffett-Smith's book, Paul Janzen, IX 1998 00148 * 00149 */ 00150 00151 #include <ctype.h> 00152 #ifdef HAVE_ERR_H 00153 #include <err.h> 00154 #endif 00155 #include <math.h> 00156 #include <string.h> 00157 #include <stdlib.h> 00158 #include <time.h> 00159 #include <unistd.h> 00160 00161 #ifndef PI 00162 #define PI 3.14159265358979323846 00163 #endif 00164 00165 /* 00166 * The EPOCH in the third edition of the book is 1990 Jan 0.0 TDT. 00167 * In this program, we do not bother to correct for the differences 00168 * between UTC (as shown by the UNIX clock) and TDT. (TDT = TAI + 32.184s; 00169 * TAI-UTC = 32s in Jan 1999.) 00170 */ 00171 #define EPOCH_MINUS_1970 (20 * 365 + 5 - 1) /* 20 years, 5 leaps, back 1 day to Jan 0 */ 00172 #define EPSILONg 279.403303 /* solar ecliptic long at EPOCH */ 00173 #define RHOg 282.768422 /* solar ecliptic long of perigee at EPOCH */ 00174 #define ECCEN 0.016713 /* solar orbit eccentricity */ 00175 #define lzero 318.351648 /* lunar mean long at EPOCH */ 00176 #define Pzero 36.340410 /* lunar mean long of perigee at EPOCH */ 00177 #define Nzero 318.510107 /* lunar mean long of node at EPOCH */ 00178 00179 /* 00180 * percentFull -- 00181 * return phase of the moon as a percentage of full 00182 */ 00183 static double percentFull( uint tmpt ) 00184 { 00185 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 00186 double A4, lprime, V, ldprime, D, Nm; 00187 00188 double days; 00189 days = ( tmpt - EPOCH_MINUS_1970 * 86400 ) / 86400.0; 00190 00191 N = 360 * days / 365.242191; /* sec 46 #3 */ 00192 adj360( &N ); 00193 Msol = N + EPSILONg - RHOg; /* sec 46 #4 */ 00194 adj360( &Msol ); 00195 Ec = 360 / PI * ECCEN * sin( degreesToRadians( Msol ) ); /* sec 46 #5 */ 00196 LambdaSol = N + Ec + EPSILONg; /* sec 46 #6 */ 00197 adj360( &LambdaSol ); 00198 l = 13.1763966 * days + lzero; /* sec 65 #4 */ 00199 adj360( &l ); 00200 Mm = l - ( 0.1114041 * days ) - Pzero; /* sec 65 #5 */ 00201 adj360( &Mm ); 00202 Nm = Nzero - ( 0.0529539 * days ); /* sec 65 #6 */ 00203 adj360( &Nm ); 00204 Ev = 1.2739 * sin( degreesToRadians( 2 * ( l - LambdaSol ) - Mm ) ); /* sec 65 #7 */ 00205 Ac = 0.1858 * sin( degreesToRadians( Msol ) ); /* sec 65 #8 */ 00206 A3 = 0.37 * sin( degreesToRadians( Msol ) ); 00207 Mmprime = Mm + Ev - Ac - A3; /* sec 65 #9 */ 00208 Ec = 6.2886 * sin( degreesToRadians( Mmprime ) ); /* sec 65 #10 */ 00209 A4 = 0.214 * sin( degreesToRadians( 2 * Mmprime ) ); /* sec 65 #11 */ 00210 lprime = l + Ev + Ec - Ac + A4; /* sec 65 #12 */ 00211 V = 0.6583 * sin( degreesToRadians( 2 * ( lprime - LambdaSol ) ) );/* sec 65 #13 */ 00212 ldprime = lprime + V; /* sec 65 #14 */ 00213 D = ldprime - LambdaSol; /* sec 67 #2 */ 00214 D = 50.0 * ( 1 - cos( degreesToRadians( D ) ) ); /* sec 67 #3 */ 00215 return D; 00216 } 00217 00218 /* 00219 * degreesToRadians -- 00220 * convert degrees to radians 00221 */ 00222 static double degreesToRadians( double degree ) 00223 { 00224 return ( degree * PI ) / 180.00; 00225 } 00226 00227 /* 00228 * adj360 -- 00229 * adjust value so 0 <= degree <= 360 00230 */ 00231 static void adj360( double *degree ) 00232 { 00233 for ( ;; ) { 00234 if ( *degree < 0 ) { 00235 *degree += 360; 00236 } else if ( *degree > 360 ) { 00237 *degree -= 360; 00238 } else { 00239 break; 00240 } 00241 } 00242 }