/[zanavi_public1]/navit/navit/android/src/bpi/sdbm/illuminance/SolarPosition.java
ZANavi

Contents of /navit/navit/android/src/bpi/sdbm/illuminance/SolarPosition.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: 7186 byte(s)
new map version, lots of fixes and experimental new features
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 }

   
Visit the ZANavi Wiki