Add error calculation in flat-field corr. procedure

function metis_flat_field, data, header, cal_pack, history = history

function metis_flat_field, data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history

	journal, 'Flat-field correction:'

		ff_info = cal_pack.vl_channel.flat_field

		ff_image = readfits(cal_pack.path + ff_info.file_name, /silent)

		ff_error = readfits(cal_pack.path + ff_info.file_name, exten_no = 1, /silent)

		ff_file = ff_info.file_name

		first_row = header.ref_rows eq 'Included' ? 2 : 0

		detector_size = cal_pack.vl_channel.detector_size.value

		ff_image = ff_image[*, first_row : first_row + detector_size - 1]

		ff_error = ff_error[*, first_row : first_row + detector_size - 1]

		journal, '  VL ref. rows = ' + header.ref_rows.tolower()

	endif

	if header.filter.contains('UV', /fold) then begin

		ff_info = cal_pack.uv_channel.flat_field

		nbin = sqrt(header.nbin)

		i = where(ff_info.nbin eq nbin)

		i = where(ff_info.nbin eq sqrt(header.nbin))

		ff_image = readfits(cal_pack.path + ff_info[i].file_name, /silent)

		; WARN - flat field for the uv channel must still be provided

		ff_error = ff_image * 0.

		ff_file = ff_info[i].file_name

	endif

	; WARN - error for binning must be checked

	ff_nobin = ff_image

	ff_image = rebin(ff_image, header.naxis1, header.naxis2)

	ff_error = ff_image * sqrt(rebin((ff_error/ff_nobin)^2, header.naxis1, header.naxis2)/header.nbin)

	mask = where(ff_image le 0.)

	ff_image[mask] = 1.

	data = data/ff_image

	data[mask] = 0.

	if isa(error) then error += (ff_error/ff_image)^2

	if ~ isa(history) then history = !null

	history = [history, 'Flat-field correction: ', '  ' + ff_file]