wave_cal_slit.c

00001 /*----------------------------------------------------------------------------
00002  
00003    File name    :       wave_cal_slit.c
00004    Author       :   A. Modigliani
00005    Created on   :   Sep 17, 2003
00006    Description  : 
00007 
00008  ---------------------------------------------------------------------------*/
00009 
00010 /*----------------------------------------------------------------------------
00011                                 Includes
00012  ---------------------------------------------------------------------------*/
00013 #include "wave_cal_slit.h"
00014 #include "sinfoni_pro_save.h"
00015 #include "wavecal_ini_by_cpl.h"
00016 #include "sinfoni_wcal_functions.h"
00017 /*----------------------------------------------------------------------------
00018                                 Defines
00019  ---------------------------------------------------------------------------*/
00020 
00021 /*----------------------------------------------------------------------------
00022                              Function Definitions
00023  ---------------------------------------------------------------------------*/
00024 int dumpFitParamsToTbl ( FitParams ** params, char * filename );
00025 
00026 /*----------------------------------------------------------------------------
00027    Function     :       wave_cal_slit()
00028    In           :       ini_file: file name of according .ini file
00029    Out          :       integer (0 if it worked, -1 if it doesn't) 
00030    Job          :
00031 
00032 
00033   Normal method:
00034 
00035   does the wavelength calibration and the fitting of the slitlet edge 
00036   positions (ASCII file 32 x 2 values) if wished
00037   produces an array of the bcoefs and of the fit parameters if wished and a 
00038   wavelength calibration map input is an emission line frame and a line list
00039 
00040 
00041   o searching for lines by cross correlation with a line list
00042   o Gaussian fitting of emission lines in each column->positions of the lines->
00043     resulting fit parameters can be stored in an ASCII file
00044   o Fitting of a polynomial to the line positions for each column
00045   o Smoothing: fitting of each polynomial coefficient by another polynomial
00046     across the whole frame -> resulting polynomial coefficients can be stored 
00047     in an ASCII file.
00048   o Wavelength calibration map (micron value for each frame pixel) can be
00049     produced by using these coefficients and a cross correlation to the
00050     original frame
00051 
00052   o The slitlet edge positions can be fitted:
00053     1) Automatically (not really stable) or by using guess edge positions
00054     2) By using a Boltzmann or a linear slope function
00055 
00056   Slit method:
00057 
00058   does the wavelength calibration and the fitting of the slitlet edge 
00059   positions (ASCII file 32 x 2 values) if wished produces a list of the fit 
00060   parameters and of the smoothed coefficients if wished and a wavelength 
00061   calibration map input is an emission line frame and a line list
00062 
00063   o Does the same as other method but smoothes the found polynomial 
00064     coefficients within each slitlet and not over the whole frame.
00065 
00066   o Produces always a wavelength calibration map and does not crosscorrelate.
00067 
00068  ---------------------------------------------------------------------------*/
00069 
00070 
00071 int wave_cal_slit (const char* plugin_id,cpl_parameterlist* config, cpl_frameset* sof)
00072 {
00073     const char* _id = "wave_cal_slit";
00074     wave_config * cfg =NULL;
00075 
00076 
00077 
00078 
00079     char col_name[MAX_NAME_SIZE];
00080     char tbl_name[MAX_NAME_SIZE];
00081 
00082     char tbl_fitpar_name[MAX_NAME_SIZE];
00083     char tbl_line_list_name[MAX_NAME_SIZE];
00084     char tbl_slitpos_guess_name[MAX_NAME_SIZE];
00085     char* col=NULL;
00086     char key_value[MAX_NAME_SIZE];
00087     char key_name[MAX_NAME_SIZE];
00088 
00089     int check = 0;
00090     int lx = 0;
00091     int ly = 0;
00092     int n_lines=0;
00093     int i = 0;
00094     int j = 0;
00095     int n = 0;
00096     int readsum=0;
00097     int sum=0;
00098     int fit=0;
00099     int sw=0;
00100     int trow=0;
00101 
00102     int* status=NULL;
00103     int* n_found_lines=NULL;
00104     int* sum_pointer=NULL;
00105     int** row_clean;
00106 
00107     float a=0;
00108     float shift=0;
00109     float val_x=0;
00110     float val_y=0;
00111 
00112     float* wave=NULL;
00113     float* intens=NULL;
00114 
00115     float** acoefs;
00116     float** wavelength_clean;
00117     float** slit_pos;
00118 
00119     double fwhm_med=0;
00120     double fwhm_avg=0;
00121     double coef_med=0;
00122     double coef_avg=0;
00123  
00124     OneImage *  map=NULL ;
00125     OneImage * im=NULL ;
00126 
00127     FitParams** par;
00128 
00129     cpl_table* tbl_wcal=NULL;
00130     cpl_table* tbl_spos=NULL;
00131     cpl_table* tbl_fitpar = NULL;
00132     cpl_table* tbl_line_list = NULL;
00133     cpl_table* tbl_slitpos_guess=NULL;
00134     cpl_table * tbl_fp =NULL;
00135     cpl_table* qclog_tbl=NULL;
00136 
00137     cpl_image* map_img=NULL;
00138     cpl_image* map_imgw=NULL;
00139 
00140     cpl_frameset* raw=NULL;
00141     cpl_parameter* p=NULL;
00142 
00143     struct qc_wcal qc;
00144     /*        -----------------------------------------------------------------
00145        1) parse the file names and parameters to the ns_config data 
00146           structure cfg
00147        -----------------------------------------------------------------
00148      */
00149 
00150     cpl_msg_info(_id,"Parsing cpl input");
00151     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00152       cpl_msg_error(_id, (char* ) cpl_error_get_message());
00153       return -1;
00154     }
00155     raw=cpl_frameset_new();
00156 
00157     cfg = parse_cpl_input_wave(config,sof,&raw) ;
00158 
00159     if (cfg == NULL)
00160     {
00161         cpl_msg_error(_id,"could not parse cpl input!\n") ;
00162         cpl_frameset_delete(raw);
00163         return -1 ;
00164     }
00165 
00166     p = cpl_parameterlist_find(config,"sinfoni.wavecal.slitpos_boostrap");
00167     sw=cpl_parameter_get_bool(p);
00168    
00169 
00170     if(sw==1) {
00171        cfg->nslitlets=32;
00172        cfg->calibIndicator=1;
00173        cfg->wavemapInd=0;
00174        cfg->slitposIndicator=1;
00175        cpl_msg_info(_id,"***********************************");
00176        cpl_msg_info(_id,"parameter setting for %s",PRO_WAVE_SLITPOS_STACKED);
00177        cpl_msg_info(_id,"***********************************");
00178     }
00179 
00180     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00181       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00182       return -1;
00183     }
00184    
00185 
00186     if(is_fits_file(cfg->inFrame) != 1) {
00187       cpl_msg_error(_id,"Input file cfg->inFrame %s is not FITS",cfg->inFrame);
00188       cpl_frameset_delete(raw);
00189       wavecal_free(cfg);
00190       return -1;
00191     }
00192 
00193 
00194     if (cfg->slitposIndicator == 1 && cfg->estimateIndicator == 1) {
00195       if (is_fits_file(cfg->slitposGuessName) != 1) {
00196         cpl_msg_error(_id,"slitlet position guess list not given!");
00197         cpl_frameset_delete(raw);
00198         wavecal_free(cfg);
00199         return -1;
00200       }
00201     }
00202 
00203     if (cfg->calibIndicator == 0 && cfg->wavemapInd == 1) {
00204       if (is_fits_file(cfg->coeffsName) != 1) {
00205         cpl_msg_error(_id,"coefficients list not given!");
00206         cpl_frameset_delete(raw);
00207         wavecal_free(cfg);
00208         return -1;
00209       }
00210     }
00211 
00212 
00213     if (cfg->slitposIndicator == 1) {
00214       if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00215        
00216         if (is_fits_file(cfg->paramsList) != 1) {
00217       cpl_msg_error(_id,"parameter list not given!");
00218           cpl_frameset_delete(raw);
00219           wavecal_free(cfg);
00220       return -1;
00221     }
00222     
00223       }
00224     }
00225 
00226 
00227 /*---load the emission line frame---*/
00228     im = load_image(cfg->inFrame);
00229 
00230     if (im == NULL) {
00231       cpl_msg_error(_id,"could not load image\n");
00232       cpl_frameset_delete(raw);
00233       wavecal_free(cfg);
00234       return -1;
00235     }
00236 
00237 
00238     lx = im->lx;
00239     ly = im->ly;
00240 
00241 
00242 
00243 
00244 if (cfg->calibIndicator == 1 || cfg->wavemapInd == 1) {
00245     /*---open the line list and read the number of lines---*/
00246 
00247 
00248     strcpy(tbl_line_list_name,cfg->lineList);
00249     tbl_line_list = cpl_table_load(tbl_line_list_name,1,0);
00250     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00251       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00252       destroy_image(im);
00253       cpl_frameset_delete(raw);
00254       wavecal_free(cfg);
00255       return -1;
00256     }
00257     n = cpl_table_get_nrow(tbl_line_list);
00258     n_lines = n;
00259     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00260       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00261       cpl_table_delete(tbl_line_list);
00262       destroy_image(im);
00263       cpl_frameset_delete(raw);
00264       wavecal_free(cfg);
00265       return -1;
00266     }
00267 
00268     /* THIS ORIGINATES A MEMORY LEAK 
00269     wave   = new_float_array (n);
00270     intens = new_float_array (n);
00271     if (wave == NULL || intens == NULL) {
00272       cpl_msg_error(_id,"could not allocate memory for the line list values\n" );
00273       return -1;
00274     }
00275     */
00276 
00277     wave   = cpl_table_get_data_float(tbl_line_list,"wave");
00278     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00279       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00280       cpl_table_delete(tbl_line_list);
00281       destroy_image(im);
00282       cpl_frameset_delete(raw);
00283       wavecal_free(cfg);
00284       return -1;
00285     }
00286 
00287     intens = cpl_table_get_data_float(tbl_line_list,"int");
00288     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00289       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00290       cpl_table_delete(tbl_line_list);
00291       destroy_image(im);
00292       cpl_frameset_delete(raw);
00293       wavecal_free(cfg);
00294       return -1;
00295     }
00296 
00297 }
00298 
00299 
00300 /*
00301  ----------------------------------------------------------------------
00302  ---------------------------FINDLINES----------------------------------
00303  ----------------------------------------------------------------------
00304  */
00305 
00306 /*if not yet done:
00307   do the wavelength calibration, that means: 
00308   find the dispersion relation and parameterize its coefficients
00309  */
00310  cpl_msg_info(_id,"guessBeginWave=%g",cfg->guessBeginWavelength);
00311  cpl_msg_info(_id,"guessDisp1=%g",cfg->guessDispersion1);
00312  cpl_msg_info(_id,"guessDisp2=%g",cfg->guessDispersion2);
00313 
00314 
00315 
00316 
00317 if (cfg->calibIndicator == 1 && cfg->wavemapInd == 0) {
00318    cpl_msg_info(_id,"Findlines");
00319    acoefs  = new_2Dfloat_array(cfg->nrDispCoefficients, lx);
00320    /*allocate memory*/
00321    n_found_lines    = new_int_array(lx); 
00322    row_clean        = new_2Dint_array(lx, n_lines);
00323    wavelength_clean = new_2Dfloat_array(lx, n_lines);
00324    sum_pointer      = new_int_array(1) ;
00325 
00326 
00327 
00328    /*find the emission lines in each image column*/
00329    intarray_set_value(sum_pointer, 0, 0);
00330    check = findLines(im, wave, intens, n_lines, row_clean, 
00331                       wavelength_clean, cfg->guessBeginWavelength, 
00332              cfg->guessDispersion1, cfg->guessDispersion2,
00333                       cfg->mindiff, cfg->halfWidth, 
00334                       n_found_lines, cfg->sigma, sum_pointer );
00335 
00336    if (-1 == check) {
00337      cpl_msg_error(_id,"findLines failed!\n");
00338      destroy_2Darray ( wavelength_clean, lx );
00339      destroy_2Dintarray (row_clean, lx);
00340      destroy_intarray ( n_found_lines );
00341      destroy_intarray ( sum_pointer );
00342      destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00343      cpl_table_delete(tbl_line_list);
00344      destroy_image(im);
00345      cpl_frameset_delete(raw);
00346      wavecal_free(cfg);
00347      return -1;
00348    }
00349 
00350 
00351 
00352       
00353 
00354 /*---------------------------------------------------------------------
00355  *-----------------------WAVE_CALIB-------------------------------------
00356  *---------------------------------------------------------------------
00357  */
00358 
00359      
00360    cpl_msg_info(_id,"Wave Calibration");
00361    sum = intarray_get_value(sum_pointer,0);
00362    /* allocate memory for the fit parameters */
00363    par = newFitParams( sum );
00364 
00365    if (par == NULL) {
00366         cpl_msg_error(_id,"newFitParams failed!\n");
00367     destroy_2Darray ( wavelength_clean, lx );
00368     destroy_2Dintarray (row_clean, lx);
00369     destroy_intarray ( n_found_lines );
00370     destroy_intarray ( sum_pointer );
00371     destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00372     destroyFitParams(par);
00373     cpl_table_delete(tbl_line_list);
00374     destroy_image(im);
00375     cpl_frameset_delete(raw);
00376     wavecal_free(cfg);
00377         return -1;
00378    }
00379 
00380   /*
00381    fit each line, make a polynomial fit and fit the resulting fit 
00382    coefficients across the columns of the slitlet
00383    */
00384    map = waveCal(im, 
00385                  par, 
00386                  acoefs, 
00387                  cfg->nslitlets, 
00388                  row_clean, 
00389                  wavelength_clean, 
00390                  n_found_lines, 
00391                  cfg->guessDispersion1, 
00392                  cfg->halfWidth, 
00393                  cfg->minAmplitude, 
00394                  cfg->maxResidual, 
00395                  cfg->fwhm, 
00396                  cfg->nrDispCoefficients, 
00397                  cfg->nrCoefCoefficients, 
00398                  cfg->sigmaFactor, 
00399                  cfg->pixeldist, 
00400                  cfg->pixel_tolerance);
00401 
00402 
00403    if (map == NULL ) { 
00404      cpl_msg_error(_id,"wave_cal failed!\n");
00405      destroy_2Darray ( wavelength_clean, lx );
00406      destroy_2Dintarray (row_clean, lx);
00407      destroy_intarray ( n_found_lines );
00408      destroy_intarray ( sum_pointer );
00409      destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00410      destroyFitParams(par);
00411      cpl_table_delete(tbl_line_list);
00412      destroy_image(im);
00413      cpl_frameset_delete(raw);
00414      wavecal_free(cfg);
00415      return -1;
00416    }
00417    cpl_msg_info(_id,"Check line positions");
00418   
00419    shift = checkLinePositions (im, acoefs, cfg->nrDispCoefficients, par);
00420    if (FLAG == shift){ 
00421       cpl_msg_error(_id,"checkForLinePositions failed!\n");
00422    }
00423 
00424    map_imgw=cpl_image_wrap_float(map->lx,map->ly,map->data);
00425    map_img=cpl_image_duplicate(map_imgw);
00426    cpl_image_unwrap(map_imgw);
00427 
00428    det_ncounts(raw, cfg->qc_thresh_max, &qc);
00429 
00430     qclog_tbl = sinfoni_qclog_init(4);
00431      sprintf(key_value,"%d",n_lines);
00432     sinfoni_qclog_add(qclog_tbl,0,"QC WAVE ALL","CPL_TYPE_INT",
00433                       key_value,"Number of found lines");
00434 
00435     sprintf(key_value,"%d",qc.nsat);
00436     sinfoni_qclog_add(qclog_tbl,1,"QC WAVE NPIXSAT","CPL_TYPE_INT",
00437                       key_value,"Number of saturated pixels");
00438 
00439     sprintf(key_value,"%g",qc.max_di);
00440     sinfoni_qclog_add(qclog_tbl,2,"QC WAVE MAXFLUX","CPL_TYPE_DOUBLE",
00441                       key_value,"Max int off-lamp subracted frm");
00442 
00443     sprintf(key_value,"%g",shift);
00444     sinfoni_qclog_add(qclog_tbl,3,"QC WAVE POSERR","CPL_TYPE_DOUBLE",
00445                       key_value,"Overall positioning error in microns");
00446 
00447 
00448 
00449 
00450 
00451    if(-1 == sinfoni_pro_save_ima(map_img,raw,sof,cfg->outName,
00452          PRO_WAVE_MAP,qclog_tbl,plugin_id,config)) {
00453          cpl_msg_error(_id,"cannot save ima %s", cfg->outName);
00454      destroy_2Darray ( wavelength_clean, lx );
00455      destroy_2Dintarray (row_clean, lx);
00456      destroy_intarray ( n_found_lines );
00457      destroy_intarray ( sum_pointer );
00458      destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00459      destroyFitParams(par);
00460      cpl_table_delete(tbl_line_list);
00461      destroy_image(im);
00462      destroy_image(map);
00463      cpl_frameset_delete(raw);
00464      wavecal_free(cfg);
00465      cpl_table_delete(qclog_tbl);
00466      cpl_image_delete(map_img);
00467          return -1;
00468     }
00469    cpl_table_delete(qclog_tbl);
00470    cpl_image_delete(map_img);
00471 
00472     /* destroy_image(map); */
00473 
00474    /*
00475     #store the resulting polynomial fit coefficients in an 
00476      ASCII file if wished
00477     */
00478 
00479 
00480 
00481    if (cfg->writeCoeffsInd == 1) {
00482          tbl_wcal = cpl_table_new(lx);
00483          for (i=0; i< cfg->nrDispCoefficients; i++) {
00484              sprintf(col_name,"%s%d","coeff",i);
00485              cpl_table_new_column(tbl_wcal,col_name, CPL_TYPE_DOUBLE);
00486      }
00487 
00488 
00489          qclog_tbl = sinfoni_qclog_init(1+2*cfg->nrDispCoefficients);
00490          sprintf(key_value,"%d",n_lines);
00491          sinfoni_qclog_add(qclog_tbl,0,"QC WAVE ALL","CPL_TYPE_INT",
00492                key_value,"Number of found lines");
00493 
00494 
00495 
00496          for (j=0; j< lx; j++) { 
00497         for (i=0; i< cfg->nrDispCoefficients; i++) {
00498                 sprintf(col_name,"%s%d","coeff",i);
00499             a = array2D_get_value(acoefs, i, j);
00500         /* fprintf(acoefs_file, "%15.13g ", a) ; */
00501                 cpl_table_set_double(tbl_wcal,col_name,j,a);
00502         }
00503 
00504 
00505      }
00506 
00507 
00508      for (i=0; i< cfg->nrDispCoefficients; i++) {
00509             sprintf(col_name,"%s%d","coeff",i);
00510             coef_avg=cpl_table_get_column_mean(tbl_wcal,col_name);
00511             coef_med=cpl_table_get_column_median(tbl_wcal,col_name);
00512 
00513             trow=1+i;
00514             sprintf(key_name,"%s%d%s","QC COEF",i," AVG");
00515             sprintf(key_value,"%g",coef_avg);
00516 
00517             sinfoni_qclog_add(qclog_tbl,trow,key_name,"CPL_TYPE_DOUBLE",
00518                key_value,"Average wavecal Coef");
00519 
00520             trow=1+i+cfg->nrDispCoefficients;
00521             sprintf(key_name,"%s%d%s","QC COEF",i," MED");
00522             sprintf(key_value,"%g",coef_med);
00523             sinfoni_qclog_add(qclog_tbl,trow,key_name,"CPL_TYPE_DOUBLE",
00524                key_value,"Median wavecal Coef");
00525 
00526      }
00527 
00528      /*
00529          fclose(acoefs_file);
00530      */
00531          strcpy(tbl_name,cfg->coeffsName);
00532          if(-1 == sinfoni_pro_save_tbl(tbl_wcal,raw,sof,tbl_name,
00533              PRO_WAVE_COEF_SLIT,qclog_tbl,plugin_id,config)) {
00534          cpl_msg_error(_id,"cannot save tbl %s", tbl_name);
00535 
00536          cpl_table_delete(tbl_wcal);
00537          cpl_table_delete(qclog_tbl);
00538              wavecal_free(cfg);
00539              cpl_frameset_delete(raw);
00540 
00541          destroy_2Darray ( wavelength_clean, lx );
00542          destroy_2Dintarray (row_clean, lx);
00543          destroy_intarray ( n_found_lines );
00544          destroy_intarray ( sum_pointer );
00545          destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00546          cpl_table_delete(tbl_line_list);
00547 
00548 
00549              if ( (cfg->slitposIndicator == 1 && 
00550                    cfg->estimateIndicator != 1) ||
00551                   (cfg->calibIndicator == 1)  || 
00552                   (cfg->wavemapInd == 1) ){
00553                   destroyFitParams(par);
00554          }
00555 
00556          destroy_image ( im );
00557          destroy_image ( map );
00558 
00559          return -1;
00560          }
00561          cpl_table_delete(tbl_wcal);
00562          cpl_table_delete(qclog_tbl);
00563 
00564    }
00565 
00566    /*
00567     #store the resulting Gaussian fit parameters in an ASCII file if wished
00568    */
00569    if (cfg->writeParInd == 1) {
00570 
00571 
00572       dumpFitParamsToAscii(par, cfg->paramsList);
00573  
00574       if ( NULL == par )
00575       {
00576          cpl_msg_error (_id,"no fit parameters available!") ;
00577          return -1;
00578       }
00579 
00580       if ( NULL == cfg->paramsList )
00581       {
00582          cpl_msg_error (_id,"no filename available!") ;
00583          return -1;
00584       }
00585 
00586       tbl_fp = cpl_table_new(par[0] -> n_params);
00587       cpl_table_new_column(tbl_fp,"n_params", CPL_TYPE_INT);
00588       cpl_table_new_column(tbl_fp,"column", CPL_TYPE_INT);
00589       cpl_table_new_column(tbl_fp,"line", CPL_TYPE_INT);
00590       col = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00591 
00592       for(j=0;j<4;j++) {
00593          sprintf(col,"%s%d","fpar",j);
00594          cpl_table_new_column(tbl_fp,col, CPL_TYPE_DOUBLE);
00595          sprintf(col,"%s%d","dpar",j);
00596          cpl_table_new_column(tbl_fp,col, CPL_TYPE_DOUBLE);
00597       }
00598 
00599 
00600 
00601     qclog_tbl = sinfoni_qclog_init(3);
00602     sprintf(key_value,"%d",n_lines);
00603     sinfoni_qclog_add(qclog_tbl,0,"QC NLINES","CPL_TYPE_INT",
00604               key_value,"Number of Found lines");
00605 
00606       for ( i = 0 ; i < par[0] -> n_params ; i++ )
00607       {
00608          cpl_table_set_int(tbl_fp,"n_params",i,par[i]->n_params);
00609          cpl_table_set_int(tbl_fp,"column",i,par[i]->column);
00610          cpl_table_set_int(tbl_fp,"line",i,par[i]->line);
00611 
00612 
00613          for(j=0;j<4;j++) {
00614        sprintf(col,"%s%d","fpar",j);
00615        if(qfits_isnan(par[i]->fit_par[j])) {
00616          cpl_table_set_invalid(tbl_fp,col,i);
00617        } else {
00618          cpl_table_set_double(tbl_fp,col,i,par[i]->fit_par[j]);
00619        }
00620        sprintf(col,"%s%d","dpar",j);
00621        if(qfits_isnan(par[i]->derv_par[j])) {
00622          cpl_table_set_invalid(tbl_fp,col,i);
00623        } else {
00624          cpl_table_set_double(tbl_fp,col,i,par[i]->derv_par[j]);
00625        }
00626      }
00627       }
00628 
00629       fwhm_avg = cpl_table_get_column_mean(tbl_fp,"fpar1");
00630       fwhm_med = cpl_table_get_column_median(tbl_fp,"fpar1");
00631       sprintf(key_value,"%f",fwhm_med);
00632  
00633       sinfoni_qclog_add(qclog_tbl,1,"QC FWHM MED","CPL_TYPE_DOUBLE",
00634               key_value,"Median FWHM of found lines");
00635 
00636       sprintf(key_value,"%f",fwhm_avg);
00637       sinfoni_qclog_add(qclog_tbl,2,"QC FWHM AVG","CPL_TYPE_DOUBLE",
00638               key_value,"Average FWHM of found lines");
00639 
00640       if(-1 == sinfoni_pro_save_tbl(tbl_fp,raw,sof,cfg->paramsList,
00641          PRO_WAVE_PAR_LIST,qclog_tbl,plugin_id,config)) {
00642          cpl_msg_error(_id,"cannot save tbl %s", cfg->paramsList);
00643      cpl_table_delete(qclog_tbl);
00644      cpl_table_delete(tbl_fp) ;
00645      cpl_free(col);
00646      destroy_2Darray ( wavelength_clean, lx );
00647      destroy_2Dintarray (row_clean, lx);
00648      destroy_intarray ( n_found_lines );
00649      destroy_intarray ( sum_pointer );
00650      destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00651      destroyFitParams(par);
00652      cpl_table_delete(tbl_line_list);
00653      destroy_image(im);
00654      destroy_image(map);
00655      cpl_frameset_delete(raw);
00656      wavecal_free(cfg);
00657      return -1;
00658       }
00659 
00660       cpl_table_delete(qclog_tbl);
00661        cpl_table_delete(tbl_fp) ;
00662       cpl_free(col);
00663 
00664    }
00665    /* free memory */
00666    destroy_2Darray ( wavelength_clean, lx );
00667    destroy_2Dintarray (row_clean, lx);
00668    destroy_intarray ( n_found_lines );
00669    destroy_intarray ( sum_pointer );
00670    destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00671 
00672    cpl_table_delete(tbl_line_list);
00673 
00674 
00675 /*----------------------------------------------------------------------
00676  *-------------------WAVEMAP--------------------------------------------
00677  *----------------------------------------------------------------------
00678  */
00679 
00680 /*
00681 #---now do the cross correlation and produce a wavelength map---
00682  */
00683 } else if (cfg->wavemapInd == 1 && cfg->calibIndicator == 0) { 
00684   cpl_msg_info(_id,"Wavemap");
00685   acoefs = new_2Dfloat_array ( cfg->nrDispCoefficients, lx);
00686    /* #read the parameterized dispersion relation */
00687 
00688    strcpy(tbl_name,cfg->coeffsName);
00689    tbl_wcal = cpl_table_load(tbl_name,1,0);
00690     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00691       cpl_msg_info(_id,"cannot load table %s",tbl_name);
00692       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00693       destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00694       cpl_table_delete(tbl_wcal);
00695       cpl_table_delete(tbl_line_list);
00696       destroy_image(im);
00697       cpl_frameset_delete(raw);
00698       wavecal_free(cfg);
00699       return -1;
00700     }
00701    for (i =0; i < lx; i++) {
00702       for (j = 0; j< cfg->nrDispCoefficients; j++) {
00703           sprintf(col_name,"%s%d","coeff",j);
00704             acoefs[j][i]=cpl_table_get_double(tbl_wcal,col_name,i,status);
00705       }
00706    }
00707     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00708       cpl_msg_info(_id,"cannot read table %s",tbl_name);
00709       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00710       destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00711       cpl_table_delete(tbl_wcal);
00712       cpl_table_delete(tbl_line_list);
00713       destroy_image(im);
00714       cpl_frameset_delete(raw);
00715       wavecal_free(cfg);
00716       return -1;
00717     }
00718    cpl_table_delete(tbl_wcal);
00719 
00720    map = createShiftedSlitWavemap2 ( im, 
00721                                           acoefs, 
00722                                           cfg->nrDispCoefficients,
00723                                           wave, 
00724                                           intens, 
00725                                           n_lines, 
00726                                           cfg->magFactor, 
00727                       cfg->guessDispersion1, 
00728                                           cfg->pixeldist );
00729    if (map == NULL) {
00730            cpl_msg_error(_id,"createShiftedSlitWavemap2 failed!\n");
00731            return -1;
00732    }
00733 
00734    par = newFitParams(15*n_lines);
00735    cpl_msg_info(_id,"Check shifts");
00736 
00737    shift = checkCorrelatedLinePositions ( im, acoefs, 
00738                                            cfg->nrDispCoefficients, 
00739                                            wave, 
00740                                            intens, 
00741                                            n_lines, 
00742                                            cfg->fwhm, 
00743                                            cfg->halfWidth, 
00744                                            cfg->minAmplitude, 
00745                                            cfg->guessDispersion1, 
00746                                            par );
00747 
00748 
00749    if (FLAG == shift){
00750       cpl_msg_error(_id,"checkCorrelatedLinePositions failed!\n");
00751    }
00752 
00753    map_imgw=cpl_image_wrap_float(map->lx,map->ly,map->data);
00754    map_img=cpl_image_duplicate(map_img);
00755    cpl_image_unwrap(map_imgw);
00756 
00757     det_ncounts(raw, cfg->qc_thresh_max,&qc);
00758     qclog_tbl = sinfoni_qclog_init(3);
00759 
00760     sprintf(key_value,"%d",n_lines);
00761     sinfoni_qclog_add(qclog_tbl,0,"QC NLINES","CPL_TYPE_INT",
00762               key_value,"Number of found lines");
00763 
00764 
00765 
00766     fwhm_avg = cpl_table_get_column_mean(tbl_fp,"fpar1");
00767     fwhm_med = cpl_table_get_column_median(tbl_fp,"fpar1");
00768 
00769     sprintf(key_value,"%f",fwhm_med);
00770 
00771     sinfoni_qclog_add(qclog_tbl,1,"QC FWHM MED","CPL_TYPE_DOUBLE",
00772               key_value,"Median FWHM of found lines");
00773 
00774     sprintf(key_value,"%f",fwhm_avg);
00775     sinfoni_qclog_add(qclog_tbl,2,"QC FWHM AVG","CPL_TYPE_DOUBLE",
00776               key_value,"Average FWHM of found lines");
00777 
00778 
00779     if(-1 == sinfoni_pro_save_ima(map_img,raw,sof,cfg->outName,
00780          PRO_WAVE_MAP,qclog_tbl,plugin_id,config)) {
00781          cpl_msg_error(_id,"cannot save ima %s", cfg->outName);
00782          destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00783      cpl_table_delete(qclog_tbl);
00784      cpl_image_delete(map_img);
00785          cpl_frameset_delete(raw);
00786          wavecal_free(cfg);
00787          return -1;
00788     }
00789 
00790     cpl_table_delete(qclog_tbl);
00791     cpl_image_delete(map_img);
00792 
00793 
00794    /* # ---free memory--- */
00795    /* destroy_image(map); */
00796    destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00797    /* To fix a memory bug we comment the following. But check! */
00798    /* cpl_free ( wave ); */
00799    /* cpl_free ( intens ); */
00800 
00801 } else if (cfg->wavemapInd == 1 && cfg->calibIndicator == 1) {
00802    cpl_msg_error(_id,"give either wavemapIndicator = yes and calibIndicator = no \
00803                 or wavemapIndicator = no and calibIndicator = yes") ;
00804 
00805    cpl_table_delete(tbl_line_list);
00806    destroy_image(im);
00807    cpl_frameset_delete(raw);
00808    wavecal_free(cfg);
00809    return -1;
00810 }
00811 
00812 /*-------------------------------------------------------------------
00813  *-------------------SLITFITS----------------------------------------
00814  *-------------------------------------------------------------------
00815  #--fit the slitlet edge positions if desired--
00816  */
00817 if (cfg->slitposIndicator == 1) {
00818   cpl_msg_info(_id,"fit the slitlet edge positions");
00819    if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00820 
00821    /* #read the first integer to determine the number of fit 
00822       parameters to allocate
00823     */
00824 
00825      if(is_fits_file(cfg->paramsList) !=1 ) {
00826        cpl_msg_error(_id,"cannot read FITS file %s ",cfg->paramsList);
00827        return -1; 
00828      }
00829          strcpy(tbl_fitpar_name,cfg->paramsList);
00830          tbl_fitpar = cpl_table_load(tbl_fitpar_name,1,0);
00831          if(cpl_error_get_code() != CPL_ERROR_NONE) {
00832 
00833             cpl_msg_error(_id, "error loading tbl %s ",
00834                  tbl_fitpar_name);
00835             cpl_msg_error(_id, (char* ) cpl_error_get_message());
00836             destroyFitParams(par);
00837             cpl_table_delete(tbl_fitpar);
00838             cpl_table_delete(tbl_line_list);
00839             destroy_image(im);
00840             destroy_image(map);
00841             cpl_frameset_delete(raw);
00842             wavecal_free(cfg);
00843             return -1;
00844      }
00845          readsum=cpl_table_get_int(tbl_fitpar,"n_params",1,status);
00846           if(cpl_error_get_code() != CPL_ERROR_NONE) {
00847 
00848             cpl_msg_error(_id, "error getting readsum from tbl %s ",
00849                  tbl_fitpar_name);
00850             cpl_msg_error(_id, (char* ) cpl_error_get_message());
00851             return -1;
00852      }
00853          
00854          cpl_table_delete(tbl_fitpar);
00855  
00856 
00857      /*
00858       if ( NULL == (par_list = fopen (cfg->paramsList, "r" ) ) )
00859       {
00860          cpl_msg_error(_id,"cannot open %s\n", cfg->paramsList) ;
00861          return -1 ;
00862       }
00863      */
00864       /* read 5 elements keep last */
00865       /* fscanf( par_list, "%i", &readsum); */
00866 
00867       /* readsum = string.atoi(par_list.read(5)); */
00868       /* fclose(par_list); */
00869       /* #allocate memory for the fit parameters and read 
00870          them from an ASCII file
00871        */
00872        par = newFitParams( readsum );
00873        if (par == NULL) {
00874           cpl_msg_error(_id,"newFitParams failed!\n");
00875           return -1;
00876        }
00877        /*
00878        dumpAsciiToFitParams(par, cfg->paramsList);
00879        */
00880       if( -1 == dumpTblToFitParams(par, cfg->paramsList) ) {
00881      cpl_msg_error(_id,"reading tbl %s ", cfg->paramsList);
00882          return -1;
00883        }
00884  
00885    }
00886 
00887    /* #allocate memory for the slitlet position array */
00888    slit_pos = new_2Dfloat_array(32,2);
00889    /* #now decide which fit model function you want to use by 
00890       reading a given character
00891     */
00892 
00893    cpl_msg_info(_id,"cfg->fitBoltzIndicator=%d\n",cfg->fitBoltzIndicator);
00894    cpl_msg_info(_id,"cfg->estimateIndicator=%d\n",cfg->estimateIndicator);
00895    cpl_msg_info(_id,"cfg->fitEdgeIndicator=%d\n",cfg->fitEdgeIndicator);
00896 
00897    if (cfg->fitBoltzIndicator == 1) {
00898 
00899       if (cfg->estimateIndicator == 1) {
00900       /* #open the ASCII list of the slitlet positions--- */
00901 
00902 
00903      /*READ TFITS TABLE*/
00904          strcpy(tbl_slitpos_guess_name,cfg->slitposGuessName);
00905          tbl_slitpos_guess = cpl_table_load(tbl_slitpos_guess_name,1,0);
00906          if(cpl_error_get_code() != CPL_ERROR_NONE) {
00907 
00908             cpl_msg_error(_id, "error loading tbl %s ",
00909                  tbl_slitpos_guess_name);
00910             cpl_msg_error(_id, (char* ) cpl_error_get_message());
00911             return -1;
00912      }
00913          n = cpl_table_get_nrow(tbl_slitpos_guess);
00914          if(cpl_error_get_code() != CPL_ERROR_NONE) {
00915 
00916             cpl_msg_error(_id, "error getting ncol from tbl %s ",
00917                  tbl_slitpos_guess_name);
00918             cpl_msg_error(_id, (char* ) cpl_error_get_message());
00919             return -1;
00920      }
00921 
00922          for (i =0 ; i< 32; i++){
00923            val_x=cpl_table_get_double(tbl_slitpos_guess,"pos1",i,status);
00924            val_y=cpl_table_get_double(tbl_slitpos_guess,"pos2",i,status);
00925        array2D_set_value(slit_pos,val_x,i,0);
00926        array2D_set_value(slit_pos,val_y,i,1);
00927      }
00928          if(cpl_error_get_code() != CPL_ERROR_NONE) {
00929 
00930             cpl_msg_error(_id, "error reading tbl %s ",
00931                  tbl_slitpos_guess_name);
00932             cpl_msg_error(_id, (char* ) cpl_error_get_message());
00933             return -1;
00934      }
00935          cpl_table_delete(tbl_slitpos_guess);
00936 
00937          cpl_msg_info(_id,"fitSlitsBoltzWithEstimate");
00938          fit = fitSlitsBoltzWithEstimate( im, slit_pos, 
00939                            cfg->boxLength, cfg->yBox, 
00940                                        cfg->diffTol, 
00941                                        cfg->loPos, 
00942                                        cfg->hiPos );
00943          if (fit < 0){
00944             cpl_msg_error(_id,"fitSlitsBoltzWithEstimate failed" );
00945         destroy_2Darray ( slit_pos, 32 );
00946         destroy_image ( im );
00947         destroy_image ( map );
00948         wavecal_free(cfg);
00949         cpl_frameset_delete(raw);
00950                 destroyFitParams(par);
00951             return -1;
00952      }
00953       } else {
00954               cpl_msg_info(_id,"fitSlitsBoltz");
00955               fit = fitSlitsBoltz( im, par, slit_pos, cfg->boxLength, 
00956                                  cfg->yBox, cfg->diffTol );
00957 
00958           if (fit < 0) {
00959                 cpl_msg_error (_id, "fitSlitsBoltz failed" );
00960         destroy_2Darray ( slit_pos, 32 );
00961         destroy_image ( im );
00962         destroy_image ( map );
00963         wavecal_free(cfg);
00964                 destroyFitParams(par);
00965         cpl_frameset_delete(raw);
00966                 return -1;
00967           }
00968 
00969       } 
00970 
00971    } else if (cfg->fitEdgeIndicator == 1) {
00972 
00973       if (cfg->estimateIndicator == 1){
00974      /*READ TFITS TABLE*/
00975          strcpy(tbl_slitpos_guess_name,cfg->slitposGuessName);
00976          tbl_slitpos_guess = cpl_table_load(tbl_slitpos_guess_name,1,0);
00977          n = cpl_table_get_nrow(tbl_slitpos_guess);
00978 
00979          for (i =0 ; i< 32; i++){
00980            val_x=cpl_table_get_double(tbl_slitpos_guess,"pos1",i,status);
00981            val_y=cpl_table_get_double(tbl_slitpos_guess,"pos2",i,status);
00982        array2D_set_value(slit_pos,val_x,i,0);
00983        array2D_set_value(slit_pos,val_y,i,1);
00984      }
00985          cpl_table_delete(tbl_slitpos_guess);
00986 
00987          cpl_msg_info(_id,"fitSlitsEdgeWithEstimate");
00988          fit = fitSlitsEdgeWithEstimate( im, slit_pos,
00989                                            cfg->boxLength, cfg->yBox, 
00990                                            cfg->diffTol, cfg->loPos, 
00991                         cfg->hiPos );
00992          if (fit < 0) {
00993                cpl_msg_error(_id, "fitSlitsBoltz failed" );
00994            destroy_2Darray ( slit_pos, 32 );
00995            destroy_image ( im );
00996            destroy_image ( map );
00997            wavecal_free(cfg);
00998            destroyFitParams(par);
00999            cpl_frameset_delete(raw);
01000                return -1;
01001      }
01002       } else {
01003          cpl_msg_info(_id,"fitSlitsEdge");
01004          fit = fitSlitsEdge( im, par, slit_pos, 
01005                                 cfg->boxLength, cfg->yBox, cfg->diffTol );
01006          if (fit < 0) {
01007                 cpl_msg_error(_id,"fitSlitsEdge failed" );
01008             destroy_2Darray ( slit_pos, 32 );
01009             destroy_image ( im );
01010             destroy_image ( map );
01011             wavecal_free(cfg);
01012                     destroyFitParams(par);
01013             cpl_frameset_delete(raw);
01014                 return -1;
01015      }
01016       }
01017    } else {
01018          cpl_msg_error(_id,"no indication of desired fit function given" );
01019          return -1;
01020    }
01021 
01022     /* #store the resulting sitlet positions in an TFITS table */
01023    tbl_spos = cpl_table_new(32);
01024    cpl_table_new_column(tbl_spos,"pos1", CPL_TYPE_DOUBLE);
01025    cpl_table_new_column(tbl_spos,"pos2", CPL_TYPE_DOUBLE);
01026    cpl_table_set_column_format(tbl_spos,"pos1", "15.9f");
01027    cpl_table_set_column_format(tbl_spos,"pos2", "15.9f");
01028 
01029     for (i =0; i< 32; i++) {
01030      /*
01031      fprintf( slitpos_file, "%15.9f %15.9f \n",
01032       array2D_get_value(slit_pos,i,0),array2D_get_value(slit_pos,i,1));
01033      */
01034       cpl_table_set_double(tbl_spos,"pos1",i,array2D_get_value(slit_pos,i,0));
01035       cpl_table_set_double(tbl_spos,"pos2",i,array2D_get_value(slit_pos,i,1));
01036      
01037    }
01038     if(sw == 1) {
01039        strcpy(tbl_name,"out_slitpos_guess.tfits");
01040 
01041        if(-1 == sinfoni_pro_dump_tbl(tbl_spos,raw,sof,tbl_name,
01042                      PRO_SLIT_POS_GUESS,plugin_id,config)) {
01043          cpl_msg_error(_id,"cannot dump tbl %s", tbl_name);
01044      cpl_table_delete(tbl_spos);
01045      destroy_2Darray ( slit_pos, 32 );
01046          if ( (cfg->slitposIndicator == 1 && cfg->estimateIndicator != 1) ||
01047               (cfg->calibIndicator == 1)  || (cfg->wavemapInd == 1) ){
01048               destroyFitParams(par);
01049      }
01050          destroy_image ( im );
01051          destroy_image ( map );
01052          wavecal_free(cfg);
01053          cpl_frameset_delete(raw);
01054          return -1;
01055        }
01056 
01057     } else {
01058        strcpy(tbl_name,cfg->slitposName);
01059 
01060        if(-1 == sinfoni_pro_dump_tbl(tbl_spos,raw,sof,tbl_name,
01061                      PRO_SLIT_POS,plugin_id,config)) {
01062          cpl_msg_error(_id,"cannot dump tbl %s", tbl_name);
01063      cpl_table_delete(tbl_spos);
01064      destroy_2Darray ( slit_pos, 32 );
01065          if ( (cfg->slitposIndicator == 1 && cfg->estimateIndicator != 1) ||
01066               (cfg->calibIndicator == 1)  || (cfg->wavemapInd == 1) ){
01067               destroyFitParams(par);
01068      }
01069          destroy_image ( im );
01070          destroy_image ( map );
01071          wavecal_free(cfg);
01072          cpl_frameset_delete(raw);
01073          return -1;
01074 
01075        }
01076 
01077 
01078     }
01079    
01080     cpl_table_delete(tbl_spos);
01081    /*# free memory*/
01082    destroy_2Darray ( slit_pos, 32 );
01083    
01084    
01085 }
01086 
01087 /* #-----free the rest memory--*/
01088 if ( (cfg->slitposIndicator == 1 && cfg->estimateIndicator != 1) || 
01089      (cfg->calibIndicator == 1)  || (cfg->wavemapInd == 1) ){
01090      destroyFitParams(par);
01091 }
01092 
01093 destroy_image ( im );
01094 destroy_image ( map );
01095 wavecal_free(cfg);
01096 cpl_frameset_delete(raw);
01097 
01098  return 0;
01099 
01100 
01101 
01102 }
01103 

Generated on Wed Oct 26 13:08:56 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001