line of sight intensity image calculator

def calc_emissivity(point):

    x = point[0]
    y = point[1]
    z = point[2]

    if x**2 + y**2 + z**2 < 0.7**2:
        return 1
    else:
        return 0

def calc_emissivity_2(point):

    x = float(point[0])
    y = float(point[1])
    z = float(point[2][0])

    r2 = x**2 + y**2 + z**2
    if r2< 0.7**2:
        return r2**2
    else:
        return 0

def los_integrate(emiss_func, x, y, zmin, zmax):

    import scipy.integrate

    return (scipy.integrate.quadrature
            (lambda z: emiss_func([x,y,z]),
             zmin, zmax)[0])

class RichBremsstrahlung:

    def __init__(self, fname, inclination=0):

        import h5py
        import numpy
        import scipy.interpolate

        self.q = inclination
        with h5py.File(fname,'r+') as f:
            x_list = numpy.array(f['x_coordinate'])
            y_list = numpy.array(f['y_coordinate'])
            t_list = numpy.array(f['temperature'])
            d_list = numpy.array(f['density'])
        self.z_min = numpy.min(y_list)
        self.z_max = numpy.max(y_list)
        self.temperature_func = (
            scipy.interpolate.LinearNDInterpolator
            (numpy.array([x_list, y_list]).T,
             t_list,
             fill_value=0))
        self.density_func = (
            scipy.interpolate.LinearNDInterpolator
            (numpy.array([x_list,y_list]).T,
             d_list,
             fill_value=0))

    def __call__(self, p):

        import numpy

        x = p[0]
        y = p[1]*numpy.cos(self.q)+p[2]*numpy.sin(self.q)
        z = p[2]*numpy.cos(self.q)-p[1]*numpy.sin(self.q)

        r = numpy.sqrt(x**2 + y**2)
        t = self.temperature_func(r,z)
        d = self.temperature_func(r,z)
        return numpy.sqrt(t)*d**2

def main():

    import numpy
    import pylab
    import matplotlib as plt
    import scipy.integrate

    q = numpy.pi/3
    x_list = numpy.linspace(-10,10,num=100)
    y_list = numpy.linspace(-10,10,num=101)
    rbs = RichBremsstrahlung('snapshot_575.h5',q)
    I_list = [[los_integrate(rbs,x,y,rbs.z_min,rbs.z_max)
               for x in x_list]
              for y in y_list]
    X_list, Y_list = numpy.meshgrid(x_list, y_list)
    pylab.contourf(X_list, Y_list, I_list, 20,
                   cmap=plt.cm.bone)
    pylab.axis('equal')
    pylab.xlabel('x [pc]')
    pylab.ylabel('y [pc]')
    pylab.title(r'$i = '+str(180*q/numpy.pi)+'^O$')
    pylab.show()

if __name__ == '__main__':

    main()