lunarphase.cpp
00001 /* 00002 This file is part of libkholidays. 00003 Copyright (c) 2004 Allen Winter <winter@kde.org> 00004 00005 Copyright (c) 1989, 1993 00006 The Regents of the University of California. All rights reserved. 00007 00008 This library is free software; you can redistribute it and/or 00009 modify it under the terms of the GNU General Public License as 00010 published by the Free Software Foundation; either version 2 of the 00011 License, or (at your option) any later version. 00012 00013 This library is distributed in the hope that it will be useful, 00014 but WITHOUT ANY WARRANTY; without even the implied warranty of 00015 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00016 General Public License for more details. 00017 00018 You should have received a copy of the GNU General Public License 00019 along with this program; if not, write to the Free Software 00020 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00021 00022 In addition, as a special exception, the copyright holders give 00023 permission to link the code of this program with any edition of 00024 the TQt library by Trolltech AS, Norway (or with modified versions 00025 of TQt that use the same license as TQt), and distribute linked 00026 combinations including the two. You must obey the GNU General 00027 Public License in all respects for all of the code used other than 00028 TQt. If you modify this file, you may extend this exception to 00029 your version of the file, but you are not obligated to do so. If 00030 you do not wish to do so, delete this exception statement from 00031 your version. 00032 */ 00033 00034 #include "config.h" 00035 00036 #include <kglobal.h> 00037 #include <klocale.h> 00038 #include <kdebug.h> 00039 00040 #include "lunarphase.h" 00041 00042 LunarPhase::LunarPhase( Hemisphere hemisphere ) 00043 { 00044 mHemisphere = hemisphere; 00045 } 00046 00047 LunarPhase::~LunarPhase() 00048 { 00049 } 00050 00051 void LunarPhase::setHemisphere( Hemisphere hemisphere ) 00052 { 00053 mHemisphere = hemisphere; 00054 } 00055 00056 LunarPhase::Hemisphere LunarPhase::hemisphere() const 00057 { 00058 return( mHemisphere ); 00059 } 00060 00061 TQString LunarPhase::hemisphereStr() const 00062 { 00063 return hemisphereName( mHemisphere ); 00064 } 00065 00066 TQString LunarPhase::hemisphereName( LunarPhase::Hemisphere hemisphere ) 00067 { 00068 switch( hemisphere ) { 00069 case Northern: 00070 default: 00071 return( i18n( "Northern" ) ); 00072 break; 00073 case Southern: 00074 return( i18n( "Southern" ) ); 00075 break; 00076 } 00077 } 00078 00079 TQString LunarPhase::phaseStr( const TQDate &date ) const 00080 { 00081 return phaseName( phase( date ) ); 00082 } 00083 00084 TQString LunarPhase::phaseName( LunarPhase::Phase phase ) 00085 { 00086 switch ( phase ) { 00087 case New: 00088 return( i18n( "New Moon" ) ); 00089 break; 00090 case Full: 00091 return( i18n( "Full Moon" ) ); 00092 break; 00093 case FirstQ: 00094 return( i18n( "First Quarter Moon" ) ); 00095 break; 00096 case LastQ: 00097 return( i18n( "Last Quarter Moon" ) ); 00098 break; 00099 default: 00100 case None: 00101 return( TQString() ); 00102 break; 00103 } 00104 } 00105 00106 LunarPhase::Phase LunarPhase::phase( const TQDate &date ) const 00107 { 00108 Phase retPhase = None; 00109 00110 // compute percent-full for the middle of today and yesterday. 00111 TQTime noontime( 12, 0, 0 ); 00112 TQDateTime today( date, noontime ); 00113 double todayPer = percentFull( today.toTime_t() ); 00114 TQDateTime yesterday( date.addDays(-1), noontime ); 00115 double yesterdayPer = percentFull( yesterday.toTime_t() ); 00116 00117 if ( ( todayPer < 0.50 ) && ( yesterdayPer > 0.50 ) ) { 00118 retPhase = New; 00119 } else if ( ( todayPer > 99.50 ) && ( yesterdayPer < 99.50 ) ) { 00120 retPhase = Full; 00121 } else { 00122 // compute percent-full for the start of today. 00123 TQTime sqt( 0, 0, 0 ); 00124 TQDateTime start( date, sqt ); 00125 double startPer = percentFull( start.toTime_t() ); 00126 // compute percent-full for the end of today. 00127 TQTime eqt( 23, 59, 59 ); 00128 TQDateTime end( date, eqt ); 00129 double endPer = percentFull( end.toTime_t() ); 00130 00131 if ( ( startPer <= 50 ) && ( endPer > 50 ) ) { 00132 if ( mHemisphere == Northern ) { 00133 retPhase = FirstQ; 00134 } else { 00135 retPhase = LastQ; 00136 } 00137 } 00138 if ( ( endPer <= 50 ) && ( startPer > 50 ) ) { 00139 if ( mHemisphere == Northern ) { 00140 retPhase = LastQ; 00141 } else { 00142 retPhase = FirstQ; 00143 } 00144 } 00145 // Note: if you want to support crescent and gibbous phases then please 00146 // read the source for the original BSD 'pom' program. 00147 } 00148 return( retPhase ); 00149 } 00150 00151 /* 00152 * Copyright (c) 1989, 1993 00153 * The Regents of the University of California. All rights reserved. 00154 * 00155 * This code is derived from software posted to USENET. 00156 * 00157 * Redistribution and use in source and binary forms, with or without 00158 * modification, are permitted provided that the following conditions 00159 * are met: 00160 * 1. Redistributions of source code must retain the above copyright 00161 * notice, this list of conditions and the following disclaimer. 00162 * 2. Redistributions in binary form must reproduce the above copyright 00163 * notice, this list of conditions and the following disclaimer in the 00164 * documentation and/or other materials provided with the distribution. 00165 * 3. All advertising materials mentioning features or use of this software 00166 * must display the following acknowledgement: 00167 * This product includes software developed by the University of 00168 * California, Berkeley and its contributors. 00169 * 4. Neither the name of the University nor the names of its contributors 00170 * may be used to endorse or promote products derived from this software 00171 * without specific prior written permission. 00172 * 00173 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 00174 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00175 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 00176 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 00177 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 00178 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 00179 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 00180 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 00181 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 00182 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 00183 * SUCH DAMAGE. 00184 */ 00185 00186 #if HAVE_SYS_CDEFS_H 00187 #include <sys/cdefs.h> 00188 #endif 00189 00190 /* 00191 * Phase of the Moon. Calculates the current phase of the moon. 00192 * Based on routines from `Practical Astronomy with Your Calculator', 00193 * by Duffett-Smith. Comments give the section from the book that 00194 * particular piece of code was adapted from. 00195 * 00196 * -- Keith E. Brandt VIII 1984 00197 * 00198 * Updated to the Third Edition of Duffett-Smith's book, Paul Janzen, IX 1998 00199 * 00200 */ 00201 00202 00203 #include <ctype.h> 00204 #if HAVE_ERR_H 00205 #include <err.h> 00206 #endif 00207 #include <math.h> 00208 #include <string.h> 00209 #include <stdlib.h> 00210 #include <time.h> 00211 #include <unistd.h> 00212 00213 #ifndef PI 00214 #define PI 3.14159265358979323846 00215 #endif 00216 00217 /* 00218 * The EPOCH in the third edition of the book is 1990 Jan 0.0 TDT. 00219 * In this program, we do not bother to correct for the differences 00220 * between UTC (as shown by the UNIX clock) and TDT. (TDT = TAI + 32.184s; 00221 * TAI-UTC = 32s in Jan 1999.) 00222 */ 00223 #define EPOCH_MINUS_1970 (20 * 365 + 5 - 1) /* 20 years, 5 leaps, back 1 day to Jan 0 */ 00224 #define EPSILONg 279.403303 /* solar ecliptic long at EPOCH */ 00225 #define RHOg 282.768422 /* solar ecliptic long of perigee at EPOCH */ 00226 #define ECCEN 0.016713 /* solar orbit eccentricity */ 00227 #define lzero 318.351648 /* lunar mean long at EPOCH */ 00228 #define Pzero 36.340410 /* lunar mean long of perigee at EPOCH */ 00229 #define Nzero 318.510107 /* lunar mean long of node at EPOCH */ 00230 00231 /* 00232 * percentFull -- 00233 * return phase of the moon as a percentage of full 00234 */ 00235 double LunarPhase::percentFull( uint tmpt ) const 00236 { 00237 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 00238 double A4, lprime, V, ldprime, D, Nm; 00239 00240 double days; 00241 days = ( tmpt - EPOCH_MINUS_1970 * 86400 ) / 86400.0; 00242 00243 N = 360 * days / 365.242191; /* sec 46 #3 */ 00244 adj360(&N); 00245 Msol = N + EPSILONg - RHOg; /* sec 46 #4 */ 00246 adj360(&Msol); 00247 Ec = 360 / PI * ECCEN * sin(degreesToRadians(Msol)); /* sec 46 #5 */ 00248 LambdaSol = N + Ec + EPSILONg; /* sec 46 #6 */ 00249 adj360(&LambdaSol); 00250 l = 13.1763966 * days + lzero; /* sec 65 #4 */ 00251 adj360(&l); 00252 Mm = l - (0.1114041 * days) - Pzero; /* sec 65 #5 */ 00253 adj360(&Mm); 00254 Nm = Nzero - (0.0529539 * days); /* sec 65 #6 */ 00255 adj360(&Nm); 00256 Ev = 1.2739 * sin(degreesToRadians(2*(l - LambdaSol) - Mm)); /* sec 65 #7 */ 00257 Ac = 0.1858 * sin(degreesToRadians(Msol)); /* sec 65 #8 */ 00258 A3 = 0.37 * sin(degreesToRadians(Msol)); 00259 Mmprime = Mm + Ev - Ac - A3; /* sec 65 #9 */ 00260 Ec = 6.2886 * sin(degreesToRadians(Mmprime)); /* sec 65 #10 */ 00261 A4 = 0.214 * sin(degreesToRadians(2 * Mmprime)); /* sec 65 #11 */ 00262 lprime = l + Ev + Ec - Ac + A4; /* sec 65 #12 */ 00263 V = 0.6583 * sin(degreesToRadians(2 * (lprime - LambdaSol))); /* sec 65 #13 */ 00264 ldprime = lprime + V; /* sec 65 #14 */ 00265 D = ldprime - LambdaSol; /* sec 67 #2 */ 00266 return(50.0 * (1 - cos(degreesToRadians(D)))); /* sec 67 #3 */ 00267 } 00268 00269 /* 00270 * degreesToRadians -- 00271 * convert degrees to radians 00272 */ 00273 double LunarPhase::degreesToRadians( double degree ) const 00274 { 00275 return( degree * PI / 180 ); 00276 } 00277 00278 /* 00279 * adj360 -- 00280 * adjust value so 0 <= degree <= 360 00281 */ 00282 void LunarPhase::adj360( double *degree ) const 00283 { 00284 for( ;; ) 00285 if( *degree < 0 ) 00286 *degree += 360; 00287 else if( *degree > 360 ) 00288 *degree -= 360; 00289 else 00290 break; 00291 }