Version 2.2

decode_obt.pro

 function decode_obt, t_coarse, t_fine, to_double = to_double, from_double = from_double
-	if keyword_set(from_double) then begin
-		obt = t_coarse
-		t_coarse = ulong(obt)
-		t_fine = uint((obt - t_coarse) * 65536.0D)
-	endif
-	
-	if keyword_set(to_double) then begin
-		return, double(t_coarse) + double(t_fine)/65536.0D
-	endif
-
-	return, t_coarse.tostring() + ':' + t_fine.tostring()
-end
\ No newline at end of file
+	if keyword_set(to_double) then return, double(t_coarse) + double(t_fine)/65536.0D
+	if keyword_set(from_double) then obt = [ulong(t_coarse), uint((t_coarse - ulong(t_coarse)) * 65536)] else obt = [t_coarse, t_fine]
+	return, (obt.tostring()).join(':')
+end

interpol_param.pro

 function interpol_param, table, par_name, date, empty_params = empty_params
-	i = where(table.par_name eq par_name, count)
-	if count lt 1 then return, 0.0 ; this should not be necessary
-	par_date = dblarr(count)
-	par_val = fltarr(count)
-	for j = 0, count - 1 do begin
-		k = i[j]
-		if table.eng_val[k] eq 'N/A' then begin
-			table.eng_val[k] = table.raw_val[k]
-			if ~ isa(empty_params) then $
-				empty_params = table.par_name[k] else $
-					empty_params = [empty_params, table.par_name[k]]
-		endif
-		par_date[j] = date_conv(table.gen_time[k], 'JULIAN')
-		par_val[j] = float(table.eng_val[k])
-	endfor 
+	sample = where(table.par_name eq par_name, n)
+	par_date = dblarr(n)
+	for k = 0, n - 1 do par_date[k] = date_conv(table.gen_time[sample[k]], 'JULIAN')
+	if total(table.eng_val[sample].contains('N/A')) then begin
+		table.eng_val[sample] = table.raw_val[sample]
+		empty_params = [empty_params, par_name]
+	endif
+	par_val = float(table.eng_val[sample])
+	s = sort(par_date)
+	par_val = par_val[s]
+	par_date = par_date[s]
 	value = interpol(par_val, par_date, date_conv(date, 'JULIAN'))
 	if finite(value) then return, value else return, 0.0
-end
\ No newline at end of file
+end

make_bin_table.pro

+function make_bin_table, table
+    n = n_elements(table.gen_time)
+    gen_time = dblarr(n)
+    for k = 0, n - 1 do gen_time[k] = date_conv(table.gen_time[k], 'JULIAN')
+
+    w = where(table.eng_val.contains('N/A'))
+    table.eng_val[w] = table.raw_val[w]
+    par_names = table.par_name
+    s = sort(par_names)
+    par_names = par_names[s]
+    u = uniq(par_names)
+    par_names = par_names[u]
+
+    bin_table = !null
+    foreach par_name, par_names do begin
+        sample = where(table.par_name eq par_name)
+
+        min_time = min(gen_time[sample], i)
+        j = sample[i]
+        param = create_struct( $
+            'PAR_NAME', par_name, $
+            'GEN_TIME', date_conv(gen_time[j], 'FITS'), $
+            'RAW_VAL', table.raw_val[j], $
+            'ENG_VAL', table.eng_val[j], $
+            'UNIT', table.unit[j], $
+            'DESCR', table.desc[j])
+        bin_table = [bin_table, param]
+
+        max_time = max(gen_time[sample], i)
+        j = sample[i]
+        param = create_struct( $
+            'PAR_NAME', par_name, $
+            'GEN_TIME', date_conv(gen_time[j], 'FITS'), $
+            'RAW_VAL', table.raw_val[j], $
+            'ENG_VAL', table.eng_val[j], $
+            'UNIT', table.unit[j], $
+            'DESCR', table.desc[j])
+        bin_table = [bin_table, param]
+
+        avg_time = (min_time + max_time)/2
+        i = max(where(gen_time[sample] le avg_time))
+        j = sample[i]
+        param = create_struct( $
+            'PAR_NAME', par_name, $
+            'GEN_TIME', date_conv(gen_time[j], 'FITS'), $
+            'RAW_VAL', table.raw_val[j], $
+            'ENG_VAL', table.eng_val[j], $
+            'UNIT', table.unit[j], $
+            'DESCR', table.desc[j])
+        bin_table = [bin_table, param]
+
+        if product(table.eng_val[sample].matches('^-?[0-9]+')) then begin
+            min_val = min(float(table.eng_val[sample]), i)
+            j = sample[i]
+            param = create_struct( $
+                'PAR_NAME', par_name, $
+                'GEN_TIME', date_conv(gen_time[j], 'FITS'), $
+                'RAW_VAL', table.raw_val[j], $
+                'ENG_VAL', table.eng_val[j], $
+                'UNIT', table.unit[j], $
+                'DESCR', 'Min. value of ' + table.desc[j])
+            bin_table = [bin_table, param]
+
+            max_val = max(float(table.eng_val[sample]), i)
+            j = sample[i]
+            param = create_struct( $
+                'PAR_NAME', par_name, $
+                'GEN_TIME', date_conv(gen_time[j], 'FITS'), $
+                'RAW_VAL', table.raw_val[j], $
+                'ENG_VAL', table.eng_val[j], $
+                'UNIT', table.unit[j], $
+                'DESCR', 'Max. value of ' + table.desc[j])
+            bin_table = [bin_table, param]
+
+            param = create_struct( $
+                'PAR_NAME', par_name, $
+                'GEN_TIME', 'N/A', $
+                'RAW_VAL', string(mean(fix(table.raw_val[sample])), format = '(I0)'), $
+                'ENG_VAL', string(mean(double(table.eng_val[sample])), format = '(D0.16)'), $
+                'UNIT', (table.unit[sample])[0], $
+                'DESCR', 'Average value of ' + (table.desc[sample])[0])
+            bin_table = [bin_table, param]
+        endif
+    endforeach
+
+    return, bin_table
+end

