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