<s> from SimPEG import * <EOL> from scipy . special import ellipk , ellipe <EOL> from scipy . constants import mu_0 , pi <EOL> def MagneticDipoleVectorPotential ( srcLoc , obsLoc , component , moment = <NUM_LIT:1.> , dipoleMoment = ( <NUM_LIT:0.> , <NUM_LIT:0.> , <NUM_LIT:1.> ) , mu = mu_0 ) : <EOL> """<STR_LIT>""" <EOL> if type ( component ) in [ list , tuple ] : <EOL> out = range ( len ( component ) ) <EOL> for i , comp in enumerate ( component ) : <EOL> out [ i ] = MagneticDipoleVectorPotential ( srcLoc , obsLoc , comp , dipoleMoment = dipoleMoment ) <EOL> return np . concatenate ( out ) <EOL> if isinstance ( obsLoc , Mesh . BaseMesh ) : <EOL> mesh = obsLoc <EOL> assert component in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] , "<STR_LIT>" <EOL> return MagneticDipoleVectorPotential ( srcLoc , getattr ( mesh , '<STR_LIT>' + component ) , component [ <NUM_LIT:1> ] , dipoleMoment = dipoleMoment ) <EOL> if component == '<STR_LIT:x>' : <EOL> dimInd = <NUM_LIT:0> <EOL> elif component == '<STR_LIT:y>' : <EOL> dimInd = <NUM_LIT:1> <EOL> elif component == '<STR_LIT:z>' : <EOL> dimInd = <NUM_LIT:2> <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> srcLoc = np . atleast_2d ( srcLoc ) <EOL> obsLoc = np . atleast_2d ( obsLoc ) <EOL> dipoleMoment = np . atleast_2d ( dipoleMoment ) <EOL> nEdges = obsLoc . shape [ <NUM_LIT:0> ] <EOL> nSrc = srcLoc . shape [ <NUM_LIT:0> ] <EOL> m = np . array ( dipoleMoment ) . repeat ( nEdges , axis = <NUM_LIT:0> ) <EOL> A = np . empty ( ( nEdges , nSrc ) ) <EOL> for i in range ( nSrc ) : <EOL> dR = obsLoc - srcLoc [ i , np . newaxis ] . repeat ( nEdges , axis = <NUM_LIT:0> ) <EOL> mCr = np . cross ( m , dR ) <EOL> r = np . sqrt ( ( dR ** <NUM_LIT:2> ) . sum ( axis = <NUM_LIT:1> ) ) <EOL> A [ : , i ] = + ( mu / ( <NUM_LIT:4> * pi ) ) * mCr [ : , dimInd ] / ( r ** <NUM_LIT:3> ) <EOL> if nSrc == <NUM_LIT:1> : <EOL> return A . flatten ( ) <EOL> return A <EOL> def MagneticDipoleFields ( srcLoc , obsLoc , component , moment = <NUM_LIT:1.> , mu = mu_0 ) : <EOL> """<STR_LIT>""" <EOL> if component == '<STR_LIT:x>' : <EOL> dimInd = <NUM_LIT:0> <EOL> elif component == '<STR_LIT:y>' : <EOL> dimInd = <NUM_LIT:1> <EOL> elif component == '<STR_LIT:z>' : <EOL> dimInd = <NUM_LIT:2> <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> srcLoc = np . atleast_2d ( srcLoc ) <EOL> obsLoc = np . atleast_2d ( obsLoc ) <EOL> moment = np . atleast_2d ( moment ) <EOL> nFaces = obsLoc . shape [ <NUM_LIT:0> ] <EOL> nSrc = srcLoc . shape [ <NUM_LIT:0> ] <EOL> m = np . array ( moment ) . repeat ( nFaces , axis = <NUM_LIT:0> ) <EOL> B = np . empty ( ( nFaces , nSrc ) ) <EOL> for i in range ( nSrc ) : <EOL> dR = obsLoc - srcLoc [ i , np . newaxis ] . repeat ( nFaces , axis = <NUM_LIT:0> ) <EOL> r = np . sqrt ( ( dR ** <NUM_LIT:2> ) . sum ( axis = <NUM_LIT:1> ) ) <EOL> if dimInd == <NUM_LIT:0> : <EOL> B [ : , i ] = + ( mu / ( <NUM_LIT:4> * pi ) ) / ( r ** <NUM_LIT:3> ) * ( <NUM_LIT:3> * dR [ : , <NUM_LIT:2> ] * dR [ : , <NUM_LIT:0> ] / r ** <NUM_LIT:2> ) <EOL> elif dimInd == <NUM_LIT:1> : <EOL> B [ : , i ] = + ( mu / ( <NUM_LIT:4> * pi ) ) / ( r ** <NUM_LIT:3> ) * ( <NUM_LIT:3> * dR [ : , <NUM_LIT:2> ] * dR [ : , <NUM_LIT:1> ] / r ** <NUM_LIT:2> ) <EOL> elif dimInd == <NUM_LIT:2> : <EOL> B [ : , i ] = + ( mu / ( <NUM_LIT:4> * pi ) ) / ( r ** <NUM_LIT:3> ) * ( <NUM_LIT:3> * dR [ : , <NUM_LIT:2> ] ** <NUM_LIT:2> / r ** <NUM_LIT:2> - <NUM_LIT:1> ) <EOL> else : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> if nSrc == <NUM_LIT:1> : <EOL> return B . flatten ( ) <EOL> return B <EOL> def MagneticLoopVectorPotential ( srcLoc , obsLoc , component , radius , mu = mu_0 ) : <EOL> """<STR_LIT>""" <EOL> if type ( component ) in [ list , tuple ] : <EOL> out = range ( len ( component ) ) <EOL> for i , comp in enumerate ( component ) : <EOL> out [ i ] = MagneticLoopVectorPotential ( srcLoc , obsLoc , comp , radius , mu ) <EOL> return np . concatenate ( out ) <EOL> if isinstance ( obsLoc , Mesh . BaseMesh ) : <EOL> mesh = obsLoc <EOL> assert component in [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] , "<STR_LIT>" <EOL> return MagneticLoopVectorPotential ( srcLoc , getattr ( mesh , '<STR_LIT>' + component ) , component [ <NUM_LIT:1> ] , radius , mu ) <EOL> srcLoc = np . atleast_2d ( srcLoc ) <EOL> obsLoc = np . atleast_2d ( obsLoc ) <EOL> n = obsLoc . shape [ <NUM_LIT:0> ] <EOL> nSrc = srcLoc . shape [ <NUM_LIT:0> ] <EOL> if component == '<STR_LIT:z>' : <EOL> A = np . zeros ( ( n , nSrc ) ) <EOL> if nSrc == <NUM_LIT:1> : <EOL> return A . flatten ( ) <EOL> return A <EOL> else : <EOL> A = np . zeros ( ( n , nSrc ) ) <EOL> for i in range ( nSrc ) : <EOL> x = obsLoc [ : , <NUM_LIT:0> ] - srcLoc [ i , <NUM_LIT:0> ] <EOL> y = obsLoc [ : , <NUM_LIT:1> ] - srcLoc [ i , <NUM_LIT:1> ] <EOL> z = obsLoc [ : , <NUM_LIT:2> ] - srcLoc [ i , <NUM_LIT:2> ] <EOL> r = np . sqrt ( x ** <NUM_LIT:2> + y ** <NUM_LIT:2> ) <EOL> m = ( <NUM_LIT:4> * radius * r ) / ( ( radius + r ) ** <NUM_LIT:2> + z ** <NUM_LIT:2> ) <EOL> m [ m > <NUM_LIT:1.> ] = <NUM_LIT:1.> <EOL> K = ellipk ( m ) <EOL> E = ellipe ( m ) <EOL> ind = ( r > <NUM_LIT:0> ) & ( m < <NUM_LIT:1> ) <EOL> Aphi = np . zeros ( n ) <EOL> Aphi [ ind ] = <NUM_LIT> / np . sqrt ( m [ ind ] ) * np . sqrt ( radius / r [ ind ] ) * ( ( <NUM_LIT:1.> - m [ ind ] / <NUM_LIT> ) * K [ ind ] - E [ ind ] ) <EOL> if component == '<STR_LIT:x>' : <EOL> A [ ind , i ] = Aphi [ ind ] * ( - y [ ind ] / r [ ind ] ) <EOL> elif component == '<STR_LIT:y>' : <EOL> A [ ind , i ] = Aphi [ ind ] * ( x [ ind ] / r [ ind ] ) <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> if nSrc == <NUM_LIT:1> : <EOL> return A . flatten ( ) <EOL> return A <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> from SimPEG import Mesh <EOL> import matplotlib . pyplot as plt <EOL> cs = <NUM_LIT:20> <EOL> ncx , ncy , ncz = <NUM_LIT> , <NUM_LIT> , <NUM_LIT> <EOL> hx = np . ones ( ncx ) * cs <EOL> hy = np . ones ( ncy ) * cs <EOL> hz = np . ones ( ncz ) * cs <EOL> mesh = Mesh . TensorMesh ( [ hx , hy , hz ] , '<STR_LIT>' ) <EOL> srcLoc = np . r_ [ <NUM_LIT:0.> , <NUM_LIT:0.> , <NUM_LIT:0.> ] <EOL> Ax = MagneticLoopVectorPotential ( srcLoc , mesh . gridEx , '<STR_LIT:x>' , <NUM_LIT:200> ) <EOL> Ay = MagneticLoopVectorPotential ( srcLoc , mesh . gridEy , '<STR_LIT:y>' , <NUM_LIT:200> ) <EOL> Az = MagneticLoopVectorPotential ( srcLoc , mesh . gridEz , '<STR_LIT:z>' , <NUM_LIT:200> ) <EOL> A = np . r_ [ Ax , Ay , Az ] <EOL> B0 = mesh . edgeCurl * A <EOL> J0 = mesh . edgeCurl . T * B0 <EOL> mesh . plotImage ( B0 , vType = '<STR_LIT>' ) <EOL> mesh . plotImage ( B0 , vType = '<STR_LIT>' ) <EOL> mesh . plotImage ( B0 , vType = '<STR_LIT>' ) <EOL> plt . show ( ) </s>