metis_l1_prep.pro

 pro metis_l1_prep
 
-	; start the log 
+	; start the log
 
 	journal, 'output/metis_l1_prep_log.txt'
 
 	; some definitions
 
 	metis_datatype = list( $
-		'vl-image', $				; 0
-		'uv-image', $				; 1
-		'pcu-accumul', $			; 2
-		'vl-temp-matrix', $			; 3
-		'uv-temp-matrix', $			; 4
-		'vl-c-ray-mat', $			; 5
-		'uv-c-ray-mat', $			; 6
-		'pcu-event-list', $			; 7
-		'pcu-test-event-list', $	; 8
-		'light-curve' $				; 9
+		'VL image', $						; 0
+		'UV image', $						; 1
+		'PCU accumulation matrix', $		; 2
+		'VL temporal standard deviation', $	; 3
+		'UV temporal standard deviation', $	; 4
+		'VL cosmic-ray matrix', $			; 5
+		'UV cosmic-ray matrix', $			; 6
+		'PCU event list', $					; 7
+		'PCU test event list', $			; 8
+		'VL light-curve' $					; 9
 	)
 
 	metis_obj_size = list( $
@@ -55,7 +55,7 @@ pro metis_l1_prep
 	fits_info, input.file_name, n_ext = n_ext, extname = extname, /silent
 
 	; read the primary hdu
-	
+
 	ext_no = 0
 	data = mrdfits(input.file_name, ext_no, primary_header, /silent)
 
@@ -120,25 +120,26 @@ pro metis_l1_prep
 	date_end = solo_obt2utc(decode_obt(obt_end, /from_double))
 	date_avg = solo_obt2utc(decode_obt(obt_avg, /from_double))
 
-	date_beg_string = strreplace(strreplace(strmid(date_beg, 0, 19), '-', ''), ':', '')
+	date_beg_string = strmid(date_beg, 0, 19)
+	date_beg_string = date_beg_string.replace('-', '')
+	date_beg_string = date_beg_string.replace(':', '')
 
 	journal, 'UTC time of start acquisition = ' + date_beg
 
 	; name of the fits file
-	
+
 	file_name_fields = strsplit(fxpar(primary_header, 'FILENAME'), '_', /extract)
-	
+
 	file_name = 'solo_l1_' + file_name_fields[2] + '_' + date_beg_string + '_v' + version + '.fits'
 	out_file_name = 'output/' + file_name
 
 	; exposure times
 
-	dit = float(fxpar(metadata_extension_header, 'DIT'))
-	telapse = float(obt_end - obt_beg)
+	dit = fxpar(metadata_extension_header, 'DIT')
+	telapse = obt_end - obt_beg
 	xposure = dit/1000.
 
 	; instrument keywords
-	; TODO - complete with filter information
 
 	if datatype eq 0 or datatype eq 3 or datatype eq 5 or datatype eq 9 then begin
 		filter = 'VL'
@@ -175,15 +176,15 @@ pro metis_l1_prep
 		metadata_extension_header[j + 1 : k - 1], $
 		primary_header[i : *] $
 	]
