Improve history handling in polariz. pipeline

    ; calibration block

    data = !null

    history = !null

    abs_error = !null

    data_header = !null

    data_subdark = !null

    for k = 0, 3 do begin

        ; read the input image

        image = mrdfits(input.file_name[k], 0, primary_header, /silent)

        image = mrdfits(input.file_name[k], 0, header, /silent)

        ; patch to fix the lack of the keyword SEQ_NUM

        ; fix the lack of the keyword SEQ_NUM

        seq_num = fxpar(primary_header, 'SEQ_NUM', missing = 0)

        seq_num = fxpar(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)

            obj_cnt = fxpar(header, 'OBJ_CNT')

            n_pol = fxpar(header, 'N_POL', missing = 4)

            seq_num = ((obj_cnt - 1) / n_pol + 1)

            fxaddpar, primary_header, 'SEQ_NUM', seq_num, before = 'POL_ID'

            fxaddpar, header, 'SEQ_NUM', seq_num, before = 'POL_ID'

        endif

        if k eq 0 then comment = fxpar(primary_header, 'COMMENT') else begin

            sxdelpar, primary_header, 'COMMENT'

            foreach line, comment do $

                fxaddpar, primary_header, 'COMMENT', line

        endelse

        if k eq 0 then primary_header = header

        header = fits_hdr2struct(primary_header)

        ; clean header from comments and history

        sxdelpar, header, 'COMMENT'

        sxdelpar, header, 'HISTORY'

        ; convert header to idl structure

        header = fits_hdr2struct(header)

        ; ======================================================================

        ; for old l1 data

        ; apply dark correction to compute stokes i and total brightness

        tb_history = !null

        image_subdark = metis_dark_vlda(image, header, error = image_error, quality_matrix = image_quality_matrix, cal_pack, history = tb_history)

        image_subdark = metis_dark_vlda(image, header, error = image_error, quality_matrix = image_quality_matrix, cal_pack, history = k eq 0 ? history : !null)

        ; ======================================================================

        ; for data already subtracted of dark

    demod_history = 'Demodulation performed for angles ' + string(angles, format = '(3(f0.1, ", "), f0.1)') + ' deg'

    tb_history = [tb_history, demod_history]

    pb_history = demod_history

    i = reform(stokes_subdark[*, *, 0])

    i_abs_error = reform(stokes_abs_error[*, *, 0])

    q = reform(stokes[*, *, 1])

    ; compute the tb from the dark-subtracted stokes i and apply other calibrations

    cal_history = !null

    journal, 'Calibrating total brightness...'

    tb_image = i

    tb_error = (i_abs_error / i) ^ 2

    tb_image = metis_flat_field(tb_image, header, error = tb_error, quality_matrix = quality_matrix, cal_pack, history = tb_history)

    tb_image = metis_vignetting(tb_image, header, error = tb_error, quality_matrix = quality_matrix, cal_pack, history = tb_history)

    tb_image = metis_rad_cal(tb_image, header, error = tb_error, quality_matrix = quality_matrix, cal_pack, history = tb_history, /polarimetric)

    tb_image = metis_flat_field(tb_image, header, error = tb_error, quality_matrix = quality_matrix, cal_pack, history = cal_history)

    tb_image = metis_vignetting(tb_image, header, error = tb_error, quality_matrix = quality_matrix, cal_pack, history = cal_history)

    tb_image = metis_rad_cal(tb_image, header, error = tb_error, quality_matrix = quality_matrix, cal_pack, history = cal_history, /polarimetric)

    ; compute the pb from the stokes q and u and apply other calibrations

    pb_image = sqrt(q ^ 2 + u ^ 2)

    pb_error = (q ^ 2 * q_abs_error ^ 2 + u ^ 2 * u_abs_error ^ 2) / pb_image ^ 4

    pb_image = metis_flat_field(pb_image, header, error = pb_error, cal_pack, history = pb_history)

    pb_image = metis_vignetting(pb_image, header, error = pb_error, cal_pack, history = pb_history)

    pb_image = metis_rad_cal(pb_image, header, error = pb_error, cal_pack, history = pb_history, /polarimetric)

    pb_image = metis_flat_field(pb_image, header, error = pb_error, cal_pack)

    pb_image = metis_vignetting(pb_image, header, error = pb_error, cal_pack)

    pb_image = metis_rad_cal(pb_image, header, error = pb_error, cal_pack, /polarimetric)

    ; ==========================================================================

    ; for simple radiometric calibration

    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'

    wcs_history = ['Update WCS and solar ephemeris:', '  SKD version = ' + kernel_version]

    ; append solar ephemeris keywords

    ephemeris = solo_get_ephemeris(header, cal_pack)

    foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'

    ; update the comment and history keywords

    ; update the comment keyword

    fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'

    fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'

    history = ['Update WCS and solar ephemeris:', '  SKD version = ' + kernel_version]

    tb_history = [tb_history, history]

    pb_history = [pb_history, history]

    ; delete useless keywords

    sxdelpar, primary_header, 'BLANK'

    ; add the history keyword

    pb_history = [demod_history, cal_history, wcs_history]

    for k = 0, n_elements(pb_history) - 1 do fxaddpar, primary_pb_header, 'HISTORY', pb_history[k]

    fxaddpar, primary_pb_header, 'HISTORY', 'L2 FITS file created on ' + date

    ; write the fits file

    fxaddpar, primary_tb_header, 'DATAMAX', max(tb_image, /nan)

    ; add the history keyword

    tb_history = [history, demod_history, cal_history, wcs_history]

    for k = 0, n_elements(tb_history) - 1 do fxaddpar, primary_tb_header, 'HISTORY', tb_history[k]

    fxaddpar, primary_tb_header, 'HISTORY', 'L2 FITS file created on ' + date

    ; write the fits file

    ; add the history keyword

    for k = 0, n_elements(pb_history) - 1 do fxaddpar, primary_polangle_header, 'HISTORY', pb_history[k]

    fxaddpar, primary_polangle_header, 'HISTORY', 'L2 FITS file created on ' + date

    ; write the fits file

    ; add the history keyword

    for k = 0, n_elements(tb_history) - 1 do fxaddpar, primary_stokes_header, 'HISTORY', tb_history[k]

    fxaddpar, primary_stokes_header, 'HISTORY', 'L2 FITS file created on ' + date

    ; write the fits file