Commit 33c16356 authored by bepinat's avatar bepinat

Add an option (ZRMIN, ZRMAX) in both create_config, camel and samel to...

Add an option (ZRMIN, ZRMAX) in both create_config, camel and samel to distinguish between the range for searching the redshift of the spectrum and the spectral range considered around lines, currently in the _test version
parent c08b2f43
......@@ -80,13 +80,13 @@ def compute_fwhm(dz, z=None, zcosmo=0):
def readconfig(filename, config):
"""This function reads the configuration file and returns a configuration dictionnary
"""This function reads the configuration file and returns a configuration dictionary
Parameters
----------
filename: string
name of the configuration file
config: dictionnary
config: dictionary
contains options with priority on the configuration file
"""
......@@ -135,7 +135,7 @@ def readconfig(filename, config):
else:
conf[key] = value
# add input options to configuration dictionnary + check needed keywords
# add input options to configuration dictionary + check needed keywords
if config is not None:
for key in config:
if config[key] is not None:
......@@ -166,7 +166,7 @@ def readcubes(conf, multi_ext=False):
Parameters
----------
conf: configuration dictionnary
conf: configuration dictionary
multi_ext: bool
keyword to indicate if input fits has multiple extensions
......@@ -302,7 +302,7 @@ def cutcube(cube, hdr, conf):
input data cube or variance array
hdr: fits header
input data cube header
conf: configuration dictionnary
conf: configuration dictionary
"""
......@@ -472,7 +472,7 @@ class line:
class lines:
"""This class enables to deal with lines. A dictionnary stored in lines will contain informations on each lines.
"""This class enables to deal with lines. A dictionary stored in lines will contain informations on each lines.
"""
......@@ -600,6 +600,7 @@ def build_parinfo_lines(els, linename, conf, free_ratio, comind, multi_comp=1):
zcosmo = conf['REDSHIFT']
zmax = conf['REDMAX']
zmin = conf['REDMIN']
#xxx ICI redefinir zmin et zmax avec le nouveau mot clé pour définir les limites et idem utiliser ce nouveau zmin, zmax pour explorer en z les parametres initiaux
parbase = {'value': 0., 'fixed': 1, 'limited': [0, 0], 'limits': [0., 0.], 'tied': '', 'mpmaxstep': 0, 'mpminstep': 0, 'parname': ''}
parinfo = []
......@@ -742,6 +743,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Condition on the wavelength to be within the redshift cuts
#XXX try zmin > zmax?
#xxx CES zmin, zmax sont bien ceux à utiliser pour les bornes en lbda
zcosmo = conf['REDSHIFT']
zmax = conf['REDMAX']
zmin = conf['REDMIN']
......@@ -824,6 +826,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
zmean = (zmax + zmin) / 2.
zzmax = zmean + np.floor((zmax - zmean) / zstep + 1) * zstep
zzmin = zmean + np.ceil((zmin - zmean) / zstep) * zstep
#xxx ici utiliser le nouveau zmin, zmax (nouveau mot clé), car cela définit les bornes de l'exploration en z
#Parameters contraints and informations: each line has z, sigma, intensity and ratio constraint
......@@ -1202,8 +1205,8 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
keyword to indicate the number of components to be adjusted
multi_ext: bool
keyword to indicate if input fits has multiple extensions
config: dictionnary
optional dictionnary that contains options usually stored in the configuration file
config: dictionary
optional dictionary that contains options usually stored in the configuration file
"""
# XXX gérer message d'erreur si on trouve pas le fichier d'input, s'il n'y a aucune raie, si un paramètre a un problème (e.g. de formatage)...
......@@ -1326,8 +1329,8 @@ def usage():
keyword to allow unconstrained line ratio
multi_ext: bool
keyword to indicate if input fits has multiple extensions
config: dictionnary
optional dictionnary that contains options usually stored in the configuration file
config: dictionary
optional dictionary that contains options usually stored in the configuration file
#######################
......@@ -1417,7 +1420,7 @@ def main(argv):
if np.size(argv) < 1:
sys.exit(2)
opts = vars(options) # this is a dictionnary
opts = vars(options) # this is a dictionary
opts['OII3726'] = opts['OII3729']
camel(options.filename, plot=options.plot, debug=options.debug, free_ratio=options.free_ratio, multi_comp=options.multi_comp, multi_ext=options.multi_ext, ftol=options.ftol, xtol=options.xtol, gtol=options.gtol, maxiter=options.maxiter, factor=options.factor, config=opts, cut=options.cut, clip=options.clip, smooth=options.smooth)
return
......
#!/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 time
import os
import sys
import copy
import numpy as np
#import pyfits as pf
import astropy.io.fits as fits
from cap_mpfit import mpfit # python 2.7 & 3!
#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 argparse
import logging
logging.basicConfig(level=logging.DEBUG)
logger = logging.getLogger('CAMEL')
def dz_from_dv(dv, z):
'''
This function enables to compute the redshift interval corresponding to a given velocity interval
Parameters
----------
dv: float
velocity interval in km/s
z: float
mean redshift of the source
Returns
-------
dz: float
corresponding redshift interval
'''
dz = dv * 1e3 / ct.c * (1 + z)
return dz
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 is 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 is 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 dictionary
Parameters
----------
filename: string
name of the configuration file
config: dictionary
contains options with priority on the configuration file
"""
conf = {}
if filename is not None:
try:
data = open(filename)
except:
logger.info('Camel is not able to open ' + filename)
data = open(filename)
count_extral = 0
for raw in data:
keyvalcom = raw.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.count('.') == 1):
try:
value = float(value)
except:
value = value
#elif value.isdigit():
#value = int(value)
elif (value.count('.') == 0):
try:
value = int(value)
except:
value = value
if (key == 'EXTRAL') & (value is not None):
if count_extral == 0:
conf[key] = list()
conf[key].append([value])
count_extral += 1
else:
conf[key] = value
# add input options to configuration dictionary + check needed keywords
if config is not None:
for key in config:
if config[key] is not None:
conf[key] = config[key]
elif not(key in conf.keys()):
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', 'NII6548', 'NII6583', 'SII6716', 'SII6731', 'OIII4959', 'OIII5007', 'HBETA', 'OII3729', 'OII3726', 'EXTRAL', 'COMMW', 'REDSHIFT', 'REDMIN', 'REDMAX', 'ZRMIN', 'ZRMAX', 'INITW', 'WMIN', 'WMAX', 'DFIT', 'DGCTNUM', 'MFIT', 'SCLIP', 'XYCLIP', 'NCLIP', 'SPSF', 'WSMOOTH', 'SSMOOTH', 'THRES', 'MEDIAN', 'FITSPSF', 'XMIN', 'YMIN', 'ZMIN', 'XMAX', 'YMAX', 'ZMAX']
for key in needed_keys:
if not(key in conf.keys()):
conf[key] = None
if 'OII' in conf.keys():
conf['OII3729'] = conf['OII3726'] = conf.pop('OII')
if conf['SPSF'] is None:
conf['SPSF'] = 0.
if conf['COMMW'] is None:
conf['COMMW'] = False
# XXX or give an error message but in that case, must test that the keywords have correct values
if conf['ZRMIN'] is None:
deltared = conf['REDMAX'] - conf['REDMIN']
conf['ZRMIN'] = conf['REDMIN'] - deltared / 2.
if conf['ZRMAX'] is None:
deltared = conf['REDMAX'] - conf['REDMIN']
conf['ZRMAX'] = conf['REDMAX'] + deltared / 2.
return conf
def readcubes(conf, multi_ext=False):
"""This function reads the input cubes (cube and variance)
Parameters
----------
conf: configuration dictionary
multi_ext: bool
keyword to indicate if input fits has multiple extensions
"""
if multi_ext:
hdul = fits.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 = None
varhdr = ''
else:
hdul = fits.open(conf['SKYFILE'])
var = hdul[0].data
varhdr = hdul[0].header
else:
hdul = fits.open(conf['FITSFILE'])
cube = hdul[0].data
hdr = hdul[0].header
if conf['SKYFILE'] == conf['FITSFILE']:
var = np.zeros(cube.shape, dtype='>f4')
if conf['XMIN'] is None:
conf['XMIN'] = 0
if conf['YMIN'] is None:
conf['YMIN'] = 0
if conf['XMAX'] is None:
conf['XMAX'] = cube.shape[2] - 1
if conf['YMAX'] is None:
conf['YMAX'] = cube.shape[1] - 1
for i in np.arange(cube.shape[0]):
var[i, :, :] = np.var(cube[i, conf['YMIN']:conf['YMAX'], conf['XMIN']:conf['XMAX']])
varhdr = hdul[0].header
var = writedata(var, varhdr, conf['OUTPUT'] + '_variance.fits')
elif (not conf['SKYFILE']):
var = None
varhdr = ''
else:
hdul = fits.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 is None:
return
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 = fits.PrimaryHDU(data, hdr)
hdulist = fits.HDUList(hdu)
hdulist.writeto(filename, overwrite=True, output_verify='fix')
return hdulist[0].header
def muse_whitelight_image(cube, hdr, filename):
"""This function creates and writes a white light image
Parameters
----------
data: ndarray
cube to be collapsed
hdr: fits header
fits header of the cube
filename: string
name of the output
"""
data = np.sum(cube, axis=0)
hdrw = copy.deepcopy(hdr)
if 'CTYPE3' in hdr.keys():
del hdrw['CTYPE3']
if 'CRVAL3' in hdr.keys():
del hdrw['CRVAL3']
if 'CRPIX3' in hdr.keys():
del hdrw['CRPIX3']
if 'CUNIT3' in hdr.keys():
del hdrw['CUNIT3']
if 'CD3_3' in hdr.keys():
del hdrw['CD3_3']
if 'CD3_2' in hdr.keys():
del hdrw['CD3_2']
if 'CD3_1' in hdr.keys():
del hdrw['CD3_1']
if 'CD2_3' in hdr.keys():
del hdrw['CD2_3']
if 'CD1_3' in hdr.keys():
del hdrw['CD1_3']
wlh = writedata(data, hdrw, filename)
return wlh
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 dictionary
"""
if conf['XMIN'] is None:
conf['XMIN'] = 0
if conf['YMIN'] is None:
conf['YMIN'] = 0
if conf['ZMIN'] is None:
conf['ZMIN'] = 0
if conf['XMAX'] is None:
conf['XMAX'] = cube.shape[2] - 1
if conf['YMAX'] is None:
conf['YMAX'] = cube.shape[1] - 1
if conf['ZMAX'] is None:
conf['ZMAX'] = cube.shape[0] - 1
if cube.ndim == 3:
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 == 1:
conf['ZMIN'] = np.max([conf['ZMIN'], 0], axis=0)
conf['ZMAX'] = np.min([conf['ZMAX'], np.array(cube.shape) - 1], axis=0)
conf['ZMIN'] = np.min([conf['ZMIN'], conf['ZMAX']], axis=0)
conf['ZMAX'] = np.max([conf['ZMAX'], conf['ZMIN']], axis=0)
if cube.ndim == 3:
cubecut = cube[int(conf['ZMIN']):int(conf['ZMAX']) + 1, int(conf['YMIN']):int(conf['YMAX']) + 1, int(conf['XMIN']):int(conf['XMAX']) + 1]
hdr['CRPIX1'] -= conf['XMIN']
hdr['CRPIX2'] -= conf['YMIN']
hdr['CRPIX3'] -= conf['ZMIN']
elif cube.ndim == 1:
cubecut = cube[int(conf['ZMIN']):int(conf['ZMAX']) + 1]
hdr['CRPIX1'] -= conf['ZMIN']
return cubecut, hdr
def remove_nanvar(var):
"""This function removes nan in variance cube if any and replaces them by median variance (spectrally)
Parameters
----------
var: array
input variance cube
Returns
-------
corrected variance
"""
cnan = np.isnan(var)
nanmap = cnan.prod(axis=0)
indnan = np.where(nanmap == 1)
indnum = np.where(nanmap == 0)
qvar = var[:, indnum[0], indnum[1]].reshape(var.shape[0], np.size(indnum[0]))
medianvar = np.median(qvar, axis=1)
if np.size(indnan[0]) == 0:
logger.info('Variance cube contains NaN and has been corrected')
for ind in range(np.size(indnan[0])):
var[:, indnan[0][ind], indnan[1][ind]] = medianvar
return var
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, dtype='>f4') # 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, dtype='>f4') # 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, dtype='>f4') # 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 is 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 dictionary 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.801, ref='HBETA', low=2.75, th=2.85))
self.append(line('HBETA', 4861.363))
self.append(line('HGAMMA', 4340.47, ref='HBETA', low=0.44, up=0.5, th=0.468))
self.append(line('HDELTA', 4101.73, ref='HBETA', low=0.23, up=0.29, th=0.259))
self.append(line('HEPS', 3970.07, ref='HBETA', low=0.13, up=0.19, th=0.159))
self.append(line('NII6583', 6583.45, ref='NII6548', low=2.7, up=3.3, th=3.))
#self.append(line('NII6583', 6583., ref='NII6548', low=2.7, up=3.3, th=3.))
self.append(line('NII6548', 6548.05))
self.append(line('SII6731', 6730.82))
self.append(line('SII6716', 6716.44, ref='SII6731', low=0.45, up=1.45, th=1.))
self.append(line('OIII5007', 5006.843, ref='OIII4959', low=2.7, up=3.3, th=3.))
#self.append(line('OIII5007', 5006.843, ref='OIII4959', low=2.9, up=3.1, th=3.))
self.append(line('OIII4959', 4958.911))
self.append(line('OIII4363', 4363.21, ref='OIII4959')) # XXX low=, up=, th=
self.append(line('OII3729', 3728.80, ref='OII3726', low=0.35, up=1.5, th=1.))
#self.append(line('OII3729', 3728.80, ref='OII3726', low=0.7, up=1.5, th=1.))
#self.append(line('OII3729', 3728.80, ref='OII3726', low=0.95, up=1.05, th=1.))
#self.append(line('OII3729', 3728.80, ref='OII3726', low=1.3, up=1.5, th=1.4))
self.append(line('OII3726', 3726.04))
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.3))
self.append(line('NEIII3868', 3868.))
#self['Ha'] =
#self.names = ['Ha', 'Hb', 'Hga', 'Hd', 'Heps', 'NII6583', 'NII6548', 'SII6731', 'SII6716', '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.,'SII6716':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','SII6716':'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', 'micron', 'micrometers', 'mum'}:
lconvfac = 1e4
if cunit.strip().lower() in {'angstroms', 'angstrom'}:
lconvfac = 1
if cunit.strip().lower() in {'nanometers', 'nm'}:
lconvfac = 1e1
if cunit.strip().lower() in {'centimeters', 'cm'}:
lconvfac = 1e8
if 'CD3_3' in hdr.keys():
sclp = hdr['CD3_3']
elif 'CDELT3' in hdr.keys():
sclp = hdr['CDELT3']
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, zind=None, sigind=None, intind=None, ratioind=None, contind=None):
wavel = wave.reshape(wave.size, 1)
deg = np.arange(degc + 1)
coefs = p[contind]
continuum = (coefs * (wavel - wave[0]) ** deg).sum(axis=1)
wavel1 = (lines * (p[zind] + 1))
dlines = (p[sigind] * lines)
coefs = (p[intind] * p[ratioind])
slines = (coefs * np.exp(-0.5 * (wavel - wavel1) ** 2 / (dlines ** 2 + psf ** 2))).sum(axis=1)
#wavel = wave.reshape(wave.size, 1)
#deg = np.arange(degc + 1)