maps.py 文件源码

def make_emission_measure_map(time: u.s, field, instr, temperature_bin_edges=None, **kwargs):
    """
    Return a cube of maps showing the true emission meausure in each pixel
    as a function of electron temperature.
    """
    plot_settings = {'cmap': cm.get_cmap('magma'),
                     'norm': colors.SymLogNorm(1, vmin=1e25, vmax=1e29)}
    plot_settings.update(kwargs.get('plot_settings', {}))

    # read unbinned temperature and density
    with h5py.File(instr.counts_file, 'r') as hf:
        try:
            i_time = np.where(np.array(hf['time'])*u.Unit(hf['time'].attrs['units']) == time)[0][0]
        except IndexError:
            raise IndexError('{} is not a valid time in observing time for {}'.format(time, instr.name))
        unbinned_temperature = np.array(hf['electron_temperature'][i_time,:])
        temperature_unit = u.Unit(hf['electron_temperature'].attrs['units'])
        unbinned_density = np.array(hf['density'][i_time,:])
        density_unit = u.Unit(hf['density'].attrs['units'])

    # setup bin edges and weights
    if temperature_bin_edges is None:
        temperature_bin_edges = 10.**(np.arange(5.5, 7.5, 0.1))*u.K
    x_bin_edges = np.diff(instr.bin_range.x)/instr.bins.x*np.arange(instr.bins.x+1) + instr.bin_range.x[0]
    y_bin_edges = np.diff(instr.bin_range.y)/instr.bins.y*np.arange(instr.bins.y+1) + instr.bin_range.y[0]
    z_bin_edges = np.diff(instr.bin_range.z)/instr.bins.z*np.arange(instr.bins.z+1) + instr.bin_range.z[0]
    z_bin_indices = np.digitize(instr.total_coordinates.value[:,2], z_bin_edges, right=True)
    dh = np.diff(z_bin_edges)[z_bin_indices - 1]
    emission_measure_weights = (unbinned_density**2)*dh
    # bin in x,y,T space with emission measure weights
    xyT_coordinates = np.append(instr.total_coordinates.value[:,:2], 
                                unbinned_temperature[:,np.newaxis], axis=1)
    hist, _ = np.histogramdd(xyT_coordinates, 
                             bins=[x_bin_edges, y_bin_edges, temperature_bin_edges.value], 
                             weights=emission_measure_weights)

    meta_base = instr.make_fits_header(field, instr.channels[0])
    del meta_base['wavelnth']
    del meta_base['waveunit']
    meta_base['detector'] = r'$\mathrm{EM}(T)$'
    meta_base['comment'] = 'LOS Emission Measure distribution'
    data = np.transpose(hist, (1,0,2))*density_unit*density_unit*instr.total_coordinates.unit

    return EMCube(data, meta_base, temperature_bin_edges, plot_settings=plot_settings)