In [ ]:

    

#%pdb
%matplotlib inline
import matplotlib
matplotlib.rcParams['font.family'] = 'stixgeneral'
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt
import os
import yt
yt.mylog.setLevel("ERROR")
import numpy as np
from synchrotron.yt_synchrotron_emissivity import StokesFieldName, synchrotron_fits_filename
import pyfits


    

In [ ]:

    

dir = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/'
ds = yt.load(dir+'data/MHD_Jet_10Myr_hdf5_plt_cnt_0100')
ptype = 'lobe'
proj_axis = 'x'
fitsname = synchrotron_fits_filename(ds, dir, ptype, proj_axis)
print(fitsname)
#hdulist = pyfits.open(fitsname)
times = {100:'1.6 Myr', 200:'3.2 Myr', 600:'9.5 Myr',
         910:'57 Myr', 1050:'101 Myr', 1380:'300 Myr', 1700:'500 Myr'}


    

In [ ]:

    

nus = np.array([100,150,300,600,1400,8000])
fitsname = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/cos_synchrotron_QU_nn_lobe/fits/synchrotron_lobe_%04d.fits'
for fno in [100, 200, 600, 910, 1050, 1380, 1700]:
    hdulist = pyfits.open(fitsname % fno)
    fluxes = []
    for nu in [(nu, 'MHz') for nu in nus]:
        stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
        field = stokes.I[1]
        fluxes.append(np.sum(hdulist[field].data))
    plt.plot(nus, fluxes, 'x-', label=times[fno])
    if fno in [1700]:
        powerlaw = fluxes[0]*(nus/100)**(-0.5)
        plt.plot(nus, powerlaw, ':', lw=1, label=r'$\nu^{-0.5}$')
        powerlaw = fluxes[0]*(nus/100)**(-1)
        plt.plot(nus, powerlaw, ':', lw=1, label=r'$\nu^{-1}$')

plt.xlabel('Frequency (MHz)')
plt.semilogx()
plt.ylabel('Flux')
plt.semilogy()
plt.legend()
plt.grid()
plt.show()


    

In [ ]:

    

nus = np.array([100,150,300,600,1400,8000])
dirs = ['fits_linear', 'fits_log', 'fits_emis', 'fits_gamc']

fitsname = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/cos_synchrotron_QU_nn_lobe/%s/synchrotron_lobe_0910.fits'
for dir in dirs:
    hdulist = pyfits.open(fitsname % dir)
    fluxes = []
    for nu in [(nu, 'MHz') for nu in nus]:
        stokes = StokesFieldName(ptype, nu, proj_axis, field_type='flash')
        field = stokes.I[1]
        mask = hdulist[field].data < 1.0E0
        fluxes.append(np.sum(hdulist[field].data[mask]))
    fluxes = np.array(fluxes)
    fluxes /= fluxes[0]
    plt.plot(nus, fluxes, 'x-', label=dir.split('_')[-1])

powerlaw = fluxes[0]*(nus/100)**(-0.5)
plt.plot(nus, powerlaw, '--', lw=1, label=r'$\nu^{-0.5}$')
powerlaw = fluxes[0]*(nus/100)**(-1)
plt.plot(nus, powerlaw, '-.', lw=1, label=r'$\nu^{-1}$')
powerlaw = fluxes[0]*(nus/100)**(-1.5)
plt.plot(nus, powerlaw, ':', lw=1, label=r'$\nu^{-1.5}$')

plt.xlabel('Frequency (MHz)')
plt.semilogx()
plt.ylabel('Flux')
plt.semilogy()
plt.legend()
plt.grid()
plt.show()


    

In [ ]:

    

nus = np.array([100,150,300,600,1400,8000])
fitsname = '/d/d5/ychen/2015_production_runs/1022_h1_10Myr/cos_synchrotron_QU_nn_lobe/fits_linear/synchrotron_lobe_1_0_2_%04d.fits'
for fno in [100, 200, 600, 910, 1050, 1380, 1700]:
    hdulist = pyfits.open(fitsname % fno)
    fluxes = []
    for nu in [(nu, 'MHz') for nu in nus]:
        stokes = StokesFieldName(ptype, nu, [1,0,2], field_type='flash')
        field = stokes.I[1]
        fluxes.append(np.sum(hdulist[field].data))
    fluxes = np.array(fluxes)
    fluxes /= fluxes[0]
    plt.plot(nus, fluxes, 'x-', label=times[fno])
    if fno in [1700]:
        powerlaw = fluxes[0]*(nus/100)**(-0.5)
        plt.plot(nus, powerlaw, ':', lw=1, label=r'$\nu^{-0.5}$')
        powerlaw = fluxes[0]*(nus/100)**(-1)
        plt.plot(nus, powerlaw, ':', lw=1, label=r'$\nu^{-1}$')

plt.xlabel('Frequency (MHz)')
plt.semilogx()
plt.ylabel('Flux')
plt.semilogy()
plt.legend()
plt.grid()
plt.show()


    

In [ ]: