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

Further update of error calculation procedure

parent f5f45c1c

calculate_error.pro

+29 −17
Original line number Diff line number Diff line 
function calculate_error, cal_pack, err_info, detector

    compile_opt idl2



    if detector.toupper() eq 'UV' then channel = cal_pack.uv_channel

    if detector.toupper() eq 'VL' then channel = cal_pack.vl_channel
@@ -15,30 +17,40 @@ function calculate_error, cal_pack, err_info, detector 
    data = err_info.data.image

    s = where(data eq 0., count)



    ; calculate the initial error from dark/bias subtraction (relative, squared)

    ; calculate the initial relative error after dark/bias subtraction



    if detector.toupper() eq 'UV' then $

        error = (data * radiometry.iaps_gain.value * radiometry.f_noise.value ^ 2 + $

            2. * err_info.dark.error ^ 2) / data ^ 2

        start_error = sqrt((data * radiometry.iaps_gain.value * radiometry.f_noise.value ^ 2 + $

            2. * err_info.dark.error ^ 2)) / data



    if detector.toupper() eq 'VL' then $

        error = (data * radiometry.gain.value + 2. * err_info.bias.error ^ 2 + $

            2. * err_info.dark.error ^ 2) / data ^ 2

        start_error = sqrt((data * radiometry.gain.value + 2. * err_info.bias.error ^ 2 + $

            2. * err_info.dark.error ^ 2)) / data



    ; add error for flat-field correction

    error += (err_info.flat.error / err_info.flat.image) ^ 2

    ; error for flat-field correction



    ; add error for vignetting correction

    error += (err_info.vign.error / err_info.vign.image) ^ 2

    flat_error = err_info.flat.error / err_info.flat.image



    ; add error for radiometric calibration

    error += (radiometry.rad_response.error / radiometry.rad_response.value) ^ 2 $

        + (cal_pack.instrument.pupil_area.error / cal_pack.instrument.pupil_area.value) ^ 2 $

        + (channel.angular_pixel.error / channel.angular_pixel.value) ^ 2 $

        + (unit_error / unit_factor) ^ 2

    ; error for vignetting correction



    if detector.toupper() eq 'UV' then $

        error += (err_info.eff.error / err_info.eff.image) ^ 2

    vign_error = err_info.vign.error / err_info.vign.image



    ; error for radiometric calibration



    rad_error = radiometry.rad_response.error / radiometry.rad_response.value $

        + cal_pack.instrument.pupil_area.error / cal_pack.instrument.pupil_area.value $

        + channel.angular_pixel.error / channel.angular_pixel.value $

        + unit_error / unit_factor



    ; error specific for efficiency (only for the uv channel)



    if detector.toupper() eq 'UV' then begin

        eff_error = err_info.eff.error / err_info.eff.image

    endif else eff_error = data * 0.



    ; total relative error



    error = sqrt(start_error ^ 2 + flat_error ^ 2) ;+ vign_error + eff_error + rad_error



    if count gt 0 then error[s] = 0.

calculate_error_polariz.pro

+47 −64
Original line number Diff line number Diff line 
function calculate_error_polariz, cal_pack, err_info

    compile_opt idl2



    channel = cal_pack.vl_channel

    radiometry = channel.radiometry
