camel.py 49.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
#!/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, sys
import copy
import numpy as np
import pyfits as pf
from mpfit import mpfit  # python 2.7...
bepinat's avatar
bepinat committed
28
import ipdb
29 30 31 32
#from matplotlib import pyplot as plt
from scipy import ndimage
from scipy import constants as ct

bepinat's avatar
bepinat committed
33
import argparse
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85

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)
bepinat's avatar
bepinat committed
86 87 88
        count_extral = 0
        for raw in data:
            keyvalcom = raw.split('= ')
89
            key = keyvalcom[0].rstrip()   # on vire les espaces de fin de keyword
bepinat's avatar
bepinat committed
90
            value = keyvalcom[1].split('/ ')[0].rstrip('\t').rstrip()  # on vire les tabultation et les espaces de fin de keyword
91 92 93 94 95 96 97 98 99 100
            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
bepinat's avatar
bepinat committed
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
            #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 != None):
                if count_extral == 0:
                    conf[key] = np.array([value])
                else:
                    conf[key] = np.append(conf[key], value)
                count_extral += 1
            else:
                conf[key] = value
123 124 125 126 127 128 129 130 131 132 133 134 135

    # 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?
bepinat's avatar
bepinat committed
136
    needed_keys = ['FITSFILE', 'OUTPUT', 'SKYFILE', 'HALPHA', 'NII6548', 'NII6583', 'SII6716', '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']
137 138 139 140 141 142
    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'):
bepinat's avatar
bepinat committed
143
        conf['OII3729'] = conf['OII3726'] = conf.pop('OII')
144
    if conf['SPSF'] == None: conf['SPSF'] = 0.
bepinat's avatar
bepinat committed
145
    if conf['COMMW'] == None: conf['COMMW'] = False
146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
            # 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']):
bepinat's avatar
bepinat committed
168
            var = None
169 170 171 172 173 174
            varhdr = ''
        else:
            hdul = pf.open(conf['SKYFILE'])
            var = hdul[0].data
            varhdr = hdul[0].header
    else:
bepinat's avatar
bepinat committed
175 176 177
        hdul = pf.open(conf['FITSFILE'])
        cube = hdul[0].data
        hdr = hdul[0].header
178
        if (not conf['SKYFILE']):
bepinat's avatar
bepinat committed
179
            var = None
180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
            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
    
    """
bepinat's avatar
bepinat committed
198 199
    if var == None:
        return
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
    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

bepinat's avatar
bepinat committed
266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291
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:
        print('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

292 293 294 295 296 297 298 299 300 301 302 303 304 305 306
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
bepinat's avatar
bepinat committed
307
    
308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
    """

    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 = {}
bepinat's avatar
bepinat committed
404 405
        self.append(line('HALPHA', 6562.801, ref='HBETA', low=2.75, th=2.85))
        self.append(line('HBETA', 4861.363))
bepinat's avatar
bepinat committed
406
        self.append(line('HGAMMA', 4340.47, ref='HBETA', low=0.44, up=0.5, th=0.468))
bepinat's avatar
bepinat committed
407 408
        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))
bepinat's avatar
bepinat committed
409
        self.append(line('NII6583', 6583.45, ref='NII6548', low=2.7, up=3.3, th=3.))
bepinat's avatar
bepinat committed
410 411 412 413
        #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.))
bepinat's avatar
bepinat committed
414 415
        self.append(line('OIII5007', 5006.843, ref='OIII4959', low=2.7, up=3.3, th=3.))
        self.append(line('OIII4959', 4958.911))
416
        self.append(line('OIII4363', 4363., ref='OIII4959'))  # XXX low=, up=, th=
bepinat's avatar
bepinat committed
417 418
        self.append(line('OII3729', 3728.80, ref='OII3726', low=0.3, up=1.5, th=1.))
        self.append(line('OII3726', 3726.04))
419 420 421 422 423 424 425
        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'] = 
bepinat's avatar
bepinat committed
426
        #self.names = ['Ha', 'Hb', 'Hga', 'Hd', 'Heps', 'NII6583', 'NII6548', 'SII6731', 'SII6716', 'OIII5007', 'OIII4959', 'OIII4363', 'OII3729', 'OII3726', 'HeI4471', 'HeI5876', 'HeII4686', 'OI6300', 'NeIII3868']
