/[aagtl_public1]/src/bpi/sdbm/illuminance/SolarPosition.java
aagtl

Contents of /src/bpi/sdbm/illuminance/SolarPosition.java

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (show annotations) (download)
Sun Aug 5 14:00:28 2012 UTC (11 years, 7 months ago) by zoffadmin
File size: 7342 byte(s)
license text correction
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 <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 1999–2015."
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 }

   
Visit the aagtl Website