diff --git a/interpol_param.pro b/interpol_param.pro index 100da790e639e7be0243e693da9e08efb514a64e..88dcb5e228ba112ca7c82e0988af210c46070199 100644 --- a/interpol_param.pro +++ b/interpol_param.pro @@ -1,6 +1,6 @@ function interpol_param, table, par_name, date, empty_params = empty_params i = where(table.par_name eq par_name, count) - if count lt 1 then return, fix(32768) ; this should not be necessary + if count lt 1 then return, 0.0 ; this should not be necessary par_date = dblarr(count) par_val = fltarr(count) for j = 0, count - 1 do begin @@ -15,5 +15,5 @@ function interpol_param, table, par_name, date, empty_params = empty_params par_val[j] = float(table.eng_val[k]) endfor value = interpol(par_val, par_date, date_conv(date, 'JULIAN')) - return, value + if finite(value) then return, value else return, 0.0 end \ No newline at end of file diff --git a/json_write.pro b/json_write.pro new file mode 100644 index 0000000000000000000000000000000000000000..e5bfb459c80850ac46ebaf8732e99154bf362927 --- /dev/null +++ b/json_write.pro @@ -0,0 +1,5 @@ +pro json_write, struct, filename + openw, unit, filename, /get_lun + printf, unit, struct, /implied_print + free_lun, unit +end diff --git a/metis_l1_prep.pro b/metis_l1_prep.pro index a43a8fa7496cb5daefa99c6eefb201a5fcbe8d23..62af567827af0953cff0c5c485346e1c36d582fb 100644 --- a/metis_l1_prep.pro +++ b/metis_l1_prep.pro @@ -2,7 +2,7 @@ pro metis_l1_prep ; start the log - journal,'output/metis_l1_prep_log.txt' + journal, 'output/metis_l1_prep_log.txt' ; some definitions @@ -33,6 +33,8 @@ pro metis_l1_prep input = json_parse('input/contents.json', /toarray, /tostruct) + journal, 'File ' + input.file_name + ; prepare the spice kernels kernels = [ $ @@ -44,6 +46,10 @@ pro metis_l1_prep cspice_furnsh, kernels + journal, 'SPICE kernel files correctly loaded:' + journal, ' sclk kernel = ' + input.tsc_file_name + journal, ' leap seconds kernel = ' + input.tls_file_name + ; read l0 fits structure fits_info, input.file_name, n_ext = n_ext, extname = extname, /silent @@ -51,7 +57,7 @@ pro metis_l1_prep ; read the primary hdu ext_no = 0 - image = mrdfits(input.file_name, ext_no, primary_header, /silent) + data = mrdfits(input.file_name, ext_no, primary_header, /silent) ; if the data is not an image, read the data binary-table extension @@ -73,14 +79,9 @@ pro metis_l1_prep naxis = 2 naxis1 = metis_obj_size[datatype] naxis2 = naxis1 - bitpix = 16 - data_size = naxis1 * naxis2 - data_volume = (data_size * bitpix) / 8. endif else begin naxis = 0 bitpix = 8 - data_size = 0 - data_volume = 0 endelse ; if the data product is an accumulation matrix, look for the accumulation vector extension and read it if it exists @@ -92,6 +93,10 @@ pro metis_l1_prep endif else accumul_vector = !null endif + journal, 'L0 FITS file correctly read:' + journal, ' datatype = ' + string(datatype, format = '(I0)') + journal, ' size = ' + string(naxis1, format = '(I0)') + '×' + string(naxis2, format = '(I0)') + ; NOTE - the radialization extension is ignored ; read the planning data @@ -117,6 +122,8 @@ pro metis_l1_prep date_beg_string = strreplace(strreplace(strmid(date_beg, 0, 19), '-', ''), ':', '') + journal, 'UTC time of start acquisition = ' + date_beg + ; name of the fits file file_name_fields = strsplit(fxpar(primary_header, 'FILENAME'), '_', /extract) @@ -126,8 +133,8 @@ pro metis_l1_prep ; exposure times - dit = fxpar(metadata_extension_header, 'DIT') - telapse = obt_end - obt_beg + dit = float(fxpar(metadata_extension_header, 'DIT')) + telapse = float(obt_end - obt_beg) xposure = dit/1000. ; instrument keywords @@ -135,18 +142,18 @@ pro metis_l1_prep if datatype eq 0 or datatype eq 3 or datatype eq 5 or datatype eq 9 then begin filter = 'VL' - wavelnth = 0.0 - wavemin = 0.0 - wavemax = 0.0 + wavelnth = 610.0 + wavemin = 580.0 + wavemax = 640.0 waveband = 'Visible light 580-640 nm' endif else begin filter = 'UV' - wavelnth = 1215.67 - wavemin = 0.0 - wavemax = 0.0 + wavelnth = 121.6 + wavemin = 111.6 + wavemax = 131.6 waveband = 'H I Lyman-alpha 121.6 nm' endelse - detector = filter + 'DA' + detector = filter + 'D' telescope = 'SOLO/Metis/' + detector ; campaign keywords @@ -169,6 +176,27 @@ pro metis_l1_prep primary_header[i : *] $ ] + ; rebin the image if binning was applied during the acquisition + + if fxpar(primary_header, 'COMPR') then bin_type = fxpar(primary_header, 'B0_BIN') else bin_type = 0 + + fpfilt = fxpar(primary_header, 'VLFPFILT', count = count) + if count gt 0 and fpfilt eq 0 and bin_type eq 0 then bin_type = 1 + + bin_fact = 2. ^ bin_type + if bin_fact gt 1. then begin + data = rebin(data, naxis1 / bin_fact, naxis2 / bin_fact, /sample) + if ~ isa(comment) then comment = !null + comment = [comment, 'Image was rebinned according to the commanded binning factor.'] + + journal, 'Data was binned during acquisition:' + journal, ' binning = ' + string(bin_fact, format = '(I1)') + '×' + string(bin_fact, format = '(I1)') + journal, ' data was correctly rebinned:' + journal, ' new size = ' + string(naxis1 / bin_fact, format = '(I0)') + '×' + string(naxis1 / bin_fact, format = '(I0)') + endif else begin + journal, 'Data was not binned during acquistion.' + endelse + ; adjust the primary header (it is almost the same for all data product types) fxaddpar, primary_header, 'FILENAME', file_name @@ -202,7 +230,7 @@ pro metis_l1_prep fxaddpar, primary_header, 'OBS_ID', planning_data.obs_id, 'Unique identifier for the observation that is associated with the data acquisition', before = 'DATAMIN' fxaddpar, primary_header, 'TARGET', 'TBD', 'Taget as defined in the SOOP description', before = 'DATAMIN' fxaddpar, primary_header, 'BSCALE', 1.0, before = 'DATAMIN' - fxaddpar, primary_header, 'BZERO', 0, before = 'DATAMIN' + fxaddpar, primary_header, 'BZERO', 0.0, before = 'DATAMIN' fxaddpar, primary_header, 'BTYPE', strupcase(metis_datatype[datatype]), 'Science data object type', before = 'DATAMIN' fxaddpar, primary_header, 'BUNIT', 'DN', 'Units of physical value, after application of BSCALE and BZERO', before = 'DATAMIN' if datatype le 6 then begin @@ -211,9 +239,9 @@ pro metis_l1_prep fxaddpar, primary_header, 'PXEND1', naxis1, 'Last pixel that has been read out in dimension 1', before = 'COMPRESS' fxaddpar, primary_header, 'PXEND2', naxis2, 'Last pixel that has been read out in dimension 2', before = 'COMPRESS' endif - fxaddpar, primary_header, 'NBIN1', 1, 'Binning factor in the dimension 1', before = 'COMPRESS' - fxaddpar, primary_header, 'NBIN2', 1, 'Binning factor in the dimension 2', before = 'COMPRESS' - fxaddpar, primary_header, 'NBIN', 1, 'Product of all NBIN values above', before = 'COMPRESS' + fxaddpar, primary_header, 'NBIN1', bin_fact, 'Binning factor in the dimension 1', before = 'COMPRESS' + fxaddpar, primary_header, 'NBIN2', bin_fact, 'Binning factor in the dimension 2', before = 'COMPRESS' + fxaddpar, primary_header, 'NBIN', bin_fact * bin_fact, 'Product of all NBIN values above', before = 'COMPRESS' fxaddpar, primary_header, 'IDB_VERS', input.idb_version, '', before = 'HDR_VERS' fxaddpar, primary_header, 'INFO_URL', 'http://metis.oato.inaf.it', 'Link to more information on the instrument data', before = 'HISTORY' @@ -222,7 +250,6 @@ pro metis_l1_prep hk_table = json_parse(input.hk_file_name, /toarray, /tostruct) ; replace raw values with calibrated values in the primary header - ; WARN - must be done here? if datatype eq 0 or datatype eq 3 or datatype eq 5 then begin @@ -239,28 +266,76 @@ pro metis_l1_prep fxaddpar, primary_header, 'TSENSOR ', interpol_param(hk_table, 'NIT0E0E0', date_avg, empty_params = empty_params) fxaddpar, primary_header, 'PMPTEMP ', interpol_param(hk_table, 'NIT0L00D', date_avg, empty_params = empty_params) + + journal, 'Header keywords were calibrated using HK parameters:' + journal, ' TSENSOR = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)') + journal, ' PMPTEMP = ' + string(fxpar(primary_header, 'PMPTEMP'), format = '(F0)') endif if datatype eq 1 or datatype eq 4 or datatype eq 6 then begin fxaddpar, primary_header, 'HVU_SCR ', interpol_param(hk_table, 'NIT0E070', date_avg, empty_params = empty_params) fxaddpar, primary_header, 'HVU_MCP ', interpol_param(hk_table, 'NIT0E071', date_avg, empty_params = empty_params) fxaddpar, primary_header, 'TSENSOR ', interpol_param(hk_table, 'NIT0E050', date_avg, empty_params = empty_params) + + journal, 'Header keywords were calibrated using HK parameters:' + journal, ' HVU_SCR = ' + string(fxpar(primary_header, 'HVU_SCR'), format = '(F0)') + journal, ' HVU_MCP = ' + string(fxpar(primary_header, 'HVU_MCP'), format = '(F0)') + journal, ' TSENSOR = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)') endif if ~ isa(empty_params) then empty_params = '' - ; add keywords for file history + ; replace verbose keyword values + + key = fxpar(primary_header, 'MEASKIND', count = count) + if count gt 0 then begin + pol_id = fxpar(primary_header, 'POL_ID') + if key eq 0 and pol_id eq 0 then fxaddpar, primary_header, 'MEASKIND', 'Fixed pol.' + if key eq 0 and pol_id ne 0 then fxaddpar, primary_header, 'MEASKIND', 'pB' + if key eq 1 then fxaddpar, primary_header, 'MEASKIND', 'tB' + endif + + key = fxpar(primary_header, 'FRAMEMOD', count = count) + if count gt 0 then begin + if key eq 0 then fxaddpar, primary_header, 'FRAMEMOD', 'Single' + if key eq 1 then fxaddpar, primary_header, 'FRAMEMOD', 'Multiple' + endif + + key = fxpar(primary_header, 'VLFPFILT', count = count) + if count gt 0 then begin + pol_id = fxpar(primary_header, 'POL_ID') + if pol_id ne 0 then fxaddpar, primary_header, 'VLFPFILT', 'None' else begin + if key eq 0 then fxaddpar, primary_header, 'VLFPFILT', 'Binning' + if key eq 1 then fxaddpar, primary_header, 'VLFPFILT', 'Masking' + if key eq 2 then fxaddpar, primary_header, 'VLFPFILT', 'No filter' + endelse + endif + + key = fxpar(primary_header, 'PMPSTAB', count = count) + if count gt 0 then begin + if key eq 0 then fxaddpar, primary_header, 'PMPSTAB', 'Not stable' + if key eq 1 then fxaddpar, primary_header, 'PMPSTAB', 'Stable' + endif + + key = fxpar(primary_header, 'REF_ROWS', count = count) + if count gt 0 then begin + if key eq 0 then fxaddpar, primary_header, 'REF_ROWS', 'Not included' + if key eq 1 then fxaddpar, primary_header, 'REF_ROWS', 'Included' + endif - fxaddpar, primary_header, 'HISTORY', '' + keys = ['PRESUM', 'CME_OBS', 'SUNDISK', 'CR_SEP', 'SP_NOISE', 'COMPR', 'MASKING', 'RADIAL'] + foreach key_name, keys do begin + key = fxpar(primary_header, key_name, count = count) + if count gt 0 then begin + if key eq 0 then fxaddpar, primary_header, key_name, 'Disabled' + if key eq 1 then fxaddpar, primary_header, key_name, 'Enabled' + endif + endforeach ; remove unused keywords - sxdelpar, primary_header, 'CAD_BEG' - sxdelpar, primary_header, 'CAD_END' sxdelpar, primary_header, 'WIDTH' sxdelpar, primary_header, 'HEIGHT' - sxdelpar, primary_header, 'NB_IMG' - sxdelpar, primary_header, 'N' sxdelpar, primary_header, 'X_SIZE' sxdelpar, primary_header, 'Y_SIZE' sxdelpar, primary_header, 'Z_SIZE' @@ -268,14 +343,45 @@ pro metis_l1_prep sxdelpar, primary_header, 'N_BANDS' sxdelpar, primary_header, 'ORIG_X' sxdelpar, primary_header, 'ORIG_Y' - sxdelpar, primary_header, 'FIRSTROW' + + ; modify keywords for file history + + date = date_conv(systime(/julian, /utc), 'FITS') + + old_history = fxpar(primary_header, 'HISTORY') + sxdelpar, primary_header, 'HISTORY' + history = ['L1 FITS file created on ' + date, old_history] + for k = 0, n_elements(history) - 1 do $ + fxaddpar, primary_header, 'HISTORY', history[k] + + ; modify keywords for comments + + old_comment = fxpar(primary_header, 'COMMENT', count = count) + if count eq 0 then old_comment = '' + sxdelpar, primary_header, 'COMMENT' + if max(old_comment.startswith('Warning: correction of the OBT for pB acquisitions was applied.')) then begin + if ~ isa(comment) then comment = !null + comment = ['Correction of the OBT for pB acquisitions was applied at L0 level.', comment] + endif + if max(old_comment.startswith('Warning: the OBT of the acquisition is missing.')) then begin + if ~ isa(comment) then comment = !null + comment = ['Warning: the OBT of the acquisition is missing.',' The generation time of the first science packet was used', ' to compute OBT_BEG. Polarimetric keywords DAC*POL*',' may be wrong.', comment] + endif + if isa(comment) then begin + for k = 0, n_elements(comment) - 1 do $ + fxaddpar, primary_header, 'COMMENT', comment[k] + endif ; add checksum and datasum to the fits header ; WARN - should this be done at the end? i don't know fits_add_checksum, primary_header, image - mwrfits, image, out_file_name, primary_header, /no_comment, /create, /silent + mwrfits, data, out_file_name, primary_header, /no_comment, /create, /silent + + journal, 'Fits file created:' + journal, ' file name = ' + out_file_name + ; if applicable, save the data binary-table extension as it is @@ -309,6 +415,14 @@ pro metis_l1_prep mwrfits, hk_bin_table, out_file_name, hk_extension_header, /no_comment, /silent + journal, ' HK binary-table extension correctly added' + + ; unload the spice kernels + + cspice_unload, kernels + + journal, 'SPICE kernel files unloaded.' + ; write the auxiliary information file output = { $ @@ -318,15 +432,18 @@ pro metis_l1_prep empty_params : empty_params[uniq(empty_params)] $ } - openw, unit, 'output/contents.json', /get_lun - printf, unit, output, /implied_print - free_lun, unit + json_write, output, 'output/contents.json' - ; unload the spice kernels + journal, 'Output saved.' - cspice_unload, kernels + ;close the log - ; close the log + journal, 'Exiting without errors.' + journal + exit, status = 0 + error_handling: + journal, 'Errors occurred while processing.', /continue journal + exit, status = 1 end \ No newline at end of file