diff --git a/metis_dark_uvda.pro b/metis_dark_uvda.pro index b62e7490ab1e14a8ff934bae7eb7b58d234242bc..21d72510f9a22f984d41afb2810aefcd54930044 100644 --- a/metis_dark_uvda.pro +++ b/metis_dark_uvda.pro @@ -1,92 +1,97 @@ -function metis_dark_uvda, data, header, cal_pack, history = history +function metis_dark_uvda, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history dark = cal_pack.uv_channel.dark - dit = header.dit/1000.D0 + dit = header.dit + nbin = header.nbin ndit1 = header.ndit1 ndit2 = header.ndit2 - nbin = sqrt(header.nbin) - tsensor = header.tsensor obj_cnt = header.obj_cnt obt_beg = header.obt_beg - + masking = header.masking + tsensor = header.tsensor + ; WARN - temporary patch to handle local l1 fits files - ; if tsensor eq 0. then tsensor = -25. + if tsensor eq 0. then tsensor = -25. - for i = 0, n_elements(dark) - 1 do begin - if (dark[i].dit eq dit) and (dark[i].ndit1 eq ndit1) and (dark[i].ndit2 eq ndit2) and (abs(dark[i].tsensor - tsensor) lt 5.) and (abs(obt_beg - dark[i].obt_beg) lt (3600. * 24. * 10.)) then begin - if obj_cnt eq 1 then begin - dark_file = dark[i].file_name.cnt_1 - dark_image = float(readfits(cal_pack.path + dark_file, /silent)) - dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2 - - mask = where(data eq 0., count) - data = data - dark_image * ndit1 * ndit2 - data[mask] = 0. - - goto, jump - endif else begin - dark_file = dark[i].file_name.cnt_2 - dark_image = float(readfits(cal_pack.path + dark_file, /silent)) - dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2 - - mask = where(data eq 0., count) - data = data - dark_image * ndit1 * ndit2 - data[mask] = 0. - - goto, jump - endelse - endif - endfor + flag_use_extrapolated = 1 + flag_use_notfullset = 1 + flag_normalize_dark = 0 obt_available = !null + + ; check if one or more dark matrices with same dit, ndit1, ndit2, and tsensor exist for i = 0, n_elements(dark) - 1 do begin - if (dark[i].dit eq dit) and (abs(dark[i].tsensor - tsensor) lt 5.) then obt_available = [obt_available, dark[i].obt_beg] + if $ + (dark[i].dit eq dit) and $ + (abs(dark[i].tsensor - tsensor) lt 5) and $ + (dark[i].ndit1 eq ndit1) and (dark[i].ndit2 eq ndit2) and $ + flag_use_notfullset and flag_use_extrapolated then $ + obt_available = [obt_available, dark[i].obt_beg] endfor - if not isa(obt_available) then begin + ; if not, exit with error code + if ~ isa(obt_available) then begin journal, 'Error 12: UV applicable dark file not found.' journal exit, status = 12 endif + ; if yes, chose the temporally closest delta_obt = min(abs(obt_beg - obt_available), j) i = where(dark.obt_beg eq obt_available[j]) - if obj_cnt eq 1 then begin - if (ndit1 eq dark[i].ndit1) and (ndit2 eq dark[i].ndit2) then begin - dark_file = dark[i].file_name.cnt_1 - dark_image = float(readfits(cal_pack.path + dark_file, /silent)) - endif else begin - q = ndit2/dark[i].ndit2 - dark_file = dark[i].file_name.cnt_1 + ' - ' + dark[i].file_name.cnt_2 - dark_1 = float(readfits(cal_pack.path + dark[i].file_name.cnt_1, /silent)) - dark_2 = float(readfits(cal_pack.path + dark[i].file_name.cnt_2, /silent)) - dark_image = (dark_1 + (q - 1) * dark_2)/q - endelse - endif else begin - dark_file = dark[i].file_name.cnt_2 + 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 + + ; apply normalization if required + if flag_normalize_dark and masking.contains('disabled', /fold) then $ + dark_image = correzione_uvda(data, cal_pack.path + dark_file) else $ dark_image = float(readfits(cal_pack.path + dark_file, /silent)) - endelse - dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2 + dark_nobin = dark_image + + ; 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 mask = where(data eq 0., count) - data = data - dark_image * ndit1 * ndit2 - data[mask] = 0. + 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(dark_error_file) then $ + dark_error = float(readfits(cal_pack.path + dark_error_file, /silent)) else $ + dark_error = fltarr(1024, 1024) + + ; rebin it and take into account the exposure time correction + dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin) + + ; calculate the total uncertainty + if isa(error) then error += (data/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2 jump: - + s = where(data le 0., count) + quality_matrix[s] = 0 + if ~ isa(history) then history = !null history = [history, 'Dark-current correction: ', ' ' + dark_file] journal, 'Dark-current correction:' - journal, ' dit = ' + string(dit, format = '(F0)') + ' s' + 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' journal, ' dark file = ' + dark_file + journal, ' extrapolated dark = ' + dark[i].extrapol.tolower() return, data end