00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "slit_pos.h"
00014 #include "sinfoni_pro_save.h"
00015 #include "wavecal_ini_by_cpl.h"
00016 #include "sinfoni_wcal_functions.h"
00017
00018
00019
00020
00021
00022
00023
00024 int dumpFitParamsToTbl ( FitParams ** params, char * filename );
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 int slit_pos (cpl_parameterlist* config, cpl_frameset* sof)
00073 {
00074 const char* _id = "slit_pos";
00075 wave_config * cfg =NULL;
00076 int check = 0;
00077 int lx = 0;
00078 int ly = 0;
00079 int n_lines=0;
00080 int i = 0;
00081 int j = 0;
00082 int n = 0;
00083
00084 int sum=0;
00085
00086 int* n_found_lines=NULL;
00087 int* sum_pointer=NULL;
00088 int** row_clean;
00089
00090 float a=0;
00091 float shift=0;
00092 float* wave=NULL;
00093 float* intens=NULL;
00094
00095 float** acoefs;
00096 float** wavelength_clean;
00097
00098 float** slit_pos;
00099
00100 OneImage * map=NULL ;
00101 OneImage * im=NULL ;
00102
00103 FitParams** par;
00104
00105 cpl_table* tbl_wcal=NULL;
00106 cpl_table* tbl_spos=NULL;
00107
00108 char* col_name=NULL;
00109 char* tbl_name=NULL;
00110
00111 char* tbl_line_list_name=NULL;
00112 cpl_table* tbl_line_list = NULL;
00113 int* status=NULL;
00114 cpl_image* map_img=NULL;
00115
00116 cpl_frameset* raw=NULL;
00117
00118
00119 cpl_table * tbl_fp =NULL;
00120 char* col=NULL;
00121 cpl_table* qclog_tbl=NULL;
00122 char* key_value=NULL;
00123 char* key_name=NULL;
00124 double fwhm_med=0;
00125 double fwhm_avg=0;
00126 double coef_med=0;
00127 double coef_avg=0;
00128 int trow=0;
00129 struct qc_wcal qc;
00130
00131
00132
00133
00134
00135
00136 cpl_msg_info(_id,"Parsing cpl input");
00137 cfg = parse_cpl_input_wave(config,sof,&raw) ;
00138
00139
00140 cfg->nslitlets=32;
00141 cfg->calibIndicator=1;
00142 cfg->wavemapInd=0;
00143 cfg->slitposIndicator=1;
00144
00145 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00146 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00147 return -1;
00148 }
00149
00150 if (cfg == NULL)
00151 {
00152 cpl_msg_error(_id,"could not parse cpl input!\n") ;
00153 return -1 ;
00154 }
00155 if(is_fits_file(cfg->inFrame) != 1) {
00156 cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrame);
00157 return -1;
00158 }
00159
00160
00161 if (cfg->slitposIndicator == 1 && cfg->estimateIndicator == 1) {
00162 if (is_fits_file(cfg->slitposGuessName) != 1) {
00163 cpl_msg_error(_id,"slitlet position guess list not given!");
00164 return -1;
00165 }
00166 }
00167
00168 if (cfg->calibIndicator == 0 && cfg->wavemapInd == 1) {
00169 if (is_fits_file(cfg->coeffsName) != 1) {
00170 cpl_msg_error(_id,"coefficients list not given!");
00171 return -1;
00172 }
00173 }
00174
00175
00176 if (cfg->slitposIndicator == 1) {
00177 if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00178
00179 if (is_fits_file(cfg->paramsList) != 1) {
00180 cpl_msg_error(_id,"parameter list not given!");
00181 return -1;
00182 }
00183
00184 }
00185 }
00186
00187
00188 im = load_image(cfg->inFrame);
00189 if (im == NULL) {
00190 cpl_msg_error(_id,"could not load image\n");
00191 return -1;
00192 }
00193
00194
00195 lx = im->lx;
00196 ly = im->ly;
00197
00198
00199
00200 if (cfg->calibIndicator == 1 || cfg->wavemapInd == 1) {
00201
00202
00203 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00204 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00205 return -1;
00206 }
00207
00208 tbl_line_list_name=cfg->lineList;
00209 tbl_line_list = cpl_table_load(tbl_line_list_name,1,0);
00210 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00211 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00212 return -1;
00213 }
00214 n = cpl_table_get_nrow(tbl_line_list);
00215 n_lines = n;
00216 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00217 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00218 return -1;
00219 }
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230 wave = cpl_table_get_data_float(tbl_line_list,"wave");
00231 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00232 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00233 return -1;
00234 }
00235
00236 intens = cpl_table_get_data_float(tbl_line_list,"int");
00237 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00238 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00239 return -1;
00240 }
00241
00242 }
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256 cpl_msg_info(_id,"guessBeginWave=%g",cfg->guessBeginWavelength);
00257 cpl_msg_info(_id,"guessDisp1=%g",cfg->guessDispersion1);
00258 cpl_msg_info(_id,"guessDisp2=%g",cfg->guessDispersion2);
00259
00260 if (cfg->calibIndicator == 1 && cfg->wavemapInd == 0) {
00261 cpl_msg_info(_id,"Findlines");
00262 acoefs = new_2Dfloat_array(cfg->nrDispCoefficients, lx);
00263
00264 n_found_lines = new_int_array(lx);
00265 row_clean = new_2Dint_array(lx, n_lines);
00266 wavelength_clean = new_2Dfloat_array(lx, n_lines);
00267 sum_pointer = new_int_array(1) ;
00268
00269 intarray_set_value(sum_pointer, 0, 0);
00270 check = findLines(im, wave, intens, n_lines, row_clean,
00271 wavelength_clean, cfg->guessBeginWavelength,
00272 cfg->guessDispersion1, cfg->guessDispersion2,
00273 cfg->mindiff, cfg->halfWidth,
00274 n_found_lines, cfg->sigma, sum_pointer );
00275 if (-1 == check) {
00276 cpl_msg_error(_id,"findLines failed!\n");
00277 return -1;
00278 }
00279
00280
00281
00282
00283
00284
00285
00286
00287 cpl_msg_info(_id,"Wave Calibration");
00288 sum = intarray_get_value(sum_pointer,0);
00289
00290 par = newFitParams( sum );
00291 if (par == NULL) {
00292 cpl_msg_error(_id,"newFitParams failed!\n");
00293 return -1;
00294 }
00295
00296
00297
00298
00299
00300 slit_pos = new_2Dfloat_array(32,2);
00301
00302 map = spred_waveCal(im,
00303 par,
00304 acoefs,
00305 cfg->nslitlets,
00306 row_clean,
00307 wavelength_clean,
00308 n_found_lines,
00309 cfg->guessDispersion1,
00310 cfg->halfWidth,
00311 cfg->minAmplitude,
00312 cfg->maxResidual,
00313 cfg->fwhm,
00314 cfg->nrDispCoefficients,
00315 cfg->nrCoefCoefficients,
00316 cfg->sigmaFactor,
00317 cfg->pixeldist,
00318 cfg->pixel_tolerance,
00319 slit_pos);
00320
00321
00322 if (map == NULL ) {
00323 cpl_msg_error(_id,"wave_cal failed!\n");
00324 return -1;
00325 }
00326 cpl_msg_info(_id,"Check line positions");
00327
00328 shift = checkLinePositions (im, acoefs, cfg->nrDispCoefficients, par);
00329 if (FLAG == shift){
00330 cpl_msg_error(_id,"checkForLinePositions failed!\n");
00331 }
00332
00333 map_img=cpl_image_wrap_float(map->lx,map->ly,map->data);
00334
00335 det_ncounts(raw, cfg->qc_thresh_max,&qc);
00336 qclog_tbl = sinfoni_qclog_init(3);
00337 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00338 sprintf(key_value,"%d",n_lines);
00339 sinfoni_qclog_add(qclog_tbl,0,"QC WAVE ALL","CPL_TYPE_INT",
00340 key_value,"Number of found lines");
00341
00342 sprintf(key_value,"%d",qc.nsat);
00343 sinfoni_qclog_add(qclog_tbl,1,"QC WAVE NPIXSAT","CPL_TYPE_INT",
00344 key_value,"Number of saturated pixels");
00345
00346 sprintf(key_value,"%g",qc.max_di);
00347 sinfoni_qclog_add(qclog_tbl,2,"QC WAVE MAXFLUX","CPL_TYPE_INT",
00348 key_value,"Max int off-lamp subracted frm");
00349
00350 if(-1 == sinfoni_pro_save_ima(map_img,raw,sof,cfg->outName,
00351 PRO_WAVE_MAP,qclog_tbl,_id,config)) {
00352 cpl_msg_error(_id,"cannot save ima %s", cfg->outName);
00353 }
00354 cpl_table_delete(qclog_tbl);
00355 cpl_free(key_value);
00356
00357
00358
00359
00360
00361
00362
00363
00364 if (cfg->writeCoeffsInd == 1) {
00365 col_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00366 tbl_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00367 tbl_wcal = cpl_table_new(lx);
00368 for (i=0; i< cfg->nrDispCoefficients; i++) {
00369 sprintf(col_name,"%s%d","coeff",i);
00370 cpl_table_new_column(tbl_wcal,col_name, CPL_TYPE_DOUBLE);
00371 }
00372
00373
00374 qclog_tbl = sinfoni_qclog_init(1+2*cfg->nrDispCoefficients);
00375 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00376 key_name = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00377 sprintf(key_value,"%d",n_lines);
00378 sinfoni_qclog_add(qclog_tbl,0,"QC WAVE ALL","CPL_TYPE_INT",
00379 key_value,"Number of found lines");
00380
00381
00382
00383 for (j=0; j< lx; j++) {
00384 for (i=0; i< cfg->nrDispCoefficients; i++) {
00385 sprintf(col_name,"%s%d","coeff",i);
00386 a = array2D_get_value(acoefs, i, j);
00387
00388 cpl_table_set_double(tbl_wcal,col_name,j,a);
00389 }
00390
00391
00392 }
00393
00394
00395 for (i=0; i< cfg->nrDispCoefficients; i++) {
00396 sprintf(col_name,"%s%d","coeff",i);
00397 coef_avg=cpl_table_get_column_mean(tbl_wcal,col_name);
00398 coef_med=cpl_table_get_column_median(tbl_wcal,col_name);
00399
00400 trow=1+i;
00401 sprintf(key_name,"%s%d%s","QC COEF",i," AVG");
00402 sprintf(key_value,"%g",coef_avg);
00403
00404 sinfoni_qclog_add(qclog_tbl,trow,key_name,"CPL_TYPE_DOUBLE",
00405 key_value,"Average wavecal Coef");
00406
00407 trow=1+i+cfg->nrDispCoefficients;
00408 sprintf(key_name,"%s%d%s","QC COEF",i," MED");
00409 sprintf(key_value,"%g",coef_med);
00410 sinfoni_qclog_add(qclog_tbl,trow,key_name,"CPL_TYPE_DOUBLE",
00411 key_value,"Median wavecal Coef");
00412
00413 }
00414
00415
00416
00417
00418 strcpy(tbl_name,cfg->coeffsName);
00419
00420 if(-1 == sinfoni_pro_save_tbl(tbl_wcal,raw,sof,tbl_name,
00421 PRO_WAVE_COEF_SLIT,qclog_tbl,_id,config)) {
00422 cpl_msg_error(_id,"cannot save tbl %s", tbl_name);
00423 }
00424 cpl_table_delete(tbl_wcal);
00425 cpl_table_delete(qclog_tbl);
00426 cpl_free(key_value);
00427 cpl_free(key_name);
00428
00429 cpl_free(col_name);
00430 cpl_free(tbl_name);
00431
00432 }
00433
00434
00435
00436
00437
00438 if (cfg->writeParInd == 1) {
00439
00440
00441 dumpFitParamsToAscii(par, cfg->paramsList);
00442
00443 if ( NULL == par )
00444 {
00445 cpl_msg_error (_id,"no fit parameters available!") ;
00446 return -1;
00447 }
00448
00449 if ( NULL == cfg->paramsList )
00450 {
00451 cpl_msg_error (_id,"no filename available!") ;
00452 return -1;
00453 }
00454
00455 tbl_fp = cpl_table_new(par[0] -> n_params);
00456 cpl_table_new_column(tbl_fp,"n_params", CPL_TYPE_INT);
00457 cpl_table_new_column(tbl_fp,"column", CPL_TYPE_INT);
00458 cpl_table_new_column(tbl_fp,"line", CPL_TYPE_INT);
00459 col = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00460
00461 for(j=0;j<4;j++) {
00462 sprintf(col,"%s%d","fpar",j);
00463 cpl_table_new_column(tbl_fp,col, CPL_TYPE_DOUBLE);
00464 sprintf(col,"%s%d","dpar",j);
00465 cpl_table_new_column(tbl_fp,col, CPL_TYPE_DOUBLE);
00466 }
00467
00468
00469
00470 qclog_tbl = sinfoni_qclog_init(3);
00471 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00472 sprintf(key_value,"%d",n_lines);
00473 sinfoni_qclog_add(qclog_tbl,0,"QC NLINES","CPL_TYPE_INT",
00474 key_value,"Number of Found lines");
00475
00476 for ( i = 0 ; i < par[0] -> n_params ; i++ )
00477 {
00478 cpl_table_set_int(tbl_fp,"n_params",i,par[i]->n_params);
00479 cpl_table_set_int(tbl_fp,"column",i,par[i]->column);
00480 cpl_table_set_int(tbl_fp,"line",i,par[i]->line);
00481
00482
00483 for(j=0;j<4;j++) {
00484 sprintf(col,"%s%d","fpar",j);
00485 cpl_table_set_double(tbl_fp,col,i,par[i]->fit_par[j]);
00486 sprintf(col,"%s%d","dpar",j);
00487 cpl_table_set_double(tbl_fp,col,i,par[i]->derv_par[j]);
00488 }
00489 }
00490
00491 fwhm_avg = cpl_table_get_column_mean(tbl_fp,"fpar1");
00492 fwhm_med = cpl_table_get_column_median(tbl_fp,"fpar1");
00493 sprintf(key_value,"%f",fwhm_med);
00494
00495 sinfoni_qclog_add(qclog_tbl,1,"QC FWHM MED","CPL_TYPE_DOUBLE",
00496 key_value,"Median FWHM of found lines");
00497
00498 sprintf(key_value,"%f",fwhm_avg);
00499 sinfoni_qclog_add(qclog_tbl,2,"QC FWHM AVG","CPL_TYPE_DOUBLE",
00500 key_value,"Average FWHM of found lines");
00501
00502
00503 if(-1 == sinfoni_pro_save_tbl(tbl_fp,raw,sof,cfg->paramsList,
00504 PRO_WAVE_PAR_LIST,qclog_tbl,_id,config)) {
00505 cpl_msg_error(_id,"cannot save tbl %s", cfg->paramsList);
00506 }
00507 cpl_table_delete(qclog_tbl);
00508 cpl_free(key_value);
00509
00510 cpl_table_delete(tbl_fp) ;
00511 cpl_free(col);
00512
00513 }
00514
00515 destroy_2Darray ( wavelength_clean, lx );
00516 destroy_2Dintarray (row_clean, lx);
00517 destroy_intarray ( n_found_lines );
00518 destroy_intarray ( sum_pointer );
00519 destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529 } else if (cfg->wavemapInd == 1 && cfg->calibIndicator == 0) {
00530 cpl_msg_info(_id,"Wavemap");
00531 acoefs = new_2Dfloat_array ( cfg->nrDispCoefficients, lx);
00532
00533
00534 tbl_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00535 col_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00536 strcpy(tbl_name,cfg->coeffsName);
00537 tbl_wcal = cpl_table_load(tbl_name,1,0);
00538 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00539 cpl_msg_info(_id,"cannot load table %s",tbl_name);
00540 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00541 return -1;
00542 }
00543 for (i =0; i < lx; i++) {
00544 for (j = 0; j< cfg->nrDispCoefficients; j++) {
00545 sprintf(col_name,"%s%d","coeff",j);
00546 acoefs[j][i]=cpl_table_get_double(tbl_wcal,col_name,i,status);
00547 }
00548 }
00549 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00550 cpl_msg_info(_id,"cannot read table %s",tbl_name);
00551 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00552 return -1;
00553 }
00554 cpl_free(col_name);
00555 cpl_free(tbl_name);
00556 cpl_table_delete(tbl_wcal);
00557
00558 map = createShiftedSlitWavemap2 ( im,
00559 acoefs,
00560 cfg->nrDispCoefficients,
00561 wave,
00562 intens,
00563 n_lines,
00564 cfg->magFactor,
00565 cfg->guessDispersion1,
00566 cfg->pixeldist );
00567 if (map == NULL) {
00568 cpl_msg_error(_id,"createShiftedSlitWavemap2 failed!\n");
00569 return -1;
00570 }
00571
00572 par = newFitParams(15*n_lines);
00573 cpl_msg_info(_id,"Check shifts");
00574
00575 shift = checkCorrelatedLinePositions ( im, acoefs,
00576 cfg->nrDispCoefficients,
00577 wave,
00578 intens,
00579 n_lines,
00580 cfg->fwhm,
00581 cfg->halfWidth,
00582 cfg->minAmplitude,
00583 cfg->guessDispersion1,
00584 par );
00585
00586
00587 if (FLAG == shift){
00588 cpl_msg_error(_id,"checkCorrelatedLinePositions failed!\n");
00589 }
00590
00591 map_img=cpl_image_wrap_float(map->lx,map->ly,map->data);
00592
00593 det_ncounts(raw, cfg->qc_thresh_max,&qc);
00594 qclog_tbl = sinfoni_qclog_init(3);
00595
00596 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00597 sprintf(key_value,"%d",n_lines);
00598 sinfoni_qclog_add(qclog_tbl,0,"QC NLINES","CPL_TYPE_INT",
00599 key_value,"Number of found lines");
00600
00601
00602
00603 fwhm_avg = cpl_table_get_column_mean(tbl_fp,"fpar1");
00604 fwhm_med = cpl_table_get_column_median(tbl_fp,"fpar1");
00605
00606 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00607 sprintf(key_value,"%f",fwhm_med);
00608
00609 sinfoni_qclog_add(qclog_tbl,1,"QC FWHM MED","CPL_TYPE_DOUBLE",
00610 key_value,"Median FWHM of found lines");
00611
00612 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00613 sprintf(key_value,"%f",fwhm_avg);
00614 sinfoni_qclog_add(qclog_tbl,2,"QC FWHM AVG","CPL_TYPE_DOUBLE",
00615 key_value,"Average FWHM of found lines");
00616
00617
00618 if(-1 == sinfoni_pro_save_ima(map_img,raw,sof,cfg->outName,
00619 PRO_WAVE_MAP,qclog_tbl,_id,config)) {
00620 cpl_msg_error(_id,"cannot save ima %s", cfg->outName);
00621 }
00622 cpl_table_delete(qclog_tbl);
00623 cpl_free(key_value);
00624
00625
00626
00627
00628 destroy_2Darray ( acoefs, cfg->nrDispCoefficients );
00629
00630
00631
00632
00633 } else if (cfg->wavemapInd == 1 && cfg->calibIndicator == 1) {
00634 cpl_msg_error(_id,"give either wavemapIndicator = yes and calibIndicator = no \
00635 or wavemapIndicator = no and calibIndicator = yes") ;
00636 return -1;
00637 }
00638
00639
00640
00641
00642
00643
00644
00645 if (cfg->slitposIndicator == 1) {
00646 cpl_msg_info(_id,"fit the slitlet edge positions");
00647
00648
00649
00650 tbl_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00651 tbl_spos = cpl_table_new(32);
00652 cpl_table_new_column(tbl_spos,"pos1", CPL_TYPE_DOUBLE);
00653 cpl_table_new_column(tbl_spos,"pos2", CPL_TYPE_DOUBLE);
00654 cpl_table_set_column_format(tbl_spos,"pos1", "15.9f");
00655 cpl_table_set_column_format(tbl_spos,"pos2", "15.9f");
00656
00657 for (i =0; i< 32; i++) {
00658
00659
00660
00661
00662 cpl_table_set_double(tbl_spos,"pos1",i,array2D_get_value(slit_pos,i,0));
00663 cpl_table_set_double(tbl_spos,"pos2",i,array2D_get_value(slit_pos,i,1));
00664
00665 }
00666
00667 strcpy(tbl_name,"out_guess_slit_pos.tfits");
00668 if(-1 == sinfoni_pro_dump_tbl(tbl_spos,raw,sof,tbl_name,
00669 PRO_SLIT_POS,_id,config)) {
00670 cpl_msg_error(_id,"cannot dump tbl %s", tbl_name);
00671 }
00672 cpl_table_delete(tbl_spos);
00673 cpl_free(tbl_name);
00674
00675 destroy_2Darray ( slit_pos, 32 );
00676
00677
00678 }
00679
00680
00681 if ( (cfg->slitposIndicator == 1 && cfg->estimateIndicator != 1) ||
00682 (cfg->calibIndicator == 1) || (cfg->wavemapInd == 1) ){
00683 destroyFitParams(par);
00684 }
00685 destroy_image ( im );
00686 destroy_image ( map );
00687 wavecal_free(cfg);
00688
00689 return 0;
00690
00691
00692
00693 }
00694
00695
00696