In [1]:

    

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(context="talk")

from bokeh import plotting
plotting.output_notebook()

import pandas as pd

from astropy.modeling import models, fitting


    

    





    
        
        
        
    

        
        BokehJS successfully loaded.
    

In [84]:

    

nstep = 10000


    

In [3]:

    

%%bash -s "$nstep"
make
./sedov3 $1


    

    




make: `sedov3' is up to date.
        1000
 spherical_standard_omega0p00_nstep_1000.dat                                     
 standard
 xgeom =  3.000000E+00 eblast=  1.000000E+00 omega =  0.000000E+00 alpha =  4.935902E-01 j1    =  2.221251E-02 j2    =  1.878160E-02 
 r2    =  1.151666E+00 rho2  =  4.000000E+00 u2    =  3.454999E-01 e2    =  5.968510E-02 p2    =  1.591603E-01 cs2   =  2.575204E-01 

In [4]:

    

i, x, den, energy, pressure, velocity, cs = np.loadtxt('spherical_standard_omega0p00_nstep_' + str(nstep).zfill(4) + '.dat', skiprows=2, unpack=True)


    

In [85]:

    

filename = 'spherical_standard_omega0p00_nstep_' + str(nstep).zfill(5) + '.dat'

df = pd.read_csv(filename, skiprows=1, sep=r"\s+", index_col=0, usecols=np.arange(6)+1)
(df / df.max()).plot()


    

    
Out[85]:





<matplotlib.axes._subplots.AxesSubplot at 0x10ec86ba8>






    

In [81]:

    

power_law_model = models.PowerLaw1D(alpha=-4.5)
fitter = fitting.LevMarLSQFitter()

fitting_range = (df.index<.6)
fitted = fitter(power_law_model, df.index[fitting_range], df.den[fitting_range])
print(fitted)


    

    




Model: PowerLaw1D
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Parameters:
      amplitude        x_0          alpha     
    ------------- ------------- --------------
    2.07893518303 1.25337004394 -4.53217010146

In [82]:

    

plt.plot(df.index, df.den, label="data")
plt.plot(df.index, fitted(df.index), label="fit")
plt.xscale("log")
plt.yscale("log")
plt.legend(loc="best")


    

    
Out[82]:





<matplotlib.legend.Legend at 0x10f314cc0>






    

In [83]:

    

p = plotting.figure(x_axis_type = "log", y_axis_type = "log", tools="crosshair, wheel_zoom")

p.line(df.index, df.den, legend="density: data" )

fit_by_hand = models.PowerLaw1D(amplitude=5.199e-4, x_0 =.2, alpha=-4.5)
p.line(df.index, fit_by_hand(df.index), legend="density: fit")

plotting.show(p)