1/*
2 This file is part of the kholidays library.
3
4 Copyright (c) 2004,2007,2009 Allen Winter <winter@kde.org>
5
6 Copyright (c) 1989, 1993 //krazy:exclude=copyright
7 The Regents of the University of California. All rights reserved.
8
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Library General Public
11 License as published by the Free Software Foundation; either
12 version 2 of the License, or (at your option) any later version.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU Library General Public License for more details.
18
19 You should have received a copy of the GNU Library General Public License
20 along with this library; see the file COPYING.LIB. If not, write to the
21 Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 Boston, MA 02110-1301, USA.
23*/
24
25#include "lunarphase.h"
26#include <config-kholidays.h>
27
28#include <KDebug>
29#include <KGlobal>
30#include <KLocalizedString>
31
32#include <QtCore/QDateTime>
33
34using namespace KHolidays;
35
36static double percentFull( uint tmpt );
37static double degreesToRadians( double degree );
38static void adj360( double *degree );
39
40QString LunarPhase::phaseNameAtDate( const QDate &date )
41{
42 return phaseName( phaseAtDate( date ) );
43}
44
45QString LunarPhase::phaseName( LunarPhase::Phase phase )
46{
47 switch ( phase ) {
48 case NewMoon:
49 return( i18n( "New Moon" ) );
50 case FullMoon:
51 return( i18n( "Full Moon" ) );
52 case FirstQuarter:
53 return( i18n( "First Quarter Moon" ) );
54 case LastQuarter:
55 return( i18n( "Last Quarter Moon" ) );
56 default:
57 case None:
58 return QString();
59 }
60}
61
62LunarPhase::Phase LunarPhase::phaseAtDate( const QDate &date )
63{
64 Phase retPhase = None;
65
66 // compute percent-full for the middle of today and yesterday.
67 const QTime anytime( 12, 0, 0 );
68 const QDateTime today( date, anytime, Qt::UTC );
69 const double todayPer = percentFull( today.toTime_t() ) + 0.5;
70
71 const QDateTime tomorrow( date.addDays( 1 ), anytime, Qt::UTC );
72 const double tomorrowPer = percentFull( tomorrow.toTime_t() ) + 0.5;
73
74 if ( static_cast<int>( todayPer ) == 100 &&
75 static_cast<int>( tomorrowPer ) != 100 ) {
76 retPhase = FullMoon;
77 } else if ( static_cast<int>( todayPer ) == 0 &&
78 static_cast<int>( tomorrowPer ) != 0 ) {
79 retPhase = NewMoon;
80 } else {
81 if ( todayPer > 50 && tomorrowPer < 50 ) {
82 retPhase = LastQuarter;
83 }
84 if ( todayPer < 50 && tomorrowPer > 50 ) {
85 retPhase = FirstQuarter;
86 }
87
88 // Note: if you want to support crescent and gibbous phases then please
89 // read the source for the original BSD 'pom' program.
90 }
91
92 return( retPhase );
93}
94
95/*
96 * Copyright (c) 1989, 1993
97 * The Regents of the University of California. All rights reserved.
98 *
99 * This code is derived from software posted to USENET.
100 *
101 * Redistribution and use in source and binary forms, with or without
102 * modification, are permitted provided that the following conditions
103 * are met:
104 * 1. Redistributions of source code must retain the above copyright
105 * notice, this list of conditions and the following disclaimer.
106 * 2. Redistributions in binary form must reproduce the above copyright
107 * notice, this list of conditions and the following disclaimer in the
108 * documentation and/or other materials provided with the distribution.
109 * 3. All advertising materials mentioning features or use of this software
110 * must display the following acknowledgement:
111 * This product includes software developed by the University of
112 * California, Berkeley and its contributors.
113 * 4. Neither the name of the University nor the names of its contributors
114 * may be used to endorse or promote products derived from this software
115 * without specific prior written permission.
116 *
117 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
118 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
119 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
120 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
121 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
122 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
123 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
124 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
125 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
126 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
127 * SUCH DAMAGE.
128 */
129
130#ifdef HAVE_SYS_CDEFS_H
131#include <sys/cdefs.h>
132#endif
133
134/*
135 * Phase of the Moon. Calculates the current phase of the moon.
136 * Based on routines from `Practical Astronomy with Your Calculator',
137 * by Duffett-Smith. Comments give the section from the book that
138 * particular piece of code was adapted from.
139 *
140 * -- Keith E. Brandt VIII 1984
141 *
142 * Updated to the Third Edition of Duffett-Smith's book, Paul Janzen, IX 1998
143 *
144 */
145
146#include <ctype.h>
147#ifdef HAVE_ERR_H
148#include <err.h>
149#endif
150#include <math.h>
151#include <string.h>
152#include <stdlib.h>
153#include <time.h>
154#include <unistd.h>
155
156static double PI = 3.14159265358979323846;
157
158/*
159 * The EPOCH in the third edition of the book is 1990 Jan 0.0 TDT.
160 * In this program, we do not bother to correct for the differences
161 * between UTC (as shown by the UNIX clock) and TDT. (TDT = TAI + 32.184s;
162 * TAI-UTC = 32s in Jan 1999.)
163 */
164static int EPOCH_MINUS_1970 = ( 20 * 365 + 5 - 1 ); /* 20 years, 5 leaps, back 1 day to Jan 0 */
165static double EPSILONg = 279.403303; /* solar ecliptic long at EPOCH */
166static double RHOg = 282.768422; /* solar ecliptic long of perigee at EPOCH */
167static double ECCEN = 0.016713; /* solar orbit eccentricity */
168static double lzero = 318.351648; /* lunar mean long at EPOCH */
169static double Pzero = 36.340410; /* lunar mean long of perigee at EPOCH */
170static double Nzero = 318.510107; /* lunar mean long of node at EPOCH */
171
172/*
173 * percentFull --
174 * return phase of the moon as a percentage of full
175 */
176static double percentFull( uint tmpt )
177{
178 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
179 double A4, lprime, V, ldprime, D, Nm;
180
181 double days;
182 days = ( tmpt - EPOCH_MINUS_1970 * 86400 ) / 86400.0;
183
184 N = 360 * days / 365.242191; /* sec 46 #3 */
185 adj360( &N );
186 Msol = N + EPSILONg - RHOg; /* sec 46 #4 */
187 adj360( &Msol );
188 Ec = 360 / PI * ECCEN * sin( degreesToRadians( Msol ) ); /* sec 46 #5 */
189 LambdaSol = N + Ec + EPSILONg; /* sec 46 #6 */
190 adj360( &LambdaSol );
191 l = 13.1763966 * days + lzero; /* sec 65 #4 */
192 adj360( &l );
193 Mm = l - ( 0.1114041 * days ) - Pzero; /* sec 65 #5 */
194 adj360( &Mm );
195 Nm = Nzero - ( 0.0529539 * days ); /* sec 65 #6 */
196 adj360( &Nm );
197 Ev = 1.2739 * sin( degreesToRadians( 2 * ( l - LambdaSol ) - Mm ) ); /* sec 65 #7 */
198 Ac = 0.1858 * sin( degreesToRadians( Msol ) ); /* sec 65 #8 */
199 A3 = 0.37 * sin( degreesToRadians( Msol ) );
200 Mmprime = Mm + Ev - Ac - A3; /* sec 65 #9 */
201 Ec = 6.2886 * sin( degreesToRadians( Mmprime ) ); /* sec 65 #10 */
202 A4 = 0.214 * sin( degreesToRadians( 2 * Mmprime ) ); /* sec 65 #11 */
203 lprime = l + Ev + Ec - Ac + A4; /* sec 65 #12 */
204 V = 0.6583 * sin( degreesToRadians( 2 * ( lprime - LambdaSol ) ) );/* sec 65 #13 */
205 ldprime = lprime + V; /* sec 65 #14 */
206 D = ldprime - LambdaSol; /* sec 67 #2 */
207 D = 50.0 * ( 1 - cos( degreesToRadians( D ) ) ); /* sec 67 #3 */
208 return D;
209}
210
211/*
212 * degreesToRadians --
213 * convert degrees to radians
214 */
215static double degreesToRadians( double degree )
216{
217 return ( degree * PI ) / 180.00;
218}
219
220/*
221 * adj360 --
222 * adjust value so 0 <= degree <= 360
223 */
224static void adj360( double *degree )
225{
226 for ( ;; ) {
227 if ( *degree < 0 ) {
228 *degree += 360;
229 } else if ( *degree > 360 ) {
230 *degree -= 360;
231 } else {
232 break;
233 }
234 }
235}
236