Unexpected values through a 24 hour cycle

[JamesLawrenceGSI]

Using the following code to show a 24 hour cycle at 0.1 hour intervals using the pysolar library, I get very unusual outputs which I can't explain:

import matplotlib.pyplot as plt
from pysolar.solar import *

SI=[]
date = datetime.datetime(2010, 1, 1, 0, 0, 0, 130320)
end_date = datetime.datetime(2010, 1, 1, 23, 0, 0, 130320)
long = -71.4
lat = 42.3

while (date <= end_date):
    date = date + datetime.timedelta(hours=0.1)
    altitude_deg = get_altitude(lat, long, date)
    azimuth_deg = get_azimuth(lat, long, date)
    SI.append( radiation.get_radiation_direct(date, altitude_deg) )

plt.plot(SI)
plt.show()

The plot from this is:

The spikes occur when the 'altitude' parameter tends to 0, but not necessarily when it is negative, for example, at 0.564 degrees, the solar irradiance is computed at 15953163.

If I only capture the solar irradiance when the altitude > 1 degree, I get this (as expected):

The code for this is:

import matplotlib.pyplot as plt
SI=[]
date = datetime.datetime(2010, 1, 1, 0, 0, 0, 130320)
end_date = datetime.datetime(2010, 1, 1, 23, 0, 0, 130320)
long = -71.4
lat = 42.3

while (date <= end_date):
    date = date + datetime.timedelta(hours=temporal_resolution)
    altitude_deg = get_altitude(lat, long, date)
    azimuth_deg = get_azimuth(lat, long, date)
    if altitude_deg > 1:
        SI.append( radiation.get_radiation_direct(date, altitude_deg) )

plt.plot(SI)
plt.show()

Also, shouldn't altitude be called elevation? In physics, the angle subtended from the horizon to some object is generally referred to as the elevation angle, or am I misunderstanding its meaning?

[pingswept]

This looks like a bug to me. We don't have any good validation of radiation.get_radiation_direct, but the spike you're seeing at low altitude is obviously insane. The function that is going awry is here: https://github.com/pingswept/pysolar/blob/master/pysolar/radiation.py#L43

A quick fix would be to just wipe out all radiation below 1 degree, given that it's close to zero anyway.

If anyone can figure out what's going wrong and a good way to fix it, patches are welcome.

[louis-red]

I can't reproduce this issue with my pysolar version (0.7, python 3.5.5 64-bit). As the get_radiation_direct function only takes account of the when parameter for the tm_yday (days since January 1) attribute, the function get_air_mass_ratio that uniquely uses the altitude_deg parameter should be at fault. But the only formula in there is

result = 1 / math.sin(math.radians(altitude_deg))

which tends near zero to +inf. But the result of the function is then used here :

air_mass_ratio = get_air_mass_ratio(altitude_deg)
return  flux * math.exp(-1 * optical_depth * air_mass_ratio)  

and the result tends to zero as air_mass_ratio tends to +infinity. So there is no bug. My best guess would be that you're on a 32 bits platform, and the result overflew to the negatives.

[louis-red]

@JamesLawrenceGSI would you have by any means executed the script on a 32-bit platform ?

[adbuerger]

Hello @JacquotLeHaricot I wanted to let you know that I can reproduce the error reported by @JamesLawrenceGSI and according to this post I'm running a 64 bit version of Python 3:

Python 3.5.2 (default, Nov 23 2017, 16:37:01) 
[GCC 5.4.0 20160609] on linux

>>> import struct
>>> print(struct.calcsize("P") * 8)
64

Can you think of another reason why the error is not showing for you?

[AmitAronovitch]

@adbuerger @JamesLawrenceGSI can you test with master? I think this was solved a year ago, in this commit: https://github.com/pingswept/pysolar/commit/d94fa5549e4f1435a947fdd7a97e13f374818f56

However, this was after the release of 0.7

(to test with master, clone the repo, pip uninstall pysolar and pip install -e . inside the project tree)

[adbuerger]

@AmitAronovitch thanks for pointing this out! Yes, the peaks are gone, see below the output when using the version from master:

import matplotlib.pyplot as plt
from pysolar.solar import *
import pytz

cet = pytz.timezone("Europe/Berlin")

SI=[]
date = cet.localize(datetime.datetime(2010, 1, 1, 0, 0, 0, 130320))
end_date = cet.localize(datetime.datetime(2010, 1, 1, 23, 0, 0, 130320))
long = -71.4
lat = 42.3

while (date <= end_date):
    date = date + datetime.timedelta(hours=0.1)
    altitude_deg = get_altitude(lat, long, date)
    azimuth_deg = get_azimuth(lat, long, date)
    SI.append( radiation.get_radiation_direct(date, altitude_deg) )

plt.plot(SI)
plt.show()

[pingswept]

Fix released in 0.8, so I think we're done with this.