Development of Synchrotron model and fitting for MWA data

By Kamen Kozarev



In [3]:

    
from matplotlib import pyplot as plt
%matplotlib inline
import math
import scipy.integrate as integrate
import numpy as np
import scipy.special as special

Formulae for asymptotic synchrotron spectrum values



In [97]:

    
#Asymptotic synchrotron values
x=np.arange(1000)/250.
f=np.zeros(1000)
f[0:200]=4*math.pi/(math.sqrt(3)*special.gamma(1./3.))*pow(x[0:200]/2.,1./3.)
f[201:]=math.sqrt(math.pi/2.)*np.exp(-x[201:])*pow(x[201:],1./2.)
plt.plot(x,f)









    Out[97]:





[<matplotlib.lines.Line2D at 0x1513e15810>]

Algorithm for the F(x) function, which approximates the synchrotron spectral shape.



In [1]:

    
#The F(x) function
def fnu(values):
    result=[]
    for x in values:
        integral = integrate.quad(lambda x: special.kv(5./3.,x), x, np.inf)
        result.append(x * integral[0])
    return np.array(result)

def fnu_single(x):
    integral = integrate.quad(lambda x: special.kv(5./3.,x), x, np.inf)
    result= x * integral[0]
    return np.array(result)

A generic model for testing the algorithm



In [9]:

    
#Define the constants and arrays
nu_c=0.5
nfreq=100
x=np.arange(nfreq)/(nfreq/4.)
Ptot=np.zeros(nfreq)
p=3.
#Pitch angle
a=45*math.pi/180.

FF=fnu(x)
plt.plot(x,FF)









    Out[9]:





[<matplotlib.lines.Line2D at 0x110049690>]

A more realistic model, using proper frequencies.



In [45]:

    
#start and end frequency, MHz
nu_start=1.e6
nu_end=300.e6
#Critical frequency, MHz
nu_c=120.e6
#electron mass, g
me=9.10938e-28
#electron charge Statcoulomb
q=4.80326e-10
#speed of light [cm/s]
c=2.99792458e10
#Bmag [Gauss]
bmag=10.

#Pitch angle
a=45.*math.pi/180.
#Frequency array, MHz
nus=np.linspace(nu_start,nu_end,num=nfreq)



In [43]:

    
factors=math.sqrt(3)*q**3*bmag*math.sin(a)/(2.*math.pi*me*c**2)
#Array to hold total power
Ptot=[] 
for nu in nus:
    ff=fnu_single(nu/nu_c)
    ff.shape
    Ptot.append(factors*ff)
Ptot=np.array(Ptot)
plt.plot(nus/1.e6,Ptot)
plt.xlim(1,300)
plt.title("Electron Synchrotron spectrum")
plt.xscale("log")
plt.xlabel("f [MHz]")
plt.yscale("log")
plt.ylabel("Power")









    Out[43]:





<matplotlib.text.Text at 0x115583f10>

The Synchrotron spectrum for a power-law electron distribution, including self-absorption



In [189]:

    
#CONSTANTS:
#electron mass, g
me=9.10938e-28
#electron charge Statcoulomb
q=4.80326e-10
#speed of light [cm/s]
c=2.99792458e10

#PARAMETERS:
#Constant of proportionality
const=10.
#Power law index of electron spectrum
p=3.5
#Pitch angle
a=80.*math.pi/180.
#Bmag [Gauss]
bmag=10.
#Critical frequency, MHz
nu_c=100.e6

#start and end frequency, MHz
nu_start=80.e6
nu_end=300.e6
#Frequency array, MHz
nus=np.linspace(nu_start,nu_end,num=nfreq)#/(nfreq/4.)

#Array to hold total power
Ppow=[]
Ssyn=[]

jfactor1=math.sqrt(3)*q**3*const*bmag*math.sin(a)/(2.*math.pi*me*c**2*(p+1))
jfactor2=special.gamma(p/4. + 19./12.)*special.gamma(p/4. - 1./12.)

#afactor1=math.sqrt(3)*q**3/(8*math.pi*me) * const*(bmag*math.sin(a))**(0.5*p+1)
#afactor2=special.gamma((3.*p+2.)/12.)*special.gamma((3*p+22.)/12.)

for nu in nus:
    jfactor3=pow(((me*c*(nu/nu_c)*2*math.pi)/(3*q*bmag*math.sin(a))),(0.5*(1-p)))
    Pemit=jfactor1*jfactor2*jfactor3
    #Abs=afactor1*afactor2*pow(nu,-0.5*(p+4.))
    #Ppow.append(Pemit/(Abs*4*math.pi))
    Ppow.append(jfactor3*pow(nu/nu_c,5./2))
    #Ssyn.append(pow(nu/nu_c,5./2)*(1-math.exp(-pow(nu/nu_c,-0.5*(p+4)))))
    Ssyn.append(pow(nu/nu_c,5./2)*(1-math.exp(-pow(nu/nu_c,-0.5*(p+3)))))
Ppow=np.array(Ppow)
Ssyn=np.array(Ssyn)

plt.plot(nus/1.e6,Ssyn)
#plt.plot(nus/1.e6,Ppow)
#plt.xlim(30,300)
plt.title("Electron Power Law Synchrotron spectrum")
plt.xscale("log")
plt.xlabel("f [MHz]")
plt.yscale("log")
plt.ylabel("Power")









    Out[189]:





<matplotlib.text.Text at 0x152ec190d0>



In [ ]: