Checking calculation of Radiation at Polar Region

Setup environment for tutorials

import os
import sys
import numpy as np
import math
from matplotlib import pyplot as plt
from datetime import datetime
import shyft.hydrology as api
from shyft.time_series import (Calendar,deltahours,deltaminutes,UtcPeriod,TimeSeries,TimeAxis,DoubleVector,point_interpretation_policy)
#import seaborn as sns
albedo=0.5
turbidity=0.5

radparam = api.RadiationParameter(albedo, turbidity)
radcal = api.RadiationCalculator(radparam)
radres = api.RadiationResponse()

latitude_deg = 76.8 # positive in the northern hemisphere
longitude_deg = 0 # negative reckoning west from prime meridian in Greenwich, England
elevation = 0
# ds0 = utc.time(2012, 8, 21, 0)

utc = Calendar()
ds0 = utc.time(2017, 1, 14, 0, 0, 0)

res = []
t = []
for i in range(0, 24*60*360 , 60):
    ds = ds0 + deltaminutes(i)
    t.append(int(ds))
    # adcal.net_radiation_step(radres, latitude_deg, ds, api.deltaminutes(10), 0, 0, 0, 60, elevation, 0)
    radcal.net_radiation_inst(radres, latitude_deg, ds, 0, 0, 0, 60, elevation, 0)
    res.append(radres.sw_cs_p)

fig, ax = plt.subplots(figsize=(18, 6))
fig.autofmt_xdate()
t = np.array(t).astype('datetime64[s]')
ax.plot(t, res, label='radiation @ {}'.format(latitude_deg))
ax.xaxis_date()
plt.autoscale(enable=True, axis='x', tight=True)
plt.xlabel('hour', fontsize=15)
plt.ylabel('W/m2', fontsize=15)
plt.title('Clear sky radiation on normal plane', fontsize=18)
plt.legend()
plt.grid(True)
plt.show()

ds0 = utc.time(2017, 6, 14, 0)

res = []
t = []
for i in range(0, 23):
    ds0 = utc.time(2017, 6, 14, i, 0, 0)
    ds1 = utc.time(2017, 6, 14, i+1, 0, 0)
    t.append(int(ds0))
    # adcal.net_radiation_step(radres, latitude_deg, ds, api.deltaminutes(10), 0, 0, 0, 60, elevation, 0)
    radcal.net_radiation_step(radres, latitude_deg, ds0,deltahours(1), 0, 0, 0, 60, elevation, 0)
    res.append(radres.sw_cs_p)

fig, ax = plt.subplots(figsize=(18, 6))
fig.autofmt_xdate()
t = np.array(t).astype('datetime64[s]')
ax.plot(t, res, label='radiation @ {}'.format(latitude_deg))
ax.xaxis_date()
plt.autoscale(enable=True, axis='x', tight=True)
plt.xlabel('hour', fontsize=15)
plt.ylabel('W/m2', fontsize=15)
plt.title('Clear sky radiation on normal plane', fontsize=18)
plt.legend()
plt.grid(True)
plt.show()

ds0 = utc.time(2017, 6, 14, 0)

res = []
t = []
for i in range(0, 23):
    ds0 = utc.time(2017, 6, 14, i, 0, 0)
    ds1 = utc.time(2017, 6, 14, i+1, 0, 0)
    t.append(int(ds0))
    # adcal.net_radiation_step(radres, latitude_deg, ds, api.deltaminutes(10), 0, 0, 0, 60, elevation, 0)
    radcal.net_radiation_step(radres, latitude_deg, ds0,deltahours(1), 0, 0, 0, 60, elevation, 0)
    res.append(radres.sw_cs_p)

fig, ax = plt.subplots(figsize=(18, 6))
fig.autofmt_xdate()
t = np.array(t).astype('datetime64[s]')
ax.plot(t, res, label='radiation @ {}'.format(latitude_deg))
ax.xaxis_date()
plt.autoscale(enable=True, axis='x', tight=True)
plt.xlabel('hour', fontsize=15)
plt.ylabel('W/m2', fontsize=15)
plt.title('Clear sky radiation on normal plane', fontsize=18)
plt.legend()
plt.grid(True)
plt.show()