427 428 429
        #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']
        
bepinat's avatar
bepinat committed
430
        #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.}
431
        ##Initialisation of groups of lines
bepinat's avatar
bepinat committed
432
        #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'}
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473

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()
bepinat's avatar
bepinat committed
474
    for lline in els.lines:
475 476
        #check the line name is set in the configuration file
        #if (not conf.keys().isdisjoint([line])): # c'est du python3
bepinat's avatar
bepinat committed
477 478
        if conf.has_key(lline): # python 2.7
            els.lines[lline].fit = conf[lline]
479 480 481 482 483 484 485 486 487 488 489 490 491

    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']
bepinat's avatar
bepinat committed
492
    zrange = zmax - zmin
493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547
    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]
    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
bepinat's avatar
bepinat committed
548 549
    #if var == None: var = np.ones([indok.size, cube.shape[1], cube.shape[2]])
    if var == None: var = np.ones(cube.shape)
550 551 552 553 554 555
    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
bepinat's avatar
bepinat committed
556 557 558 559 560

    if zstep >= (zmax - zmin) / 2: zstep = (zmax - zmin) / 2.
    zmean = (zmax + zmin) / 2.
    zzmax = zmean + (np.round(zmean / zstep) + 1) * zstep
    zzmin = zmean - np.round(zmean / zstep) * zstep
561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580
    
    #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)
bepinat's avatar
bepinat committed
581

582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627
    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:
bepinat's avatar
bepinat committed
628 629
                parinfo[2 * els.index + 1 + i]['value'] = els.lines[name].th
                params[2 * els.index + 1 + i] = els.lines[name].th
630
            else:
bepinat's avatar
bepinat committed
631
                parinfo[2 * els.index + 1 + i]['value'] = 1
632 633 634 635 636 637 638 639 640 641 642
                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]
bepinat's avatar
bepinat committed
643 644
            parinfo[2 * els.index + 1 + i]['value'] = 1
            params[2 * els.index + 1 + i] = 1
645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667
            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
bepinat's avatar
bepinat committed
668 669
    #ipdb.set_trace()

670 671 672 673 674 675 676
    #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:]))
bepinat's avatar
bepinat committed
677
    statusmap = np.zeros(cube.shape[1:])
678 679 680 681 682 683 684 685 686
    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
bepinat's avatar
bepinat committed
687
    counter = 0
bepinat's avatar
bepinat committed
688 689
    fac_prev = -1
    number_of_pixels = np.float(cube.shape[2] * cube.shape[1])
690 691
    for x in range(cube.shape[2]):
        for y in range(cube.shape[1]):
bepinat's avatar
bepinat committed
692 693 694 695 696 697
            if np.int((counter / number_of_pixels) * 100) != fac_prev:
                fac_prev = np.int((counter / number_of_pixels) * 100)
                sys.stdout.write("Progress => {:3d}%\r".format(fac_prev))
                sys.stdout.flush()
            #if (counter / 10.) == np.int(counter / 10.): print('pixel %i / %i'%(counter, number_of_pixels))
            
bepinat's avatar
bepinat committed
698
            counter += 1
699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717
            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'])
bepinat's avatar
bepinat committed
718
            for z in np.arange(zzmin, zzmax, zstep):
