Skip to content
Snippets Groups Projects
Select Git revision
  • e70511b7e96a94dd21186aec98a02aa589c4f92f
  • main default
  • develop
  • feature/arch_support
  • feature/global_refactoring
  • v1.0.0
6 results

build

Blame
  • metis_dark_uvda.pro 3.21 KiB
    function metis_dark_uvda, 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
    	
    	; WARN - temporary patch to handle local l1 fits files
    	if tsensor eq 0. then tsensor = -25.
    
    	flag_use_extrapolated = 1
    	flag_use_notfullset = 1
    	flag_normalize_dark = 0
    
    	obt_available = !null
    
    	; check if one or more dark matrices with same dit, ndit1, ndit2, and tsensor exist
    	for i = 0, n_elements(dark) - 1 do begin
    		if $
    			(dark[i].dit eq dit) and $
    			(abs(dark[i].tsensor - tsensor) lt 5) and $
    			(dark[i].ndit1 eq ndit1) and (dark[i].ndit2 eq ndit2) and $
    			flag_use_notfullset and flag_use_extrapolated  then $
    			obt_available = [obt_available, dark[i].obt_beg]
    	endfor
    
    	; if not, exit with error code
    	if ~ isa(obt_available) then begin
    		journal, 'Error 12: UV applicable dark file not found.'
    		journal
    		exit, status = 12
    	endif
    
    	; if yes, chose the temporally closest
    	delta_obt = 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
    
    	; apply normalization if required
    	if flag_normalize_dark and masking.contains('disabled', /fold) then $
    		dark_image = correzione_uvda(data, cal_pack.path + dark_file) else $
    		dark_image = float(readfits(cal_pack.path + dark_file, /silent))
    
    	dark_nobin = dark_image
    
    	; 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
    
    	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(dark_error_file) then $
    		dark_error = float(readfits(cal_pack.path + dark_error_file, /silent)) else $
    		dark_error = fltarr(1024, 1024)
    
    	; rebin it and take into account the exposure time correction
    	dark_error = dark_image * sqrt(rebin((dark_error/dark_nobin)^2, header.naxis1, header.naxis2)/nbin)
    	
    	; calculate the total uncertainty
    	if isa(error) then error += (data/cal_pack.uv_channel.radiometry.iaps_gain.value * cal_pack.uv_channel.radiometry.f_noise.value + 2. * dark_error^2)/data^2
    
    	jump:
    	s = where(data le 0., count)
    	quality_matrix[s] = 0
    	
    	if ~ isa(history) then history = !null
    	history = [history, 'Dark-current correction: ', '  ' + dark_file]
    
    	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'
    	journal, '  dark file = ' + dark_file
    	journal, '  extrapolated dark = ' + dark[i].extrapol.tolower()
    
    	return, data
    end