Commit 6291b19e authored by Roberto Susino's avatar Roberto Susino
Browse files

Merge branch 'test' 

# Conflicts:
#	metis_l2_prep_vl_polariz.pro
parents 9e472063 3156e371

calculate_error.pro

0 → 100644
+58 −0
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



    radiometry = channel.radiometry



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

        unit_factor = radiometry.msb.value

        unit_error = radiometry.msb.error

    endif else begin

        unit_factor = 1.

        unit_error = 0.

    endelse



    data = err_info.data.image

    s = where(data eq 0., count)



    ; calculate the initial relative error after dark/bias subtraction



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

        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 $

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

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



    ; error for flat-field correction



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



    ; error for vignetting correction



    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.



    return, error

end

calculate_error_polariz.pro

0 → 100644
+64 −0
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



    data = err_info.data.image

    error = data * 0.

    i_error = 0.

    q_error = 0.

    u_error = 0.



    ; initial errors over stokes parameters



    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)



    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] ^ 2) / err_info.stokes.image[*, *, 0] ^ 2



    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)



    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] ^ 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



    ; relative errors



    i_error = sqrt(i_error)

    q_error = sqrt(q_error)

    u_error = sqrt(u_error)



    ; error for flat-field correction



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



    ; error for vignetting correction



    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 $

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



    ; total relative errors



    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_tot_error]], [[q_tot_error]], [[u_tot_error]]]

end

check_compression.pro

0 → 100644
+44 −0
Original line number Diff line number Diff line 
pro check_compression, image, header

    ; Check compression



    compr = fxpar(header, 'COMPR', missing = '')

    delta_q = fxpar(header, 'B0_DQ', missing = -1)



    if compr.contains('enabled', /fold) then fxaddpar, header, 'COMPRESS', 'Lossless'

    if delta_q ge 0 then if delta_q gt 0 then fxaddpar, header, 'COMPRESS', 'Lossy'



    ; Check binning



    bin_type = fxpar(header, 'BIN_TYPE')

    bin_size = 2 ^ bin_type



    if fxpar(header, 'NBIN') ne bin_size ^ 2 then begin

        naxis1 = fxpar(header, 'NAXIS1') / bin_size

        naxis2 = fxpar(header, 'NAXIS2') / bin_size

        image = rebin(image, naxis1, naxis2)

        image = long(image) * bin_size ^ 2



        fxaddpar, header, 'NAXIS1', naxis1

        fxaddpar, header, 'NAXIS2', naxis2

        fxaddpar, header, 'NBIN1', bin_size

        fxaddpar, header, 'NBIN2', bin_size

        fxaddpar, header, 'NBIN', bin_size ^ 2

        fxaddpar, header, 'PXEND1', naxis1

        fxaddpar, header, 'PXEND2', naxis2

        fxaddpar, header, 'FIRSTROW', 0, before = 'CHECKSUM'

        fxaddpar, header, 'B0_BIN', bin_type, before = 'CHECKSUM'

        fxaddpar, header, 'B0_DQ', 0, before = 'CHECKSUM'

        fxaddpar, header, 'B0_STOP', 2048, before = 'CHECKSUM'

        fxaddpar, header, 'B1_BIN', 0, before = 'CHECKSUM'

        fxaddpar, header, 'B1_DQ', 0, before = 'CHECKSUM'

        fxaddpar, header, 'B1_STOP', 0, before = 'CHECKSUM'

        fxaddpar, header, 'B2_BIN', 0, before = 'CHECKSUM'

        fxaddpar, header, 'B2_DQ', 0, before = 'CHECKSUM'

        fxaddpar, header, 'B2_STOP', 0, before = 'CHECKSUM'



        comment = fxpar(header, 'COMMENT')

        comment = ['Image was rebinned according to the commanded binning factor.', comment]

        sxdelpar, header, 'COMMENT'

        foreach line, comment do fxaddpar, header, 'COMMENT', line

    endif

end

compute_target.pro

0 → 100644
+19 −0
Original line number Diff line number Diff line 
function compute_target, header

    p1 = fxpar(header, 'SC_YAW')

    p2 = fxpar(header, 'SC_PITCH')

    r = fxpar(header, 'RSUN_ARC')



    rho = sqrt(p1 ^ 2 + p2 ^ 2) / r

    psi = (atan(p2, p1) * !radeg + 270) mod 360



    if rho le 0.1 then target = 'On disc, disc centre'

    if rho gt 0.1 and rho lt 0.9 then target = 'On disc'

    if rho ge 0.9 then begin

        target = 'Limb'

        if psi le 1 or psi ge 359 then target += ', North Pole'

        if psi gt 1 and psi lt 179 then target += ', East'

        if psi ge 179 and psi le 181 then target += ', South Pole'

        if psi gt 181 and psi lt 359 then target += ', West'

    endif

    return, target

end

fits_hdr2struct.pro

+12 −11
Original line number Diff line number Diff line 
function fits_hdr2struct, string_hdr

    compile_opt idl2



    struct = !null

    n = n_elements(string_hdr)