/*---------------------------------------------------------------------------- File name : lw_spjitter.c Author : Y. Jung Created on : March 2001 Description : ISAAC spectro jitter recipe in LW mode ---------------------------------------------------------------------------*/ /* $Id: lw_spjitter.c,v 1.34 2002/03/04 14:03:40 yjung Exp $ $Author: yjung $ $Date: 2002/03/04 14:03:40 $ $Revision: 1.34 $ */ /*---------------------------------------------------------------------------- Includes ---------------------------------------------------------------------------*/ #include "eclipse.h" #include "isaacp_lib.h" /*---------------------------------------------------------------------------- Defines ---------------------------------------------------------------------------*/ #define DISXMIN 1 #define DISYMIN 50 #define DISXMAX 1024 #define DISYMAX 975 #define WC_DISHIGH 50 #define WC_DISLOW 50 #define SPECTRACT_BADLEFT 50 #define SPECTRACT_BADRIGHT 50 #define SPECTRACT_BADTOP 0 #define SPECTRACT_BADBOT 0 /*--------------------------------------------------------------------------- New types ---------------------------------------------------------------------------*/ typedef struct _lw_spjitter_BB_ { /* Frames */ char frames_filelist[FILENAMESZ]; /* Input ascii file */ int nframes ; /* Nb of frames */ char ** framelist ; /* List of input images */ int nb_hc_free ; char ** framelist_hc_free ; int nb_clean_frames ; char ** framelist_clean ; int nb_hc1 ; /* Nb of half cycle frames */ char ** hc1_list ; int nb_hc2 ; char ** hc2_list ; int spectrum_length ; /* Calibration frames */ int cal_arc_active ; /* Arc table available */ int cal_startrace_active ; /* Startrace table available */ int cal_spflat_active ; /* Flatfield available */ char cal_arc_name[FILENAMESZ]; char cal_startrace_name[FILENAMESZ]; char cal_spflat_name[FILENAMESZ]; /* Classification */ int divided_by_flat ; /* Flag for flatfield division */ int offsets_source ; /* OFFSOURCE_HEADER or OFFSOURCE_FILE */ char offsets_sourcefile[FILENAMESZ]; char keyword_used_for_classif[ASCIILINESZ] ; double * offsets_total_list ; int * cube_id ; int nb_classified_cubes ; /* Array of nframes int whose values go from 1 to */ /* nb_classified_cubes */ /* Wavelength calibration */ int wavecal_active ; int wavecal_arc_active ; char wavecal_arcfile[FILENAMESZ] ; int wavecal_discard_hi ; int wavecal_discard_lo ; int wavecal_discard_le ; int wavecal_discard_ri ; double * wavecal_disprel ; int wavecal_nb_coeff ; /* Correction of the distortion */ int distortion_active ; int auto_dark_subtraction ; int distor_xmin ; int distor_ymin ; /* Zone to consider for distortion estimation */ int distor_xmax ; int distor_ymax ; /* Combination */ double * main_offset_diff ; int circular_shift ; char final_comb_method[32] ; double average_hi_rejection ; double average_lo_rejection ; /* Extraction of the brightest spectrum */ int spectrum_extr_active ; int detect_bad_left ; int detect_bad_right ; /* number of pixels to reject */ int detect_bad_top ; int detect_bad_bot ; int spectrum_detected ; int spectrum_position ; int spectrum_width ; int res_sky_hi_width ; int res_sky_lo_width ; int res_sky_hi_dist ; int res_sky_lo_dist ; int apply_filter ; /* to apply a median filter before extraction */ int spectrum_extracted ; double * extracted_values ; double * extr_x_coordinate ; /* Post-processing features */ int postproc_active ; int postproc_startviewer ; char postproc_viewer[ASCIILINESZ]; int postproc_gnuplot ; int postproc_statusreport ; /* Saving the output */ int output_averaged_frames ; int output_difference_frames ; char output_basename[FILENAMESZ]; /* Algorithm status */ int status_classification ; int status_average ; int status_disto_slit_curv ; int status_disto_startrace ; int status_wavecal_sky ; int status_wavecal_arc ; int status_wavecal_done ; int status_combination ; int status_extraction ; int status_save ; } lw_spjitter_bb ; /*--------------------------------------------------------------------------- Function prototypes ---------------------------------------------------------------------------*/ /* Data reduction */ static long lw_spjitter_engine(char *, char *, char *, char *) ; static image_t * lw_spjit_comb_create(lw_spjitter_bb *, image_t **) ; static int lw_spjit_save_results(lw_spjitter_bb *, image_t *) ; /* Blackboard handling */ static lw_spjitter_bb * lw_spjitter_bb_create(void) ; static void lw_spjitter_bb_destroy(lw_spjitter_bb *) ; static void lw_spjitter_bb_dump(lw_spjitter_bb *, char *) ; /* Ini file generation and parsing */ static void parse_section_general(dictionary *, lw_spjitter_bb *, int *); static void parse_section_calibration(dictionary*, lw_spjitter_bb *, int*); static void parse_section_classification(dictionary *, lw_spjitter_bb *, int *) ; static void parse_section_distortion(dictionary *, lw_spjitter_bb *, int *) ; static void parse_section_wlcalib(dictionary *, lw_spjitter_bb *, int *) ; static void parse_section_combination(dictionary *, lw_spjitter_bb *, int *) ; static void parse_section_extraction(dictionary *, lw_spjitter_bb *, int *) ; static void parse_section_postproc(dictionary *, lw_spjitter_bb *, int *) ; static void parse_section_output(dictionary*, lw_spjitter_bb *, int*); static void gen_ini(char*, char*) ; static lw_spjitter_bb * lw_spjitter_ini_parse(char *, char *, char *, char *) ; /*---------------------------------------------------------------------------- Main code ---------------------------------------------------------------------------*/ int isaac_lw_spjitter_main(void * dict) { dictionary * d ; int generate_flag ; int timing_flag ; char * ini_name ; char * calib_name ; long nbp ; char argname[10] ; char * name_i ; char * name_o ; int nfiles ; int errors ; int i ; /* Initialize eclipse environment */ hello_world("lw_spjitter", "Spectroscopic Jitter reduction in LW mode"); d = (dictionary*)dict ; /* Get options */ generate_flag = dictionary_getint(d, "arg.generate", 0) ; timing_flag = dictionary_getint(d, "arg.timing", 0) ; calib_name = dictionary_get(d, "arg.calname") ; ini_name = dictionary_get(d, "arg.ininame") ; if (ini_name == NULL) ini_name = strdup("lw_spjitter.ini") ; else ini_name = strdup(ini_name) ; /* IN INI FILE GENERATION MODE */ if (generate_flag) { /* Get input/output file names */ nfiles = dictionary_getint(d, "arg.n", -1) ; if (nfiles != 2) { e_error("expecting one and only one INI file name: aborting") ; errors = -1 ; } else { free(ini_name) ; ini_name = strdup(dictionary_get(d, "arg.1")) ; gen_ini(ini_name, get_eclipse_version()) ; e_comment(0, "spjitter ini file generated [%s]\n", ini_name) ; errors = 0 ; } /* IN REDUCTION MODE */ } else { /* Get input/output file names */ nfiles = dictionary_getint(d, "arg.n", -1) ; if (nfiles < 2) { e_error("missing input file name(s): aborting"); errors = -1 ; } else { /* Loop on input file names */ errors = 0 ; for (i=1 ; i START lw_spjitter engine\n") ; time(&local_t) ; (void)fprintf(stdout, "%s", ctime(&local_t)) ; (void)fprintf(stdout, "pid is %ld\n", (long)getpid()) ; /* CLASSIFICATION OF INPUT FRAMES */ e_comment(0, "---> part 1 of 6: classification of frames") ; /* Remove HC[12] frames from bb->framelist and store them in bb */ /* Get hc1 frames */ if ((bb->hc1_list = get_hc_frames(bb->framelist, bb->nframes, &(bb->nb_hc1), 1)) == NULL) e_warning("cannot get HC1 frames") ; /* Get hc2 frames */ if ((bb->hc2_list = get_hc_frames(bb->framelist, bb->nframes, &(bb->nb_hc2), 2)) == NULL) e_warning("cannot get HC2 frames") ; /* Get rid of HC frames */ if ((bb->framelist_hc_free = reject_hc_frames(bb->framelist, bb->nframes, &(bb->nb_hc_free))) == NULL) { e_error("cannot get rid of HC frames") ; lw_spjitter_bb_destroy(bb) ; return -1 ; } /* Get rid of non-science frames */ if ((bb->framelist_clean = reject_non_science_frames(bb->framelist_hc_free, bb->nb_hc_free, &(bb->nb_clean_frames))) == NULL) { e_error("cannot clean frame list") ; lw_spjitter_bb_destroy(bb) ; return -1 ; } /* Classify the files */ if (spectro_jitter_classification(bb->framelist_clean, bb->nb_clean_frames, bb->offsets_source, bb->keyword_used_for_classif, bb->offsets_sourcefile, &(bb->offsets_total_list), &(bb->nb_classified_cubes), &(bb->cube_id), &(bb->status_classification)) == -1) { e_error("cannot classify the frames") ; lw_spjitter_bb_destroy(bb) ; return -1 ; } /* Load the cubes */ if ((spjit_classified_cubes=spectro_jitter_framelist_to_cubes(bb->framelist_clean, bb->nb_clean_frames, bb->cube_id, bb->nb_classified_cubes)) == NULL) { e_error("cannot create classified cubes") ; lw_spjitter_bb_destroy(bb) ; return -1 ; } /* Nb of computed pixels */ for (i=0 ; inb_classified_cubes ; i++) { total_npix += spjit_classified_cubes[i]->lx * spjit_classified_cubes[i]->ly * spjit_classified_cubes[i]->np ; } /* Apply the flatfield if required */ if (bb->cal_spflat_active) { for (i=0 ; inb_classified_cubes ; i++) { isaac_ff_dark_badpix_handling(&(spjit_classified_cubes[i]), bb->cal_spflat_name, NULL, NULL) ; } bb->divided_by_flat = 1 ; } /* SHIFT & AVERAGE CUBES TO IMAGES */ e_comment(0, "---> part 2 of 6: averaging") ; spjit_averaged_images = spectro_jitter_average(spjit_classified_cubes, bb->nb_classified_cubes, bb->offsets_total_list, bb->cube_id, bb->nb_clean_frames, &(bb->status_average)) ; for (i=0 ; inb_classified_cubes ; i++) { cube_del(spjit_classified_cubes[i]) ; } free(spjit_classified_cubes) ; if (spjit_averaged_images == NULL) { e_error("in distortion correction or averaging") ; lw_spjitter_bb_destroy(bb) ; return -1 ; } /* WAVELENGTH CALIBRATION */ e_comment(0, "---> part 3 of 6: wavelength calibration") ; if (bb->wavecal_active == 1) { if (bb->nb_hc1 > 0) sky_frame = strdup(bb->hc1_list[0]) ; else if (bb->nb_hc2 > 0) sky_frame = strdup(bb->hc2_list[0]) ; else sky_frame = NULL ; if ((bb->wavecal_disprel = spectro_jitter_wlcalib(bb->wavecal_nb_coeff, bb->wavecal_arc_active, bb->wavecal_arcfile, bb->framelist_clean, bb->nb_clean_frames, bb->spectrum_length, sky_frame, 1, bb->wavecal_discard_lo, bb->wavecal_discard_hi, bb->wavecal_discard_le, bb->wavecal_discard_ri, &(bb->status_wavecal_arc), &(bb->status_wavecal_sky), &(bb->status_wavecal_done))) == NULL) { e_error("no wavelength calibration computed: skipping") ; } if (sky_frame != NULL) free(sky_frame) ; } else { e_comment(0, "skipping the wavelength calibration") ; bb->status_wavecal_done = ALGO_SKIPPED ; } /* CREATION OF THE COMBINED IMAGE */ e_comment(0, "---> part 4 of 6: combination") ; spjit_combined_image = lw_spjit_comb_create(bb, spjit_averaged_images) ; for (i=0 ; inb_classified_cubes ; i++) { image_del(spjit_averaged_images[i]) ; } free(spjit_averaged_images) ; if (spjit_combined_image == NULL) { e_error("cannot create the compute image: aborting") ; lw_spjitter_bb_destroy(bb) ; return -1 ; } /* EXTRACTION OF THE BRIGHTEST SPECTRUM */ e_comment(0, "---> part 5 of 6: brightest spectrum extraction") ; if (bb->spectrum_extr_active == 1) { if ((extraction_array = spectro_jitter_extract(spjit_combined_image, bb->spectrum_width, bb->res_sky_lo_dist, bb->res_sky_hi_dist, bb->res_sky_lo_width, bb->res_sky_hi_width, bb->apply_filter, NULL, bb->main_offset_diff, (int)(bb->nb_classified_cubes/2), bb->status_wavecal_done, bb->wavecal_disprel, bb->wavecal_nb_coeff, &(bb->spectrum_position), &(bb->spectrum_detected), &(bb->spectrum_extracted), &(bb->status_extraction))) == NULL) { e_warning("cannot extract the brightest spectrum: skipping") ; } else { bb->extr_x_coordinate = extraction_array->x ; extraction_array->x = NULL ; bb->extracted_values = extraction_array->y ; extraction_array->y = NULL ; double3_del(extraction_array) ; } } else { e_comment(0, "skipping the spectrum extraction") ; bb->spectrum_detected = 0 ; bb->spectrum_extracted = 0 ; bb->status_extraction = ALGO_SKIPPED ; } /* SAVE RESULTS */ e_comment(0, "---> part 6 of 6: saving output data") ; lw_spjit_save_results(bb, spjit_combined_image) ; image_del(spjit_combined_image) ; e_comment(0, "--> STOP lw_spjitter engine\n"); time(&local_t) ; (void)fprintf(stdout, "%s", ctime(&local_t)) ; /* If requested: launch an image viewer on the result */ if ((bb->postproc_startviewer == 1) && (bb->postproc_active == 1)) { sprintf(filename, "%s.fits", bb->output_basename); show_image(filename, bb->postproc_viewer) ; } /* If requested and if the spectrum has been extracted: plot */ /* the spectrum */ if ((bb->postproc_gnuplot == 1) && (bb->postproc_active == 1) && (bb->spectrum_extracted == 1)) { plot_signal(bb->extr_x_coordinate, bb->extracted_values, bb->spectrum_length, "Wavelength (in angstroms)", "Extracted spectrum (in ADU)") ; } /* If requested, output a status file as basename_status.ascii */ if ((bb->postproc_statusreport == 1) && (bb->postproc_active == 1)) { sprintf(filename, "%s_status.ascii", bb->output_basename); lw_spjitter_bb_dump(bb, filename); } lw_spjitter_bb_destroy(bb); return total_npix ; } /*-------------------------------------------------------------------------*/ /** @name lw_spjit_comb_create @memo Create the combined image @param bb blackboard @param averaged averaged images @return Combined image, and image used for spectrum detection */ /*--------------------------------------------------------------------------*/ static image_t * lw_spjit_comb_create( lw_spjitter_bb * bb, image_t ** averaged) { image_t ** difference ; image_t ** corrected ; image_t * tmp_image ; char * estim_frame ; cube_t * combined_shifted_list ; image_t * combined ; char out_aver[FILENAMESZ] ; char out_diff[FILENAMESZ] ; qfits_header* hdr ; double off_val1 ; double off_val2 ; double off_val3 ; double off_val4 ; double offset ; int lo_rej, hi_rej; int i, j ; /* Output or not the averaged images */ if (bb->output_averaged_frames == 1) { /* Read the fits header of the first input frame */ if ((hdr = qfits_header_read(bb->framelist_clean[0])) == NULL) { e_error("cannot output averaged frames") ; } else { sprintf(out_aver, "%s_aver_XX.fits", bb->output_basename) ; isaac_header_for_image(hdr) ; isaac_pro_fits(hdr, out_aver, "REDUCED", NULL, isaac_spec_lw_jitter_aver, "OK", "lw_spec_obs_chopnod", bb->nframes, NULL, NULL) ; /* For each averaged frame */ for (i=0 ; inb_classified_cubes ; i++) { /* Name of the out files */ sprintf(out_aver, "%s_aver_%d.fits", bb->output_basename, i+1) ; /* Write averaged frames on disk */ image_save_fits_hdrdump(averaged[i], out_aver, hdr, BPP_DEFAULT) ; } qfits_header_destroy(hdr) ; } } /* Compute the difference between the averaged frames */ e_comment(1, "computing frame differences"); /* Allocate the difference frames */ difference = malloc(bb->nb_classified_cubes*sizeof(image_t*)/2) ; off_val1 = get_offset_in_list(bb->offsets_total_list, bb->nframes, bb->cube_id, 0, 0) ; off_val2 = get_offset_in_list(bb->offsets_total_list, bb->nframes, bb->cube_id, 1, 0) ; for (i=0 ; inb_classified_cubes/2 ; i++) { off_val3 = get_offset_in_list(bb->offsets_total_list, bb->nframes, bb->cube_id, 2*i, 0) ; off_val4 = get_offset_in_list(bb->offsets_total_list, bb->nframes, bb->cube_id, 2*i+1, 0) ; /* Subtraction */ /* Test the offsets order to choose between a-b or b-a */ if ((off_val2 - off_val1) * (off_val4 - off_val3) > 0) { tmp_image = image_sub(averaged[2*i], averaged [2*i+1]) ; } else { tmp_image = image_sub(averaged[2*i+1], averaged [2*i]) ; } /* Division by 2 */ difference[i] = image_cst_op(tmp_image, 2.0, '/') ; image_del(tmp_image) ; if (difference[i] == NULL) { e_error("subtraction failed - aborting") ; for (j=0 ; jstatus_combination = ALGO_FAILED ; return NULL ; } } /* Correction of the distortion */ if (bb->distortion_active) { if (bb->nb_hc1 > 0) estim_frame = strdup(bb->hc1_list[0]) ; else if (bb->nb_hc2 > 0) estim_frame = strdup(bb->hc2_list[0]) ; else estim_frame = NULL ; if ((corrected = isaac_distortion_estim_corr(difference, bb->nb_classified_cubes, estim_frame, bb->cal_arc_active, bb->cal_arc_name, bb->cal_startrace_active, bb->cal_startrace_name, bb->distor_xmin, bb->distor_ymin, bb->distor_xmax, bb->distor_ymax, bb->auto_dark_subtraction, &(bb->status_disto_slit_curv), &(bb->status_disto_startrace))) != NULL) { /* Free the difference images */ for (i=0 ; inb_classified_cubes ; i++) { image_del(difference[i]) ; } free(difference) ; } else { e_warning("Distortion correction failed") ; corrected = difference ; } if (estim_frame != NULL) free(estim_frame) ; } else { corrected = difference ; } /* Output or not the difference corrected images */ if (bb->output_difference_frames == 1) { /* Read the fits header of the first input frame */ if ((hdr = qfits_header_read(bb->framelist_clean[0])) == NULL) { e_error("cannot output difference frames") ; } else { sprintf(out_diff, "%s_diff_XX.fits", bb->output_basename) ; isaac_header_for_image(hdr) ; isaac_pro_fits(hdr, out_diff, "REDUCED", NULL, isaac_spec_lw_jitter_diff, "OK", "lw_spec_obs_chopnod", bb->nframes, NULL, NULL) ; /* For each averaged frame */ for (i=0 ; inb_classified_cubes/2 ; i++) { /* Name of the out files */ sprintf(out_diff, "%s_diff_%d.fits", bb->output_basename, i+1) ; /* Write averaged frames on disk */ image_save_fits_hdrdump(corrected[i], out_diff, hdr, BPP_DEFAULT) ; } qfits_header_destroy(hdr) ; } } /* Shift the combined (difference - distortion corrected) images */ /* Allocate the combined list */ combined_shifted_list = cube_new(corrected[0]->lx, corrected[0]->ly, bb->nb_classified_cubes/2) ; /* Allocate the offsets array */ bb->main_offset_diff = malloc((bb->nb_classified_cubes/2)*sizeof(double)) ; off_val1 = get_offset_in_list(bb->offsets_total_list, bb->nframes, bb->cube_id, 0, 0) ; off_val2 = get_offset_in_list(bb->offsets_total_list, bb->nframes, bb->cube_id, 1, 0) ; for (i=0 ; inb_classified_cubes/2 ; i++) { off_val3 = get_offset_in_list(bb->offsets_total_list, bb->nframes, bb->cube_id, 2*i, 0) ; off_val4 = get_offset_in_list(bb->offsets_total_list, bb->nframes, bb->cube_id, 2*i+1, 0) ; /* Test the offsets order to choose between a-b or b-a */ if ((off_val2-off_val1)*(off_val4-off_val3) > 0) { bb->main_offset_diff[i] = off_val4 - off_val3 ; offset = off_val3 - off_val1 ; } else { bb->main_offset_diff[i] = off_val3 - off_val4 ; offset = off_val4 - off_val1 ; } combined_shifted_list->plane[i] = image_shift( corrected[i], 0, offset, (double*)NULL) ; } /* Free corrected */ for (i=0 ; inb_classified_cubes/2 ; i++) { image_del(corrected[i]) ; } free(corrected) ; /* Final combination of the combined shifted images */ /* If nb of planes to combine is < 3 -> linear average */ if (combined_shifted_list->np < 3) { bb->final_comb_method[0] = 'l' ; } if (bb->final_comb_method[0] == 'r') { lo_rej = (int)(bb->average_lo_rejection * combined_shifted_list->np); hi_rej = (int)(bb->average_hi_rejection * combined_shifted_list->np); combined = cube_avg_reject(combined_shifted_list, lo_rej, hi_rej); } else if (bb->final_comb_method[0] == 'l') { combined = cube_avg_linear(combined_shifted_list) ; } else if (bb->final_comb_method[0] == 'm') { combined = cube_avg_median(combined_shifted_list) ; } else { e_warning("final combination method not recognized - use median") ; combined = cube_avg_median(combined_shifted_list) ; } cube_del(combined_shifted_list) ; if (combined == NULL) { e_error("in averaging the combined images") ; free(combined) ; bb->status_combination = ALGO_FAILED ; return NULL ; } bb->status_combination = ALGO_OK ; /* Free and return */ return combined ; } /*-------------------------------------------------------------------------*/ /** @name lw_spjit_save_results @memo write outpu files @param bb blackboard @param combined combined image @return -1 in error case */ /*--------------------------------------------------------------------------*/ static int lw_spjit_save_results( lw_spjitter_bb * bb, image_t * combined) { qfits_header * fh ; char outname[1024] ; qfits_table * table ; qfits_col * col ; double * out_table[2] ; FILE * paf ; char pafname[FILENAMESZ] ; char * mjd_obs ; char wl_method[FILENAMESZ] ; char * strvar ; char cval[80] ; framelist * lnames ; int i ; sprintf(outname, "%s.fits", bb->output_basename); /* Create the framelist object */ lnames = framelist_new(bb->nframes) ; for (i=0 ; inframes ; i++) { lnames->name[i] = strdup(bb->framelist[i]) ; } /* WRITE THE OUTPUT FITS FILE (combined image) */ /* Read FITS header from the first frame */ fh = qfits_header_read(bb->framelist_clean[0]); /* Prepare it for image output */ if (isaac_header_for_image(fh) == -1) { e_error("in writing the output fits file") ; qfits_header_destroy(fh) ; framelist_del(lnames) ; return -1 ; } /* Update FITS header with PRO keywords for the combined image */ isaac_pro_fits(fh, outname, "REDUCED", NULL, isaac_spec_lw_jitter_comb, "OK", "lw_spec_obs_chopnod", lnames->n, lnames, NULL); /* Write HISTORY keywords in the header */ if (isaac_add_files_history(fh, lnames) == -1) { e_warning("cannot write HISTORY keywords in out file") ; } if (bb->wavecal_disprel != NULL) { sprintf(cval, "DISPCOE1= %g", bb->wavecal_disprel[0]) ; qfits_header_add(fh, "HISTORY", cval, NULL, NULL) ; sprintf(cval, "DISPCOE2= %g", bb->wavecal_disprel[1]) ; qfits_header_add(fh, "HISTORY", cval, NULL, NULL) ; } /* Write the file */ image_save_fits_hdrdump(combined, outname, fh, BPP_DEFAULT) ; qfits_header_destroy(fh) ; e_comment(0, "combined image produced: [%s]", outname) ; /* WRITE THE FITS OUT TABLE (extracted spectrum) */ if (bb->status_extraction == ALGO_OK) { sprintf(outname, "%s.tfits", bb->output_basename); /* Table informations */ table = qfits_table_new(outname, QFITS_BINTABLE, -1, 2, bb->spectrum_length) ; col = table->col ; for (i=0 ; inc ; i++) { qfits_col_fill(col, 1, sizeof(double), TFITS_BIN_TYPE_D, " ", " ", " ", " ", 0, 0.0, 0, 1.0, i*sizeof(double)) ; col++ ; } col = table->col ; sprintf(col->tlabel, "X_coordinate") ; col++ ; sprintf(col->tlabel, "Extracted_spectrum_value") ; out_table[0] = bb->extr_x_coordinate ; out_table[1] = bb->extracted_values ; /* Read the input header */ if ((fh = qfits_header_read(bb->framelist_clean[0])) == NULL) { e_error("in writing the output fits file") ; qfits_table_close(table) ; framelist_del(lnames) ; return -1 ; } /* Prepare it for table output */ if (isaac_header_for_table(fh) == -1) { e_error("in writing the output fits file") ; qfits_header_destroy(fh) ; qfits_table_close(table) ; framelist_del(lnames) ; return -1 ; } /* Write the PRO keywords in the header */ if (isaac_pro_fits(fh, outname, "REDUCED", NULL, isaac_spec_lw_jitter_extr, "OK", "lw_spec_obs_chopnod", lnames->n, lnames, NULL) == -1) { e_error("in writing PRO keywords in output file") ; qfits_header_destroy(fh) ; framelist_del(lnames) ; qfits_table_close(table) ; return -1 ; } /* Write the HISTORY keywords with the input file names */ if (isaac_add_files_history(fh, lnames) == -1) { e_warning("cannot write HISTORY keywords in out file") ; } /* Write the file on disk */ if (qfits_save_table_hdrdump((void**)out_table, table, fh) == -1) { e_error("cannot write file: %s", outname) ; qfits_header_destroy(fh) ; qfits_table_close(table) ; framelist_del(lnames) ; return -1 ; } qfits_table_close(table) ; qfits_header_destroy(fh) ; e_comment(0, "extracted spectrum table produced: [%s]", outname) ; } framelist_del(lnames) ; /* Write the PAF file */ sprintf(pafname, "%s.paf", get_rootname(outname)) ; paf = paf_print_header( pafname, "ISAAC/lw_spjitter", "lw_spjitter recipe results") ; if (paf == NULL) { e_warning("cannot output PAF file") ; } else { fprintf(paf, "\n"); /* ARCFILE */ strvar = isaac_get_arcfile(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "ARCFILE \"%s\" \n", strvar) ; } /* MJD-OBS */ mjd_obs = isaac_get_mjdobs(bb->framelist_clean[0]) ; if (mjd_obs!=NULL) { fprintf(paf, "MJD-OBS %s; # Obs start\n\n", mjd_obs); } else { fprintf(paf, "MJD-OBS 0.0; # Obs start unknown\n\n"); } /* INSTRUME keyword */ strvar = isaac_get_instrument(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "INSTRUME \"%s\" \n", strvar) ; } /* TPL.ID */ strvar = isaac_get_templateid(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "TPL.ID \"%s\" \n", strvar) ; } /* TPL.NEXP */ strvar = isaac_get_numbexp(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "TPL.NEXP %s \n", strvar) ; } /* DPR.CATG */ strvar = isaac_get_dpr_catg(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "DPR.CATG \"%s\" \n", strvar) ; } /* DPR.TYPE */ strvar = isaac_get_dpr_type(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "DPR.TYPE \"%s\" \n", strvar) ; } /* DPR.TECH */ strvar = isaac_get_dpr_tech(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "DPR.TECH \"%s\" \n", strvar) ; } /* Add PRO.CATG */ fprintf(paf, "PRO.CATG \"%s\" ;# Product category\n", isaac_get_pro_catg_value(isaac_spec_lw_jitter_qc)) ; /* Add the date */ fprintf(paf, "DATE-OBS \"%s\" ;# Date\n", isaac_get_date_obs(bb->framelist_clean[0])) ; /* INS.GRAT.NAME */ strvar = isaac_get_resolution(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "INS.GRAT.NAME \"%s\" \n", strvar) ; } /* INS.GRAT.WLEN */ fprintf(paf, "INS.GRAT.WLEN %g \n", isaac_get_central_wavelength(bb->framelist_clean[0])) ; /* QC.STDNAME */ strvar = isaac_get_obs_targ_name(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "QC.STDNAME \"%s\" \n", strvar) ; } /* INS.FILTER.ID */ strvar = isaac_get_filter(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "INS.FILTER.ID \"%s\" \n", strvar) ; } /* INS.OPTI1.ID */ strvar = isaac_get_optical_id(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "INS.OPTI1.ID \"%s\" \n", strvar) ; } /* OBS.ID */ strvar = isaac_get_obs_id(bb->framelist_clean[0]) ; if (strvar != NULL) { fprintf(paf, "OBS.ID \"%s\" \n", strvar) ; } if (bb->status_wavecal_done==ALGO_OK) { if (bb->status_wavecal_sky==ALGO_OK) { sprintf(wl_method, "sky lines") ; } else if (bb->status_wavecal_arc==ALGO_OK) { sprintf(wl_method, "arc file") ; } else { sprintf(wl_method, "physical model") ; } /* QC.CENWL */ fprintf(paf, "QC.WLEN %g \n", (double)(bb->wavecal_disprel[0]+bb->wavecal_disprel[1]*(1024/2))); /* QC.DISPCO1 */ fprintf(paf, "QC.DISPCO1 %g \n", (double)bb->wavecal_disprel[0]); /* QC.DISPCO2 */ fprintf(paf, "QC.DISPCO2 %g \n", (double)bb->wavecal_disprel[1]); } else { sprintf(wl_method, "none") ; } /* QC.WLMETTHOD */ fprintf(paf, "QC.WLMETHOD \"%s\" \n", wl_method) ; fclose(paf) ; } bb->status_save = ALGO_OK ; return 0 ; } /*-------------------------------------------------------------------------*/ /** @name lw_spjitter_bb_create @memo allocate and initialize with zeros a blackboard struct @return newly allocated blackboard structure @doc deallocation MUST be performed by lw_spjitter_bb_destroy() */ /*--------------------------------------------------------------------------*/ static lw_spjitter_bb * lw_spjitter_bb_create(void) { lw_spjitter_bb * bb ; bb = calloc(1, sizeof(lw_spjitter_bb)); return bb ; } /*-------------------------------------------------------------------------*/ /** @name lw_spjitter_bb_destroy @memo deallocate a blackboard structure allocated previously by lw_spjitter_bb_create() @param bb blackboard @return void */ /*--------------------------------------------------------------------------*/ static void lw_spjitter_bb_destroy(lw_spjitter_bb * bb) { int i ; /* Deallocate frame names */ if (bb->nframes>0) { if (bb->framelist!=NULL) { for (i=0 ; inframes ; i++) { if (bb->framelist[i] != NULL) free(bb->framelist[i]); } free(bb->framelist); } } /* Deallocate HC free frame names */ if (bb->nb_hc_free>0) { if (bb->framelist_hc_free != NULL) { for (i=0 ; inb_hc_free ; i++) { if (bb->framelist_hc_free[i] != NULL) free(bb->framelist_hc_free[i]); } free(bb->framelist_hc_free); } } /* Deallocate clean frame names */ if (bb->nb_clean_frames>0) { if (bb->framelist_clean != NULL) { for (i=0 ; inb_clean_frames ; i++) { if (bb->framelist_clean[i] != NULL) free(bb->framelist_clean[i]); } free(bb->framelist_clean); } } /* Deallocate half cycle frame names */ if (bb->nb_hc1>0) { if (bb->hc1_list!=NULL) { for (i=0 ; inb_hc1 ; i++) { if (bb->hc1_list[i] != NULL) free(bb->hc1_list[i]); } free(bb->hc1_list); } } /* Deallocate half cycle frame names */ if (bb->nb_hc2>0) { if (bb->hc2_list!=NULL) { for (i=0 ; inb_hc2 ; i++) { if (bb->hc2_list[i] != NULL) free(bb->hc2_list[i]); } free(bb->hc2_list); } } /* Deallocate offsets lists */ if (bb->cube_id != NULL) free(bb->cube_id) ; if (bb->offsets_total_list != NULL) free(bb->offsets_total_list) ; if (bb->main_offset_diff != NULL) free(bb->main_offset_diff) ; /* Deallocate wavelength calibration */ if (bb->wavecal_disprel != NULL) free(bb->wavecal_disprel) ; /* Deallocate extracted values */ if (bb->extracted_values != NULL) free(bb->extracted_values) ; if (bb->extr_x_coordinate != NULL) free(bb->extr_x_coordinate) ; free(bb); return ; } /*-------------------------------------------------------------------------*/ /** @name lw_spjitter_bb_dump @memo dump out contents of a blackboard structure to stdout @param bb pointer to allocated blackboard structure @param filename file name @return void */ /*--------------------------------------------------------------------------*/ #define IBOOL(i) ((i)? "yes" : "no") static void lw_spjitter_bb_dump(lw_spjitter_bb * bb, char * filename) { FILE * dm ; int i ; if (bb == NULL) return ; if (filename == NULL) dm = stderr ; else { /* Open the file for dumping */ dm = fopen(filename, "w") ; if (dm==NULL) { e_error("cannot create file [%s]: ", filename); e_error("dumping config to stderr"); dm = stderr ; } } fprintf(dm, "#\n\ # lw_spjitter_config from pid [%ld]\n", (long)getpid()); fprintf(dm, "#\n"); fprintf(dm, "\n\ # The configuration: \n\ # WLCalibrationSky != Ok ;\n\ # WLCalibrationArc != Ok ;\n\ # WLCalibrationDone = Ok ;\n\ # means that the wavelength calibration has been done with the\n\ # physical model.\n\n") ; fprintf(dm, "[AlgorithmStatus]\n"); fprintf(dm, "Classification = %s ;\n", algo_stat(bb->status_classification)) ; fprintf(dm, "Averaging = %s ;\n", algo_stat(bb->status_average)) ; fprintf(dm, "SlitCurvDist = %s ;\n", algo_stat(bb->status_disto_slit_curv)) ; fprintf(dm, "StartraceDist = %s ;\n", algo_stat(bb->status_disto_startrace)) ; fprintf(dm, "WLCalibSky = %s ;\n", algo_stat(bb->status_wavecal_sky)) ; fprintf(dm, "WLCalibArc = %s ;\n", algo_stat(bb->status_wavecal_arc)) ; fprintf(dm, "WLCalibDone = %s ;\n", algo_stat(bb->status_wavecal_done)) ; fprintf(dm, "Combination = %s ;\n", algo_stat(bb->status_combination)) ; fprintf(dm, "Extraction = %s ;\n", algo_stat(bb->status_extraction)) ; fprintf(dm, "Save = %s ;\n\n\n", algo_stat(bb->status_save)) ; fprintf(dm, "[Frames]\n\n"); fprintf(dm, "FileList = %s ;\n", bb->frames_filelist) ; fprintf(dm, "NFrames = %d ;\n\n", bb->nframes); for (i=0 ; inframes ; i++) fprintf(dm, "Frame:%03d = %s ;\n", i+1, bb->framelist[i]) ; fprintf(dm, "\n"); fprintf(dm, "NHalfCycleFrames1 = %d ;\n", bb->nb_hc1) ; for (i=0 ; inb_hc1 ; i++) fprintf(dm, "HC1Frame:%03d = %s ;\n", i+1, bb->hc1_list[i]) ; fprintf(dm, "\n"); fprintf(dm, "NHalfCycleFrames2 = %d ;\n", bb->nb_hc2) ; for (i=0 ; inb_hc2 ; i++) fprintf(dm, "HC2Frame:%03d = %s ;\n", i+1, bb->hc2_list[i]) ; fprintf(dm, "\n"); fprintf(dm, "SpectrumLength = %d ;\n\n\n", bb->spectrum_length) ; fprintf(dm, "[Calibration]\n\n") ; fprintf(dm, "CalArcActive = %s ;\n", IBOOL(bb->cal_arc_active)) ; fprintf(dm, "CalStartraceActive = %s ;\n", IBOOL(bb->cal_startrace_active)) ; fprintf(dm, "CalFlatActive = %s ;\n\n", IBOOL(bb->cal_spflat_active)) ; fprintf(dm, "ArcTable = %s ;\n", bb->cal_arc_name) ; fprintf(dm, "StarTraceTable = %s ;\n", bb->cal_startrace_name) ; fprintf(dm, "MasterSpFlat = %s ;\n\n\n", bb->cal_spflat_name) ; fprintf(dm, "[Classification]\n\n\n") ; if (bb->offsets_source == OFFSOURCE_HEADER) fprintf(dm, "Source = header ;\n") ; else if (bb->offsets_source == OFFSOURCE_FILE) fprintf(dm, "Source = file ;\n") ; fprintf(dm, "OffsetKeyword = %s ;\n", bb->keyword_used_for_classif) ; fprintf(dm, "OffsetFile = %s ;\n\n", bb->offsets_sourcefile) ; for (i=0 ; inb_hc_free ; i++) fprintf(dm, "OffsetsTotalList:%03d = %g ;\n", i+1, bb->offsets_total_list[i]) ; fprintf(dm, "\n"); fprintf(dm, "NbClassifiedCubes = %d ;\n\n\n", bb->nb_classified_cubes) ; fprintf(dm, "[Distortion]\n\n") ; fprintf(dm, "DistoActive = %s ;\n\n", IBOOL(bb->distortion_active)); fprintf(dm, "AutoDarkSubtraction = %s ;\n\n", IBOOL(bb->auto_dark_subtraction)) ; fprintf(dm, "XMin = %d ;\n", bb->distor_xmin) ; fprintf(dm, "YMin = %d ;\n", bb->distor_ymin) ; fprintf(dm, "XMax = %d ;\n", bb->distor_xmax) ; fprintf(dm, "YMax = %d ;\n\n", bb->distor_ymax) ; fprintf(dm, "[WavelengthCalibration]\n\n") ; fprintf(dm, "WaveCalActive = %s ;\n\n", IBOOL(bb->wavecal_active)) ; fprintf(dm, "WaveArcTable = %s ;\n\n", bb->wavecal_arcfile) ; fprintf(dm, "DiscardHigh = %d ;\n", bb->wavecal_discard_hi) ; fprintf(dm, "DiscardLow = %d ;\n", bb->wavecal_discard_lo) ; fprintf(dm, "DiscardLeft = %d ;\n", bb->wavecal_discard_le) ; fprintf(dm, "DiscardRight = %d ;\n\n", bb->wavecal_discard_ri) ; if (bb->wavecal_disprel != NULL) for (i=0 ; iwavecal_nb_coeff ; i++) fprintf(dm, "DispCoeff%d = %g ;\n", i+1, bb->wavecal_disprel[i]) ; fprintf(dm, "\n\n"); fprintf(dm, "[Combination]\n\n") ; fprintf(dm, "DividedByFlat = %s ;\n\n", IBOOL(bb->divided_by_flat)) ; for (i=0 ; inb_classified_cubes/2 ; i++) fprintf(dm, "MainOffsetDiff%d = %g ;\n", i+1, bb->main_offset_diff[i]) ; fprintf(dm, "\nCircularShift = %s ;\n\n", IBOOL(bb->circular_shift)) ; if (bb->final_comb_method[0] == 'l') fprintf(dm, "FinalCombMethod = linear ;\n") ; else if (bb->final_comb_method[0] == 'r') fprintf(dm, "FinalCombMethod = rejection ;\n") ; else if (bb->final_comb_method[0] == 'm') fprintf(dm, "FinalCombMethod = median ;\n") ; fprintf(dm, "AverageHiRejection = %g ;\n", bb->average_hi_rejection) ; fprintf(dm, "AverageLoRejection = %g ;\n\n\n", bb->average_lo_rejection) ; fprintf(dm, "[SpectrumExtraction]\n\n") ; fprintf(dm, "SpecExtrActive = %s ;\n\n", IBOOL(bb->spectrum_extr_active)) ; fprintf(dm, "BadLeft = %d ;\n", bb->detect_bad_left) ; fprintf(dm, "BadRight = %d ;\n", bb->detect_bad_right) ; fprintf(dm, "BadTop = %d ;\n", bb->detect_bad_top) ; fprintf(dm, "BadBot = %d ;\n\n", bb->detect_bad_bot) ; fprintf(dm, "SpectrumDetected = %s ;\n\n", IBOOL(bb->spectrum_detected)) ; fprintf(dm, "SpectrumPosition = %d ;\n", bb->spectrum_position) ; fprintf(dm, "SpectrumWidth = %d ;\n", bb->spectrum_width) ; fprintf(dm, "ResSkyHiWidth = %d ;\n", bb->res_sky_hi_width) ; fprintf(dm, "ResSkyLoWidth = %d ;\n", bb->res_sky_lo_width) ; fprintf(dm, "ResSkyHiDist = %d ;\n", bb->res_sky_hi_dist) ; fprintf(dm, "ResSkyLoDist = %d ;\n", bb->res_sky_lo_dist) ; fprintf(dm, "ApplyFilter = %d ;\n", bb->apply_filter) ; fprintf(dm, "\nSpectrumExtracted = %s ;\n\n", IBOOL(bb->spectrum_extracted)) ; fprintf(dm, "[PostProcessing]\n\n") ; fprintf(dm, "PostProcActive = %s ;\n\n", IBOOL(bb->postproc_active)) ; fprintf(dm, "ProduceStatusReport = %s ;\n", IBOOL(bb->postproc_statusreport)) ; fprintf(dm, "StartViewer = %s ;\n", IBOOL(bb->postproc_startviewer)) ; fprintf(dm, "StartCommand = %s ;\n", bb->postproc_viewer) ; fprintf(dm, "PlotSpectrum = %s ;\n\n", IBOOL(bb->postproc_gnuplot)) ; fprintf(dm, "[Output]\n\n") ; fprintf(dm, "BaseName = %s ;\n\n", bb->output_basename) ; fprintf(dm, "OutAveraged = %s ;\n", IBOOL(bb->output_averaged_frames)) ; fprintf(dm, "OutDifference = %s ;\n\n", IBOOL(bb->output_difference_frames)) ; fprintf(dm, "# end of file\n"); if (dm!=stderr) fclose(dm); return ; } #undef IBOOL /*-------------------------------------------------------------------------*/ /** @name gen_ini @memo Generate an ini file @param filename INI file name @param ecl_version eclipse version @return nothing */ /*--------------------------------------------------------------------------*/ static void gen_ini( char * filename, char * ecl_version) { FILE * f ; f=fopen(filename,"w"); fprintf(f, "#\n" "# Configuration file for lw_spjitter\n" "# %s\n" "#\n" "\n", create_timestamp()); fprintf(f, "#\n" "# Please check out the WWW documentation at:\n" "# http://www.eso.org/eclipse/\n" "# \n" "# ----- A note about output names\n" "# \n" "# All files created by 'lw_spjitter' will be named according to the\n" "# following convention:\n" "# basename_type.fits or basename_type.tfits or basename_type.ascii\n" "# where basename is declared in the output option\n" "# and type depends on what kind of data are contained in the file.\n" "# \n" "# The following output files are possible (if requested):\n" "# \n" "# * basename.fits is the combined image\n" "# \n" "# * basename.tfits is a 2 columns fits table contaninig the extracted\n" "# spectrum (only if the spectrum has been extracted)\n" "# * basename_aver[X].fits are the averaged frames\n" "# * basename_diff[X].fits are the difference frames\n" "# \n" "# \n" "# \n" "# ----- A note about FITS keywords\n" "# \n" "# Either provide the name of a FITS keyword when\n" "# requested, or provide the string 'isaac-default'\n" "# that means that 'lw_spjitter' should try to look for the\n" "# information with ISAAC conventions.\n" "# \n" "# The 'lw_spjitter' code knows about all current and past\n" "# conventions for ISAAC keywords and will try to find\n" "# the corresponding value no matter what the file version\n" "# is.\n" "# Example: in the first six months, ISAAC produces the\n" "# RA information in a keyword named 'RA' but later in a keyword\n" "# named 'ISAACRA'. If you indicate 'isaac-default', lw_spjitter will\n" "# first look for a valid keyword named 'RA' and if it does not\n" "# find it, will look for 'ISAACRA'.\n" "# If getting the information requires some more complex search\n" "# in the header, it will be implemented in the code itself.\n" "\n" "# Input and Calibration files\n" "# \n" "# Input file names are stored in a separate file.\n" "# The input frames list is an ASCII file containing\n" "# the files name in first column and the optionally the frames type\n" "# in second column.\n" "# \n" "# See the -i option to change that name when you generate\n" "# this file (with the -g or --generate option).\n" "# put file names are stored in a separate file.\n" "# \n" "[Eclipse]\n" "Version = %s\n", ecl_version); fprintf(f, "# ----- Section [Calibration]\n" "# Provide in this section names for the Arc table, the star trace table,\n" "# and the spectroscopic master flat you want to use.\n" "# These three files come respectively from the calibration recipes\n" "# lw_is_spec_arc\n" "# lw_is_spec_startrace (name_coef.tfits)\n" "# is_spec_flat.\n" "# \n" "# If you do not have such files, specify 'none' as filename\n" "# \n" "# See the -c option to generate a default lw_spjitter.ini\n" "# containing calibration file names as produced by the\n" "# Reduction Block Scheduler in the DataFlow environment.\n" "# If these files are given through this calibration ASCII file\n" "# (with the -c option), their type has to be specified in the\n" "# second column:\n" "# LW_ARC_COEF for the arc file\n" "# LW_STAR_TRACE_COEF for the star trace file\n" "# LW_MASTER_SP_FLAT for the master flat field\n" "# The wavelength calibration is computed by using the ARC calibration\n" "# table (produced by lw_is_spec_arc) specified here by WavecalArcFile\n" "# and NOT by ArcTable in the calibration section.\n" "# ArcTable can only be used for the distortion correction!\n" "\n" "[Calibration]\n" "ArcTable = none ; arc table name\n" "StarTraceTable = none ; StarTrace table name\n" "MasterSpFlatTable = none ; Master flat name\n") ; fprintf(f, "\n" "\n" "# ----- Section [Classification]\n" "# Classification of the input files\n" "# \n" "# In the input file list, only the existing files whose DPR type is\n" "# \"OBJECT\" or \"STD\" (checked in the headers) are taken into account\n" "# (others are rejected).\n" "# \n" "# The input files are classified in series of frames.\n" "# The classification is done according the offsets read in the\n" "# header of the input files or in a file provided by the user.\n" "# This offsets file is a 1 column ASCII file with on each line the offset\n" "# value.\n" "#\n" "# The classification is done according the offsets order encountered.\n" "# EXPLANATIONS .............\n" "#\n" "# OffsetKeyword is the keyword to look for in the headers to get\n" "# the offsets.\n" "\n" "[Classification]\n" "\n" "# The classification is done according the offsets read in the\n" "# header of the input files or in a file provided by the user\n" "[Classification/Source]\n" "Select = header\n" "# OffsetKeyword is the keyword to look for in the headers to get\n" "# the offsets.\n" "OffsetKeyword = isaac-default\n" "# OffsetFile is the ascii file name with offsets or none\n" "OffsetFile = none\n" "\n" "\n" "\n"); fprintf(f, "# ----- Section [WavelengthCalibration]\n" "# The wavelength calibration can be switched off (Select to no).\n" "# \n" "# The wavelength calibration is computed by using the ARC calibration\n" "# table (produced by lw_is_spec_arc) specified here by WavecalArcFile\n" "# and NOT by ArcTable in the calibration section.\n" "# ArcTable can only be used for the distortion correction!\n" "#\n" "# If not available, a wavelength calibration using the sky\n" "# lines in the input frames is computed.\n" "# If this is also not possible (no sky lines), a rough\n" "# wavelength calibration is done using a physical model,\n" "# using the central wavelength, resolution and slit used\n" "# (header infos)).\n" "# \n" "# The discarded pixels concern only the calibration using the sky lines.\n" "# They are discarded to improve the sky lines detection.\n" "\n" "[WavelengthCalibration]\n" "Select = yes\n" "WavecalArcFile = none\n" "DiscardHigh = 50 ; number of pixels to discard at the top and\n" "DiscardLow = 50 ; the bottom of the image used for calibration\n" "DiscardLeft = -1 ; left columns set to 0 before lines matching\n" "DiscardRight = -1 ; same as left. -1 for automatic mode\n" "WavecalNbCoeff = 2 ; nb of coeff of the dispersion polynomial\n" "\n" "\n"); fprintf(f, "# ----- Section [Combination]\n" "# DISTORTION CORRECTION NOT IMPLEMENTED YET - SWITCHED OFF\n" "# For each combined image that has to be created, the two corresponding\n" "# averaged images are subtracted (1-2 and 2-1), the difference images \n" "# are distortion-corrected, and then, 2-1 is shifted and added to 1-2.\n" "#\n" "# It is possible to switch off the distortion correction (DistoActive \n" "# set to no).\n" "#\n" "# If the ARC calibration table is provided, the distortion correction is\n" "# computed with the polynomial stored in it.\n" "# If not, the distortion polynomial is found using the sky lines in the \n" "# input frames. If there are not enough sky lines the distortion correction\n" "# is skipped.\n" "# The AutoDarkSubtraction is a method used to clean the raw image before \n" "# the sky lines detection.\n" "#\n" "# The bad zone specified is rejected only for the distortion polynomial \n" "# estimation.\n" "#\n" "# The shift can be circular or not.\n" "# \n" "# The obtained combined images are shifted according the offsets\n" "# differences, and combined together with a linear average, a rejection\n" "# average or a median. The rejection coefficients only concern the\n" "# rejection average method.\n" "\n" "[Combination]\n" "CircularShift = no\n" "\n" "[Combination/FinalCombMethod]\n" "Select = median ; median - rejection - linear\n" "AverageHiRejection = 0.1 ; high rejection rate for averaging\n" "AverageLoRejection = 0.1 ; low rejection rate for averaging\n" "\n"); fprintf(f, "\n" "[Combination/Distortion]\n" "Select = no ; activate the distortion correction\n" "AutoDarkSubtraction = yes ; auto. dark subt. before sky lines detection.\n" "XMin = 1\n" "YMin = 50\n" "XMax = 1024\n" "YMax = 975\n" "\n"); fprintf(f, "\n" "\n" "# ----- Section [SpectrumExtraction]\n" "# The spectrum extraction can be switched off (SpecExtrActive set to no).\n" "#\n" "# The bad region specified is only used for the spectrum detection.\n" "# If the spectrum position specified is negative, 'spjitter' tries\n" "# to detect the brightest spectrum in the combined image. If not,\n" "# the spectrum to extract is assumed to stand at the specified position.\n" "# The residual sky coefficients are set to default values if negative in\n" "# this INI file.\n" "\n" "[SpectrumExtraction]\n" "Select = yes\n" "SpectrumWidth = 10\n" "BadTop = 0\n" "BadLeft = 50\n" "BadRight = 50\n" "BadBot = 0\n" "ResSkyHiWidth = -1 ; residual sky width above the spectrum \n" "ResSkyHiDist = -1 ; residual sky width above the spectrum \n" "ResSkyLoWidth = -1 ; residual sky width below the spectrum \n" "ResSkyLoDist = -1 ; residual sky width below the spectrum \n" "ApplyFilter = no ; to apply a median filter before extraction\n" "\n" "# This is the detection method\n" "\n" "[SpectrumExtraction/method]\n" "Select = automatic ; automatic or manual\n" "SpectrumPosition = -1\n" "\n" "\n" "\n"); fprintf(f, "# ----- Section [Output]\n" "# The final result base name is declared here.\n" "# See the -o option to change that name when generating the\n" "# file (with the -g or --generate option).\n" "# The flags OutAveraged and OutDifference allow to have the\n" "# corresponding frames written on disk.\n" "\n" "[Output]\n" "OutAveraged = no ; write averaged frames on disk\n" "OutDifference = no ; write difference frames on disk\n\n") ; fprintf(f, "# ----- Section [PostProcessing]\n" "# Examples:\n" "# StartCommand = \"saoimage -fits %%s\" ;\n" "# StartCommand = \"rtd %%s\" ;\n" "# StartCommand = \"xv %%s\" ;\n" "\n" "[PostProcessing]\n" "Select = yes\n" "ProduceStatusReport = yes ; to produce global status report\n" "PlotSpectrum = no ; to plot the extracted spectrum\n" "\n" "[PostProcessing/StartViewer]\n" "Select = no ; to launch a viewer when finished\n" "StartCommand = ""saoimage -fits %%s""\n" "\n" "\n" "\n"); fclose(f); } /* end */ /*-------------------------------------------------------------------------*/ /** @name lw_spjitter_ini_parse @memo parse an ini file and return a filled blackboard structure @param ini_name name of ini file to parse @param name_i input ascii file @param name_o output basename @param name_c calibration ascii file (NULL if none) @return allocated blackboard structure @doc returns NULL in case of error. The returned blackboard must be freed by the caller. */ /*--------------------------------------------------------------------------*/ static lw_spjitter_bb * lw_spjitter_ini_parse( char * ini_name, char * name_i, char * name_o, char * name_c) { FILE * calib ; char frame_name[LINE_MAX+1] ; char frame_type[LINE_MAX+1] ; char line[LINE_MAX+1] ; lw_spjitter_bb * bb ; dictionary * sym ; int status ; charmatrix * charm ; char ** framelist ; int nval ; char * name ; char * cval ; int i, j ; /* Create output table to fill in */ bb = lw_spjitter_bb_create() ; if (bb==NULL) { e_error("allocating memory: aborting\n"); return NULL ; } /* FILL THE FIRST FIELDS OF THE BLACKBOARD */ strcpy(bb->frames_filelist, name_i) ; /* Read input char matrix */ charm = charmatrix_read(name_i); if (charm == NULL) { e_error("parsing input list [%s]", name_i); lw_spjitter_bb_destroy(bb) ; return NULL ; } /* Check input matrix */ nval = charm->ly ; for (i=0 ; ily ; i++) { /* Check file existence */ name = charmatrix_elem(charm, 0, i) ; if (file_exists(name) != 1 ) { e_warning("file [%s] declared in list does not exist", name) ; nval -- ; } } if (nval < 1) { e_error("no valid plane found in list [%s]", name_i) ; charmatrix_del(charm) ; lw_spjitter_bb_destroy(bb) ; return NULL ; } /* Allocate structures to go into the blackboard */ framelist = malloc(nval * sizeof(char*)); /* Browse through the charmatrix to get names */ i=0 ; for (j=0 ; jly ; j++) { name = charmatrix_elem(charm, 0, j) ; if (file_exists(name) == 1) framelist[i++] = strdup(name) ; } charmatrix_del(charm) ; bb->framelist = framelist ; bb->nframes = nval ; /* Output base name */ strcpy(bb->output_basename, name_o); /* Find out if named file really exists */ if (!file_exists(ini_name)) { e_error("cannot find ini file [%s]", ini_name); lw_spjitter_bb_destroy(bb) ; return NULL ; } /* Load ini file into symbolic table */ sym = iniparser_load(ini_name) ; if (sym == NULL) { e_error("parsing ini file [%s]: aborting", ini_name) ; lw_spjitter_bb_destroy(bb) ; return NULL ; } /* Parse symbolic table and fill up blackboard as necessary */ status = 0 ; parse_section_general(sym, bb, &status) ; parse_section_calibration(sym, bb, &status) ; parse_section_classification(sym, bb, &status) ; parse_section_distortion(sym, bb, &status) ; parse_section_wlcalib(sym, bb, &status) ; parse_section_combination(sym, bb, &status) ; parse_section_extraction(sym, bb, &status) ; parse_section_postproc(sym, bb, &status) ; parse_section_output(sym, bb, &status) ; /* Deallocate symbolic table */ iniparser_freedict(sym) ; if (status > 0) { e_error("%d error(s) in ini file [%s]", status, ini_name); lw_spjitter_bb_destroy(bb) ; bb = NULL ; } /* Write the spectrum length in the bb */ bb->spectrum_length = -1 ; cval = qfits_query_hdr(bb->framelist[0], "NAXIS1"); if (cval != NULL) bb->spectrum_length = (int)atoi(cval) ; /* Calibration files */ /* If a file name was given */ if (name_c != NULL) /* And it corresponds to a readable file */ if ((calib = fopen(name_c, "r")) != NULL) { /* Get one more line from the file */ while (fgets(line, LINE_MAX, calib) != NULL) /* If it is not a comment */ if (line[0] != '#') { nval = sscanf(line, "%s %s", frame_name, frame_type) ; /* If two string values can be read from the line */ if (nval == 2) { if (!strcmp(frame_type, "LW_ARC_COEF")) { strcpy(bb->cal_arc_name, frame_name) ; bb->cal_arc_active = 1 ; } if (!strcmp(frame_type, "LW_STAR_TRACE_COEF")) { strcpy(bb->cal_startrace_name, frame_name) ; bb->cal_startrace_active = 1 ; } if (!strcmp(frame_type, "LW_MASTER_SP_FLAT")) { strcpy(bb->cal_spflat_name, frame_name) ; bb->cal_spflat_active = 1 ; } } } fclose(calib); } return bb ; } /*-------------------------------------------------------------------------*/ /** @name parse_section_??? @memo update a blackboard from what can be found in the ini file @param sym symbolic table read from ini file @param bb blackboard @param status nb of errors counter @return void @doc all of these functions update a status integer to indicate if an error occurred, or leave it as it is if everything went Ok. parse_section_general() parse_section_calibration() parse_section_classification() parse_section_distortion() parse_section_wlcalib() parse_section_combination() parse_section_extraction() parse_section_postproc() parse_section_output() */ /*--------------------------------------------------------------------------*/ static void parse_section_general( dictionary * sym, lw_spjitter_bb * bb, int * status) { char * cval ; /* Check eclipse version */ cval = iniparser_getstr(sym, "eclipse:versionnumber") ; if (cval != NULL) if (strcmp(cval, get_eclipse_version())) { e_warning("this ini file produced by version %s", cval) ; e_warning("you are running version %s", get_eclipse_version()) ; } bb->framelist_hc_free = NULL ; bb->nb_hc_free = 0 ; bb->framelist_clean = NULL ; bb->nb_clean_frames = 0 ; bb->nb_hc1 = bb->nb_hc2 = 0 ; bb->hc1_list = bb->hc2_list = NULL ; return ; } static void parse_section_calibration( dictionary * sym, lw_spjitter_bb * bb, int * status) { char * cval ; cval = iniparser_getstr(sym, "calibration:arctable"); if (cval!=NULL) { if (strcmp(cval, "none")) { bb->cal_arc_active = 1 ; strcpy(bb->cal_arc_name, cval); } else { strcpy(bb->cal_arc_name, "none") ; } } cval = iniparser_getstr(sym, "calibration:startracetable"); if (cval!=NULL) { if (strcmp(cval, "none")) { bb->cal_startrace_active = 1 ; strcpy(bb->cal_startrace_name, cval); } else { strcpy(bb->cal_startrace_name, "none") ; } } cval = iniparser_getstr(sym, "calibration:masterspflattable"); if (cval!=NULL) { if (strcmp(cval, "none")) { bb->cal_spflat_active = 1 ; strcpy(bb->cal_spflat_name, cval); } else { strcpy(bb->cal_spflat_name, "none") ; } } return ; } static void parse_section_classification( dictionary * sym, lw_spjitter_bb * bb, int * status) { char * cval ; char * cval2 ; cval = iniparser_getstr(sym, "classification/source:select") ; if (cval==NULL) { e_warning("no source specified: switching to header classification"); bb->offsets_source = OFFSOURCE_HEADER ; cval2 = iniparser_getstr(sym, "classification/source:offsetkeyword") ; if (cval2==NULL) { e_error("no offset keyword provided") ; (*status)++ ; return ; } else { strcpy(bb->keyword_used_for_classif, cval2) ; } } else { if (!strcmp(cval, "header")) { bb->offsets_source = OFFSOURCE_HEADER ; cval2 = iniparser_getstr(sym, "classification/source:offsetkeyword") ; if (cval2==NULL) { e_error("no offset keyword provided") ; (*status)++ ; return ; } else { strcpy(bb->keyword_used_for_classif, cval2) ; } } else if (!strcmp(cval, "file")) { bb->offsets_source = OFFSOURCE_FILE ; cval2 = iniparser_getstr(sym, "classification/source:offsetfile"); if (cval2==NULL) { e_error("classification source is file"); e_error("but no file name provided"); (*status)++; } else if (!strcmp(cval2, "none")) { e_error("classification source is file"); e_error("but no file name provided"); (*status)++; } else { strcpy(bb->offsets_sourcefile, cval2); } } } return ; } static void parse_section_distortion( dictionary * sym, lw_spjitter_bb * bb, int * status) { int ival ; if (iniparser_getboolean(sym,"combination/distortion:select",0) == 1) { bb->distortion_active = 1 ; } if (iniparser_getboolean(sym, "combination/distortion:autodarksubtraction",0) == 1) { bb->auto_dark_subtraction = 1 ; } ival = iniparser_getint(sym, "combination/distortion:xmin", -1) ; if (ival < 0) { e_warning("illegal or no value for [Combination/Distortion]:XMin") ; e_warning("using default XMin [%d]", DISXMIN) ; ival = DISXMIN ; } bb->distor_xmin = ival ; ival = iniparser_getint(sym, "combination/distortion:ymin", -1) ; if (ival < 0) { e_warning("illegal or no value for [Combination/Distortion]:YMin") ; e_warning("using default YMin [%d]", DISYMIN) ; ival = DISYMIN ; } bb->distor_ymin = ival ; ival = iniparser_getint(sym, "combination/distortion:xmax", -1) ; if (ival < 0) { e_warning("illegal or no value for [Combination/Distortion]:XMax") ; e_warning("using default XMax [%d]", DISXMAX) ; ival = DISXMAX ; } bb->distor_xmax = ival ; ival = iniparser_getint(sym, "combination/distortion:ymax", -1) ; if (ival < 0) { e_warning("illegal or no value for [Combination/Distortion]:YMax") ; e_warning("using default YMax [%d]", DISYMAX) ; ival = DISYMAX ; } bb->distor_ymax = ival ; return ; } static void parse_section_wlcalib( dictionary * sym, lw_spjitter_bb * bb, int * status) { int ival ; char * cval ; cval = iniparser_getstr(sym, "wavelengthcalibration:wavecalarcfile"); if (cval!=NULL) { if (strcmp(cval, "none")) { bb->wavecal_arc_active = 1 ; strcpy(bb->wavecal_arcfile, cval); } else { strcpy(bb->wavecal_arcfile, "none") ; } } if (iniparser_getboolean(sym, "wavelengthcalibration:select",0) == 1) { bb->wavecal_active = 1 ; } ival = iniparser_getint(sym, "wavelengthcalibration:discardhigh", -1) ; if (ival < 0) { e_warning("illegal [WavelengthCalibration]:DiscardHigh") ; e_warning("using default DiscardHigh [%d]", WC_DISHIGH) ; ival = WC_DISHIGH ; } bb->wavecal_discard_hi = ival ; ival = iniparser_getint(sym, "wavelengthcalibration:discardlow", -1) ; if (ival < 0) { e_warning("illegal [WavelengthCalibration]:DiscardLow") ; e_warning("using default DiscardLow [%d]", WC_DISLOW) ; ival = WC_DISLOW ; } bb->wavecal_discard_lo = ival ; bb->wavecal_discard_le=iniparser_getint(sym, "wavelengthcalibration:discardleft", -1) ; bb->wavecal_discard_ri=iniparser_getint(sym, "wavelengthcalibration:discardright", -1) ; ival = iniparser_getint(sym, "wavelengthcalibration:wavecalnbcoeff", -1) ; if (ival <= 0) { e_warning("illegal [WavelengthCalibration]:WavecalNbCoeff - use 2") ; ival = 2 ; } bb->wavecal_nb_coeff = ival ; return ; } static void parse_section_combination( dictionary * sym, lw_spjitter_bb * bb, int * status) { double dval ; char * cval ; bb->divided_by_flat = -1 ; if (iniparser_getboolean(sym, "combination:circularshift", 0) == 1) { bb->circular_shift = 1 ; } /* Get the final combination method */ cval = iniparser_getstr(sym, "combination/finalcombmethod:select") ; if (cval==NULL) { e_warning("default final combination method used: [median]"); strcpy(bb->final_comb_method, "median"); } else { strcpy(bb->final_comb_method, cval); } /* Get the high rejected rate for averaging */ dval = iniparser_getdouble(sym, "combination/finalcombmethod:averagehirejection", -1.0) ; if (dval<= 0.0) { bb->average_hi_rejection = 0.0 ; } else { bb->average_hi_rejection = dval ; } /* Get the low rejected rate for averaging */ dval = iniparser_getdouble(sym, "combination/finalcombmethod:averagelorejection", -1.0) ; if (dval<= 0.0) { bb->average_lo_rejection = 0.0 ; } else { bb->average_lo_rejection = dval ; } return ; } static void parse_section_extraction( dictionary * sym, lw_spjitter_bb * bb, int * status) { int ival ; char * cval ; ival = iniparser_getint(sym, "spectrumextraction:badleft", -1) ; if (ival < 0) { e_warning("illegal [SpectrumExtraction]:BadLeft") ; e_warning("using default BadLeft [%d]", SPECTRACT_BADLEFT) ; ival = SPECTRACT_BADLEFT ; } bb->detect_bad_left = ival ; ival = iniparser_getint(sym, "spectrumextraction:badright", -1) ; if (ival < 0) { e_warning("illegal [SpectrumExtraction]:BadRight") ; e_warning("using default BadRight [%d]", SPECTRACT_BADRIGHT) ; ival = SPECTRACT_BADRIGHT ; } bb->detect_bad_right = ival ; ival = iniparser_getint(sym, "spectrumextraction:badtop", -1) ; if (ival < 0) { e_warning("illegal [SpectrumExtraction]:BadTop") ; e_warning("using default BadTop [%d]", SPECTRACT_BADTOP) ; ival = SPECTRACT_BADTOP ; } bb->detect_bad_top = ival ; ival = iniparser_getint(sym, "spectrumextraction:badbot", -1) ; if (ival < 0) { e_warning("illegal [SpectrumExtraction]:BadBot") ; e_warning("using default BadBot [%d]", SPECTRACT_BADBOT) ; ival = SPECTRACT_BADBOT ; } bb->detect_bad_bot = ival ; if (iniparser_getboolean(sym, "spectrumextraction:select",0) == 1) { bb->spectrum_extr_active = 1 ; } cval = iniparser_getstr(sym, "spectrumextraction/method:select"); if (cval!=NULL) { if (strcmp(cval, "manual")) { ival = iniparser_getint(sym, "spectrumextraction/method:spectrumposition", -1); if (ival < 0) ival = -1 ; bb->spectrum_position = ival ; } else if (strcmp(cval, "automatic")) { bb->spectrum_position = -1 ; } else { bb->spectrum_position = -1 ; e_warning("Detection method not recognized - use automatic method") ; } } ival = iniparser_getint(sym, "spectrumextraction:spectrumwidth", -1) ; if (ival < 0) { e_warning("illegal [SpectrumExtraction]:SpectrumWidth") ; e_warning("using SpectrumWidth = 10") ; ival = 10 ; } bb->spectrum_width = ival ; ival = iniparser_getint(sym, "spectrumextraction:resskyhiwidth", -1) ; if (ival < 0) ival = -1 ; bb->res_sky_hi_width = ival ; ival = iniparser_getint(sym, "spectrumextraction:resskylowidth", -1) ; if (ival < 0) ival = -1 ; bb->res_sky_lo_width = ival ; ival = iniparser_getint(sym, "spectrumextraction:resskyhidist", -1) ; if (ival < 0) ival = -1 ; bb->res_sky_hi_dist = ival ; ival = iniparser_getint(sym, "spectrumextraction:resskylodist", -1) ; if (ival < 0) ival = -1 ; bb->res_sky_lo_dist = ival ; if (iniparser_getboolean(sym, "spectrumextraction:applyfilter",0) == 1) { bb->apply_filter = 1 ; } bb->spectrum_detected = -1 ; bb->spectrum_extracted = -1 ; return ; } static void parse_section_postproc( dictionary * sym, lw_spjitter_bb * bb, int * status) { char * cval ; if (iniparser_getboolean(sym,"postprocessing:select",0)==1){ bb->postproc_active = 1 ; if (iniparser_getboolean(sym,"postprocessing/startviewer:select",0)==1) { bb->postproc_startviewer = 1 ; cval = iniparser_getstr(sym, "postprocessing/startviewer:startcommand") ; if (cval == NULL) { e_error("viewer start requested but no viewer provided"); bb->postproc_startviewer = 0 ; (*status)++ ; } else { strcpy(bb->postproc_viewer, cval); } } /* Plot the extracted spectrum? */ if (iniparser_getboolean(sym,"postprocessing:plotspectrum",0)==1){ bb->postproc_gnuplot = 1 ; } /* Produce final status report? */ if (iniparser_getboolean(sym,"postprocessing:producestatusreport",1)==1){ bb->postproc_statusreport = 1 ; } } return ; } static void parse_section_output( dictionary * sym, lw_spjitter_bb * bb, int * status) { if (iniparser_getboolean(sym, "output:outaveraged", 0)==1) { bb->output_averaged_frames = 1 ; } if (iniparser_getboolean(sym, "output:outdifference", 0)==1) { bb->output_difference_frames = 1 ; } return ; }