Clean up repo.

2026-07-19 11:28:00 +08:00 · 2014-05-16 15:51:57 -07:00
parent 95a837b1fc
commit 48ad667c1f
4 changed files with 0 additions and 1615 deletions
@@ -0,0 +1,72 @@
+import numpy as np
+
+def ismember(a, b):
+    tf = np.array([i in b for i in a])
+    return tf
+
+def path2edgeModel(mesh, pts):
+    edm_x = np.zeros(np.prod(mesh.vnEx))
+    edm_y = np.zeros(np.prod(mesh.vnEy))
+    edm_z = np.zeros(np.prod(mesh.vnEz))
+
+    for ii in range (pts.shape[0]-1):
+        pt1 = pts[ii,:]
+        pt2 = pts[ii+1,:]
+        delta = pt2 - pt1
+        deltaDim = np.argwhere(delta)
+        #assert(np.size(deltaDim)==1), "Path must be orthoginal to mesh"
+        if deltaDim == 0:
+            xLoc = mesh.vectorCCx[(min(pt1[0],pt2[0]) < mesh.vectorCCx ) & (mesh.vectorCCx < max(pt1[0],pt2[0]))]
+            yLoc = pts[ii,1]
+            zLoc = pts[ii,2]
+            delDir = np.sign(pt2[0]-pt1[0])
+            xyz = np.c_[xLoc, np.ones(np.size(xLoc))*yLoc, np.ones(np.size(xLoc))*zLoc]
+            edgeInd=ismember(map(tuple,mesh.gridEx),map(tuple,xyz))
+            edm_x[edgeInd] = delDir
+            # print '>> x-direction', ii
+            # print mesh.gridEx[edgeInd]
+        if deltaDim == 1:
+            xLoc = pts[ii,0]
+            yLoc = mesh.vectorCCy[(min(pt1[1],pt2[1]) < mesh.vectorCCy ) & (mesh.vectorCCy < max(pt1[1],pt2[1]))]
+            zLoc = pts[ii,2]
+            delDir = np.sign(pt2[1]-pt1[1])
+            xyz = np.c_[np.ones(np.size(yLoc))*xLoc, yLoc, np.ones(np.size(yLoc))*zLoc]
+            edgeInd=ismember(map(tuple,mesh.gridEy),map(tuple,xyz))
+            edm_y[edgeInd] = delDir
+            # print '>> y-direction', ii
+            # print mesh.gridEy[edgeInd]
+        if deltaDim == 2:
+            xLoc = pts[ii,0]
+            yLoc = pts[ii,1]
+            zLoc = mesh.vectorCCz[(min(pt1[2],pt2[2]) < mesh.vectorCCz ) & (mesh.vectorCCz < max(pt1[2],pt2[2]))]
+            delDir = np.sign(pt2[2]-pt1[2])
+            xyz = np.c_[np.ones(np.size(zLoc))*xLoc, np.ones(np.size(zLoc))*yLoc, zLoc]
+            edgeInd=ismember(map(tuple,mesh.gridEz),map(tuple,xyz))
+            edm_z[edgeInd] = delDir
+            # print '>> z-direction', ii
+            # print mesh.gridEz[edgeInd]
+
+
+    edgeModel = np.r_[edm_x, edm_y, edm_z]
+    return edgeModel
+
+def MMRhalf(loc1, loc2, x, y):
+    """ Anaytic function for MMR response (B^{1D})
+            - loc1=(x1,y1): x, y location for (+) charge
+            - loc2=(x2,y2): x, y1
+            - x : observation points in x-direction
+            - y : observation points in y-direction
+    """
+    x1=loc1[0]
+    x2=loc2[0]
+    y1=loc1[1]
+    y2=loc2[1]
+    mu0 = 4*np.pi*1e-7
+    I = 1
+
+    rho  = lambda x1, y1, x, y: np.sqrt((x-x1)**2+(y-y1)**2)
+
+    By =mu0*I/(4*np.pi)*np.array((x-x1)/rho(x1,y1,x,y)**2-(x-x2)/rho(x2,y2,x,y)**2)
+    Bx =mu0*I/(4*np.pi)*np.array(-(y-y1)/rho(x1,y1,x,y)**2+(y-y2)/rho(x2,y2,x,y)**2)
+
+    return Bx, By