1 |
/*
|
2 |
* Copyright 2008-2009 Mike Reedell / LuckyCatLabs.
|
3 |
*
|
4 |
* Licensed under the Apache License, Version 2.0 (the "License");
|
5 |
* you may not use this file except in compliance with the License.
|
6 |
* You may obtain a copy of the License at
|
7 |
*
|
8 |
* http://www.apache.org/licenses/LICENSE-2.0
|
9 |
*
|
10 |
* Unless required by applicable law or agreed to in writing, software
|
11 |
* distributed under the License is distributed on an "AS IS" BASIS,
|
12 |
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
13 |
* See the License for the specific language governing permissions and
|
14 |
* limitations under the License.
|
15 |
*/
|
16 |
|
17 |
package com.luckycatlabs.sunrisesunset.calculator;
|
18 |
|
19 |
import java.math.BigDecimal;
|
20 |
import java.math.RoundingMode;
|
21 |
import java.util.Calendar;
|
22 |
import java.util.TimeZone;
|
23 |
|
24 |
import com.luckycatlabs.sunrisesunset.Zenith;
|
25 |
import com.luckycatlabs.sunrisesunset.dto.Location2;
|
26 |
|
27 |
/**
|
28 |
* Parent class of the Sunrise and Sunset calculator classes.
|
29 |
*/
|
30 |
public class SolarEventCalculator
|
31 |
{
|
32 |
private Location2 location;
|
33 |
private TimeZone timeZone;
|
34 |
|
35 |
/**
|
36 |
* Constructs a new <code>SolarEventCalculator</code> using the given parameters.
|
37 |
*
|
38 |
* @param location
|
39 |
* <code>Location</code> of the place where the solar event should be calculated from.
|
40 |
* @param timeZoneIdentifier
|
41 |
* time zone identifier of the timezone of the location parameter. For example,
|
42 |
* "America/New_York".
|
43 |
*/
|
44 |
public SolarEventCalculator(Location2 location, String timeZoneIdentifier)
|
45 |
{
|
46 |
this.location = location;
|
47 |
this.timeZone = TimeZone.getTimeZone(timeZoneIdentifier);
|
48 |
}
|
49 |
|
50 |
public double div(double a, double b)
|
51 |
{
|
52 |
int x;
|
53 |
int y;
|
54 |
x = (int) a;
|
55 |
y = (int) b;
|
56 |
return (double) (x / y);
|
57 |
}
|
58 |
|
59 |
public double Rev(double input)
|
60 |
{
|
61 |
double x;
|
62 |
x = input - Math.floor(input / 360.0) * 360;
|
63 |
return (x);
|
64 |
}
|
65 |
|
66 |
public double Radians(double input)
|
67 |
{
|
68 |
return Math.toRadians(input);
|
69 |
}
|
70 |
|
71 |
public double Deg(double input)
|
72 |
{
|
73 |
return Math.toDegrees(input);
|
74 |
}
|
75 |
|
76 |
public double ElevationRefraction(double El_geometric)
|
77 |
{
|
78 |
double El_observed;
|
79 |
double x, a0, a1, a2, a3, a4;
|
80 |
a0 = 0.58804392;
|
81 |
a1 = -0.17941557;
|
82 |
a2 = 0.29906946E-1;
|
83 |
a3 = -0.25187400E-2;
|
84 |
a4 = 0.82622101E-4;
|
85 |
El_observed = El_geometric;
|
86 |
x = Math.abs(El_geometric + 0.589);
|
87 |
double refraction = Math.abs(a0 + a1 * x + a2 * x * x + a3 * x * x * x + a4 * x * x * x * x);
|
88 |
|
89 |
if (El_geometric > 10.2)
|
90 |
{
|
91 |
El_observed = El_geometric + 0.01617 * (Math.cos(Radians(Math.abs(El_geometric))) / Math.sin(Radians(Math.abs(El_geometric))));
|
92 |
}
|
93 |
else
|
94 |
{
|
95 |
El_observed = El_geometric + refraction;
|
96 |
|
97 |
}
|
98 |
return El_observed;
|
99 |
}
|
100 |
|
101 |
public double CalcJD(double day, double month, double year)
|
102 |
{
|
103 |
double jd = 2415020.5 - 64; // 1.1.1900 - correction of algorithm
|
104 |
if (month <= 2)
|
105 |
{
|
106 |
year = year - 1;
|
107 |
month += 12;
|
108 |
}
|
109 |
jd += (double) ((int) ((year - 1900) * 365.25));
|
110 |
jd += (double) ((int) (30.6001 * (1 + month)));
|
111 |
return (jd + day);
|
112 |
}
|
113 |
|
114 |
public double frac(double x)
|
115 |
{
|
116 |
return (x - Math.floor(x));
|
117 |
}
|
118 |
|
119 |
public double Mod(double a, double b)
|
120 |
{
|
121 |
return (a - Math.floor(a / b) * b);
|
122 |
}
|
123 |
|
124 |
public double GMST(double JD)
|
125 |
{
|
126 |
double UT = frac(JD - 0.5) * 24.; // UT in hours
|
127 |
JD = Math.floor(JD - 0.5) + 0.5; // JD at 0 hours UT
|
128 |
double T = (JD - 2451545.0) / 36525.0;
|
129 |
double T0 = 6.697374558 + T * (2400.051336 + T * 0.000025862);
|
130 |
return (Mod(T0 + UT * 1.002737909, 24.));
|
131 |
}
|
132 |
|
133 |
public double GMST2LMST(double gmst, double lon)
|
134 |
{
|
135 |
double RAD = 180.0 / Math.PI;
|
136 |
double lmst = Mod(gmst + RAD * lon / 15, 24.);
|
137 |
return (lmst);
|
138 |
}
|
139 |
|
140 |
public class cart_ret
|
141 |
{
|
142 |
double x;
|
143 |
double y;
|
144 |
double z;
|
145 |
double radius;
|
146 |
double lat;
|
147 |
double lon;
|
148 |
}
|
149 |
|
150 |
public cart_ret EquPolar2Cart(double lon, double lat, double distance)
|
151 |
{
|
152 |
cart_ret cart = new cart_ret();
|
153 |
double rcd = Math.cos(lat) * distance;
|
154 |
cart.x = rcd * Math.cos(lon);
|
155 |
cart.y = rcd * Math.sin(lon);
|
156 |
cart.z = distance * Math.sin(lat);
|
157 |
return (cart);
|
158 |
}
|
159 |
|
160 |
public cart_ret Observer2EquCart(double lon, double lat, double height, double gmst)
|
161 |
{
|
162 |
double flat = 298.257223563; // WGS84 flatening of earth
|
163 |
double aearth = 6378.137; // GRS80/WGS84 semi major axis of earth ellipsoid
|
164 |
cart_ret cart = new cart_ret();
|
165 |
// Calculate geocentric latitude from geodetic latitude
|
166 |
double co = Math.cos(lat);
|
167 |
double si = Math.sin(lat);
|
168 |
double fl = 1.0 - 1.0 / flat;
|
169 |
fl = fl * fl;
|
170 |
si = si * si;
|
171 |
double u = 1.0 / Math.sqrt(co * co + fl * si);
|
172 |
double a = aearth * u + height;
|
173 |
double b = aearth * fl * u + height;
|
174 |
double radius = Math.sqrt(a * a * co * co + b * b * si); // geocentric distance from earth center
|
175 |
cart.y = Math.acos(a * co / radius); // geocentric latitude, rad
|
176 |
cart.x = lon; // longitude stays the same
|
177 |
if (lat < 0.0)
|
178 |
{
|
179 |
cart.y = -cart.y;
|
180 |
} // adjust sign
|
181 |
cart = EquPolar2Cart(cart.x, cart.y, radius); // convert from geocentric polar to geocentric cartesian, with regard to Greenwich
|
182 |
// rotate around earth's polar axis to align coordinate system from Greenwich to vernal equinox
|
183 |
double x = cart.x;
|
184 |
double y = cart.y;
|
185 |
double rotangle = gmst / 24 * 2 * Math.PI; // sideral time gmst given in hours. Convert to radians
|
186 |
cart.x = x * Math.cos(rotangle) - y * Math.sin(rotangle);
|
187 |
cart.y = x * Math.sin(rotangle) + y * Math.cos(rotangle);
|
188 |
cart.radius = radius;
|
189 |
cart.lon = lon;
|
190 |
cart.lat = lat;
|
191 |
return (cart);
|
192 |
}
|
193 |
|
194 |
public class moonCoor_ret
|
195 |
{
|
196 |
public double lon;
|
197 |
public double lat;
|
198 |
public double orbitLon;
|
199 |
public double distance;
|
200 |
public double diameter;
|
201 |
public double parallax;
|
202 |
public double raGeocentric;
|
203 |
public double decGeocentric;
|
204 |
public double ra;
|
205 |
public double dec;
|
206 |
public double raTopocentric;
|
207 |
public double decTopocentric;
|
208 |
public double distanceTopocentric;
|
209 |
public double moonAge;
|
210 |
public double phase;
|
211 |
public double az;
|
212 |
public double alt;
|
213 |
public String moonPhase;
|
214 |
public String sign;
|
215 |
}
|
216 |
|
217 |
public double Mod2Pi(double x)
|
218 |
{
|
219 |
x = Mod(x, 2. * Math.PI);
|
220 |
return (x);
|
221 |
}
|
222 |
|
223 |
public moonCoor_ret Ecl2Equ(moonCoor_ret coor, double TDT)
|
224 |
{
|
225 |
double pi = Math.PI;
|
226 |
double DEG = pi / 180.0;
|
227 |
|
228 |
double T = (TDT - 2451545.0) / 36525.; // Epoch 2000 January 1.5
|
229 |
double eps = (23. + (26 + 21.45 / 60.) / 60. + T * (-46.815 + T * (-0.0006 + T * 0.00181)) / 3600.) * DEG;
|
230 |
double coseps = Math.cos(eps);
|
231 |
double sineps = Math.sin(eps);
|
232 |
|
233 |
double sinlon = Math.sin(coor.lon);
|
234 |
coor.ra = Mod2Pi(Math.atan2((sinlon * coseps - Math.tan(coor.lat) * sineps), Math.cos(coor.lon)));
|
235 |
coor.dec = Math.asin(Math.sin(coor.lat) * coseps + Math.cos(coor.lat) * sineps * sinlon);
|
236 |
|
237 |
return coor;
|
238 |
}
|
239 |
|
240 |
public moonCoor_ret GeoEqu2TopoEqu(moonCoor_ret coor, cart_ret observer, double lmst)
|
241 |
{
|
242 |
double cosdec = Math.cos(coor.dec);
|
243 |
double sindec = Math.sin(coor.dec);
|
244 |
double coslst = Math.cos(lmst);
|
245 |
double sinlst = Math.sin(lmst);
|
246 |
double coslat = Math.cos(observer.lat); // we should use geocentric latitude, not geodetic latitude
|
247 |
double sinlat = Math.sin(observer.lat);
|
248 |
double rho = observer.radius; // observer-geocenter in Kilometer
|
249 |
|
250 |
double x = coor.distance * cosdec * Math.cos(coor.ra) - rho * coslat * coslst;
|
251 |
double y = coor.distance * cosdec * Math.sin(coor.ra) - rho * coslat * sinlst;
|
252 |
double z = coor.distance * sindec - rho * sinlat;
|
253 |
|
254 |
coor.distanceTopocentric = Math.sqrt(x * x + y * y + z * z);
|
255 |
coor.decTopocentric = Math.asin(z / coor.distanceTopocentric);
|
256 |
coor.raTopocentric = Mod2Pi(Math.atan2(y, x));
|
257 |
|
258 |
return coor;
|
259 |
}
|
260 |
|
261 |
public moonCoor_ret Equ2Altaz(moonCoor_ret coor, double TDT, double geolat, double lmst)
|
262 |
{
|
263 |
double cosdec = Math.cos(coor.dec);
|
264 |
double sindec = Math.sin(coor.dec);
|
265 |
double lha = lmst - coor.ra;
|
266 |
double coslha = Math.cos(lha);
|
267 |
double sinlha = Math.sin(lha);
|
268 |
double coslat = Math.cos(geolat);
|
269 |
double sinlat = Math.sin(geolat);
|
270 |
|
271 |
double N = -cosdec * sinlha;
|
272 |
double D = sindec * coslat - cosdec * coslha * sinlat;
|
273 |
coor.az = Mod2Pi(Math.atan2(N, D));
|
274 |
coor.alt = Math.asin(sindec * sinlat + cosdec * coslha * coslat);
|
275 |
|
276 |
return coor;
|
277 |
}
|
278 |
|
279 |
public class sunCoor_ret
|
280 |
{
|
281 |
double lon;
|
282 |
double anomalyMean;
|
283 |
}
|
284 |
|
285 |
public moonCoor_ret MoonPosition(sunCoor_ret sunCoor, double TDT, cart_ret observer, double lmst)
|
286 |
{
|
287 |
double D = TDT - 2447891.5;
|
288 |
|
289 |
double pi = Math.PI;
|
290 |
double DEG = pi / 180.0;
|
291 |
|
292 |
// Mean Moon orbit elements as of 1990.0
|
293 |
double l0 = 318.351648 * DEG;
|
294 |
double P0 = 36.340410 * DEG;
|
295 |
double N0 = 318.510107 * DEG;
|
296 |
double i = 5.145396 * DEG;
|
297 |
double e = 0.054900;
|
298 |
double a = 384401; // km
|
299 |
double diameter0 = 0.5181 * DEG; // angular diameter of Moon at a distance
|
300 |
double parallax0 = 0.9507 * DEG; // parallax at distance a
|
301 |
|
302 |
double l = 13.1763966 * DEG * D + l0;
|
303 |
double MMoon = l - 0.1114041 * DEG * D - P0; // Moon's mean anomaly M
|
304 |
double N = N0 - 0.0529539 * DEG * D; // Moon's mean ascending node longitude
|
305 |
double C = l - sunCoor.lon;
|
306 |
double Ev = 1.2739 * DEG * Math.sin(2 * C - MMoon);
|
307 |
double Ae = 0.1858 * DEG * Math.sin(sunCoor.anomalyMean);
|
308 |
double A3 = 0.37 * DEG * Math.sin(sunCoor.anomalyMean);
|
309 |
double MMoon2 = MMoon + Ev - Ae - A3; // corrected Moon anomaly
|
310 |
double Ec = 6.2886 * DEG * Math.sin(MMoon2); // equation of centre
|
311 |
double A4 = 0.214 * DEG * Math.sin(2 * MMoon2);
|
312 |
double l2 = l + Ev + Ec - Ae + A4; // corrected Moon's longitude
|
313 |
double V = 0.6583 * DEG * Math.sin(2 * (l2 - sunCoor.lon));
|
314 |
double l3 = l2 + V; // true orbital longitude;
|
315 |
|
316 |
double N2 = N - 0.16 * DEG * Math.sin(sunCoor.anomalyMean);
|
317 |
|
318 |
moonCoor_ret moonCoor = new moonCoor_ret();
|
319 |
moonCoor.lon = Mod2Pi(N2 + Math.atan2(Math.sin(l3 - N2) * Math.cos(i), Math.cos(l3 - N2)));
|
320 |
moonCoor.lat = Math.asin(Math.sin(l3 - N2) * Math.sin(i));
|
321 |
moonCoor.orbitLon = l3;
|
322 |
|
323 |
moonCoor = Ecl2Equ(moonCoor, TDT);
|
324 |
// relative distance to semi mayor axis of lunar oribt
|
325 |
moonCoor.distance = (1 - Math.sqrt(e)) / (1 + e * Math.cos(MMoon2 + Ec));
|
326 |
moonCoor.diameter = diameter0 / moonCoor.distance; // angular diameter in radians
|
327 |
moonCoor.parallax = parallax0 / moonCoor.distance; // horizontal parallax in radians
|
328 |
moonCoor.distance *= a; // distance in km
|
329 |
|
330 |
// Calculate horizonal coordinates of sun, if geographic positions is given
|
331 |
if ((observer != null) && (lmst != 0.0))
|
332 |
{
|
333 |
// transform geocentric coordinates into topocentric (==observer based) coordinates
|
334 |
moonCoor = GeoEqu2TopoEqu(moonCoor, observer, lmst);
|
335 |
moonCoor.raGeocentric = moonCoor.ra; // backup geocentric coordinates
|
336 |
moonCoor.decGeocentric = moonCoor.dec;
|
337 |
moonCoor.ra = moonCoor.raTopocentric;
|
338 |
moonCoor.dec = moonCoor.decTopocentric;
|
339 |
moonCoor = Equ2Altaz(moonCoor, TDT, observer.lat, lmst); // now ra and dec are topocentric
|
340 |
}
|
341 |
|
342 |
// Age of Moon in radians since New Moon (0) - Full Moon (pi)
|
343 |
moonCoor.moonAge = Mod2Pi(l3 - sunCoor.lon);
|
344 |
moonCoor.phase = 0.5 * (1 - Math.cos(moonCoor.moonAge)); // Moon phase, 0-1
|
345 |
|
346 |
//String[] phases = {"Neumond", "Zunehmende Sichel", "Erstes Viertel", "Zunehmender Mond",
|
347 |
// "Vollmond", "Abnehmender Mond", "Letztes Viertel", "Abnehmende Sichel", "Neumond"};
|
348 |
double mainPhase = 1. / 29.53 * 360 * DEG; // show 'Newmoon, 'Quarter' for +/-1 day arond the actual event
|
349 |
double p = Mod(moonCoor.moonAge, 90. * DEG);
|
350 |
if (p < mainPhase || p > 90 * DEG - mainPhase)
|
351 |
p = 2 * Math.round(moonCoor.moonAge / (90. * DEG));
|
352 |
else
|
353 |
p = 2 * Math.floor(moonCoor.moonAge / (90. * DEG)) + 1;
|
354 |
// moonCoor.moonPhase = phases[p];
|
355 |
|
356 |
// moonCoor.sign = Sign(moonCoor.lon);
|
357 |
|
358 |
return (moonCoor);
|
359 |
}
|
360 |
|
361 |
public double Refraction(double alt)
|
362 |
{
|
363 |
double pi = Math.PI;
|
364 |
double DEG = pi / 180.0;
|
365 |
double RAD = 180.0 / pi;
|
366 |
|
367 |
double altdeg = alt * RAD;
|
368 |
if (altdeg < -2 || altdeg >= 90) return (0);
|
369 |
|
370 |
double pressure = 1015;
|
371 |
double temperature = 10;
|
372 |
if (altdeg > 15) return (0.00452 * pressure / ((273 + temperature) * Math.tan(alt)));
|
373 |
|
374 |
double y = alt;
|
375 |
double D = 0.0;
|
376 |
double P = (pressure - 80.) / 930.;
|
377 |
double Q = 0.0048 * (temperature - 10.);
|
378 |
double y0 = y;
|
379 |
double D0 = D;
|
380 |
|
381 |
for (int i = 0; i < 3; i++)
|
382 |
{
|
383 |
double N = y + (7.31 / (y + 4.4));
|
384 |
N = 1. / Math.tan(N * DEG);
|
385 |
D = N * P / (60. + Q * (N + 39.));
|
386 |
N = y - y0;
|
387 |
y0 = D - D0 - N;
|
388 |
if ((N != 0.) && (y0 != 0.))
|
389 |
{
|
390 |
N = y - N * (alt + D - y) / y0;
|
391 |
}
|
392 |
else
|
393 |
{
|
394 |
N = alt + D;
|
395 |
}
|
396 |
y0 = y;
|
397 |
D0 = D;
|
398 |
y = N;
|
399 |
}
|
400 |
return (D); // Hebung durch Refraktion in radians
|
401 |
}
|
402 |
|
403 |
public moonCoor_ret computeMoonStats(Calendar date)
|
404 |
{
|
405 |
double pi = Math.PI;
|
406 |
double DEG = pi / 180.0;
|
407 |
double RAD = 180.0 / pi;
|
408 |
|
409 |
date.setTimeZone(TimeZone.getTimeZone("UTC"));
|
410 |
|
411 |
double JD0 = CalcJD((double) date.get(Calendar.DAY_OF_MONTH), (double) date.get(Calendar.MONTH) + 1, (double) date.get(Calendar.YEAR));
|
412 |
double JD = JD0 + ((double) date.get(Calendar.HOUR_OF_DAY) + (double) date.get(Calendar.MINUTE) / 60 + (double) date.get(Calendar.SECOND) / 3600) / 24;
|
413 |
// UTC
|
414 |
double TDT = JD;
|
415 |
//System.out.println("TDT=" + TDT);
|
416 |
|
417 |
double lat = this.location.getLatitude().doubleValue() * DEG; // geodetic latitude of observer on WGS84
|
418 |
double lon = this.location.getLongitude().doubleValue() * DEG; // latitude of observer
|
419 |
double height = 0 * 0.001; // altiude of observer in meters above WGS84 ellipsoid (and converted to kilometers)
|
420 |
double gmst = GMST(JD);
|
421 |
double lmst = GMST2LMST(gmst, lon);
|
422 |
|
423 |
cart_ret observerCart = Observer2EquCart(lon, lat, height, gmst); // geocentric cartesian coordinates of observer
|
424 |
|
425 |
// sunCoor = SunPosition(TDT, lat, lmst*15.*DEG); // Calculate data for the Sun at given time
|
426 |
sunCoor_ret sunCoor = new sunCoor_ret();
|
427 |
BigDecimal longitudeHour = getLongitudeHour(date, true);
|
428 |
BigDecimal meanAnomaly = getMeanAnomaly(longitudeHour);
|
429 |
sunCoor.lon = longitudeHour.doubleValue();
|
430 |
sunCoor.anomalyMean = meanAnomaly.doubleValue();
|
431 |
moonCoor_ret moonCoor = MoonPosition(sunCoor, TDT, observerCart, lmst * 15. * DEG); // Calculate data for the Moon at given time
|
432 |
moonCoor.az = moonCoor.az * RAD;
|
433 |
moonCoor.alt = moonCoor.alt * RAD + Refraction(moonCoor.alt);
|
434 |
//System.out.println("moon azimuth=" + moonCoor.az);
|
435 |
//System.out.println("moon elevation=" + moonCoor.alt);
|
436 |
|
437 |
return moonCoor;
|
438 |
}
|
439 |
|
440 |
public void computeMoonStats2(Calendar date)
|
441 |
{
|
442 |
date.setTimeZone(TimeZone.getTimeZone("UTC"));
|
443 |
|
444 |
double Year = date.get(Calendar.YEAR);
|
445 |
double Month = date.get(Calendar.MONTH) + 1; // month starts from "0" zero !!
|
446 |
double Day = date.get(Calendar.DAY_OF_MONTH);
|
447 |
double Hour = date.get(Calendar.HOUR_OF_DAY);
|
448 |
double Minute = date.get(Calendar.MINUTE);
|
449 |
double Second = date.get(Calendar.SECOND);
|
450 |
double d = 367 * Year - div((7 * (Year + (div((Month + 9), 12)))), 4) + div((275 * Month), 9) + Day - 730530;
|
451 |
//System.out.println("dd=" + d);
|
452 |
d = d + Hour / 24 + Minute / (60 * 24) + Second / (24 * 60 * 60); // OK
|
453 |
|
454 |
// System.out.println("y=" + Year);
|
455 |
// System.out.println("m=" + Month);
|
456 |
// System.out.println("d=" + Day);
|
457 |
// System.out.println("h=" + Hour);
|
458 |
// System.out.println("m=" + Minute);
|
459 |
// System.out.println("s=" + Second);
|
460 |
// System.out.println("dd=" + d);
|
461 |
// System.out.println("c=" + date.getTimeInMillis());
|
462 |
|
463 |
double N = 125.1228 - 0.0529538083 * d;
|
464 |
double i = 5.1454;
|
465 |
|
466 |
double w = 318.0634 + 0.1643573223 * d; //OK
|
467 |
double a = 60.2666;
|
468 |
//a=6.6107940559473451507806351067866;
|
469 |
//a=0;
|
470 |
|
471 |
//a=149476000; //km average distance
|
472 |
double e = 0.054900;
|
473 |
double M = 115.3654 + 13.0649929509 * d;
|
474 |
|
475 |
w = Rev(w);
|
476 |
M = Rev(M);
|
477 |
N = Rev(N);
|
478 |
|
479 |
double E = M + (180 / Math.PI) * e * Math.sin(Math.toRadians(M)) * (1 + e * Math.cos(Math.toRadians(M)));
|
480 |
E = Rev(E); // OK
|
481 |
|
482 |
double Ebeforeit = E;
|
483 |
// now iterate until difference between E0 and E1 is less than 0.005_deg
|
484 |
// use E0, calculate E1
|
485 |
|
486 |
int Iterations = 0;
|
487 |
double E_error = 9;
|
488 |
double E0;
|
489 |
double E1;
|
490 |
//double Eafterit;
|
491 |
//double E_ErrorBefore;
|
492 |
|
493 |
while ((E_error > 0.0005) && (Iterations < 20)) // ok - itererer korrekt
|
494 |
{
|
495 |
Iterations = Iterations + 1;
|
496 |
E0 = E;
|
497 |
E1 = E0 - (E0 - (180 / Math.PI) * e * Math.sin(Math.toRadians(E0)) - M) / (1 - e * Math.cos(Math.toRadians(E0)));
|
498 |
//alert('1 E0='+E0+'\nNew E1='+E1+'\nE='+E+'\Diff='+Rev(E0-E1));
|
499 |
E = Rev(E1);
|
500 |
//alert(Math.abs(E-E0));
|
501 |
|
502 |
//Eafterit = E;
|
503 |
|
504 |
if (E < E0)
|
505 |
{
|
506 |
E_error = E0 - E;
|
507 |
}
|
508 |
else
|
509 |
{
|
510 |
E_error = E - E0;
|
511 |
}
|
512 |
|
513 |
// if (E < Ebeforeit)
|
514 |
// {
|
515 |
// //E_ErrorBefore = Ebeforeit - E;
|
516 |
// }
|
517 |
//
|
518 |
// else
|
519 |
// {
|
520 |
// //E_ErrorBefore = E - Ebeforeit;
|
521 |
// }
|
522 |
|
523 |
//System.out.println("(loop) E=" + E);
|
524 |
}
|
525 |
|
526 |
double x = a * (Math.cos(Math.toRadians(E)) - e);
|
527 |
double y = a * Math.sin(Math.toRadians(Rev(E))) * Math.sqrt(1 - e * e);
|
528 |
double r = Math.sqrt(x * x + y * y);
|
529 |
double v = Math.toDegrees(Math.atan2(y, x));
|
530 |
|
531 |
//System.out.println("E=" + E);
|
532 |
|
533 |
x = a * (Math.cos(Math.toRadians(E)) - e);
|
534 |
y = a * Math.sin(Math.toRadians(Rev(E))) * Math.sqrt(1 - e * e);
|
535 |
r = Math.sqrt(x * x + y * y);
|
536 |
v = Math.toDegrees(Math.atan2(y, x));
|
537 |
|
538 |
//alert('E='+E);
|
539 |
|
540 |
// ok s� langt
|
541 |
|
542 |
double sunlon = Rev(v + w); // trolig ok
|
543 |
|
544 |
x = r * Math.cos(Math.toRadians(sunlon));
|
545 |
y = r * Math.sin(Math.toRadians(sunlon));
|
546 |
double z = 0;
|
547 |
|
548 |
double xeclip = r * (Math.cos(Math.toRadians(N)) * Math.cos(Math.toRadians(v + w)) - Math.sin(Math.toRadians(N)) * Math.sin(Math.toRadians(v + w)) * Math.cos(Math.toRadians(i)));
|
549 |
double yeclip = r * (Math.sin(Math.toRadians(N)) * Math.cos(Math.toRadians(v + w)) + Math.cos(Math.toRadians(N)) * Math.sin(Math.toRadians(v + w)) * Math.cos(Math.toRadians(i)));
|
550 |
double zeclip = r * Math.sin(Math.toRadians(v + w)) * Math.sin(Math.toRadians(i));
|
551 |
|
552 |
double moon_longitude = Rev(Math.toDegrees(Math.atan2(yeclip, xeclip))); // OK
|
553 |
double moon_latitude = Math.toDegrees(Math.atan2(zeclip, Math.sqrt(xeclip * xeclip + yeclip * yeclip))); // trolig OK
|
554 |
|
555 |
// ----------- SUN -----------
|
556 |
// ----------- SUN -----------
|
557 |
// date.setTimeZone(this.timeZone);
|
558 |
BigDecimal longitudeHour = getLongitudeHour(date, true);
|
559 |
BigDecimal meanAnomaly = getMeanAnomaly(longitudeHour);
|
560 |
//BigDecimal sunTrueLong = getSunTrueLongitude(meanAnomaly);
|
561 |
//BigDecimal cosineSunLocalHour = getCosineSunLocalHour(sunTrueLong, Zenith.OFFICIAL);
|
562 |
|
563 |
// System.out.println("Sun MA=" + meanAnomaly);
|
564 |
|
565 |
// gesch�tzt!!!!!
|
566 |
// see -> http://www.obliquity.com/info/meaning.html
|
567 |
double sun_Obliquity = 23.45;
|
568 |
// gesch�tzt!!!!!
|
569 |
|
570 |
// sunangles[11] ????
|
571 |
double w_S = 282.9404 + 4.70935E-5 * d; //OK
|
572 |
//double a_S = 1;
|
573 |
//a=6.6107940559473451507806351067866;
|
574 |
//a=0;
|
575 |
|
576 |
//a=149476000; //km average distance
|
577 |
//double e_S = 0.016709 - 1.151E-9 * d;
|
578 |
double M_S = 356.0470 + 0.9856002585 * d;
|
579 |
double oblecl_S = 23.4393 - 3.563E-7 * d;
|
580 |
double L_S = w_S + Rev(M_S);
|
581 |
double GMST0_sun = (L_S + 180);
|
582 |
|
583 |
L_S = Rev(L_S);
|
584 |
sun_Obliquity = oblecl_S;
|
585 |
|
586 |
// System.out.println("GMST0_sun=" + GMST0_sun);
|
587 |
// System.out.println("L_S=" + L_S);
|
588 |
// System.out.println("oblecl_S=" + oblecl_S);
|
589 |
// sunangles[11] ????
|
590 |
// ----------- SUN -----------
|
591 |
// ----------- SUN -----------
|
592 |
|
593 |
double Mm = Rev(M); // Moons mean anomaly
|
594 |
double Lm = Rev(N + w + M); // moon mean longitude
|
595 |
double Ms = meanAnomaly.doubleValue(); // sun mean anomaly
|
596 |
//double Ls = sunTrueLong.doubleValue(); // sun mean longtitude
|
597 |
double Ls = L_S;
|
598 |
double D = Rev(Lm - Ls); //Moons mean elongation
|
599 |
double F = Rev(Lm - N); //Moons argument of latitude
|
600 |
|
601 |
// Perbutations Moons Longitude
|
602 |
|
603 |
double P_lon1 = -1.274 * Math.sin(Radians(Mm - 2 * D)); // (Evection)
|
604 |
double P_lon2 = +0.658 * Math.sin(Radians(2 * D)); // (Variation)
|
605 |
double P_lon3 = -0.186 * Math.sin(Radians(Ms)); // (Yearly equation)
|
606 |
double P_lon4 = -0.059 * Math.sin(Radians(2 * Mm - 2 * D));
|
607 |
double P_lon5 = -0.057 * Math.sin(Radians(Mm - 2 * D + Ms));
|
608 |
double P_lon6 = +0.053 * Math.sin(Radians(Mm + 2 * D));
|
609 |
double P_lon7 = +0.046 * Math.sin(Radians(2 * D - Ms));
|
610 |
double P_lon8 = +0.041 * Math.sin(Radians(Mm - Ms));
|
611 |
double P_lon9 = -0.035 * Math.sin(Radians(D)); // (Parallactic equation)
|
612 |
double P_lon10 = -0.031 * Math.sin(Radians(Mm + Ms));
|
613 |
double P_lon11 = -0.015 * Math.sin(Radians(2 * F - 2 * D));
|
614 |
double P_lon12 = +0.011 * Math.sin(Radians(Mm - 4 * D));
|
615 |
// Perbutations Moons Latitude
|
616 |
|
617 |
double P_lat1 = -0.173 * Math.sin(Radians(F - 2 * D));
|
618 |
double P_lat2 = -0.055 * Math.sin(Radians(Mm - F - 2 * D));
|
619 |
double P_lat3 = -0.046 * Math.sin(Radians(Mm + F - 2 * D));
|
620 |
double P_lat4 = +0.033 * Math.sin(Radians(F + 2 * D));
|
621 |
double P_lat5 = +0.017 * Math.sin(Radians(2 * Mm + F));
|
622 |
|
623 |
double P_lon = P_lon1 + P_lon2 + P_lon3 + P_lon4 + P_lon5 + P_lon6 + P_lon7 + P_lon8 + P_lon9 + P_lon10 + P_lon11 + P_lon12;
|
624 |
double P_lat = P_lat1 + P_lat2 + P_lat3 + P_lat4 + P_lat5;
|
625 |
double P_moondistance = -0.58 * Math.cos(Radians(Mm - 2 * D)) - 0.46 * Math.cos(Radians(2 * D));
|
626 |
|
627 |
//alert('P_lon='+P_lon+'\nP_lat='+P_lat+'\nP_moondistance='+P_moondistance);
|
628 |
|
629 |
moon_longitude = moon_longitude + P_lon;
|
630 |
moon_latitude = moon_latitude + P_lat;
|
631 |
r = r + P_moondistance;
|
632 |
|
633 |
// OK so far
|
634 |
// now calculate RA & Decl
|
635 |
// get the Eliptic coordinates
|
636 |
|
637 |
double xh = r * Math.cos(Radians(moon_longitude)) * Math.cos(Radians(moon_latitude));
|
638 |
double yh = r * Math.sin(Radians(moon_longitude)) * Math.cos(Radians(moon_latitude));
|
639 |
double zh = r * Math.sin(Radians(moon_latitude));
|
640 |
// rotate to rectangular equatorial coordinates
|
641 |
double xequat = xh;
|
642 |
|
643 |
double yequat = yh * Math.cos(Radians(sun_Obliquity)) - zh * Math.sin(Radians(sun_Obliquity));
|
644 |
double zequat = yh * Math.sin(Radians(sun_Obliquity)) + zh * Math.cos(Radians(sun_Obliquity));
|
645 |
double Moon_RA = Rev(Deg(Math.atan2(yh, xh))); // OK
|
646 |
double Moon_Decl = Deg(Math.atan2(zh, Math.sqrt(xh * xh + yh * yh))); // trolig OK
|
647 |
|
648 |
Moon_RA = Rev(Deg(Math.atan2(yequat, xequat))); // OK
|
649 |
Moon_Decl = Deg(Math.atan2(zequat, Math.sqrt(xequat * xequat + yequat * yequat))); // trolig OK
|
650 |
|
651 |
// System.out.println("Moon Ra=" + Moon_RA);
|
652 |
|
653 |
// System.out.println("Ls=" + Ls);
|
654 |
// war "+180" mit "-180" funkts aber besser :-)
|
655 |
double GMST0 = (Ls - 180);
|
656 |
// System.out.println("GMST0=" + GMST0);
|
657 |
|
658 |
//*********CALCULATE TIME *********************
|
659 |
|
660 |
// System.out.println("d1=" + d);
|
661 |
double UT = d - Math.floor(d);
|
662 |
//UT = 0.9;
|
663 |
// System.out.println("d1=" + UT);
|
664 |
//alert("UT="+UT);
|
665 |
|
666 |
// ???????????????????
|
667 |
// ???????????????????
|
668 |
double SiteLon = this.location.getLatitude().doubleValue();
|
669 |
double SiteLat = this.location.getLongitude().doubleValue();
|
670 |
// ???????????????????
|
671 |
// ???????????????????
|
672 |
|
673 |
double SIDEREALTIME = GMST0 + UT * 360 + SiteLon; // ok
|
674 |
double HourAngle = SIDEREALTIME - Moon_RA; // trolig ok
|
675 |
// System.out.println("GMST0 + UT * 360 + SiteLon=" + GMST0 + " " + UT + " " + SiteLon);
|
676 |
// System.out.println("SIDEREALTIME - Moon_RA=" + SIDEREALTIME + " " + Moon_RA);
|
677 |
// System.out.println("SIDEREALTIME=" + SIDEREALTIME);
|
678 |
// System.out.println("HourAngle=" + HourAngle);
|
679 |
|
680 |
// make things easier!!
|
681 |
double pi = Math.PI;
|
682 |
|
683 |
x = Math.cos(HourAngle * Math.PI / 180) * Math.cos(Moon_Decl * Math.PI / 180);
|
684 |
y = Math.sin(HourAngle * Math.PI / 180) * Math.cos(Moon_Decl * Math.PI / 180);
|
685 |
z = Math.sin(Moon_Decl * Math.PI / 180);
|
686 |
|
687 |
double xhor = x * Math.sin(SiteLat * pi / 180) - z * Math.cos(SiteLat * pi / 180);
|
688 |
//alert('sitelat='+SiteLat+'\nsitelon='+SiteLon);
|
689 |
double yhor = y;
|
690 |
double zhor = x * Math.cos(SiteLat * pi / 180) + z * Math.sin(SiteLat * pi / 180);
|
691 |
|
692 |
double MoonElevation = Deg(Math.asin(zhor)); // ok regner ikke m�ne elevation helt riktig...
|
693 |
//System.out.println("MoonElevation=" + MoonElevation);
|
694 |
|
695 |
MoonElevation = MoonElevation - Deg(Math.asin(1 / r * Math.cos(Radians(MoonElevation))));
|
696 |
//System.out.println("MoonElevation=" + MoonElevation);
|
697 |
|
698 |
double GeometricElevation = MoonElevation;
|
699 |
MoonElevation = ElevationRefraction(MoonElevation); // atmospheric refraction
|
700 |
|
701 |
double MoonAzimuth = Deg(Math.atan2(yhor, xhor));
|
702 |
|
703 |
// System.out.println("MoonElevation=" + MoonElevation);
|
704 |
// System.out.println("MoonAzimuth=" + MoonAzimuth);
|
705 |
// System.out.println("GeometricElevation=" + GeometricElevation);
|
706 |
// System.out.println("Moon_Decl=" + Moon_Decl);
|
707 |
// System.out.println("moon_longitude=" + moon_longitude);
|
708 |
// System.out.println("moon_latitude=" + moon_latitude);
|
709 |
// System.out.println("P_moondistance" + P_moondistance);
|
710 |
// System.out.println("r=" + r);
|
711 |
}
|
712 |
|
713 |
/**
|
714 |
* Computes the sunrise time for the given zenith at the given date.
|
715 |
*
|
716 |
* @param solarZenith
|
717 |
* <code>Zenith</code> enum corresponding to the type of sunrise to compute.
|
718 |
* @param date
|
719 |
* <code>Calendar</code> object representing the date to compute the sunrise for.
|
720 |
* @return the sunrise time, in HH:MM format (24-hour clock), 00:00 if the sun does not rise on the given
|
721 |
* date.
|
722 |
*/
|
723 |
public String computeSunriseTime(Zenith solarZenith, Calendar date)
|
724 |
{
|
725 |
return computeSolarEventTime(solarZenith, date, true);
|
726 |
}
|
727 |
|
728 |
/**
|
729 |
* Computes the sunset time for the given zenith at the given date.
|
730 |
*
|
731 |
* @param solarZenith
|
732 |
* <code>Zenith</code> enum corresponding to the type of sunset to compute.
|
733 |
* @param date
|
734 |
* <code>Calendar</code> object representing the date to compute the sunset for.
|
735 |
* @return the sunset time, in HH:MM format (24-hour clock), 00:00 if the sun does not set on the given
|
736 |
* date.
|
737 |
*/
|
738 |
public String computeSunsetTime(Zenith solarZenith, Calendar date)
|
739 |
{
|
740 |
return computeSolarEventTime(solarZenith, date, false);
|
741 |
}
|
742 |
|
743 |
private String computeSolarEventTime(Zenith solarZenith, Calendar date, boolean isSunrise)
|
744 |
{
|
745 |
date.setTimeZone(this.timeZone);
|
746 |
BigDecimal longitudeHour = getLongitudeHour(date, isSunrise);
|
747 |
|
748 |
BigDecimal meanAnomaly = getMeanAnomaly(longitudeHour);
|
749 |
BigDecimal sunTrueLong = getSunTrueLongitude(meanAnomaly);
|
750 |
BigDecimal cosineSunLocalHour = getCosineSunLocalHour(sunTrueLong, solarZenith);
|
751 |
if ((cosineSunLocalHour.doubleValue() < -1.0) || (cosineSunLocalHour.doubleValue() > 1.0))
|
752 |
{
|
753 |
return "99:99";
|
754 |
}
|
755 |
|
756 |
BigDecimal sunLocalHour = getSunLocalHour(cosineSunLocalHour, isSunrise);
|
757 |
BigDecimal localMeanTime = getLocalMeanTime(sunTrueLong, longitudeHour, sunLocalHour);
|
758 |
BigDecimal localTime = getLocalTime(localMeanTime, date);
|
759 |
return getLocalTimeAsString(localTime);
|
760 |
}
|
761 |
|
762 |
/**
|
763 |
* Computes the base longitude hour, lngHour in the algorithm.
|
764 |
*
|
765 |
* @return the longitude of the location of the solar event divided by 15 (deg/hour), in <code>BigDecimal</code> form.
|
766 |
*/
|
767 |
private BigDecimal getBaseLongitudeHour()
|
768 |
{
|
769 |
return divideBy(location.getLongitude(), BigDecimal.valueOf(15));
|
770 |
}
|
771 |
|
772 |
/**
|
773 |
* Computes the longitude time, t in the algorithm.
|
774 |
*
|
775 |
* @return longitudinal time in <code>BigDecimal</code> form.
|
776 |
*/
|
777 |
private BigDecimal getLongitudeHour(Calendar date, Boolean isSunrise)
|
778 |
{
|
779 |
int offset = 18;
|
780 |
if (isSunrise)
|
781 |
{
|
782 |
offset = 6;
|
783 |
}
|
784 |
BigDecimal dividend = BigDecimal.valueOf(offset).subtract(getBaseLongitudeHour());
|
785 |
BigDecimal addend = divideBy(dividend, BigDecimal.valueOf(24));
|
786 |
BigDecimal longHour = getDayOfYear(date).add(addend);
|
787 |
return setScale(longHour);
|
788 |
}
|
789 |
|
790 |
/**
|
791 |
* Computes the mean anomaly of the Sun, M in the algorithm.
|
792 |
*
|
793 |
* @return the suns mean anomaly, M, in <code>BigDecimal</code> form.
|
794 |
*/
|
795 |
private BigDecimal getMeanAnomaly(BigDecimal longitudeHour)
|
796 |
{
|
797 |
BigDecimal meanAnomaly = multiplyBy(new BigDecimal("0.9856"), longitudeHour).subtract(new BigDecimal("3.289"));
|
798 |
return setScale(meanAnomaly);
|
799 |
}
|
800 |
|
801 |
/**
|
802 |
* Computes the true longitude of the sun, L in the algorithm, at the given location, adjusted to fit in
|
803 |
* the range [0-360].
|
804 |
*
|
805 |
* @param meanAnomaly
|
806 |
* the suns mean anomaly.
|
807 |
* @return the suns true longitude, in <code>BigDecimal</code> form.
|
808 |
*/
|
809 |
private BigDecimal getSunTrueLongitude(BigDecimal meanAnomaly)
|
810 |
{
|
811 |
BigDecimal sinMeanAnomaly = new BigDecimal(Math.sin(convertDegreesToRadians(meanAnomaly).doubleValue()));
|
812 |
BigDecimal sinDoubleMeanAnomaly = new BigDecimal(Math.sin(multiplyBy(convertDegreesToRadians(meanAnomaly), BigDecimal.valueOf(2)).doubleValue()));
|
813 |
|
814 |
BigDecimal firstPart = meanAnomaly.add(multiplyBy(sinMeanAnomaly, new BigDecimal("1.916")));
|
815 |
BigDecimal secondPart = multiplyBy(sinDoubleMeanAnomaly, new BigDecimal("0.020")).add(new BigDecimal("282.634"));
|
816 |
BigDecimal trueLongitude = firstPart.add(secondPart);
|
817 |
|
818 |
if (trueLongitude.doubleValue() > 360)
|
819 |
{
|
820 |
trueLongitude = trueLongitude.subtract(BigDecimal.valueOf(360));
|
821 |
}
|
822 |
return setScale(trueLongitude);
|
823 |
}
|
824 |
|
825 |
/**
|
826 |
* Computes the suns right ascension, RA in the algorithm, adjusting for the quadrant of L and turning it
|
827 |
* into degree-hours. Will be in the range [0,360].
|
828 |
*
|
829 |
* @param sunTrueLong
|
830 |
* Suns true longitude, in <code>BigDecimal</code>
|
831 |
* @return suns right ascension in degree-hours, in <code>BigDecimal</code> form.
|
832 |
*/
|
833 |
private BigDecimal getRightAscension(BigDecimal sunTrueLong)
|
834 |
{
|
835 |
BigDecimal tanL = new BigDecimal(Math.tan(convertDegreesToRadians(sunTrueLong).doubleValue()));
|
836 |
|
837 |
BigDecimal innerParens = multiplyBy(convertRadiansToDegrees(tanL), new BigDecimal("0.91764"));
|
838 |
BigDecimal rightAscension = new BigDecimal(Math.atan(convertDegreesToRadians(innerParens).doubleValue()));
|
839 |
rightAscension = setScale(convertRadiansToDegrees(rightAscension));
|
840 |
|
841 |
if (rightAscension.doubleValue() < 0)
|
842 |
{
|
843 |
rightAscension = rightAscension.add(BigDecimal.valueOf(360));
|
844 |
}
|
845 |
else if (rightAscension.doubleValue() > 360)
|
846 |
{
|
847 |
rightAscension = rightAscension.subtract(BigDecimal.valueOf(360));
|
848 |
}
|
849 |
|
850 |
BigDecimal ninety = BigDecimal.valueOf(90);
|
851 |
BigDecimal longitudeQuadrant = sunTrueLong.divide(ninety, 0, RoundingMode.FLOOR);
|
852 |
longitudeQuadrant = longitudeQuadrant.multiply(ninety);
|
853 |
|
854 |
BigDecimal rightAscensionQuadrant = rightAscension.divide(ninety, 0, RoundingMode.FLOOR);
|
855 |
rightAscensionQuadrant = rightAscensionQuadrant.multiply(ninety);
|
856 |
|
857 |
BigDecimal augend = longitudeQuadrant.subtract(rightAscensionQuadrant);
|
858 |
return divideBy(rightAscension.add(augend), BigDecimal.valueOf(15));
|
859 |
}
|
860 |
|
861 |
private BigDecimal getCosineSunLocalHour(BigDecimal sunTrueLong, Zenith zenith)
|
862 |
{
|
863 |
BigDecimal sinSunDeclination = getSinOfSunDeclination(sunTrueLong);
|
864 |
BigDecimal cosineSunDeclination = getCosineOfSunDeclination(sinSunDeclination);
|
865 |
|
866 |
BigDecimal zenithInRads = convertDegreesToRadians(zenith.degrees());
|
867 |
BigDecimal cosineZenith = BigDecimal.valueOf(Math.cos(zenithInRads.doubleValue()));
|
868 |
BigDecimal sinLatitude = BigDecimal.valueOf(Math.sin(convertDegreesToRadians(location.getLatitude()).doubleValue()));
|
869 |
BigDecimal cosLatitude = BigDecimal.valueOf(Math.cos(convertDegreesToRadians(location.getLatitude()).doubleValue()));
|
870 |
|
871 |
BigDecimal sinDeclinationTimesSinLat = sinSunDeclination.multiply(sinLatitude);
|
872 |
BigDecimal dividend = cosineZenith.subtract(sinDeclinationTimesSinLat);
|
873 |
BigDecimal divisor = cosineSunDeclination.multiply(cosLatitude);
|
874 |
|
875 |
return setScale(divideBy(dividend, divisor));
|
876 |
}
|
877 |
|
878 |
private BigDecimal getSinOfSunDeclination(BigDecimal sunTrueLong)
|
879 |
{
|
880 |
BigDecimal sinTrueLongitude = BigDecimal.valueOf(Math.sin(convertDegreesToRadians(sunTrueLong).doubleValue()));
|
881 |
BigDecimal sinOfDeclination = sinTrueLongitude.multiply(new BigDecimal("0.39782"));
|
882 |
return setScale(sinOfDeclination);
|
883 |
}
|
884 |
|
885 |
private BigDecimal getCosineOfSunDeclination(BigDecimal sinSunDeclination)
|
886 |
{
|
887 |
BigDecimal arcSinOfSinDeclination = BigDecimal.valueOf(Math.asin(sinSunDeclination.doubleValue()));
|
888 |
BigDecimal cosDeclination = BigDecimal.valueOf(Math.cos(arcSinOfSinDeclination.doubleValue()));
|
889 |
return setScale(cosDeclination);
|
890 |
}
|
891 |
|
892 |
private BigDecimal getSunLocalHour(BigDecimal cosineSunLocalHour, Boolean isSunrise)
|
893 |
{
|
894 |
BigDecimal arcCosineOfCosineHourAngle = getArcCosineFor(cosineSunLocalHour);
|
895 |
BigDecimal localHour = convertRadiansToDegrees(arcCosineOfCosineHourAngle);
|
896 |
if (isSunrise)
|
897 |
{
|
898 |
localHour = BigDecimal.valueOf(360).subtract(localHour);
|
899 |
}
|
900 |
return divideBy(localHour, BigDecimal.valueOf(15));
|
901 |
}
|
902 |
|
903 |
private BigDecimal getLocalMeanTime(BigDecimal sunTrueLong, BigDecimal longitudeHour, BigDecimal sunLocalHour)
|
904 |
{
|
905 |
BigDecimal rightAscension = this.getRightAscension(sunTrueLong);
|
906 |
BigDecimal innerParens = longitudeHour.multiply(new BigDecimal("0.06571"));
|
907 |
BigDecimal localMeanTime = sunLocalHour.add(rightAscension).subtract(innerParens);
|
908 |
localMeanTime = localMeanTime.subtract(new BigDecimal("6.622"));
|
909 |
|
910 |
if (localMeanTime.doubleValue() < 0)
|
911 |
{
|
912 |
localMeanTime = localMeanTime.add(BigDecimal.valueOf(24));
|
913 |
}
|
914 |
else if (localMeanTime.doubleValue() > 24)
|
915 |
{
|
916 |
localMeanTime = localMeanTime.subtract(BigDecimal.valueOf(24));
|
917 |
}
|
918 |
return setScale(localMeanTime);
|
919 |
}
|
920 |
|
921 |
private BigDecimal getLocalTime(BigDecimal localMeanTime, Calendar date)
|
922 |
{
|
923 |
BigDecimal utcTime = localMeanTime.subtract(getBaseLongitudeHour());
|
924 |
BigDecimal utcOffSet = getUTCOffSet(date);
|
925 |
BigDecimal utcOffSetTime = utcTime.add(utcOffSet);
|
926 |
return adjustForDST(utcOffSetTime, date);
|
927 |
}
|
928 |
|
929 |
private BigDecimal adjustForDST(BigDecimal localMeanTime, Calendar date)
|
930 |
{
|
931 |
BigDecimal localTime = localMeanTime;
|
932 |
if (timeZone.inDaylightTime(date.getTime()))
|
933 |
{
|
934 |
localTime = localTime.add(BigDecimal.ONE);
|
935 |
}
|
936 |
if (localTime.doubleValue() > 24.0)
|
937 |
{
|
938 |
localTime = localTime.subtract(BigDecimal.valueOf(24));
|
939 |
}
|
940 |
return localTime;
|
941 |
}
|
942 |
|
943 |
/**
|
944 |
* Returns the local rise/set time in the form HH:MM.
|
945 |
*
|
946 |
* @param localTime
|
947 |
* <code>BigDecimal</code> representation of the local rise/set time.
|
948 |
* @return <code>String</code> representation of the local rise/set time in HH:MM format.
|
949 |
*/
|
950 |
private String getLocalTimeAsString(BigDecimal localTime)
|
951 |
{
|
952 |
String[] timeComponents = localTime.toPlainString().split("\\.");
|
953 |
int hour = Integer.parseInt(timeComponents[0]);
|
954 |
|
955 |
BigDecimal minutes = new BigDecimal("0." + timeComponents[1]);
|
956 |
minutes = minutes.multiply(BigDecimal.valueOf(60)).setScale(0, RoundingMode.HALF_EVEN);
|
957 |
if (minutes.intValue() == 60)
|
958 |
{
|
959 |
minutes = BigDecimal.ZERO;
|
960 |
hour += 1;
|
961 |
}
|
962 |
|
963 |
String minuteString = minutes.intValue() < 10 ? "0" + minutes.toPlainString() : minutes.toPlainString();
|
964 |
String hourString = (hour < 10) ? "0" + String.valueOf(hour) : String.valueOf(hour);
|
965 |
return hourString + ":" + minuteString;
|
966 |
}
|
967 |
|
968 |
/** ******* UTILITY METHODS (Should probably go somewhere else. ***************** */
|
969 |
|
970 |
private BigDecimal getDayOfYear(Calendar date)
|
971 |
{
|
972 |
return new BigDecimal(date.get(Calendar.DAY_OF_YEAR));
|
973 |
}
|
974 |
|
975 |
private BigDecimal getUTCOffSet(Calendar date)
|
976 |
{
|
977 |
int offSetInMillis = date.get(Calendar.ZONE_OFFSET);
|
978 |
BigDecimal offSet = new BigDecimal(offSetInMillis / 3600000);
|
979 |
return offSet.setScale(0, RoundingMode.HALF_EVEN);
|
980 |
}
|
981 |
|
982 |
private BigDecimal getArcCosineFor(BigDecimal radians)
|
983 |
{
|
984 |
BigDecimal arcCosine = BigDecimal.valueOf(Math.acos(radians.doubleValue()));
|
985 |
return setScale(arcCosine);
|
986 |
}
|
987 |
|
988 |
private BigDecimal convertRadiansToDegrees(BigDecimal radians)
|
989 |
{
|
990 |
return multiplyBy(radians, new BigDecimal(180 / Math.PI));
|
991 |
}
|
992 |
|
993 |
private BigDecimal convertDegreesToRadians(BigDecimal degrees)
|
994 |
{
|
995 |
return multiplyBy(degrees, BigDecimal.valueOf(Math.PI / 180.0));
|
996 |
}
|
997 |
|
998 |
private BigDecimal multiplyBy(BigDecimal multiplicand, BigDecimal multiplier)
|
999 |
{
|
1000 |
return setScale(multiplicand.multiply(multiplier));
|
1001 |
}
|
1002 |
|
1003 |
private BigDecimal divideBy(BigDecimal dividend, BigDecimal divisor)
|
1004 |
{
|
1005 |
return dividend.divide(divisor, 4, RoundingMode.HALF_EVEN);
|
1006 |
}
|
1007 |
|
1008 |
private BigDecimal setScale(BigDecimal number)
|
1009 |
{
|
1010 |
return number.setScale(4, RoundingMode.HALF_EVEN);
|
1011 |
}
|
1012 |
}
|