def tmy_plane_of_array(dni_list, dhi_list, latitude=0, longitude=0,
    timezone=0, plane_azim=180, plane_tilt=0):
    """
    from TMY data calculate the yearly plane-of-array energy yield 

    the function accepts one tmy list for dni (beam) and one tmy list
    for dhi insolation values, the azimuth, tilt values for an oriented
    surface for solar energy collection and a location latitude. It
    returns the PlaneOfArray (POA) irradiation value on an hourly
    basis for the whole TMY year.

    results with reference to the unit surface (1 m2)
    
    input:
    - dni_list, list with 8760 DNI values from TMY dataset, integer/float
    - dhi_list, list with 8760 DHI values from TMY dataset, integer/float
    - latitude, latitude of the location (positive North of equator),
    float, default =0
    - longitude, longitude of the location (positive North of equator),
    float, default =0
    - timezone, positive East of Greenwich, integer/float 
    - plane_azim, azimuth of the oriented plane, in degrees
    (South is 180), float, default = 180
    - plane_tilt, tilt of the oriented plane in degrees
    (horizontal tilt is 0), float, default =0

    output:
    - poa_list, POA 8760 hourly values, float
    - dir_list, direct radiation on plane, 8760 hourly values, float 
    - dif_list, diffuse radiation on plane, 8760 hourly values, float
    - incid_list, beam angle of incidence on plane, 8760 hourly values,
    float

    """
    import math
    import solprim.solartimeposition as stp
    latitude_rad = math.radians(latitude)
    plane_azim_rad = math.radians(plane_azim)
    plane_tilt_sin = math.sin(math.radians(plane_tilt))
    plane_tilt_cos = math.cos(math.radians(plane_tilt))
    dhi_coeff = 0.5 * (1 + plane_tilt_cos)
    year = 2015
    poa_list = []
    dir_list = []
    dif_list = []
    incid_list = []
    month_len = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    month_ptr = [0, 744, 1416, 2160, 2880, 3624, 4344, 5088, 5832, 6552, 
        7296, 8016]
    day_of_year = 0
    for m in range(0, 12):
        m_ptr = month_ptr[m]
        for d in range(0, month_len[m]):
            d_ptr = m_ptr + d * 24
            day_of_year = day_of_year + 1
            sol_decl = stp.solar_declination(day_of_year)
            sol_decl_rad = math.radians(sol_decl)
            eot = stp.equation_of_time(day_of_year)
            for h in range(0, 24):
                h_ptr = d_ptr + h
                hour_poa = 0
                plane_dir = 0
                plane_dif = 0
                incid_angle = 0
                dni = dni_list[h_ptr]
                dhi = dhi_list[h_ptr]
                if dni > 0 or dhi > 0:
                    hra = stp.hour_angle_decimal(h + 0.5, eot, longitude,
                        timezone)
                    hra_rad = math.radians(hra)
                    sol_pos = stp.solar_azim_elev(sol_decl_rad, hra_rad,
                        latitude_rad)
                    sol_elev_sin = sol_pos[4]
                    sol_elev_cos = sol_pos[5]
                    sol_azim_rad = sol_pos[3]
                    cos_theta = stp.incidence_factor(sol_elev_sin,
                        sol_elev_cos, plane_tilt_sin, plane_tilt_cos,
                        sol_azim_rad, plane_azim_rad)
                    incid_angle = math.degrees(math.acos(cos_theta))
                    plane_dir = dni * cos_theta
                    plane_dif = dhi * dhi_coeff
                    hour_poa = plane_dir + plane_dif
                poa_list.append(hour_poa)
                dir_list.append(plane_dir)
                dif_list.append(plane_dif)
                incid_list.append(round(incid_angle, 2))
    return poa_list, dir_list, dif_list, incid_list
