From 1c47a0c4f9fe449c6cec31158e7284caa233ce59 Mon Sep 17 00:00:00 2001
From: Roberto <roberto.susino@inaf.it>
Date: Mon, 13 Nov 2023 09:21:55 +0100
Subject: [PATCH] Manage rad. calibration of UV images with OBSW 2.2
---
metis_dark_uvda_new.pro | 165 ++++++++++++++++++++++++++++++++++++++++
metis_l2_prep_uv.pro | 6 +-
metis_rad_cal.pro | 1 -
3 files changed, 170 insertions(+), 2 deletions(-)
create mode 100644 metis_dark_uvda_new.pro
diff --git a/metis_dark_uvda_new.pro b/metis_dark_uvda_new.pro
new file mode 100644
index 0000000..6a815ad
--- /dev/null
+++ b/metis_dark_uvda_new.pro
@@ -0,0 +1,165 @@
+function m_angle_uv, matrix
+
+ a = mean(matrix[8:43, 8:43])
+ b = mean(matrix[980:1015, 8:43])
+ c = mean(matrix[8:43, 980:1015])
+ d = mean(matrix[980:1015, 980:1015])
+
+ return, [a, b, c, d]
+end
+
+function m_center_uv, matrix
+
+ a = mean(matrix[452:571, 460:579])
+
+ return, a
+end
+
+function correct_dark, data, dark, offset
+
+ dark_bin = rebin(dark, 1024, 1024, /sample)
+ data_bin = rebin(data, 1024, 1024, /sample)
+
+ x = [m_angle_uv(data_bin)]
+ y = [m_angle_uv(dark_bin)]
+
+ offset = (y[0] + y[2])/2. - (x[0] + x[2])/2.
+
+ dark_corr = dark - offset
+
+ return, dark_corr
+end
+
+function metis_dark_uvda_new, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history
+
+ dark = cal_pack.uv_channel.dark
+
+ dit = header.dit
+ nbin = header.nbin
+ ndit1 = header.ndit1
+ ndit2 = header.ndit2
+ obj_cnt = header.obj_cnt
+ obt_beg = header.obt_beg
+ masking = header.masking
+ tsensor = header.tsensor
+
+ journal, 'Dark-current correction:'
+ journal, ' dit = ' + string(dit, format = '(I0)') + ' ms'
+ journal, ' ndit1 = ' + string(ndit1, format = '(I0)')
+ journal, ' ndit2 = ' + string(ndit2, format = '(I0)')
+ journal, ' obj_cnt = ' + string(obj_cnt, format = '(I0)')
+ journal, ' tsensor = ' + string(tsensor, format = '(F0)') + ' degC'
+
+ flag_use_extrapolated = boolean(0)
+ flag_use_notfullset = boolean(1)
+ flag_correct_dark = boolean(1)
+
+ obt_available = !null
+
+ ; check if one or more dark matrices with same parameters exist
+ for i = 0, n_elements(dark) - 1 do begin
+ if $
+ (dark[i].extrapol eq 'False') and $
+ (dark[i].dit eq dit) and $
+ (abs(dark[i].tsensor - tsensor) lt 5) and $
+ (flag_use_notfullset or dark[i].fullset) then $
+ obt_available = [obt_available, dark[i].obt_beg]
+ endfor
+
+ ; if not, search for extrapolated dark matrices with same parameters
+ if ~ isa(obt_available) and flag_use_extrapolated then begin
+ for i = 0, n_elements(dark) - 1 do begin
+ if $
+ (dark[i].extrapol eq 'True') and $
+ (dark[i].dit eq dit) and $
+ (abs(dark[i].tsensor - tsensor) lt 5) and $
+ (flag_use_notfullset or dark[i].fullset) then $
+ obt_available = [obt_available, dark[i].obt_beg]
+ endfor
+ endif
+
+ ; if not, search for dark matrices with no constraint on tsensor
+ if ~ isa(obt_available) then begin
+ for i = 0, n_elements(dark) - 1 do begin
+ if $
+ (dark[i].extrapol eq 'False') and $
+ (dark[i].dit eq dit) and $
+ (flag_use_notfullset or dark[i].fullset) then $
+ obt_available = [obt_available, dark[i].obt_beg]
+ endfor
+ endif
+
+ ; if not, exit with error code
+ if ~ isa(obt_available) then begin
+ journal, 'Error 02: UV applicable dark file not found.'
+ journal
+ exit, status = 2
+ endif
+
+ ; if yes, chose the temporally closest
+ delta = min(abs(obt_beg - obt_available), j)
+ i = where(dark.obt_beg eq obt_available[j])
+
+ transient = where(dark[i].filename.cnt_1 ne '', n)
+
+ ; select the dark matrix corresponding to the correct transient phase
+ if obj_cnt le n then $
+ dark_file = dark[i].filename.cnt_1[obj_cnt - 1] else $
+ dark_file = dark[i].filename.cnt_2
+
+ dark_image = float(readfits(cal_pack.path + dark_file, /silent))
+
+ ; rebin the dark and take into account the exposure time correction of l1 images
+ dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin
+ dark_image = dark_image * ndit1 * ndit2
+
+ ; apply correction if required
+ if flag_correct_dark and masking.contains('disabled', /fold) then $
+ dark_image = correct_dark(data, dark_image, offset) else $
+ offset = 0
+
+ offset = offset/(nbin * ndit1 * ndit2)
+
+ mask = where(data eq 0., count)
+ data = data - dark_image
+
+ ; set values in masked areas equal to zero
+ if count gt 0 then data[mask] = 0.
+
+ ; read the dark error matrix
+ dark_error_file = dark[i].error
+
+ if file_test(cal_pack.path + dark_error_file) then $
+ dark_error = float(readfits(cal_pack.path + dark_error_file, hdr, /silent)) else $
+ dark_error = fltarr(1024, 1024)
+
+ ; rebin the error matrix and take into account the exposure time correction of l1 images
+ dark_error = sqrt(rebin(dark_error^2, header.naxis1, header.naxis2) * nbin)/nbin
+ dark_error = dark_error * ndit1 * ndit2
+
+ if not isa(error) then error = 0.
+
+ error += ((data > 0)/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2
+
+ ; set error values in masked areas equal to zero
+ if count gt 0 then error[mask] = 0.
+
+ jump:
+ s = where(data le 0., count)
+ if isa(quality_matrix) and count gt 0 then quality_matrix[s] = 0
+
+ if ~ isa(history) then history = !null
+ history = [history, $
+ 'Dark-current correction: ', $
+ ' ' + dark_file, $
+ ' extrapolated dark = ' + dark[i].extrapol.tolower(), $
+ ' corrected dark = ' + json_serialize(flag_correct_dark), $
+ ' offset = ' + string(offset, format = '(f0.1, " DN")')]
+
+ journal, ' dark file = ' + dark_file
+ journal, ' extrapolated dark = ' + dark[i].extrapol.tolower()
+ journal, ' corrected dark = ' + json_serialize(flag_correct_dark)
+ journal, ' offset = ' + string(offset, format = '(f0.1, " DN")')
+
+ return, data
+end
diff --git a/metis_l2_prep_uv.pro b/metis_l2_prep_uv.pro
index bc2457b..387bb4d 100644
--- a/metis_l2_prep_uv.pro
+++ b/metis_l2_prep_uv.pro
@@ -67,7 +67,11 @@ pro metis_l2_prep_uv
; header.nsumexp = header.ndit1 * header.ndit2
; ====================================
- data = metis_dark_uvda(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ if fix(header.hdr_vers) le 4 then begin
+ data = metis_dark_uvda(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ endif else begin
+ data = metis_dark_uvda_new(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
+ endelse
; ====================================
; for data already subtracted of dark
diff --git a/metis_rad_cal.pro b/metis_rad_cal.pro
index 26cb407..017b1b4 100644
--- a/metis_rad_cal.pro
+++ b/metis_rad_cal.pro
@@ -33,7 +33,6 @@ function metis_rad_cal, data, header, cal_pack, polarimetric = polarimetric, nor
eff_file = radiometry.opt_efficiency.file_name
eff = float(readfits(cal_pack.path + eff_file, /silent))
-
eff = rebin(eff, header.naxis1, header.naxis2)
eff_error_file = radiometry.opt_efficiency.error
--
GitLab