00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "cube_create_from_resampled.h"
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 void change_plist_cube_create_from_resampled (cpl_propertylist * plist,
00033 char * ini_name,
00034 char * name,
00035 char * outName,
00036 float cenLambda,
00037 float dispersion,
00038 int cenpix,
00039 float center_x,
00040 float center_y )
00041 {
00042
00043 float pixelscale ;
00044 float ra ;
00045 float dec ;
00046 float angle ;
00047 float radangle ;
00048 float cd1_1, cd1_2, cd2_1, cd2_2 ;
00049 char firsttext[2*FILENAMESZ] ;
00050 char * date ;
00051
00052 strcpy(firsttext, "spiffi_cube_create_from_resampled -f \0") ;
00053
00054 pixelscale = spiffi_get_pixelscale(name) ;
00055 ra = spiffi_get_RA(name) ;
00056 dec = spiffi_get_DEC(name) ;
00057 angle = spiffi_get_posangle(name) ;
00058 radangle = angle * PI_NUMB / 180. ;
00059 cd1_1 = cos(radangle) * pixelscale / 3600. ;
00060 cd1_2 = -sin(radangle) * pixelscale / 3600. ;
00061 cd2_1 = sin(radangle) * pixelscale / 3600. ;
00062 cd2_2 = cos(radangle) * pixelscale / 3600. ;
00063
00064
00065
00066 cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "RA---TAN" ) ;
00067 cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00068
00069 cpl_propertylist_insert_after_float(plist, "CTYPE1", "CRPIX1", center_x ) ;
00070 cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00071
00072
00073
00074 cpl_propertylist_insert_after_float(plist, "CRPIX1", "CRVAL1", ra ) ;
00075 cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00076
00077 cpl_propertylist_insert_after_float(plist, "CRVAL1", "CDELT1", -pixelscale/3600. ) ;
00078 cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00079
00080 cpl_propertylist_insert_after_string(plist, "CDELT1", "CUNIT1", "DEGREE" ) ;
00081 cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00082
00083 cpl_propertylist_insert_after_string(plist, "CUNIT1", "CTYPE2", "DEC--TAN" ) ;
00084 cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00085
00086
00087
00088
00089 cpl_propertylist_insert_after_float(plist, "CTYPE2", "CRPIX2", center_y ) ;
00090 cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00091
00092 cpl_propertylist_insert_after_float(plist, "CRPIX2", "CRVAL2", dec) ;
00093 cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00094
00095 cpl_propertylist_insert_after_float(plist, "CRVAL2", "CDELT2", pixelscale/3600. ) ;
00096 cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00097
00098 cpl_propertylist_insert_after_string(plist, "CDELT2", "CUNIT2", "DEGREE" ) ;
00099 cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00100
00101 cpl_propertylist_insert_after_string(plist, "CUNIT2", "CTYPE3", "WAVE" ) ;
00102 cpl_propertylist_set_comment(plist, "CTYPE3", "wavelength axis in microns") ;
00103
00104
00105 cpl_propertylist_insert_after_int(plist, "CTYPE3", "CRPIX3", cenpix+1 ) ;
00106 cpl_propertylist_set_comment(plist, "CRPIX3", "Reference pixel in z") ;
00107
00108 cpl_propertylist_insert_after_float(plist, "CRPIX3", "CRVAL3", cenLambda ) ;
00109 cpl_propertylist_set_comment(plist, "CRVAL3", "central wavelength") ;
00110
00111 cpl_propertylist_insert_after_float(plist, "CRVAL3", "CDELT3", dispersion ) ;
00112 cpl_propertylist_set_comment(plist, "CDELT3", "microns per pixel") ;
00113
00114 cpl_propertylist_insert_after_string(plist, "CDELT3", "CUNIT3", "MICRON" ) ;
00115 cpl_propertylist_set_comment(plist, "CUNIT3", "spectral unit" ) ;
00116
00117
00118
00119 cpl_propertylist_insert_after_float(plist, "CUNIT3", "CD1_1", cd1_1 ) ;
00120 cpl_propertylist_set_comment(plist, "CD1_1", "CD rotation matrix" ) ;
00121
00122 cpl_propertylist_insert_after_float(plist, "CD1_1", "CD1_2", cd1_2 ) ;
00123 cpl_propertylist_set_comment(plist, "CD1_2", "CD rotation matrix" ) ;
00124
00125 cpl_propertylist_insert_after_float(plist, "CD1_2", "CD2_1", cd2_1 ) ;
00126 cpl_propertylist_set_comment(plist, "CD2_1", "CD rotation matrix" ) ;
00127
00128 cpl_propertylist_insert_after_float(plist, "CD2_1", "CD2_2", cd2_2 ) ;
00129 cpl_propertylist_set_comment(plist, "CD2_2", "CD rotation matrix" ) ;
00130
00131 date = get_datetime_iso8601() ;
00132 strcat (firsttext , ini_name) ;
00133
00134
00135 cpl_propertylist_set_string(plist,KEY_NAME_PIPEFILE, outName ) ;
00136 cpl_propertylist_set_string(plist,KEY_NAME_HPRO_TYPE, "REDUCED" ) ;
00137 cpl_propertylist_set_string(plist,KEY_NAME_HPRO_CATG, "Science Observation" ) ;
00138 cpl_propertylist_set_string(plist,KEY_NAME_HPRO_STATUS, "OK" ) ;
00139 cpl_propertylist_set_string(plist,KEY_NAME_HPRO_DATE, date ) ;
00140 cpl_propertylist_set_string(plist,KEY_NAME_HPRO_RECID, firsttext ) ;
00141 cpl_propertylist_set_string(plist,KEY_NAME_HPRO_DRSID, PACKAGE_VERSION ) ;
00142
00143
00144 }
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 int cube_create_from_resampled (cpl_parameterlist* config, cpl_frameset* sof)
00163 {
00164
00165 const char* _id = "cube_create_from_resampled";
00166
00167 cube_config * cfg ;
00168 OneImage * im ;
00169 OneCube * cube ;
00170 OneCube * outcube ;
00171 fits_header * head ;
00172 float fcol ;
00173 FILE * slitlist ;
00174 FILE * list ;
00175 FILE * fc_list ;
00176 float ** slit_edges ;
00177 float * correct_dist ;
00178 float * distances ;
00179 int i ;
00180 cpl_frameset* raw=NULL;
00181
00182
00183 cfg = parse_cpl_input_cube(config,sof,&raw) ;
00184
00185
00186 if (cfg == NULL)
00187 {
00188 cpl_msg_error(_id," could not parse .ini file!") ;
00189 return -1 ;
00190 }
00191
00192 if(is_fits_file(cfg->inFrame) != 1 ){
00193 cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrame);
00194 return -1;
00195 }
00196
00197 if (cfg->halocorrectInd == 1)
00198 {
00199 if(is_fits_file(cfg->halospectrum) != 1) {
00200 cpl_msg_error(_id,"Input file %s is not FITS",cfg->halospectrum);
00201 return -1;
00202 }
00203 }
00204 if (is_fits_file(cfg->poslist) != 1)
00205 {
00206 cpl_msg_error(_id,"Input file %s is not FITS",cfg->poslist) ;
00207 return -1 ;
00208 }
00209
00210 im = load_image(cfg->inFrame) ;
00211 if (im == NULL)
00212 {
00213 cpl_msg_error(_id, " could not load inFrame" ) ;
00214 return -1 ;
00215 }
00216 head = fits_read_header(cfg->inFrame) ;
00217
00218
00219
00220
00221
00222
00223 if (NULL == (correct_dist = (float*) cpl_calloc(cfg->nslits, sizeof (float))) )
00224 {
00225 cpl_msg_error(_id," could not allocate memory!") ;
00226 return -1 ;
00227 }
00228
00229 if (cfg->northsouthInd == 0)
00230 {
00231 if (NULL == (slit_edges = (float**) cpl_calloc (cfg->nslits, sizeof(float*)) ) )
00232 {
00233 cpl_msg_error(_id," could not allocate memory!") ;
00234 return -1 ;
00235 }
00236 for ( i = 0 ; i < cfg->nslits ; i++ )
00237 {
00238 if (NULL == (slit_edges[i] = (float*) cpl_calloc (2, sizeof(float))))
00239 {
00240 cpl_msg_error(_id," could not allocate memory!") ;
00241 return -1 ;
00242 }
00243 }
00244
00245 if ( NULL == (slitlist = fopen(cfg->poslist, "r" )) )
00246 {
00247 cpl_msg_error(_id," could not open poslist!") ;
00248 return -1 ;
00249 }
00250
00251 i = 0 ;
00252 while (EOF != fscanf(slitlist, "%f %f\n", &slit_edges[i][0] , &slit_edges[i][1]) )
00253 {
00254 i++ ;
00255 }
00256 if ( i != cfg->nslits)
00257 {
00258 cpl_msg_error(_id," wrong file format!") ;
00259 fclose(slitlist) ;
00260 return -1 ;
00261 }
00262 fclose(slitlist) ;
00263
00264
00265 cube = makeCubeSpi( im, slit_edges, correct_dist ) ;
00266 if (cube == NULL)
00267 {
00268 cpl_msg_error(_id," could not construct data cube!") ;
00269 return -1 ;
00270 }
00271
00272 for ( i = 0 ; i < cfg->nslits ; i++ )
00273 {
00274 cpl_free (slit_edges[i]) ;
00275 }
00276 cpl_free (slit_edges) ;
00277 }
00278 else
00279 {
00280 if ( NULL == (distances = (float*) cpl_calloc (cfg->nslits - 1, sizeof (float))) )
00281 {
00282 cpl_msg_error(_id," could allocate memory!") ;
00283 return -1 ;
00284 }
00285
00286 if ( NULL == (list = fopen(cfg->poslist, "r" ) ) )
00287 {
00288 cpl_msg_error(_id," could not open poslist!") ;
00289 return -1 ;
00290 }
00291
00292 i = 0 ;
00293 while (EOF != fscanf(list, "%f\n", &distances[i] ))
00294 {
00295 i++ ;
00296 }
00297
00298 if ( i != cfg->nslits - 1)
00299 {
00300 cpl_msg_error(_id," wrong file format!") ;
00301 fclose(list) ;
00302 return -1 ;
00303 }
00304
00305 fclose(list) ;
00306
00307 if ( NULL == (fc_list = fopen(cfg->firstCol, "r" ) ) )
00308 {
00309 cpl_msg_error(_id," could not open firstCol list!") ;
00310 return -1 ;
00311 }
00312 fscanf(fc_list, "%f\n", &fcol) ;
00313 fclose(fc_list) ;
00314
00315
00316 cube = makeCubeDist( im, fcol, distances, correct_dist ) ;
00317 if (cube == NULL)
00318 {
00319 cpl_msg_error(_id," could not construct a data cube!") ;
00320 return -1 ;
00321 }
00322 cpl_free(distances) ;
00323 }
00324
00325
00326
00327
00328
00329
00330
00331
00332 if (cfg->method == 'P')
00333 {
00334 outcube = fineTuneCube( cube, correct_dist, cfg->order ) ;
00335 if (outcube == NULL)
00336 {
00337 cpl_msg_error(_id," could not fine tune the data cube") ;
00338 return -1 ;
00339 }
00340 }
00341 else if (cfg->method == 'F')
00342 {
00343 for ( i = 0 ; i < cfg->nslits ; i++ )
00344 {
00345 correct_dist[i] = -correct_dist[i] ;
00346 }
00347 outcube = fineTuneCubeByFFT( cube, correct_dist ) ;
00348 if ( outcube == NULL )
00349 {
00350 cpl_msg_error(_id," could not fine tune the data cube") ;
00351 return -1 ;
00352 }
00353 }
00354 else if (cfg->method == 'S')
00355 {
00356 outcube = fineTuneCubeBySpline( cube, correct_dist ) ;
00357 if ( outcube == NULL )
00358 {
00359 cpl_msg_error(_id," could not fine tune the data cube") ;
00360 return -1 ;
00361 }
00362 }
00363 else
00364 {
00365 cpl_msg_error(_id," wrong method indicator given!") ;
00366 return -1 ;
00367 }
00368
00369 save_cube_to_fits_hdr_dump( outcube, cfg -> outName, head ) ;
00370
00371
00372 fits_header_destroy(head) ;
00373 cpl_free(correct_dist) ;
00374 destroy_image(im) ;
00375 destroy_cube (cube) ;
00376 destroy_cube (outcube) ;
00377 cube_cfg_destroy (cfg) ;
00378 return 0 ;
00379 }
00380
00381