offset_supernova_gravity_sink_analysis

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_triangle_area(x1,x2,x3,y1,y2,y3):

    dx_a = x2 - x1
    dx_b = x3 - x1
    dy_a = y2 - y1
    dy_b = y3 - y1
    return abs(0.5*(dx_a*dy_b-dx_b*dy_a))

def calc_polygon_area(x_list, y_list):

    res = 0
    for i in range(2,len(x_list)):
        res += calc_triangle_area(x_list[0],
                                  x_list[i-1],
                                  x_list[i],
                                  y_list[0],
                                  y_list[i-1],
                                  y_list[i])
    return res

def calc_mass(fname):

    import numpy

    data = consolidate(fname)
    vert_x_raw = data['x position of vertices']
    vert_y_raw = data['y position of vertices']
    vert_n_raw = data['Number of vertices in cell']
    vert_idx_list = numpy.concatenate([[0],numpy.cumsum(vert_n_raw)])
    areas = []
    for i in range(len(data['Number of vertices in cell'])):
        upbound = vert_idx_list[i+1]
        lowbound = vert_idx_list[i]
        areas.append(calc_polygon_area(vert_x_raw[lowbound:upbound],
                                       vert_y_raw[lowbound:upbound]))
    areas = numpy.array(areas)
    masses = areas*data['density']

    return sum([m for x,y,m in zip(data['x_coordinate'],
                                   data['y_coordinate'],
                                   masses) if x**2+y**2<500**2])

def mass_history():

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

    time_list, mass_list, x_mom_list, y_mom_list, energy_list, tracer = numpy.loadtxt('conserved.txt',unpack=True)
    pylab.plot(time_list, mass_list[0] - mass_list)
    pylab.xlabel('Time [year]')
    pylab.ylabel(r'Mass [$M_{\odot}$]')
    pylab.show()

def main():

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

    data = consolidate('snapshot_60.h5')
    g = 5./3.
    data['abs_velocity'] = numpy.sqrt(data['x_velocity']**2+
                                      data['y_velocity']**2)
    data['sound_speed'] = numpy.sqrt(g*data['pressure']/data['density'])
    data['mach_number'] = data['abs_velocity']/data['sound_speed']
    pylab.tricontourf(data['x_coordinate'],
                      data['y_coordinate'],
                      data['y_velocity'],50)
    pylab.colorbar()
    pylab.axis('equal')
    pylab.show()

if __name__ == '__main__':

    mass_history()