planar_relativistic_breakout_analysis

def main():

    import matplotlib
    matplotlib.use('Qt4Agg')
    import pylab
    import h5py
    import numpy

    initial = h5py.File('initial.h5')
    final = h5py.File('final.h5')

    for n, key in enumerate(['density','pressure','celerity']):
        pylab.subplot(3,1,n+1)
        pylab.plot(initial['position'],initial[key])
        pylab.plot(final['position'],final[key])
    pylab.subplot(311)
    pylab.title('t = '+str(numpy.array(final['time'])[0]))
    pylab.show()

def asymptotic():

    import matplotlib
    matplotlib.use('Qt4Agg')
    import pylab
    import h5py
    import numpy
    import math

    final = h5py.File('final.h5')
    initial = h5py.File('initial.h5')

    gamma_list = numpy.sqrt(numpy.array(final['celerity'])**2+1)
    filtered = {}
    filtered['position'] = [x for x,w in zip(final['position'],
                                             gamma_list)
                            if w>5 and x<2*numpy.max(initial['position'])]
    assert(len(filtered['position'])>0)
    filtered['lorentz'] = [w for x,w in zip(final['position'],
                                             gamma_list)
                           if w>5 and x<2*numpy.max(initial['position'])]
    assert(numpy.max(filtered['position'])<0)
    fit_data = numpy.polyfit(numpy.log(-numpy.array(filtered['position'])),
                             numpy.log(filtered['lorentz']),1)
    pylab.loglog(-numpy.array(final['position']),
                 gamma_list)
    pylab.loglog(-numpy.array(final['position']),
                 numpy.exp(numpy.polyval(fit_data,
                                         numpy.log(-numpy.array(final['position'])))))
    print(fit_data)
    pylab.show()

def consolidate(fname):

    import h5py
    import numpy

    with h5py.File(fname) as f:
        return dict((key,numpy.array(f[key])) for key in f)

def calc_celerity_history():

    import glob
    import h5py
    import numpy

    all_snapshot_files = glob.glob('snapshot_*.h5')
    all_snapshot_files = sorted(all_snapshot_files,
                                key=lambda x:int(x.replace('snapshot_','').replace('.h5','')))

    celerity_histories = []
    time_list = []
    for fname in all_snapshot_files:
        print(fname)
        with h5py.File(fname) as f:
            celerity_histories.append(numpy.array(f['celerity']))
            time_list.append(numpy.array(f['time'])[0])
    celerity_histories = numpy.array(celerity_histories)
    return time_list, celerity_histories

def midvals(ar):

    import numpy

    return 0.5*(numpy.array(ar[1:])+numpy.array(ar[:-1]))

def velocity_histories():

    import matplotlib
    matplotlib.use('Qt4Agg')
    import pylab
    import numpy
    from joblib import Memory

    mem = Memory(cachedir='/home/almog/joblib')
    cached_calc_celerity_history = mem.cache(calc_celerity_history)
    time_list, celerity_histories = cached_calc_celerity_history()
    pylab.subplot(211)
    pylab.plot(time_list, celerity_histories[:,500])
    pylab.subplot(212)
    pylab.plot(midvals(time_list), numpy.diff(celerity_histories[:,500])/numpy.diff(time_list))
    pylab.show()

def asymptotic_2():

    import matplotlib
    matplotlib.use('Qt4Agg')
    import pylab
    import h5py
    import numpy
    import math
    import glob

    mid = consolidate('mid.h5')
    finals = [consolidate(fname) for fname in glob.glob('milestone_*.h5')]
    mid['lorentz factor'] = numpy.sqrt(mid['celerity']**2+1)
    for itm in finals:
        itm['lorentz factor'] = numpy.sqrt(itm['celerity']**2+1)
        itm['x_list'] = 0.5*(numpy.log10(itm['lorentz factor'])[1:]+
                             numpy.log10(itm['lorentz factor'])[:-1])
        itm['y_list'] = (numpy.diff(numpy.log10(mid['lorentz factor']))/
                         numpy.diff(numpy.log10(itm['lorentz factor'])))

    for itm in finals:
        pylab.subplot(211)
        pylab.loglog(itm['lorentz factor'], mid['lorentz factor'],
                     basex=2, basey=2)
        pylab.subplot(212)
        pylab.plot(itm['x_list'], itm['y_list'])
    pylab.plot(itm['x_list'], [1./(1.+numpy.sqrt(3.)) for x in itm['x_list']],'k')
    pylab.ylim((0,1))
    pylab.show()

def show_single_snapshot():

    import matplotlib
    matplotlib.use('Qt4Agg')
    import pylab
    import h5py
    import numpy

    data = consolidate('mid.h5')

    for n, key in enumerate(['density','pressure','celerity']):
        pylab.subplot(3,1,1+n)
        pylab.plot(data['position'],data[key])
    pylab.subplot(311)
    pylab.title('t = '+str(data['time'][0]))
    pylab.show()

def asymptotic_3():

    import matplotlib
    matplotlib.use('Qt4Agg')
    import pylab
    import h5py
    import numpy
    import math
    import glob

    mid = consolidate('mid.h5')
    final = consolidate('milestone_10000.h5')
    mid['lorentz_factor'] = numpy.sqrt(mid['celerity']**2+1)
    final['lorentz_factor'] = numpy.sqrt(final['celerity']**2+1)
    #x_list = [numpy.log(gf) for gf,gb in zip(final['lorentz_factor']
    x_list = [numpy.log(gf) for gf,gb in zip(final['lorentz_factor'],
                                             mid['lorentz_factor'])
              if 60>gf>2]
    y_list = [numpy.log(gb) for gf,gb in zip(final['lorentz_factor'],
                                             mid['lorentz_factor'])
              if 60>gf>2]
    fit_data = numpy.polyfit(x_list, y_list, 1)
    print(fit_data)
    print(1./(1+numpy.sqrt(3.)))
    pylab.loglog(final['lorentz_factor'],
                 mid['lorentz_factor'],
                 basex=2, basey=2)
    pylab.loglog(final['lorentz_factor'],
                 numpy.exp(fit_data[1])*final['lorentz_factor']**fit_data[0])
    pylab.show()

if __name__ == '__main__':

    #main()
    asymptotic_3()
    #show_single_snapshot()