Commit 448d192a authored by bepinat's avatar bepinat

bug in NII line fit

parent 1937494c
......@@ -25,7 +25,7 @@ import copy
import numpy as np
import pyfits as pf
from mpfit import mpfit # python 2.7...
#import ipdb
import ipdb
#from matplotlib import pyplot as plt
from scipy import ndimage
from scipy import constants as ct
......@@ -359,8 +359,9 @@ class lines:
self.append(line('HGAMMA', 4340., ref='HBETA', low=0.44, up=0.5, th=0.468))
self.append(line('HDELTA', 4101., ref='HBETA', low=0.23, up=0.29, th=0.259))
self.append(line('HEPS', 3968., ref='HBETA', low=0.13, up=0.19, th=0.159))
self.append(line('NII6583', 6583., ref='NII6548', low=2.7, up=3.3, th=3.))
self.append(line('NII6548', 6548.))
self.append(line('NII6584', 6583., ref='NII6545', low=2.7, up=3.3, th=3.))
#self.append(line('NII6583', 6583., ref='NII6545', low=2.7, up=3.3, th=3.))
self.append(line('NII6545', 6545.))
self.append(line('SII6731', 6731.))
self.append(line('SII6717', 6717., ref='SII6731', low=0.45, up=1.45, th=1.))
self.append(line('OIII5007', 5007., ref='OIII4959', low=2.7, up=3.3, th=3.))
......@@ -447,7 +448,6 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
linefit = np.array([els.lines[i].fit for i in els.lines])[argsorted]
linewave = np.array([els.lines[i].wave for i in els.lines])[argsorted]
linename = np.array([els.lines[i].name for i in els.lines])[argsorted]
#ipdb.set_trace()
indok = np.zeros(0, dtype = int)
firstind = np.zeros(linefit.size)
lastind = np.zeros(linefit.size)
......@@ -572,8 +572,10 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
parinfo[2 * els.index + 1 + i]['limited'][1] = 1
parinfo[2 * els.index + 1 + i]['limits'][1] = els.lines[name].up
if els.lines[name].th != None:
parinfo[2 * els.index + 1 + i]['value'] = params[2 * els.index + 1 + i] = els.lines[name].th
parinfo[2 * els.index + 1 + i]['value'] = els.lines[name].th
params[2 * els.index + 1 + i] = els.lines[name].th
else:
parinfo[2 * els.index + 1 + i]['value'] = 1
params[2 * els.index + 1 + i] = 1
#Free parameters
......@@ -585,7 +587,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
if (not free_ratio) & (i == j) & (els.lines[name].fit is True):
#When line is reference line, the ratio is fixed equal to 1
parinfo[2 * els.index + 1 + i]['limited'] = [0, 0]
parinfo[2 * els.index + 1 + i]['value'] = params[2 * els.index + 1 + i] = 1
parinfo[2 * els.index + 1 + i]['value'] = 1
params[2 * els.index + 1 + i] = 1
parinfo[els.index + 1 + i]['fixed'] = 0
#XXX When line ratio is fixed: all ratios are fixed equal to theoretical values and intensities are tied (the ref intensity is free)
......@@ -609,7 +612,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#parinfo[2 * els.index + 1 + i]['limited'] = [0, 0]
#parinfo[2 * els.index + 1 + i]['value'] = params[2 * els.index + 1 + i] = 1
#parinfo[els.index + 1 + i]['fixed'] = 0
#ipdb.set_trace()
#Continuum parametres
for i in range(conf['DGCTNUM'] + 1):
parinfo[1 + els.index * 3 + i]['fixed'] = 0
......@@ -617,7 +621,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Output initialisation to be filled
paramscube = np.zeros(np.append(len(params), cube.shape[1:]))
perrorcube = np.zeros(np.append(len(params), cube.shape[1:]))
statusmap = np.zeros(cube.shape)
statusmap = np.zeros(cube.shape[1:])
dofmap = np.zeros(cube.shape[1:])
fnormmap = np.zeros(cube.shape[1:])
modcube = np.zeros(cube.shape)
......@@ -627,8 +631,12 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
fa = {'wave': wave, 'spectrum': None, 'err': None, 'lines': linewave, 'psf': conf['SPSF'], 'degc': conf['DGCTNUM']}
#Loop on pixels
counter = 0
number_of_pixels = cube.shape[2] * cube.shape[1]
for x in range(cube.shape[2]):
for y in range(cube.shape[1]):
if (counter / 10.) == np.int(counter / 10.): print('pixel %i / %i'%(counter, number_of_pixels))
counter += 1
fa['spectrum'] = cube[indok,y,x] # Spectrum
fa['err'] = np.sqrt(var[indok,y,x]) # Error
#fa['weight'] = np.sqrt(weight[indok,y,x]) # Weight
......@@ -1006,6 +1014,5 @@ def main(argv):
camel(options.filename, plot=options.plot, debug=options.debug, free_ratio=options.free_ratio, multi_ext=options.multi_ext, config=opts)
return
if __name__ == "__main__":
main(sys.argv[1:])
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment