notebook.community

Edit and run



In [ ]:

    
import matplotlib.pyplot as plt
import numpy as np
e = 4.803E-10 #esu
me = 9.109E-28 #g
c = 2.998E10 #cm/s



In [ ]:

    
xmin, xmax = -10, 10
ymin, ymax = -10, 10
nx, ny = 41, 41
x = np.linspace(xmin,xmax,nx)
y = np.linspace(ymin,ymax,ny)
xx, yy = np.meshgrid(x,y)

plt.imshow(xx, extent=(xmin,xmax,ymin,ymax), origin='lower')
plt.colorbar()



In [ ]:

    
B0 = 1E-5
radius = 5
B = B0*np.exp(-(xx**2+yy**2)/radius**2)

fig = plt.figure(figsize=(16,6))
ax1 = fig.add_subplot(121)
im1 = ax1.imshow(B, extent=(xmin,xmax,ymin,ymax))
cb1 = plt.colorbar(im1)
ax1.set_title('B')

P0 = 5E-10
P = P0*np.power(10,np.cos(xx/5*np.pi))

ax2 = fig.add_subplot(122)

im2 = ax2.imshow(np.log10(P), extent=(xmin,xmax,ymin,ymax))
cb2 = plt.colorbar(im2)
cb2.set_label('log P')
ax2.set_title('P')



In [ ]:

    
plt.figure(figsize=(8,6))
plt.imshow(np.log10(P*B**1.5), extent=(xmin,xmax,ymin,ymax))
plt.colorbar()
plt.title('P*B^1.5')



In [ ]:

    
def nn_deposit(positions, deposit_field):
    grid = np.zeros((ny,nx))
    distfield = np.inf*np.ones((ny,nx))
    def process(pos, field_val):
        for ix in range(nx):
            for iy in range(ny):
                dist = (x[ix]-pos[0])**2 + (y[iy]-pos[1])**2
                if dist < distfield[iy,ix]:
                    distfield[iy,ix] = dist
                    grid[iy,ix] = field_val
    
    for ip, pos in enumerate(positions):
        process(pos, deposit_field[ip])
    return grid, distfield



In [ ]:

    
def nnw_deposit(positions, deposit_field, power):
    nwpart = 10
    nnfield = np.zeros((ny,nx,nwpart))
    distfield = np.zeros((ny,nx,nwpart))
    distarg = np.zeros((ny,nx,nwpart), dtype=int)
    def process(pos, field_val):
        for ix in range(nx):
            for iy in range(ny):
                wdist = 1.0/((x[ix]-pos[0])**power + (y[iy]-pos[1])**power)
                if wdist < distfield[iy,ix,distarg[iy,ix,0]]:
                    continue
                for p in distarg[iy,ix,:]:
                    if wdist > distfield[iy,ix,p]:
                        distfield[iy,ix,p] = wdist
                        nnfield[iy,ix,p] = field_val
                        distarg[iy,ix,:] = np.argsort(distfield[iy,ix,:])
                        break
    
    for ip, pos in enumerate(positions):
        process(pos, deposit_field[ip])
    nnw = np.sum(nnfield[:,:,:]*distfield[:,:,:], axis=-1)
    distw = np.sum(distfield[:,:,:], axis=-1)

    return nnw/distw



In [ ]:

    
vmin, vmax = -50, -34
vmin2, vmax2 = -38, -34

