Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- ##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
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement