wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-11 20:44:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #20 from simpeg/PropMap

Browse Source 
Prop map
This commit is contained in:
Lindsey
2015-06-02 14:15:22 -07:00
parent 3df2140a88 66b26a1b5e
commit c28e2c10de
5 changed files with 144 additions and 221 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/Base.py
+50 -113
View File
@@ -1,6 +1,14 @@
from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver
from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
class EMPropMap(Maps.PropMap):
sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap))
mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))
class BaseEMProblem(Problem.BaseProblem):
def __init__(self, mesh, **kwargs):
@@ -9,6 +17,8 @@ class BaseEMProblem(Problem.BaseProblem):
surveyPair = Survey.BaseSurvey
dataPair = Survey.Data
PropMap = EMPropMap
Solver = SimpegSolver
solverOpts = {}
@@ -24,94 +34,20 @@ class BaseEMProblem(Problem.BaseProblem):
self.__makeASymmetric = True
return self.__makeASymmetric
####################################################
# Phys Props
####################################################
# Mu
@property
def mu(self):
if getattr(self, '_mu', None) is None:
# if getattr(self, '_mui', None) is not None:
# self._mu = sel
self._mu = mu_0
return self._mu
@mu.setter
def mu(self, value):
if getattr(self, '_MfMui', None) is not None:
del self._MfMui
if getattr(self, '_MfMuiI', None) is not None:
del self._MfMuiI
if getattr(self, '_MeMu', None) is not None:
del delf._MeMu
if getattr(self, '_MeMuI', None) is not None:
del self._MeMuI
self._mu = value
# TODO: hardcoded to assume diagonal mu
@property
def mui(self):
if getattr(self, '_mui', None) is None:
self._mui = 1./mu_0
return self._mui
@mui.setter
def mui(self, value):
if getattr(self, '_MfMui', None) is not None:
del self._MfMui
if getattr(self, '_MfMuiI', None) is not None:
del self._MfMuiI
if getattr(self, '_MeMu', None) is not None:
del delf._MeMu
if getattr(self, '_MeMuI', None) is not None:
del self._MeMuI
self._mui = value
# Sigma
# deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI']
@property
def sigma(self):
if getattr(self, '_sigma', None) is None:
self._sigma = self.curModel.transform
return self._sigma
@sigma.setter
def sigma(self, value):
if getattr(self, '_MeSigma', None) is not None:
del self._MeSigma
if getattr(self, '_MeSigmaI', None) is not None:
del self._MeSigmaI
if getattr(self, '_MfSigmai', None) is not None:
del delf._MfSigmai
if getattr(self, '_MfSigmaiI', None) is not None:
del self._MfSigmaiI
self._sigma = value
# def dsigma_dm(self):
# return self.curModel.transformDeriv
# TODO: hardcoded to assume diagonal sigma
@property
def sigmai(self):
if getattr(self, '_sigmai', None) is None:
self._sigmai = 1./self.curModel.transform
return self._sigmai
@sigmai.setter
def sigmai(self, value):
if getattr(self, '_MeSigma', None) is not None:
del self._MeSigma
if getattr(self, '_MeSigmaI', None) is not None:
del self._MeSigmaI
if getattr(self, '_MfSigmai', None) is not None:
del delf._MfSigmai
if getattr(self, '_MfSigmaiI', None) is not None:
del self._MfSigmaiI
self._sigma = value
####################################################
# Mass Matrices
####################################################
@property
def deleteTheseOnModelUpdate(self):
toDelete = []
if self.mapping.sigmaMap is not None or self.mapping.rhoMap is not None:
toDelete += ['_MeSigma', '_MeSigmaI','_MfRho','_MfRhoI']
if self.mapping.muMap is not None or self.mapping.muiMap is not None:
toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI']
return toDelete
@property
def Me(self):
if getattr(self, '_Me', None) is None:
@@ -129,25 +65,25 @@ class BaseEMProblem(Problem.BaseProblem):
@property
def MfMui(self):
if getattr(self, '_MfMui', None) is None:
self._MfMui = self.mesh.getFaceInnerProduct(self.mui)
self._MfMui = self.mesh.getFaceInnerProduct(self.curModel.mui)
return self._MfMui
@property
def MfMuiI(self):
if getattr(self, '_MfMuiI', None) is None:
self._MfMuiI = self.mesh.getFaceInnerProduct(self.mui, invMat=True)
self._MfMuiI = self.mesh.getFaceInnerProduct(self.curModel.mui, invMat=True)
return self._MfMuiI
@property
def MeMu(self):
if getattr(self, '_MeMu', None) is None:
self._MeMu = self.mesh.getEdgeInnerProduct(self.mu)
self._MeMu = self.mesh.getEdgeInnerProduct(self.curModel.mu)
return self._MeMu
@property
def MeMuI(self):
if getattr(self, '_MeMuI', None) is None:
self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True)
self._MeMuI = self.mesh.getEdgeInnerProduct(self.curModel.mu, invMat=True)
return self._MeMuI
# ----- Electrical Conductivity ----- #
@@ -155,42 +91,43 @@ class BaseEMProblem(Problem.BaseProblem):
@property
def MeSigma(self):
if getattr(self, '_MeSigma', None) is None:
self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma)
self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma)
return self._MeSigma
# def dMeSigma_dsigma(self, u):
# return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u)
def MeSigmaDeriv(self, u):
"""
Deriv of MeSigma wrt sigma
"""
return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u)
@property
def MeSigmaI(self):
if getattr(self, '_MeSigmaI', None) is None:
self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.sigma, invMat=True)
self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True)
return self._MeSigmaI
# def dMeSigmaI_dsigma(self,u)
def MeSigmaIDeriv(self, u):
"""
Deriv of MeSigma wrt sigma
"""
return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
@property
def MfSigmai(self):
if getattr(self, '_MfSigmai', None) is None:
self._MfSigmai = self.mesh.getFaceInnerProduct(self.sigmai)
return self._MfSigmai
def MfRho(self):
if getattr(self, '_MfRho', None) is None:
self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho)
return self._MfRho
# def dMfSigmai_dsigmai(self,u)
def MfRhoDeriv(self,u):
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
@property
def MfSigmaiI(self):
if getattr(self, '_MfSigmaiI', None) is None:
self._MfSigmaiI = self.mesh.getFaceInnerProduct(self.sigmai, invMat=True)
return self._MfSigmaiI
def MfRhoI(self):
if getattr(self, '_MfRhoI', None) is None:
self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True)
return self._MfRhoI
# def dMfSigmaiI(self,u)
####################################################
# Fields
####################################################
def fields(self, m):
self.curModel = m
F = self.forward(m)
return F
def dMfRhoIDeriv(self,u):
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)
simpegEM/FDEM/FDEM.py
+19 -17
View File
@@ -19,8 +19,8 @@ class BaseFDEMProblem(BaseEMProblem):
surveyPair = SurveyFDEM
fieldsPair = FieldsFDEM
def forward(self, m):
def fields(self, m):
self.curModel = m
F = self.fieldsPair(self.mesh, self.survey)
for freq in self.survey.freqs:
@@ -175,9 +175,11 @@ class ProblemFDEM_e(BaseFDEMProblem):
def getADeriv(self, freq, u, v, adjoint=False):
sig = self.sigma
dsig_dm = self.curModel.transformDeriv
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u)
# sig = self.sigma
# dsig_dm = self.curModel.transformDeriv
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u)
dsig_dm = self.curModel.sigmaDeriv
dMe_dsig = self.MeSigmaDeriv(u)
if adjoint:
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
@@ -327,14 +329,14 @@ class ProblemFDEM_j(BaseFDEMProblem):
"""
MeMuI = self.MeMuI
MfSigmai = self.MfSigmai
MfRho = self.MfRho
C = self.mesh.edgeCurl
iomega = 1j * omega(freq) * sp.eye(self.mesh.nF)
A = C * MeMuI * C.T * MfSigmai + iomega
A = C * MeMuI * C.T * MfRho + iomega
if self._makeASymmetric is True:
return MfSigmai.T*A
return MfRho.T*A
return A
@@ -347,7 +349,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
"""
MeMuI = self.MeMuI
MfSigmai = self.MfSigmai
MfRho = self.MfRho
C = self.mesh.edgeCurl
sigi = self.sigmai
dsig_dm = self.curModel.transformDeriv
@@ -356,11 +358,11 @@ class ProblemFDEM_j(BaseFDEMProblem):
if adjoint:
if self._makeASymmetric is True:
v = MfSigmai * v
v = MfRho * v
return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) )
if self._makeASymmetric is True:
return MfSigmai.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) )
return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) )
return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) )
@@ -378,8 +380,8 @@ class ProblemFDEM_j(BaseFDEMProblem):
RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e
if self._makeASymmetric is True:
MfSigmai = self.MfSigmai
return MfSigmai.T*RHS
MfRho = self.MfRho
return MfRho.T*RHS
return RHS
@@ -429,10 +431,10 @@ class ProblemFDEM_h(BaseFDEMProblem):
"""
MeMu = self.MeMu
MfSigmai = self.MfSigmai
MfRho = self.MfRho
C = self.mesh.edgeCurl
return C.T * MfSigmai * C + 1j*omega(freq)*MeMu
return C.T * MfRho * C + 1j*omega(freq)*MeMu
def getADeriv(self, freq, u, v, adjoint=False):
@@ -458,9 +460,9 @@ class ProblemFDEM_h(BaseFDEMProblem):
S_m, S_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MfSigmai = self.MfSigmai
MfRho = self.MfRho
RHS = S_m + C.T * ( MfSigmai * S_e )
RHS = S_m + C.T * ( MfRho * S_e )
return RHS
simpegEM/FDEM/FieldsFDEM.py
+33 -33
View File
@@ -24,9 +24,9 @@ class FieldsFDEM_e(FieldsFDEM):
def _e(self, e_sol, srcList):
e = e_sol
for i, src in enumerate(srcList):
e_p = src.e_p(self.survey.prob)
if e_p is not None:
e[:,i] += e_p
ePrimary = src.ePrimary(self.survey.prob)
if ePrimary is not None:
e[:,i] += ePrimary
return e
def _b(self, e_sol, srcList):
@@ -38,9 +38,9 @@ class FieldsFDEM_e(FieldsFDEM):
if S_m is not None:
b[:,i] += 1./(1j*omega(src.freq)) * S_m
b_p = src.b_p(self.survey.prob)
if b_p is not None:
b[:,i] += b_p
bPrimary = src.bPrimary(self.survey.prob)
if bPrimary is not None:
b[:,i] += bPrimary
return b
@@ -72,9 +72,9 @@ class FieldsFDEM_b(FieldsFDEM):
b = b_sol
for i, src in enumerate(srcList):
b_p = src.b_p(self.survey.prob)
if b_p is not None:
b[:,i] += b_p
bPrimary = src.bPrimary(self.survey.prob)
if bPrimary is not None:
b[:,i] += bPrimary
return b
def _e(self, b_sol, srcList):
@@ -85,9 +85,9 @@ class FieldsFDEM_b(FieldsFDEM):
if S_e is not None:
e += -self._MeSigmaI*S_e
e_p = src.e_p(self.survey.prob)
if e_p is not None:
e[:,i] += e_p
ePrimary = src.ePrimary(self.survey.prob)
if ePrimary is not None:
e[:,i] += ePrimary
return e
@@ -114,23 +114,23 @@ class FieldsFDEM_j(FieldsFDEM):
def startup(self):
self._edgeCurl = self.survey.prob.mesh.edgeCurl
self._MeMuI = self.survey.prob.MeMuI
self._MfSigmai = self.survey.prob.MfSigmai
self._MfRho = self.survey.prob.MfRho
self._curModel = self.survey.prob.curModel
def _j(self, j_sol, srcList):
j = j_sol
for i, src in enumerate(srcList):
j_p = src.j_p(self.survey.prob)
if j_p is not None:
j[:,i] += j_p
jPrimary = src.jPrimary(self.survey.prob)
if jPrimary is not None:
j[:,i] += jPrimary
return j
def _h(self, j_sol, srcList):
MeMuI = self._MeMuI
C = self._edgeCurl
MfSigmai = self._MfSigmai
MfRho = self._MfRho
h = MeMuI * (C.T * (MfSigmai * j_sol) )
h = MeMuI * (C.T * (MfRho * j_sol) )
for i, src in enumerate(srcList):
h[:,i] *= -1./(1j*omega(src.freq))
@@ -138,9 +138,9 @@ class FieldsFDEM_j(FieldsFDEM):
if S_m is not None:
h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m
h_p = src.h_p(self.survey.prob)
if h_p is not None:
h[:,i] += h_p
hPrimary = src.hPrimary(self.survey.prob)
if hPrimary is not None:
h[:,i] += hPrimary
return h
@@ -151,7 +151,7 @@ class FieldsFDEM_j(FieldsFDEM):
dsig_dm = self._curModel.transformDeriv
dsigi_dsig = -Utils.sdiag(sigi)**2
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j)
sigi = self._MfSigmai
sigi = self._MfRho
S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint)
@@ -177,14 +177,14 @@ class FieldsFDEM_h(FieldsFDEM):
def startup(self):
self._edgeCurl = self.survey.prob.mesh.edgeCurl
self._MeMuI = self.survey.prob.MeMuI
self._MfSigmai = self.survey.prob.MfSigmai
self._MfRho = self.survey.prob.MfRho
def _h(self, h_sol, srcList):
h = h_sol
for i, src in enumerate(srcList):
h_p = src.h_p(self.survey.prob)
if h_p is not None:
h[:,i] += h_p
hPrimary = src.hPrimary(self.survey.prob)
if hPrimary is not None:
h[:,i] += hPrimary
return h
def _j(self, h_sol, srcList):
@@ -194,9 +194,9 @@ class FieldsFDEM_h(FieldsFDEM):
if S_e is not None:
j[:,i] += -S_e
j_p = src.j_p(self.survey.prob)
if j_p is not None:
j[:,i] += j_p
jPrimary = src.jPrimary(self.survey.prob)
if jPrimary is not None:
j[:,i] += jPrimary
return j
def _jDeriv(self, h_sol, srcList, v, adjoint=False):
@@ -215,11 +215,11 @@ class FieldsFDEM_h(FieldsFDEM):
# elif fieldType == 'h':
# MeMuI = self._MeMuI
# C = self.mesh.edgeCurl
# MfSigmai = self._MfSigmai
# MfRho = self._MfRho
# if not adjoint:
# h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigmai * j ) )
# h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfRho * j ) )
# else:
# h = -(1./(1j*omega(freq))) * MfSigmai.T * ( C * ( MeMuI.T * j ) )
# h = -(1./(1j*omega(freq))) * MfRho.T * ( C * ( MeMuI.T * j ) )
# return h
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
@@ -235,7 +235,7 @@ class FieldsFDEM_h(FieldsFDEM):
# dsig_dm = self._curModel.transformDeriv
# dsigi_dsig = -Utils.sdiag(sigi)**2
# dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j)
# sigi = self._MfSigmai
# sigi = self._MfRho
# if not adjoint:
# return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) )
# else:
simpegEM/FDEM/SurveyFDEM.py
+39 -55
View File
@@ -96,51 +96,35 @@ class SrcFDEM(Survey.BaseSrc):
rxPair = RxFDEM
def eval(self, prob):
S_m = self._getS_m(prob)
S_e = self._getS_e(prob)
S_m = self.S_m(prob)
S_e = self.S_e(prob)
return S_m, S_e
def evalDeriv(self, prob, v, adjoint=None):
return self._getS_mDeriv(prob,v,adjoint), self._getS_eDeriv(prob,v,adjoint)
return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint)
def b_p(self,prob):
b_p = self._getb_p(prob)
return b_p
def bPrimary(self,prob):
return None
def h_p(self,prob):
h_p = self._geth_p(prob)
return h_p
def e_p(self,prob):
e_p = self._gete_p(prob)
return e_p
def j_p(self,prob):
j_p = self._getj_p(prob)
return j_p
def _getb_p(self,prob):
def hPrimary(self,prob):
return None
def _geth_p(self,prob):
def ePrimary(self,prob):
return None
def _gete_p(self,prob):
def jPrimary(self,prob):
return None
def _getj_p(self,prob):
def S_m(self,prob):
return None
def _getS_m(self,prob):
def S_e(self,prob):
return None
def _getS_e(self,prob):
def S_mDeriv(self, prob, v, adjoint = False):
return None
def _getS_mDeriv(self, prob, v, adjoint = False):
return None
def _getS_eDeriv(self, prob, v, adjoint = False):
def S_eDeriv(self, prob, v, adjoint = False):
return None
@@ -154,12 +138,12 @@ class SrcFDEM_RawVec_e(SrcFDEM):
"""
def __init__(self, rxList, freq, S_e):
self.S_e = np.array(S_e,dtype=float)
self._S_e = np.array(S_e,dtype=float)
self.freq = float(freq)
SrcFDEM.__init__(self, rxList)
def _getS_e(self, prob):
return self.S_e
def S_e(self, prob):
return self._S_e
class SrcFDEM_RawVec_m(SrcFDEM):
@@ -172,12 +156,12 @@ class SrcFDEM_RawVec_m(SrcFDEM):
"""
def __init__(self, rxList, freq, S_m):
self.S_m = np.array(S_m,dtype=float)
self._S_m = np.array(S_m,dtype=float)
self.freq = float(freq)
SrcFDEM.__init__(self, rxList)
def _getS_m(self, prob):
return self.S_m
def S_m(self, prob):
return self._S_m
class SrcFDEM_RawVec(SrcFDEM):
@@ -190,16 +174,16 @@ class SrcFDEM_RawVec(SrcFDEM):
:param rxList: receiver list
"""
def __init__(self, rxList, freq, S_m, S_e):
self.S_m = np.array(S_m,dtype=float)
self.S_e = np.array(S_e,dtype=float)
self._S_m = np.array(S_m,dtype=float)
self._S_e = np.array(S_e,dtype=float)
self.freq = float(freq)
SrcFDEM.__init__(self, rxList)
def _getS_m(self,prob):
return self.S_m
def S_m(self,prob):
return self._S_m
def _getS_e(self,prob):
return self.S_e
def S_e(self,prob):
return self._S_e
class SrcFDEM_MagDipole(SrcFDEM):
@@ -212,7 +196,7 @@ class SrcFDEM_MagDipole(SrcFDEM):
self.moment = moment
SrcFDEM.__init__(self, rxList)
def _getb_p(self,prob):
def bPrimary(self,prob):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
@@ -243,12 +227,12 @@ class SrcFDEM_MagDipole(SrcFDEM):
return C*a
def _geth_p(self,prob):
b = self._getb_p(prob)
def hPrimary(self,prob):
b = self.bPrimary(prob)
return h_from_b(prob,b)
def _getS_m(self,prob):
b_p = self._getb_p(prob)
def S_m(self,prob):
b_p = self.bPrimary(prob)
return -1j*omega(self.freq)*b_p
@@ -264,7 +248,7 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM):
self.moment = moment
SrcFDEM.__init__(self, rxList)
def _getb_p(self,prob):
def bPrimary(self,prob):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
@@ -295,12 +279,12 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM):
return b
def _geth_p(self,prob):
b = self._getb_p(prob)
def hPrimary(self,prob):
b = self.bPrimary(prob)
return h_from_b(prob, b)
def _getS_m(self,prob):
b = self._getb_p(prob)
def S_m(self,prob):
b = self.bPrimary(prob)
return -1j*omega(self.freq)*b
@@ -313,7 +297,7 @@ class SrcFDEM_CircularLoop(SrcFDEM):
self.radius = radius
SrcFDEM.__init__(self, rxList)
def _getb_p(self,prob):
def bPrimary(self,prob):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
@@ -343,12 +327,12 @@ class SrcFDEM_CircularLoop(SrcFDEM):
return C*a
def _geth_p(self,prob):
b = self._getb_p(prob)
def hPrimary(self,prob):
b = self.bPrimary(prob)
return h_from_b
def _getS_m(self, prob):
b = self._getb_p(prob)
def S_m(self, prob):
b = self.bPrimary(prob)
return -1j*omega(self.freq)*b
simpegEM/Utils/EMUtils.py
+3 -3
View File
@@ -16,7 +16,7 @@ def e_from_j(prob,j):
if eqLocs is 'FE':
MSigmaI = prob.MeSigmaI
elif eqLocs is 'EF':
MSigmaI = prob.MfSigmai
MSigmaI = prob.MfRho
return MSigmaI*j
def j_from_e(prob,e):
@@ -24,13 +24,13 @@ def j_from_e(prob,e):
if eqLocs is 'FE':
MSigma = prob.MeSigma
elif eqLocs is 'EF':
MSigma = prob.MfSigma
MSigma = prob.MfRhoI
return MSigma*e
def b_from_h(prob,h):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MMu = prob.MfMu
MMu = prob.MfMuiI
elif eqLocs is 'EF':
MMu = prob.MeMu
return MMu*h