import math
import cmath

def main():
    wdir = [150,150,160,150,150,160,200,270,290,310,290,320,310,330,310,310,320,320,330,310,310,220,280,260]
    wspd = [24,17, 8,10, 9, 8, 9,14,14,30,34,33,30,24,33,20,17,17,16, 8, 6, 5, 6, 6]
    rwd = ""
    rws = ""

    wdir_nws = []

    for i in range(0,len(wdir)):
        wdir_nws.append((wdir[i] + 180) * (3.14156/180)) # Dir converted to radians

    rwd, rws = r_wind(len(wdir_nws), wdir_nws, wspd)

    print 'RWD =', rwd
    print 'RWS =', rws

    rwd, rws = r_wind_new(len(wdir_nws), wdir_nws, wspd)

    print 'RWD =', rwd
    print 'RWS =', rws

def r_wind(n_rwd, rrdd, rwnd):
    """
    Takes an array of houry wind directions and speeds and returns the
    resultant average wind speed and direction, or None for speed and direction
    if there are fewer than 16 hourly observations.

    n_rwd -- number of hourly wind observations
    rrdd  -- array of wind directions plus 180 degrees, in radians
    rwnd  -- array of wind speeds

    Returns a tuple of resultant wind direction and wind speed.  Both are
    returned as strings.  Wind direction is in tens of degrees.
    """
    rwdir = rwspd = None
    totrx = 0.0
    totry = 0.0

    if n_rwd > 15:
        x = 0
        while x < n_rwd:
            totrx = (rwnd[x]*(math.sin(rrdd[x]))) + totrx # tot x & y coor.
            totry = (rwnd[x]*(math.cos(rrdd[x]))) + totry
            x = x + 1

        if totry == 0.0:
            r_dir = 0
        else:
            r_dir = math.atan(totrx/totry) # avg r_wind dir
        r_dir = r_dir / (3.14156/180) # convert radians back to degrees

        if (totrx*totry) < 0:
            if totrx < 0:
                r_dir = 180 + r_dir
            else:
                r_dir = 360 + r_dir
        else:
            if totrx > 0:
                r_dir = 180 + r_dir

        r_spd = ((totrx**2 + totry**2) / n_rwd**2)**.5

        if r_dir == 0 and r_spd != 0:
            r_dir = 360

        if r_spd == 0:
            r_dir = 0

        wc1 = '%#4.0f' % r_dir
        wi1 = int(wc1[0:2])
        wi2 = int(wc1[2:3])

        if wi2 > 4:
            wi1 = wi1 + 1

        rwdir = '%2i' % wi1
        rwspd = '%#5.1f' % r_spd

    return rwdir, rwspd

def r_wind_new(n_rwd, rrdd, rwnd):
    wind_vector_sum = None

    if n_rwd > 15:
        for i in range(0, n_rwd):
            wind_polar = cmath.rect(rwnd[i], rrdd[i] - math.pi)
            if wind_vector_sum is None:
                wind_vector_sum = wind_polar
            else:
                wind_vector_sum += wind_polar

        r, phi = cmath.polar(wind_vector_sum / n_rwd)

        rwdir = int(round(int(round(math.degrees(phi) % 360)) / 10.0))
        rwspd = int(round(r*10))/10.0
    else:
        return None, None

    return '%2i' % rwdir, '%#5.1f' % rwspd

main()