Discontinuity just after sunset

[martinmCGG]

PySolar seems to produce discontinuous sun elevation/altitude around sunset (and also near sunrise): there is a ~0.6deg jump in the following test case. May be related to #115.

Minimal example showing the issue:

import pysolar.solar as solar
import datetime
import matplotlib.pyplot as plt
import numpy as np

lat, long = 37.249870621332086, -115.81458511013524 # a randomly picked position
strange_timestamp = 1565491275.92175
epsilon = 420 #3600*12 #0.00008 # <- scale
x = np.linspace(strange_timestamp - epsilon, strange_timestamp + epsilon, 100)
y = [solar.get_altitude(lat, long, datetime.datetime.fromtimestamp(t).astimezone(datetime.timezone.utc)) for t in x]
plt.scatter(x, y, label='pysolar')

# optional, for comparison (pip install pysolar suncalc)
import suncalc
y2 = [np.rad2deg(suncalc.get_position(datetime.datetime.fromtimestamp(t).astimezone(datetime.timezone.utc), long, lat)['altitude']) for t in x]
plt.scatter(x, y2, label='suncalc')
plt.legend()
plt.show()

[pingswept]

Interesting. That certainly looks like a bug. Thanks for the test case.

My guess is that some correction factor is blowing up as altitude nears zero.

I’d welcome more details from anyone who wanted to dig into the details.

[st-bender]

Hi,
I believe this is caused by hard-limiting the refraction correction with where in

pysolar/pysolar/solar.py

Lines 338 to 342 in eabf856

	a = pressure * 2.830 * 1.02
	b = 1010.0 * temperature * 60.0 * math.tan(math.radians(tea + (10.3/(tea + 5.11))))
	del_e = math.where(tea >= -1.0*(sun_radius + atmos_refract),
	a / b, 0.)

I was trying to figure out the difference between pysolar and the algorithm by Roberto Grena in 2012 (https://doi.org/10.1016/j.solener.2012.01.024) as implemented in (https://github.com/david-salac/Fast-SZA-and-SAA-computation). The refraction correction is not limited there, giving a smoother transition in the problematic range (Fig. 1). (The refraction formula is the same as in pysolar, just converted to radians everywhere.)

However, as discussed in #14 and #23 this leads to other problems (Fig. 2).

The problem starts when tea approaches -5.11 degrees. Then the denominator in the tangent argument approaches zero, and the argument of tan() grows beyond limits. Being 2-pi periodic, tan() then occasionally gives very large positive or negative numbers.

For testing, I have increased in the limit in where to -5°, which has a smoother transition in the "jump" region and avoids the "blow-up" angles at the expense of a little kink at -5°, see Figs. 3 and 4.