-	
+
 	; rebin the image if binning was applied during the acquisition
 
 	if fxpar(primary_header, 'COMPR') then bin_type = fxpar(primary_header, 'B0_BIN') else bin_type = 0
-	
+
 	pol_id = fxpar(primary_header, 'POL_ID')
 	fpfilt = fxpar(primary_header, 'VLFPFILT', count = count)
 	if count gt 0 and pol_id eq 0 and fpfilt eq 0 and bin_type eq 0 then bin_type = 1
-	
+
 	bin_fact = 2. ^ bin_type
 	if bin_fact gt 1. then begin
 		data = rebin(data, naxis1 / bin_fact, naxis2 / bin_fact, /sample)
@@ -208,31 +209,31 @@ pro metis_l1_prep
 	fxaddpar, primary_header, 'DATE-AVG', date_avg, 'Average time of observation', before = 'OBT_BEG'
 	fxaddpar, primary_header, 'DATE-END', date_end, 'End time of observation', before = 'OBT_BEG'
 	fxaddpar, primary_header, 'TIMESYS', 'UTC', 'System used for time keywords', before = 'OBT_BEG'
-	fxaddpar, primary_header, 'TIMRDER', 0.0, 'Estimated random error in time values', before = 'OBT_BEG'
-	fxaddpar, primary_header, 'TIMSYER', 0.0, 'Estimated systematic error in time values', before = 'OBT_BEG'
+	fxaddpar, primary_header, 'TIMRDER', 0.0, '[s] Estimated random error in time values', before = 'OBT_BEG'
+	fxaddpar, primary_header, 'TIMSYER', 0.0, '[s] Estimated systematic error in time values', before = 'OBT_BEG'
 	fxaddpar, primary_header, 'LEVEL', 'L1'
 	fxaddpar, primary_header, 'CREATOR', 'metis_l1_prep.pro'
 	fxaddpar, primary_header, 'OBSRVTRY', 'Solar Orbiter', 'Satellite name', before = 'INSTRUME'
 	fxaddpar, primary_header, 'TELESCOP', telescope, 'Telescope that took the measurement', before = 'INSTRUME'
 	fxaddpar, primary_header, 'DETECTOR', detector, 'Subunit/sensor', before = 'DATAMIN'
 	fxaddpar, primary_header, 'OBJECT', 'TBD', 'The use of the keyword OBJECT is [TBD]', before = 'DATAMIN'
-	fxaddpar, primary_header, 'OBS_MODE', planning_data.obs_mode, 'Observation mode or study that has been used to acquire this image', before = 'DATAMIN'
-	fxaddpar, primary_header, 'OBS_TYPE', obs_type, 'Subfield of OBS_ID that only contains an encoded version of OBS_MODE', before = 'DATAMIN'
+	fxaddpar, primary_header, 'OBS_MODE', planning_data.obs_mode, 'Observation mode', before = 'DATAMIN'
+	fxaddpar, primary_header, 'OBS_TYPE', obs_type, 'Encoded version of OBS_MODE', before = 'DATAMIN'
 	fxaddpar, primary_header, 'FILTER', filter, 'Filter used to acquire this image', before = 'DATAMIN'
