import math


def sun_NR(doy,lat):
    """
    Function to calculate the maximum sunshine duration N and incoming radiation
    at the top of the atmosphere from day of year and latitude.

    NOTE: Only valid for latitudes between 0 and 67 degrees (tropics and
    temperate zone)

    Input:
        - doy: (array of) day of year
        - lat: latitude in degrees, negative for southern hemisphere

    Output:
        - N: (array of) maximum sunshine hours [h]
        - Rext: (array of) extraterrestrial radiation [J/day]
    """

    # Set solar constant [W/m2]
    S = 1367.0 # [W/m2]
    # Convert latitude [degrees] to radians
    latrad = lat * math.pi / 180.0
    # Determine length of doy array
    #l = scipy.size(doy)
    # Check if we have a single value or an array
    #if l < 2:   # Dealing with single value...
        # calculate solar declination dt [radians]
    dt = 0.409 * math.sin(2 * math.pi / 365 * doy - 1.39)
    # calculate sunset hour angle [radians]
    ws = math.acos(-math.tan(latrad) * math.tan(dt))
    # Calculate sunshine duration N [h]
    N = 24 / math.pi * ws
    # Calculate day angle j [radians]
    j = 2 * math.pi / 365.25 * doy
    # Calculate relative distance to sun
    dr = 1.0 + 0.03344 * math.cos(j - 0.048869)
    # Calculate Rext
    Rext = S * 86400 / math.pi * dr * (ws * math.sin(latrad) * math.sin(dt) + math.sin(ws) * math.cos(latrad) * math.cos(dt))
    return Rext