Manage rad. calibration of UV images with OBSW 2.2

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

		; header.nsumexp = header.ndit1 * header.ndit2

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

		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

		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