-	fxaddpar, primary_header, 'WAVELNTH', wavelnth, 'Characteristic wavelength at which the observation was taken', before = 'DATAMIN'
-	fxaddpar, primary_header, 'WAVEMIN', wavemin, 'The shortest wavelength at which the net (approximate) response function becomes 0.05 times the maximum response. ', before = 'DATAMIN'
-	fxaddpar, primary_header, 'WAVEMAX', wavemax, 'The longest wavelength at which the net (approximate) response function becomes 0.05 times the maximum response', before = 'DATAMIN'
-	fxaddpar, primary_header, 'WAVEBAND', waveband, 'Description of the wavelength band', before = 'DATAMIN'
-	fxaddpar, primary_header, 'XPOSURE', xposure, 'Total effective exposure time of the observation, in seconds', before = 'DATAMIN'
-	fxaddpar, primary_header, 'NSUMEXP', 1, 'Number of images summed together to form the observation', before = 'DATAMIN'
-	fxaddpar, primary_header, 'TELAPSE', telapse, 'Total elapsed time between the beginning and end of the complete observation in seconds, including any dead times between exposures', before = 'DATAMIN'
-	fxaddpar, primary_header, 'SOOPNAME', planning_data.soop_name, 'SOOP(s) that this observation belongs to', before = 'DATAMIN'
-	fxaddpar, primary_header, 'SOOPTYPE', soop_type, before = 'DATAMIN'
-	fxaddpar, primary_header, 'OBS_ID', planning_data.obs_id, 'Unique identifier for the observation that is associated with the data acquisition', before = 'DATAMIN' 
-	fxaddpar, primary_header, 'TARGET', 'TBD', 'Taget as defined in the SOOP description', before = 'DATAMIN'
-	fxaddpar, primary_header, 'BSCALE', 1.0, before = 'DATAMIN'
-	fxaddpar, primary_header, 'BZERO', 0.0, before = 'DATAMIN'
-	fxaddpar, primary_header, 'BTYPE', strupcase(metis_datatype[datatype]), 'Science data object type', before = 'DATAMIN'
+	fxaddpar, primary_header, 'WAVELNTH', wavelnth, '[nm] Characteristic wavelength of observation', before = 'DATAMIN'
+	fxaddpar, primary_header, 'WAVEMIN', wavemin, '[nm] Min. wavelength where response > 0.05 of max.', before = 'DATAMIN'
+	fxaddpar, primary_header, 'WAVEMAX', wavemax, '[nm] Max. wavelength where response > 0.05 of max.', before = 'DATAMIN'
+	fxaddpar, primary_header, 'WAVEBAND', waveband, 'Bandpass description', before = 'DATAMIN'
+	fxaddpar, primary_header, 'XPOSURE', xposure, '[s] Total effective exposure time', before = 'DATAMIN'
+	fxaddpar, primary_header, 'NSUMEXP', 1, 'Number of detector readouts summed together', before = 'DATAMIN'
+	fxaddpar, primary_header, 'TELAPSE', telapse, '[s] Elapsed time between beginning and end of observation', before = 'DATAMIN'
+	fxaddpar, primary_header, 'SOOPNAME', planning_data.soop_name, 'Name of the SOOP campaign that the data belong to', before = 'DATAMIN'
+	fxaddpar, primary_header, 'SOOPTYPE', soop_type, 'Campaign ID(s) that the data belong to', before = 'DATAMIN'
+	fxaddpar, primary_header, 'OBS_ID', planning_data.obs_id, 'Unique ID of the individual observation', before = 'DATAMIN'
+	fxaddpar, primary_header, 'TARGET', 'TBD', 'Type of target from planning', before = 'DATAMIN'
+	fxaddpar, primary_header, 'BSCALE', 1., 'Ratio of physical to array value at 0 offset', before = 'DATAMIN'
+	fxaddpar, primary_header, 'BZERO', 0., 'Physical value for the array value 0', before = 'DATAMIN'
+	fxaddpar, primary_header, 'BTYPE', metis_datatype[datatype], 'Science data object type', before = 'DATAMIN'
 	fxaddpar, primary_header, 'BUNIT', 'DN', 'Units of physical value, after application of BSCALE and BZERO', before = 'DATAMIN'
 	if datatype le 6 then begin
 		fxaddpar, primary_header, 'PXBEG1', 1, 'First pixel that has been read out in dimension 1', before = 'COMPRESS'
@@ -252,19 +253,21 @@ pro metis_l1_prep
 
 	; replace raw values with calibrated values in the primary header
 
+	empty_params = !null
+
 	if datatype eq 0 or datatype eq 3 or datatype eq 5 then begin
