1 |
/**
|
2 |
* aagtl Advanced Geocaching Tool for Android
|
3 |
* loosely based on agtl by Daniel Fett <fett@danielfett.de>
|
4 |
* Copyright (C) 2010 - 2012 Zoff.cc <aagtl@work.zoff.cc>
|
5 |
*
|
6 |
* This program is free software; you can redistribute it and/or
|
7 |
* modify it under the terms of the GNU General Public License
|
8 |
* version 2 as published by the Free Software Foundation.
|
9 |
*
|
10 |
* This program is distributed in the hope that it will be useful,
|
11 |
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
12 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
13 |
* GNU General Public License for more details.
|
14 |
*
|
15 |
* You should have received a copy of the GNU General Public License
|
16 |
* along with this program; if not, write to the
|
17 |
* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
|
18 |
* Boston, MA 02110-1301, USA.
|
19 |
*/
|
20 |
|
21 |
package bpi.sdbm.illuminance;
|
22 |
|
23 |
import java.util.Calendar;
|
24 |
import java.util.Date;
|
25 |
import java.util.GregorianCalendar;
|
26 |
import java.util.TimeZone;
|
27 |
|
28 |
/**
|
29 |
*
|
30 |
* Compute sun position for a given date/time and longitude/latitude.
|
31 |
*
|
32 |
* This is a simple Java port of the "PSA" solar positioning algorithm, as
|
33 |
* documented in:
|
34 |
*
|
35 |
* Blanco-Muriel et al.: Computing the Solar Vector. Solar Energy Vol 70 No 5
|
36 |
* pp 431-441. http://dx.doi.org/10.1016/S0038-092X(00)00156-0
|
37 |
*
|
38 |
* According to the paper, "The algorithm allows .. the true solar vector to
|
39 |
* be determined with an accuracy of 0.5 minutes of arc for the period 19992015."
|
40 |
*
|
41 |
* @author Klaus A. Brunner
|
42 |
*
|
43 |
*/
|
44 |
|
45 |
public class SolarPosition
|
46 |
{
|
47 |
|
48 |
private static final double dEarthMeanRadius = 6371.01; // in km
|
49 |
private static final double dAstronomicalUnit = 149597890; // in km
|
50 |
|
51 |
private static final double pi = Math.PI;
|
52 |
private static final double twopi = (2 * pi);
|
53 |
private static final double rad = (pi / 180);
|
54 |
|
55 |
|
56 |
/** result wrapper class */
|
57 |
public static class SunCoordinates
|
58 |
{
|
59 |
/**
|
60 |
* zenith angle, in degrees
|
61 |
*/
|
62 |
public double zenithAngle;
|
63 |
|
64 |
/**
|
65 |
* azimuth, in degrees, 0 - 360 degrees clockwise from north
|
66 |
*/
|
67 |
public double azimuth;
|
68 |
}
|
69 |
|
70 |
/**
|
71 |
* Calculate sun position for a given time and location.
|
72 |
*
|
73 |
* @param time
|
74 |
* Note that it's unclear how well the algorithm performs before the year 1990 or after the year 2015.
|
75 |
* @param latitude
|
76 |
* (positive east of Greenwich)
|
77 |
* @param longitude
|
78 |
* (positive north of equator)
|
79 |
* @return
|
80 |
*/
|
81 |
public static SunCoordinates getSunPosition(final Date time, double latitude, double longitude)
|
82 |
{
|
83 |
|
84 |
SunCoordinates retval = new SunCoordinates();
|
85 |
|
86 |
Calendar utcTime = new GregorianCalendar(TimeZone.getTimeZone("GMT"));
|
87 |
utcTime.setTimeInMillis(time.getTime());
|
88 |
|
89 |
// Main variables
|
90 |
double dElapsedJulianDays;
|
91 |
double dDecimalHours;
|
92 |
double dEclipticLongitude;
|
93 |
double dEclipticObliquity;
|
94 |
double dRightAscension;
|
95 |
double dDeclination;
|
96 |
|
97 |
// Auxiliary variables
|
98 |
double dY;
|
99 |
double dX;
|
100 |
|
101 |
// Calculate difference in days between the current Julian Day
|
102 |
// and JD 2451545.0, which is noon 1 January 2000 Universal Time
|
103 |
|
104 |
{
|
105 |
long liAux1;
|
106 |
long liAux2;
|
107 |
double dJulianDate;
|
108 |
// Calculate time of the day in UT decimal hours
|
109 |
dDecimalHours = utcTime.get(Calendar.HOUR_OF_DAY)
|
110 |
+ (utcTime.get(Calendar.MINUTE) + utcTime.get(Calendar.SECOND) / 60.0) / 60.0;
|
111 |
// Calculate current Julian Day
|
112 |
liAux1 = (utcTime.get(Calendar.MONTH) + 1 - 14) / 12;
|
113 |
liAux2 = (1461 * (utcTime.get(Calendar.YEAR) + 4800 + liAux1)) / 4
|
114 |
+ (367 * (utcTime.get(Calendar.MONTH) + 1 - 2 - 12 * liAux1)) / 12
|
115 |
- (3 * ((utcTime.get(Calendar.YEAR) + 4900 + liAux1) / 100)) / 4
|
116 |
+ utcTime.get(Calendar.DAY_OF_MONTH) - 32075;
|
117 |
dJulianDate = (double) (liAux2) - 0.5 + dDecimalHours / 24.0;
|
118 |
// Calculate difference between current Julian Day and JD 2451545.0
|
119 |
dElapsedJulianDays = dJulianDate - 2451545.0;
|
120 |
}
|
121 |
|
122 |
// System.err.println("elapsed julian " + dElapsedJulianDays);
|
123 |
// System.err.println("decimal hrs " + dDecimalHours);
|
124 |
|
125 |
// Calculate ecliptic coordinates (ecliptic longitude and obliquity of
|
126 |
// the
|
127 |
// ecliptic in radians but without limiting the angle to be less than
|
128 |
// 2*Pi
|
129 |
// (i.e., the result may be greater than 2*Pi)
|
130 |
{
|
131 |
double dMeanLongitude;
|
132 |
double dMeanAnomaly;
|
133 |
double dOmega;
|
134 |
dOmega = 2.1429 - 0.0010394594 * dElapsedJulianDays;
|
135 |
dMeanLongitude = 4.8950630 + 0.017202791698 * dElapsedJulianDays; // Radians
|
136 |
dMeanAnomaly = 6.2400600 + 0.0172019699 * dElapsedJulianDays;
|
137 |
dEclipticLongitude = dMeanLongitude + 0.03341607 * Math.sin(dMeanAnomaly) + 0.00034894
|
138 |
* Math.sin(2 * dMeanAnomaly) - 0.0001134 - 0.0000203 * Math.sin(dOmega);
|
139 |
dEclipticObliquity = 0.4090928 - 6.2140e-9 * dElapsedJulianDays + 0.0000396
|
140 |
* Math.cos(dOmega);
|
141 |
}
|
142 |
|
143 |
// System.err.println("ecl. longi. " + dEclipticLongitude);
|
144 |
// System.err.println("ecl. obliq. " + dEclipticObliquity);
|
145 |
|
146 |
|
147 |
// Calculate celestial coordinates ( right ascension and declination )
|
148 |
// in radians
|
149 |
// but without limiting the angle to be less than 2*Pi (i.e., the result
|
150 |
// may be
|
151 |
// greater than 2*Pi)
|
152 |
{
|
153 |
double dSin_EclipticLongitude;
|
154 |
dSin_EclipticLongitude = Math.sin(dEclipticLongitude);
|
155 |
dY = Math.cos(dEclipticObliquity) * dSin_EclipticLongitude;
|
156 |
dX = Math.cos(dEclipticLongitude);
|
157 |
dRightAscension = Math.atan2(dY, dX);
|
158 |
if (dRightAscension < 0.0) dRightAscension = dRightAscension + 2 * Math.PI;
|
159 |
dDeclination = Math.asin(Math.sin(dEclipticObliquity) * dSin_EclipticLongitude);
|
160 |
}
|
161 |
|
162 |
// System.err.println("right asc " + dRightAscension);
|
163 |
// System.err.println("decl. " + dDeclination);
|
164 |
|
165 |
// Calculate local coordinates ( azimuth and zenith angle ) in degrees
|
166 |
{
|
167 |
double dGreenwichMeanSiderealTime;
|
168 |
double dLocalMeanSiderealTime;
|
169 |
double dLatitudeInRadians;
|
170 |
double dHourAngle;
|
171 |
double dCos_Latitude;
|
172 |
double dSin_Latitude;
|
173 |
double dCos_HourAngle;
|
174 |
double dParallax;
|
175 |
dGreenwichMeanSiderealTime = 6.6974243242 + 0.0657098283 * dElapsedJulianDays
|
176 |
+ dDecimalHours;
|
177 |
dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime * 15 + longitude) * rad;
|
178 |
dHourAngle = dLocalMeanSiderealTime - dRightAscension;
|
179 |
dLatitudeInRadians = latitude * rad;
|
180 |
dCos_Latitude = Math.cos(dLatitudeInRadians);
|
181 |
dSin_Latitude = Math.sin(dLatitudeInRadians);
|
182 |
dCos_HourAngle = Math.cos(dHourAngle);
|
183 |
retval.zenithAngle = (Math.acos(dCos_Latitude * dCos_HourAngle * Math.cos(dDeclination)
|
184 |
+ Math.sin(dDeclination) * dSin_Latitude));
|
185 |
dY = -Math.sin(dHourAngle);
|
186 |
dX = Math.tan(dDeclination) * dCos_Latitude - dSin_Latitude * dCos_HourAngle;
|
187 |
retval.azimuth = Math.atan2(dY, dX);
|
188 |
if (retval.azimuth < 0.0) retval.azimuth = retval.azimuth + twopi;
|
189 |
retval.azimuth = retval.azimuth / rad;
|
190 |
// Parallax Correction
|
191 |
dParallax = (dEarthMeanRadius / dAstronomicalUnit) * Math.sin(retval.zenithAngle);
|
192 |
retval.zenithAngle = (retval.zenithAngle + dParallax) / rad;
|
193 |
}
|
194 |
|
195 |
return retval;
|
196 |
}
|
197 |
|
198 |
|
199 |
public static void main(String[] args)
|
200 |
{
|
201 |
|
202 |
SunCoordinates sc = SolarPosition.getSunPosition(new Date(), 48.2, 16.4);
|
203 |
|
204 |
//System.out.println("current azimuth " + sc.azimuth);
|
205 |
//System.out.println("current zenith angle " + sc.zenithAngle);
|
206 |
|
207 |
}
|
208 |
}
|