=================================================
Checking calculation of Radiation at Polar Region
=================================================

:ref:`setup-shyft`

.. code-block:: python

    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()

.. image:: clear_sky_radiation.png

.. code-block:: python

    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()

.. image:: clear_sky_radiation_normal.png

.. code-block:: python

    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()

.. image:: clear_sky_radiation_normal2.png