-		
-		; WARN - DACPOL parameters are not calibrated since a calibration curve does not exist in the IDB. Their calibration in physical units (e.g., voltages or angles) should be done later
-
-		; fxaddpar, primary_header, 'DAC1POL1', interpol_param(hk_table, 'NIT0E061', date_avg, empty_param = empty_param)
-		; fxaddpar, primary_header, 'DAC2POL1', interpol_param(hk_table, 'NIT0E062', date_avg, empty_param = empty_param)
-		; fxaddpar, primary_header, 'DAC1POL2', interpol_param(hk_table, 'NIT0E064', date_avg, empty_param = empty_param)
-		; fxaddpar, primary_header, 'DAC2POL2', interpol_param(hk_table, 'NIT0E065', date_avg, empty_param = empty_param)
-		; fxaddpar, primary_header, 'DAC1POL3', interpol_param(hk_table, 'NIT0E067', date_avg, empty_param = empty_param)
-		; fxaddpar, primary_header, 'DAC2POL3', interpol_param(hk_table, 'NIT0E068', date_avg, empty_param = empty_param)
-		; fxaddpar, primary_header, 'DAC1POL4', interpol_param(hk_table, 'NIT0E06A', date_avg, empty_param = empty_param)
-		; fxaddpar, primary_header, 'DAC2POL4', interpol_param(hk_table, 'NIT0E06B', date_avg, empty_param = empty_param)
-		
+
+		; NOTE - DACPOL parameters are not calibrated since a calibration curve does not exist in the IDB. Their calibration in physical units (e.g., voltages or angles) should be done later
+
+		; fxaddpar, primary_header, 'DAC1POL1', interpol_param(hk_table, 'NIT0E061', date_avg, empty_params = empty_params)
+		; fxaddpar, primary_header, 'DAC2POL1', interpol_param(hk_table, 'NIT0E062', date_avg, empty_params = empty_params)
+		; fxaddpar, primary_header, 'DAC1POL2', interpol_param(hk_table, 'NIT0E064', date_avg, empty_params = empty_params)
+		; fxaddpar, primary_header, 'DAC2POL2', interpol_param(hk_table, 'NIT0E065', date_avg, empty_params = empty_params)
+		; fxaddpar, primary_header, 'DAC1POL3', interpol_param(hk_table, 'NIT0E067', date_avg, empty_params = empty_params)
+		; fxaddpar, primary_header, 'DAC2POL3', interpol_param(hk_table, 'NIT0E068', date_avg, empty_params = empty_params)
+		; fxaddpar, primary_header, 'DAC1POL4', interpol_param(hk_table, 'NIT0E06A', date_avg, empty_params = empty_params)
+		; fxaddpar, primary_header, 'DAC2POL4', interpol_param(hk_table, 'NIT0E06B', date_avg, empty_params = empty_params)
+
 		fxaddpar, primary_header, 'TSENSOR ', interpol_param(hk_table, 'NIT0E0E0', date_avg, empty_params = empty_params)
 		fxaddpar, primary_header, 'PMPTEMP ', interpol_param(hk_table, 'NIT0L00D', date_avg, empty_params = empty_params)
 
@@ -276,18 +279,24 @@ pro metis_l1_prep
 	if datatype eq 1 or datatype eq 4 or datatype eq 6 then begin
 		fxaddpar, primary_header, 'HVU_SCR ', interpol_param(hk_table, 'NIT0E070', date_avg, empty_params = empty_params)
 		fxaddpar, primary_header, 'HVU_MCP ', interpol_param(hk_table, 'NIT0E071', date_avg, empty_params = empty_params)
+		fxaddpar, primary_header, 'HV_SCR_V', interpol_param(hk_table, 'NIT0E0B7', date_avg, empty_params = empty_params), after = 'HVU_MCP'
+		fxaddpar, primary_header, 'HV_MCP_V', interpol_param(hk_table, 'NIT0E0B6', date_avg, empty_params = empty_params), after = 'HV_SCR_V'
+		fxaddpar, primary_header, 'HV_MCP_I', interpol_param(hk_table, 'NIT0E0BF', date_avg, empty_params = empty_params), after = 'HV_MCP_V'
 		fxaddpar, primary_header, 'TSENSOR ', interpol_param(hk_table, 'NIT0E050', date_avg, empty_params = empty_params)
 
 		journal, 'Header keywords were calibrated using HK parameters:'
