Discontinuity of first derivative in `solar.get_altitude()`

[citypilgrim]

I live in Singapore, and am studying the altitude of the sun across the day, for different times of the year. I realised that the output generated by 'solar.get_azimuth()' is discontinuous in the first derivative. This is the code that I have used:

# imports
import datetime as dt

import matplotlib.pyplot as plt
import matplotlib.basic_units as pbu
import numpy as np
import pysolar.solar as pssl

# params
lt, lg = 1.299119, 103.771232
ele = 0

year = 2019
mon_start, mon_end, mon_step = 1, 12, 1 # step !< 1
day_start, day_end, day_step = 1, 1, 1
hr_start, hr_end, hr_step = 0, 23, 1
mn_start, mn_end, mn_step = 0, 59, 1

# generating data
when_ara = np.array([[year, mon_step, day_step, hr_step, mn_step]]) #1st ent tb rm ltr
alt_ara = np.array([])

for mon in range(mon_start, mon_end + 1, mon_step):
    print('retrieving for month {}...'.format(mon))
    for day in range(day_start, day_end + 1, day_step):
        for hr in range(hr_start, hr_end + 1, hr_step):
            for mn in range(mn_start, mn_end + 1, mn_step):
                when = dt.datetime(year, mon, day, hr, mn
                                   , tzinfo=dt.timezone(dt.timedelta(hours=8)))
                when_ara = np.append(when_ara, [[year, mon, day, hr, mn]], axis=0)

                # TOGGLE BETWEEN FAST AND SLOW HERE
                alt = np.deg2rad(pssl.get_altitude(lt, lg, when, elevation= ele))
                # alt = np.deg2rad(pssl.get_altitude_fast(lt, lg, when))#, elevation= ele))                
                alt_ara = np.append(alt_ara, alt)

when_ara = when_ara[1:]
when_ara = when_ara.reshape((mon_end - mon_start) // mon_step + 1
                            , (day_end - day_start) // day_step + 1
                            , (hr_end - hr_start) // hr_step + 1
                            , (mn_end - mn_start) // mn_step + 1
                            , when_ara.shape[-1])
alt_ara = alt_ara.reshape((mon_end - mon_start) // mon_step + 1
                                , (day_end - day_start) // day_step + 1
                                , (hr_end - hr_start) // hr_step + 1
                                , (mn_end - mn_start) // mn_step + 1)

# adjusting data for plot
when_ara = when_ara.reshape(np.prod(when_ara.shape[:-3])
                            , np.prod(when_ara.shape[-3:-1])
                            , when_ara.shape[-1]) # when_ara axis=-1 include yr,...
when_ara = when_ara[..., -2] + when_ara[..., -1]/60
alt_ara = alt_ara.reshape(np.prod(alt_ara.shape[:-2])
                              , np.prod(alt_ara.shape[-2:]))
# taking first derivative
der1alt_ara = np.gradient(alt_ara, when_ara[0][1] - when_ara[0][0], axis=-1)

# plot
fig = plt.figure()
ax = fig.add_subplot(111)
ax1 = ax.twinx()

for i in range(len(when_ara)):
    plot, = ax.plot(when_ara[i], alt_ara[i], yunits=pbu.degrees
                    , label='month {}'.format(i+1), marker='o')

    ax1.plot(when_ara[i], der1alt_ara[i], yunits=pbu.degrees
             , color=plot.get_color())

ax.legend(fontsize='x-small')
ax.set_xlabel('local time [24hr]')
ax.set_ylabel('alt [deg]')
ax1.set_ylabel(r'$\frac{d alt}{dt}$ [deg/hr]')

plt.show()

This is the output I got by using solar.get_altitude(), the 'o' markers represent altitude, '-' lines represent the first derivative.

this is the output I got by using solar.get_altitude_fast()

[pingswept]

Your original function looks continuous, so the problem is probably in your derivative calculation. I'm closing this; feel free to re-open if it's a persistent problem.

[citypilgrim]

The problem still persists, the function I am using to get the gradient is np.gradient, there should not be any problem with that. Also, the bump exists within the original array, it is small but it is there, I have zoomed in for you in the plot below. I have condensed the code to a single date range across one day.

# imports
import datetime as dt

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pysolar.solar as pssl

# params
lt, lg = 1.299119, 103.771232
ele = 0

# getting time array
start = pd.Timestamp('2020-02-10')
end = pd.Timestamp('2020-02-11')
when_ara = pd.date_range(start, end, freq='min', tz='Asia/Singapore')

# calculating altitude
alt_ara = np.array([
    np.deg2rad(pssl.get_altitude(lt, lg, when, elevation=ele))\
    for when in when_ara
])

# taking first derivative
der1alt_ara = np.gradient(alt_ara, when_ara[1].value - when_ara[0].value)

# plot
fig = plt.figure()
ax = fig.add_subplot(111)
ax1 = ax.twinx()

ax.plot_date(when_ara, alt_ara, '-', color='C0')
ax1.plot_date(when_ara, der1alt_ara, '-', color='C1')

ax.legend(fontsize='x-small')
ax.set_xlabel('local time [24hr]')
ax.set_ylabel('alt [deg]')
ax1.set_ylabel(r'$\frac{d alt}{dt}$ [a.u.]')

plt.show()

[pingswept]

Ah, it does appear that there is a bump in the altitude-- about a 0.03 degree change in around 1 minute, which looks like about twice the rate of change before and after the bump.

To figure out why, I would look at each of the components used to calculate the altitude, and see which one of them is changing quickly. My guess is that something is getting weird because it's in denominator and near zero.

I'm not sure what to tell you about np.gradient. In the blue line, it looks like the slope roughly doubles, rather than spiking violently.

[citypilgrim]

Thanks, yes the slope indeed does roughly double, the spike in the np.gradient plot is due to me plotting as a continuous line. matplotlib performs a linear interpolation between discrete points

[pingswept]

Closing this for now. Please reopen if I’m in error here. See also #155.