Commit f1ef23fb authored by bepinat's avatar bepinat

Initial commit with contributors, code and licence

parents
This program has been developped from the CAMEL IDL code written in the frame of MASSIV and MUSE.
This programs is based on Levendberg-Marquart minimisation technique and uses the mpfit library that can be found here:
https://code.google.com/p/astrolibpy/downloads/detail?name=mpfit.tar.gz&can=2&q=
It runs with python 2.7 and uses standard python libraries:
time, os, sys, copy, numpy, pyfits, scipy, optparse
\ No newline at end of file
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
Created on 12/2014
@author: Benoît Epinat
"""
#import getopt
import time
import os, sys
import copy
import numpy as np
import pyfits as pf
from mpfit import mpfit # python 2.7...
#import ipdb
#from matplotlib import pyplot as plt
from scipy import ndimage
from scipy import constants as ct
import optparse
def compute_rv(z, zcosmo=None):
"""This function computes the velocity from local and global redshifts
Parameters
----------
z: float (or numpy array of floats)
local redshift (local velocity)
zcosmo: float
cosmological redshift (sytemic velocity)
"""
if zcosmo == None:
v = ((z + 1) ** 2 - 1) / ((z + 1) ** 2 + 1) * ct.c * 1e-3
else:
v = (z - zcosmo) / (1 + z) * ct.c * 1e-3
return v
def compute_fwhm(dz, z=None, zcosmo=0):
"""This function computes the velocity variation
Parameters
----------
dz: float (or numpy array of floats)
local redshift variation (local velocity dispersion)
z: float (or numpy array of floats)
local redshift (local velocity)
zcosmo: float
cosmological redshift (sytemic velocity)
"""
if z == None: z = zcosmo
dv = dz * (1 + zcosmo) / (1 + z) ** 2 * ct.c *1e-3
return dv
def readconfig(filename, config):
"""This function reads the configuration file and returns a configuration dictionnary
Parameters
----------
filename: string
name of the configuration file
config: dictionnary
contains options with priority on the configuration file
"""
conf = {}
if filename != None:
data = open(filename)
for line in data:
keyvalcom = line.split('= ')
key = keyvalcom[0].rstrip() # on vire les espaces de fin de keyword
value = keyvalcom[1].split('/')[0].rstrip('\t').rstrip() # on vire les tabultation et les espaces de fin de keyword
value = value.replace("'", '')
value = value.replace('"', '')
if value == '':
value = None
elif value.upper() == 'NONE':
value = None
elif value.upper() == 'FALSE':
value = False
elif value.upper() == 'TRUE':
value = True
elif (value.count('.') == 1) & (value.replace('.', '').isdigit()):
value = float(value)
elif value.isdigit():
value = int(value)
conf[key] = value
# add input options to configuration dictionnary + check needed keywords
if config != None:
keys = conf.keys()
for key in config:
if config[key] != None:
conf[key] = config[key]
#elif keys.isdisjoint([key]): # c'est du python3
elif (not conf.has_key(key)): # c'est du python 2.7
conf[key] = None
# XXX or give an error message but in that case, must test that the keywords have correct values
## XXX more keywords? add default values?
needed_keys = ['FITSFILE', 'OUTPUT', 'SKYFILE', 'HALPHA', 'NII6545', 'NII6584', 'SII6717', 'SII6731', 'OIII4959', 'OIII5007', 'HBETA', 'OII3729', 'OII3726', 'COMMW', 'REDSHIFT', 'REDMIN', 'REDMAX', 'INITW', 'WMIN', 'WMAX', 'DFIT', 'DGCTNUM', 'MFIT', 'SCLIP', 'XYCLIP', 'NCLIP', 'SPSF', 'WSMOOTH', 'SSMOOTH', 'THRES', 'MEDIAN', 'FITSPSF', 'XMIN', 'YMIN', 'ZMIN', 'XMAX', 'YMAX', 'ZMAX']
keys = conf.keys()
for key in needed_keys:
#if keys.isdisjoint([key]): # c'est du python3
if (not conf.has_key(key)): # c'est du python 2.7
conf[key] = None
if conf.has_key('OII'):
conf['OII3729'] = conf['OII3726'] = conf.pop('OII')
if conf['SPSF'] == None: conf['SPSF'] = 0.
# XXX or give an error message but in that case, must test that the keywords have correct values
return conf
def readcubes(conf, multi_ext=False):
"""This function reads the input cubes (cube and variance)
Parameters
----------
conf: configuration dictionnary
multi_ext: bool
keyword to indicate if input fits has multiple extensions
"""
if multi_ext:
hdul = pf.open(conf['FITSFILE'])
cube = hdul[1].data
hdr = hdul[1].header
if conf['SKYFILE'] == conf['FITSFILE']:
var = hdul[2].data
varhdr = hdul[2].header
elif (not conf['SKYFILE']):
var = 0
varhdr = ''
else:
hdul = pf.open(conf['SKYFILE'])
var = hdul[0].data
varhdr = hdul[0].header
else:
hdul=pf.open(conf['FITSFILE'])
cube=hdul[0].data
hdr=hdul[0].header
if (not conf['SKYFILE']):
var = 0
varhdr = ''
else:
hdul = pf.open(conf['SKYFILE'])
var = hdul[0].data
varhdr = hdul[0].header
return cube, hdr, var, varhdr
def testcubes(cube, var):
"""This function tests the sizes of the variance data
Parameters
----------
cube: ndarray
input data cube
var: ndarray
input variance data
"""
if var.ndim == 3:
if var.shape != cube.shape:
sys.exit(2)
if var.ndim == 1:
if var.shape[0] != cube.shape[0]:
sys.exit(2)
return
def writedata(data, hdr, filename):
"""This function writes the output data
Parameters
----------
data: ndarray
data to be written
hdr: fits header
fits header to be written
filename: string
name of the output
"""
hdu = pf.PrimaryHDU(data, hdr)
hdulist = pf.HDUList(hdu)
hdulist.writeto(filename, clobber=True)
return hdulist[0].header
def cutcube(cube, hdr, conf):
"""This function cuts the cube according to the limits requested
Parameters
----------
cube: array
input data cube or variance array
hdr: fits header
input data cube header
conf: configuration dictionnary
"""
if conf['XMIN'] == None:
conf['XMIN'] = 0
if conf['YMIN'] == None:
conf['YMIN'] = 0
if conf['ZMIN'] == None:
conf['ZMIN'] = 0
if conf['XMAX'] == None:
conf['XMAX'] = cube.shape[2] - 1
if conf['YMAX'] == None:
conf['YMAX'] = cube.shape[1] - 1
if conf['ZMAX'] == None:
conf['ZMAX'] = cube.shape[0] - 1
conf['ZMIN'], conf['YMIN'], conf['XMIN'] = np.max([[conf['ZMIN'], conf['YMIN'], conf['XMIN']], [0, 0, 0]],axis=0)
conf['ZMAX'], conf['YMAX'], conf['XMAX'] = np.min([[conf['ZMAX'], conf['YMAX'], conf['XMAX']], np.array(cube.shape) - 1],axis=0)
conf['ZMIN'], conf['YMIN'], conf['XMIN'] = np.min([[conf['ZMIN'], conf['YMIN'], conf['XMIN']], [conf['ZMAX'], conf['YMAX'], conf['XMAX']]],axis=0)
conf['ZMAX'], conf['YMAX'], conf['XMAX'] = np.max([[conf['ZMAX'], conf['YMAX'], conf['XMAX']], [conf['ZMIN'], conf['YMIN'], conf['XMIN']]],axis=0)
if cube.ndim == 3:
cubecut = cube[conf['ZMIN']:conf['ZMAX'] + 1, conf['YMIN']:conf['YMAX'] + 1, conf['XMIN']:conf['XMAX'] + 1]
hdr['CRPIX1'] -= conf['XMIN']
hdr['CRPIX2'] -= conf['YMIN']
hdr['CRPIX3'] -= conf['ZMIN']
elif cube.ndim == 1:
cubecut = cube[conf['ZMIN']:conf['ZMAX'] + 1]
hdr['CRPIX1'] -= conf['ZMIN']
return cubecut, hdr
def clipcube(cube, hdr, xy=3, sclip=3, npass=10):
"""This function performs a sigma clipping on the cube
Parameters
----------
cube: array
input datacube
hdr: fits header
input data cube header
xy: integer
size of the box in pixels in which the clipping is done
sclip: float
value of the clipping (n sigma)
npass: integer
number of pass to perform the clipping
"""
cube1 = np.zeros(cube.shape) # initialization of output cube
for z in range(cube.shape[0]): # loop on the z range
im = cube[z,:]
count = 0
ok = False
while count < npass | ok:
count += 1
im1 = ndimage.median_filter(im, xy)
diff = im - im1
bad = np.abs(diff) > (sclip * np.std(diff))
im[bad] = im1[bad]
if np.size(bad) == 0:
ok=True
cube1[z,:]=im
hdr.append(('XYCLIP', xy, 'Size of the clipping box'))
hdr.append(('SCLIP', sclip, 'n-sigma clipping'))
hdr.append(('NCLIP', npass, 'Number of pass of the clipping process'))
return cube1, hdr
def spatialsmooth(cube, hdr, fwhm):
"""This function performs a Gaussian spatial smoothing on the cube
Parameters
----------
cube: array
input datacube
hdr: fits header
input data cube header
fwhm: float
full width half maximum of the 2D Gaussian kernel
"""
cube1 = np.zeros(cube.shape) # initialization of output cube
sigma = fwhm / (2. * np.sqrt(2. * np.log(2.)))
for z in range(cube.shape[0]): # loop on the z range
cube1[z,:] = ndimage.gaussian_filter(cube[z,:], sigma)
hdr.append(('SSMOOTH', fwhm, 'FWHM in pixels of 2D Gaussian spatial smoothing'))
return cube1, hdr
def spectralsmooth(cube, hdr, fwhm):
"""This function performs a Gaussian spectral smoothing on the cube
Parameters
----------
cube: array
input datacube
hdr: fits header
input data cube header
fwhm: float
full width half maximum of the 1D Gaussian kernel
"""
cube1 = np.zeros(cube.shape) # initialization of output cube
sigma = fwhm / (2. * np.sqrt(2. * np.log(2.)))
for x in range(cube.shape[2]): # loop on the x range
for y in range(cube.shape[1]): # loop on the y range
cube1[:, y, x] = ndimage.gaussian_filter(cube[:, y, x], sigma)
hdr.append(('WSMOOTH', fwhm, 'FWHM in pixels of Gaussian spectral smoothing'))
# XXX add Hanning, Gauss, ...
return cube1, hdr
class line:
"""This class is the basic for defining a line. A line is defined by its name, its wavelength and the reference line to which it is attached if the ratio has to be constrained.
"""
def __init__(self, name, wave, ref=None, fit=False, low=0, up=None, th=None, index=None):
self.index = index
self.name = name
self.wave = wave
self.fit = fit
if ref == None: ref = name
self.ref = ref
if ref == name:
self.low = 1
self.up = 1
self.th = 1
else:
self.low = low
self.up = up
self.th = th
class lines:
"""This class enables to deal with lines. A dictionnary stored in lines will contain informations on each lines.
"""
def append(self, line):
self.lines[line.name] = line
self.lines[line.name].index = self.index
self.index += 1
def __init__(self):
self.index = 0
self.lines = {}
self.append(line('HALPHA', 6562.8, ref='HBETA', low=2.75, th=2.85))
self.append(line('HBETA', 4861.))
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('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.))
self.append(line('OIII4959', 4959.))
self.append(line('OIII4363', 4363., ref='OIII4959')) # XXX low=, up=, th=
self.append(line('OII3729', 3729., ref='OII3726', low=0.3, up=1.5, th=1.))
self.append(line('OII3726', 3726.))
self.append(line('HEI4471', 4471.))
self.append(line('HEI5876', 5876., ref='HEI4471', low=2.5, th=2.5))
self.append(line('HEII4686', 4686.))
self.append(line('OI6300', 6300.))
self.append(line('NEIII3868', 3868.))
#self['Ha'] =
#self.names = ['Ha', 'Hb', 'Hga', 'Hd', 'Heps', 'NII6583', 'NII6548', 'SII6731', 'SII6717', 'OIII5007', 'OIII4959', 'OIII4363', 'OII3729', 'OII3726', 'HeI4471', 'HeI5876', 'HeII4686', 'OI6300', 'NeIII3868']
#self.waves = [6562.8, 4861., 4340., 4101., 3968., 6583., 6548., 6731., 6717., 5007., 4959., 4363., 3729., 3726., 4471., 5876., 4686., 6300., 3868.]
#self.ref = ['Hb', 'Hb', 'Hb', 'Hb', 'Hb', 'NII6548', 'NII6548', 'SII6731', 'SII6731', 'OIII4959', 'OIII4959', 'OIII4959', 'OII3726', 'OII3726', 'HeI4471', 'HeI4471', 'HeII4686', 'OI6300', 'NeIII3868']
#lines={'Ha':6562.8,'Hb':4861.,'Hga':4340.,'Hd':4101.,'Heps':3968.,'NII6583':6583.,'NII6548':6548.,'SII6731':6731.,'SII6717':6717.,'OIII5007':5007.,'OIII4959':4959.,'OIII4363':4363.,'OII3729':3729.,'OII3726':3726.,'HeI4471':4471.,'HeI5876':5876.,'HeII4686':4686.,'OI6300':6300.,'NeIII3868':3868.}
##Initialisation of groups of lines
#sgroup={'Ha':'Hb','Hb':'Hb','Hga':'Hb','Hd':'Hb','Heps':'Hb','NII6583':'NII6548','NII6548':'NII6548','SII6731':'SII6731','SII6717':'SII6731','OIII5007':'OIII4959','OIII4959':'OIII4959','OIII4363':'OIII4959','OII3729':'OII3726','OII3726':'OII3726','HeI4471':'HeI4471','HeI5876':'HeI4471','HeII4686':'HeII4686','OI6300':'OI6300','NeIII3868':'NeIII3868'}
def waveindgen(hdr):
'''
This function returns the wavelength index in angstroms from the header information and the conversion factor from header data to angstroms.
'''
cunit = hdr['CUNIT3']
if cunit.strip().lower() in {'microns', 'micrometers', 'mum'}: lconvfac=1e4
if cunit.strip().lower() in {'angstroms', 'angstrom'}: lconvfac=1
if cunit.strip().lower() in {'nanometers', 'nm'}: lconvfac=1e3
if 'CDELT3' in hdr.ascard.keys():
sclp = hdr['CDELT3']
elif 'CD3_3' in hdr.ascard.keys():
sclp = hdr['CD3_3']
wave = lconvfac * (sclp * (np.arange(hdr['NAXIS3']) - hdr['CRPIX3'] + 1) + hdr['CRVAL3'])
return wave, lconvfac, sclp
def elspectrum(p, wave=None, lines=None, psf=None, degc=None):
wavel = wave.reshape(1, wave.size)
deg = np.arange(degc + 1).reshape(degc + 1, 1)
coefs = p[lines.size * 3 + 1:lines.size * 3 + 1 + degc + 1].reshape(degc + 1, 1)
continuum = (coefs * (wavel - wave[0]) ** deg).sum(axis=0)
wlines = lines.reshape(lines.size, 1)
dlines = (p[1:lines.size + 1] * lines).reshape(lines.size, 1)
coefs = (p[lines.size + 1:lines.size * 2 + 1] * p[lines.size *2 + 1:lines.size * 3 + 1]).reshape(lines.size, 1)
slines = (coefs * np.exp( -0.5 * (wavel - wlines * (p[0] + 1)) ** 2 / (dlines ** 2 + psf ** 2))).sum(axis=0)
return (continuum + slines)
def myelspectrum(p, fjac=None, wave=None, spectrum= None, err=None, lines=None, psf=None, degc=None):
model = elspectrum(p, wave=wave, lines=lines, psf=psf, degc=degc)
status = 0
return [status, (spectrum - model) / err]
def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=False):
"""
"""
#Initialisation of lines
els = lines()
for line in els.lines:
#check the line name is set in the configuration file
#if (not conf.keys().isdisjoint([line])): # c'est du python3
if conf.has_key(line): # python 2.7
els.lines[line].fit = conf[line]
if conf['EXTRAL'] != None:
for i in np.arange(np.shape(conf['EXTRAL'])[0]):
els.append(line('EXTRAL%i'%(i+1), conf['EXTRAL'][i], fit=True))
#Construction of wavelength index
wave0, lconvfac, sclp = waveindgen(hdr) # wavelength array in Angstrom and convertion factor
#Condition on the wavelength to be within the redshift cuts
#XXX try zmin > zmax?
zcosmo = conf['REDSHIFT']
zmax = conf['REDMAX']
zmin = conf['REDMIN']
zrange = zmax-zmin
argsorted = np.argsort(np.array([els.lines[i].index for i in els.lines]))
lineindex = np.array([els.lines[i].index for i in els.lines])[argsorted]
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)
firstind0 = np.zeros(linefit.size)
lastind0 = np.zeros(linefit.size)
for i in range(linefit.size):
if (not linefit[i]): continue
linf = np.max([(zmin - zrange / 2. + 1) * linewave[i] , np.min(wave0)])
lsup = np.min([(zmax + zrange / 2. + 1) * linewave[i] , np.max(wave0)])
condo = np.where((wave0 >= linf) & (wave0 <= lsup))
if np.size(condo) != 0:
firstind0[i] = condo[0][0]
lastind0[i] = condo[0][np.size(condo) - 1]
indok = np.append(indok, condo)
indok = np.unique(indok)
wave = wave0[indok]
for i in range(linefit.size):
if (not linefit[i]): continue
linf = np.max([(zmin - zrange / 2. + 1) * linewave[i] , np.min(wave0)])
lsup = np.min([(zmax + zrange / 2. + 1) * linewave[i] , np.max(wave0)])
condo = np.where((wave >= linf) & (wave <= lsup))
if np.size(condo) != 0:
firstind[i] = condo[0][0]
lastind[i] = condo[0][np.size(condo) - 1]
#firstind = np.zeros(np.size(cond))
#lastind = np.zeros(np.size(cond))
#firstind0 = np.zeros(np.size(cond))
#lastind0 = np.zeros(np.size(cond))
#for i in range(np.size(cond)):
#linf = np.max([(zmin - zrange / 2. + 1) * linewave[cond[0][i]] , np.min(wave0)])
#lsup = np.min([(zmax + zrange / 2. + 1) * linewave[cond[0][i]] , np.max(wave0)])
#condo = np.where((wave0 >= linf) & (wave0 <= lsup))
#if np.size(condo) != 0:
#firstind0[i] = condo[0][0]
#lastind0[i] = condo[0][np.size(condo) - 1]
#indok = np.append(indok, condo)
#indok = np.unique(indok)
#wave = wave0[indok]
#for i in range(np.size(cond)):
#linf = np.max([(zmin - zrange / 2. + 1) * linewave[cond[0][i]] , np.min(wave0)])
#lsup = np.min([(zmax + zrange / 2. + 1) * linewave[cond[0][i]] , np.max(wave0)])
#condo = np.where((wave >= linf) & (wave <= lsup))
#if np.size(condo) != 0:
#firstind[i] = condo[0][0]
#lastind[i] = condo[0][np.size(condo) - 1]
#Weights
if var == None: var = np.ones([indok.size, cube.shape[1], cube.shape[2]])
if var.ndim == 1: var = np.tile(var.reshape(var.size, 1, 1), (1, cube.shape[1], cube.shape[2]))
#if var.ndim == 3: weight = ((1. / var[indok,:,:]) / np.sum(1. / var[indok,:,:], axis=0)) * indok.size
if var.ndim == 3: weight = ((1. / var) / np.sum(1. / var, axis=0)) * var.shape[0]
#Redshift step
zstep = (conf['DFIT'] * 1e3 / ct.c + zcosmo) / (1 - conf['DFIT'] * 1e3 / ct.c) - zcosmo
if zstep >= (zmax - zmin): zstep = (zmax - zmin) * 0.99
#Parameters contraints and informations
params = np.zeros(1 + els.index * 3 + conf['DGCTNUM'] + 1)
parbase = {'value':0., 'fixed':1, 'limited':[0,0], 'limits':[0.,0.], 'tied':'', 'mpmaxstep':0, 'mpminstep':0, 'parname':''}
parinfo = []
for i in range(len(params)):
parinfo.append(copy.deepcopy(parbase))
#Redshift range
parinfo[0]['limited'] = [1, 1]
parinfo[0]['limits'] = [zmin, zmax]
parinfo[0]['fixed'] = 0
#Line parameters
#When common dispersion, we determine the reference line for the width and set the dispersion as free parameter
commwindex = -1
if conf['COMMW']:
commwindex = np.min(lineindex[np.where(linefit)])
parinfo[1 + commwindex]['fixed'] = 0
linew = np.zeros(els.index)
for i in range(els.index):
name = linename[i]
j = els.lines[els.lines[name].ref].index
#Dispersion
parinfo[1 + i]['limited'] = [1, 1]
parinfo[1 + i]['limits'] = [(conf['WMIN'] * 1e3 / ct.c + zcosmo) / (1 - conf['WMIN'] * 1e3 / ct.c) - zcosmo, (conf['WMAX'] * 1e3 / ct.c + zcosmo) / (1 - conf['WMAX'] * 1e3 / ct.c) - zcosmo] #in redshift unit
parinfo[1 + i]['value'] = params[1 + i] = (conf['INITW'] * 1e3 / ct.c + zcosmo) / (1 - conf['INITW'] * 1e3 / ct.c) - zcosmo
#When dispersion is common for all lines, widths are tied to the first free line (the first free width is free)
if (conf['COMMW']) & (i != commwindex):
parinfo[1 + i]['tied'] = 'p[%i]'%(1 + commwindex)
#When dispersion is independant for all lines (common by groups), widths are tied to the reference line and the ref width is free
if (not conf['COMMW']) & (i != j):
parinfo[1 + i]['tied'] = 'p[%i]'%(1 + j)
if els.lines[name].fit:
parinfo[1 + j]['fixed'] = 0 # the width of the reference line is free
if (not conf['COMMW']) & (i == j) & (els.lines[name].fit is True): parinfo[1 + i]['fixed'] = 0 # the width of the reference line is free
#XXX When dispersion is completely independant, widths are not tied and are free
#if (not conf['COMMW']):
#if els.lines[name].fit: parinfo[1 + i]['fixed'] = 0 # the width of the reference line is free
#Intensity
parinfo[els.index + 1 + i]['limited'][0] = 1
parinfo[els.index + 1 + i]['limits'][0] = 0.
#When line ratio is free: all ratios are fixed (= 1) and lines intensities are not tied and free parameters
if free_ratio:
parinfo[2 * els.index + 1 + i]['value'] = 1
if els.lines[name].fit is True:
parinfo[els.index + 1 + i]['fixed'] = 0
#When line ratio is constrained: line intensities are tied and ratios are free, except for the reference line (the ref intensity is free, the ref ratio is fixed = 1)
if (not free_ratio) & (i != j):
parinfo[els.index + 1 + i]['tied'] = 'p[%i]'%(els.index + 1 + j) # line intensity is tied to reference
#Ratio limits
parinfo[2 * els.index + 1 + i]['limited'][0] = 1
if els.lines[name].low != None:
parinfo[2 * els.index + 1 + i]['limits'][0] = els.lines[name].low
if els.lines[name].up != None:
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
else:
params[2 * els.index + 1 + i] = 1
#Free parameters
if els.lines[name].fit is True:
parinfo[els.index + 1 + j]['fixed'] = 0 # the intensity of the reference line is free
parinfo[2 * els.index + 1 + i]['fixed'] = 0 # line ratio is free
#XXX When reference line is not fitted, we have to find a solution if several lines of a group are fitted
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[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)
#if (not free_ratio) & (i!=j):
#parinfo[els.index + 1 + i]['tied'] = 'p[%i]'%(els.index + 1 + j)
#parinfo[2 * els.index + 1 + i]['limited'][0] = 1
#if els.lines[name].low != None:
#parinfo[2 * els.index + 1 + i]['limits'][0] = els.lines[name].low
#if els.lines[name].up != None:
#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
#else:
#params[2 * els.index + 1 + i] = 1
#if els.lines[name].fit is True:
#parinfo[els.index + 1 + j]['fixed'] = 0 # the intensity of the reference line is free
#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[els.index + 1 + i]['fixed'] = 0
#Continuum parametres
for i in range(conf['DGCTNUM'] + 1):
parinfo[1 + els.index * 3 + i]['fixed'] = 0
#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)
dofmap = np.zeros(cube.shape[1:])
fnormmap = np.zeros(cube.shape[1:])
modcube = np.zeros(cube.shape)
#List of extra parameters
#xxx Ajuster wave et spectrum pour ne garder que ce qui nous intéresse?
fa = {'wave': wave, 'spectrum': None, 'err': None, 'lines': linewave, 'psf': conf['SPSF'], 'degc': conf['DGCTNUM']}
#Loop on pixels
for x in range(cube.shape[2]):
for y in range(cube.shape[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
if (np.min(fa['spectrum']) == np.max(fa['spectrum'])) & (np.min(fa['spectrum']) == 0): continue # No need for fitting
minfnorm = np.infty
#index = None
std = np.std(fa['spectrum'])
#Parameters initialisation
p = np.copy(params)
#Line intensity initisalisation
for i in range(els.index):
p[els.index + 1 + i] = 2 * std
#Inititalisation of continuum (degree 0)
p[1 + els.index * 3] = np.median(fa['spectrum'])
for z in np.arange(zmin, zmax, zstep):
#print(z)
p[0] = z
pstart = np.copy(p)
fit0 = mpfit.mpfit(myelspectrum, xall=pstart, functkw=fa, quiet=True, parinfo=parinfo)
if fit0.status == 0: continue
if fit0.status == -16: continue
if fit0.fnorm <= minfnorm:
minfnorm = fit0.fnorm
fit = copy.copy(fit0)
if minfnorm == np.infty: continue
#print(fit.params[2 * els.index + 1 + 12], parinfo[2 * els.index + 1 + 12]['limits'], fit.params[2 * els.index + 1 + 13])
paramscube[:, y, x] = fit.params[:]
perrorcube[:, y, x] = fit.perror[:]
fnormmap[y, x] = fit.fnorm
dofmap[y, x] = fit.dof
statusmap[y, x] = fit.status
modcube[indok, y, x] = elspectrum(fit.params, wave=wave, lines=linewave, psf=conf['SPSF'], degc=conf['DGCTNUM'])
#Creating outputs
#Error normalisation by reduced chi2
perrorcube = perrorcube * np.sqrt(fnormmap / dofmap)
#Continuum
#if we want the coefficient, we must multiply them by factor + give lambda[0]
contcube = np.zeros(cube.shape)
econtcube = np.zeros(cube.shape)
for i in range(conf['DGCTNUM'] + 1):
contcube += paramscube[1 + els.index * 3 + i,:,:] * (wave0 - wave[0]).reshape(wave0.size, 1, 1) ** i
econtcube += perrorcube[1 + els.index * 3 + i,:,:] * (wave0 - wave[0]).reshape(wave0.size, 1, 1) ** i
#contmap
#; continuum is computed at the longest wavelength, not necessarily optimal
#cont[m,n]+=pfinal[2*n_elements(l)+1+i]*double(lambda[n_elements(lambda)-1]-lambda[0])^i/double(lambda[n_elements(lambda)-1]-lambda[0])
#econt[m,n]+=(paramerr[2*n_elements(l)+1+i]*double(lambda[n_elements(lambda)-1]-lambda[0])^i/double(lambda[n_elements(lambda)-1]-lambda[0]))^2
#Residual cube
#residual