00001 #include "utilities_scired.h"
00002
00025
00026 OneCube * cube_getvig(
00027 OneCube * cube_in,
00028 int loleft_x,
00029 int loleft_y,
00030 int upright_x,
00031 int upright_y)
00032 {
00033 const char * fctid="cube_getvig";
00034 OneCube * cube_out ;
00035 int i ;
00036 int outlx,
00037 outly ;
00038
00039 if (cube_in==NULL) return NULL ;
00040
00041 if ((loleft_x>upright_x) ||
00042 (loleft_y>upright_y)) {
00043 cpl_msg_error(fctid,"ill-defined slit for extraction: aborting");
00044 return NULL ;
00045 }
00046
00047
00048 outlx = upright_x - loleft_x + 1 ;
00049 outly = upright_y - loleft_y + 1 ;
00050
00051 cube_out = new_cube(outlx, outly, cube_in->np) ;
00052
00053 for (i=0 ; i<cube_in->np ; i++) {
00054
00055
00056
00057
00058
00059 cube_out->plane[i] =
00060 image_getvig(cube_in->plane[i],
00061 loleft_x, loleft_y,
00062 upright_x, upright_y) ;
00063 }
00064 return cube_out ;
00065 }
00066
00083
00084
00085 OneImage * image_getvig(
00086 OneImage * image_in,
00087 int loleft_x,
00088 int loleft_y,
00089 int upright_x,
00090 int upright_y)
00091 {
00092 OneImage * slit_img ;
00093 int i, j ;
00094 register
00095 pixelvalue * inpt,
00096 * outpt ;
00097 int outlx, outly ;
00098 const char* fctid="image_getvig";
00099 if (image_in==NULL) return NULL ;
00100
00101 if ((loleft_x<1) || (loleft_x>image_in->lx) ||
00102 (loleft_y<1) || (loleft_y>image_in->ly) ||
00103 (upright_x<1) || (upright_x>image_in->lx) ||
00104 (upright_y<1) || (upright_y>image_in->ly) ||
00105 (loleft_x>upright_x) || (loleft_y>upright_y)) {
00106 cpl_msg_error(fctid,"extraction zone is [%d %d] [%d %d]\n"
00107 "cannot extract such zone: aborting slit extraction",
00108 loleft_x, loleft_y, upright_x, upright_y) ;
00109 return NULL ;
00110 }
00111
00112 outlx = upright_x - loleft_x + 1 ;
00113 outly = upright_y - loleft_y + 1 ;
00114 slit_img = new_image(outlx, outly) ;
00115
00116 for (j=0 ; j<outly ; j++) {
00117 inpt = image_in->data+loleft_x-1 + (j+loleft_y-1)*image_in->lx ;
00118 outpt = slit_img->data + j*outlx ;
00119 for (i=0 ; i<outlx ; i++) {
00120 *outpt++ = *inpt++ ;
00121 }
00122 }
00123 return slit_img ;
00124 }
00125
00126
00127 int
00128 set_wcs_cube(cpl_imagelist* cub, char* name, double clambda,
00129 double dis, double cpix, double cx, double cy)
00130 {
00131 const char* _id="set_wcs_cube";
00132 cpl_propertylist* plist=NULL;
00133
00134 if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00135 cpl_msg_error(_id, "getting header from frame %s",name);
00136 cpl_propertylist_delete(plist) ;
00137 return -1 ;
00138 }
00139
00140 change_plist_cube(plist, name, clambda, dis, cpix, cx, cy) ;
00141
00142 if (cpl_imagelist_save(cub, name, CPL_BPP_DEFAULT,plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00143 cpl_msg_error(_id, "Cannot save the product %s",name);
00144 cpl_propertylist_delete(plist) ;
00145 return -1 ;
00146 }
00147 cpl_propertylist_delete(plist) ;
00148 return 0;
00149
00150 }
00151
00152 int
00153 set_wcs_image(cpl_image* img, char* name, double cx, double cy)
00154 {
00155 const char* _id="set_wcs_image";
00156 cpl_propertylist* plist=NULL;
00157
00158 if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00159 cpl_msg_error(_id, "getting header from frame %s",name);
00160 cpl_propertylist_delete(plist) ;
00161 return -1 ;
00162 }
00163 change_plist_image(plist, name, cx, cy) ;
00164
00165 if (cpl_image_save(img, name, CPL_BPP_DEFAULT,plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00166 cpl_msg_error(_id, "Cannot save the product %s",name);
00167 cpl_propertylist_delete(plist) ;
00168 return -1 ;
00169 }
00170 cpl_propertylist_delete(plist) ;
00171 return 0;
00172
00173 }
00174
00175
00176
00177 int
00178 set_wcs_spectrum(cpl_image* img, char* name, double clambda,
00179 double dis, double cpix)
00180 {
00181 const char* _id="set_wcs_image";
00182 cpl_propertylist* plist=NULL;
00183
00184 if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00185 cpl_msg_error(_id, "getting header from frame %s",name);
00186 cpl_propertylist_delete(plist) ;
00187 return -1 ;
00188 }
00189 change_plist_spectrum(plist, clambda, dis,cpix) ;
00190
00191 if (cpl_image_save(img, name, CPL_BPP_DEFAULT,plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00192 cpl_msg_error(_id, "Cannot save the product %s",name);
00193 cpl_propertylist_delete(plist) ;
00194 return -1 ;
00195 }
00196 cpl_propertylist_delete(plist) ;
00197 return 0;
00198
00199 }
00200
00201 void change_plist_cube (cpl_propertylist * plist,
00202 char * name,
00203 float cenLambda,
00204 float dispersion,
00205 int cenpix,
00206 float center_x,
00207 float center_y )
00208 {
00209
00210 float pixelscale ;
00211 float ra ;
00212 float dec ;
00213 float angle ;
00214 float radangle ;
00215 float cd1_1, cd1_2, cd2_1, cd2_2 ;
00216 char firsttext[2*FILENAMESZ] ;
00217
00218 strcpy(firsttext, "spiffi_objnod -f \0") ;
00219
00220 pixelscale = spiffi_get_pixelscale(name)/2. ;
00221 ra = spiffi_get_RA(name) ;
00222 dec = spiffi_get_DEC(name) ;
00223 angle = spiffi_get_posangle(name) ;
00224 radangle = angle * PI_NUMB / 180. ;
00225 cd1_1 = cos(radangle) * pixelscale / 3600. ;
00226 cd1_2 = sin(radangle) * pixelscale / 3600. ;
00227 cd2_1 = -sin(radangle) * pixelscale / 3600. ;
00228 cd2_2 = cos(radangle) * pixelscale / 3600. ;
00229 cpl_propertylist_erase_regexp(plist, "^CTYPE1");
00230 cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "RA---TAN" ) ;
00231 cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00232
00233 cpl_propertylist_erase_regexp(plist, "^CRPIX1");
00234 cpl_propertylist_insert_after_float(plist, "CTYPE1", "CRPIX1", center_x ) ;
00235 cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00236
00237 cpl_propertylist_erase_regexp(plist, "^CRVAL1");
00238 cpl_propertylist_insert_after_float(plist, "CRPIX1", "CRVAL1", ra ) ;
00239 cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00240
00241 cpl_propertylist_erase_regexp(plist, "^CDELT1");
00242 cpl_propertylist_insert_after_float(plist, "CRVAL1", "CDELT1", -pixelscale/3600. ) ;
00243 cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00244
00245 cpl_propertylist_erase_regexp(plist, "^CUNIT1");
00246 cpl_propertylist_insert_after_string(plist, "CDELT1", "CUNIT1", "DEGREE" ) ;
00247 cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00248
00249 cpl_propertylist_erase_regexp(plist, "^CTYPE2");
00250 cpl_propertylist_insert_after_string(plist, "CUNIT1", "CTYPE2", "DEC--TAN" ) ;
00251 cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00252
00253 cpl_propertylist_erase_regexp(plist, "^CRPIX2");
00254 cpl_propertylist_insert_after_float(plist, "CTYPE2", "CRPIX2", center_y ) ;
00255 cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00256
00257 cpl_propertylist_erase_regexp(plist, "^CRVAL2");
00258 cpl_propertylist_insert_after_float(plist, "CRPIX2", "CRVAL2", dec) ;
00259 cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00260
00261 cpl_propertylist_erase_regexp(plist, "^CDELT2");
00262 cpl_propertylist_insert_after_float(plist, "CRVAL2", "CDELT2", pixelscale/3600. ) ;
00263 cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00264
00265 cpl_propertylist_erase_regexp(plist, "^CUNIT2");
00266 cpl_propertylist_insert_after_string(plist, "CDELT2", "CUNIT2", "DEGREE" ) ;
00267 cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00268
00269 cpl_propertylist_erase_regexp(plist, "^CTYPE3");
00270 cpl_propertylist_insert_after_string(plist, "CUNIT2", "CTYPE3", "WAVE" ) ;
00271 cpl_propertylist_set_comment(plist, "CTYPE3", "wavelength axis in microns") ;
00272
00273 cpl_propertylist_erase_regexp(plist, "^CRPIX3");
00274 cpl_propertylist_insert_after_int(plist, "CTYPE3", "CRPIX3", cenpix+1 ) ;
00275 cpl_propertylist_set_comment(plist, "CRPIX3", "Reference pixel in z") ;
00276
00277 cpl_propertylist_erase_regexp(plist, "^CRVAL3");
00278 cpl_propertylist_insert_after_float(plist, "CRPIX3", "CRVAL3", cenLambda ) ;
00279 cpl_propertylist_set_comment(plist, "CRVAL3", "central wavelength") ;
00280
00281 cpl_propertylist_erase_regexp(plist, "^CDELT3");
00282 cpl_propertylist_insert_after_float(plist, "CRVAL3", "CDELT3", dispersion ) ;
00283 cpl_propertylist_set_comment(plist, "CDELT3", "microns per pixel") ;
00284
00285 cpl_propertylist_erase_regexp(plist, "^CUNIT3");
00286 cpl_propertylist_insert_after_string(plist, "CDELT3", "CUNIT3", "MICRON" ) ;
00287 cpl_propertylist_set_comment(plist, "CUNIT3", "spectral unit" ) ;
00288
00289 cpl_propertylist_erase_regexp(plist, "^CD1_1");
00290 cpl_propertylist_insert_after_float(plist, "CUNIT3", "CD1_1", cd1_1 ) ;
00291 cpl_propertylist_set_comment(plist, "CD1_1", "CD rotation matrix" ) ;
00292
00293 cpl_propertylist_erase_regexp(plist, "^CD1_2");
00294 cpl_propertylist_insert_after_float(plist, "CD1_1", "CD1_2", cd1_2 ) ;
00295 cpl_propertylist_set_comment(plist, "CD1_2", "CD rotation matrix" ) ;
00296
00297 cpl_propertylist_erase_regexp(plist, "^CD2_1");
00298 cpl_propertylist_insert_after_float(plist, "CD1_2", "CD2_1", cd2_1 ) ;
00299 cpl_propertylist_set_comment(plist, "CD2_1", "CD rotation matrix" ) ;
00300
00301 cpl_propertylist_erase_regexp(plist, "^CD2_2");
00302 cpl_propertylist_insert_after_float(plist, "CD2_1", "CD2_2", cd2_2 ) ;
00303 cpl_propertylist_set_comment(plist, "CD2_2", "CD rotation matrix" ) ;
00304
00305
00306 }
00307
00308
00309 void change_plist_spectrum (cpl_propertylist * plist,
00310 double cenLambda,
00311 double dispersion,
00312 int cenpix)
00313 {
00314
00315 cpl_propertylist_erase_regexp(plist, "^CTYPE3");
00316 cpl_propertylist_erase_regexp(plist, "^CRPIX3");
00317 cpl_propertylist_erase_regexp(plist, "^CRVAL3");
00318 cpl_propertylist_erase_regexp(plist, "^CDELT3");
00319 cpl_propertylist_erase_regexp(plist, "^CUNIT3");
00320
00321 cpl_propertylist_erase_regexp(plist, "^CTYPE2");
00322 cpl_propertylist_erase_regexp(plist, "^CRPIX2");
00323 cpl_propertylist_erase_regexp(plist, "^CRVAL2");
00324 cpl_propertylist_erase_regexp(plist, "^CDELT2");
00325 cpl_propertylist_erase_regexp(plist, "^CUNIT2");
00326
00327
00328 cpl_propertylist_erase_regexp(plist, "^CD1_1");
00329 cpl_propertylist_erase_regexp(plist, "^CD1_2");
00330 cpl_propertylist_erase_regexp(plist, "^CD2_1");
00331 cpl_propertylist_erase_regexp(plist, "^CD2_2");
00332
00333
00334 cpl_propertylist_erase_regexp(plist, "^CTYPE1");
00335 cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "PIXEL");
00336 cpl_propertylist_set_comment(plist, "CTYPE1", "Pixel coordinate system.");
00337
00338 cpl_propertylist_erase_regexp(plist, "^CRPIX1");
00339 cpl_propertylist_insert_after_int(plist, "CTYPE1", "CRPIX1", 1 ) ;
00340 cpl_propertylist_set_comment(plist, "CRPIX1", "Reference pixel in x") ;
00341
00342 cpl_propertylist_erase_regexp(plist, "^CRVAL1");
00343 cpl_propertylist_insert_after_double(plist, "CRPIX1", "CRVAL1", 1 ) ;
00344 cpl_propertylist_set_comment(plist, "CRVAL1", "value of ref pixel.") ;
00345
00346 cpl_propertylist_erase_regexp(plist, "^CDELT1");
00347 cpl_propertylist_insert_after_double(plist, "CRVAL1", "CDELT1", 1 ) ;
00348 cpl_propertylist_set_comment(plist, "CDELT1", "pixel scale") ;
00349
00350
00351
00352
00353
00354 cpl_propertylist_erase_regexp(plist, "^CUNIT1");
00355 cpl_propertylist_insert_after_string(plist, "CDELT1", "CUNIT1", "Pixel" );
00356 cpl_propertylist_set_comment(plist, "CUNIT1", "spectral unit" );
00357
00358
00359 cpl_propertylist_erase_regexp(plist, "^CTYPE2");
00360 cpl_propertylist_insert_after_string(plist, "EXPTIME","CTYPE2","WAVE" );
00361 cpl_propertylist_set_comment(plist, "CTYPE2","wavelength axis in microns");
00362
00363
00364 cpl_propertylist_erase_regexp(plist, "^CRPIX2");
00365 cpl_propertylist_insert_after_int(plist, "CTYPE2", "CRPIX2",cenpix+1 ) ;
00366 cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in x") ;
00367
00368 cpl_propertylist_erase_regexp(plist, "^CRVAL2");
00369 cpl_propertylist_insert_after_double(plist, "CRPIX2","CRVAL2",cenLambda ) ;
00370 cpl_propertylist_set_comment(plist, "CRVAL2", "central wavelength") ;
00371
00372 cpl_propertylist_erase_regexp(plist, "^CDELT2");
00373 cpl_propertylist_insert_after_double(plist, "CRVAL2", "CDELT2",dispersion);
00374 cpl_propertylist_set_comment(plist, "CDELT2", "microns per pixel");
00375
00376 cpl_propertylist_erase_regexp(plist, "^CUNIT2");
00377 cpl_propertylist_insert_after_string(plist, "CDELT2", "CUNIT2", "MICRON");
00378 cpl_propertylist_set_comment(plist, "CUNIT2", "spectral unit" );
00379
00380
00381
00382 }
00383
00384
00385 void change_plist_image (cpl_propertylist * plist,
00386 char * name,
00387 float center_x,
00388 float center_y )
00389 {
00390
00391 float pixelscale ;
00392 float ra ;
00393 float dec ;
00394 float angle ;
00395 float radangle ;
00396 float cd1_1, cd1_2, cd2_1, cd2_2 ;
00397 char firsttext[2*FILENAMESZ] ;
00398
00399 strcpy(firsttext, "spiffi_objnod -f \0") ;
00400
00401 pixelscale = spiffi_get_pixelscale(name)/2. ;
00402 ra = spiffi_get_RA(name) ;
00403 dec = spiffi_get_DEC(name) ;
00404 angle = spiffi_get_posangle(name) ;
00405 radangle = angle * PI_NUMB / 180. ;
00406 cd1_1 = cos(radangle) * pixelscale / 3600. ;
00407 cd1_2 = sin(radangle) * pixelscale / 3600. ;
00408 cd2_1 = -sin(radangle) * pixelscale / 3600. ;
00409 cd2_2 = cos(radangle) * pixelscale / 3600. ;
00410
00411 cpl_propertylist_erase_regexp(plist, "^CTYPE1");
00412 cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "RA---TAN" ) ;
00413 cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00414
00415 cpl_propertylist_erase_regexp(plist, "^CRPIX1");
00416 cpl_propertylist_insert_after_float(plist, "CTYPE1", "CRPIX1", center_x ) ;
00417 cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00418
00419 cpl_propertylist_erase_regexp(plist, "^CRVAL1");
00420 cpl_propertylist_insert_after_float(plist, "CRPIX1", "CRVAL1", ra ) ;
00421 cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00422
00423 cpl_propertylist_erase_regexp(plist, "^CDELT1");
00424 cpl_propertylist_insert_after_float(plist, "CRVAL1", "CDELT1", -pixelscale/3600. ) ;
00425 cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00426
00427 cpl_propertylist_erase_regexp(plist, "^CUNIT1");
00428 cpl_propertylist_insert_after_string(plist, "CDELT1", "CUNIT1", "DEGREE" ) ;
00429 cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00430
00431 cpl_propertylist_erase_regexp(plist, "^CTYPE2");
00432 cpl_propertylist_insert_after_string(plist, "CUNIT1", "CTYPE2", "DEC--TAN" ) ;
00433 cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00434
00435 cpl_propertylist_erase_regexp(plist, "^CRPIX2");
00436 cpl_propertylist_insert_after_float(plist, "CTYPE2", "CRPIX2", center_y ) ;
00437 cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00438
00439 cpl_propertylist_erase_regexp(plist, "^CRVAL2");
00440 cpl_propertylist_insert_after_float(plist, "CRPIX2", "CRVAL2", dec) ;
00441 cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00442
00443 cpl_propertylist_erase_regexp(plist, "^CDELT2");
00444 cpl_propertylist_insert_after_float(plist, "CRVAL2", "CDELT2", pixelscale/3600. ) ;
00445 cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00446
00447 cpl_propertylist_erase_regexp(plist, "^CUNIT2");
00448 cpl_propertylist_insert_after_string(plist, "CDELT2", "CUNIT2", "DEGREE" ) ;
00449 cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00450
00451 cpl_propertylist_erase_regexp(plist, "^CD1_1");
00452 cpl_propertylist_insert_after_float(plist, "CUNIT2", "CD1_1", cd1_1 ) ;
00453 cpl_propertylist_set_comment(plist, "CD1_1", "CD rotation matrix" ) ;
00454
00455 cpl_propertylist_erase_regexp(plist, "^CD1_2");
00456 cpl_propertylist_insert_after_float(plist, "CD1_1", "CD1_2", cd1_2 ) ;
00457 cpl_propertylist_set_comment(plist, "CD1_2", "CD rotation matrix" ) ;
00458
00459 cpl_propertylist_erase_regexp(plist, "^CD2_1");
00460 cpl_propertylist_insert_after_float(plist, "CD1_2", "CD2_1", cd2_1 ) ;
00461 cpl_propertylist_set_comment(plist, "CD2_1", "CD rotation matrix" ) ;
00462
00463 cpl_propertylist_erase_regexp(plist, "^CD2_2");
00464 cpl_propertylist_insert_after_float(plist, "CD2_1", "CD2_2", cd2_2 ) ;
00465 cpl_propertylist_set_comment(plist, "CD2_2", "CD rotation matrix" ) ;
00466
00467
00468 }
00469
00470
00471 cpl_imagelist*
00472 sinfoni_cube2imglist(OneCube* cube)
00473 {
00474 cpl_imagelist* ims=NULL;
00475 int i=0;
00476 OneImage* eimg=NULL;
00477 cpl_image* cimg=NULL;
00478 cpl_image* wimg=NULL;
00479 ims=cpl_imagelist_new(cube->lx,cube->ly,cube->np,CPL_TYPE_FLOAT);
00480
00481 for(i=0;i<cube->np;i++) {
00482 eimg=cube->plane[i];
00483 wimg=cpl_image_wrap_float(eimg->lx,eimg->ly,eimg->data);
00484 cimg=cpl_image_duplicate(wimg);
00485 cpl_image_unwrap(wimg);
00486 cpl_imagelist_set(ims,cimg,i);
00487 }
00488
00489 return ims;
00490 }
00491
00492 OneCube** sinfoni_correct_median(OneCube** cubes, const int n_cubes)
00493 {
00494 const char* fctid="sinfoni_correct_median";
00495 int i=0;
00496 OneCube** cubes_cor=NULL;
00497 double median=0;
00498 pixelvalue pv=0;
00499 int z=0;
00500 if ( cubes == NULL ) {
00501 cpl_msg_error (fctid,"no cube list given!") ;
00502 return NULL ;
00503 }
00504 if ( n_cubes <= 0 ) {
00505 cpl_msg_error (fctid,"wrong number of data cubes in list!") ;
00506 return NULL ;
00507 }
00508
00509 cubes_cor = (OneCube**) cpl_calloc (n_cubes, sizeof (OneCube*)) ;
00510
00511 for ( i = 0 ; i < n_cubes ; i++ ) {
00512 cubes_cor[i] = newCube(cubes[i]->lx,cubes[i]->ly,cubes[i]->np);
00513 for(z=0;z< cubes[i]->np; z++) {
00514 pv=get_median_pixelvalue_image(cubes[i]->plane[z]);
00515 median=(double)pv;
00516 if(!qfits_isnan(median)) {
00517 cubes_cor[i]->plane[z]=cst_op_image(cubes[i]->plane[z],median,'-');
00518 }
00519 }
00520 }
00521
00522 return cubes_cor;
00523 }
00524
00525 int sinfoni_correct_median_it(OneCube** inp)
00526 {
00527
00528 const char* fctid="sinfoni_correct_median";
00529 double median=0;
00530 int z=0;
00531
00532 for(z=0;z< (*inp)->np; z++) {
00533 median=median_image((*inp)->plane[z]);
00534 if(!qfits_isnan(median)) {
00535 cst_op_image_local((*inp)->plane[z],median,'-');
00536 } else {
00537 cpl_msg_error(fctid,"median is NAN");
00538 }
00539 }
00540
00541 return 0;
00542 }
00543
00544
00545
00546 OneCube** sinfoni_correct_sky(OneCube** cubes,
00547 const int nc,
00548 OneCube* sky_cube)
00549
00550 {
00551 const char* fctid="sinfoni_correct_sky";
00552 OneCube** cubes_sky=NULL;
00553 int x=0;
00554 int y=0;
00555 int z=0;
00556 int i=0;
00557 float k=0.5;
00558 int ovr=0;
00559 int ks=0;
00560 int nclip=0;
00561 double med=0;
00562 double avg=0;
00563 double sig=0;
00564 int msk_sum=0;
00565 double val_msk_sum=0;
00566 cpl_vector* val=NULL;
00567 cpl_vector* msk=NULL;
00568
00569
00570 if ( cubes == NULL ) {
00571 cpl_msg_error (fctid,"no cube list given!") ;
00572 return NULL ;
00573 }
00574
00575 if ( nc <= 0 ) {
00576 cpl_msg_error (fctid,"wrong number of data cubes in list!") ;
00577 return NULL ;
00578 }
00579
00580 cubes_sky = (OneCube**) cpl_calloc (nc, sizeof (OneCube*)) ;
00581 for(z=0;z< cubes[0]->np; z++) {
00582 for(y=0;y< cubes[0]->ly; y++) {
00583 for(x=0;x< cubes[0]->lx; x++) {
00584
00585 msk=cpl_vector_new(nc);
00586 for (i=0;i<nc;i++) {
00587 cpl_vector_set(msk,i,1);
00588 }
00589 nclip=0;
00590 for (ks=0;ks<nc;ks++) {
00591 sig=0;
00592 med=0;
00593 ovr=0;
00594 val=cpl_vector_new(nc-nclip);
00595
00596 for ( i = 0 ; i < nc ; i++ ) {
00597 if ((!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) &&
00598 (cpl_vector_get(msk,i) != 0) ) {
00599 cpl_vector_set(val,ovr,(double)cubes[i]->plane[z]->data[x+y*cubes[i]->lx]);
00600 ovr++;
00601 }
00602 }
00603
00604 if(ovr>0) {
00605 avg=cpl_vector_get_mean(val);
00606 med=cpl_vector_get_median(val);
00607 if(ovr>1) {
00608 sig=cpl_vector_get_stdev(val);
00609 } else {
00610 sig=0;
00611 }
00612 } else {
00613 avg=cpl_vector_get(val,0);
00614 med=avg;
00615 sig=0;
00616 }
00617
00618 cpl_vector_delete(val);
00619 for ( i = 0 ; i < nc ; i++ ) {
00620
00621 if ((!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) &&
00622 (cpl_vector_get(msk,i) != 0)) {
00623 if(abs((cubes[i]->plane[z]->data[x+y*cubes[i]->lx]-med))> k*sig) {
00624 cpl_vector_set(msk,i,0);
00625 nclip++;
00626 }
00627 }
00628 }
00629 }
00630 msk_sum=0;
00631 val_msk_sum=0;
00632 for ( i = 0 ; i < nc ; i++ ) {
00633
00634 if (!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) {
00635 msk_sum+=cpl_vector_get(msk,i);
00636 val_msk_sum+=cubes[i]->plane[z]->data[x+y*cubes[i]->lx]*
00637 cpl_vector_get(msk,i);
00638 }
00639 }
00640 sky_cube->plane[z]->data[x+y*sky_cube->lx]=val_msk_sum/msk_sum;
00641 cpl_vector_delete(msk);
00642 }
00643 }
00644 }
00645 for ( i = 0 ; i < nc ; i++ ) {
00646 cubes_sky[i]=copy_cube(cubes[i]);
00647
00648 op_cube(&(cubes_sky[i]),sky_cube,'-');
00649
00650 }
00651
00652
00653 return cubes_sky;
00654 }
00655
00656
00657 OneCube** sinfoni_correct_sky2(OneCube** cubes,
00658 const int nc,
00659 OneCube* sky_cube,
00660 OneCube* med_cube,
00661 OneCube* msk_cube,
00662 OneCube* avg_cube,
00663 OneCube* sig_cube,
00664 OneCube* ovr_cube)
00665 {
00666 const char* fctid="sinfoni_correct_sky";
00667 OneCube** cubes_sky=NULL;
00668 int x=0;
00669 int y=0;
00670 int z=0;
00671 int i=0;
00672 float k=0.5;
00673 int ovr=0;
00674 int ks=0;
00675 int nclip=0;
00676 double med=0;
00677 double avg=0;
00678 double sig=0;
00679 int msk_sum=0;
00680 double val_msk_sum=0;
00681 cpl_vector* val=NULL;
00682 cpl_vector* msk=NULL;
00683
00684
00685 if ( cubes == NULL ) {
00686 cpl_msg_error (fctid,"no cube list given!") ;
00687 return NULL ;
00688 }
00689
00690 if ( nc <= 0 ) {
00691 cpl_msg_error (fctid,"wrong number of data cubes in list!") ;
00692 return NULL ;
00693 }
00694
00695 cubes_sky = (OneCube**) cpl_calloc (nc, sizeof (OneCube*)) ;
00696 for(z=0;z< cubes[0]->np; z++) {
00697 for(y=0;y< cubes[0]->ly; y++) {
00698 for(x=0;x< cubes[0]->lx; x++) {
00699
00700 msk=cpl_vector_new(nc);
00701 for (i=0;i<nc;i++) {
00702 cpl_vector_set(msk,i,1);
00703 }
00704 ovr_cube->plane[z]->data[x+y*ovr_cube->lx]=nc;
00705 msk_cube->plane[z]->data[x+y*msk_cube->lx]=nc;
00706 nclip=0;
00707 for (ks=0;ks<nc;ks++) {
00708 sig=0;
00709 med=0;
00710 ovr=0;
00711 val=cpl_vector_new(nc-nclip);
00712
00713 for ( i = 0 ; i < nc ; i++ ) {
00714 if ((!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) &&
00715 (cpl_vector_get(msk,i) != 0) ) {
00716 cpl_vector_set(val,ovr,(double)cubes[i]->plane[z]->data[x+y*cubes[i]->lx]);
00717 ovr++;
00718 }
00719 }
00720
00721 if(ovr>1) {
00722 avg=cpl_vector_get_mean(val);
00723 med=cpl_vector_get_median(val);
00724 sig=cpl_vector_get_stdev(val);
00725 } else {
00726 avg=cpl_vector_get(val,0);
00727 med=avg;
00728 sig=0;
00729 }
00730
00731 med_cube->plane[z]->data[x+y*med_cube->lx]=med;
00732 avg_cube->plane[z]->data[x+y*avg_cube->lx]=avg;
00733 sig_cube->plane[z]->data[x+y*sig_cube->lx]=sig;
00734 cpl_vector_delete(val);
00735 for ( i = 0 ; i < nc ; i++ ) {
00736
00737 if ((!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) &&
00738 (cpl_vector_get(msk,i) != 0)) {
00739 if(abs((cubes[i]->plane[z]->data[x+y*cubes[i]->lx]-med))> k*sig) {
00740
00741 msk_cube->plane[z]->data[x+y*msk_cube->lx]-=1;
00742 cpl_vector_set(msk,i,0);
00743 nclip++;
00744 }
00745 }
00746 }
00747 }
00748 msk_sum=0;
00749 val_msk_sum=0;
00750 for ( i = 0 ; i < nc ; i++ ) {
00751
00752 if (!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) {
00753
00754
00755
00756
00757
00758 msk_sum+=cpl_vector_get(msk,i);
00759 val_msk_sum+=cubes[i]->plane[z]->data[x+y*cubes[i]->lx]*
00760 cpl_vector_get(msk,i);
00761 }
00762 }
00763 sky_cube->plane[z]->data[x+y*sky_cube->lx]=val_msk_sum/msk_sum;
00764 cpl_vector_delete(msk);
00765 }
00766 }
00767 }
00768 for ( i = 0 ; i < nc ; i++ ) {
00769 cubes_sky[i]=copy_cube(cubes[i]);
00770
00771 op_cube(&(cubes_sky[i]),sky_cube,'-');
00772 }
00773
00774
00775 return cubes_sky;
00776 }
00777
00778 int
00779 extract_ao_info(cpl_frameset* obj_set,cpl_frameset** set,const char* recid,cpl_parameterlist* config)
00780 {
00781
00782 const char* fctid=NULL;
00783 cpl_frame* frm=NULL;
00784 cpl_table* tbl=NULL;
00785 cpl_frame* product_frame=NULL;
00786 cpl_propertylist* plist=NULL;
00787 cpl_image* image=NULL;
00788
00789 int nobj=0;
00790 int i=0;
00791 const char* name_o="out_ao_info.tfits";
00792
00793 name_o="out_ao_info.tfits";
00794 nobj=cpl_frameset_get_size(obj_set);
00795 for(i=0;i<nobj;i++) {
00796 frm=cpl_frameset_get_frame(obj_set,i);
00797 tbl=cpl_table_load(cpl_frame_get_filename(frm),1,1);
00798 plist=cpl_propertylist_load(cpl_frame_get_filename(frm),0);
00799 if(i==0) {
00800 cpl_table_save(tbl,NULL,plist,name_o,CPL_IO_DEFAULT);
00801 } else {
00802 cpl_table_save(tbl,NULL,plist,name_o,CPL_IO_EXTEND);
00803 }
00804 }
00805 image=cpl_image_new(0,0,CPL_TYPE_FLOAT);
00806
00807 product_frame = cpl_frame_new();
00808 cpl_frame_set_filename(product_frame, name_o) ;
00809 cpl_frame_set_tag(product_frame, PRO_AO_INFO) ;
00810 cpl_frame_set_type(product_frame, CPL_FRAME_TYPE_IMAGE) ;
00811 cpl_frame_set_group(product_frame, CPL_FRAME_GROUP_PRODUCT) ;
00812 cpl_frame_set_level(product_frame, CPL_FRAME_LEVEL_FINAL) ;
00813 if (cpl_error_get_code()) {
00814 cpl_msg_error(fctid, "Error while initialising the product frame") ;
00815 cpl_propertylist_delete(plist) ;
00816 cpl_frame_delete(product_frame) ;
00817 cpl_image_delete(image) ;
00818 return -1 ;
00819 }
00820
00821
00822 if (cpl_dfs_setup_product_header(plist, product_frame, *set, config,
00823 recid, "SINFONI", "?Dictionary?") != CPL_ERROR_NONE) {
00824 cpl_msg_error(fctid, "Problem in the product DFS-compliance") ;
00825 cpl_propertylist_delete(plist) ;
00826 cpl_frame_delete(product_frame) ;
00827 cpl_image_delete(image) ;
00828 return -1 ;
00829 }
00830
00831
00832 if (cpl_image_save(image, name_o, CPL_BPP_DEFAULT, plist,
00833 CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
00834 cpl_msg_error(fctid, "Could not save product");
00835 cpl_propertylist_delete(plist) ;
00836 cpl_frame_delete(product_frame) ;
00837 cpl_image_delete(image) ;
00838 return -1 ;
00839 }
00840 cpl_propertylist_delete(plist) ;
00841 cpl_image_delete(image) ;
00842
00843
00844 cpl_frameset_insert(*set, product_frame) ;
00845
00846 return 0;
00847
00848 }
00849
00850 int
00851 assign_offset(const int n,char* name,float* offsetx,float* offsety,
00852 const float ref_offx,const float ref_offy)
00853 {
00854
00855 const char * _id = "assign_offset";
00856 float offx=0;
00857 float offy=0;
00858 double mjd_obs=0;
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871 cpl_propertylist * plist=NULL;
00872 cpl_msg_info(_id,"Assign offsets");
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912 offx = spiffi_get_cumoffsetx(name) - ref_offx ;
00913 if (offx == FLAG) {
00914 cpl_msg_warning(_id," could not read fits header keyword cummoffsetx!") ;
00915 cpl_msg_warning(_id," Set relative offset to 0 - %f!",ref_offx) ;
00916 offx = - ref_offx;
00917
00918 }
00919
00920 offy = spiffi_get_cumoffsety(name) - ref_offy ;
00921 if (offy == FLAG ) {
00922 cpl_msg_warning(_id," could not read fits header keyword! cumoffsety") ;
00923 cpl_msg_warning(_id," Set relative offset to 0 - %f!",ref_offx) ;
00924 offy = - ref_offy;
00925
00926 }
00927 cpl_msg_info(_id,"offx=%f offy=%f",offx,offy);
00928
00929
00930
00931
00932 if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00933 cpl_msg_error(_id, "getting header from reference frame %s",name);
00934 cpl_propertylist_delete(plist) ;
00935 cpl_free(name);
00936 return -1 ;
00937 }
00938
00939 if (cpl_propertylist_contains(plist, KEY_NAME_MJD_OBS)) {
00940 mjd_obs=cpl_propertylist_get_double(plist, KEY_NAME_MJD_OBS);
00941 } else {
00942 cpl_msg_error(_id,"keyword %s does not exist",KEY_NAME_MJD_OBS);
00943 cpl_propertylist_delete(plist) ;
00944 return -1;
00945 }
00946 cpl_propertylist_delete(plist) ;
00947
00948 if (mjd_obs > 53421.58210082 ) {
00949
00950
00951
00952
00953
00954 array_set_value(offsetx,-2*offx,n);
00955 array_set_value(offsety,2*offy,n);
00956 } else {
00957
00958
00959
00960
00961
00962 array_set_value(offsetx,2*offx,n);
00963 array_set_value(offsety,-2*offy,n);
00964 }
00965
00966 return 0;
00967
00968
00969 }
00970
00971 int
00972 object_assign_offset(char* name, const int n, double* ref_offx, double* ref_offy, float** offsetx, float** offsety)
00973 {
00974
00975 const char * _id = "object_assign_offset";
00976 float offx=0;
00977 float offy=0;
00978 double mjd_obs=0;
00979 cpl_propertylist * plist=NULL;
00980 cpl_msg_info(_id,"Assign offsets");
00981 if ( n == 0 )
00982 {
00983
00984
00985 *ref_offx = spiffi_get_cumoffsetx(name) ;
00986 if (*ref_offx == FLAG)
00987 {
00988 cpl_msg_error(_id," could not read fits header keyword cummoffsetx!") ;
00989
00990 }
00991
00992 *ref_offy = spiffi_get_cumoffsety(name) ;
00993 if (*ref_offy == FLAG )
00994 {
00995 cpl_msg_error(_id," could not read fits header keyword! cumoffsety") ;
00996
00997 }
00998 cpl_msg_info(_id,"Reference offx=%f offy=%f",*ref_offx,*ref_offy);
00999
01000 offx = 0. ;
01001 offy = 0. ;
01002
01003 } else {
01004
01005 offx = spiffi_get_cumoffsetx(name) - *ref_offx ;
01006 if (offx == FLAG)
01007 {
01008 cpl_msg_error(_id," could not read fits header keyword cummoffsetx!") ;
01009
01010 }
01011
01012 offy = spiffi_get_cumoffsety(name) - *ref_offy ;
01013 if (offy == FLAG )
01014 {
01015 cpl_msg_error(_id," could not read fits header keyword! cumoffsety") ;
01016
01017 }
01018 cpl_msg_info(_id,"offx=%f offy=%f",offx,offy);
01019
01020 }
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032 if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
01033 cpl_msg_error(_id, "getting header from reference frame %s",name);
01034 cpl_propertylist_delete(plist) ;
01035 cpl_free(name);
01036 return -1 ;
01037 }
01038
01039 if (cpl_propertylist_contains(plist, KEY_NAME_MJD_OBS)) {
01040 mjd_obs=cpl_propertylist_get_double(plist, KEY_NAME_MJD_OBS);
01041 } else {
01042 cpl_msg_error(_id,"keyword %s does not exist",KEY_NAME_MJD_OBS);
01043 cpl_propertylist_delete(plist) ;
01044 return -1;
01045 }
01046 cpl_propertylist_delete(plist) ;
01047
01048 if (mjd_obs > 53421.58210082 ) {
01049
01050 array_set_value(*offsetx,-offx*2,n);
01051 array_set_value(*offsety,+offy*2,n);
01052 } else {
01053
01054 array_set_value(*offsetx,+offx*2,n);
01055 array_set_value(*offsety,-offy*2,n);
01056 }
01057
01058 return 0;
01059 }
01060
01061
01062
01063
01064 OneCube*
01065 fine_tune(OneCube* cube,float* correct_dist, const char* method, const int order, const int nslits) {
01066 const char * _id= "fine_tune";
01067 int i =0;
01068 OneCube* outcube2=NULL;
01069 float* neg_dist=NULL;
01070 cpl_msg_info(_id,"Finetuning, method=%s",method);
01071
01072 if (strcmp(method,"P")==0)
01073 {
01074 outcube2 = fineTuneCube( cube, correct_dist, order ) ;
01075 if (outcube2 == NULL)
01076 {
01077 cpl_msg_error (_id," could not fine tune the data cube\n") ;
01078 return NULL ;
01079 }
01080 }
01081 else if (strcmp(method,"F")==0)
01082 {
01083 neg_dist=cpl_calloc(nslits,sizeof(float));
01084 for ( i = 0 ; i < nslits ; i++ )
01085 {
01086 neg_dist[i] = -correct_dist[i] ;
01087 }
01088 outcube2 = fineTuneCubeByFFT( cube, neg_dist ) ;
01089 cpl_free(neg_dist);
01090 if ( outcube2 == NULL )
01091 {
01092 cpl_msg_error (_id," could not fine tune the data cube\n") ;
01093 return NULL ;
01094 }
01095 }
01096 else if (strcmp(method,"S")==0)
01097 {
01098 outcube2 = fineTuneCubeBySpline( cube, correct_dist ) ;
01099 if ( outcube2 == NULL )
01100 {
01101 cpl_msg_error (_id," could not fine tune the data cube\n") ;
01102 return NULL ;
01103 }
01104 }
01105 else
01106 {
01107 cpl_msg_error (_id," wrong method indicator given!") ;
01108 return NULL ;
01109 }
01110
01111
01112
01113 return outcube2;
01114
01115 }