Update UVDA dark selection procedure

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.

	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.

	if tsensor eq 0. then tsensor = -25.

				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