719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792
                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
    residualcube = np.zeros(cube.shape)
    residualcube[indok,:,:] = cube[indok,:,:] - modcube[indok,:,:]
    #dispcont
    dispcontmap = np.std(cube[indok,:,:] - modcube[indok,:,:], axis=0)
    
    #Chi2
    chi2map = fnormmap / dofmap
    
    #Velocity field
    zmap = paramscube[0,:,:]
    ezmap = perrorcube[0,:,:]
    rvmap = compute_rv(zmap, zcosmo=zcosmo)
    ervmap = compute_fwhm(ezmap, z=zmap, zcosmo=zcosmo)
    
    #Lines
    wavemaps = np.zeros(np.append(els.index, cube.shape[1:]))
    ewavemaps = np.zeros(np.append(els.index, cube.shape[1:]))
    wavedispmaps = np.zeros(np.append(els.index, cube.shape[1:]))
    ewavedispmaps = np.zeros(np.append(els.index, cube.shape[1:]))
    dispmaps = np.zeros(np.append(els.index, cube.shape[1:]))
    edispmaps = np.zeros(np.append(els.index, cube.shape[1:]))
    intmaps = np.zeros(np.append(els.index, cube.shape[1:]))
    eintmaps = np.zeros(np.append(els.index, cube.shape[1:]))
    fluxmaps = np.zeros(np.append(els.index, cube.shape[1:]))
    efluxmaps = np.zeros(np.append(els.index, cube.shape[1:]))    
    snrmaps = np.zeros(np.append(els.index, cube.shape[1:]))
    
    # dzmap and edzmap are intermediate products
    if conf['COMMW']:
        dzmap = paramscube[1 + commwindex,:,:]
        edzmap = perrorcube[1 + commwindex,:,:]
    
    for i in range(els.index):
        if (not linefit[i]): continue
        name = linename[i]
        j = els.lines[els.lines[name].ref].index
        
bepinat's avatar
bepinat committed
793 794 795 796
        #if (not conf['COMMW']) & (i != j):
            #dzmap = paramscube[1 + j,:,:]
            #edzmap = perrorcube[1 + j,:,:]
        if (not conf['COMMW']):
