In [ ]:

    

import numpy as np
import matplotlib.pyplot as plt

r1 = 5.0E20
r2 = 1.5E21
bf = 1.5E21
rout = r2+bf
rpol = rout + 2.0*r2
bz = 1.86E-6
PI = np.pi


    

In [ ]:

    

def Aphi(r):
    if r < r2:
        return 0.5*bz*r
    # Current Setting
    if False:
        const = r2*(0.5*r2-4.0/PI**2*bf**2/r2)*bz
        if r < rpol:
            return (bz*2.0/PI**2*bf*(2*bf/r*np.sin(PI/2/bf*(rout-r))+PI*np.cos(PI/2/bf*(rout-r))) + const/r)*\
                    0.5*(1.0+cos(PI*(max(min(r,rpol),rout)-rout)/(rpol-rout)))
        else:
            return 0.0
        
    if False:
        const = r2*(0.5*r2-4.0/PI**2*bf**2/r2)*bz
        if r < rout+2*bf:
            return bz*2.0/PI**2*bf*(2*bf/r*np.sin(PI/2/bf*(rout-r))+PI*np.cos(PI/2/bf*(rout-r))) + const/r
        else:
            return (-bz/PI*2.0*bf*(rout+2*bf)+const)/r

    if False:
        const = r2*(0.5*r2-4.0/PI**2*bf**2/r2)*bz
        if r < rpol:
            return 0.5*bz*r2*(1.0+bf*2/PI/r2*sin(0.5*PI*(r-r2)/bf))*(rpol-r)/(rpol-r2)
        else:
            return 0.0
    if False:
        if r < rpol:
            return 0.5*bz*r2*(1.0+bf*2/PI/r2*sin(0.5*PI*(r-r2)/bf))*0.5*(1.0+cos(PI*(max(r,rout)-rout)/(rpol-rout)))
        else:
            return 0.0
    if True:
        bz1 = 0.3*1.86E-6
        rpol2 = r2*sqrt((bz+bz1)/bz1)
        if r < rpol2:
            return -0.5*bz1*r+0.5*(bz1+bz)*r2**2/r
        else:
            return 0
    
    
def Bz(rr, Aphis):
    Bz = np.zeros(len(rr))
    for i in range(len(rr)-1):
        Bz[i] = 1/rr[i]*(rr[i+1]*Aphis[i+1]-rr[i]*Aphis[i])/(rr[i+1]-rr[i])
        #Bz[i] = (Aphis[i+1]-Aphis[i])/(rr[i+1]-rr[i])
    return Bz


    

In [ ]:

    

rr = np.arange(0.0, rpol+bf, (rpol+bf)/1000)
Aphis = np.array([Aphi(r) for r in rr])/1E-10

plt.plot(rr,Aphis, label='Aphi')
plt.plot(rr,Bz(rr,Aphis)*r1, label='Bz')
plt.axvline(r2, color='black')
plt.axvline(rout, color='black')
plt.axvline(rpol, color='black')
plt.axhline(0.0, color='black')
plt.legend(loc=0)
#plt.xlim(28,30)
#plt.ylim(-1E-6, 1E-6)
plt.show()


    

In [ ]:

    

Rs = np.genfromtxt('/home/ychen/d9/FLASH4/MHD_Jet_3D/0829_test/efield.txt', usecols=0)
Ex = np.genfromtxt('/home/ychen/d9/FLASH4/MHD_Jet_3D/0829_test/efield.txt', usecols=1)
plt.scatter(Rs, Ex)
Bzs = Bz(Rs,Ex)
#plt.plot(rr,Aphis, label='Aphi')
plt.ylim(-3E-6, 3E-6)


    

In [ ]:

    

t = []
Bphi = []
Bzs = []

for file in files[0:10]:
    #print file
    ds = yt.load(file.fullpath)
    t.append(ds.current_time)
    Bx = ds.find_field_values_at_point('magnetic_field_x', posb)
    print Bx
    By = ds.find_field_values_at_point('magnetic_field_y', posb)
    Bz = ds.find_field_values_at_point('magnetic_field_z', posb)
    Bzs.append(Bz)
    Bphi.append(sqrt(Bx**2+By**2))
    #Bphi.append(-Bx)

plt.plot(t,Bzs, 'o-', label='Bz')
plt.plot(t,Bphi, 'o-', label='Bphi')
plt.axhline(y=5.2548E-6*np.sqrt(4*np.pi))
plt.legend(loc=0)
plt.xlabel('t')
plt.ylabel('B')
print Bzs[-1]
print Bphi[-1]