@@ -309,7 +309,7 @@ def ephemeris(time, location, pressure=101325, temperature=12):
309
309
UnivDate = DayOfYear
310
310
UnivHr = DecHours
311
311
312
- Yr = Year - 1900
312
+ Yr = time_utc . year - 1900
313
313
YrBegin = 365 * Yr + np .floor ((Yr - 1 ) / 4. ) - 0.5
314
314
315
315
Ezero = YrBegin + UnivDate
@@ -320,16 +320,19 @@ def ephemeris(time, location, pressure=101325, temperature=12):
320
320
45.836 + 8640184.542 * T + 0.0929 * T ** 2 ) / 86400.
321
321
GMST0 = 360 * (GMST0 - np .floor (GMST0 ))
322
322
GMSTi = np .mod (GMST0 + 360 * (1.0027379093 * UnivHr / 24. ), 360 )
323
-
323
+
324
+ # Local apparent sidereal time
324
325
LocAST = np .mod ((360 + GMSTi - Longitude ), 360 )
326
+
325
327
EpochDate = Ezero + UnivHr / 24.
326
328
T1 = EpochDate / 36525.
329
+
327
330
ObliquityR = np .radians (
328
331
23.452294 - 0.0130125 * T1 - 1.64e-06 * T1 ** 2 + 5.03e-07 * T1 ** 3 )
329
332
MlPerigee = 281.22083 + 4.70684e-05 * EpochDate + 0.000453 * T1 ** 2 + (
330
333
3e-06 * T1 ** 3 )
331
334
MeanAnom = np .mod ((358.47583 + 0.985600267 * EpochDate - 0.00015 *
332
- T1 ** 2 - 3e-06 * T1 ** 3 ), 360 )
335
+ T1 ** 2 - 3e-06 * T1 ** 3 ), 360 )
333
336
Eccen = 0.01675104 - 4.18e-05 * T1 - 1.26e-07 * T1 ** 2
334
337
EccenAnom = MeanAnom
335
338
E = 0
@@ -340,7 +343,7 @@ def ephemeris(time, location, pressure=101325, temperature=12):
340
343
341
344
TrueAnom = (
342
345
2 * np .mod (np .degrees (np .arctan2 (((1 + Eccen ) / (1 - Eccen )) ** 0.5 *
343
- np .tan (np .radians (EccenAnom ) / float ( 2 ) ), 1 )), 360 ))
346
+ np .tan (np .radians (EccenAnom ) / 2. ), 1 )), 360 ))
344
347
EcLon = np .mod (MlPerigee + TrueAnom , 360 ) - Abber
345
348
EcLonR = np .radians (EcLon )
346
349
DecR = np .arcsin (np .sin (ObliquityR )* (np .sin (EcLonR )))
@@ -350,16 +353,18 @@ def ephemeris(time, location, pressure=101325, temperature=12):
350
353
351
354
HrAngle = LocAST - RtAscen
352
355
HrAngleR = np .radians (HrAngle )
353
-
354
356
HrAngle = HrAngle - (360 * ((abs (HrAngle ) > 180 )))
357
+
355
358
SunAz = np .degrees (np .arctan2 (- 1 * np .sin (HrAngleR ), np .cos (LatR ) *
356
359
(np .tan (DecR )) - np .sin (LatR )* (np .cos (HrAngleR ))))
357
- SunAz = SunAz + ( SunAz < 0 ) * 360
360
+ SunAz [ SunAz < 0 ] += 360
358
361
# potential error in the following line
359
362
SunEl = np .degrees (np .arcsin ((np .cos (LatR ) * (np .cos (DecR )) *
360
363
(np .cos (HrAngleR )) + np .sin (LatR )* (np .sin (DecR )))))
364
+
361
365
SolarTime = (180 + HrAngle ) / 15.
362
366
367
+ # Calculate refraction correction
363
368
# replace with conditional array assignment
364
369
Refract = []
365
370
for Elevation in SunEl :
0 commit comments