MagneticFieldSimulation_version4.py

"""
========================================
======Magnetic field simulation=========
========================================

Created by R. van Manen and William Rijksen
Version 4

Magnetic field simulation is based on:
Modeling a cylindrical permanent magnet with a surface charge of magnetic monopoles
(http://jimhawley.ca/downloads/Permanent_Magnet_With_Magnetic_Monopole_Surface_Charge.pdf)

Structure of this python script:
    (1) Initialization
        (a) loading libraries
        (b) closing plot windows
    (2) Constants
        (a) sensor related constants
        (b) magnet related constants
        (c) axis limits of the field plot
        (d) resolutions
        (e) trajectory
    (3) Variables
    (4) Pre processing
        (a) trajectory start and end points
        (b) sensor orientation vector format
        (c) sensor limits
    (5) Function definitions
        (a) magnetic Field Simulation
        (b) trajectory length
        (c) all inclusive length function
    (6) Processing
        (a) field calculations
        (b) trajectory calculations
        (c) brute force trajectory length
    (7) Post processing
        (a) plot window
        (b) field plot
        (c) trajectory plot

"""

"============(1) Initialization==============================================="
#Loading libraries
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from scipy import integrate
from sys import platform

#closing plot windows
plt.close("all")

"============(2) Constants===================================================="
#sensor related constants
S_sat       = 37.5               #Sensor saturation [mT]
S_sen       = 40                 #Sensor sensitivity [mV/mT]
S_vrf       = 3.3                #Sensor reference voltage [V]
ADC_res     = 12                 #ADconverter resolution [Bit]
S_res       = 0.5                #Target resolution [mm/bit]
#magnet related constants
M_h         = 5e-3               #Magnet height [m]
M_r         = 2.5e-3             #Magnet radius [m]
Br          = 1.45               #Remanence field [T] 
#axis limits of the field plot
F_Zmax      = 50e-3              #[m]
F_Zmin      = -50e-3             #[m]
F_Rmax      = 40e-3              #[m]
F_Rmin      = -40e-3             #[m]
#resolutions
F_res       = 25                 #Resolution of Field plot [-]
T_res       = 100                #Resolution of Trajectory plot [-]
#trajectory
T_L         = 40e-3              #Trajectory length [m]

"============(3) Variables===================================================="
T_dis       = 0e-3               #Trajectory distance from center [m]
M_or        = 0/180*np.pi      #Trajectory orientation [rad]
S_OR        = 0/180*np.pi

"============(4) Pre processing==============================================="
#trajectory start and end points
T_R0        = T_dis*np.cos(M_or)                    #[m]
T_R1        = T_dis*np.cos(M_or)+np.sin(M_or)*T_L   #[m]
T_Z0        = -T_dis*np.sin(M_or)                   #[m]
T_Z1        = -T_dis*np.sin(M_or)+np.cos(M_or)*T_L  #[m]
T_lin       = np.linspace(0,T_L,T_res)
#sensor orientation vector format
M_OR        = [np.sin(M_or),np.cos(M_or)]

#sensor limits
B_max       = S_sat/10**3                           #Maximum magnetic field [T]
B_min       = -S_sat/10**3                          #Minimum magnetic field [T]
dB_min      = S_vrf*10**3/(2**ADC_res*S_res*S_sen)  #Minimal magnetic field gradient [T/m]
dB_max      = -dB_min                               #Maximum magnetic field gradient [T/m]

"============(5) Function definitions========================================="
#magnetic Field Function
def B_r(z,f,P_z,P_r):
    return Br*M_r/(4*np.pi)*((P_z-z)*np.cos(f))/(P_r**2+M_r**2+(P_z-z)**2-2*M_r*P_r*np.cos(f))**(3/2)
def B_z(z,f,P_z,P_r):
    return Br*M_r/(4*np.pi)*(M_r-P_r*np.cos(f))/(P_r**2+M_r**2+(P_z-z)**2-2*M_r*P_r*np.cos(f))**(3/2)
def lb(z):
    return -M_h/2
def ub(z):
    return M_h/2
def MFF(P_z,P_r):
    (R,abserr) = integrate.dblquad(B_r,0,2*np.pi,lb,ub,args=(P_z,P_r))
    (Z,abserr) = integrate.dblquad(B_z,0,2*np.pi,lb,ub,args=(P_z,P_r))
    
    return (R,Z)
MFF         = np.vectorize(MFF)

#trajectory length
def Length(B,dB,Bmin,Bmax,dBmin,dBmax):
    T_linR = np.linspace(0,T_L,T_res*15)
    B      = np.interp(T_linR,T_lin,B)
    dB     = np.interp(T_linR,T_lin,dB)
    vector = np.invert(np.all([dB<dBmin,dB>dBmax],axis=0))
    vector = np.all([vector,B<Bmax,B>Bmin],axis=0)
    temp   = np.sign(dB)
    vector[np.where(temp[1:]!=temp[:-1])]=0
    
    T_length = 0
    while any(vector):
        T_length += 1
        vector = [vector[i] and vector[i] == vector[i+1] for i in range(len(vector)-1)]

    T_length = T_length*T_L/(T_res*15-1)
    return T_length
#Length      = np.vectorize(Length)

