Read MDM Spectrum

read_mdm.py

def read_mdm(filename,to_plot=False,return_header=False):
    with fits.open(filename) as spec:
#             print(spec.info())
#             print(spec[0].data)
        
            flux = spec[0].data[0][0]

            w0 = np.float(spec[0].header["CRVAL1"])
            wi = np.int(spec[0].header["CRPIX1"])
            wstep = np.float(spec[0].header["CD1_1"])
            w00 = (0-wi)*wstep + w0
            lf = len(flux)
            wavelength = np.arange(w00,w00+(lf-1)*wstep,wstep)
            while len(wavelength)<lf:
                wavelength = np.append(wavelength,
                                       wavelength[-1]+wstep)

            var = spec[0].data[3][0]
#             print(len(wavelength), len(flux), len(var))

            print("sec(z) =",spec[0].header["AIRMASS"],
                  "HA =",spec[0].header["HA"])
            
            if to_plot is True:
                plt.figure(figsize=(10,7))
                for i in range(4):
                    plt.plot(wavelength,spec[0].data[i][0],
                             label="row {0}".format(i))
                    plt.yscale("log")
                    plt.ylim(1e-16,1e-12)
                    plt.xlim(6000,8000)
                    plt.legend(loc="best")
                    
            if return_header is True:
                header = spec[0].header
                
    if return_header is True:
        return wavelength, flux, var, header
    else:
        return wavelength, flux, var