@@ -267,6 +267,10 @@ def ephemeris(time, location, pressure=101325, temperature=12):
267
267
268
268
# Added by Rob Andrews (@Calama-Consulting), Calama Consulting, 2014
269
269
# Edited by Will Holmgren (@wholmgren), University of Arizona, 2014
270
+
271
+ # Most comments in this function are from PVLIB_MATLAB or from
272
+ # pvlib-python's attempt to understand and fix problems with the
273
+ # algorithm. The comments are *not* based on the reference material.
270
274
271
275
pvl_logger .warning ('Using solarposition.ephemeris is discouraged. '
272
276
'solarposition.pyephem and solarposition.spa are '
@@ -275,55 +279,50 @@ def ephemeris(time, location, pressure=101325, temperature=12):
275
279
pvl_logger .debug ('location={}, temperature={}, pressure={}' .format (
276
280
location , temperature , pressure ))
277
281
278
- Latitude = location .latitude
279
- ''' the inversion of longitude is due to the fact that this code was
280
- originally written for the convention that positive longitude were for
281
- locations west of the prime meridian. However, the correct convention (as
282
- of 2009) is to use negative longitudes for locations west of the prime
283
- meridian. Therefore, the user should input longitude values under the
284
- correct convention (e.g. Albuquerque is at -106 longitude), but it needs
285
- to be inverted for use in the code.
286
- '''
282
+ # the inversion of longitude is due to the fact that this code was
283
+ # originally written for the convention that positive longitude were for
284
+ # locations west of the prime meridian. However, the correct convention (as
285
+ # of 2009) is to use negative longitudes for locations west of the prime
286
+ # meridian. Therefore, the user should input longitude values under the
287
+ # correct convention (e.g. Albuquerque is at -106 longitude), but it needs
288
+ # to be inverted for use in the code.
289
+
287
290
Latitude = location .latitude
288
291
Longitude = 1 * location .longitude
289
- Year = time .year
290
- Hour = time .hour
291
- Minute = time .minute
292
- Second = time .second
293
- DayOfYear = time .dayofyear
294
-
295
- DecHours = Hour + Minute / float (60 ) + Second / float (3600 )
296
-
297
- Abber = 20 / float (3600 )
292
+
293
+ Abber = 20 / 3600.
298
294
LatR = np .radians (Latitude )
299
-
300
- # this code is needed to handle the new location.tz format string
301
- try :
302
- utc_offset = pytz .timezone (location .tz ).utcoffset (
303
- time [0 ]).total_seconds () / 3600.
304
- except ValueError :
305
- utc_offset = time .tzinfo .utcoffset (time [0 ]).total_seconds () / 3600.
306
- pvl_logger .debug ('utc_offset={}' .format (utc_offset ))
295
+
296
+ # the SPA algorithm needs time to be expressed in terms of
297
+ # decimal UTC hours of the day of the year.
298
+
299
+ # first convert to utc
300
+ time_utc = localize_to_utc (time , location )
301
+
302
+ # strip out the day of the year and calculate the decimal hour
303
+ DayOfYear = time_utc .dayofyear
304
+ DecHours = (time_utc .hour + time_utc .minute / 60. + time_utc .second / 3600.
305
+ time_utc .microsecond / 3600.e6 )
307
306
308
307
UnivDate = DayOfYear
309
-
310
- # ToDo: Will H: surprised to see the 0.5 here, but moving on...
311
- UnivHr = DecHours + utc_offset - .5
308
+ UnivHr = DecHours
312
309
313
310
Yr = Year - 1900
314
311
315
312
YrBegin = 365 * Yr + np .floor ((Yr - 1 ) / float (4 )) - 0.5
316
313
317
314
Ezero = YrBegin + UnivDate
318
- T = Ezero / float (36525 )
319
- GMST0 = 6 / float (24 ) + 38 / float (1440 ) + (
320
- 45.836 + 8640184.542 * T + 0.0929 * T ** 2 ) / float (86400 )
315
+ T = Ezero / 36525.
316
+
317
+ # Calculate Greenwich Mean Sidereal Time (GMST)
318
+ GMST0 = 6 / 24. + 38 / 1440. + (
319
+ 45.836 + 8640184.542 * T + 0.0929 * T ** 2 ) / 86400.
321
320
GMST0 = 360 * (GMST0 - np .floor (GMST0 ))
322
321
GMSTi = np .mod (GMST0 + 360 * (1.0027379093 * UnivHr / float (24 )), 360 )
323
322
324
323
LocAST = np .mod ((360 + GMSTi - Longitude ), 360 )
325
- EpochDate = Ezero + UnivHr / float ( 24 )
326
- T1 = EpochDate / float ( 36525 )
324
+ EpochDate = Ezero + UnivHr / 24.
325
+ T1 = EpochDate / 36525.
327
326
ObliquityR = np .radians (
328
327
23.452294 - 0.0130125 * T1 - 1.64e-06 * T1 ** 2 + 5.03e-07 * T1 ** 3 )
329
328
MlPerigee = 281.22083 + 4.70684e-05 * EpochDate + 0.000453 * T1 ** 2 + (
@@ -338,7 +337,6 @@ def ephemeris(time, location, pressure=101325, temperature=12):
338
337
E = EccenAnom
339
338
EccenAnom = MeanAnom + np .degrees (Eccen ) * (np .sin (np .radians (E )))
340
339
341
- # pdb.set_trace()
342
340
TrueAnom = (
343
341
2 * np .mod (np .degrees (np .arctan2 (((1 + Eccen ) / (1 - Eccen )) ** 0.5 *
344
342
np .tan (np .radians (EccenAnom ) / float (2 )), 1 )), 360 ))
@@ -361,8 +359,8 @@ def ephemeris(time, location, pressure=101325, temperature=12):
361
359
(np .cos (HrAngleR )) + np .sin (LatR )* (np .sin (DecR )))))
362
360
SolarTime = (180 + HrAngle ) / float (15 )
363
361
362
+ # replace with conditional array assignment
364
363
Refract = []
365
-
366
364
for Elevation in SunEl :
367
365
TanEl = np .tan (np .radians (Elevation ))
368
366
if Elevation > 5 and Elevation <= 85 :
@@ -378,7 +376,7 @@ def ephemeris(time, location, pressure=101325, temperature=12):
378
376
Refract .append (0 )
379
377
380
378
Refract = (np .array (Refract ) * ((283 / float (273 + temperature ))) *
381
- pressure / float ( 101325 ) / float ( 3600 ) )
379
+ pressure / 101325. / 3600. )
382
380
383
381
SunZen = 90 - SunEl
384
382
0 commit comments