In [ ]:
#%pdb
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.special import kv
from scipy.special import gamma as gamma_func
from scipy.integrate import quad
from scipy.optimize import curve_fit
In [ ]:
me = 9.109E-28 # electron mass in g
c = 2.998E10 # speed of light in cm/s
e = 4.803E-10 # electron charge in esu
B53 = lambda x: kv(5./3., x)
def F(x):
x = min(x, 1E3)
integral, err = quad(B53, x, np.Inf)#, epsabs=0, epsrel=1.49e-08)
#print x, np.Inf, integral
ub = 1E9*x
while integral < 0.0 or np.isnan(integral):
integral, err = quad(B53, x, ub)
#print x, ub, integral
ub /= 10
if ub <= x: raise ArithmeticError
return x*integral
def G(x):
return x*kv(2./3., x)
In [ ]:
xx=np.logspace(-6, np.log10(4), 1000)
Fxx = np.array([F(x) for x in xx])
Gxx = np.array([G(x) for x in xx])
plt.plot(xx,Fxx, label='F(x)')
plt.plot(xx,Gxx, label='G(x)')
plt.plot(xx,Gxx/Fxx, ':', label='G(x)/F(x)')
plt.xlabel('x')
plt.legend()
plt.show()
#plt.ylabel('F(x)')
In [ ]:
xx=np.logspace(-6, 1, 1000)
Fxx = np.array([F(x) for x in xx])
Gxx = np.array([G(x) for x in xx])
plt.figure(facecolor='w')
# Perpendicular component F + G
plt.plot(xx,Fxx+Gxx, '--', label=r'$\perp$')
# Parallel component F - G
plt.plot(xx,Fxx-Gxx, ':', label=r'$\parallel$')
# Summation of the two
plt.plot(xx,2*Fxx, color='k', label='Total')
plt.xlabel(r'$x=\nu/\nu_c$')
plt.semilogx()
plt.legend(loc=2)
#plt.ylabel('F(x)')
In [ ]:
# exponential + power law fitting function
def F_fit(x, A, pind):
return A*x**pind*np.exp(-x)
xx=np.logspace(-6,2, 1000)
Fxx = np.array([F(x) for x in xx])
par, cov = curve_fit(F_fit, xx, Fxx, p0=[2.0, 0.33])
A, pind = par
print(par)
print(cov)
# fitting function in log space
def F_fit_log10(x, A, pind):
return np.log10(A*x**pind*np.exp(-x))
par, cov = curve_fit(F_fit_log10, xx, np.log10(Fxx), p0=[2.0, 0.33])
A2, pind2 = par
print(par)
print(cov)
In [ ]:
xx=np.logspace(-6, 2, 100)
Fxx = np.array([F(x) for x in xx])
#plt.plot(xx,[F(x) for x in xx])
#plt.xlim(0,4.0)
#plt.show()
fig = plt.figure(figsize=(8,6))
ax1=fig.add_axes((.1,.5,.8,.6))
ax1.plot(xx,Fxx, c='r', label='F(x)')
ax1.plot(xx,F_fit(xx, A, pind), c='b', ls=':', alpha=0.5, label=r'$%.3f\,x^{%.3f}e^{-x}$ (linear fit)' % (A,pind))
ax1.plot(xx,F_fit(xx, A2, pind2), c='g', ls='--', alpha=0.5, label=r'$%.3f\,x^{%.3f}e^{-x}$ (log fit)' % (A2,pind2))
ax1.semilogy()
ax1.semilogx()
ax1.set_ylim(1E-2,1.0)
#plt.xlim(0.01, 10)
ax1.legend(loc=2)
ax1.xaxis.set_ticks_position('top')
#ax1.set_xticklabels([]) #Remove x-tic labels for the first frame
ax1.grid()
ax2=fig.add_axes((.1,.3,.8,.2))
ax2.plot(xx,F_fit(xx, A, pind)-Fxx, c='b', ls=':')
ax2.plot(xx,F_fit(xx, A2, pind2)-Fxx, c='g', ls='--')
#plt.semilogy()
plt.semilogx()
ax2.set_ylabel('residue')
#resmax = np.abs(Fxx-F_fit(xx, A2, pind2)).max()
ax2.set_ylim(-0.1, 0.1)
ax2.set_yticks([-0.1, -0.05, 0.00, 0.05, 0.1])
ax2.set_yticklabels([-0.1, -0.05, 0.00, 0.05])
ax2.hlines(0, xx.min(), xx.max(), lw=0.2, alpha=0.5)
ax2.set_xticklabels([])
ax3=fig.add_axes((.1,.1,.8,.2))
ax3.plot(xx,(F_fit(xx, A, pind)-Fxx)/Fxx*100, c='b', ls=':')
ax3.plot(xx,(F_fit(xx, A2, pind2)-Fxx)/Fxx*100, c='g', ls='--')
#plt.semilogy()
ax3.semilogx()
ax3.set_ylabel('% error')
ax3.set_ylim(-40,40)
ax3.set_yticks([-40, -20, 0, 20, 40])
ax3.set_yticklabels([-40, -20, 0, 20])
ax3.hlines(0, xx.min(), xx.max(), lw=0.2, alpha=0.5)
plt.show()
In [ ]:
# Integrate over gamma
def func(gamma, B, nu):
#return 3**0.5*B*e**3/c**2/me*F(4*np.pi*me*c/3.0/gamma**2/e/B*nu)*gamma**(-2)
return F(4*np.pi*me*c/3.0/gamma**2/e/B*nu)*gamma**(-2)
def J(gmax, gmin, B, nu):
#print '(gmin,gmax) = ', gmin, gmax
#print '(xmin,xmax) = ', 4*np.pi*me*c/3.0/gmax**2/e/B*nu, 4*np.pi*me*c/3.0/gmin**2/e/B*nu
integral, err = quad(func, gmin, gmax, args=(B, nu))
#print integral, err
return 3**0.5*B*e**3/c**2/me*integral
In [ ]:
gmax=1E4
nu=1E9
B = 1E-5
#print 'nc = %e' % (1.0/(4.0*np.pi*me*c/3.0/gmax**2/e/B))
nc = 1.0/(4.0*np.pi*me*c/3.0/gmax**2/e/B)
print ('nc = %e' % nc)
print ('nu/nc = %e' % (nu/nc))
#print 'F(x(gmax)) =', F(4*np.pi*me*c/3.0/e/B/gmax**2*nu)
#print 'F(x(gmin)) =', F(4*np.pi*me*c/3.0/e/B/1.0**2*nu)
#print 3**0.5*B*e**3/c**2/me*gamma**(-2)
print (J(1E5, 1.0, 1E-5, 1E6))
print (J(1E4, 1.0, 1E-5, 1E6))
In [ ]:
nus=np.logspace(-5, 15, 50)
color = { 1E3: 'r', 1E4: 'g', 1E5: 'b'}
lw = {1E-5: 1, 1E-6: 0.5}
t0 = time.time()
for B in [1E-5, 1E-6]:
for gmax in [1E3, 1E4, 1E5]:
nc = 3.0*gmax**2*e*B/(4.0*np.pi*me*c)
print B, gmax, 'nc = %e' % nc
Js = np.array([J(gmax, 1.0, B, nu) for nu in nus])
plt.plot(nus, Js, '-', lw=lw[B], c=color[gmax])
plt.vlines(nc, 1E-34, 1E-26, lw=lw[B], linestyles=':', colors=color[gmax])
#plt.plot(nus, 5E-19*B**(1.5)*nus**(-0.5)*exp(-nus/nc), '--', lw=lw[B],
#c=color[gmax], alpha=0.7)
print ('Elapsed time:', time.time()-t0)
plt.semilogx()
plt.semilogy()
plt.ylim(1E-34, 1E-26)
plt.xlabel('Hz')
In [ ]:
# Integrate over x
def func(x):
#print 'x =', x, 'F(x) = ', F(x)
#xx.append(x)
#FF.append(x**(-0.5)*F(x))
return x**(-0.5)*F(x)
def J(gmax, gmin, B, nu):
xmin = 4.0*np.pi*me*c/3.0/gmax**2/e/B*nu
xmax = min(4.0*np.pi*me*c/3.0/gmin**2/e/B*nu, 1E3)
#print '(gmin,gmax) = ', gmin, gmax
#print '(xmin,xmax) = ', xmin, xmax
integral, err = quad(func, xmin, xmax)
#print integral, err
return 3.0/8.0*B**(1.5)*e**3.5/c**2.5/me**1.5/(np.pi)**0.5/nu**0.5*integral
#return 3**0.5*B*e**3/c**2/me**2*0.5*(3.0*B*e/(4*np.pi*me*c*nu))**0.5*integral
#return 0.5*(3.0*B*e/(4*np.pi*me*c*nu))**0.5*integral
#return 1.0/nu**0.5*integral
In [ ]:
def perp(x):
#print 'x =', x, 'F(x) = ', F(x)
#xx.append(x)
#FF.append(x**(-0.5)*F(x))
return x**(-0.5)*(F(x)+G(x))
def para(x):
return x**(-0.5)*(F(x)-G(x))
def J_perp(gmax, gmin, B, nu):
xmin = 4.0*np.pi*me*c/3.0/gmax**2/e/B*nu
xmax = min(4.0*np.pi*me*c/3.0/gmin**2/e/B*nu, 1E3)
#print '(gmin,gmax) = ', gmin, gmax
#print '(xmin,xmax) = ', xmin, xmax
integral, err = quad(perp, xmin, xmax)
#print integral, err
return 3.0/8.0*B**(1.5)*e**3.5/c**2.5/me**1.5/(np.pi)**0.5/nu**0.5*integral
#return 3**0.5*B*e**3/c**2/me**2*0.5*(3.0*B*e/(4*np.pi*me*c*nu))**0.5*integral
#return 0.5*(3.0*B*e/(4*np.pi*me*c*nu))**0.5*integral
#return 1.0/nu**0.5*integral
def J_para(gmax, gmin, B, nu):
xmin = 4.0*np.pi*me*c/3.0/gmax**2/e/B*nu
xmax = min(4.0*np.pi*me*c/3.0/gmin**2/e/B*nu, 1E3)
integral, err = quad(para, xmin, xmax)
return 3.0/8.0*B**(1.5)*e**3.5/c**2.5/me**1.5/(np.pi)**0.5/nu**0.5*integral
In [ ]:
xx = []
FF = []
print (J(1E5, 1.0, 1E-5, 1E2))
#plt.scatter(np.log10(xx), np.log10(FF), s=1, lw=0, c='r')
xx = []
FF = []
print (J(1E4, 1.0, 1E-5, 1E2))
#plt.scatter(np.log10(xx), np.log10(FF), s=1, lw=0, c='b')
#plt.ylim(-2,2)
In [ ]:
nus=np.logspace(-5, 15, 50)
color = {1E0: 'r', 1E1: 'g', 1E2: 'b', 1E3: 'r', 1E4: 'g', 1E5: 'b'}
lw = {1E-5: 2, 1E-6: 0.5}
gmax = 1E5
gmin = 1.0
t0 = time.time()
for B in [1E-5, 1E-6]:
#for gmin in [1E0, 1E1, 1E2]:
for gmax in [1E3, 1E4, 1E5]:
nuc = 3.0*gmax**2*e*B/(4.0*np.pi*me*c)
nucmin = 3.0*gmin**2*e*B/(4.0*np.pi*me*c)
print (B, gmax, gmin, 'nuc = %e' % nuc)
Js = np.array([J_perp(gmax, gmin, B, nu) for nu in nus])
plt.plot(nus, Js, '-', lw=lw[B], c=color[gmax])
plt.vlines(nuc, 1E-34, 1E-26, lw=lw[B], linestyles=':', colors=color[gmax])
#plt.vlines(nucmin, 1E-34, 1E-26, lw=lw[B], linestyles=':', colors=color[gmin])
#plt.plot(nus, 5E-19*B**(1.5)*nus**(-0.5)*exp(-nus/nc),
#'--', lw=lw[B], c=color[gmax], alpha=0.7)
print ('Elapsed time:', time.time()-t0)
#plt.grid()
plt.semilogx()
plt.semilogy()
plt.xlabel('Hz')
plt.ylim(1E-34, 1E-26)
In [ ]:
gmax=1E6
gmin=1E2
nu=1.5E9
xmin = 4.0*np.pi*me*c/3.0/gmax**2/e/B*nu
xmax = 4.0*np.pi*me*c/3.0/gmin**2/e/B*nu
print (xmin, F(xmin))
print (xmax, F(xmax))
In [ ]:
class Synchrotron():
B = 0.0
gmax = 0.0
gmin = 0.0
nus = 0.0
Js = 0.0
nucmin = 0.0
nuc = 0.0
def gen_sync(self, B, gmax, gmin=1.0, lognumin=-5, lognumax=13, n_nu=51, pol='para'):
self.B = B
self.gmax = gmax
self.gmin = gmin
self.nucmin = 3.0*gmin**2*e*B/(4.0*np.pi*me*c)
self.nuc = 3.0*gmax**2*e*B/(4.0*np.pi*me*c)
t0 = time.time()
nus = np.logspace(lognumin, lognumax, n_nu)
if pol == 'para':
Js = np.array([J_para(gmax, gmin, B, nu) for nu in nus])
elif pol == 'perp':
Js = np.array([J_perp(gmax, gmin, B, nu) for nu in nus])
self.nus = nus[Js.nonzero()]
self.Js = Js[Js.nonzero()]
print ('Synchrotron Spectra created with B = %.2e gauss, (gmax, gmin) = (%.1e, %.1e)' % (B, gmax, gmin))
print ('Elapsed time:', time.time()-t0)
__init__ = gen_sync
def f1(self, nus, A1, pind):
return A1*nus**pind
def f2(self, nus, A2):
return A2*nus**(-0.5)*np.exp(-nus/self.nuc)
def f(self, nus, A1, pind1, A2):
norm = 3.0/8.0*self.B**1.5*e**3.5/(c**2.5*me**1.5*(np.pi)**0.5)
return norm/(1./self.f1(nus, A1, pind1) + 1./self.f2(nus, A2))
def f_log(self, nus, A1, pind1, A2):
return np.log10(self.f(nus, A1, pind1, A2))
par3 = [1., 1./3., 1.]
def f_log_Aonly(self, nus, A1, A2):
return np.log10(self.f(nus, A1, self.par3[1], A2))
par1 = [1., 1./3., 1.]
par2 = [1., 1./3., 1.]
def fit_sync(self):
self.par1, cov1 = curve_fit(self.f, self.nus, self.Js, p0=self.par1)
#print par1
#print cov1
self.par2, cov2 = curve_fit(self.f_log, self.nus, np.log10(self.Js), p0=self.par2)
#print par2
#print cov2
(self.par3[0], self.par3[2]), cov3 = curve_fit(self.f_log_Aonly, self.nus, np.log10(self.Js), p0=[1, 1])
#print par3
#print cov3
def plot_sync(self, show=True):
nus = self.nus
Js = self.Js
par1 = self.par1
par2 = self.par2
par3 = self.par3
f = self.f
fig = plt.figure(figsize=(8,6))
ax1=fig.add_axes((.1,.4,.8,.5))
if pol == 'para':
label = u'$J_\parallel(x)$'
elif pol == 'perp':
label = u'$J_\perp(x)$'
ax1.plot(nus[::2], Js[::2], 'x', label=label, c='k', ms=4)
#ax1.plot(nus, f1(nus, par1[0], par1[1]))
#ax1.plot(nus, f2(nus, par1[2], par1[3]))
label = r'$(\frac{1}{%.3f\nu^{%.3f}}+\frac{1}{%.3f\nu^{-0.5}e^{-x}})^{-1}$'
ax1.plot(nus, f(nus, par1[0], par1[1], par1[2]), c='g', ls=':', alpha=0.5,
label=label % tuple(par1) +' (linear fit)')
ax1.plot(nus, f(nus, par2[0], par2[1], par2[2]), c='b', ls='--', alpha=0.5,
label=label % tuple(par2) +' (log fit)')
ax1.plot(nus, f(nus, par3[0], par3[1], par3[2]), c='r', ls='-.', alpha=0.5,
label=label % tuple(par3) +' (log fit)')
ax1.vlines(self.nucmin, 1E-34, 1E-26, lw=5, colors='grey', alpha=0.2)
ax1.vlines(self.nuc, 1E-34, 1E-26, lw=5, colors='grey', alpha=0.2)
ax1.semilogx()
ax1.semilogy()
ax1.set_xlim(1E-5, 1E13)
ax1.set_ylim(1E-34, 1E-26)
ax1.legend(loc=3)
txt = ('B = %.0E G\n'+r'$\gamma_{max}$ = %.0E'+'\n'+r'$\gamma_{min}$ = %.0E') % (self.B, self.gmax, self.gmin)
ax1.annotate(txt, (0.8,0.8), xycoords='axes fraction')
ax1.xaxis.set_ticks_position('top')
ax1.xaxis.set_label_position('top')
ax1.set_xlabel(r'$\nu$ (Hz)')
ax1.set_ylabel(r'$J(\nu)\,(\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{Hz}^{-1})$')
ax1.grid()
res_scale = 1E-28
ax2=fig.add_axes((.1,.25,.8,.15))
ax2.plot(nus, (f(nus, par1[0], par1[1], par1[2])-Js)/res_scale, c='g', ls=':')
ax2.plot(nus, (f(nus, par2[0], par2[1], par2[2])-Js)/res_scale, c='b', ls='--')
ax2.plot(nus, (f(nus, par3[0], par3[1], par3[2])-Js)/res_scale, c='r', ls='-.')
ax2.semilogx()
ax2.set_xlim(1E-5, 1E13)
ax2.set_ylabel('residue\n(%.0e)' % res_scale)
ax2.set_ylim(-5, 5)
ax2.set_yticks(np.array([-4, -2, 0., 2, 4]))
ax2.set_yticklabels([-4, -2, 0., 2])
ax2.hlines(0, nus.min(), nus.max(), lw=0.2, alpha=0.5)
ax2.set_xticklabels([])
ax3=fig.add_axes((.1,.1,.8,.15))
ax3.plot(nus, (f(nus, par1[0], par1[1], par1[2])-Js)/Js*100, c='g', ls=':')
ax3.plot(nus, (f(nus, par2[0], par2[1], par2[2])-Js)/Js*100, c='b', ls='--')
ax3.plot(nus, (f(nus, par3[0], par3[1], par3[2])-Js)/Js*100, c='r', ls='-.')
#plt.semilogy()
ax3.semilogx()
ax3.set_xlim(1E-5, 1E13)
ax3.set_ylabel('rel. error (%)')
ax3.set_ylim(-100,100)
ax3.set_yticks([-100, -50, 0, 50, 100])
ax3.set_yticklabels([-100, -50, 0, 50])
ax3.hlines(0, nus.min(), nus.max(), lw=0.2, alpha=0.5)
if pol == 'para':
fname = './synchrotron_spectra_fitting/Para_B%i_gmax%i_gmin%i.png'
elif pol == 'perp':
fname = './synchrotron_spectra_fitting/Perp_B%i_gmax%i_gmin%i.png'
plt.savefig(fname % (np.log10(self.B), np.log10(self.gmax), np.log10(self.gmin)))
if show:
plt.show()
return fig
In [ ]:
for pol in ['perp', 'para']:
for B in [1E-5, 1E-6]:
for gmin in [1E0, 1E1]:
for gmax in [1E2, 1E5, 1E8]:
sync = Synchrotron(B, gmax, gmin, pol=pol)
sync.fit_sync()
fig = sync.plot_sync(show=True)
In [ ]:
p = 2.
# Ribicki & Lightman eq. (6.36), mu = (p-3)/2
f = 2**((p+1)/2)/(p+1)*gamma_func(p/4 + 19./12)*gamma_func(p/4 - 1./12)
g = 2**((p-3)/2)*gamma_func(p/4 + 7./12)*gamma_func(p/4 - 1./12)
print (f, g)
print (g/f)
print (2*f)
perp_const = f+g
para_const = f-g
print (perp_const, para_const)
#print 1./(2.*np.pi)**1.5/(p+1)*gamma_func(p/4 + 19./12)*gamma_func(p/4 - 1./12)
#print 1./2./(4.*np.pi)**0.5*5.8
In [ ]:
# Ribicki & Lightman eq. (6.36)
mu = -0.5
f = 2**(mu+1)/(mu+2)*gamma_func(mu/2+7/3)*gamma_func(mu/2+2/3)
g = 2**mu*gamma_func(mu/2+4/3)*gamma_func(mu/2+2/3)
print (f, g)
In [ ]:
xx=np.logspace(-6, 2, 100)
xmax = 1E2
Js = np.array([quad(func, x, xmax)[0] for x in xx])
#print(Js)
plt.plot(xx, [func(x) for x in xx], label=u'$x^{-1/2}F(x)$')
plt.plot(xx, f*np.exp(-xx), label=u'$e^{-x}$')
plt.plot(xx, Js, label=u'$\int_x^{\infty} y^{-1/2}F(y) dy$')
plt.semilogx()
plt.semilogy()
plt.ylim(1E-10, 1E1)
plt.xlim(1E-1,1E2)
plt.legend()
In [ ]:
print ('B = %.0e, gmax = %.0e' % (sync.B, sync.gmax))
print ('(nu_min, nu_max) =', sync.nus.min(), sync.nus.max())
print (sync.par1)
print (sync.par2)
print (sync.par3)
In [ ]:
#for x << 1, F(x) ~ x**(1./3.)
print (4*np.pi/np.sqrt(3)/gamma_func(1./3.)/2.**(1./3.))
#for x >> 1, F(x) ~ e**(-x)*x**(1./2)
print (np.sqrt(np.pi/2.))
In [ ]: