gigrating.c

/*
 * This file is part of the GIRAFFE Pipeline
 * Copyright (C) 2002-2019 European Southern Observatory
 *
 * 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 2 of the License, or
 * (at your option) 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, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */
 
#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif
 
#include <math.h>
 
#include <cxmemory.h>
#include <cxmessages.h>
 
#include "gimacros.h"
#include "gialias.h"
#include "gigrating.h"
 
#define GIFITS_KEYWORD_MISSING_MSG "FITS KEYWORD [%s] not found!! Aborting..."
 


 
static const cxdouble GI_WAVELENGTH_EPSILON = 1.e-8;  /* 1.e-4 */
 
 
inline static cxint
_giraffe_grating_copy(GiGrating *self, cpl_table *grating, cxint row,
                      GiInstrumentMode mode)
{
 
    const cxchar *c_resolution = NULL;
 
 
    cx_assert(self != NULL);
    cx_assert(grating != NULL);
 
 
    if (row >= cpl_table_get_nrow(grating)) {
        return 1;
    }
 
 
    if (!cpl_table_has_column(grating, "ORDER")) {
        return 2;
    }
    else {
        self->order = cpl_table_get_int(grating, "ORDER", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "WLEN0")) {
        return 2;
    }
    else {
        self->wlen0 = cpl_table_get_double(grating, "WLEN0", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "WLMIN")) {
        return 2;
    }
    else {
        self->wlenmin = cpl_table_get_double(grating, "WLMIN", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "WLMAX")) {
        return 2;
    }
    else {
        self->wlenmax = cpl_table_get_double(grating, "WLMAX", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "BAND")) {
        return 2;
    }
    else {
        self->band = cpl_table_get_double(grating, "BAND", row, NULL);
    }
 
 
    switch (mode) {
        case GIMODE_MEDUSA:
            c_resolution = "RMED";
            break;
 
        case GIMODE_IFU:
        case GIMODE_ARGUS:
            c_resolution = "RIFA";
            break;
 
        default:
            return 3;
            break;
    }
 
    if (!cpl_table_has_column(grating, c_resolution)) {
        return 2;
    }
    else {
        self->resol = cpl_table_get_int(grating, c_resolution, row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "THETA")) {
        return 2;
    }
    else {
        self->theta = cpl_table_get_double(grating, "THETA", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "FCOLL")) {
        return 2;
    }
    else {
        self->fcoll = cpl_table_get_double(grating, "FCOLL", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "GCAM")) {
        return 2;
    }
    else {
        self->gcam = cpl_table_get_double(grating, "GCAM", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "SDX")) {
        return 2;
    }
    else {
        self->sdx = cpl_table_get_double(grating, "SDX", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "SDY")) {
        return 2;
    }
    else {
        self->sdy = cpl_table_get_double(grating, "SDY", row, NULL);
    }
 
 
    if (!cpl_table_has_column(grating, "SPHI")) {
        return 2;
    }
    else {
        self->sphi = cpl_table_get_double(grating, "SPHI", row, NULL);
    }
 
    return 0;
 
}
 

 
GiGrating *
giraffe_grating_new(void)
{
 
    GiGrating *grating = (GiGrating *) cx_calloc(1, sizeof *grating);
 
    grating->name   = cx_string_new();
    grating->setup  = cx_string_new();
    grating->filter = cx_string_new();
    grating->slit   = cx_string_new();
 
    return grating;
 
}
 

 
GiGrating *
giraffe_grating_create(const GiImage *spectra, const GiTable *grating)
{
 
    const cxchar *const c_setup = "SETUP";
 
    const cxchar *setup = NULL;
 
    cxint i;
    cxint order;
    cxint row = -1;
    cxint status = 0;
 
    cxdouble wlen0;
 
    cpl_propertylist *properties = NULL;
 
    cpl_table *_grating = NULL;
 
    GiGrating *self = NULL;
 
    GiInstrumentMode mode;
 
 
 
    if (spectra == NULL || grating == NULL) {
        return NULL;
    }
 
    properties = giraffe_image_get_properties(spectra);
 
    if (properties == NULL) {
        return NULL;
    }
 
    _grating = giraffe_table_get(grating);
 
    if (_grating == NULL) {
        return NULL;
    }
 
 
    /*
     * Get instrument setup information from the reference image
     * and fill the grating setup structure with the appropriate
     * data from the grating table. The central wavelength and the
     * order information provided by the reference image is used
     * to select the grating data.
     */
 
    self = giraffe_grating_new();
 
    /* Grating name */
 
    if (!cpl_propertylist_has(properties, GIALIAS_GRATNAME)) {
        giraffe_grating_delete(self);
        return NULL;
    }
    else {
        cx_string_set(self->name,
                      cpl_propertylist_get_string(properties,
                                                  GIALIAS_GRATNAME));
    }
 
    /* Order sorting filter name */
 
    if (!cpl_propertylist_has(properties, GIALIAS_FILTNAME)) {
        giraffe_grating_delete(self);
        return NULL;
    }
    else {
        cx_string_set(self->filter,
                      cpl_propertylist_get_string(properties,
                                                  GIALIAS_FILTNAME));
    }
 
    /* Slit name */
 
    if (!cpl_propertylist_has(properties, GIALIAS_SLITNAME)) {
        giraffe_grating_delete(self);
        return NULL;
    }
    else {
        cx_string_set(self->slit,
                      cpl_propertylist_get_string(properties,
                                                  GIALIAS_SLITNAME));
    }
 
    /* Spacing of grating grooves */
 
    if (!cpl_propertylist_has(properties, GIALIAS_GRATGRV)) {
        giraffe_grating_delete(self);
        return NULL;
    }
    else {
 
        cxdouble grooves = cpl_propertylist_get_double(properties,
                                                       GIALIAS_GRATGRV);
        self->space = 1. / fabs(GI_MM_TO_NM * grooves);
 
    }
 
    /* Instrument mode */
 
    mode = giraffe_get_mode(properties);
 
 
    /*
     * Select the grating data for the current setup
     */
 
    if (!cpl_table_has_column(_grating, "ORDER") ||
        !cpl_table_has_column(_grating, "WLEN0")) {
 
        giraffe_grating_delete(self);
        return NULL;
 
    }
 
    if (!cpl_propertylist_has(properties, GIALIAS_GRATWLEN)) {
        giraffe_grating_delete(self);
        return NULL;
    }
    else {
        wlen0 = cpl_propertylist_get_double(properties, GIALIAS_GRATWLEN);
    }
 
    if (!cpl_propertylist_has(properties, GIALIAS_GRATORDER)) {
        giraffe_grating_delete(self);
        return NULL;
    }
    else {
        order = cpl_propertylist_get_int(properties, GIALIAS_GRATORDER);
    }
 
 
    for (i = 0; i < cpl_table_get_nrow(_grating); i++) {
 
        cxint _order = cpl_table_get_int(_grating, "ORDER", i, NULL);
 
        if (order == _order) {
 
            cxdouble _wlen0 = cpl_table_get_double(_grating, "WLEN0",
                                                   i, NULL);
 
            if (fabs(wlen0 - _wlen0) < GI_WAVELENGTH_EPSILON) {
                row = i;
                break;
            }
 
        }
 
    }
 
    if (row < 0) {
        giraffe_grating_delete(self);
        return NULL;
    }
 
    /*
     * Try to figure out the setup identifier.
     */
 
    /* FIXME: Check whether we can live with an empty setup string
     *        at later stages of the wavelength calibration.
     */
 
    if (cpl_propertylist_has(properties, GIALIAS_SETUPNAME)) {
        setup = cpl_propertylist_get_string(properties, GIALIAS_SETUPNAME);
        cx_string_set(self->setup, setup);
    }
    else {
        if (cpl_table_has_column(_grating, c_setup)) {
            setup = cpl_table_get_string(_grating, c_setup, row);
            cx_string_set(self->setup, setup);
        }
    }
 
 
    status = _giraffe_grating_copy(self, _grating, row, mode);
 
    if (status != 0) {
        giraffe_grating_delete(self);
        return NULL;
    }
 
    return self;
 
}
 

 
void
giraffe_grating_delete(GiGrating *self)
{
 
    if (self == NULL) {
        return;
    }
 
    if (self->name != NULL) {
        cx_string_delete(self->name);
    }
 
    if (self->setup != NULL) {
        cx_string_delete(self->setup);
    }
 
    if (self->filter != NULL) {
        cx_string_delete(self->filter);
    }
 
    if (self->slit != NULL) {
        cx_string_delete(self->slit);
    }
 
    cx_free(self);
 
    return;
 
}
 

 
cxint
giraffe_grating_setup(GiTable *grating_table, GiImage *spectra,
                      GiGrating *grating_setup)
{
 
    const cxchar* const fctid = "giraffe_grating_setup";
 
    const cxchar* c_name_setup  = "SETUP";
    const cxchar* c_name_order  = "ORDER";
    const cxchar* c_name_wl0    = "WLEN0";
    const cxchar* c_name_wlmin  = "WLMIN";
    const cxchar* c_name_wlmax  = "WLMAX";
    const cxchar* c_name_band   = "BAND";
    const cxchar* c_name_theta  = "THETA";
    const cxchar* c_name_fcoll  = "FCOLL";
    const cxchar* c_name_gcam   = "GCAM";
    const cxchar* c_name_sdx    = "SDX";
    const cxchar* c_name_sdy    = "SDY";
    const cxchar* c_name_sdphi  = "SPHI";
    const cxchar* c_name_rmed   = "RMED";
    const cxchar* c_name_rifa   = "RIFA";
 
    cxint i = 0;
    cxint row_match = 0;
    cxint row_nulls = 0;
 
    cxdouble wlen_match = 0.;
    cxdouble wlen = 0.;
    cxdouble tmp_gratgrv = 0.;
 
 
    cx_string* slit_name = NULL;
 
    cpl_propertylist* ref_plimg = NULL;
 
    cpl_table* ref_gtable = NULL;
 
    GiInstrumentMode  instrument_mode;
 
 
    /************************************************************************
                                    Preprocessing
    ************************************************************************/
 
    if (grating_table  == NULL) {
        return 1;
    }
 
    if (spectra == NULL) {
        return 1;
    }
 
    if (grating_setup == NULL) {
        return 1;
    }
 
    ref_plimg = giraffe_image_get_properties(spectra);
    if (ref_plimg == NULL) {
        return 128;
    }
 
    ref_gtable = giraffe_table_get(grating_table);
    if (ref_gtable == NULL) {
        return 128;
    }
 
    slit_name = cx_string_new();
 
    /************************************************************************
                                     Processing
    ************************************************************************/
 
    /*
     *  Retrieve Grating information from associated image...
     */
 
    if (cpl_propertylist_has(ref_plimg, GIALIAS_GRATWLEN) == FALSE) {
        cpl_msg_error(fctid, GIFITS_KEYWORD_MISSING_MSG, GIALIAS_GRATWLEN);
        cx_string_delete(slit_name);
        return 2;
    }
    else {
        grating_setup->wlen0 = cpl_propertylist_get_double(ref_plimg,
                                                           GIALIAS_GRATWLEN);
    }
 
    if (cpl_propertylist_has(ref_plimg, GIALIAS_SLITNAME) == FALSE) {
        cpl_msg_error(fctid, GIFITS_KEYWORD_MISSING_MSG, GIALIAS_SLITNAME);
        cx_string_delete(slit_name);
        return 2;
    }
    else {
        cx_string_set(slit_name,
                      cpl_propertylist_get_string(ref_plimg,
                                                   GIALIAS_SLITNAME));
    }
 
    if (cpl_propertylist_has(ref_plimg, GIALIAS_GRATGRV) == FALSE) {
        cpl_msg_error(fctid, GIFITS_KEYWORD_MISSING_MSG, GIALIAS_GRATGRV);
        cx_string_delete(slit_name);
        return 2;
    }
    else {
        tmp_gratgrv = cpl_propertylist_get_double(ref_plimg,
                                                  GIALIAS_GRATGRV);
    }
 
    if (cpl_propertylist_has(ref_plimg, GIALIAS_GRATNAME) == FALSE) {
        cpl_msg_error(fctid, GIFITS_KEYWORD_MISSING_MSG, GIALIAS_GRATNAME);
        cx_string_delete(slit_name);
        return 2;
    }
    else {
        cx_string_set(grating_setup->name,
                      cpl_propertylist_get_string(ref_plimg,
                                                  GIALIAS_GRATNAME));
    }
 
    if (cpl_propertylist_has(ref_plimg, GIALIAS_FILTNAME) == FALSE) {
        cpl_msg_error(fctid, GIFITS_KEYWORD_MISSING_MSG, GIALIAS_FILTNAME);
        cx_string_delete(slit_name);
        return 2;
    }
    else {
        cx_string_set(grating_setup->filter,
                      cpl_propertylist_get_string(ref_plimg,
                                                  GIALIAS_FILTNAME));
    }
 
 
    /*
     *  Find wavelength nearest to central wavelength...
     */
 
    for (i = 0; i < cpl_table_get_nrow(ref_gtable); i++) {
 
        wlen = cpl_table_get(ref_gtable, c_name_wl0, i, &row_nulls);
 
        if (fabs(wlen - grating_setup->wlen0) <
            fabs(wlen_match - grating_setup->wlen0)) {
            wlen_match = wlen;
            row_match  = i;
        }
 
    }
 
 
    /*
     *  Have we found a match?...
     */
 
    if (fabs(wlen_match - grating_setup->wlen0) > GI_WAVELENGTH_EPSILON) {
 
        cpl_msg_error(fctid, "Central wavelength [%f] nout found in grating "
                      "table, aborting...", grating_setup->wlen0);
        cx_string_delete(slit_name);
        return 3;
    }
    else {
        cpl_msg_debug(fctid, "Found wlen0 in grating table at position %d",
                      row_match);
    }
 
 
    /*
     * Retrieve values associated to matched wavelength from
     * grating table...
     */
 
    cx_string_set(grating_setup->setup,
                  (cxchar*) cpl_table_get_string(ref_gtable, c_name_setup,
                                                 row_match));
 
    cx_string_set(grating_setup->slit, cx_string_get(slit_name));
 
    grating_setup->order = cpl_table_get(ref_gtable, c_name_order,
                                         row_match, &row_nulls);
 
    grating_setup->wlenmin = cpl_table_get(ref_gtable, c_name_wlmin,
                                           row_match, &row_nulls);
 
    grating_setup->wlenmax = cpl_table_get(ref_gtable, c_name_wlmax,
                                           row_match, &row_nulls);
 
    grating_setup->band = cpl_table_get(ref_gtable, c_name_band,
                                        row_match, &row_nulls);
 
    grating_setup->theta = cpl_table_get(ref_gtable, c_name_theta,
                                         row_match, &row_nulls);
 
    grating_setup->space = 1.0 / fabs(GI_MM_TO_NM * tmp_gratgrv);
 
 
    instrument_mode = giraffe_get_mode(ref_plimg);
 
    switch (instrument_mode) {
        case GIMODE_MEDUSA:
            grating_setup->resol = cpl_table_get(ref_gtable, c_name_rmed,
                                                 row_match, &row_nulls);
            break;
 
        case GIMODE_IFU:
            grating_setup->resol = cpl_table_get(ref_gtable, c_name_rifa,
                                                 row_match, &row_nulls);
            break;
 
        case GIMODE_ARGUS:
            grating_setup->resol = cpl_table_get(ref_gtable, c_name_rifa,
                                                 row_match, &row_nulls);
            break;
 
        default:
            grating_setup->resol = -1.0;
            break;
    }
 
    grating_setup->fcoll   =
        cpl_table_get(ref_gtable, c_name_fcoll, row_match, &row_nulls);
 
    grating_setup->gcam    =
        cpl_table_get(ref_gtable, c_name_gcam,  row_match, &row_nulls);
 
    grating_setup->sdx  =
        cpl_table_get(ref_gtable, c_name_sdx,   row_match, &row_nulls);
 
    grating_setup->sdy  =
        cpl_table_get(ref_gtable, c_name_sdy,   row_match, &row_nulls);
 
    grating_setup->sphi =
        cpl_table_get(ref_gtable, c_name_sdphi, row_match, &row_nulls);
 
    cx_string_delete(slit_name);
 
    return 0;
 
}
 

 
void
giraffe_grating_dump(const GiGrating *grating)
{
 
    const cxchar *fctid = "giraffe_grating_dump";
 
    if (grating == NULL) {
        return;
    }
 
    cpl_msg_debug(fctid, "---- GiGrating -------------------------");
 
    cpl_msg_debug(fctid, "Grating Name        : %s",
                  cx_string_get(grating->name));
    cpl_msg_debug(fctid, "Grating Filter Name : %s",
                  cx_string_get(grating->filter));
    cpl_msg_debug(fctid, "Grating Setup Name  : %s",
                  cx_string_get(grating->setup));
    cpl_msg_debug(fctid, "Grating Order       : %12d", grating->order);
    cpl_msg_debug(fctid, "Grating Wlen0       : %12.6f", grating->wlen0);
    cpl_msg_debug(fctid, "Grating Wlen Min    : %12.6f", grating->wlenmin);
    cpl_msg_debug(fctid, "Grating Wlen Max    : %12.6f", grating->wlenmax);
    cpl_msg_debug(fctid, "Grating Band        : %12.6f", grating->band);
    cpl_msg_debug(fctid, "Grating Resol       : %12d", grating->resol);
    cpl_msg_debug(fctid, "Grating Space       : %12.6f", grating->space);
    cpl_msg_debug(fctid, "Grating Theta       : %12.6f", grating->theta);
    cpl_msg_debug(fctid, "Grating FColl       : %12.6f", grating->fcoll);
    cpl_msg_debug(fctid, "Grating GCam        : %12.6f", grating->gcam);
    cpl_msg_debug(fctid, "Grating SlitDx      : %12.6f", grating->sdx);
    cpl_msg_debug(fctid, "Grating SlitDy      : %12.6f", grating->sdy);
    cpl_msg_debug(fctid, "Grating SlitPhi     : %12.6f", grating->sphi);
 
    return;
 
}