-		journal, '  HVU_SCR = ' + string(fxpar(primary_header, 'HVU_SCR'), format = '(F0)')
-		journal, '  HVU_MCP = ' + string(fxpar(primary_header, 'HVU_MCP'), format = '(F0)')
-		journal, '  TSENSOR = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)')
+		journal, '  HVU_SCR  = ' + string(fxpar(primary_header, 'HVU_SCR'), format = '(F0)')
+		journal, '  HVU_MCP  = ' + string(fxpar(primary_header, 'HVU_MCP'), format = '(F0)')
+		journal, '  HV_SCR_V = ' + string(fxpar(primary_header, 'HV_SCR_V'), format = '(F0)')
+		journal, '  HV_MCP_V = ' + string(fxpar(primary_header, 'HV_MCP_V'), format = '(F0)')
+		journal, '  HV_MCP_I = ' + string(fxpar(primary_header, 'HV_MCP_I'), format = '(F0)')
+		journal, '  TSENSOR  = ' + string(fxpar(primary_header, 'TSENSOR'), format = '(F0)')
 	endif
 
-	if ~ isa(empty_params) then empty_params = ''
+	if empty_params eq !null then empty_params = ''
 
 	; replace verbose keyword values
-	
+
 	key = fxpar(primary_header, 'MEASKIND', count = count)
 	if count gt 0 then begin
 		pol_id = fxpar(primary_header, 'POL_ID')
@@ -348,12 +357,7 @@ pro metis_l1_prep
 	; modify keywords for file history
 
 	date = date_conv(systime(/julian, /utc), 'FITS')
-	
-	old_history = fxpar(primary_header, 'HISTORY')
-	sxdelpar, primary_header, 'HISTORY'
-	history = ['L1 FITS file created on ' + date, old_history]
-	for k = 0, n_elements(history) - 1 do $
-		fxaddpar, primary_header, 'HISTORY', history[k]
+	fxaddpar, primary_header, 'HISTORY', 'L1 FITS file created on ' + date
 
 	; modify keywords for comments
 
@@ -374,7 +378,6 @@ pro metis_l1_prep
 	endif
 
 	; add checksum and datasum to the fits header
-	; WARN - should this be done at the end? i don't know
 
 	fits_add_checksum, primary_header, image
 
@@ -383,7 +386,6 @@ pro metis_l1_prep
 	journal, 'Fits file created:'
 	journal, '  file name = ' + out_file_name
 
-
 	; if applicable, save the data binary-table extension as it is
 
 	if isa(data_bin_table) then mwrfits, data_bin_table, 'output/' + file_name, data_extension_header, /no_comment, /silent
@@ -400,19 +402,7 @@ pro metis_l1_prep
 	fxaddpar, hk_extension_header, 'GCOUNT', 1, 'Group count'
 	fxaddpar, hk_extension_header, 'EXTNAME', 'House-keeping', 'Extension name'
 
-	for i = 0, n_elements(hk_table.par_name) - 1 do begin
-		param = create_struct( $
-			'PAR_NAME', hk_table.par_name[i], $
-			'PACKET', hk_table.packet[i], $
-			'GEN_TIME', hk_table.gen_time[i], $
-			'REC_TIME', hk_table.rec_time[i], $
-			'RAW_VAL', hk_table.raw_val[i], $
-			'ENG_VAL', hk_table.eng_val[i], $
-			'UNIT', hk_table.unit[i], $
-			'DESCR', hk_table.desc[i] $
-		)
-		if ~ isa(hk_bin_table) then hk_bin_table = param else hk_bin_table = [hk_bin_table, param]
-	endfor
+	hk_bin_table = make_bin_table(hk_table)
 
 	mwrfits, hk_bin_table, out_file_name, hk_extension_header, /no_comment, /silent
 
@@ -432,7 +422,7 @@ pro metis_l1_prep
 		log_file_name: 'output/metis_l1_prep_log.txt', $
 		empty_params : empty_params[uniq(empty_params)]  $
  	}
-	
+
 	json_write, output, 'output/contents.json'
 
 	journal, 'Output saved.'
@@ -447,4 +437,4 @@ pro metis_l1_prep
 	journal, 'Errors occurred while processing.'
 	journal
 	exit, status = 1
-end
\ No newline at end of file
+end

strreplace.pro

-function strreplace, input, substring, replacement
-	return, strjoin(strsplit(input, substring, /extract, /preserve), replacement)
-end