Commit a33897b7 authored by Roberto Susino's avatar Roberto Susino
Browse files

Add error calculation in vignetting corr. procedure

parent 89cc7f42
Loading
Loading
Loading
Loading
+19 −3
Original line number Diff line number Diff line
function metis_vignetting, data, header, cal_pack, history = history
function metis_vignetting, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history

	journal, 'Vignetting correction:'

@@ -6,6 +6,7 @@ function metis_vignetting, data, header, cal_pack, history = history
	if header.filter.contains('UV', /fold) then vig_info = cal_pack.uv_channel.vignetting

	vig_image = readfits(cal_pack.path + vig_info.file_name, /silent)
	vig_error = readfits(cal_pack.path + vig_info.error, /silent)
	vig_image = median(vig_image, 5)
	vig_file = vig_info.file_name

@@ -32,18 +33,33 @@ function metis_vignetting, data, header, cal_pack, history = history
		vig_image = rotate(vig_image, 3)
		vig_image = reverse(vig_image, 1)

		; WARN - error for the shift of the uv vignetting function must be included

		vig_error = rotate(vig_error, 1)
		vig_error = shift(vig_error, dx, dy)
		vig_error[2047 - abs(dx) : 2047, *] = 0.
		vig_error[*, 0 : abs(dy) - 1] = 0.
		vig_error = rotate(vig_error, 3)
		vig_error = reverse(vig_error, 1)

		journal, '  UV shift applied = [' + string(dx, format = '(F0.1)') + ', ' + string(dy, format = '(F0.1)') + '] pixel'
	endif

	; WARN - error for binning must be checked

	vig_nobin = vig_image
	vig_image = rebin(vig_image, header.naxis1, header.naxis2)
	vig_image = (vig_image > 0.) < 1.
	vig_error = vig_image * sqrt(rebin((vig_error/vig_nobin)^2, header.naxis1, header.naxis2)/header.nbin)

	vig_image = (vig_image > 0.) < 1.
	mask = where(vig_image eq 0.)
	vig_image[mask] = 1.

	data = data/vig_image
	data[mask] = 0.

	if isa(error) then error += (vig_error/vig_image)^2

	if ~ isa(history) then history = !null
	history = [history, 'Vignetting correction: ', '  ' + vig_file + ' shifted by [' + string(dx, format = '(F0.1)') + ', ' + string(dy, format = '(F0.1)') + '] pixel']