/[zanavi_public1]/navit/navit/android/src/com/luckycatlabs/sunrisesunset/calculator/SolarEventCalculator.java
ZANavi

Contents of /navit/navit/android/src/com/luckycatlabs/sunrisesunset/calculator/SolarEventCalculator.java

Parent Directory Parent Directory | Revision Log Revision Log


Revision 31 - (show annotations) (download)
Mon Feb 4 17:41:59 2013 UTC (11 years, 1 month ago) by zoff99
File size: 36223 byte(s)
new map version, lots of fixes and experimental new features
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 }

   
Visit the ZANavi Wiki