def plot_emission(field, particle_dtau):
    fig = plt.figure(figsize=(16,6))
    ax1 = fig.add_subplot(121)
    im = ax1.imshow(np.log10(field), extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=vmin, vmax=vmax)
    cb1 = plt.colorbar(im)
    sc1 = ax1.scatter(posx, posy, c=np.log10(particle_dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
    ax1.set_xlim(xmin, xmax)
    ax1.set_ylim(ymin, ymax)

    ax2 = fig.add_subplot(122)
    im2 = ax2.imshow(np.log10(field), extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=vmin2, vmax=vmax2)
    cb2 = plt.colorbar(im2)
    sc2 = ax2.scatter(posx, posy, c=np.log10(particle_dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
    ax2.set_xlim(xmin, xmax)
    ax2.set_ylim(ymin, ymax)

    return fig

Generate Particles dtau values



In [ ]:

    
npart = 50
vminp, vmaxp = -6, -2
posx = xmin + (xmax-xmin)*np.random.rand(npart)
posy = ymin + (ymax-ymin)*np.random.rand(npart)
#pfield = np.random.lognormal(mean=-30, sigma=50, size=npart)
dtau = np.random.lognormal(mean=-8.6, sigma=2, size=npart)
n, bins, patches = plt.hist(np.log10(dtau), range=(vminp,vmaxp), bins=20)

bin_centers = 0.5 * (bins[:-1] + bins[1:])
# scale values to interval [0,1]
cols = bin_centers - min(bin_centers)
cols /= max(cols)

cm = plt.cm.get_cmap('jet_r')
for col, p in zip(cols, patches):
    plt.setp(p, 'facecolor', cm(col))

plt.show()



In [ ]:

    
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
im = ax.imshow(np.log10(P*B**1.5), extent=(xmin,xmax,ymin,ymax), origin='lower')
cb1 = plt.colorbar(im)
cb1.set_label(u'grid $\log_{10}{PB^{3/2}}$')
sc = ax.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
cb2 = plt.colorbar(sc)
cb2.set_label(u'particle $\log_{10}{dtau}}$')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

Interpolate B and P values to particle positions



In [ ]:

    
from scipy.interpolate import interp2d
B_inter_f = interp2d(x, y, B, kind='cubic')
B_part = np.array([B_inter_f(px, py)[0] for (px, py) in zip(posx, posy)])
P_inter_f = interp2d(x, y, P, kind='cubic')
P_part = np.array([P_inter_f(px, py)[0] for (px, py) in zip(posx, posy)])
#B_part = B_inter_f(posx, posy).diagonal() #not working



In [ ]:

    
fig = plt.figure(figsize=(16,6))
ax1 = fig.add_subplot(121)
im1 = ax1.imshow(B, extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=0, vmax=1E-5)
sc1 = ax1.scatter(posx, posy, c=B_part, edgecolors='w', vmin=0, vmax=1E-5)
cb1 = plt.colorbar(im1)
ax1.set_xlim(xmin, xmax)
ax1.set_ylim(ymin, ymax)
ax1.set_title('B')

ax2 = fig.add_subplot(122)
im2 = ax2.imshow(P, extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=0, vmax=5E-9)
sc2 = ax2.scatter(posx, posy, c=P_part, edgecolors='w', vmin=0, vmax=5E-9)
cb2 = plt.colorbar(im2)
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(ymin, ymax)
ax2.set_title('P')

Distance weighted - spectra



In [ ]:

    
nu = 1.5E8
gamma_min = 10

gamc = 1/dtau
nuc = 3*gamc**2*e*B_part/(4*np.pi*me*c)
norm = 3/8*e**3.5/(c**2.5*me**1.5*np.pi**0.5)
N0 = 3/me/c/c/(np.log(gamc/gamma_min))/4/np.pi
sync_spec_part = N0*norm*nu**(-0.5)*np.exp(-nu/nuc)
deposited_nnw = P*B**1.5*nnw_deposit(zip(posx, posy), sync_spec_part, 2)
fig = plot_emission(deposited_nnw, dtau)

Nearest Neighbor - dtau



In [ ]:

    
deposited, distfield = nn_deposit(zip(posx, posy), dtau)

fig = plt.figure(figsize=(16,6))
ax = fig.add_subplot(121)
im = ax.imshow(np.log10(deposited), extent=(xmin,xmax,ymin,ymax), cmap='jet_r', origin='lower', vmin=vminp, vmax=vmaxp)
cb1 = plt.colorbar(im)
sc = ax.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_title('dtau')

ax2 = fig.add_subplot(122)
im2 = ax2.imshow(np.log10(1/deposited), extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=1, vmax=6)
cb2 = plt.colorbar(im2)
sc2 = ax2.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(ymin, ymax)
ax2.set_title('gamc')

Distance weighted - dtau



In [ ]:

    
deposited_nnw = nnw_deposit(zip(posx, posy), dtau, 2)

fig = plt.figure(figsize=(16,6))
ax = fig.add_subplot(121)
im = ax.imshow(np.log10(deposited_nnw), cmap='jet_r', extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=vminp, vmax=vmaxp)
cb1 = plt.colorbar(im)
cb1.set_label('log dtau')
sc = ax.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
#cb2 = plt.colorbar(sc)
#cb2.set_label(u'$\log_{10}{gmax}}$')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_title('dtau')

gamc = 1/deposited_nnw

ax2 = fig.add_subplot(122)
im2 = ax2.imshow(np.log10(gamc), cmap='jet', extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=1, vmax=6)
cb2 = plt.colorbar(im2)
cb2.set_label('log gamc')
sc2 = ax2.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
#cb2 = plt.colorbar(sc)
#cb2.set_label(u'$\log_{10}{gmax}}$')
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(ymin, ymax)
ax2.set_title('gamc')

Emission calculated with distance weighted dtau



In [ ]:

    
nu = 1.5E8

gamc = 1/deposited_nnw
nuc = 3*gamc**2*e*B/(4*np.pi*me*c)
fig = plt.figure(figsize=(16,5))
ax = fig.add_subplot(121)
im = ax.imshow(np.log10(nuc), cmap='jet', extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=6, vmax=10)
cb1 = plt.colorbar(im)
cb1.set_label(r'grid $\log \nu_c$')
sc1 = ax.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
cb2 = plt.colorbar(sc)
cb2.set_label(u'particle $\log_{10}{dtau}}$')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

norm = 3/8*B**1.5*e**3.5/(c**2.5*me**1.5*np.pi**0.5)
N0 = P*3/me/c/c/(np.log(gamc/gamma_min))/4/np.pi
sync_spec = N0*norm*nu**(-0.5)*np.exp(-nu/nuc)
ax2 = fig.add_subplot(122)
im2 = ax2.imshow(np.log10(sync_spec), extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=vmin2, vmax=vmax2)
cb2 = plt.colorbar(im2)
cb2.set_label('emission at %.1e Hz' % nu)
sc2 = ax2.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(ymin, ymax)



In [ ]:

    
fig = plot_emission(sync_spec, dtau)



In [ ]:

    
plt.hist(np.log10(sync_spec.flatten()+1E-100), range=(-80, -30))

Calculate Emission from Particles



In [ ]:

    
gamc = 1/dtau
nuc = 3*gamc**2*e*B_part/(4*np.pi*me*c)
norm = 3/8*B_part**1.5*e**3.5/(c**2.5*me**1.5*np.pi**0.5)
N0 = P_part*3/me/c/c/(np.log(gamc/gamma_min))/4/np.pi
sync_spec_part = N0*norm*nu**(-0.5)*np.exp(-nu/nuc)

Nearest neighbor - emission



In [ ]:

    
deposited, distfield = nn_deposit(zip(posx, posy), sync_spec_part)
fig = plot_emission(deposited, dtau)

Distance weighted - emission



In [ ]:

    
deposited_nnw = nnw_deposit(zip(posx, posy), sync_spec_part, 2)
fig = plot_emission(deposited_nnw, dtau)

Distance weighted - log emission



In [ ]:

    
deposited_nnw = nnw_deposit(zip(posx, posy), np.log10(sync_spec_part+np.finfo(np.float64).tiny), 2)
deposited_nnw = np.power(10, deposited_nnw)
fig = plot_emission(deposited_nnw, dtau)



In [ ]:

Distance weighted - log dtau



In [ ]:

    
deposited_nnw = nnw_deposit(zip(posx, posy), np.log10(dtau), 2)
deposited_nnw = np.power(10, deposited_nnw)

fig = plt.figure(figsize=(16,6))
ax = fig.add_subplot(121)
im = ax.imshow(np.log10(deposited_nnw), cmap='jet_r', extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=vminp, vmax=vmaxp)
cb1 = plt.colorbar(im)
cb1.set_label('log dtau')
sc = ax.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
#cb2 = plt.colorbar(sc)
#cb2.set_label(u'$\log_{10}{gmax}}$')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_title('dtau')

gamc = 1/deposited_nnw

ax2 = fig.add_subplot(122)
im2 = ax2.imshow(np.log10(gamc), cmap='jet', extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=1, vmax=6)
cb2 = plt.colorbar(im2)
cb2.set_label('log gamc')
sc2 = ax2.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
#cb2 = plt.colorbar(sc)
#cb2.set_label(u'$\log_{10}{gmax}}$')
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(ymin, ymax)
ax2.set_title('gamc')

Emission calculated with distance weighted log dtau



In [ ]:

    
gamc = 1/deposited_nnw
nuc = 3*gamc**2*e*B/(4*np.pi*me*c)
fig = plt.figure(figsize=(16,5))
ax = fig.add_subplot(121)
im = ax.imshow(np.log10(nuc), extent=(xmin,xmax,ymin,ymax), cmap='jet', origin='lower', vmin=6, vmax=10)
cb1 = plt.colorbar(im)
cb1.set_label(r'$\log \nu_c$')
sc1 = ax.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
cb2 = plt.colorbar(sc)
cb2.set_label(u'$\log_{10}{dtau}}$')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

norm = 3/8*B**1.5*e**3.5/(c**2.5*me**1.5*np.pi**0.5)
N0 = P*3/me/c/c/(np.log(gamc/gamma_min))/4/np.pi
sync_spec = N0*norm*nu**(-0.5)*np.exp(-nu/nuc)
ax2 = fig.add_subplot(122)
im2 = ax2.imshow(np.log10(sync_spec), extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=vmin2, vmax=vmax2)
cb2 = plt.colorbar(im2)
cb2.set_label('emission at %.1e Hz' % nu)
sc2 = ax2.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(ymin, ymax)



In [ ]:

    
fig = plot_emission(sync_spec, dtau)



In [ ]:

Distance weighted - log gamc



In [ ]:

    
gamc = 1/dtau
deposited_nnw = nnw_deposit(zip(posx, posy), np.log10(gamc), 2)
deposited_nnw = np.power(10, deposited_nnw)

fig = plt.figure(figsize=(16,5))
ax1 = fig.add_subplot(121)
im1 = ax1.imshow(np.log10(deposited_nnw), cmap='jet', extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=2, vmax=5)
cb1 = plt.colorbar(im1)
sc1 = ax1.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
ax1.set_xlim(xmin, xmax)
ax1.set_ylim(ymin, ymax)
ax1.set_title('gamc')

nuc = 3*deposited_nnw**2*e*B/(4*np.pi*me*c)

ax2 = fig.add_subplot(122)
im2 = ax2.imshow(np.log10(nuc), cmap='jet', extent=(xmin,xmax,ymin,ymax), origin='lower', vmin=6, vmax=10)
cb2 = plt.colorbar(im)
cb2.set_label(r'grid $\log \nu_c$')
sc2 = ax2.scatter(posx, posy, c=np.log10(dtau), cmap='jet_r', edgecolors='w', vmin=vminp, vmax=vmaxp)
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(ymin, ymax)
ax2.set_title('nuc')

Emission calculated with distance weighted log gamc



In [ ]:

    
sync_spec = N0*norm*nu**(-0.5)*np.exp(-nu/nuc)
fig = plot_emission(sync_spec, dtau)



In [ ]: