Skip to content
Snippets Groups Projects
Select Git revision
  • 7e071883b973c0116bf538c92fb5dacd71d8049a
  • master default protected
  • test
  • Version-2.1.1
  • Version-2.1.0
  • Version-2.0.8
  • Version-2.0.7
  • Version-2.0.6
  • Version-2.0.5
  • Version-2.0.2
  • Version-2.0.4
  • Version-2.0.3
  • Version-2.0.0
  • Version-1.0.2
  • Version-1.0.1
  • Version-1.0.0
  • Version-0.0.0
17 results

metis_l2_prep_uv.pro

Blame
  • metis_l2_prep_uv.pro 9.57 KiB
    pro metis_l2_prep_uv
    
    	; keyword defining if the detector reference frame must be used for the output
    
    	ref_detector = 0
    
    	; start the log
    
    	journal,'output/metis_l2_prep_log.txt'
    
    	; read the auxiliary file - here we have all the inputs we need
    
    	input = json_parse('input/contents.json', /toarray, /tostruct)
    
    	journal, 'File ' + file_basename(input.file_name)
    
    	; load the spice kernels
    
    	load_spice_kernels, input.spice_kernels, kernel_list = kernel_list, kernel_version = kernel_version
    
    	journal, 'SPICE kernel files correctly loaded:'
    	journal, '  SDK version= ' + kernel_version
    
    	; read the calibration package index
    
    	cal_pack = json_parse(input.cal_pack_path + '/index.json', /toarray, /tostruct)
    
    	; include the calibration package path into the cal_pack structure
    
    	cal_pack = create_struct('path', input.cal_pack_path + path_sep(), cal_pack)
    
    	journal, 'Calibration package correctly imported:'
    	journal, '  version = ' + cal_pack.version
    	journal, '  validity range = ' + string(cal_pack.validity_range.start, cal_pack.validity_range._end, format = '(A, "-", A)')
    
    	; read the primary hdu
    
    	data = mrdfits(input.file_name, 0, primary_header, /silent)
    
    	; read the quality matrix
    
    	quality_matrix = mrdfits(input.file_name, 'quality matrix', /silent)
    
    	; convert the string header into an idl structure
    
    	header = fits_hdr2struct(primary_header)
    
    	journal, 'L1 FITS file correctly read:'
    	journal, '  datatype = ' + string(header.datatype, format = '(I0)')
    	journal, '  sess_num = ' + header.sess_num
    	journal, '  nbin = ' + string(sqrt(header.nbin), format = '(I0)')
    
    	; error is the quadratic relative error
    	
    	error = 0.
    
    	; calibration block
    
    	case header.datatype of
    	1 : begin
    		data = double(data)
    
    		; ====================================
    		; for old l1 data
    		; data = data * header.nbin * header.ndit1 * header.ndit2
    		; header.xposure = header.dit/1000. * header.ndit1 * header.ndit2
    		; header.nsumexp = header.ndit1 * header.ndit2
    		; ====================================
    
    		data = metis_dark_uvda(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
    
    		; ====================================
    		; for data already subtracted of dark
    		; file = file_basename(header.parent, '.fits') + '_subdark.fits'
    		; data = mrdfits('UV subdark/' + file, 0, /silent)
    		; data = rebin(data, header.naxis1, header.naxis2)
    		; data = data * header.nbin * header.ndit1 * header.ndit2
    		; history = ['Dark correction: ', '  ' + file]
    		; ====================================
    
    		data = metis_flat_field(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
    		data = metis_vignetting(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
    		data = metis_rad_cal(data, header, cal_pack, error = error, quality_matrix = quality_matrix, history = history)
    
    		; ====================================
    		; for simple radiometric calibration
    		; cal_pack.uv_channel.cal_units = 'DN/s'
    		; ====================================
    
    		; ====================================
    		; temporary correction based on alfa-leo data
    		; data = data * 3.4
    		; history = [history, '  additional cal. factor from stars = 3.4']
    		; ====================================
    		
    		btype = 'UV Lyman-alpha intensity'
    		bunit = cal_pack.uv_channel.cal_units
    	end
    
    	4 : begin
    		; calibration of temporal noise images
    		btype = 'UV temporal standard deviation'
    		bunit = 'DN'
    		history = !null
    	end
    
    	6 : begin
    		; calibration of cr/sep log matrices
    		btype = 'UV cosmic-ray matrix'
    		bunit = 'DN'
    		history = !null
    	end
    
    	else : begin
    		journal, 'Error 01: wrong input data product (expected data types 1, 4, or 6).'
    		journal
    		exit, status = 1
    	end
    	endcase
    
    	; definitions for the primary header
    	; version of the fits file
    
    	version = string(input.l2_version, format = '(I02)')
    
    	; fits creation date
    
    	date = date_conv(systime(/julian, /utc), 'FITS')
    
    	; name of the fits file
    
    	file_name_fields = strsplit(header.filename, '_', /extract)
    	file_name = 'solo_L2_' + file_name_fields[2] + '_' + file_name_fields[3] + '_V' + version + '.fits'
    	out_file_name = 'output/' + file_name
    
    	; adjust the primary header
    
    	fxaddpar, primary_header, 'FILENAME', file_name
    	fxaddpar, primary_header, 'PARENT', file_basename(input.file_name)
    	fxaddpar, primary_header, 'LEVEL', 'L2'
    	fxaddpar, primary_header, 'CREATOR', 'metis_l2_prep_uv.pro'
    	fxaddpar, primary_header, 'VERS_SW', input.sw_version
    	fxaddpar, primary_header, 'VERS_CAL', cal_pack.version
    	fxaddpar, primary_header, 'VERSION', version
    	fxaddpar, primary_header, 'BTYPE', btype
    	fxaddpar, primary_header, 'BUNIT', bunit
    	fxaddpar, primary_header, 'DATAMIN', min(data, /nan)
    	fxaddpar, primary_header, 'DATAMAX', max(data, /nan)
    	fxaddpar, primary_header, 'WAVEBAND', cal_pack.uv_channel.name
    	fxaddpar, primary_header, 'XPOSURE', header.xposure
    	fxaddpar, primary_header, 'NSUMEXP', header.nsumexp
    
    	sxdelpar, primary_header, 'BLANK'
    	
    	; fix planning info keywords
    
    	if header.soopname.startswith('unknown') then soopname = 'none' else soopname = header.soopname
    	if header.obs_mode.startswith('unknown') then obs_mode = 'METIS_GENERIC_OBS' else obs_mode = header.obs_mode
    	if soopname eq 'none' then sooptype = 'none' else sooptype = header.sooptype
    	if obs_mode eq 'none' then obs_type = 'none' else obs_type = header.obs_type
    	if soopname eq 'none' and obs_mode eq 'none' then obs_id = 'none' else obs_id = header.obs_id
    
    	fxaddpar, primary_header, 'SOOPNAME', soopname
    	fxaddpar, primary_header, 'SOOPTYPE', sooptype
    	fxaddpar, primary_header, 'OBS_MODE', obs_mode
    	fxaddpar, primary_header, 'OBS_TYPE', obs_type
    	fxaddpar, primary_header, 'OBS_ID',  obs_id
    
    	; append wcs keywords
    
    	wcs = metis_wcs(header, cal_pack, ref_detector = ref_detector)
    	foreach element, wcs do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
    
    	; append solar ephemeris keywords
    
    	ephemeris = solo_get_ephemeris(header, cal_pack)
    	foreach element, ephemeris do fxaddpar, primary_header, element.name, element.value, element.comment, before = 'DATATYPE'
    
    	history = [history, 'Update WCS and solar ephemeris:', '  SKD version = ' + kernel_version]
    
    	; update the comment and history keywords
    
    	fxaddpar, primary_header, 'COMMENT', 'WARNING: UV radiometric calibration is still preliminary.'
    	fxaddpar, primary_header, 'COMMENT', 'Uncertainty matrix in the FITS extension is preliminary.'
    	
    	if ref_detector then begin
    		fxaddpar, primary_header, 'COMMENT', 'Flip vertically and rotate CROTA degrees counter-clockwise'
    		fxaddpar, primary_header, 'COMMENT', '  to have Solar North up.' 
    	endif else $
    		fxaddpar, primary_header, 'COMMENT', 'Rotate CROTA degrees counter-clockwise to have Solar North up.'
    
    	for k = 0, n_elements(history) - 1 do $
    	fxaddpar, primary_header, 'HISTORY', history[k]
    	fxaddpar, primary_header, 'HISTORY', 'L2 FITS file created on ' + date
    
    	if not ref_detector then data = metis_rectify(data, 'UV')
    	fits_add_checksum, primary_header, data
    	mwrfits, float(data), out_file_name, primary_header, /no_comment, /create, /silent
    
    	journal, 'Fits file created:'
    	journal, '  file name = ' + file_basename(out_file_name)
    
    	; add the extension with the quality matrix
    
    	base_header = primary_header
    	sxdelpar, base_header, 'SIMPLE'
    	sxdelpar, base_header, 'EXTEND'
    	sxdelpar, base_header, 'DATASUM'
    	sxdelpar, base_header, 'CHECKSUM'
    	sxdelpar, base_header, 'COMMENT'
    	sxdelpar, base_header, 'HISTORY'
    	
    	extension_header = base_header
    	fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
    	fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
    	fxaddpar, extension_header, 'EXTNAME', 'Quality matrix', 'extension name', before = 'LONGSTRN'
    	fxaddpar, extension_header, 'BTYPE', 'Pixel quality'
    	fxaddpar, extension_header, 'BUNIT', 'None'
    	fxaddpar, extension_header, 'DATAMIN', min(quality_matrix, /nan)
    	fxaddpar, extension_header, 'DATAMAX', max(quality_matrix, /nan)
    	fxaddpar, extension_header, 'COMMENT', 'Quality matrix values:'
    	fxaddpar, extension_header, 'COMMENT', '  NaN = saturated or null L0 pixel counts'
    	fxaddpar, extension_header, 'COMMENT', '  0   = unreliable pixel value'
    	fxaddpar, extension_header, 'COMMENT', '  1   = good pixel'
    	if not ref_detector then quality_matrix = metis_rectify(quality_matrix, 'UV')
    	fits_add_checksum, extension_header, quality_matrix
    	mwrfits, float(quality_matrix), out_file_name, extension_header, /no_comment, /silent
    
    	journal, 'Quality-matrix extension correctly added.'
    
    	; add the extension with the error matrix
    
    	if not ref_detector then error = metis_rectify(error, 'UV')
    	error_matrix = data * sqrt(error)
    	extension_header = base_header
    	fxaddpar, extension_header, 'PCOUNT', 0, 'parameter count', before = 'LONGSTRN'
    	fxaddpar, extension_header, 'GCOUNT', 1, 'group count', before = 'LONGSTRN'
    	fxaddpar, extension_header, 'EXTNAME', 'Error matrix', 'extension name', before = 'LONGSTRN'
    	fxaddpar, extension_header, 'BTYPE', 'Absolute error'
    	fxaddpar, extension_header, 'DATAMIN', min(error_matrix, /nan)
    	fxaddpar, extension_header, 'DATAMAX', max(error_matrix, /nan)
    	
    	fits_add_checksum, extension_header, float(error_matrix)
    	mwrfits, float(error_matrix), out_file_name, extension_header, /no_comment, /silent
    
    	journal, 'Error-matrix extension correctly added.'
    
    	; write the auxiliary information file
    
    	output = { $
    		file_name: out_file_name, $
    		l1_file_name: input.file_name, $
    		log_file_name: 'output/metis_l2_prep_log.txt' $
    	}
    
    	json_write, output, 'output/contents.json'
    
    	; unload the spice kernels
    
    	load_spice_kernels, kernel_list = kernel_list, /unload
    
    	journal, 'SPICE kernel files unloaded.'
    
    	; close the log
    
    	journal, 'Exiting without errors.'
    	journal
    	
    	exit, status = 0
    end