import numpy as np
from scipy import ndimage
#import itertools
from numba import jit, njit
# It has to be converted to py (simply copy this text in a py file named MultiProcFuncAtomMagnetometry.py)
# Multiprocessing functions
######################################
StudyProcessMemory=False
"""
# To study processes
StudyProcessMemory=True
import sys,psutil
try:
RunCodeOnline=True
from google.colab import files
from google.colab import drive
drive.mount('/content/drive', force_remount=True)
pathScriptData='/content/drive/My Drive/ResearchInnovation/Consultancy/JofreFamily_EMFluidInterface/Calculations/AtomMagnetic/'
sys.path.append(pathScriptData)
!cp -r "/content/drive/My Drive/ResearchInnovation/Consultancy/JofreFamily_EMFluidInterface/Calculations/AtomMagnetic/MultiProcFuncAtomMagnetometry.ipynb" '/content/'
except:
RunCodeOnline=False
# To run locally it is feasible using jupiter notebook from winPython and following the instructions: https://research.google.com/colaboratory/local-runtimes.html
# Launching the Jupyter Notebook from winpython cmd: jupyter notebook --NotebookApp.allow_origin='https://colab.research.google.com' --port=8888 --NotebookApp.port_retries=0
# Also, winpython cmd allows to install packages need by using pip install command.
# The generated (highly time computation) files can then be easily updated to the specific google drive for hosted runtime (online) execution.
pathScriptData='C:/GeneratedDataPythonJupyter/'
sys.path.append(pathScriptData)
sys.path.append('./Scripts/')
#######################################################################
"""
"""
def initializerProc(init_data):
#global initializerAux_data
#initializerAux_data=np.zeros_like(init_data,dtype=object)
initializerAux_data = np.asarray(init_data,dtype=object)
"""
class initializerDataClass:
initializerAux_data=np.asarray([],dtype=object)
InitData=initializerDataClass()
def initializerProc(init_data):
InitData.initializerAux_data=np.asarray(init_data,dtype=object)
def ScalarConeCylCoordCalcIntegralsProc_star(a_b):
"""Convert `f([1,2])` to `f(1,2)` call."""
return ScalarConeCylCoordCalcIntegralsProc(*a_b)
# It is not njit supported the np.nansum (but it is already very efficient). jit for np.sum is not more efficient
def DoubleSumProc4D(MatrixLarge,d1Value):
# Zeros in the Distancediff might be present and nansum deals with nans
# It is better to be a sum with the differentiator sinceit is regarding the scan in all positions of rho and phi
#return np.sum(np.sum(MatrixLarge,axis=0),axis=0)
#MatrixShape=MatrixLarge.shape
#NumbaMeshgrid4D(d1Value*np.ones(shape=MatrixShape[0],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[1],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[2],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[3],dtype=MatrixLarge.dtype))[0]
#NumbaMeshgrid3D(d2Value*np.ones(shape=MatrixShape[1],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[2],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[3],dtype=MatrixLarge.dtype))[0]
#return np.sum(NumbaMeshgrid3D(d2Value*np.ones(shape=MatrixShape[1],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[2],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[3],dtype=MatrixLarge.dtype))[0]*np.sum(NumbaMeshgrid4D(d1Value*np.ones(shape=MatrixShape[0],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[1],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[2],dtype=MatrixLarge.dtype), np.ones(shape=MatrixShape[3],dtype=MatrixLarge.dtype))[0]*MatrixLarge,axis=0),axis=0)
return d1Value*np.nansum(np.nansum(MatrixLarge,axis=0),axis=0)
@jit
def NumbaMeshgrid2D(x, y):
xx = np.empty(shape=(x.size, y.size), dtype=x.dtype)
yy = np.empty(shape=(x.size, y.size), dtype=y.dtype)
#for i, j in itertools.product(np.arange(0,x.size,1), np.arange(0,y.size,1)):
# xx[i,j] = x[i]
# yy[i,j] = y[j]
for i in range(x.size):
for j in range(y.size):
xx[i,j] = x[i]
yy[i,j] = y[j]
return xx, yy
@jit
def NumbaMeshgrid3D(x, y, z):
xx = np.empty(shape=(x.size, y.size, z.size), dtype=x.dtype)
yy = np.empty(shape=(x.size, y.size, z.size), dtype=y.dtype)
zz = np.empty(shape=(x.size, y.size, z.size), dtype=z.dtype)
#for i, j, k in itertools.product(np.arange(0,x.size,1), np.arange(0,y.size,1), np.arange(0,z.size,1)):
# xx[i,j,k] = x[i]
# yy[i,j,k] = y[j]
# zz[i,j,k] = z[k]
for i in range(x.size):
for j in range(y.size):
for k in range(z.size):
xx[i,j,k] = x[i]
yy[i,j,k] = y[j]
zz[i,j,k] = z[k]
return xx, yy, zz
@jit
def NumbaMeshgrid3D1axis(x, y, z):
xx = np.empty(shape=(x.size, y.size, z.size), dtype=y.dtype)
#for i, j, k in itertools.product(np.arange(0,x.size,1), np.arange(0,y.size,1), np.arange(0,z.size,1)):
# yy[i,j,k] = y[j]
for i in range(x.size):
for j in range(y.size):
for k in range(z.size):
xx[i,j,k] = x[i]
return xx
@jit
def NumbaMeshgrid4D(x, y, z, w):
xx = np.empty(shape=(x.size, y.size, z.size, w.size), dtype=x.dtype)
yy = np.empty(shape=(x.size, y.size, z.size, w.size), dtype=y.dtype)
zz = np.empty(shape=(x.size, y.size, z.size, w.size), dtype=z.dtype)
ww = np.empty(shape=(x.size, y.size, z.size, w.size), dtype=w.dtype)
#for i, j, k, l in itertools.product(np.arange(0,x.size,1), np.arange(0,y.size,1), np.arange(0,z.size,1), np.arange(0,w.size,1)):
# xx[i,j,k,l] = x[i]
# yy[i,j,k,l] = y[j]
# zz[i,j,k,l] = z[k]
# ww[i,j,k,l] = w[l]
for i in range(x.size):
for j in range(y.size):
for k in range(z.size):
for l in range(w.size):
xx[i,j,k,l] = x[i]
yy[i,j,k,l] = y[j]
zz[i,j,k,l] = z[k]
ww[i,j,k,l] = w[l]
return xx, yy, zz, ww
@jit
def MatrixZeroIndices4D(MatrixLarge):
MatrixShape=MatrixLarge.shape
MatrixComp = np.empty(shape=(MatrixShape), dtype=MatrixLarge.dtype)
#for iIter1, iIter2, iIter3, iIter4 in itertools.product(np.arange(0,MatrixShape[0],1), np.arange(0,MatrixShape[1],1), np.arange(0,MatrixShape[2],1), np.arange(0,MatrixShape[3],1)):
# if (MatrixLarge[iIter1,iIter2,iIter3,iIter4]==0.0):
# MatrixComp[iIter1,iIter2,iIter3,iIter4]=0.0
# else:
# MatrixComp[iIter1,iIter2,iIter3,iIter4]=1.0/MatrixLarge[iIter1,iIter2,iIter3,iIter4]
for iIter1 in range(0,MatrixShape[0],1):
for iIter2 in range(0,MatrixShape[1],1):
for iIter3 in range(0,MatrixShape[2],1):
for iIter4 in range(0,MatrixShape[3],1):
if (MatrixLarge[iIter1,iIter2,iIter3,iIter4]==0.0):
MatrixComp[iIter1,iIter2,iIter3,iIter4]=0.0
else:
MatrixComp[iIter1,iIter2,iIter3,iIter4]=1.0/MatrixLarge[iIter1,iIter2,iIter3,iIter4]
return MatrixComp
def ScalarConeCylCoordCalcIntegralsProc(iIterRho):
if (StudyProcessMemory==True and iIterRho==0):
np.save(pathScriptData+'MultiProcScriptProcessMemory_Scalar_MemInfoBefore.npy',psutil.Process().memory_info().rss)
[RhoZparamiIterTop,CosPhiparamTopPhiparamiIterAux1Top,zparamiIterZTop,RhoZparamiIterBottom,CosPhiparamBottomPhiparamiIterAux1Bottom,zparamiIterZBottom,RhoZparamiIterSide,RhoZparamiIterSideDiff,CosPhiparamSidePhiparamiIterAux1Side,zparamiIterZSide,rhoparam,SumdPhi,SumdZ,DimensionsParamCone,SurfChargetopP,SurfChargebottomP,SurfChargesideP,SurfChargeInsideP,PositionParam2]=InitData.initializerAux_data
if (StudyProcessMemory==True and iIterRho==0):
np.save(pathScriptData+'MultiProcScriptProcessMemory_Scalar_MemInfoMiddle.npy',psutil.Process().memory_info().rss)
# https://ieeexplore.ieee.org/document/9085339
#lenphiparam=len(phiparam)
#lenzparam=len(zparam)
#lenzparamiIterAux2=len(zparamiIterAux2)
tanpi2DimensionsParamCone2=np.tan(np.pi/2.0-DimensionsParamCone[2])
rhoparamiIterRho=rhoparam[iIterRho]
rhoparamiIterRhoSquare=rhoparamiIterRho**2
# The below line is currently computed in the main script to avoid repeating it multiple times. If discommented, then it would be computed here.
#[SurfChargesideP,SurfChargeInsideP]=NJITscalarlenzparamiIterAux2For(lenzparamiIterAux2,lenphiparam,lenzparam,nsideP,nInsidesideP,Lindex0,phiparam,MagnetizationParamRhoP,MagnetizationParamPhiP,MagnetizationParamZP,rhoparam,zparam,DimensionsParamCone,zparamiIterAux2)
SumdRho=rhoparam[1]-rhoparam[0]
#SumdPhi=phiparam[1]-phiparam[0]
#SumdZ=zparam[1]-zparam[0]
#MinGridCyl=(SumdRho**2+(SumdPhi*SumdRho)**2+SumdZ**2)/64.0 # Without np.sqrt
TopBottomSurfaceDif=SumdRho*SumdPhi# The rho multiplicaton is done in the Double SumProc
# Top
## The below line is currently computed in the main script to avoid repeating it multiple times. If discommented, then it would be computed here.
##[RhoZparamiIterTop,PhiparamiIterAux1Top,PhiparamTop,zparamiIterZTop]=NumbaMeshgrid4D(rhoparamiIterAux2rindex,phiparam,phiparam,zparam)
##[RhoparamiIterAux2,PhiparamiIterAux1,Phiparam,zparamiIterZ]=np.meshgrid(rhoparamiIterAux2rindex,phiparamiIterAux1,phiparam,zparam, indexing='ij')
#DiffDistance=rhoparamiIterRhoSquare+RhoZparamiIterTop**2+(zparamiIterZTop-(PositionParam2+DimensionsParamCone[1]))**2-2.0*rhoparamiIterRho*RhoZparamiIterTop*CosPhiparamTopPhiparamiIterAux1Top
#DiffDistanceLogic=(DiffDistance<MinGridCyl)
#if (np.sum(DiffDistanceLogic)>0):
# #print('np.sum(DiffDistanceLogic): '+str(np.sum(DiffDistanceLogic)))
# DiffDistance[DiffDistanceLogic]=MinGridCyl#/1024.0 # Without np.sqrt
#DoubleIntegralTerm=DoubleSumProc4D(RhoZparamiIterTop*SurfChargetopP/np.sqrt(DiffDistance),TopBottomSurfaceDif)
## Asuming there are no 0 distances
DoubleIntegralTerm=DoubleSumProc4D(np.nan_to_num(RhoZparamiIterTop*SurfChargetopP/np.sqrt(rhoparamiIterRhoSquare+RhoZparamiIterTop**2+(zparamiIterZTop-(PositionParam2+DimensionsParamCone[1]))**2-2.0*rhoparamiIterRho*RhoZparamiIterTop*CosPhiparamTopPhiparamiIterAux1Top), copy=False, nan=0.0, posinf=0.0, neginf=0.0),TopBottomSurfaceDif)
# Bottom
## The below line is currently computed in the main script to avoid repeating it multiple times. If discommented, then it would be computed here.
##[RhoZparamiIterBottom,PhiparamiIterAux1Bottom,PhiparamBottom,zparamiIterZBottom]=NumbaMeshgrid4D(rhoparamiIterAux2Rindex,phiparam,phiparam,zparam)
##[RhoparamiIterAux2,PhiparamiIterAux1,Phiparam,zparamiIterZ]=np.meshgrid(rhoparamiIterAux2Rindex,phiparamiIterAux1,phiparam,zparam, indexing='ij')
#DiffDistance=rhoparamiIterRhoSquare+RhoZparamiIterBottom**2+(zparamiIterZBottom-PositionParam2)**2-2.0*rhoparamiIterRho*RhoZparamiIterBottom*CosPhiparamBottomPhiparamiIterAux1Bottom
#DiffDistanceLogic=(DiffDistance<MinGridCyl)
#if (np.sum(DiffDistanceLogic)>0):
# DiffDistance[DiffDistanceLogic]=MinGridCyl#/1024.0 # Without np.sqrt
#DoubleIntegralTerm+=DoubleSumProc4D(RhoZparamiIterBottom*SurfChargebottomP/np.sqrt(DiffDistance),TopBottomSurfaceDif)
## Asuming there are no 0 distances
DoubleIntegralTerm+=DoubleSumProc4D(np.nan_to_num(RhoZparamiIterBottom*SurfChargebottomP/np.sqrt(rhoparamiIterRhoSquare+RhoZparamiIterBottom**2+(zparamiIterZBottom-PositionParam2)**2-2.0*rhoparamiIterRho*RhoZparamiIterBottom*CosPhiparamBottomPhiparamiIterAux1Bottom), copy=False, nan=0.0, posinf=0.0, neginf=0.0),TopBottomSurfaceDif)
SideInSideSurfaceDif=SumdPhi*SumdZ # The rho multiplicaton is done in the Double SumProc
# Side
## The below line is currently computed in the main script to avoid repeating it multiple times. If discommented, then it would be computed here.
##[RhoZparamiIterSide,PhiparamiIterAux1Side,PhiparamSide,zparamiIterZSide]=NumbaMeshgrid4D(zparamiIterAux2,phiparam,phiparam,zparam)
##[ZparamiIterAux2,PhiparamiIterAux1,Phiparam,zparamiIterZ]=np.meshgrid(zparamiIterAux2,phiparamiIterAux1,phiparam,zparam, indexing='ij')
DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2=DimensionsParamCone[0]-RhoZparamiIterSideDiff*tanpi2DimensionsParamCone2
#DiffDistance=rhoparamiIterRhoSquare+DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2**2+(zparamiIterZSide-RhoZparamiIterSide)**2-2.0*rhoparamiIterRho*DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2*CosPhiparamSidePhiparamiIterAux1Side
#DiffDistanceLogic=(DiffDistance<MinGridCyl)
#if (np.sum(DiffDistanceLogic)>0):
# DiffDistance[DiffDistanceLogic]=MinGridCyl#/1204.0 # Without np.sqrt
#DoubleIntegralTerm+=DoubleSumProc4D(DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2*SurfChargesideP/np.sqrt(DiffDistance),SideInSideSurfaceDif)
## Asuming there are no 0 distances
DoubleIntegralTerm+=DoubleSumProc4D(np.nan_to_num(DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2*SurfChargesideP/np.sqrt(rhoparamiIterRhoSquare+DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2**2+(zparamiIterZSide-RhoZparamiIterSide)**2-2.0*rhoparamiIterRho*DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2*CosPhiparamSidePhiparamiIterAux1Side), copy=False, nan=0.0, posinf=0.0, neginf=0.0),SideInSideSurfaceDif)
# Inside
if (DimensionsParamCone[3]>0.0):
##[ZparamiIterAux2,PhiparamiIterAux1,Phiparam,zparamiIterZ]=np.meshgrid(zparamiIterAux2,phiparamiIterAux1,phiparam,zparam, indexing='ij')
DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2=DimensionsParamCone[3]-RhoZparamiIterSideDiff*tanpi2DimensionsParamCone2
#DiffDistance=rhoparamiIterRhoSquare+DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2**2+(zparamiIterZSide-RhoZparamiIterSide)**2-2.0*rhoparamiIterRho*DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2*CosPhiparamSidePhiparamiIterAux1Side
#DiffDistanceLogic=(DiffDistance<MinGridCyl)
#if (np.sum(DiffDistanceLogic)>0):
# DiffDistance[DiffDistanceLogic]=MinGridCyl#/1024.0 # Without np.sqrt
#DoubleIntegralTerm-=DoubleSumProc4D(DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2*SurfChargeInsideP/np.sqrt(DiffDistance),SideInSideSurfaceDif)
## Asuming there are no 0 distances
DoubleIntegralTerm-=DoubleSumProc4D(np.nan_to_num(DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2*SurfChargeInsideP/np.sqrt(rhoparamiIterRhoSquare+DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2**2+(zparamiIterZSide-RhoZparamiIterSide)**2-2.0*rhoparamiIterRho*DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2*CosPhiparamSidePhiparamiIterAux1Side), copy=False, nan=0.0, posinf=0.0, neginf=0.0),SideInSideSurfaceDif)
if (StudyProcessMemory==True and iIterRho==0):
np.save(pathScriptData+'MultiProcScriptProcessMemory_Scalar_DoubleIntegralTerm.npy',DoubleIntegralTerm)
np.save(pathScriptData+'MultiProcScriptProcessMemory_Scalar_MemInfoAfter.npy',psutil.Process().memory_info().rss)
return DoubleIntegralTerm
@jit
def NJITscalarSurfaceFor(rhoparam,phiparam,zparam,MagnetizationParamRhoP,MagnetizationParamPhiP,MagnetizationParamZP,DimensionsParamCone,zparamiIterAux2,ntopP,nbottomP,nsideP,nInsidesideP,Lindex0,Lindex,rindex0,Rindex0,lenrhoparamiIterAux2rindex,lenrhoparamiIterAux2Rindex):
dRhoPart=0.0#(rhoparam[1]-rhoparam[0])/2.0
lenphiparam=len(phiparam)
lenzparam=len(zparam)
lenzparamiIterAux2=len(zparamiIterAux2)
SurfChargetopP=np.empty((lenrhoparamiIterAux2rindex,lenphiparam,lenphiparam,lenzparam),dtype=np.float32)
for iIterAux in range(0,lenrhoparamiIterAux2rindex,1):
rindex0iIterAux=rindex0+iIterAux
SurfChargetopP[iIterAux]=NumbaMeshgrid3D1axis(ntopP[0]*MagnetizationParamRhoP[rindex0iIterAux,:,Lindex]+ntopP[1]*MagnetizationParamPhiP[rindex0iIterAux,:,Lindex]+ntopP[2]*MagnetizationParamZP[rindex0iIterAux,:,Lindex],phiparam,zparam)
SurfChargebottomP=np.empty((lenrhoparamiIterAux2Rindex,lenphiparam,lenphiparam,lenzparam),dtype=np.float32)
for iIterAux in range(0,lenrhoparamiIterAux2Rindex,1):
Rindex0iIterAux=Rindex0+iIterAux
SurfChargebottomP[iIterAux]=NumbaMeshgrid3D1axis(nbottomP[0]*MagnetizationParamRhoP[Rindex0iIterAux,:,Lindex0]+nbottomP[1]*MagnetizationParamPhiP[Rindex0iIterAux,:,Lindex0]+nbottomP[2]*MagnetizationParamZP[Rindex0iIterAux,:,Lindex0],phiparam,zparam)
SurfChargesideP=np.empty((lenzparamiIterAux2,lenphiparam,lenphiparam,lenzparam),dtype=np.float32)
SurfChargeInsideP=np.empty((lenzparamiIterAux2,lenphiparam,lenphiparam,lenzparam),dtype=np.float32)
tanpi2DimensionsParamCone2=np.tan(np.pi/2.0-DimensionsParamCone[2])
for iIterAux in range(0,lenzparamiIterAux2,1):
Lindex0iIterAux=Lindex0+iIterAux
zparamiIterAux2iIterAuxzparamiIterAux20tanpi2DimensionsParamCone2=(zparamiIterAux2[iIterAux]-zparamiIterAux2[0])*tanpi2DimensionsParamCone2
# side
SurfChargesideP[iIterAux]=NumbaMeshgrid3D1axis(nsideP[0]*MagnetizationParamRhoP[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamiIterAux2iIterAuxzparamiIterAux20tanpi2DimensionsParamCone2))),:,Lindex0iIterAux]+nsideP[1]*MagnetizationParamPhiP[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamiIterAux2iIterAuxzparamiIterAux20tanpi2DimensionsParamCone2))),:,Lindex0iIterAux]+nsideP[2]*MagnetizationParamZP[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamiIterAux2iIterAuxzparamiIterAux20tanpi2DimensionsParamCone2))),:,Lindex0iIterAux],phiparam,zparam)
# inside
SurfChargeInsideP[iIterAux]=NumbaMeshgrid3D1axis(nInsidesideP[0]*MagnetizationParamRhoP[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamiIterAux2iIterAuxzparamiIterAux20tanpi2DimensionsParamCone2))),:,Lindex0iIterAux]+nInsidesideP[1]*MagnetizationParamPhiP[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamiIterAux2iIterAuxzparamiIterAux20tanpi2DimensionsParamCone2))),:,Lindex0iIterAux]+nInsidesideP[2]*MagnetizationParamZP[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamiIterAux2iIterAuxzparamiIterAux20tanpi2DimensionsParamCone2))),:,Lindex0iIterAux],phiparam,zparam)
return SurfChargetopP,SurfChargebottomP,SurfChargesideP,SurfChargeInsideP
def ScalarConeCylCoordMagScalarPotCtoCylProc_star(a_b):
"""Convert `f([1,2])` to `f(1,2)` call."""
return ScalarConeCylCoordMagScalarPotCtoCylProc(*a_b)
def ScalarConeCylCoordMagScalarPotCtoCylProc(iIterRho):
[rhoparam,phiparam,GridVolume,MagneticScalarPotentialOrigin]=InitData.initializerAux_data
dRhoPart=0.0#(rhoparam[1]-rhoparam[0])/2.0
lenGridVolume0=len(GridVolume[0])
lenGridVolume1=len(GridVolume[1])
lenGridVolume2=len(GridVolume[2])
#for iIterRho in range(0,len(rhoparam),1):
#for iIterRho in range(0,len(rhoparam),1):
# Interpolate Phi axis to have more points and populate X and Y
InterpolationOrder=3
InterpolationFactor=3
rhoparamiIterRho=rhoparam[iIterRho]+dRhoPart
phiparamInterpolated=(ndimage.zoom(phiparam,InterpolationFactor,mode = 'wrap',order=InterpolationOrder,grid_mode=False)).astype(np.float32)
MagneticScalarPotentialOriginInterpolated=(1.0/float(InterpolationFactor))*(ndimage.zoom(MagneticScalarPotentialOrigin[iIterRho],[InterpolationFactor,1],mode = 'wrap',order=InterpolationOrder,grid_mode=False)).astype(np.float32)
CosphiparamiIterPhi=np.cos(phiparamInterpolated)
SinphiparamiIterPhi=np.sin(phiparamInterpolated)
return NJITForScalarConeCylCoordMagScalarPotCtoCylProc(lenGridVolume0,lenGridVolume1,lenGridVolume2,MagneticScalarPotentialOriginInterpolated,phiparamInterpolated,rhoparamiIterRho,GridVolume[0],GridVolume[1],CosphiparamiIterPhi,SinphiparamiIterPhi)
@jit
def NJITForScalarConeCylCoordMagScalarPotCtoCylProc(lenGridVolume0,lenGridVolume1,lenGridVolume2,MagneticScalarPotentialOriginInterpolated,phiparamInterpolated,rhoparamiIterRho,GridVolume0,GridVolume1,CosphiparamiIterPhi,SinphiparamiIterPhi):
MagneticScalarPotentialOriginC=np.zeros((lenGridVolume0,lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticScalarPotentialOriginCcount=np.zeros((lenGridVolume0,lenGridVolume1,lenGridVolume2),dtype=np.float32)
OnesZ=np.ones((lenGridVolume2),dtype=np.float32)
for iIterPhi in range(0,len(phiparamInterpolated),1):
Xindex=np.argmin(np.abs(GridVolume0-rhoparamiIterRho*CosphiparamiIterPhi[iIterPhi]))
Yindex=np.argmin(np.abs(GridVolume1-rhoparamiIterRho*SinphiparamiIterPhi[iIterPhi]))
MagneticScalarPotentialOriginC[Xindex,Yindex]+=MagneticScalarPotentialOriginInterpolated[iIterPhi]
MagneticScalarPotentialOriginCcount[Xindex,Yindex]+=OnesZ
return MagneticScalarPotentialOriginC,MagneticScalarPotentialOriginCcount
def VectorConeCylCoordCalcIntegralsProc_star(a_b):
"""Convert `f([1,2])` to `f(1,2)` call."""
return VectorConeCylCoordCalcIntegralsProc(*a_b)
def VectorConeCylCoordCalcIntegralsProc(iIterRho):
[RhoZparamiIterTop,CosPhiparamTopPhiparamiIterAux1Top,zparamiIterZTop,RhoZparamiIterBottom,CosPhiparamBottomPhiparamiIterAux1Bottom,zparamiIterZBottom,RhoZparamiIterSide,RhoZparamiIterSideDiff,CosPhiparamSidePhiparamiIterAux1Side,zparamiIterZSide,rhoparam,SumdPhi,SumdZ,DimensionsParamCone,VolumeCurrenttopP,VolumeCurrentbottomP,VolumeCurrentsideP,VolumeCurrentInsideP,PositionParam2]=InitData.initializerAux_data
# https://ieeexplore.ieee.org/document/9085339
#lenphiparam=len(phiparam)
#lenzparam=len(zparam)
#lenzparamiIterAux2=len(zparamiIterAux2)
tanpi2DimensionsParamCone2=np.tan(np.pi/2.0-DimensionsParamCone[2])
rhoparamiIterRho=rhoparam[iIterRho]
rhoparamiIterRhoSquare=rhoparamiIterRho**2
# The below line is currently computed in the main script to avoid repeating it multiple times. If discommented, then it would be computed here.
#[VolumeCurrentsideP,VolumeCurrentInsideP]=NJITvectorlenzparamiIterAux2For(lenzparamiIterAux2,lenphiparam,lenzparam,nsideP,nInsidesideP,Lindex0,phiparam,rhoparam,zparamiIterAux2,DimensionsParamCone,MagnetizationParamRhoP,MagnetizationParamPhiP,MagnetizationParamZP,zparam)
SumdRho=rhoparam[1]-rhoparam[0]
#SumdPhi=phiparam[1]-phiparam[0]
#SumdZ=zparam[1]-zparam[0]
MinGridCyl=(SumdRho**2+(SumdPhi*SumdRho)**2+SumdZ**2)/64.0 # Without np.sqrt
TopBottomSurfaceDif=SumdRho*SumdPhi# The rho multiplicaton is done in the Double SumProc
# Top
## The below line is currently computed in the main script to avoid repeating it multiple times. If discommented, then it would be computed here.
##[RhoZparamiIterTop,PhiparamiIterAux1Top,PhiparamTop,zparamiIterZTop]=NumbaMeshgrid4D(rhoparamiIterAux2rindex,phiparam,phiparam,zparam)
#DiffDistance=rhoparamiIterRhoSquare+RhoZparamiIterTop**2+(zparamiIterZTop-(PositionParam2+DimensionsParamCone[1]))**2-2.0*rhoparamiIterRho*RhoZparamiIterTop*CosPhiparamTopPhiparamiIterAux1Top
#DiffDistanceLogic=(DiffDistance<MinGridCyl)
#if (np.sum(DiffDistanceLogic)>0):
# DiffDistance[DiffDistanceLogic]=MinGridCyl#/1024.0 # Without np.sqrt
#SumCommonFactor=RhoZparamiIterTop/np.sqrt(DiffDistance)
## Asuming there are no 0 distances
SumCommonFactor=np.nan_to_num(RhoZparamiIterTop/np.sqrt(rhoparamiIterRhoSquare+RhoZparamiIterTop**2+(zparamiIterZTop-(PositionParam2+DimensionsParamCone[1]))**2-2.0*rhoparamiIterRho*RhoZparamiIterTop*CosPhiparamTopPhiparamiIterAux1Top), copy=False, nan=0.0, posinf=0.0, neginf=0.0)
DoubleIntegralTermRho=DoubleSumProc4D(SumCommonFactor*VolumeCurrenttopP[0],TopBottomSurfaceDif)
DoubleIntegralTermPhi=DoubleSumProc4D(SumCommonFactor*VolumeCurrenttopP[1],TopBottomSurfaceDif)
DoubleIntegralTermZ=DoubleSumProc4D(SumCommonFactor*VolumeCurrenttopP[2],TopBottomSurfaceDif)
# Bottom
## The below line is currently computed in the main script to avoid repeating it multiple times. If discommented, then it would be computed here.
##[RhoZparamiIterBottom,PhiparamiIterAux1Bottom,PhiparamBottom,zparamiIterZBottom]=NumbaMeshgrid4D(rhoparamiIterAux2Rindex,phiparam,phiparam,zparam)
#DiffDistance=rhoparamiIterRhoSquare+RhoZparamiIterBottom**2+(zparamiIterZBottom-PositionParam2)**2-2.0*rhoparamiIterRho*RhoZparamiIterBottom*CosPhiparamBottomPhiparamiIterAux1Bottom
#DiffDistanceLogic=(DiffDistance<MinGridCyl)
#if (np.sum(DiffDistanceLogic)>0):
# DiffDistance[DiffDistanceLogic]=MinGridCyl#/1024.0 # Without np.sqrt
#SumCommonFactor=RhoZparamiIterBottom/np.sqrt(DiffDistance)
## Asuming there are no 0 distances
SumCommonFactor=np.nan_to_num(RhoZparamiIterBottom/np.sqrt(rhoparamiIterRhoSquare+RhoZparamiIterBottom**2+(zparamiIterZBottom-PositionParam2)**2-2.0*rhoparamiIterRho*RhoZparamiIterBottom*CosPhiparamBottomPhiparamiIterAux1Bottom), copy=False, nan=0.0, posinf=0.0, neginf=0.0)
DoubleIntegralTermRho+=DoubleSumProc4D(SumCommonFactor*VolumeCurrentbottomP[0],TopBottomSurfaceDif)
DoubleIntegralTermPhi+=DoubleSumProc4D(SumCommonFactor*VolumeCurrentbottomP[1],TopBottomSurfaceDif)
DoubleIntegralTermZ+=DoubleSumProc4D(SumCommonFactor*VolumeCurrentbottomP[2],TopBottomSurfaceDif)
SideInSideSurfaceDif=SumdPhi*SumdZ# The rho multiplicaton is done in the Double SumProc
# Side
## The below line is currently computed in the main script to avoid repeating it multiple times. If discommented, then it would be computed here.
##[RhoZparamiIterSide,PhiparamiIterAux1Side,PhiparamSide,zparamiIterZSide]=NumbaMeshgrid4D(zparamiIterAux2,phiparam,phiparam,zparam)
DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2=DimensionsParamCone[0]-RhoZparamiIterSideDiff*tanpi2DimensionsParamCone2
#DiffDistance=rhoparamiIterRhoSquare+DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2**2+(zparamiIterZSide-RhoZparamiIterSide)**2-2.0*rhoparamiIterRho*DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2*CosPhiparamSidePhiparamiIterAux1Side
#DiffDistanceLogic=(DiffDistance<MinGridCyl)
#if (np.sum(DiffDistanceLogic)>0):
# DiffDistance[DiffDistanceLogic]=MinGridCyl#/1024.0 # Without np.sqrt
#SumCommonFactor=DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2/np.sqrt(DiffDistance)
## Asuming there are no 0 distances
SumCommonFactor=np.nan_to_num(DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2/np.sqrt(rhoparamiIterRhoSquare+DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2**2+(zparamiIterZSide-RhoZparamiIterSide)**2-2.0*rhoparamiIterRho*DimensionsParamCone0ZparamiIterAux2tanpi2DimensionsParamCone2*CosPhiparamSidePhiparamiIterAux1Side), copy=False, nan=0.0, posinf=0.0, neginf=0.0)
DoubleIntegralTermRho+=DoubleSumProc4D(SumCommonFactor*VolumeCurrentsideP[0],SideInSideSurfaceDif)
DoubleIntegralTermPhi+=DoubleSumProc4D(SumCommonFactor*VolumeCurrentsideP[1],SideInSideSurfaceDif)
DoubleIntegralTermZ+=DoubleSumProc4D(SumCommonFactor*VolumeCurrentsideP[2],SideInSideSurfaceDif)
# Inside
if (DimensionsParamCone[3]>0.0):
DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2=DimensionsParamCone[3]-RhoZparamiIterSideDiff*tanpi2DimensionsParamCone2
#DiffDistance=rhoparamiIterRhoSquare+DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2**2+(zparamiIterZSide-RhoZparamiIterSide)**2-2.0*rhoparamiIterRho*DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2*CosPhiparamSidePhiparamiIterAux1Side
#DiffDistanceLogic=(DiffDistance<MinGridCyl)
#if (np.sum(DiffDistanceLogic)>0):
# DiffDistance[DiffDistanceLogic]=MinGridCyl#/1024.0 # Without np.sqrt
#SumCommonFactor=DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2/np.sqrt(DiffDistance)
## Asuming there are no 0 distances
SumCommonFactor=np.nan_to_num(DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2/np.sqrt(rhoparamiIterRhoSquare+DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2**2+(zparamiIterZSide-RhoZparamiIterSide)**2-2.0*rhoparamiIterRho*DimensionsParamCone3ZparamiIterAux2tanpi2DimensionsParamCone2*CosPhiparamSidePhiparamiIterAux1Side), copy=False, nan=0.0, posinf=0.0, neginf=0.0)
DoubleIntegralTermRho-=DoubleSumProc4D(SumCommonFactor*VolumeCurrentInsideP[0],SideInSideSurfaceDif)
DoubleIntegralTermPhi-=DoubleSumProc4D(SumCommonFactor*VolumeCurrentInsideP[1],SideInSideSurfaceDif)
DoubleIntegralTermZ-=DoubleSumProc4D(SumCommonFactor*VolumeCurrentInsideP[2],SideInSideSurfaceDif)
return DoubleIntegralTermRho,DoubleIntegralTermPhi,DoubleIntegralTermZ
@jit
def NJITvectorVolumesFor(rhoparam,phiparam,zparam,MagnetizationParamRhoP,MagnetizationParamPhiP,MagnetizationParamZP,zparamiIterAux2,DimensionsParamCone,ntopP,nbottomP,nsideP,nInsidesideP,Lindex0,Lindex,rindex0,Rindex0,lenrhoparamiIterAux2rindex,lenrhoparamiIterAux2Rindex):
dRhoPart=0.0#(rhoparam[1]-rhoparam[0])/2.0
lenphiparam=len(phiparam)
lenzparam=len(zparam)
lenzparamiIterAux2=len(zparamiIterAux2)
VolumeCurrenttopP=np.empty((3,lenrhoparamiIterAux2rindex,lenphiparam,lenphiparam,lenzparam),dtype=np.float32)
for iIterAux in range(0,lenrhoparamiIterAux2rindex,1):
rindex0iIterAux=rindex0+iIterAux
VolumeCurrenttopP[0][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamPhiP[rindex0iIterAux,:,Lindex]*ntopP[2]-MagnetizationParamZP[rindex0iIterAux,:,Lindex]*ntopP[1],phiparam,zparam)
VolumeCurrenttopP[1][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamZP[rindex0iIterAux,:,Lindex]*ntopP[0]-MagnetizationParamRhoP[rindex0iIterAux,:,Lindex]*ntopP[2],phiparam,zparam)
VolumeCurrenttopP[2][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamRhoP[rindex0iIterAux,:,Lindex]*ntopP[1]-MagnetizationParamPhiP[rindex0iIterAux,:,Lindex]*ntopP[0],phiparam,zparam)
VolumeCurrentbottomP=np.empty((3,lenrhoparamiIterAux2Rindex,lenphiparam,lenphiparam,lenzparam),dtype=np.float32)
for iIterAux in range(0,lenrhoparamiIterAux2Rindex,1):
Rindex0iIterAux=Rindex0+iIterAux
VolumeCurrentbottomP[0][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamPhiP[Rindex0iIterAux,:,Lindex0]*nbottomP[2]-MagnetizationParamZP[Rindex0iIterAux,:,Lindex0]*nbottomP[1],phiparam,zparam)
VolumeCurrentbottomP[1][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamZP[Rindex0iIterAux,:,Lindex0]*nbottomP[0]-MagnetizationParamRhoP[Rindex0iIterAux,:,Lindex0]*nbottomP[2],phiparam,zparam)
VolumeCurrentbottomP[2][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamRhoP[Rindex0iIterAux,:,Lindex0]*nbottomP[1]-MagnetizationParamPhiP[Rindex0iIterAux,:,Lindex0]*nbottomP[0],phiparam,zparam)
VolumeCurrentsideP=np.empty((3,lenzparamiIterAux2,lenphiparam,lenphiparam,lenzparam),dtype=np.float32)
VolumeCurrentInsideP=np.empty((3,lenzparamiIterAux2,lenphiparam,lenphiparam,lenzparam),dtype=np.float32)
tanpi2DimensionsParamCone2=np.tan(np.pi/2.0-DimensionsParamCone[2])
for iIterAux in range(0,lenzparamiIterAux2,1):
LindexCommon=Lindex0+iIterAux
IndexCalCommon=np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-(zparamiIterAux2[iIterAux]-zparamiIterAux2[0])*tanpi2DimensionsParamCone2)))
VolumeCurrentsideP[0][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamPhiP[IndexCalCommon,:,LindexCommon]*nsideP[2]-MagnetizationParamZP[IndexCalCommon,:,LindexCommon]*nsideP[1],phiparam,zparam)
VolumeCurrentsideP[1][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamZP[IndexCalCommon,:,LindexCommon]*nsideP[0]-MagnetizationParamRhoP[IndexCalCommon,:,LindexCommon]*nsideP[2],phiparam,zparam)
VolumeCurrentsideP[2][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamRhoP[IndexCalCommon,:,LindexCommon]*nsideP[1]-MagnetizationParamPhiP[IndexCalCommon,:,LindexCommon]*nsideP[0],phiparam,zparam)
# inside
IndexCalCommon=np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-(zparamiIterAux2[iIterAux]-zparamiIterAux2[0])*tanpi2DimensionsParamCone2)))
VolumeCurrentInsideP[0][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamPhiP[IndexCalCommon,:,LindexCommon]*nInsidesideP[2]-MagnetizationParamZP[IndexCalCommon,:,LindexCommon]*nInsidesideP[1],phiparam,zparam)
VolumeCurrentInsideP[1][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamZP[IndexCalCommon,:,LindexCommon]*nInsidesideP[0]-MagnetizationParamRhoP[IndexCalCommon,:,LindexCommon]*nInsidesideP[2],phiparam,zparam)
VolumeCurrentInsideP[2][iIterAux]=NumbaMeshgrid3D1axis(MagnetizationParamRhoP[IndexCalCommon,:,LindexCommon]*nInsidesideP[1]-MagnetizationParamPhiP[IndexCalCommon,:,LindexCommon]*nInsidesideP[0],phiparam,zparam)
return VolumeCurrenttopP,VolumeCurrentbottomP,VolumeCurrentsideP,VolumeCurrentInsideP
def VectorConeCylCoordMagVectorPotCtoCylProc_star(a_b):
"""Convert `f([1,2])` to `f(1,2)` call."""
return VectorConeCylCoordMagVectorPotCtoCylProc(*a_b)
def VectorConeCylCoordMagVectorPotCtoCylProc(iIterRho):
[rhoparam,phiparam,GridVolume,MagneticVectorPotentialOriginRho,MagneticVectorPotentialOriginPhi,MagneticVectorPotentialOriginZ]=InitData.initializerAux_data
dRhoPart=0.0#(rhoparam[1]-rhoparam[0])/2.0
lenGridVolume0=len(GridVolume[0])
lenGridVolume1=len(GridVolume[1])
lenGridVolume2=len(GridVolume[2])
#for iIterRho in range(0,len(rhoparam),1):
# Interpolate Phi axis to have more points and populate X and Y
InterpolationOrder=3
InterpolationFactor=3
rhoparamiIterRho=rhoparam[iIterRho]+dRhoPart
phiparamInterpolated=(ndimage.zoom(phiparam,InterpolationFactor,mode = 'wrap',order=InterpolationOrder,grid_mode=False)).astype(np.float32)
MagneticVectorPotentialOriginRhoInterpolated=(1.0/float(InterpolationFactor))*(ndimage.zoom(MagneticVectorPotentialOriginRho[iIterRho],[InterpolationFactor,1],mode = 'wrap',order=InterpolationOrder,grid_mode=False)).astype(np.float32)
MagneticVectorPotentialOriginPhiInterpolated=(1.0/float(InterpolationFactor))*(ndimage.zoom(MagneticVectorPotentialOriginPhi[iIterRho],[InterpolationFactor,1],mode = 'wrap',order=InterpolationOrder,grid_mode=False)).astype(np.float32)
MagneticVectorPotentialOriginZInterpolated=(1.0/float(InterpolationFactor))*(ndimage.zoom(MagneticVectorPotentialOriginZ[iIterRho],[InterpolationFactor,1],mode = 'wrap',order=InterpolationOrder,grid_mode=False)).astype(np.float32)
CosphiparamiIterPhi=np.cos(phiparamInterpolated)
SinphiparamiIterPhi=np.sin(phiparamInterpolated)
return NJITForVectorConeCylCoordMagVectorPotCtoCylProc(phiparamInterpolated,MagneticVectorPotentialOriginRhoInterpolated,MagneticVectorPotentialOriginPhiInterpolated,MagneticVectorPotentialOriginZInterpolated,lenGridVolume0,lenGridVolume1,lenGridVolume2,GridVolume[0],GridVolume[1],CosphiparamiIterPhi,SinphiparamiIterPhi,rhoparamiIterRho)
@jit
def NJITForVectorConeCylCoordMagVectorPotCtoCylProc(phiparamInterpolated,MagneticVectorPotentialOriginRhoInterpolated,MagneticVectorPotentialOriginPhiInterpolated,MagneticVectorPotentialOriginZInterpolated,lenGridVolume0,lenGridVolume1,lenGridVolume2,GridVolume0,GridVolume1,CosphiparamiIterPhi,SinphiparamiIterPhi,rhoparamiIterRho):
MagneticVectorPotentialOriginCx=np.zeros((lenGridVolume0,lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginCy=np.zeros((lenGridVolume0,lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginCz=np.zeros((lenGridVolume0,lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginCcount=np.zeros((lenGridVolume0,lenGridVolume1,lenGridVolume2),dtype=np.float32)
OnesZ=np.ones((lenGridVolume2),dtype=np.float32)
for iIterPhi in range(0,len(phiparamInterpolated),1):
Xindex=np.argmin(np.abs(GridVolume0-rhoparamiIterRho*CosphiparamiIterPhi[iIterPhi]))
Yindex=np.argmin(np.abs(GridVolume1-rhoparamiIterRho*SinphiparamiIterPhi[iIterPhi]))
MagneticVectorPotentialOriginCx[Xindex,Yindex]+=MagneticVectorPotentialOriginRhoInterpolated[iIterPhi]*CosphiparamiIterPhi[iIterPhi]-MagneticVectorPotentialOriginPhiInterpolated[iIterPhi]*SinphiparamiIterPhi[iIterPhi]
MagneticVectorPotentialOriginCy[Xindex,Yindex]+=MagneticVectorPotentialOriginRhoInterpolated[iIterPhi]*SinphiparamiIterPhi[iIterPhi]+MagneticVectorPotentialOriginPhiInterpolated[iIterPhi]*CosphiparamiIterPhi[iIterPhi]
MagneticVectorPotentialOriginCz[Xindex,Yindex]+=MagneticVectorPotentialOriginZInterpolated[iIterPhi]
MagneticVectorPotentialOriginCcount[Xindex,Yindex]+=OnesZ
return MagneticVectorPotentialOriginCx,MagneticVectorPotentialOriginCy,MagneticVectorPotentialOriginCz,MagneticVectorPotentialOriginCcount
def ConeCylCoordMagnetizationVectorProc_star(a_b):
"""Convert `f([1,2])` to `f(1,2)` call."""
return ConeCylCoordMagnetizationVectorProc(*a_b)
def ConeCylCoordMagnetizationVectorProc(iIterRho):
[rhoparam,phiparam,zparam,GridVolume,MagnetizationParam]=InitData.initializerAux_data
dRhoPart=0.0#(rhoparam[1]-rhoparam[0])/2.0
lenrhoparam=len(rhoparam)
lenphiparam=len(phiparam)
lenzparam=len(zparam)
MagnetizationParamRhoP=np.zeros((lenrhoparam,lenphiparam,lenzparam),dtype=np.float32)
MagnetizationParamPhiP=np.zeros((lenrhoparam,lenphiparam,lenzparam),dtype=np.float32)
MagnetizationParamZP=np.zeros((lenrhoparam,lenphiparam,lenzparam),dtype=np.float32)
MagnetizationParamPcount=np.zeros((lenrhoparam,lenphiparam,lenzparam),dtype=np.float32)
#for iIterRho in range(0,len(rhoparam),1):
rhoparamiIterRho=rhoparam[iIterRho]
CosphiparamiIterPhi=np.cos(phiparam)
SinphiparamiIterPhi=np.sin(phiparam)
OnesZ=np.ones((lenzparam),dtype=np.float32)
#ZerosZ=np.zeros((lenzparam),dtype=np.float32)
for iIterPhi in range(0,lenphiparam,1):
#for iIterZ in range(0,lenzparam,1):
#Zval=zparam[iIterZ]
Xindex=np.argmin(np.abs(GridVolume[0]-CosphiparamiIterPhi[iIterPhi]*(rhoparamiIterRho+dRhoPart)))
Yindex=np.argmin(np.abs(GridVolume[1]-SinphiparamiIterPhi[iIterPhi]*(rhoparamiIterRho+dRhoPart)))
MagnetizationParamRhoP[iIterRho,iIterPhi]=(CosphiparamiIterPhi[iIterPhi]*MagnetizationParam[0][Xindex,Yindex]+SinphiparamiIterPhi[iIterPhi]*MagnetizationParam[1][Xindex,Yindex])#CosPhiP*MagnetizationParam[0][Xindex,Yindex]+SinPhiP*MagnetizationParam[1][Xindex,Yindex]#np.sqrt(MagnetizationParam[0][Xindex,Yindex]**2+MagnetizationParam[1][Xindex,Yindex]**2)#
MagnetizationParamPhiP[iIterRho,iIterPhi]=(-SinphiparamiIterPhi[iIterPhi]*MagnetizationParam[0][Xindex,Yindex]+CosphiparamiIterPhi[iIterPhi]*MagnetizationParam[1][Xindex,Yindex])#ZerosZ#(-SinphiparamiIterPhi[iIterPhi]*MagnetizationParam[0][Xindex,Yindex]+CosphiparamiIterPhi[iIterPhi]*MagnetizationParam[1][Xindex,Yindex])
MagnetizationParamZP[iIterRho,iIterPhi]=MagnetizationParam[2][Xindex,Yindex]
MagnetizationParamPcount[iIterRho,iIterPhi]=OnesZ
return MagnetizationParamRhoP,MagnetizationParamPhiP,MagnetizationParamZP,MagnetizationParamPcount
def CanonicalSurfXMagVectorPotField3DProc_star(a_b):
"""Convert `f([1,2])` to `f(1,2)` call."""
return CanonicalSurfXMagVectorPotField3DProc(*a_b)
def CanonicalSurfXMagVectorPotField3DProc(iIterX):
[GridVolume,PositionParam,Ypos,Zpos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Yn,Zn]=InitData.initializerAux_data
return NJITCanonicalSurfXMagVectorPotField3DProc(iIterX,GridVolume[0],GridVolume[1],GridVolume[2],PositionParam[0],Ypos,Zpos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Yn,Zn)
@jit
def NJITCanonicalSurfXMagVectorPotField3DProc(iIterX,GridVolume0,GridVolume1,GridVolume2,PositionParam0,Ypos,Zpos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Yn,Zn):
lenGridVolume1=len(GridVolume1)
lenGridVolume2=len(GridVolume2)
MagneticVectorPotentialOriginxaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginyaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginzaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
PermRadFactor=(Permeability_0/(4.0*np.pi))
#for iIterX in range(0,len(GridVolume0),1):
for iIterY in range(0,lenGridVolume1,1):
for iIterZ in range(0,lenGridVolume2,1):
for iIterLine in range(0,len(Ypos),1):
Distancediff=np.sqrt((GridVolume0[iIterX]-PositionParam0)**2+(GridVolume1[iIterY]-Ypos[iIterLine])**2+(GridVolume2[iIterZ]-Zpos[iIterLine])**2)
MagneticVectorPotentialOriginxaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentX/Distancediff)
MagneticVectorPotentialOriginyaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentY/Distancediff)*Yn[iIterLine]
MagneticVectorPotentialOriginzaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentZ/Distancediff)*Zn[iIterLine]
return MagneticVectorPotentialOriginxaux,MagneticVectorPotentialOriginyaux,MagneticVectorPotentialOriginzaux
def CanonicalSurfYMagVectorPotField3DProc_star(a_b):
"""Convert `f([1,2])` to `f(1,2)` call."""
return CanonicalSurfYMagVectorPotField3DProc(*a_b)
def CanonicalSurfYMagVectorPotField3DProc(iIterX):
[GridVolume,PositionParam,Xpos,Zpos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Xn,Zn]=InitData.initializerAux_data
return NJITCanonicalSurfYMagVectorPotField3DProc(iIterX,GridVolume[0],GridVolume[1],GridVolume[2],PositionParam[1],Xpos,Zpos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Xn,Zn)
@jit
def NJITCanonicalSurfYMagVectorPotField3DProc(iIterX,GridVolume0,GridVolume1,GridVolume2,PositionParam1,Xpos,Zpos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Xn,Zn):
lenGridVolume1=len(GridVolume1)
lenGridVolume2=len(GridVolume2)
MagneticVectorPotentialOriginxaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginyaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginzaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
PermRadFactor=(Permeability_0/(4.0*np.pi))
#for iIterX in range(0,len(GridVolume0),1):
for iIterY in range(0,lenGridVolume1,1):
for iIterZ in range(0,lenGridVolume2,1):
for iIterLine in range(0,len(Xpos),1):
Distancediff=np.sqrt((GridVolume0[iIterX]-Xpos[iIterLine])**2+(GridVolume1[iIterY]-PositionParam1)**2+(GridVolume2[iIterZ]-Zpos[iIterLine])**2)
MagneticVectorPotentialOriginxaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentX/Distancediff)*Xn[iIterLine]
MagneticVectorPotentialOriginyaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentY/Distancediff)
MagneticVectorPotentialOriginzaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentZ/Distancediff)*Zn[iIterLine]
return MagneticVectorPotentialOriginxaux,MagneticVectorPotentialOriginyaux,MagneticVectorPotentialOriginzaux
def CanonicalSurfZMagVectorPotField3DProc_star(a_b):
"""Convert `f([1,2])` to `f(1,2)` call."""
return CanonicalSurfZMagVectorPotField3DProc(*a_b)
def CanonicalSurfZMagVectorPotField3DProc(iIterX):
[GridVolume,PositionParam,Xpos,Ypos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Xn,Yn]=InitData.initializerAux_data
return NJITCanonicalSurfZMagVectorPotField3DProc(iIterX,GridVolume[0],GridVolume[1],GridVolume[2],PositionParam[2],Xpos,Ypos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Xn,Yn)
@jit
def NJITCanonicalSurfZMagVectorPotField3DProc(iIterX,GridVolume0,GridVolume1,GridVolume2,PositionParam2,Xpos,Ypos,Permeability_0,SurfCurrentX,SurfCurrentY,SurfCurrentZ,Xn,Yn):
lenGridVolume1=len(GridVolume1)
lenGridVolume2=len(GridVolume2)
MagneticVectorPotentialOriginxaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginyaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
MagneticVectorPotentialOriginzaux=np.zeros((lenGridVolume1,lenGridVolume2),dtype=np.float32)
PermRadFactor=(Permeability_0/(4.0*np.pi))
#for iIterX in range(0,len(GridVolume0),1):
for iIterY in range(0,lenGridVolume1,1):
for iIterZ in range(0,lenGridVolume2,1):
for iIterLine in range(0,len(Xpos),1):
Distancediff=np.sqrt((GridVolume0[iIterX]-Xpos[iIterLine])**2+(GridVolume1[iIterY]-Ypos[iIterLine])**2+(GridVolume2[iIterZ]-PositionParam2)**2)
MagneticVectorPotentialOriginxaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentX/Distancediff)*Xn[iIterLine]
MagneticVectorPotentialOriginyaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentY/Distancediff)*Yn[iIterLine]
MagneticVectorPotentialOriginzaux[iIterY,iIterZ]+=PermRadFactor*(SurfCurrentZ/Distancediff)
return MagneticVectorPotentialOriginxaux,MagneticVectorPotentialOriginyaux,MagneticVectorPotentialOriginzaux
@jit
def NJITlogicConeVolumeIndices(rhoparam,zparam,DimensionsParamCone,MagnetizationParamRhoPshape,Lindex,Lindex0):
dRhoPart=0.0#(rhoparam[1]-rhoparam[0])/2.0
# Logic cone volume indices (insiders, not counting the surface ones)
LogicConeVolumeIndices=np.zeros((MagnetizationParamRhoPshape),dtype=np.bool_)
LogicConeVolumeIndicesAux=np.zeros((MagnetizationParamRhoPshape),dtype=np.bool_)
LogicConeVolumeIndicesTop=np.zeros((MagnetizationParamRhoPshape),dtype=np.bool_)
LogicConeVolumeIndicesTopAux=np.zeros((MagnetizationParamRhoPshape),dtype=np.bool_)
LogicConeVolumeIndicesBottom=np.zeros((MagnetizationParamRhoPshape),dtype=np.bool_)
LogicConeVolumeIndicesBottomAux=np.zeros((MagnetizationParamRhoPshape),dtype=np.bool_)
if (DimensionsParamCone[3]>0.0):
for iIterConeVolume in range(0,Lindex-Lindex0+1,1):
zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2=(zparam[Lindex0+iIterConeVolume]-zparam[Lindex0])*np.tan(np.pi/2.0-DimensionsParamCone[2])
LogicConeVolumeIndices[0:np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True # Where average or maximum is computed from
LogicConeVolumeIndicesAux[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))):np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True # Where average or maximum is applied
if (iIterConeVolume==(Lindex-Lindex0-1)):
LogicConeVolumeIndicesTop[0:np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True # Where average or maximum is computed from
LogicConeVolumeIndicesTopAux[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))):np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True # Where average or maximum is applied
if (iIterConeVolume==0):
LogicConeVolumeIndicesBottom[0:np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True # Where average or maximum is computed from
LogicConeVolumeIndicesBottomAux[np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[3]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))):np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True # Where average or maximum is applied
else:
for iIterConeVolume in range(0,Lindex-Lindex0+1,1):
zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2=(zparam[Lindex0+iIterConeVolume]-zparam[Lindex0])*np.tan(np.pi/2.0-DimensionsParamCone[2])
LogicConeVolumeIndices[0:np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True
if (iIterConeVolume==(Lindex-Lindex0-1)):
LogicConeVolumeIndicesTop[0:np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True
if (iIterConeVolume==0):
LogicConeVolumeIndicesBottom[0:np.argmin(np.abs(rhoparam-(dRhoPart+DimensionsParamCone[0]-zparamLindex0iIterConeVolumezparamLindex0tanpi2DimensionsParamCone2))),:,Lindex0+iIterConeVolume]=True
LogicConeVolumeIndicesAux=LogicConeVolumeIndices
LogicConeVolumeIndicesTopAux=LogicConeVolumeIndicesTop
LogicConeVolumeIndicesBottomAux=LogicConeVolumeIndicesBottom
return LogicConeVolumeIndices,LogicConeVolumeIndicesAux,LogicConeVolumeIndicesTop,LogicConeVolumeIndicesTopAux,LogicConeVolumeIndicesBottom,LogicConeVolumeIndicesBottomAux
def NJITCartesianlogicConeVolumeIndices(GridVolume,DimensionsParamCone,PositionParam1):
# Logic cone volume indices
dx=np.abs(GridVolume[0][1]-GridVolume[0][0])
dy=np.abs(GridVolume[1][1]-GridVolume[1][0])
dz=np.abs(GridVolume[2][1]-GridVolume[2][0])
LogicConeVolumeIndices=np.zeros((len(GridVolume[0]),len(GridVolume[1]),len(GridVolume[2])),dtype=np.bool_)
LogicConeVolumeIndicesAux=np.zeros((len(GridVolume[0]),len(GridVolume[1]),len(GridVolume[2])),dtype=np.bool_)
LogicConeVolumeIndicesTop=np.zeros((len(GridVolume[0]),len(GridVolume[1]),len(GridVolume[2])),dtype=np.bool_)
LogicConeVolumeIndicesTopAux=np.zeros((len(GridVolume[0]),len(GridVolume[1]),len(GridVolume[2])),dtype=np.bool_)
LogicConeVolumeIndicesBottom=np.zeros((len(GridVolume[0]),len(GridVolume[1]),len(GridVolume[2])),dtype=np.bool_)
LogicConeVolumeIndicesBottomAux=np.zeros((len(GridVolume[0]),len(GridVolume[1]),len(GridVolume[2])),dtype=np.bool_)
# XY axis
[iIterXaux,iIterYaux]=NumbaMeshgrid2D(GridVolume[0]-PositionParam1[0],GridVolume[1]-PositionParam1[1])
# Z axis
Lindex=np.argmin(np.abs(GridVolume[2]-np.floor((PositionParam1[2]+DimensionsParamCone[1])/dz)*dz))
Lindex0=np.argmin(np.abs(GridVolume[2]-np.floor((PositionParam1[2])/dz)*dz))
if (DimensionsParamCone[3]>0.0):
for iIterZ in range(0,Lindex-Lindex0+1,1):
dziIterZtannppi2DimensionsParamCone2=(dz*iIterZ)*np.tan(np.pi/2.0-DimensionsParamCone[2])
iIterXYlogic=np.array((iIterXaux**2+iIterYaux**2)<=(DimensionsParamCone[0]-dziIterZtannppi2DimensionsParamCone2)**2,dtype=np.bool_)
LogicConeVolumeIndices[iIterXYlogic,Lindex0+iIterZ]=True
LogicConeVolumeIndicesAux[iIterXYlogic,Lindex0+iIterZ]=True
if (iIterZ==0):
LogicConeVolumeIndicesBottom[iIterXYlogic,Lindex0+iIterZ]=True
LogicConeVolumeIndicesBottomAux[iIterXYlogic,Lindex0+iIterZ]=True
if (iIterZ==(Lindex-Lindex0-1)):
LogicConeVolumeIndicesTop[iIterXYlogic,Lindex0+iIterZ]=True
LogicConeVolumeIndicesTopAux[iIterXYlogic,Lindex0+iIterZ]=True
iIterXYlogic=np.array((iIterXaux**2+iIterYaux**2)<=(DimensionsParamCone[3]-dziIterZtannppi2DimensionsParamCone2)**2,dtype=np.bool_)
LogicConeVolumeIndicesAux[iIterXYlogic,Lindex0+iIterZ]=False
if (iIterZ==0):
LogicConeVolumeIndicesBottomAux[iIterXYlogic,Lindex0+iIterZ]=False
if (iIterZ==(Lindex-Lindex0)):
LogicConeVolumeIndicesTopAux[iIterXYlogic,Lindex0+iIterZ]=False
else:
for iIterZ in range(0,Lindex-Lindex0+1,1):
dziIterZtannppi2DimensionsParamCone2=(dz*iIterZ)*np.tan(np.pi/2.0-DimensionsParamCone[2])
iIterXYlogic=np.array((iIterXaux**2+iIterYaux**2)<=(DimensionsParamCone[0]-dziIterZtannppi2DimensionsParamCone2)**2,dtype=np.bool_)
LogicConeVolumeIndices[iIterXYlogic,Lindex0+iIterZ]=True
if (iIterZ==0):
LogicConeVolumeIndicesBottom[iIterXYlogic,Lindex0+iIterZ]=True
if (iIterZ==(Lindex-Lindex0)):
LogicConeVolumeIndicesTop[iIterXYlogic,Lindex0+iIterZ]=True
LogicConeVolumeIndicesAux=LogicConeVolumeIndices
LogicConeVolumeIndicesTopAux=LogicConeVolumeIndicesTop
LogicConeVolumeIndicesBottomAux=LogicConeVolumeIndicesBottom
return LogicConeVolumeIndices,LogicConeVolumeIndicesAux,LogicConeVolumeIndicesTop,LogicConeVolumeIndicesTopAux,LogicConeVolumeIndicesBottom,LogicConeVolumeIndicesBottomAux