… | |
… | |
42 | */
|
42 | */
|
43 |
|
43 |
|
44 | public class SolarPosition
|
44 | public class SolarPosition
|
45 | {
|
45 | {
|
46 |
|
46 |
|
47 | private static final double dEarthMeanRadius = 6371.01; // in km
|
47 | private static final double dEarthMeanRadius = 6371.01; // in km
|
48 | private static final double dAstronomicalUnit = 149597890; // in km
|
48 | private static final double dAstronomicalUnit = 149597890; // in km
|
49 |
|
49 |
|
50 | private static final double pi = Math.PI;
|
50 | private static final double pi = Math.PI;
|
51 | private static final double twopi = (2 * pi);
|
51 | private static final double twopi = (2 * pi);
|
52 | private static final double rad = (pi / 180);
|
52 | private static final double rad = (pi / 180);
|
53 |
|
|
|
54 |
|
53 |
|
55 | /** result wrapper class */
|
54 | /** result wrapper class */
|
56 | public static class SunCoordinates
|
55 | public static class SunCoordinates
|
57 | {
|
56 | {
|
58 | /**
|
57 | /**
|
59 | * zenith angle, in degrees
|
58 | * zenith angle, in degrees
|
60 | */
|
59 | */
|
61 | public double zenithAngle;
|
60 | public double zenithAngle;
|
62 |
|
61 |
|
63 | /**
|
62 | /**
|
64 | * azimuth, in degrees, 0 - 360 degrees clockwise from north
|
63 | * azimuth, in degrees, 0 - 360 degrees clockwise from north
|
65 | */
|
64 | */
|
66 | public double azimuth;
|
65 | public double azimuth;
|
67 | }
|
66 | }
|
68 |
|
67 |
|
69 | /**
|
68 | /**
|
70 | * Calculate sun position for a given time and location.
|
69 | * Calculate sun position for a given time and location.
|
71 | *
|
70 | *
|
72 | * @param time
|
71 | * @param time
|
73 | * Note that it's unclear how well the algorithm performs before the year 1990 or after the year 2015.
|
72 | * Note that it's unclear how well the algorithm performs before the year 1990 or after the year 2015.
|
74 | * @param latitude
|
73 | * @param latitude
|
75 | * (positive east of Greenwich)
|
74 | * (positive east of Greenwich)
|
76 | * @param longitude
|
75 | * @param longitude
|
77 | * (positive north of equator)
|
76 | * (positive north of equator)
|
78 | * @return
|
77 | * @return
|
79 | */
|
78 | */
|
80 | public static SunCoordinates getSunPosition(final Date time, double latitude, double longitude)
|
79 | public static SunCoordinates getSunPosition(final Date time, double latitude, double longitude)
|
81 | {
|
80 | {
|
82 |
|
81 |
|
… | |
… | |
103 | {
|
102 | {
|
104 | long liAux1;
|
103 | long liAux1;
|
105 | long liAux2;
|
104 | long liAux2;
|
106 | double dJulianDate;
|
105 | double dJulianDate;
|
107 | // Calculate time of the day in UT decimal hours
|
106 | // Calculate time of the day in UT decimal hours
|
108 | dDecimalHours = utcTime.get(Calendar.HOUR_OF_DAY)
|
|
|
109 | + (utcTime.get(Calendar.MINUTE) + utcTime.get(Calendar.SECOND) / 60.0) / 60.0;
|
107 | dDecimalHours = utcTime.get(Calendar.HOUR_OF_DAY) + (utcTime.get(Calendar.MINUTE) + utcTime.get(Calendar.SECOND) / 60.0) / 60.0;
|
110 | // Calculate current Julian Day
|
108 | // Calculate current Julian Day
|
111 | liAux1 = (utcTime.get(Calendar.MONTH) + 1 - 14) / 12;
|
109 | liAux1 = (utcTime.get(Calendar.MONTH) + 1 - 14) / 12;
|
112 | liAux2 = (1461 * (utcTime.get(Calendar.YEAR) + 4800 + liAux1)) / 4
|
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;
|
113 | + (367 * (utcTime.get(Calendar.MONTH) + 1 - 2 - 12 * liAux1)) / 12
|
|
|
114 | - (3 * ((utcTime.get(Calendar.YEAR) + 4900 + liAux1) / 100)) / 4
|
|
|
115 | + utcTime.get(Calendar.DAY_OF_MONTH) - 32075;
|
|
|
116 | dJulianDate = (double) (liAux2) - 0.5 + dDecimalHours / 24.0;
|
111 | dJulianDate = (double) (liAux2) - 0.5 + dDecimalHours / 24.0;
|
117 | // Calculate difference between current Julian Day and JD 2451545.0
|
112 | // Calculate difference between current Julian Day and JD 2451545.0
|
118 | dElapsedJulianDays = dJulianDate - 2451545.0;
|
113 | dElapsedJulianDays = dJulianDate - 2451545.0;
|
119 | }
|
114 | }
|
120 |
|
115 |
|
… | |
… | |
131 | double dMeanAnomaly;
|
126 | double dMeanAnomaly;
|
132 | double dOmega;
|
127 | double dOmega;
|
133 | dOmega = 2.1429 - 0.0010394594 * dElapsedJulianDays;
|
128 | dOmega = 2.1429 - 0.0010394594 * dElapsedJulianDays;
|
134 | dMeanLongitude = 4.8950630 + 0.017202791698 * dElapsedJulianDays; // Radians
|
129 | dMeanLongitude = 4.8950630 + 0.017202791698 * dElapsedJulianDays; // Radians
|
135 | dMeanAnomaly = 6.2400600 + 0.0172019699 * dElapsedJulianDays;
|
130 | dMeanAnomaly = 6.2400600 + 0.0172019699 * dElapsedJulianDays;
|
136 | dEclipticLongitude = dMeanLongitude + 0.03341607 * Math.sin(dMeanAnomaly) + 0.00034894
|
131 | dEclipticLongitude = dMeanLongitude + 0.03341607 * Math.sin(dMeanAnomaly) + 0.00034894 * Math.sin(2 * dMeanAnomaly) - 0.0001134 - 0.0000203 * Math.sin(dOmega);
|
137 | * Math.sin(2 * dMeanAnomaly) - 0.0001134 - 0.0000203 * Math.sin(dOmega);
|
|
|
138 | dEclipticObliquity = 0.4090928 - 6.2140e-9 * dElapsedJulianDays + 0.0000396
|
132 | dEclipticObliquity = 0.4090928 - 6.2140e-9 * dElapsedJulianDays + 0.0000396 * Math.cos(dOmega);
|
139 | * Math.cos(dOmega);
|
|
|
140 | }
|
133 | }
|
141 |
|
134 |
|
142 | // System.err.println("ecl. longi. " + dEclipticLongitude);
|
135 | // System.err.println("ecl. longi. " + dEclipticLongitude);
|
143 | // System.err.println("ecl. obliq. " + dEclipticObliquity);
|
136 | // System.err.println("ecl. obliq. " + dEclipticObliquity);
|
144 |
|
|
|
145 |
|
137 |
|
146 | // Calculate celestial coordinates ( right ascension and declination )
|
138 | // Calculate celestial coordinates ( right ascension and declination )
|
147 | // in radians
|
139 | // in radians
|
148 | // but without limiting the angle to be less than 2*Pi (i.e., the result
|
140 | // but without limiting the angle to be less than 2*Pi (i.e., the result
|
149 | // may be
|
141 | // may be
|
… | |
… | |
169 | double dHourAngle;
|
161 | double dHourAngle;
|
170 | double dCos_Latitude;
|
162 | double dCos_Latitude;
|
171 | double dSin_Latitude;
|
163 | double dSin_Latitude;
|
172 | double dCos_HourAngle;
|
164 | double dCos_HourAngle;
|
173 | double dParallax;
|
165 | double dParallax;
|
174 | dGreenwichMeanSiderealTime = 6.6974243242 + 0.0657098283 * dElapsedJulianDays
|
166 | dGreenwichMeanSiderealTime = 6.6974243242 + 0.0657098283 * dElapsedJulianDays + dDecimalHours;
|
175 | + dDecimalHours;
|
|
|
176 | dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime * 15 + longitude) * rad;
|
167 | dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime * 15 + longitude) * rad;
|
177 | dHourAngle = dLocalMeanSiderealTime - dRightAscension;
|
168 | dHourAngle = dLocalMeanSiderealTime - dRightAscension;
|
178 | dLatitudeInRadians = latitude * rad;
|
169 | dLatitudeInRadians = latitude * rad;
|
179 | dCos_Latitude = Math.cos(dLatitudeInRadians);
|
170 | dCos_Latitude = Math.cos(dLatitudeInRadians);
|
180 | dSin_Latitude = Math.sin(dLatitudeInRadians);
|
171 | dSin_Latitude = Math.sin(dLatitudeInRadians);
|
181 | dCos_HourAngle = Math.cos(dHourAngle);
|
172 | dCos_HourAngle = Math.cos(dHourAngle);
|
182 | retval.zenithAngle = (Math.acos(dCos_Latitude * dCos_HourAngle * Math.cos(dDeclination)
|
173 | retval.zenithAngle = (Math.acos(dCos_Latitude * dCos_HourAngle * Math.cos(dDeclination) + Math.sin(dDeclination) * dSin_Latitude));
|
183 | + Math.sin(dDeclination) * dSin_Latitude));
|
|
|
184 | dY = -Math.sin(dHourAngle);
|
174 | dY = -Math.sin(dHourAngle);
|
185 | dX = Math.tan(dDeclination) * dCos_Latitude - dSin_Latitude * dCos_HourAngle;
|
175 | dX = Math.tan(dDeclination) * dCos_Latitude - dSin_Latitude * dCos_HourAngle;
|
186 | retval.azimuth = Math.atan2(dY, dX);
|
176 | retval.azimuth = Math.atan2(dY, dX);
|
187 | if (retval.azimuth < 0.0) retval.azimuth = retval.azimuth + twopi;
|
177 | if (retval.azimuth < 0.0) retval.azimuth = retval.azimuth + twopi;
|
188 | retval.azimuth = retval.azimuth / rad;
|
178 | retval.azimuth = retval.azimuth / rad;
|
… | |
… | |
192 | }
|
182 | }
|
193 |
|
183 |
|
194 | return retval;
|
184 | return retval;
|
195 | }
|
185 | }
|
196 |
|
186 |
|
197 |
|
|
|
198 | public static void main(String[] args)
|
187 | public static void main(String[] args)
|
199 | {
|
188 | {
|
200 |
|
189 |
|
201 | SolarPosition.getSunPosition(new Date(), 48.2, 16.4);
|
190 | SolarPosition.getSunPosition(new Date(), 48.2, 16.4);
|
202 | //System.out.println("current azimuth " + sc.azimuth);
|
191 | //System.out.println("current azimuth " + sc.azimuth);
|