diff --git a/metis_dark_vlda.pro b/metis_dark_vlda.pro index 0a89a2866af7a1c1c0733c10151c5e1a5045a436..57b59880eba1a3c056a18e1b8346fc1aeacc2fb9 100644 --- a/metis_dark_vlda.pro +++ b/metis_dark_vlda.pro @@ -1,4 +1,4 @@ -function metis_dark_vlda, data, header, cal_pack, history = history +function metis_dark_vlda, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history bias = cal_pack.vl_channel.bias dark = cal_pack.vl_channel.dark @@ -6,24 +6,36 @@ function metis_dark_vlda, data, header, cal_pack, history = history dit = header.dit/1000.D0 ndit = header.ndit + nbin = header.nbin xposure = dit * ndit - nbin = sqrt(header.nbin) tsensor = header.tsensor obt_beg = header.obt_beg/1000.D0 ; check why obt_beg must be divided by 1000. + gain = cal_pack.vl_channel.radiometry.gain.value + ; WARN - temporary patch to handle local l1 fits files - ; if tsensor eq 0. then tsensor = -30. + if tsensor eq 0. then tsensor = -30. + ; check if one or more bias+dark matrices with same dit, binning, and tsensor exist for i = 0, n_elements(bias_dark) - 1 do begin - if bias_dark[i].dit eq dit and bias_dark[i].nbin eq nbin and abs(bias_dark[i].tsensor - tsensor) lt 5. then begin + if $ + bias_dark[i].dit eq dit and $ + bias_dark[i].nbin^2 eq nbin and $ + abs(bias_dark[i].tsensor - tsensor) lt 5. then begin + dark_file = bias_dark[i].file_name dark_image = float(readfits(cal_pack.path + dark_file, /silent)) - dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2 * ndit + dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin + dark_image = dark_image * ndit mask = where(data eq 0.) data = data - dark_image data[mask] = 0. + ; WARN - for this case the error matrix is not provided so the error is not computed + + if isa(error) then error += 0. + goto, jump endif endfor @@ -37,6 +49,7 @@ function metis_dark_vlda, data, header, cal_pack, history = history delta_obt = min(abs(obt_beg - bias_obt_available), j) i = where(bias.obt_beg eq bias_obt_available[j]) bias_image = float(readfits(cal_pack.path + bias[i].file_name, /silent)) + bias_error = float(readfits(cal_pack.path + bias[i].error, /silent)) dark_file = bias[i].file_name endif else begin journal, 'Error 02: VL applicable bias file not found.' @@ -53,6 +66,7 @@ function metis_dark_vlda, data, header, cal_pack, history = history delta_obt = min(abs(obt_beg - dark_obt_available), j) i = where(dark.obt_beg eq dark_obt_available[j]) dark_image = float(readfits(cal_pack.path + dark[i].file_name, /silent)) + dark_error = float(readfits(cal_pack.path + dark[i].error, /silent)) dark_file = [dark_file, dark[i].file_name] endif else begin journal, 'Error 03: VL applicable dark file not found.' @@ -60,14 +74,25 @@ function metis_dark_vlda, data, header, cal_pack, history = history exit, status = 3 endelse - bias_image = rebin(bias_image, header.naxis1, header.naxis2) * nbin^2 - dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin^2 + ; WARN - error for binning must be checked + + bias_nobin = bias_image + bias_image = rebin(bias_image, header.naxis1, header.naxis2) * nbin + bias_error = bias_image * sqrt(rebin((bias_error/bias_nobin)^2, header.naxis1, header.naxis2)/nbin) + + dark_nobin = dark_image + dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin + dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin) mask = where(data eq 0.) data = data - (bias_image * ndit + dark_image * xposure) data[mask] = 0. + if isa(error) then error += (data * gain + 2. * (bias_error * ndit)^2 + 2. * (dark_error * xposure)^2)/data^2 + jump: + s = where(data le 0., count) + if isa(quality_matrix) then quality_matrix[s] = 0 if ~ isa(history) then history = !null history = [history, 'Bias and dark-current corrections: ', ' ' + dark_file]