797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908
            dzmap = paramscube[1 + j,:,:]
            edzmap = perrorcube[1 + j,:,:]
        
        wavemaps[i,:,:] = (zmap + 1) * linewave[i]
        ewavemaps[i,:,:] = ezmap * linewave[i]
        wavedispmaps[i,:,:] = dzmap * linewave[i]
        ewavedispmaps[i,:,:] = edzmap * linewave[i]
        dispmaps[i,:,:] = compute_fwhm(dzmap, z=zmap, zcosmo=zcosmo)
        edispmaps[i,:,:] = compute_fwhm(edzmap, z=zmap, zcosmo=zcosmo)
        
        intmaps[i,:,:] = paramscube[1 + els.index + i,:,:] * paramscube[1 + 2 * els.index + i,:,:]
        fluxmaps[i,:,:] = np.sum(intmaps[i,:,:] * np.exp(-0.5 * (wave.reshape(wave.size, 1, 1) - linewave[i] * (1 + zmap)) ** 2 / ((dzmap * linewave[i]) ** 2 + conf['SPSF'] ** 2)) * sclp * lconvfac, axis=0)
        snrmaps[i,:,:] = np.sum(weight[firstind0[i]:lastind0[i],:,:] * (modcube[firstind0[i]:lastind0[i],:,:] - contcube[firstind0[i]:lastind0[i],:,:]), axis=0) * sclp * lconvfac / (np.sqrt(np.pi * 2) * np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2) * np.sqrt(np.sum(weight[firstind0[i]:lastind0[i],:,:] * (cube[firstind0[i]:lastind0[i],:,:] - modcube[firstind0[i]:lastind0[i],:,:]) ** 2, axis=0) / np.sum(weight[firstind0[i]:lastind0[i],:,:], axis=0))) * (lastind0[i] - firstind0[i]) / np.sum(weight[firstind0[i]:lastind0[i],:,:], axis=0)
        if free_ratio:
            eintmaps[i,:,:] = perrorcube[1 + els.index + i,:,:]
        elif i != j:
            eintmaps[i,:,:] = np.sqrt((paramscube[1 + els.index + i,:,:] * perrorcube[1 + 2 * els.index + i,:,:]) ** 2 + (perrorcube[1 + els.index + j,:,:] * paramscube[1 + 2 * els.index + i,:,:]) ** 2)
        else:
            eintmaps[i,:,:] = perrorcube[1 + els.index + i,:,:]
        
        efluxmaps[i,:,:] = np.sqrt(np.pi * 2) * np.sqrt((eintmaps[i,:,:] * np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2)) ** 2 + (linewave[i] * dzmap / np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2) * intmaps[i,:,:] * edzmap) ** 2)
    
    snrmap = np.sum(weight[indok,:,:] * (modcube[indok,:,:] - contcube[indok,:,:]), axis = 0) * sclp * lconvfac / (np.sqrt(np.pi * 2) * np.sqrt((dzmap * np.median(wave) / (1 + zmap)) ** 2 + conf['SPSF'] ** 2) * np.sqrt( np.sum(weight[indok,:,:] * (cube[indok,:,:] - modcube[indok,:,:]) ** 2, axis=0) / np.sum(weight[indok,:,:], axis=0)))
    
    #Output writing
    #Cubes
    hh = writedata(residualcube, hdr, conf['OUTPUT'] + '_residualcube.fits')
    hh = writedata(modcube, hdr, conf['OUTPUT'] + '_modcube.fits')
    hh = writedata(contcube, hdr, conf['OUTPUT'] + '_contcube.fits')
    hh = writedata(econtcube, hdr, conf['OUTPUT'] + '_econtcube.fits')
    
    #Images
    if 'CTYPE3' in hdr.ascard.keys(): del hdr['CTYPE3']
    if 'CRVAL3' in hdr.ascard.keys(): del hdr['CRVAL3']
    if 'CRPIX3' in hdr.ascard.keys(): del hdr['CRPIX3']
    if 'CUNIT3' in hdr.ascard.keys(): del hdr['CUNIT3']
    if 'CD3_3' in hdr.ascard.keys(): del hdr['CD3_3']
    if 'CD3_2' in hdr.ascard.keys(): del hdr['CD3_2']
    if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1']
    if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_3']
    if 'CD1_3' in hdr.ascard.keys(): del hdr['CD1_3']
    
    suffixe = ''
    if conf['COMMW']: 
        suffixe = '_common'
    else:
        suffixe = '_indep'
  
    #paramscube
    #perrorcube
    hh = writedata(dofmap, hdr, conf['OUTPUT'] + '_dof' + suffixe + '.fits')
    hh = writedata(fnormmap, hdr, conf['OUTPUT'] + '_fnorm' + suffixe + '.fits')
    hh = writedata(statusmap, hdr, conf['OUTPUT'] + '_status' + suffixe + '.fits')
    hh = writedata(chi2map, hdr, conf['OUTPUT'] + '_chi2' + suffixe + '.fits')
    hh = writedata(dispcontmap, hdr, conf['OUTPUT'] + '_distcont' + suffixe + '.fits')
    hh = writedata(zmap, hdr, conf['OUTPUT'] + '_z' + suffixe + '.fits')
    hh = writedata(ezmap, hdr, conf['OUTPUT'] + '_ez' + suffixe + '.fits')
    hh = writedata(rvmap, hdr, conf['OUTPUT'] + '_vel' + suffixe + '.fits')
    hh = writedata(ervmap, hdr, conf['OUTPUT'] + '_evel' + suffixe + '.fits')
    hh = writedata(snrmap, hdr, conf['OUTPUT'] + '_snr' + suffixe + '.fits')
    
    for i in range(els.index):
        if (not linefit[i]): continue
        hh = writedata(wavemaps[i,:,:], hdr, conf['OUTPUT'] + '_wave' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(ewavemaps[i,:,:], hdr, conf['OUTPUT'] + '_ewave' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(wavedispmaps[i,:,:], hdr, conf['OUTPUT'] + '_wavedisp' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(ewavedispmaps[i,:,:], hdr, conf['OUTPUT'] + '_ewavedisp' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(dispmaps[i,:,:], hdr, conf['OUTPUT'] + '_disp' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(edispmaps[i,:,:], hdr, conf['OUTPUT'] + '_edisp' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(intmaps[i,:,:], hdr, conf['OUTPUT'] + '_int' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(eintmaps[i,:,:], hdr, conf['OUTPUT'] + '_eint' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(fluxmaps[i,:,:], hdr, conf['OUTPUT'] + '_flux' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(efluxmaps[i,:,:], hdr, conf['OUTPUT'] + '_eflux' + suffixe + '_' + linename[i] + '.fits')
        hh = writedata(snrmaps[i,:,:], hdr, conf['OUTPUT'] + '_snr' + suffixe + '_' + linename[i] + '.fits')
    

    return rvmap

def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, config=None):
    """CAMEL (Cube Analysis: Moment maps of Emission Lines) computes moment maps

    Parameters
    ----------
    filename: string
        name of the configuration file
    plot: bool
        keyword to plot fits during the process (not implemented yet)
    debug: bool
        keyword to show debug information (not implemented yet)
    free_ratio: bool
        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

    """
    # 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)...
    # XXX gérer cas avec une LSF donnée en input (non bruitée)
    
    
    t0 = time.time()
    
    #Reading of configuration file
    print("Reading configuration file " + filename)
    try:
        conf = readconfig(filename, config)
    except:
        print("Configuration file problem")
        return
    
    #Reading data
bepinat's avatar
bepinat committed
909 910 911 912
    try:
        print("Reading data cube " + conf['FITSFILE'] + " and variance " + conf['SKYFILE'])
    except:
        print("Reading data cube " + conf['FITSFILE'] + ". No variance data.")
913 914 915 916 917 918 919 920 921
    try:
        cube, hdr, var, varhdr = readcubes(conf, multi_ext=multi_ext)
    except:
        print("Input cube problem")
        return
    
    try:
        testcubes(cube, var)
    except:
bepinat's avatar
bepinat committed
922
        print("Variance shape: %i; Data cube shape: %i"%(np.shape(var), cube.shape))
923 924 925 926 927 928 929 930 931 932
        print("Variance size problem")
        return
    
    #Cutting the cubes
    print("Cutting data cube")
    cube, hdr = cutcube(cube, hdr, conf)
    hdr = writedata(cube, hdr, conf['OUTPUT']+'_cube_cut.fits')
    if var != None:
        print("Cutting variance")
        var, varhdr = cutcube(var, varhdr, conf)
bepinat's avatar
bepinat committed
933 934
        if var.ndim == 3:
            var = remove_nanvar(var)
935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014
        varhdr = writedata(var, varhdr, conf['OUTPUT']+'_var_cut.fits')
    
    #Clipping the cube
    if (conf['SCLIP'] != 0) & (conf['SCLIP'] != None):
        print('Performing %3.1f-sigma clipping using box of size %i. Number of pass: %i'%(conf['SCLIP'], conf['XYCLIP'], conf['NCLIP']))
        cube, hdr = clipcube(cube, hdr, xy=conf['XYCLIP'], sclip=conf['SCLIP'], npass=conf['NCLIP'])
        hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube_cut_clip.fits')
    
    #Smoothing the cubes
    if (conf['WSMOOTH'] != 0) & (conf['WSMOOTH'] != None):
        print('Performing %3.1f pixels 1D spectral smoothing'%conf['WSMOOTH'])
        cube, hdr = spectralsmooth(cube, hdr, conf['WSMOOTH'])
        conf['OUTPUT'] = conf['OUTPUT'] + '_wsmooth'
        hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube.fits')
        if var != None:
            if var.ndim == 3:
                var, varhdr = spectralsmooth(var, varhdr, conf['WSMOOTH'])
                varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits') 
        
    if (conf['SSMOOTH'] != 0) & (conf['SSMOOTH'] != None):
        print('Performing %3.1f pixels 2D gaussian spatial smoothing'%conf['SSMOOTH'] )
        cube, hdr = spatialsmooth(cube, hdr, conf['SSMOOTH'])
        conf['OUTPUT'] = conf['OUTPUT'] + '_ssmooth'
        hdr = writedata(cube, hdr, conf['OUTPUT'] + '_cube.fits')
        if var != None:
            if var.ndim == 3:
                var, varhdr = spatialsmooth(var, varhdr, conf['SSMOOTH'])
                varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')                
        
    #Creating the maps
    rv = buildmaps(conf, cube, hdr, var=var, plot=plot, debug=debug, free_ratio=free_ratio)

    t1 = time.time()
    #print(t1-t0)
    
    return

#def usage(option, opt, value, parser):
def usage():
    """CAMEL (Cube Analysis: Moment maps of Emission Lines) can by launched from command line
    or from a python console.
    
    #########################
    #Python console solution#
    #########################
    
    Syntaxe:
    --------
    import camel
    camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, config=None)
    
    Parameters
    ----------
    filename: string
        name of the configuration file
    plot: bool
        keyword to plot fits during the process (not implemented yet)
    debug: bool
        keyword to show debug information (not implemented yet)
    free_ratio: bool
        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

    
    #######################
    #Command line solution#
    #######################
    
    Syntaxe:
    --------
    
    """
    print(usage.__doc__)

def main(argv):
    """
    """
bepinat's avatar
bepinat committed
1015 1016 1017 1018 1019 1020 1021 1022 1023
    parser = argparse.ArgumentParser(description=usage())
    parser.add_argument('--file', '-f', action="store", dest="filename", default=None, help="name of the configuration file")
    parser.add_argument('--plot', '-p', action="store_true", dest="plot", default=False, help="keyword to plot fits during the process (not implemented yet)")
    parser.add_argument('--debug', '-d', action="store_true", dest="debug", default=False, help="keyword to show debug information (not implemented yet)")
    parser.add_argument('--free_line_ratio', '-r', action="store_true", dest="free_ratio", default=False, help="keyword to allow unconstrained line ratio")
    parser.add_argument('--multi_ext', '-m', action="store_true", dest="multi_ext", default=False, help="keyword to indicate if input fits has multiple extensions")
    parser.add_argument('--FITSFILE', action="store", dest="FITSFILE", default=None, help="name of the input cube")
    parser.add_argument('--OUTPUT', action="store", dest="OUTPUT", default=None, help="generic output name")
    parser.add_argument('--SKYFILE', action="store", dest="SKYFILE", default=None, help="name of the variance cube (can be a sky cube)")
bepinat's avatar
bepinat committed
1024 1025 1026
    parser.add_argument('--HALPHA', action="store_true", dest="HALPHA", default=None, help="keyword to fit Halpha@6562.801 line")
    parser.add_argument('--HBETA', action="store_true", dest="HBETA", default=None, help="keyword to fit Hbeta@4861.363 line")
    parser.add_argument('--HGAMMA', action="store_true", dest="HGAMMA", default=None, help="keyword to fit Hgamma@4340.47 line")
bepinat's avatar
bepinat committed
1027 1028
    parser.add_argument('--HDELTA', action="store_true", dest="HDELTA", default=None, help="keyword to fit Hdelta@4101.73 line")
    parser.add_argument('--HEPS', action="store_true", dest="HEPS", default=None, help="keyword to fit Hepsilon@3970.07 line")
bepinat's avatar
bepinat committed
1029 1030 1031 1032
    parser.add_argument('--NII6548', action="store_true", dest="NII6548", default=None, help="keyword to fit NII@6548.05 line")
    parser.add_argument('--NII6583', action="store_true", dest="NII6583", default=None, help="keyword to fit NII@6583.45 line")
    parser.add_argument('--SII6716', action="store_true", dest="SII6716", default=None, help="keyword to fit SII@6716.44 line")
    parser.add_argument('--SII6731', action="store_true", dest="SII6731", default=None, help="keyword to fit SII@6730.82 line")
bepinat's avatar
bepinat committed
1033
    parser.add_argument('--OIII4363', action="store_true", dest="OIII4363", default=None, help="keyword to fit OIII@4363 line")
bepinat's avatar
bepinat committed
1034 1035 1036
    parser.add_argument('--OIII4959', action="store_true", dest="OIII4959", default=None, help="keyword to fit OIII@4958.911 line")
    parser.add_argument('--OIII5007', action="store_true", dest="OIII5007", default=None, help="keyword to fit OIII@5006.843 line")
    parser.add_argument('--OII', action="store_true", dest="OII3729", default=None, help="keyword to fit OII doublet @3726.04,3728.80") # ['OII3729','OII3726']
bepinat's avatar
bepinat committed
1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068
    parser.add_argument('--OI6300', action="store_true", dest="OI6300", default=None, help="keyword to fit OI@6300 line")
    parser.add_argument('--HEI4471', action="store_true", dest="HeI4471", default=None, help="keyword to fit HeI@4471 line")
    parser.add_argument('--HEI5876', action="store_true", dest="HeI5876", default=None, help="keyword to fit HeI@5876 line")
    parser.add_argument('--HEII4686', action="store_true", dest="HeII4686", default=None, help="keyword to fit HeII@4686 line")
    parser.add_argument('--NEIII3868', action="store_true", dest="NeIII3868", default=None, help="keyword to fit NeIII@3868 line")
    parser.add_argument('--EXTRAL', action="append", dest="EXTRAL", default=None, type=float, help="keyword to add a new line (wavelength in angstrom) - multiple extra wavelengths can be added")
    parser.add_argument('--COMMW', action="store_true", dest="COMMW", default=None, help="keyword to fix common widths (in velocity) between several fitted lines")
    parser.add_argument('--REDSHIFT', action="store", dest="REDSHIFT", default=None, type=float, help="mean redshift of the object")
    parser.add_argument('--REDMIN', action="store", dest="REDMIN", default=None, type=float, help="minimum redshift allowed for fitting")
    parser.add_argument('--REDMAX', action="store", dest="REDMAX", default=None, type=float, help="maximum redshift allowed for fitting")
    parser.add_argument('--INITW', action="store", dest="INITW", default=None, type=float, help="initial linewidth in km/s")
    parser.add_argument('--WMIN', action="store", dest="WMIN", default=None, type=float, help="minimum width allowed for fitting in km/s")
    parser.add_argument('--WMAX', action="store", dest="WMAX", default=None, type=float, help="maximum width allowed for fitting in km/s")
    parser.add_argument('--DFIT', action="store", dest="DFIT", default=None, type=float, help="bin in km/s for searching the position of the lines")
    parser.add_argument('--DGCTNUM', action="store", dest="DGCTNUM", default=None, type=int, help="continuum degree")
    parser.add_argument('--SCLIP', action="store", dest="SCLIP", default=None, type=float, help="sigma-clipping threshold to clean the cube")
    parser.add_argument('--XYCLIP', action="store", dest="XYCLIP", default=None, type=int, help="width of the box to compute the median for the sigma-clipping")
    parser.add_argument('--NCLIP', action="store", dest="NCLIP", default=None, type=int, help="number of pass for clipping")
    parser.add_argument('--LSF', action="store", dest="SPSF", default=None, type=float, help="spectral PSF width in angstroms (dispersion of the gaussian)")
    parser.add_argument('--WSMOOTH', action="store", dest="WSMOOTH", default=None, type=float, help="width for Gaussian spectral smoothing in pixel") # options? gaussian uniquement?
    parser.add_argument('--SSMOOTH', action="store", dest="SSMOOTH", default=None, type=float, help="width for Gaussian spatial smoothing in pixel")
    parser.add_argument('--XMIN', action="store", dest="XMIN", default=None, type=int, help="lower cut window in x direction of the cube in pixel")
    parser.add_argument('--YMIN', action="store", dest="YMIN", default=None, type=int, help="upper cut window in x direction of the cube in pixel")
    parser.add_argument('--ZMIN', action="store", dest="ZMIN", default=None, type=int, help="lower cut window in y direction of the cube in pixel")
    parser.add_argument('--XMAX', action="store", dest="XMAX", default=None, type=int, help="upper cut window in y direction of the cube in pixel")
    parser.add_argument('--YMAX', action="store", dest="YMAX", default=None, type=int, help="lower cut window in z direction of the cube in pixel")
    parser.add_argument('--ZMAX', action="store", dest="ZMAX", default=None, type=int, help="upper cut window in z direction of the cube in pixel")

    #parser.add_argument('--MFIT', action="store", dest="MFIT", default=None)
    #parser.add_argument('--THRES', action="store", dest="THRES", default=None, type=float)
    #parser.add_argument('--MEDIAN', action="store", dest="MEDIAN", default=None, type=float)
    #parser.add_argument('--FITSPSF', action="store", dest="FITSPSF", default=None, type=float)
1069 1070
    
    try:
bepinat's avatar
bepinat committed
1071
        options = parser.parse_args()
1072 1073 1074
    except:
        sys.exit(2)
    
bepinat's avatar
bepinat committed
1075
    if np.size(argv) < 1:
1076
        sys.exit(2)
bepinat's avatar
bepinat committed
1077
    
1078 1079 1080 1081 1082 1083 1084
    opts = vars(options)   # this is a dictionnary
    opts['OII3726'] = opts['OII3729']
    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:])