Newer
Older
; keyword defining if the detector reference frame must be used for the output
journal, 'output/metis_l1_prep_log.txt'
metis_datatype = list( $
'VL image', $ ; 0
'UV image', $ ; 1
'PCU accumulation matrix', $ ; 2
'VL temporal standard deviation', $ ; 3
'UV temporal standard deviation', $ ; 4
'VL cosmic-ray matrix', $ ; 5
'UV cosmic-ray matrix', $ ; 6
'PCU event list', $ ; 7
'PCU test event list', $ ; 8
'VL light-curve' $ ; 9
)
; read the auxiliary input file
input = json_parse('input/contents.json', /toarray, /tostruct)
journal, 'File ' + input.file_name
; load the spice kernels
load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
journal, 'SPICE kernel files correctly loaded:'
journal, ' SDK version = ' + kernel_version
; import the calibration package index file
cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
; include the calibration package path into the cal_pack structure
cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
journal, 'Calibration package correctly imported:'
journal, ' version = ' + cal_pack.version
journal, ' validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
data = mrdfits(input.file_name, 0, primary_header, /silent)
naxis = fxpar(primary_header, 'NAXIS')
naxis1 = fxpar(primary_header, 'NAXIS1', missing = 0)
naxis2 = fxpar(primary_header, 'NAXIS2', missing = 0)
; read the metadata extension
metadata_bin_table = mrdfits(input.file_name, 'metis metadata', metadata_extension_header, /silent)
; identify the data product and its nominal size
datatype = fxpar(metadata_extension_header, 'DATATYPE')
journal, 'L0 FITS file correctly read:'
journal, ' datatype = ' + string(datatype, format = '(I0)')
; if the data product is a light curve or a pcu-event list, read the data binary hk_table
if datatype gt 6 then data_bin_table = mrdfits(input.file_name, 1, data_extension_header, /silent) else data_bin_table = !null
; if the data product is an accumulation matrix, look for the accumulation vector extension and read it if it exists
if datatype eq 2 then $
if fxpar(metadata_extension_header, 'M') eq 256 then $
accumul_vector = mrdfits(input.file_name, 'accumulation vector', vector_extension_header, /silent) else accumul_vector = !null
; NOTE - the radialization extension is ignored
; read the planning data
planning_data = json_parse(input.planning_file_name, /toarray, /tostruct)
; definitions for the primary header
; version of the fits file
version = string(input.l1_version, format = '(I02)')
; creation and acquisition times in utc
date = date_conv(systime(/julian, /utc), 'FITS')
obt_beg = fxpar(primary_header, 'OBT_BEG')
obt_end = fxpar(primary_header, 'OBT_END')
obt_avg = (obt_beg + obt_end) / 2.0d
date_beg = solo_obt2utc(obt_beg)
date_end = solo_obt2utc(obt_end)
date_avg = solo_obt2utc(obt_avg)
date_beg_string = strmid(date_beg, 0, 19)
date_beg_string = date_beg_string.replace('-', '')
date_beg_string = date_beg_string.replace(':', '')
journal, 'Calculated UTC time of start acquisition:'
journal, ' time = ' + date_beg
; name of the fits file
file_name_fields = strsplit(fxpar(primary_header, 'FILENAME'), '_', /extract)
file_name = 'solo_L1_' + file_name_fields[2] + '_' + date_beg_string + '_V' + version + '.fits'
out_file_name = 'output/' + file_name
case datatype of
0: filter = 'VL'
1: filter = 'UV'
2: filter = 'UV'
3: filter = 'VL'
4: filter = 'UV'
5: filter = 'VL'
6: filter = 'UV'
7: filter = 'UV'
8: filter = 'UV'
9: filter = 'VL'
endcase
; choose the correct calibration-package structure based on the image filter
if filter eq 'VL' then channel = cal_pack.vl_channel else channel = cal_pack.uv_channel
; instrument/telescope/detector keywords
detector = filter + 'D'
telescope = 'SOLO/Metis/' + detector
wavelnth = channel.bandpass.value[1]
wavemin = channel.bandpass.value[0]
wavemax = channel.bandpass.value[2]
waveband = channel.name
if planning_data.obs_id ne 'none' then begin
obs_id_fields = strsplit(planning_data.obs_id, '_', /extract)
if obs_id_fields[2] eq '000' then sooptype = 'none' else sooptype = obs_id_fields[2]
if obs_id_fields[4] eq '000' then obs_type = 'none' else obs_type = obs_id_fields[4]
endif else begin
sooptype = 'none'
obs_type = 'none'
endelse
; build the fits file extensions
; join the metadata extension header to the primary header removing useless keywords
del_tags = list( $
'XTENSION', $
'BITPIX', $
'NAXIS', $
'NAXIS1', $
'NAXIS2', $
'PCOUNT', $
'GCOUNT', $
'TFIELDS', $
'TTYPE1', $
'TFORM1', $
'TTYPE2', $
'TFORM2', $
'EXTNAME', $
'CHECKSUM', $
'DATASUM', $
'END' $
)
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
foreach item, metadata_extension_header do begin
tag = strmid(item, 0, 8)
tag = tag.trim()
if ~isa(del_tags.where(tag)) then begin
value = fxpar(metadata_extension_header, tag, comment = comment)
fxaddpar, primary_header, tag, value, comment, before = 'CHECKSUM'
endif
endforeach
compr = fxpar(primary_header, 'COMPR', missing = -1)
xsize = fxpar(primary_header, 'X_SIZE', missing = -1)
if compr and xsize lt 0 then begin
bin_type = fxpar(primary_header, 'BIN_TYPE')
fxaddpar, primary_header, 'X_SIZE', 2048 / 2 ^ bin_type, before = 'CHECKSUM'
fxaddpar, primary_header, 'Y_SIZE', 2048 / 2 ^ bin_type, before = 'CHECKSUM'
fxaddpar, primary_header, 'Z_SIZE', 1, before = 'CHECKSUM'
fxaddpar, primary_header, 'P_BANDS', 0, before = 'CHECKSUM'
fxaddpar, primary_header, 'N_BANDS', 1, before = 'CHECKSUM'
fxaddpar, primary_header, 'ORIG_X', 2048, before = 'CHECKSUM'
fxaddpar, primary_header, 'ORIG_Y', 2048, before = 'CHECKSUM'
fxaddpar, primary_header, 'FIRSTROW', 0, before = 'CHECKSUM'
fxaddpar, primary_header, 'B0_BIN', bin_type, before = 'CHECKSUM'
fxaddpar, primary_header, 'B0_DQ', 0, before = 'CHECKSUM'
fxaddpar, primary_header, 'B0_STOP', 2048, before = 'CHECKSUM'
fxaddpar, primary_header, 'B1_BIN', 0, before = 'CHECKSUM'
fxaddpar, primary_header, 'B1_DQ', 0, before = 'CHECKSUM'
fxaddpar, primary_header, 'B1_STOP', 0, before = 'CHECKSUM'
fxaddpar, primary_header, 'B2_BIN', 0, before = 'CHECKSUM'
fxaddpar, primary_header, 'B2_DQ', 0, before = 'CHECKSUM'
fxaddpar, primary_header, 'B2_STOP', 0, before = 'CHECKSUM'
endif
; rebin the image if binning was applied during the acquisition and check for the data quality
; NOTE - this is done only for image data products
comment = !null
if datatype le 6 then begin
bin_type = fxpar(primary_header, 'B0_BIN', missing = 0)
; check for wrong binning type in fixed-polarization filtered images
pol_id = fxpar(primary_header, 'POL_ID', missing = -1)
vlfpfilt = fxpar(primary_header, 'VLFPFILT', missing = -1)
if vlfpfilt eq 0 and pol_id eq 0 and bin_type eq 0 then bin_type = 1
bin_fact = 2 ^ bin_type ; binning factor = 1, 2, or 4
if bin_fact gt 1 then begin
naxis1 = naxis1 / bin_fact
naxis2 = naxis2 / bin_fact
data = rebin(data, naxis1, naxis2)
; evaluate the quality matrix
quality_matrix = check_quality(data, cal_pack, filter)
; correct values including binning factor only for vl and uv images and accumulation matrices
if datatype le 2 then data = long(data) * bin_fact ^ 2
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 rebinned:'
journal, ' new size = ' + string(naxis1, format = '(I0)') + '×' + string(naxis2, format = '(I0)')
endif else begin
quality_matrix = check_quality(data, cal_pack, filter)
journal, 'Data was not binned during acquistion.'
endelse
; read the bad-pixel map
bad_pixel_map = readfits(cal_pack.path + channel.bad_pixel_map[0].file_name, /silent)
; adjust the map based on the presence of the reference rows
map_size = size(bad_pixel_map)
ref_rows = fxpar(primary_header, 'REF_ROWS', missing = -1)
if ref_rows ge 0 then first_row = ref_rows ? 2 : 0 else first_row = 0
bad_pixel_map = bad_pixel_map[*, first_row : first_row + map_size[1] - 1]
; rebin the bad-pixel map to the size of the image to be corrected
bad_pixel_map = rebin(bad_pixel_map, naxis1, naxis2) * bin_fact ^ 2
; correct the quality matrix according to the bad-pixel map
bad_pixels = where(bad_pixel_map gt 0, count)
if count gt 0 then quality_matrix[bad_pixels] = 0
endif else quality_matrix = !null
; correct for the effective exposure time
; NOTE - this is done only for vl and uv images and accumulation matrices
telapse = float(obt_end - obt_beg)
dit = fxpar(metadata_extension_header, 'DIT')
ndit = fxpar(metadata_extension_header, 'NDIT', missing = 1)
ndit1 = fxpar(metadata_extension_header, 'NDIT1', missing = 1)
ndit2 = fxpar(metadata_extension_header, 'NDIT2', missing = 1)
if datatype le 2 then begin
nsumexp = ndit * ndit1 * ndit2
xposure = dit / 1000.d0 * nsumexp
data = long(data) * nsumexp
if ~isa(comment) then comment = !null
comment = [comment, 'Image values were corrected for the total exposure time.']
journal, 'Image values were corrected for the total exposure time:'
journal, ' dit = ' + string(dit, format = '(I0)')
journal, ' nsumexp = ' + string(nsumexp, format = '(I0)')
endif else begin
nsumexp = 1
xposure = dit / 1000.d0
endelse
; adjust the primary header (it is almost the same for all data product types)
fxaddpar, primary_header, 'FILENAME', file_name
fxaddpar, primary_header, 'PARENT', file_basename(input.file_name), 'name of the parent file', before = 'APID'
fxaddpar, primary_header, 'DATE', date, 'date and time of FITS file creation', before = 'OBT_BEG'
fxaddpar, primary_header, 'DATE-OBS', date_beg, 'same as DATE-BEG', before = 'OBT_BEG'
fxaddpar, primary_header, 'DATE-BEG', date_beg, 'start time of observation', before = 'OBT_BEG'
fxaddpar, primary_header, 'DATE-AVG', date_avg, 'average time of observation', before = 'OBT_BEG'
fxaddpar, primary_header, 'DATE-END', date_end, 'end time of observation', before = 'OBT_BEG'
fxaddpar, primary_header, 'TIMESYS', 'UTC', 'system used for time keywords', before = 'OBT_BEG'
fxaddpar, primary_header, 'TIMRDER', 0.0, '[s] estimated random error in time values', before = 'OBT_BEG'
fxaddpar, primary_header, 'TIMSYER', 0.0, '[s] estimated systematic error in time values', before = 'OBT_BEG'
fxaddpar, primary_header, 'LEVEL', 'L1'
fxaddpar, primary_header, 'CREATOR', 'metis_l1_prep.pro'
fxaddpar, primary_header, 'VERS_SW', input.sw_version
fxaddpar, primary_header, 'VERS_CAL', cal_pack.version, 'version of the calibration package', after = 'VERS_SW'
fxaddpar, primary_header, 'VERSION', version
fxaddpar, primary_header, 'OBSRVTRY', 'Solar Orbiter', 'satellite name', before = 'INSTRUME'
fxaddpar, primary_header, 'TELESCOP', telescope, 'telescope that took the measurement', before = 'INSTRUME'
fxaddpar, primary_header, 'DETECTOR', detector, 'subunit/sensor', before = 'DATAMIN'
fxaddpar, primary_header, 'OBJECT', 'TBD', 'the use of the keyword OBJECT is [TBD]', before = 'DATAMIN'
fxaddpar, primary_header, 'OBS_MODE', planning_data.obs_mode, 'observation mode', before = 'DATAMIN'
fxaddpar, primary_header, 'OBS_TYPE', obs_type, 'encoded version of OBS_MODE', before = 'DATAMIN'
fxaddpar, primary_header, 'FILTER', filter, 'filter used to acquire this image', before = 'DATAMIN'
fxaddpar, primary_header, 'WAVELNTH', wavelnth, '[nm] characteristic wavelength of observation', before = 'DATAMIN'
fxaddpar, primary_header, 'WAVEMIN', wavemin, '[nm] min. bandpass wavelength', before = 'DATAMIN'
fxaddpar, primary_header, 'WAVEMAX', wavemax, '[nm] max. bandpass wavelength', before = 'DATAMIN'
fxaddpar, primary_header, 'WAVEBAND', waveband, 'bandpass description', before = 'DATAMIN'
fxaddpar, primary_header, 'XPOSURE', xposure, '[s] total effective exposure time', before = 'DATAMIN'
fxaddpar, primary_header, 'NSUMEXP', nsumexp, 'number of detector readouts summed together', before = 'DATAMIN'
fxaddpar, primary_header, 'TELAPSE', telapse, '[s] elapsed time during observation', before = 'DATAMIN'
fxaddpar, primary_header, 'SOOPNAME', planning_data.soop_name, 'name of the SOOP that the data belong to', before = 'DATAMIN'
fxaddpar, primary_header, 'SOOPTYPE', sooptype, 'campaign ID(s) that the data belong to', before = 'DATAMIN'
fxaddpar, primary_header, 'OBS_ID', planning_data.obs_id, 'unique ID of the individual observation', before = 'DATAMIN'
fxaddpar, primary_header, 'TARGET', 'TBD', 'type of target from planning', before = 'DATAMIN'
fxaddpar, primary_header, 'BSCALE', 1, 'ratio of physical to array value at 0 offset', before = 'DATAMIN'
fxaddpar, primary_header, 'BZERO', 0, 'physical value for the array value 0', before = 'DATAMIN'
fxaddpar, primary_header, 'BTYPE', metis_datatype[datatype], 'science data object type', before = 'DATAMIN'
fxaddpar, primary_header, 'BUNIT', 'DN', 'units of physical value', before = 'DATAMIN'
fxaddpar, primary_header, 'BLANK', 0, 'data value used to mark undefined pixels', after = 'BUNIT'
fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
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', before = 'HISTORY'
if datatype le 6 then begin
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, 'PXBEG1', 1, 'first pixel read out in dimension 1', before = 'COMPRESS'
fxaddpar, primary_header, 'PXBEG2', 1, 'first pixel read out in dimension 2', before = 'COMPRESS'
fxaddpar, primary_header, 'PXEND1', naxis1, 'last pixel read out in dimension 1', before = 'COMPRESS'
fxaddpar, primary_header, 'PXEND2', naxis2, 'last pixel read out in dimension 2', before = 'COMPRESS'
endif
; read the house-keeping telemetry
hk_table = json_parse(input.hk_file_name, /toarray, /tostruct)
; replace raw values with calibrated values in the primary header
empty_params = !null
if datatype eq 0 or datatype eq 3 or datatype eq 5 then begin
; patch to fix the lack of the keyword SEQ_NUM
pol_id = fxpar(primary_header, 'POL_ID', missing = 0)
if pol_id ge 1 and pol_id le 4 then begin
seq_num = fxpar(primary_header, 'SEQ_NUM', missing = 0)
if seq_num eq 0 then begin
obj_cnt = fxpar(primary_header, 'OBJ_CNT')
n_pol = fxpar(primary_header, 'N_POL', missing = 4)
seq_num = ((obj_cnt - 1) / n_pol + 1)
fxaddpar, primary_header, 'SEQ_NUM', seq_num, before = 'POL_ID'
endif
endif
; NOTE - DACPOL parameters are not calibrated since a calibration curve does not exist in the IDB. Their calibration in physical units (e.g., voltages or angles) should be done later
; fxaddpar, primary_header, 'DAC1POL1', interpol_param(hk_table, 'NIT0E061', date_avg, empty_params = empty_params)
; fxaddpar, primary_header, 'DAC2POL1', interpol_param(hk_table, 'NIT0E062', date_avg, empty_params = empty_params)
; fxaddpar, primary_header, 'DAC1POL2', interpol_param(hk_table, 'NIT0E064', date_avg, empty_params = empty_params)
; fxaddpar, primary_header, 'DAC2POL2', interpol_param(hk_table, 'NIT0E065', date_avg, empty_params = empty_params)
; fxaddpar, primary_header, 'DAC1POL3', interpol_param(hk_table, 'NIT0E067', date_avg, empty_params = empty_params)
; fxaddpar, primary_header, 'DAC2POL3', interpol_param(hk_table, 'NIT0E068', date_avg, empty_params = empty_params)
; fxaddpar, primary_header, 'DAC1POL4', interpol_param(hk_table, 'NIT0E06A', date_avg, empty_params = empty_params)
; fxaddpar, primary_header, 'DAC2POL4', interpol_param(hk_table, 'NIT0E06B', date_avg, empty_params = empty_params)
fxaddpar, primary_header, 'TSENSOR', interpol_param(hk_table, 'NIT0E0E0', date_avg, empty_params = empty_params), '[degC] VLDA temperature'
fxaddpar, primary_header, 'PMPTEMP', interpol_param(hk_table, 'NIT0L00D', date_avg, empty_params = empty_params), '[degC] PMP temperature'
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), '[raw] HVU Screen commanded voltage'
fxaddpar, primary_header, 'HVU_MCP', interpol_param(hk_table, 'NIT0E071', date_avg, empty_params = empty_params), '[raw] HVU MCP commanded voltage'
fxaddpar, primary_header, 'HV_SCR_V', interpol_param(hk_table, 'NIT0E0B7', date_avg, empty_params = empty_params), '[V] HVU Screen + MCP read voltage', after = 'HVU_MCP'
fxaddpar, primary_header, 'HV_MCP_V', interpol_param(hk_table, 'NIT0E0B6', date_avg, empty_params = empty_params), '[V] HVU MCP read voltage', after = 'HV_SCR_V'
fxaddpar, primary_header, 'HV_MCP_I', interpol_param(hk_table, 'NIT0E0BF', date_avg, empty_params = empty_params), '[uA] HVU MCP current', after = 'HV_MCP_V'
fxaddpar, primary_header, 'HV_TEMP', interpol_param(hk_table, 'NIT0E0B5', date_avg, empty_params = empty_params), '[degC] HVU temperature', after = 'HV_MCP_I'
fxaddpar, primary_header, 'TSENSOR', interpol_param(hk_table, 'NIT0E050', date_avg, empty_params = empty_params), '[degC] UVDA temperature'
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, ' HV_SCR_V = ' + string(fxpar(primary_header, 'HV_SCR_V'), format = '(F0)')
journal, ' HV_MCP_V = ' + string(fxpar(primary_header, 'HV_MCP_V'), format = '(F0)')
journal, ' HV_MCP_I = ' + string(fxpar(primary_header, 'HV_MCP_I'), format = '(F0)')
journal, ' HV_TEMP = ' + string(fxpar(primary_header, 'HV_TEMP'), format = '(F0)')
journal, ' TSENSOR = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)')
endif
if empty_params eq !null then empty_params = ''
; replace with text 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', 'Not applicable' 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
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
; convert the fits header into an idl structure
header = fits_hdr2struct(primary_header)
; append wcs keywords
if datatype le 6 then begin
wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
endif
; append solar keywords
ephemeris = solo_get_ephemeris(header, cal_pack)
foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
fxaddpar, primary_header, 'HISTORY', 'WCS and solar ephemeris:'
fxaddpar, primary_header, 'HISTORY', ' SKD version = ' + kernel_version
; remove redundant and misleading keywords
sxdelpar, primary_header, 'WIDTH'
sxdelpar, primary_header, 'HEIGHT'
sxdelpar, primary_header, 'X_SIZE'
sxdelpar, primary_header, 'Y_SIZE'
sxdelpar, primary_header, 'Z_SIZE'
sxdelpar, primary_header, 'P_BANDS'
sxdelpar, primary_header, 'N_BANDS'
sxdelpar, primary_header, 'ORIG_X'
sxdelpar, primary_header, 'ORIG_Y'
; modify keywords for file history
fxaddpar, primary_header, 'HISTORY', 'L1 FITS file created on ' + date
; 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
; rectify the image if requested
if not ref_detector then begin
data = metis_rectify(data, filter)
quality_matrix = metis_rectify(quality_matrix, filter)
endif
; add checksum and datasum to the fits header and save
fits_add_checksum, primary_header, data
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
if isa(data_bin_table) then begin
if datatype lt 9 then fits_add_checksum, data_extension_header, data_bin_table
mwrfits, data_bin_table, out_file_name, data_extension_header, /no_comment, /silent
endif
; save the quality matrix
if datatype le 6 and isa(quality_matrix) then begin
quality_matrix_header = !null
fxaddpar, quality_matrix_header, 'PCOUNT', 0, 'parameter count'
fxaddpar, quality_matrix_header, 'GCOUNT', 1, 'group count'
fxaddpar, quality_matrix_header, 'EXTNAME', 'Quality matrix', 'extension name'
fits_add_checksum, quality_matrix_header, quality_matrix
mwrfits, quality_matrix, out_file_name, quality_matrix_header, /no_comment, /silent
journal, 'Quality-matrix extension correctly added.'
endif
; build the telemetry extension
hk_extension_header = !null
fxaddpar, hk_extension_header, 'PCOUNT', 0, 'parameter count'
fxaddpar, hk_extension_header, 'GCOUNT', 1, 'group count'
fxaddpar, hk_extension_header, 'EXTNAME', 'House-keeping', 'extension name'
hk_bin_table = make_bin_table(hk_table)
mwrfits, hk_bin_table, out_file_name, hk_extension_header, /no_comment, /silent
journal, 'HK binary-table extension correctly added.'
; fix some header keywords
fix_fits_header, out_file_name
; unload the spice kernels
load_spice_kernels, kernel_list = kernel_list, /unload
journal, 'SPICE kernel files unloaded.'
; write the auxiliary information file
output = { $
file_name: out_file_name, $
l0_file_name: input.file_name, $
log_file_name: 'output/metis_l1_prep_log.txt', $
empty_params: empty_params[uniq(empty_params)] $
}