Closed thomasorb closed 6 years ago
This script proves that it works even with velocity calibration
from orcs.process import SpectralCube
import numpy as np
import pylab as pl
import orb.utils.io as io
ix, iy = 1878, 1363
vel = -60
cube = SpectralCube('/home/thomas/M31_SN3.merged.cm1.1.0.hdf5',debug=False)
skymap = np.ones((cube.dimx, cube.dimy), dtype=float) * -200 # fake wavelength correction of -200 km/s
vel += skymap[ix,iy] # initial velocity is corrected for wavelength correction
io.write_fits('skymap.fits', skymap, overwrite=True)
cube.correct_wavelength('skymap.fits')
axis, sky_spectrum = cube.extract_spectrum(500,500,30, mean_flux=True)
#sky_spectrum = None # uncommenting this line deactivates the sky spectrum subtraction
axis, spectrum, fit = cube.fit_lines_in_spectrum(ix, iy, 0, ('Halpha', '[NII]6583'),
fmodel='sinc',
pos_def='1',
pos_cov=vel,
snr_guess=None,
nofilter=True,
subtract_spectrum=sky_spectrum)
print fit['flux_gvar']
print fit['velocity_gvar']
mask = np.zeros((cube.dimx, cube.dimy), dtype=bool)
mask[ix-1:ix+2, iy-1:iy+2] = True
reg = np.nonzero(mask)
print reg
cube.fit_lines_in_region(reg, ('Halpha', '[NII]6583'),
fmodel='sinc',
pos_def='1',
pos_cov=vel,
snr_guess=None,
nofilter=True,
subtract_spectrum=sky_spectrum)
fluxmap = io.read_fits('./M31_SN3.1.0.ORCS/MAPS/M31_SN3.1.0.LineMaps.map.6563.1x1.flux.fits')
velmap = io.read_fits('./M31_SN3.1.0.ORCS/MAPS/M31_SN3.1.0.LineMaps.map.6563.1x1.velocity.fits')
# The results below should be exactly the same for a given set of parameters but different when sky
# subtraction is turned on or off
print 'fit lines in region map results : ----------------'
print fluxmap[ix, iy]
print velmap[ix, iy]
print 'single pixel fit results: ------------------'
print fit['flux'][0] # should give the exact same value as the region fit
print fit['velocity'][0] # should give the exact same value as the region fit
bug fixed with https://github.com/thomasorb/orcs/commit/6fd22a86f4d8b7eb1a9fc142b64bc02b8d8744a9