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