@@ -8,74 +10,55 @@ function calculate_error_polariz, cal_pack, err_info 
    q_error = 0.

    u_error = 0.



    ; error over stokes i

    ; initial errors over stokes parameters



    for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared)

        error[*, *, k] = data[*, *, k] * radiometry.gain.value $

    for k = 0, 3 do $ ; initial absolute error for i after dark/bias subtraction

        error[*, *, k] = sqrt(data[*, *, k] * radiometry.gain.value $

            + 2. * err_info.bias.error ^ 2 $

            + 2. * err_info.dark.error^2

            + 2. * err_info.dark.error ^ 2)



    for k = 0, 3 do $ ; error for demodulation (relative, squared)

        i_error += (data[*, *, k] ^ 2 * err_info.demod.error[*, *, k, 0] ^ 2 $

            + err_info.demod.image[*, *, k, 0]^2 * error[*, *, k])/err_info.stokes.image[*, *, 0]^2



    ; add error for flat-field correction

    i_error += (err_info.flat.error/err_info.flat.image)^2



    ; add error for vignetting correction

    i_error += (err_info.vign.error/err_info.vign.image)^2

            + err_info.demod.image[*, *, k, 0] ^ 2 * error[*, *, k] ^ 2) / err_info.stokes.image[*, *, 0] ^ 2



    ; add error for radiometric calibration

    i_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $

        + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $

        + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $

        + (radiometry.msb.error/radiometry.msb.value)^2



    ; error over stokes q

    

    for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared)

        error[*, *, k] = data[*, *, k] * radiometry.gain.value $

    for k = 0, 3 do $ ; initial absolute error for for q and u without dark subtraction

        error[*, *, k] = sqrt(data[*, *, k] * radiometry.gain.value $

            + err_info.bias.error ^ 2 $

            + err_info.dark.error^2

            + err_info.dark.error ^ 2)



    for k = 0, 3 do $ ; error for demodulation (relative, squared)

    for k = 0, 3 do begin ; error for demodulation (relative, squared)

        q_error += (data[*, *, k] ^ 2 * err_info.demod.error[*, *, k, 1] ^ 2 $

            + err_info.demod.image[*, *, k, 1]^2 * error[*, *, k])/err_info.stokes.image[*, *, 1]^2

            + err_info.demod.image[*, *, k, 1] ^ 2 * error[*, *, k] ^ 2) / err_info.stokes.image[*, *, 1] ^ 2

        u_error += (data[*, *, k] ^ 2 * err_info.demod.error[*, *, k, 2] ^ 2 $

            + err_info.demod.image[*, *, k, 2] ^ 2 * error[*, *, k] ^ 2) / err_info.stokes.image[*, *, 2] ^ 2

    endfor



    ; add error for flat-field correction

    q_error += (err_info.flat.error/err_info.flat.image)^2

    ; relative errors



    ; add error for vignetting correction

    q_error += (err_info.vign.error/err_info.vign.image)^2

    i_error = sqrt(i_error)

    q_error = sqrt(q_error)

    u_error = sqrt(u_error)



    ; add error for radiometric calibration

    q_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $

        + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $

        + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $

        + (radiometry.msb.error/radiometry.msb.value)^2

    ; error for flat-field correction



    ; error over stokes u

    flat_error = err_info.flat.error / err_info.flat.image



    for k = 0, 3 do $ ; initial error for dark/bias subtraction (absolute, squared)

        error[*, *, k] = data[*, *, k] * radiometry.gain.value $

            + err_info.bias.error^2 $

            + err_info.dark.error^2

    ; error for vignetting correction



    for k = 0, 3 do $ ; error for demodulation (relative, squared)

        u_error += (data[*, *, k]^2 * err_info.demod.error[*, *, k, 2]^2 $

            + err_info.demod.image[*, *, k, 2]^2 * error[*, *, k])/err_info.stokes.image[*, *, 2]^2

    vign_error = err_info.vign.error / err_info.vign.image



    ; error for radiometric calibration



    ; add error for flat-field correction

    u_error += (err_info.flat.error/err_info.flat.image)^2

    rad_error = radiometry.rad_response.error / radiometry.rad_response.value $

        + cal_pack.instrument.pupil_area.error / cal_pack.instrument.pupil_area.value $

        + channel.angular_pixel.error / channel.angular_pixel.value $

        + radiometry.msb.error / radiometry.msb.value



    ; add error for vignetting correction

    u_error += (err_info.vign.error/err_info.vign.image)^2

    ; total relative errors



    ; add error for radiometric calibration

    u_error += (radiometry.rad_response.error/radiometry.rad_response.value)^2 $

        + (cal_pack.instrument.pupil_area.error/cal_pack.instrument.pupil_area.value)^2 $

        + (channel.angular_pixel.error/channel.angular_pixel.value)^2 $

        + (radiometry.msb.error/radiometry.msb.value)^2

    i_tot_error = sqrt(i_error ^ 2 + flat_error ^ 2) ;+ vign_error + rad_error

    q_tot_error = sqrt(q_error ^ 2 + flat_error ^ 2) ;+ vign_error + rad_error

    u_tot_error = sqrt(u_error ^ 2 + flat_error ^ 2) ;+ vign_error + rad_error



    return, [[[i_error]], [[q_error]], [[u_error]]]

    return, [[[i_tot_error]], [[q_tot_error]], [[u_tot_error]]]

end

metis_dark_uvda.pro

+1 −1
Original line number Diff line number Diff line
@@ -145,7 +145,7 @@ function metis_dark_uvda, data, header, cal_pack, err_info, quality_matrix = qua 
    if isa(dark_syst_err, /string) then dark_syst_err = 0.



    ; 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 = sqrt(rebin(dark_error ^ 2, header.naxis1, header.naxis2) * nbin)

    dark_error = dark_error * ndit1 * ndit2



    ; save errors in the error info structure

metis_dark_uvda_new.pro

+1 −1
Original line number Diff line number Diff line
@@ -142,7 +142,7 @@ function metis_dark_uvda_new, data, header, cal_pack, err_info, quality_matrix = 
    if isa(dark_syst_err, /string) then dark_syst_err = 0.



    ; 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 = sqrt(rebin(dark_error ^ 2, header.naxis1, header.naxis2) * nbin)

    dark_error = dark_error * ndit1 * ndit2



    ; save errors in the error info structure

metis_dark_vlda.pro

+2 −2
Original line number Diff line number Diff line
@@ -84,10 +84,10 @@ function metis_dark_vlda, data, header, cal_pack, err_info, quality_matrix = qua 
    ; WARN - error for binning must be checked



    bias_image = rebin(bias_image, header.naxis1, header.naxis2) * nbin

    bias_error = sqrt(rebin(bias_error ^ 2, header.naxis1, header.naxis2) * nbin) / nbin

    bias_error = sqrt(rebin(bias_error ^ 2, header.naxis1, header.naxis2) * nbin)



    dark_image = rebin(dark_image, header.naxis1, header.naxis2) * nbin

    dark_error = sqrt(rebin(dark_error ^ 2, header.naxis1, header.naxis2) * nbin) / nbin

    dark_error = sqrt(rebin(dark_error ^ 2, header.naxis1, header.naxis2) * nbin)



    mask = where(data eq 0., count)

    data = data - (bias_image * ndit + dark_image * xposure)