capricorn_breakout_analysis

def whitham_mu(g):

    import math

    return 1./(2+math.sqrt(2*g/(g-1)))

def theoretical_density_powerlaw_index(omega, gamma):

    mu = whitham_mu(gamma)

    return (1.+omega*mu+omega)/(omega*mu)

def consolidate(fname):

    import h5py
    import numpy

    with h5py.File(fname) as f:
        res = dict((key,numpy.array(f[key])) for key in f)
    res['centres'] = 0.5*(res['grid'][1:]+res['grid'][:-1])
    return res

def main():

    import numpy
    import pylab

    initial = consolidate('temp_breakout/initial.h5')
    final = consolidate('temp_breakout/final.h5')

    #cd_pos = final['grid'][numpy.argmax(1/numpy.array(initial['grid']))]

    filtered = dict((key,numpy.array([v for x,v in zip(final['centres'],
                                                       final[key])
                                      if x>0.1*numpy.max(final['centres'])]))
                    for key in ['centres','pressure','velocity','density'])
    fit_data = numpy.polyfit(filtered['centres'],
                             filtered['velocity'], 1)
    d_fit = numpy.polyfit(numpy.log(filtered['centres']),
                        numpy.log(filtered['density']),1)
    print(d_fit)
    print('expected '+str(theoretical_density_powerlaw_index(1.5,5./3.)))

    for n, var_name in enumerate(['density','pressure','velocity']):
        pylab.subplot(3,1,n+1)
        if var_name == 'velocity' or var_name == 'pressure':
            #pylab.plot(initial['centres'], initial[var_name])
            pylab.plot(final['centres'], final[var_name])
        else:
            #pylab.plot(initial['centres'], numpy.log(initial[var_name]))
            pylab.semilogy(final['centres'], final[var_name])
    pylab.plot(final['centres'], numpy.polyval(fit_data, final['centres']))
    pylab.subplot(311)
    pylab.semilogy(filtered['centres'],
                   numpy.exp(numpy.polyval(d_fit,numpy.log(filtered['centres']))))
    pylab.show()

if __name__ == '__main__':

    main()