#all inclusive length function
def BF_length(M_or):
    T_R0      = T_dis*np.cos(M_or)                    #[m]
    T_R1      = T_dis*np.cos(M_or)+np.sin(M_or)*T_L   #[m]
    T_Z0        = -T_dis*np.sin(M_or)                   #[m]
    T_Z1        = -T_dis*np.sin(M_or)+np.cos(M_or)*T_L  #[m]
    M_OR      = [np.sin(M_or),np.cos(M_or)]
    T_Zlin    = np.linspace(T_Z0,T_Z1,T_res)
    T_Rlin    = np.linspace(T_R0,T_R1,T_res)
    T_B       = MFF(T_Zlin,T_Rlin)
    T_Bprj    = np.dot(M_OR,T_B)*np.cos(S_OR)
    T_dBprj   = np.gradient(T_Bprj,T_lin[1])
    T_l       = Length(T_Bprj,T_dBprj,B_min,B_max,dB_min,dB_max)
    return T_l
BF_length   = np.vectorize(BF_length)
#Optimization function
def OPT_length(x):
    M_or      = x[0]
    T_R0      = T_dis*np.cos(M_or)                    #[m]
    T_R1      = T_dis*np.cos(M_or)+np.sin(M_or)*T_L   #[m]
    T_Z0        = -T_dis*np.sin(M_or)                   #[m]
    T_Z1        = -T_dis*np.sin(M_or)+np.cos(M_or)*T_L  #[m]
    M_OR      = [np.sin(M_or),np.cos(M_or)]
    T_Zlin    = np.linspace(T_Z0,T_Z1,T_res)
    T_Rlin    = np.linspace(T_R0,T_R1,T_res)
    T_B       = MFF(T_Zlin,T_Rlin)
    T_Bprj    = np.dot(M_OR,T_B)*np.cos(S_OR)
    T_dBprj   = np.gradient(T_Bprj,T_lin[1])
    T_l       = 0.1/(Length(T_Bprj,T_dBprj,B_min,B_max,dB_min,dB_max)+1e-9)
    return T_l

"============(6) Processing==================================================="
#field calculations
F_Zlin    = np.linspace(F_Zmin,F_Zmax,F_res)
F_Rlin    = np.linspace(F_Rmin,F_Rmax,F_res)

F_B       = MFF(F_Zlin[:,None], F_Rlin[None,:])
F_Babs    = (F_B[0]**2+F_B[1]**2)**(1/2)

#Optimization
x0        = 0
#OPT       = sp.optimize.minimize(OPT_length,x0,method='Powell',
#                                 options={'disp': True})
#M_or      = OPT.x*1
T_R0      = T_dis*np.cos(M_or)                    #[m]
T_R1      = T_dis*np.cos(M_or)+np.sin(M_or)*T_L   #[m]
T_Z0        = -T_dis*np.sin(M_or)                   #[m]
T_Z1        = -T_dis*np.sin(M_or)+np.cos(M_or)*T_L  #[m]
M_OR      = [np.sin(M_or),np.cos(M_or)]
T_lin     = np.linspace(0,T_L,T_res)

#trajectory calculations
T_Zlin    = np.linspace(T_Z0,T_Z1,T_res)
T_Rlin    = np.linspace(T_R0,T_R1,T_res)

T_B       = MFF(T_Zlin,T_Rlin)
T_Babs    = (T_B[0]**2+T_B[1]**2)**(1/2)
T_Bprj    = np.dot(M_OR,T_B)*np.cos(S_OR)
T_dBprj   = np.gradient(T_Bprj,T_lin[1])

T_l       = Length(T_Bprj,T_dBprj,B_min,B_max,dB_min,dB_max)

"============(7) Post processing=============================================="
#plot window
PW        = plt.figure()
if platform == "win32":
    mngr      = plt.get_current_fig_manager()
    mngr.window.setGeometry(-1700,250,1500, 700)

#Field plot
Field     = PW.add_subplot(121, aspect='equal')
Field.add_patch(patches.Rectangle((-M_r*1e3, -M_h/2*1e3),M_r*2*1e3,M_h*1e3))    #magnet
plt.contourf(F_Rlin*1e3,F_Zlin*1e3,np.log(F_Babs),100,cmap=plt.cm.hot_r)        #magnetic field strenght color
plt.plot([T_R0*1e3,T_R1*1e3],[T_Z0*1e3,T_Z1*1e3], color='b',linewidth=3)        #trajectory

plt.streamplot(F_Rlin*1e3,F_Zlin*1e3,F_B[0],F_B[1],color='k',linewidth=0.2)     #magnetic field lines
plt.axis((F_Rmin*1e3,F_Rmax*1e3,F_Zmin*1e3,F_Zmax*1e3))
plt.title('Magnetic Field [M_or = {0: .4f}]'.format(M_or))
plt.xlabel('r [mm]')
plt.ylabel('z [mm]')

#Trajectory plot
Traject = PW.add_subplot(122)
plt.plot(-10,0,'ko',label='Length = %.1f [mm]'%(T_l*1e3))                       #length legend item
plt.plot(T_lin*1e3,T_Babs*1e3, label='|M| [mT]')                                #abs(M)
plt.plot(T_lin*1e3,T_Bprj*1e3, label='Projected |M| [mT]')                      #abs(Mprj)
Traject.fill_between(T_lin*1e3,-100,100,
                     where= np.logical_and(T_dBprj < dB_min,T_dBprj > dB_max),
                     facecolor='grey', label='unsensitive')                     #unsensitive regions
plt.title('Magnetic Field @ sensor path')
plt.xlabel('x [mm]')
plt.ylabel('|M| [mT]')
plt.axis((0,T_lin[-1]*1e3,B_min*1e3,B_max*1e3))
plt.legend()
plt.show()