nstools.py

##Surfaces is a dictionary which has keys for each neutral surface
##    at each neutral surface it has a multi dimensional array with this structure:
##    [[lon],[lat],[nsdepth],[[nsdepth],[temp],[salinity],[PV],[uz],[vz]]]
##    I plan on reworking it into a dictionary to make my code more clear soon
##Neighbors is a dictionary which has keys for each neutral surface
##    every key has a value which is another dictionary which has one key
##    which is the point at which we are calculating the derivative. The value  
##    connected to that key is an array of the indexs of the four nearest points to that point

def componentDistance(surfaces,k,i1,i2):
    #x = surfaces[k][0][i1] - surfaces[k][0][i2]
    #if abs(x) > abs(360 - (x+360)%360):
        #x = np.sign(x)*(360-(x+360)%360)
    if  (surfaces[k][0][i1]+180 ) > (surfaces[k][0][i2]+180):
        x = surfaces[k][0][i1]+180 - (surfaces[k][0][i2]+180)
        if x>180:
            x = -(360-x)
    else:
        x = surfaces[k][0][i1]+180 - (surfaces[k][0][i2]+180)
        if x < -180:
            x= -(-360-x)
    x=x*np.cos(np.deg2rad(surfaces[k][1][i2]))*111000.0
    #print(surfaces[k][0][i1],surfaces[k][0][i2],x)

    y = (surfaces[k][1][i1]-surfaces[k][1][i2])*111000.0
    return x,y


def addPrimeToSurfaces(surfaces,neighbors,debug=False):
    for k in surfaces.keys():
        surfaces[k][2].append(np.zeros(len(surfaces[k][1])))
        surfaces[k][2].append(np.zeros(len(surfaces[k][1])))
    alldxs = []
    for k in neighbors.keys():
        print("adding primes to: ",k)
        for r in neighbors[k].keys():
            #alright so here k is our NS
            #r is the index of the point for which we are calculating these prime values
            #adjacent is the list of adjacent points
            adjacent = neighbors[k][r]
            dxsum = 0
            dysum = 0
            dxs = []
            dys = []
            dhs = []
            dists = 0
            for i in adjacent:
                dx,dy = componentDistance(surfaces,k,i,r)
                dxs.append(dx)
                dys.append(dy)
            dxindexs = [np.argmin(dxs),np.argmax(dxs)]
            dyindexs = [np.argmin(dys),np.argmax(dys)]
            dxfinal,b = componentDistance(surfaces,k,adjacent[dxindexs[1]],adjacent[dxindexs[0]])
            b,dyfinal = componentDistance(surfaces,k,adjacent[dyindexs[1]],adjacent[dyindexs[0]])
            #print(r,adjacent)
            dhx = surfaces[k][2][0][adjacent[dxindexs[1]]] - surfaces[k][2][0][adjacent[dxindexs[0]]]
            #dhx = surfaces[k][2][0][adjacent[dxindexs[1]]]
            dhy = surfaces[k][2][0][adjacent[dyindexs[1]]]-surfaces[k][2][0][adjacent[dyindexs[0]]]
            #dhy = surfaces[k][2][0][adjacent[dyindexs[1]]] - surfaces[k][2][0][adjacent[dyindexs[0]]]
            dhdtheta = dhx/dxfinal
            dhdr = dhy/dyfinal
            surfaces[k][2][4][r] = dhdtheta
            surfaces[k][2][5][r] = dhdr 
    return surfaces