From 36f8eca25c04e443b1ec32059904013b1414b71f Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 23 Jun 2015 18:25:29 -0700 Subject: [PATCH] Sources can take mu values other than mu_0 --- simpegEM/Utils/SrcUtils.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/simpegEM/Utils/SrcUtils.py b/simpegEM/Utils/SrcUtils.py index ec638a0f..7c672bf4 100644 --- a/simpegEM/Utils/SrcUtils.py +++ b/simpegEM/Utils/SrcUtils.py @@ -2,7 +2,7 @@ from SimPEG import * from scipy.special import ellipk, ellipe from scipy.constants import mu_0, pi -def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)): +def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.), mu = mu_0): """ Calculate the vector potential of a set of magnetic dipoles at given locations 'ref. ' @@ -48,13 +48,13 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0 dR = obsLoc - srcLoc[i, np.newaxis].repeat(nEdges, axis=0) mCr = np.cross(m, dR) r = np.sqrt((dR**2).sum(axis=1)) - A[:, i] = +(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3) + A[:, i] = +(mu/(4*pi)) * mCr[:,dimInd]/(r**3) if nSrc == 1: return A.flatten() return A -def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.): +def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): """ Calculate the vector potential of a set of magnetic dipoles at given locations 'ref. ' @@ -89,11 +89,11 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.): dR = obsLoc - srcLoc[i, np.newaxis].repeat(nFaces, axis=0) r = np.sqrt((dR**2).sum(axis=1)) if dimInd == 0: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2) + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2) elif dimInd == 1: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2) + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2) elif dimInd == 2: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1) + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1) else: raise Exception("Not Implemented") if nSrc == 1: