In [ ]:

    

%matplotlib inline
import numpy as np
from numpy import cos, sin, sqrt
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 150


    

In [ ]:

    

r2 = 2E10
b = 1E10
rout = r2 + b
r3 = r2 + 2.*b
r4 = r3 + r2
rpol = r4 + b
bz0 = 2.0
pi = np.pi

rr = np.linspace(0, rpol+b , 1000)
delta = rr[1]-rr[0]

def bz(r):
    if r < r2:
        return bz0
    elif r < rout:
        return 0.5*bz0*(1.0+cos(pi*((r-r2)/b)))
    elif r < r3:
        return -0.5*bz0*(1.0+cos(pi*((r-r3)/b)))
    elif r < r4:
        return -bz0
    elif r < rpol:
        return -0.5*bz0*(1.0+cos(pi*((r-r4)/b)))
    else:
        return 0.0

def bzs(rr):
    return np.array([bz(r) for r in rr])



plt.plot(rr, bzs(rr))
plt.axvline(r2, ls=':')
plt.axvline(rout, ls=':')
plt.axvline(r3, ls=':')
plt.axvline(r4, ls=':')
plt.axvline(rpol, ls=':')
plt.show()


    

In [ ]:

    

print(delta)

def Aphi(r):
    #c1 = (0.5*r2+2.0*c0/r2)*pi**2*r2/b - b
    c1 = 0.5*(r2*pi)**2/b - b
    c2 = -(rout*pi)**2/b + 2*b - c1
    c3 = 0.25*r3*r3-0.5*b/pi**2*(b+c2)
    c4 = (-0.25*r4 + c3/r4)*(-2.0*pi**2*r4/b) - b
    #c0 = (0.25*rpol+0.5*b/pi**2/rpol*(-b+c4))*rpol
    c0 = 0
    
    if r<delta:
        return c0/delta+0.5*r
    elif r>=delta and r<r2:
        return 0.5*r + c0/r
    elif r>=r2 and r<rout:
        return (0.25*r + 0.5*b/pi**2/r*( b*cos(pi*((r-r2)/b))+pi*r*sin(pi*((r-r2)/b)) + c1 ) + c0/r)
    elif r>=rout and r<r3:
        return (-0.25*r - 0.5*b/pi**2/r*( b*cos(pi*((r-r3)/b))+pi*r*sin(pi*((r-r3)/b)) + c2 ) + c0/r)
    elif r>=r3 and r<r4:
        return (-0.5*r + c3/r + c0/r)
    elif r>=r4 and r<rpol:
        return (-0.25*r - 0.5*b/pi**2/r*( b*cos(pi*((r-r4)/b))+pi*r*sin(pi*((r-r4)/b)) + c4 ) + c0/r)
    else:
        return 0.0
    
    
def Aphis(rr):
    return np.array([bz0*Aphi(r) for r in rr])

def Bz_deriv(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

plt.plot(rr, Aphis(rr), label='Aphi')
plt.plot(rr, Bz_deriv(rr, Aphis(rr))*r2, label='Bz')
plt.axvline(r2, ls=':')
plt.axvline(rout, ls=':')
plt.axvline(r3, ls=':')
plt.axvline(r4, ls=':')
plt.axvline(rpol, ls=':')
plt.axhline(0, ls=':')
plt.ylim(-5E10,5E10)
#plt.xlim(0, 10)
#plt.semilogy()
plt.show()


    

In [ ]:

    

pi = np.pi

r2 = 1
bf = 0.25
b = bf
rout = r2 + bf
r3 = r2 + 2*bf


bz0=1
bz1=0.1*bz0


def Bz(r):
    if r < r2:
        return bz0
    elif r < rout:
        return 0.5*bz0*(1.0+cos(pi*((r-r2)/b)))
    elif r < r3:
        return -0.5*bz1*(1.0+cos(pi*((r-r3)/b)))
    elif r < r4:
        return -bz1
    elif r < rpol:
        return -0.5*bz1*(1.0+cos(pi*((r-r4)/b)))
    else:
        return 0.0
    
def Bzs(rr):
    return np.array([Bz(r) for r in rr])
    
def Bz_deriv(rr, Aphis):
    Bz = np.zeros(len(rr))
    for i in range(1,len(rr)-1):
        Bz[i] = 1/rr[i]*(rr[i+1]*Aphis[i+1]-rr[i-1]*Aphis[i-1])/(rr[i+1]-rr[i-1])
        #Bz[i] = (Aphis[i+1]-Aphis[i])/(rr[i+1]-rr[i])
    return Bz

#c1 = (0.5*r2+2.0*c0/r2)*pi**2*r2/b - b
c1 = 0.5*(r2*pi)**2/bf - bf
c2 = (-0.5*(bz0+bz1)*(rout*pi)**2/bf - bz0*(-bf+c1))/bz1 + bf
c3 = -0.25*r3*r3+0.5*bf/pi**2*(bf+c2)
r4 = -0.5*bf + sqrt(-2.0*c3+2.0*(bf/pi)**2-0.25*bf**2)
rpol = r4+bf
c4 = (0.25*r4**2 + c3)*(2.0*pi**2/bf) - bf
#c0 = (0.25*rpol+0.5*b/pi**2/rpol*(-b+c4))*rpol
c0 = 0

def Aphi2(r):
    if r<r2:
        return 0.5*bz0*r + c0/r
    elif r>=r2 and r<rout:
        return bz0*(0.25*r + 0.5*b/pi**2/r*( b*cos(pi*((r-r2)/b))+pi*r*sin(pi*((r-r2)/b)) + c1 ) + c0/r)
    elif r>=rout and r<r3:
        return -bz1*(0.25*r + 0.5*b/pi**2/r*( b*cos(pi*((r-r3)/b))+pi*r*sin(pi*((r-r3)/b)) + c2 ) + c0/r)
    elif r>=r3 and r<r4:
        return -bz1*(0.5*r + c3/r + c0/r)
    elif r>=r4 and r<rpol:
        return -bz1*(0.25*r + 0.5*b/pi**2/r*( b*cos(pi*((r-r4)/b))+pi*r*sin(pi*((r-r4)/b)) + c4 ) + c0/r)
    else:
        return 0.0
    
    
def Aphi2s(rr):
    return np.array([Aphi2(r) for r in rr])

rr = np.linspace(0, rpol+b , 320)
Amax = np.nan_to_num(Aphi2s(rr)).max()
plt.plot(rr, Aphi2s(rr)/Amax, label=u'$A_{\phi}$', c='navy')
plt.plot(rr, Bzs(rr), label=u'$B_z$', c='r')
#plt.plot(rr, Bz_deriv(rr, Aphi2s(rr))*r2, label='Bz(deriv)')
plt.axvline(r2, ls=':', c='gray', lw=1)
plt.axvline(rout, ls=':', c='gray', lw=1)
plt.axvline(r3, ls=':', c='gray', lw=1)
plt.axvline(r4, ls=':', c='gray', lw=1)
plt.axvline(rpol, ls=':', c='gray', lw=1)
plt.axhline(0, ls=':', c='gray', lw=1)
#plt.ylim(-5E10,5E10)
plt.ylim(-0.2, 1.1)
#plt.semilogy()
plt.xlabel(u'$r_\mathrm{nozzle}$')
plt.ylabel('arbitrary unit')
plt.legend(loc=1)

plt.show()


    

In [ ]:

    

pi = np.pi

r2 = 1
bf = 0.25
b = bf
rout = r2 + bf
r3 = r2 + 2*bf


bz0=1
bz1=0.1*bz0


def Bz(r):
    if r < r2:
        return bz0
    elif r < rout:
        return 0.5*bz0*(1.0+cos(pi*((r-r2)/b)))
    elif r < r3:
        return -0.5*bz1*(1.0+cos(pi*((r-r3)/b)))
    elif r < r4:
        return -bz1
    elif r < rpol:
        return -0.5*bz1*(1.0+cos(pi*((r-r4)/b)))
    else:
        return 0.0
    
def Bzs(rr):
    return np.array([Bz(r) for r in rr])
    
def Bz_deriv(rr, Aphis):
    Bz = np.zeros(len(rr))
    for i in range(1,len(rr)-1):
        Bz[i] = 1/rr[i]*(rr[i+1]*Aphis[i+1]-rr[i-1]*Aphis[i-1])/(rr[i+1]-rr[i-1])
        #Bz[i] = (Aphis[i+1]-Aphis[i])/(rr[i+1]-rr[i])
    return Bz

#c1 = (0.5*r2+2.0*c0/r2)*pi**2*r2/b - b
c1 = 0.5*(r2*pi)**2/bf - bf
c2 = (-0.5*(bz0+bz1)*(rout*pi)**2/bf - bz0*(-bf+c1))/bz1 + bf
c3 = -0.25*r3*r3+0.5*bf/pi**2*(bf+c2)
r4 = -0.5*bf + sqrt(-2.0*c3+2.0*(bf/pi)**2-0.25*bf**2)
rpol = r4+bf
c4 = (0.25*r4**2 + c3)*(2.0*pi**2/bf) - bf
#c0 = (0.25*rpol+0.5*b/pi**2/rpol*(-b+c4))*rpol
c0 = 0

def Aphi2(r):
    if r<r2:
        return 0.5*bz0*r + c0/r
    elif r>=r2 and r<rout:
        return bz0*(0.25*r + 0.5*b/pi**2/r*( b*cos(pi*((r-r2)/b))+pi*r*sin(pi*((r-r2)/b)) + c1 ) + c0/r)
    elif r>=rout and r<r3:
        return -bz1*(0.25*r + 0.5*b/pi**2/r*( b*cos(pi*((r-r3)/b))+pi*r*sin(pi*((r-r3)/b)) + c2 ) + c0/r)
    elif r>=r3 and r<r4:
        return -bz1*(0.5*r + c3/r + c0/r)
    elif r>=r4 and r<rpol:
        return -bz1*(0.25*r + 0.5*b/pi**2/r*( b*cos(pi*((r-r4)/b))+pi*r*sin(pi*((r-r4)/b)) + c4 ) + c0/r)
    else:
        return 0.0
    
    
def Aphi2s(rr):
    return np.array([Aphi2(r) for r in rr])



r1 = 0.5*r2

def Bphi(r):
    if r < r1:
        return -2*r**3/r1**3 + 3*r**2/r1**2
    elif r < r2:
        return 1.0
    elif r < rout:
        return -2*(rout-r)**3/(rout-r2)**3 + 3*(rout-r)**2/(rout-r2)**2
    else:
        return 0.0

def Bphis(rr):
    return np.array([Bphi(r) for r in rr])


def Az(r):
    if r < r1:
        return r**4/(2*r1**3) - r**3/r1**2 + 0.5*(-r1+rout+r2)
    elif r < r2:
        return -r + 0.5*(rout+r2)
    elif r < rout:
        return -(rout-r)**4/(2*(rout-r2)**3) + (rout-r)**3/(rout-r2)**2
    else:
        return 0.0

def Azs(rr):
    return np.array([Az(r) for r in rr])


def v(r):
    if r < r2:
        return 1.0
    elif r < rout:
        return 0.5*(1.0+np.cos(np.pi*(r-r2)/bf))
    else:
        return 0.0
    
def vs(r):
    return np.array([v(r) for r in rr])


    

In [ ]:

    

rr = np.linspace(0, rpol+b , 320)

fig = plt.figure(figsize=(5,5))

ax1 = fig.add_subplot(311)

plt.plot(rr, Bzs(rr), c='r', label=u'$B_z$')
Amax = np.nan_to_num(Aphi2s(rr)).max()
plt.plot(rr, Aphi2s(rr)/Amax, c='r', ls=':', label=u'$A_{\phi}$')
#plt.plot(rr, rr*Bzs(rr)/Bphis(rr), ls='--', label=u'$p$')

plt.axvline(r1, ls=':', c='gray', lw=1)
plt.axvline(r2, ls=':', c='gray', lw=1)
plt.axvline(rout, ls=':', c='gray', lw=1)
plt.axvline(r3, ls=':', c='gray', lw=1)
plt.axvline(r4, ls=':', c='gray', lw=1)
plt.axvline(rpol, ls=':', c='gray', lw=1)
plt.axhline(0, ls=':', c='gray', lw=1)
plt.ylim(-0.2,1.2)
#plt.xlim(0, 10)
#plt.semilogy()
plt.setp(ax1.get_xticklabels(), visible=False)
plt.ylim(-0.2,1.1)
ax1.set_yticks([-0.2,0.0,0.2,0.4,0.6,0.8,1.0])
ax1.text(r1,1.14, r'$r_1$', horizontalalignment='center')
ax1.text(r2,1.14, r'$r_2$', horizontalalignment='center')
ax1.text(rout,1.14, r'$r_{\mathrm{out}}$', horizontalalignment='center')
ax1.text(r3,1.14, r'$r_3$', horizontalalignment='center')
ax1.text(r4,1.14, r'$r_4$', horizontalalignment='center')
ax1.text(rpol,1.14, r'$r_{\mathrm{pol}}$', horizontalalignment='center')
plt.xlim(0,4.2)
plt.ylabel('arbitrary unit')
plt.yticks([0,0.5,1])
plt.legend(loc=1)

ax2 = fig.add_subplot(312, sharex=ax1)
plt.plot(rr, Bphis(rr), c='b', label=u'$B_{\phi}$')
Amax = np.nan_to_num(Azs(rr)).max()
plt.plot(rr, Azs(rr)/Amax, c='b', ls=':', label=u'$A_{z}$')

plt.axvline(r1, ls=':', c='gray', lw=1)
plt.axvline(r2, ls=':', c='gray', lw=1)
plt.axvline(rout, ls=':', c='gray', lw=1)
plt.axvline(r3, ls=':', c='gray', lw=1)
plt.axvline(r4, ls=':', c='gray', lw=1)
plt.axvline(rpol, ls=':', c='gray', lw=1)
plt.axhline(0, ls=':', c='gray', lw=1)
plt.ylim(-0.2,1.2)
plt.setp(ax2.get_xticklabels(), visible=False)
plt.ylim(-0.2,1.1)

plt.xlim(0,4.2)
plt.ylabel('arbitrary unit')
plt.legend(loc=1)



ax3 = fig.add_subplot(313, sharex=ax1)
plt.plot(rr, vs(rr), label=u'$v$', c='#2ca02c')
plt.axvline(r2, ls=':', c='gray', lw=1)
plt.axvline(rout, ls=':', c='gray', lw=1)
plt.axvline(r3, ls=':', c='gray', lw=1)
plt.axvline(r4, ls=':', c='gray', lw=1)
plt.axvline(rpol, ls=':', c='gray', lw=1)
plt.axhline(0, ls=':', c='gray', lw=1)

plt.legend(loc=1)
plt.xlabel(u'$r/r_\mathrm{jet}$')
plt.xlim(0,4.2)
plt.ylabel(u'$v/v_{jet}$')
plt.yticks([0,0.5,1])
plt.tight_layout(h_pad=0.02)
plt.savefig('B_v_profile_r.pdf')


    

In [ ]: