API Reference
Functions for reading in data files
- DiadFit.importing_data_files.add_column_name_descriptions(df)[source]
Adds a new to inputted dataframe, with a description of what the diadfit columns mean underneath
- Parameters:
df (pandas dataframe, including some columns from DiadFit (can have other columns too))
- Returns:
df
- Return type:
A dataframe with a new row, with descriptions for columns matching our reference key.
- DiadFit.importing_data_files.calculates_time_witec(*, path, filename)[source]
calculates time as seconds after midnight for non video files for WITEC files
- DiadFit.importing_data_files.check_for_duplicates(spectra_path, prefix=True, prefix_str=' ', exception=True)[source]
This function checks for duplicate filenames in a specified directory and prints the duplicates if found.
Parameters: spectra_path (str):
The path of the directory containing the files to be checked for duplicates.
- prefix (bool):
If True, the function will remove the specified prefix string from the filenames before checking for duplicates. Default is True.
- prefix_str (str:
The prefix string to be removed from filenames if ‘prefix’ is set to True. Default is a single space ‘ ‘.
Returns: file_m (numpy.ndarray): A numpy array containing the modified filenames after removing the prefix (if specified).
- DiadFit.importing_data_files.checks_if_general_witec(*, path, filename)[source]
Checks if a WITEC file is a spectra file with all the right metadata
- DiadFit.importing_data_files.checks_if_imagescan_witec(*, path, filename)[source]
Checks if a WITEC file is an imagescan (as doesnt have all metadata)
- DiadFit.importing_data_files.checks_if_video_witec(*, path, filename)[source]
Checks if a WITEC file is an image (as doesnt have all metadata)
- DiadFit.importing_data_files.convert_datastamp_to_metadata(path, filename, creation=True, modification=False)[source]
Gets file modification or creation time, outputs as metadata in the same format as for WITEC
- Parameters:
- Return type:
df of timestamp, and other columns to have the same format as the WITEC metadata output
- DiadFit.importing_data_files.extract_accumulations_witec(*, path, filename)[source]
Extracts accumulations from a WITEC file
- DiadFit.importing_data_files.extract_acq_params_horiba(path, filename)[source]
Extracts all relevant acquisition parameters from a HORIBA file, returns as a dataframe.
- DiadFit.importing_data_files.extract_acq_params_witec(*, path, filename, trupower=False)[source]
This function checks what type of file you have, and if its a spectra file, uses the functions above to extract various bits of metadata.
- DiadFit.importing_data_files.extract_date_witec(*, path, filename)[source]
Extracts date from a WITEC file
- DiadFit.importing_data_files.extract_duration_witec(*, path, filename)[source]
Extracts analysis duration from a WITEC file
- DiadFit.importing_data_files.extract_integration_time_witec(*, path, filename)[source]
Extracts Integration time from a WITEC file
- DiadFit.importing_data_files.extract_laser_power_witec(*, path, filename)[source]
Extracts laser power from a WITEC file
- DiadFit.importing_data_files.extract_objective_witec(*, path, filename)[source]
Extracts objective magnification from a WITEC file
- DiadFit.importing_data_files.extract_spectral_center_witec(*, path, filename)[source]
Extracts Spectral Center from a WITEC file
- DiadFit.importing_data_files.extract_temp_Aranet(df)[source]
Extracts temperature data from the aranet
- DiadFit.importing_data_files.extract_time_stamp_witec(*, path, filename)[source]
Extracts time stamp from a WITEC file
- DiadFit.importing_data_files.extracting_filenames_generic(*, names, prefix=False, str_prefix=None, suffix=False, CRR_filter=True, str_suffix=None, file_ext=None)[source]
Takes filenames from a panda series (e.g., a column of a dataframe of metadata), outputs a numpy array that is consistent with the same function for spectra, to allow stitching of spectra and metadata.
- Parameters:
names (Pandas.Series) – Series of sample names, e.g., from ‘filename’ column of metadata output
prefix (bool) – if True, has a number before the file name
str_prefix (str) – The string separating the prefix from the file name (e.g. if file is 01 test, str_prefix=” “)
suffix (bool) – if True, has a number or name after the filename
str_suffix (str) – The string separating the filename from the suffix
file_ext (str) – The file extension, e.g., ‘.csv’
- Return type:
np.array of names, with prefix, suffix and filetype stripped away
- DiadFit.importing_data_files.get_all_txt_files(path)[source]
This function takes a user path, and gets all the .txt. files in that path.
- Parameters:
path (str) – Folder user wishes to read data from
- DiadFit.importing_data_files.get_data(*, path=None, filename=None, Diad_files=None, filetype='Witec_ASCII')[source]
Extracts data as a np.array from user file of differen types
- DiadFit.importing_data_files.get_files(path, ID_str=None, file_ext='txt', exclude_str=None, exclude_type=None, sort=True)[source]
This function takes a user path, and extracts all files which contain the ID_str
- Parameters:
path (str) – Folder user wishes to read data from
sort (bool) – If true, sorts files alphabetically
ID_str (list) – Finds all files containing this string (e.g. [‘Ne’, ‘NE’]
exclude_str (str) – Excludes files with this string in the name
file_ext (str) – Gets all files of this format only (e.g. txt)
- Returns:
list
- Return type:
file names as a list.
- DiadFit.importing_data_files.get_settings()[source]
This function reads the settings file saved in step 1, and loads the options
- DiadFit.importing_data_files.get_video_mag(metadata_path)[source]
This function finds all the video files in a single folder, and returns a dataframe of the filename and the magnification used.
- DiadFit.importing_data_files.line_contains_video_image(line)[source]
This function returns video image information
- DiadFit.importing_data_files.loop_convert_datastamp_to_metadata(path, files, creation=True, modification=False)[source]
Loops over multiple files to get timestamp the file was created or modified using the convert_datastamp_to_metadata function.
- path: str
Path where spectra files are stored
- files: list
list of filenames
- creation: bool
If True, gets timestamp based on creation date of file
- modification: bool
If True, gets timestamp based on modification date of file
- Return type:
df of timestamp, and other columns to have the same format as the WITEC metadata output
- DiadFit.importing_data_files.read_HORIBA_to_df(*, path=None, filename)[source]
This function takes in a HORIBA .txt. file with headers with #, and looks down to the row where Data starts (no #), and saves this to a new csv file called pandas_…. old file. It exports the data as a pandas dataframe
- DiadFit.importing_data_files.read_witec_to_df(*, path=None, filename)[source]
This function takes in a WITec ASCII.txt. file with metadata mixed with data, and looks down to the row where Data starts, and saves this to a new file called pandas_…. old file. It exports the data as a pandas dataframe
- Parameters:
path (str) – Folder user wishes to read data from
filename – Specific file being read
- DiadFit.importing_data_files.save_settings(meta_path, spectra_path, spectra_filetype, prefix, prefix_str, spectra_file_ext, meta_file_ext, TruPower)[source]
This function saves settings so you can load them across multiple notebooks without repition
- Parameters:
meta_path (str) – Path where your metadata is stored
spectra_path (str) – path where your spectra is stored
spectra_filetype (str) – Style of data. Choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
spectra_file_ext (str) – Extension of spectra file and metadatafile. e.g. ‘.txt’, ‘.csv’
meta_file_ext (str) – Extension of spectra file and metadatafile. e.g. ‘.txt’, ‘.csv’
prefix (bool) – If True, removes 01, 02, from filename (WITEC problem) Also need to state prefix_str: prefix separating string (in this case, 01 Ne would be ‘ ‘
- TruPower: bool
If WITEC instrument and you have TruPower, set as True
- Return type:
file called settings.txt with these saved.
- DiadFit.importing_data_files.stitch_metadata_in_loop_horiba(Allfiles, path=None)[source]
Stitches acquisition parameters together from the function extract_acq_params_horiba for multiple files :param AllFiles: List of all file names :type AllFiles: list :param path: Path where files are found :type path: str
- Return type:
df of aquisitoin parameters
- DiadFit.importing_data_files.stitch_metadata_in_loop_witec(*, Allfiles, path, prefix=True, trupower=False, str_prefix=' ')[source]
Stitches together WITEC metadata for all files in a loop using the function extract_acq_params_witec and calculates_time_witec, exports as a dataframe
- Parameters:
- Return type:
DataFrame of metadata parameters with a row for each file.
Functions associated with fitting diads
- DiadFit.diads.add_peak(*, prefix=None, center=None, model_name='VoigtModel', min_cent=None, max_cent=None, min_sigma=None, max_sigma=None, amplitude=100, min_amplitude=None, max_amplitude=None, sigma=0.2)[source]
This function iteratively adds a peak using lmfit, used for iteratively fitting HB, C13 etc.
- Parameters:
- model_name: str ‘VoigtModel’ or ‘PseudoVoigtModel’
Type of peak which is being added.
Optional:
- min_cent, max_cent: float
Minimum and maximum x coordinate to place peak center at
- min_sigma, max_sigma: float
Minimum and maximum sigma of peak fit
- min_amplitude, max_amplitude: float
Minimum and maximum sigma of peak fit
- Returns:
peak (Fitted model)
params (Peak fit parameters)
- DiadFit.diads.assess_diad1_skewness(*, config1: diad1_fit_config = diad1_fit_config(model_name='PseudoVoigtModel', fit_peaks=2, N_poly_bck_diad1=1, lower_bck_diad1=(1180, 1220), upper_bck_diad1=(1300, 1350), fit_gauss=False, gauss_amp=1000, diad_sigma=0.2, diad_sigma_min_allowance=0.2, diad_sigma_max_allowance=5, diad_prom=100, HB_prom=20, x_range_baseline=75, y_range_baseline=100, dpi=200, x_range_residual=20, return_other_params=False), int_cut_off=0.3, path=None, filename=None, filetype=None, skewness='abs', height=1, prominence=5, width=0.5, plot_figure=True, peak_fit_routine=False, peak_pos=None, peak_height=None, dpi=200, diad_array=None)[source]
Assesses Skewness of Diad peaks. Useful for identifying mixed L + V phases (see DeVitre et al. in review)
- Parameters:
config1 (diad1_fit_config() dataclass. Import parameters for this function are.) –
- lower_bck_diad1, upper_bck_diad1: [float, float]
background position to left and right of overal diad region (inc diad+-HB+-C13)
- N_poly_bck_diad1: int
degree of polynomial to fit to the background.
- exclude_range1, exclude_range2: Tuple[float, float]
Range to exclude
int_cut_off (float) – Value of intensity at which to calculate the skewness (e.g. 0.15X the peak height).
path (str) – Folder user wishes to read data from
filename (str) – Specific file being read
- filetype: str
Identifies type of file choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
- skewness: ‘abs’ or ‘dir’
If ‘abs’, gives absolute skewness (e.g., largest possible number regardless of direction) if ‘dir’ does skewness as one side over the other
- height, prominence, width: float
Scipy find peak parameters
- plot_figure: bool
if True, plots and saves a figure.
- dpi: int
dpi to save image at
- peak_fit_routine: bool (default False) -
if True, doesnt find max of spline, but uses the peak position you have found from the diad fitting functions. Also needs peak_pos (peak center in cm-1) and peak_height
- Return type:
pd.DataFrame with information about file and skewness parameters.
- DiadFit.diads.assess_diad2_skewness(*, config1: diad2_fit_config = diad1_fit_config(model_name='PseudoVoigtModel', fit_peaks=2, N_poly_bck_diad1=1, lower_bck_diad1=(1180, 1220), upper_bck_diad1=(1300, 1350), fit_gauss=False, gauss_amp=1000, diad_sigma=0.2, diad_sigma_min_allowance=0.2, diad_sigma_max_allowance=5, diad_prom=100, HB_prom=20, x_range_baseline=75, y_range_baseline=100, dpi=200, x_range_residual=20, return_other_params=False), int_cut_off=0.3, path=None, filename=None, filetype=None, skewness='abs', height=1, prominence=5, width=0.5, plot_figure=True, dpi=200, peak_fit_routine=False, peak_pos=None, peak_height=None, diad_array=None)[source]
Assesses Skewness of Diad peaks. Useful for identifying mixed L + V phases (see DeVitre et al. in review)
- Parameters:
config1 (diad2_fit_config() dataclass. Import parameters for this function are.) –
- lower_bck_diad2, upper_bck_diad2: [float, float]
background position to left and right of overal diad region (inc diad+-HB+-C13)
- N_poly_bck_diad1: int
degree of polynomial to fit to the background.
- exclude_range1, exclude_range2: Tuple[float, float]
Range to exclude
int_cut_off (float) – Value of intensity at which to calculate the skewness (e.g. 0.15X the peak height).
path (str) – Folder user wishes to read data from
filename (str) – Specific file being read
- filetype: str
Identifies type of file choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
- skewness: ‘abs’ or ‘dir’
If ‘abs’, gives absolute skewness (e.g., largest possible number regardless of direction) if ‘dir’ does skewness as one side over the other
- height, prominence, width: float
Scipy find peak parameters
- plot_figure: bool
if True, plots and saves a figure.
- dpi: int
dpi to save image at
- peak_fit_routine: bool (default False) -
if True, doesnt find max of spline, but uses the peak position you have found from the diad fitting functions. Also needs peak_pos (peak center in cm-1) and peak_height
- Return type:
pd.DataFrame with information about file and skewness parameters.
- DiadFit.diads.calculate_HB_Diad_area_ratio(df)[source]
“ Calculates two area ratios of hotbands and diads
- DiadFit.diads.combine_diad_outputs(*, filename=None, prefix=True, Diad1_fit=None, Diad2_fit=None, to_csv=False, to_clipboard=False, path=None)[source]
This function combines the dataframes from fit_diad_2_w_bck and fit_diad_1_w_bck into a single dataframe
- Parameters:
filename (str) – filename data is for
prefix (bool) – if True, splits filename on ‘ ‘, if False, leaves as is
Diad1_fit (pd.DataFrame) – pd.DataFrame of fit from fit_diad_1_w_bck function
Diad2_fit (pd.DataFrame) – pd.DataFrame of fit from fit_diad_2_w_bck function
to_clipboard (bool) – If True, copies merged dataframe to clipboard
to_csv (bool) – if True, saves fit to Csv, by making a subfolder ‘Diad_Fits’ inside the folder ‘path’
- class DiadFit.diads.diad1_fit_config(model_name: str = 'PseudoVoigtModel', fit_peaks: int = 2, N_poly_bck_diad1: float = 1, lower_bck_diad1: Tuple[float, float] = (1180, 1220), upper_bck_diad1: Tuple[float, float] = (1300, 1350), fit_gauss: bool | None = False, gauss_amp: float | None = 1000, diad_sigma: float = 0.2, diad_sigma_min_allowance: float = 0.2, diad_sigma_max_allowance: float = 5, diad_prom: float = 100, HB_prom: float = 20, x_range_baseline: float = 75, y_range_baseline: float = 100, dpi: float = 200, x_range_residual: float = 20, return_other_params: bool = False)[source]
Configuration class for fitting Diad 1
This class stores configuration parameters for fitting Diad 1 and associated peaks. The class allows for easy customization of the fitting process.
- model_name
Default ‘PseudoVoigtModel’. Model name for peak fitting. Choose from ‘PseudoVoigtModel’, ‘VoigtModel’, ‘GaussianModel’, ‘LorentzianModel’, etc.
- Type:
- N_poly_bck_diad1
Default = 1 Degree of polynomial for background subtraction. Default 1.
- Type:
- lower_bck_diad1
Lower boundary for background fitting in cm-1. Default (1180, 1220)
- upper_bck_diad1
Upper boundary for background fitting in cm-1 Default (1300, 1350).
- fit_gauss
Default True Fit a Gaussian peak if True. Helpful for very elevated spectra
- Type:
Optional[bool]
- gauss_amp
Default = 1000 Gaussian peak amplitude if fit_gauss is True. Default 1000. We find 2X HB intensity is a good guess on many instruments
- Type:
Optional[float]
- diad_sigma_min_allowance
Default = 0.2 Tolerance on entered sigma. Means minimum allowed sigma is 0.2*sigma
- Type:
- diad_sigma_max_allowance
Default = 5 Tolerance on entered sigma. Means minimum allowed sigma is 5*sigma
- Type:
- diad_prom
Default = 100 Peak amplitude for of the diad. Suggest users obtain from the estimated peak parameter dataframe
- Type:
- HB_prom
Default = 20 Peak amplitude for HB. Suggest users obtain from the estimated peak parameter dataframe
- Type:
- HB_sigma_min_allowance
Default = 0.05 Means HB sigma cant be less than 0.05* sigma of first fitting peak (diad).
- Type:
- HB_sigma_max_allowance
Default = 3: Means HB sigma cant exceed 3* sigma of first fitting peak (diad)
- Type:
- HB_amp_min_allowance
Default =0.01 Means HB amplitude cant be less than 0.01*amplitude of first fitted peak (diad).
- Type:
- HB_amp_max_allowance
Default = 1 Means HB amplitude cant be more than 1*amplitude of first fitted peak (diad).
- Type:
- x_range_baseline
Default 75. Means that x axis of baseline shows diad position +- 75 x units.
- Type:
- y_range_baseline
Y-axis range for baseline display. Default 100. Shows a y axis scale that is 100 y units above the minimum baseline value, and 100 units above the maximum value
- Type:
- x_range_residual
Default 20. Shows x values +-20 either side of the diad peak position in the residual plot.
- Type:
- class DiadFit.diads.diad2_fit_config(model_name: str = 'PseudoVoigtModel', fit_peaks: int = 3, N_poly_bck_diad2: float = 1, lower_bck_diad2: Tuple[float, float] = (1300, 1360), upper_bck_diad2: Tuple[float, float] = (1440, 1470), fit_gauss: bool | None = False, gauss_amp: float | None = 1000, diad_sigma: float = 0.2, diad_sigma_min_allowance: float = 0.2, diad_sigma_max_allowance: float = 5, diad_prom: float = 100, HB_prom: float = 20, C13_prom: float = 10, x_range_baseline: float = 75, y_range_baseline: float = 100, plot_figure: bool = True, dpi: float = 200, x_range_residual: float = 20, return_other_params: bool = False)[source]
Configuration class for fitting Diad 2
This class stores configuration parameters for fitting Diad 2 and associated peaks. The class allows for easy customization of the fitting process.
- model_name
Default ‘PseudoVoigtModel’. Model name for peak fitting. Choose from ‘PseudoVoigtModel’, ‘VoigtModel’, ‘GaussianModel’, ‘LorentzianModel’, etc.
- Type:
- fit_peaks
Default 2 Number of peaks to fit. 3 = HB + Diad + C13, 2 = Diad + HB, 1 = Diad.
- Type:
- N_poly_bck_diad2
Default = 1 Degree of polynomial for background subtraction. Default 1.
- Type:
- lower_bck_diad2
Lower boundary for background fitting in cm-1. Default (1300, 1360)
- upper_bck_diad2
Upper boundary for background fitting in cm-1 Default (1440, 1470)
- fit_gauss
Default True Fit a Gaussian peak if True. Helpful for very elevated spectra
- Type:
Optional[bool]
- gauss_amp
Default = 1000 Gaussian peak amplitude if fit_gauss is True. Default 1000. We find 2X HB intensity is a good guess on many instruments
- Type:
Optional[float]
- diad_sigma_min_allowance
Default = 0.2 Tolerance on entered sigma. Means minimum allowed sigma is 0.2*sigma
- Type:
- diad_sigma_max_allowance
Default = 5 Tolerance on entered sigma. Means minimum allowed sigma is 5*sigma
- Type:
- diad_prom
Default = 100 Peak amplitude for of the diad. Suggest users obtain from the estimated peak parameter dataframe
- Type:
- HB_prom
Default = 20 Peak amplitude for HB. Suggest users obtain from the estimated peak parameter dataframe
- Type:
- C13_prom
Default = 10 Peak amplitude for C13. Suggest users obtain from the estimated peak parameter dataframe
- Type:
- HB_sigma_min_allowance
Default = 0.05 Means HB sigma cant be less than 0.05* sigma of first fitting peak (diad).
- Type:
- HB_sigma_max_allowance
Default = 3: Means HB sigma cant exceed 3* sigma of first fitting peak (diad)
- Type:
- HB_amp_min_allowance
Default =0.01 Means HB amplitude cant be less than 0.01*amplitude of first fitted peak (diad).
- Type:
- HB_amp_max_allowance
Default = 1 Means HB amplitude cant be more than 1*amplitude of first fitted peak (diad).
- Type:
- x_range_baseline
Default 75. Means that x axis of baseline shows diad position +- 75 x units.
- Type:
- y_range_baseline
Y-axis range for baseline display. Default 100. Shows a y axis scale that is 100 y units above the minimum baseline value, and 100 units above the maximum value
- Type:
- x_range_residual
Default 20. Shows x values +-20 either side of the diad peak position in the residual plot.
- Type:
- class DiadFit.diads.diad_id_config(height: float = 1, width: float = 0.5, prominence: float = 50, plot_figure: bool = True, exclude_range1: Tuple[float, float] | None = None, exclude_range2: Tuple[float, float] | None = None, approx_diad2_pos: Tuple[float, float] = (1379, 1395), approx_diad1_pos: Tuple[float, float] = (1275, 1290), Diad2_window: Tuple[float, float] = (1349, 1425), Diad1_window: Tuple[float, float] = (1245, 1290), approx_diad2_pos_3peaks: Tuple[float, float] = (1379, 1395, 1362))[source]
- DiadFit.diads.filter_by_string(*, fit_params, data_y_all, x_cord, str_filt=None)[source]
Filters Spectra based on approximate splitting, draws a plot showing spectra to discard and those to keep
- Parameters:
fit_params (pd.dataframe) – dataframe of fit parameters from loop_approx_diad_fits
data_y_all (np.array) – y coordinates of each spectra from loop_approx_diad_fits, used for plotting visualizatoins
x_cord (np.array) – x coordinates of 1 spectra. Assumes all x coordinates the same length
- str_filt: str
Filters just based on string in filename. Keeps files with string in, discards those without.
- Returns:
fit_params_filt (pd.DataFrame) – dataframe of fit parameters for spectra to keep
data_y_filt (np.array) – y coordinates of spectra to keep
fit_params_disc (pd.DataFrame) – dataframe of fit parameters for spectra to discard
data_y_disc (np.array) – y coordinates of spectra to discard
- DiadFit.diads.filter_splitting_prominence(*, fit_params, data_y_all, x_cord, splitting_limits=[100, 107], lower_diad1_prom=10, exclude_str=None, str_filt=None)[source]
Filters Spectra based on approximate splitting, draws a plot showing spectra to discard and those to keep
- Parameters:
fit_params (pd.dataframe) – dataframe of fit parameters from loop_approx_diad_fits
data_y_all (np.array) – y coordinates of each spectra from loop_approx_diad_fits, used for plotting visualizatoins
x_cord (np.array) – x coordinates of 1 spectra. Assumes all x coordinates the same length
splitting_limits (list (e.g. [100, 107])) – only keeps spectra with splitting in this range
lower_diad1_prom (int, float) – Only keeps spectra that meet the splitting parameter, and have an absolute diad1 prominence greater than this value (helps filter out other weird spectra)
exclude_str (str) – Excludes files with this string.
str_filt (str) – Filters just based on string in filename
- Returns:
fit_params_filt (pd.DataFrame) – dataframe of fit parameters for spectra to keep
data_y_filt (np.array) – y coordinates of spectra to keep
fit_params_disc (pd.DataFrame) – dataframe of fit parameters for spectra to discard
data_y_disc (np.array) – y coordinates of spectra to discard
- DiadFit.diads.fit_diad_1_w_bck(*, config1: diad1_fit_config = diad1_fit_config(model_name='PseudoVoigtModel', fit_peaks=2, N_poly_bck_diad1=1, lower_bck_diad1=(1180, 1220), upper_bck_diad1=(1300, 1350), fit_gauss=False, gauss_amp=1000, diad_sigma=0.2, diad_sigma_min_allowance=0.2, diad_sigma_max_allowance=5, diad_prom=100, HB_prom=20, x_range_baseline=75, y_range_baseline=100, dpi=200, x_range_residual=20, return_other_params=False), config2: diad_id_config = diad_id_config(height=1, width=0.5, prominence=50, plot_figure=True, exclude_range1=None, exclude_range2=None, approx_diad2_pos=(1379, 1395), approx_diad1_pos=(1275, 1290), Diad2_window=(1349, 1425), Diad1_window=(1245, 1290), approx_diad2_pos_3peaks=(1379, 1395, 1362)), path=None, filename=None, filetype=None, plot_figure=True, close_figure=True, Diad_pos=None, HB_pos=None, diad_array=None, minimise='weighted_least_squares')[source]
This function fits the background (using the function remove_diad_baseline) and then fits the peaks using fit_gaussian_voigt_generic_diad() It then checks if any parameters are right at the permitted edge (meaning fitting didnt converge), and tries fitting if so. Then it makes a plot of the various parts of the fitting procedure, and returns a dataframe of the fit parameters.
- Parameters:
diad_id_config() – Dataclass of parameters used to look for peaks. Uses any excluded ranges from this.
diad1_fit_config() – Dataclass of parameters used to find peaks. Use help on this dataclass to get options
- path: str
Folder user wishes to read data from
- filename: str
Specific name of file being used.
- filetype: str
Identifies type of file. choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
- plot_figure: bool
if True, saves figure
- Diad_pos: float
Estimate of diad position
- HB_pos: float
Estimate of HB position
- close_fig: bool
if True, displays figure in Notebook. If False, closes it (you can still find it saved)
- minimise: str
‘weighted_least_squares’ or ‘least_squares’
- Returns:
if return_other_params is False (default) – returns a dataframe of peak parameters
if True, also returns – result: fit parameters y_best_fit: best fit to all curves in y x_lin: linspace of best fit coordinates in x.
- DiadFit.diads.fit_diad_2_w_bck(*, config1: diad2_fit_config = diad2_fit_config(model_name='PseudoVoigtModel', fit_peaks=3, N_poly_bck_diad2=1, lower_bck_diad2=(1300, 1360), upper_bck_diad2=(1440, 1470), fit_gauss=False, gauss_amp=1000, diad_sigma=0.2, diad_sigma_min_allowance=0.2, diad_sigma_max_allowance=5, diad_prom=100, HB_prom=20, C13_prom=10, x_range_baseline=75, y_range_baseline=100, plot_figure=True, dpi=200, x_range_residual=20, return_other_params=False), config2: diad_id_config = diad_id_config(height=1, width=0.5, prominence=50, plot_figure=True, exclude_range1=None, exclude_range2=None, approx_diad2_pos=(1379, 1395), approx_diad1_pos=(1275, 1290), Diad2_window=(1349, 1425), Diad1_window=(1245, 1290), approx_diad2_pos_3peaks=(1379, 1395, 1362)), path=None, filename=None, peak_pos_voigt=None, filetype=None, diad_array=None, plot_figure=True, close_figure=False, Diad_pos=None, HB_pos=None, C13_pos=None, minimise='weighted_least_squares')[source]
This function fits the background (using the function remove_diad_baseline) and then fits the peaks using fit_gaussian_voigt_generic_diad() It then checks if any parameters are right at the permitted edge (meaning fitting didnt converge), and tries fitting if so. Then it makes a plot of the various parts of the fitting procedure, and returns a dataframe of the fit parameters.
- Parameters:
diad_id_config() – Dataclass of parameters used to look for peaks. Uses any excluded ranges from this.
diad2_fit_config() – Dataclass of parameters used to find peaks. Use help on this dataclass to get options.
path (str) – Folder user wishes to read data from
filename (str) – Specific name of file being used.
filetype (str) – Identifies type of file. choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
plot_figure (bool) – if True, saves figure
close_fig (bool) – if True, displays figure in Notebook. If False, closes it (you can still find it saved)
Diad_pos (float) – Estimate of diad peak position
HB_pos (float) – Estimate of HB peak position
C13_pos (float) – Estimate of C13 peak position
minimise (str) – ‘weighted_least_squares’ vs ‘least_squares’
- Returns:
if return_other_params is False (default) – returns a dataframe of peak parameters
if True, also returns – result: fit parameters y_best_fit: best fit to all curves in y x_lin: linspace of best fit coordinates in x.
- DiadFit.diads.fit_gaussian_voigt_generic_diad(config1, *, diad1=False, diad2=False, path=None, filename=None, xdat=None, ydat=None, span=None, plot_figure=True, dpi=200, Diad_pos=None, HB_pos=None, C13_pos=None, fit_peaks=None, block_print=True, py_baseline=None, minimise='weighted_least_squares')[source]
This function fits a voigt/Pseudovoigt curves to background subtracted data for your diads +- HBs+-C13 peaks. It can also fit a second Gaussian background iteratively with these peaks
- Parameters:
config1 (fit configuration file, from diad2_fit_config or diad1_fit_config data classes.) – Detailed descriptions of things in this file are found in docs for fit_diad_2_w_bck, and fit_diad_1_w_bck
diad1 (bool) – True if fitting Diad1, False if Diad2
diad2 (bool) – True if fitting Diad2, False if Diad1
path (str) – Path to spectra folder. Used to make a folder inside this folder to store peak fit images in
filename (str) – filename of specific diad being fitted
xdat (np.array) – x coordinates of the spectra
ydat (np.array) – y coordinates of the spectra
span (list (e.g. [1100, 1300])) – Comes from remove_diad_baseline function. Sets min and max x coordinate of new linspace for plotting the overall peak fit over the data
plot_figure (bool) – if True, plots figure from this function specifically. Even if False, if called through fit_diad_2_w_bck or fit_diad_1_w_bck, get a figure from those functions with relevant data on
dpi (int) – dpi of saved figure
Diad_pos (float) – Initial guess of Diad, HB, C13 pos
HB_pos (float) – Initial guess of Diad, HB, C13 pos
C13_pos (float) – Initial guess of Diad, HB, C13 pos
fit_peaks (int) – Number of peaks to fit. 1= just diad, 2 = diad + HB, 3 = diad+ HB +C13
block_print (bool) – If true, prints a few status updates as it goes to help debug
py_baseline – the baseline from the diadbaseline function
squares' (minimise='weighted_least_squares' or 'least) – uses a weighted least squares cost function. Prior to 1.0.19, least squares was default.
- Returns:
result, df_out, y_best_fit, x_lin, components, xdat, ydat, ax1_xlim, ax2_xlim, config1.fit_gauss, residual_diad_coords, ydat_inrange, xdat_inrange
result (lmfit.model.ModelResult (model form lmfit))
df_out (dataframe of peak fit parameters)
y_best_fit (np.array of best peak fit y parameters)
x_lin (np.array of x coordinates corresponding to y_best_fit)
components (y coordinates of the different components of the peak fit.)
xdat (input xdat)
ydat (input ydat)
ax1_xlim, ax2_xlim (xlim for plot in final figure in fit_w_bck for different fit components)
residual_diad_coordinates (np.array of y coordinate of residual plot)
ydat_inrange (np.array of data within ‘span’)
xdat_inrange (np.array of data within ‘span’)
- DiadFit.diads.fit_generic_peak(*, config: generic_peak_config = generic_peak_config(name='generic', lower_bck=(1060, 1065), upper_bck=(1120, 1130), model_name='PseudoVoigtModel', x_range_bck=10, N_poly_carb_bck=1, amplitude=1000, cent=1090, outlier_sigma=12, distance=10, prominence=5, width=6, threshold=0.1, height=100, exclude_range=None, return_other_params=False, N_peaks=1, int_cut_off=0.05), path=None, filename=None, filetype=None, plot_figure=True, dpi=200)[source]
This function fits a generic peak with various options
config: from dataclass generic_peak_config
- model_name: str, ‘GaussianModel’, ‘VoigtModel’, ‘PseudoVoigtModel’, ‘Spline’, ‘Poly’,
Fits Gaussian, Voigt, or PseudoVoigt, or fits a spline or polynomail and gets area underneath
- name: str
Name of feature you are fitting. Used to make column headings in output (e.g., if SO2, outputs are Peak_Cent_SO2, Peak_Area_SO2)
- int_cut_off: float
fits peak out to int_cut_off* peak height.
- lower_bkc: Tuple[float, float]
Upper and lower x coordinate for background to left of peak
- upper_bkc: Tuple[float, float]
Upper and lower x coordinate for background to right of peak
- x_range_bck: float
sets xlim of figure, as the peak center + x_range_bck, and peak cent - x_range_bck
- N_poly_carb_bck: int
Degree of polynomial to fit to the background positions
- amplitude: float
Approx amplitude of peak.
- cent: float
Approx cent of peak
- outlier_sigma: float
Discards points on the background outside outlier_sigma of the median value
- distance, prominence, threshold, height: floats
parameters for scipy find_peaks
- exclude_range: Optional [Tuple[float, float]]
Excludes region of spectra between these values
- dpi: float
dpi to save figure at
- return_other_params: bool
if True, returns other peak fit parameters and fits (df, xx_carb, y_carb, result0), else just returns df of peak positions, heights etc
- path: str
Path to file
- filename: str
filename with file extension
- filetype: str
choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
- plot_figure: bool
If true, plots figure and saves it in a subfolder in the same directory as the filepath specified
- dpi: int
dpi to save figure at
- Return type:
df of peak fits, unless return_other_params is True, then also returns x and y coordinate of fit, and fit params
- class DiadFit.diads.generic_peak_config(name: str = 'generic', lower_bck: Tuple[float, float] = (1060, 1065), upper_bck: Tuple[float, float] = (1120, 1130), model_name: str = 'PseudoVoigtModel', x_range_bck: float = 10, N_poly_carb_bck: float = 1, amplitude: float = 1000, cent: float = 1090, outlier_sigma: float = 12, distance: float = 10, prominence: float = 5, width: float = 6, threshold: float = 0.1, height: float = 100, exclude_range: Tuple[float, float] | None = None, return_other_params: bool = False, N_peaks: int = 1, int_cut_off: float = 0.05)[source]
Configuration class for fitting secondary peaks (e.g. SO2, Carbonate)
- name
Default ‘Generic’ Name of component you are fitting that gets stamped onto the peak parameters, e.g. if ‘SO2’, get Peak_Area_SO2’
- Type:
- lower_bck, upper_bck
Tuple[float, float] Background positions in wavenumbers
- model_name
Default ‘PseudoVoigt’ Type of peak fit. Choose from ‘GaussianModel’, ‘VoigtModel’, ‘PseudoVoigtModel’, ‘Spline’.
- Type:
- x_range_bck
float (default 10) When showing bck collected plot, shows peak center +- this amount
- Type:
- amplitude
float (default 100) Approximate amplitude of peak fit (first guess for Gaussian , voigt and Psuedovoigt needed)
- Type:
- cent
float (default 1090) Approximate peak position of species you are looking for. E.g. 1090 for CO2, 1150 for SO2
- Type:
- outlier_sigma
float (default 12) Code calculates standard deviation of background, and excludes points more than outlier_sigma*this standard deviatoin from nearby background points
- Type:
- distance, prominence, width, threshold, height
float Scipy find peak parameters to get initial fit of peak
- DiadFit.diads.identify_diad_group(*, fit_params, data_y, x_cord, filter_bool, y_fig_scale=0.1, grp_filter='Weak', str_filt=None)[source]
- Sorts diads into two groups. Those meeting the ‘filter_bool’ criteria, and those not
meeting this criteria. Ones meeting the criteria are shown on the left hand plot,
ones not meeting the criteria are shown in the right hand plot
- Parameters:
fit_params (pd.DataFrame) – dataframe of approximate peak parameters from loop_approx_diad_fits
data_y (np.array) – y coordinates of each file in fit_params
x_cord (np.array) – x coordinate of 1 spectra. Assumes all spectra have the same xcoordinate
filter_bool (Pandas.series of bool (e.g. True, False, True, False)) – This is the filter to use to subdivide spectra into two groups.
y_fig_scale (float, int (0.1)) – Determins how far apart the spectra are displayed on the plot. Larger number, the longer the plot is (and easier spectra are to see). The total height of the figure is 2+y_fig_scale*number of spectra
grp_filter (str) – The name of the group where the filter_bool is true. E.g. if you apply a filter to get the weak diads, enter ‘Weak’ here. If filter_bool is to identify Strong spectra, enter ‘Medium-Strong’ here.
- Returns:
pd.DataFrame of fit parameters where filter_bool is True
pd.DataFrame of fit parameters where filter_bool is False
np.array of y coordinates where filter_bool is True
np.array of y coordinates where filter_bool is False
- DiadFit.diads.identify_diad_peaks(*, config: diad_id_config = diad_id_config(height=1, width=0.5, prominence=50, plot_figure=True, exclude_range1=None, exclude_range2=None, approx_diad2_pos=(1379, 1395), approx_diad1_pos=(1275, 1290), Diad2_window=(1349, 1425), Diad1_window=(1245, 1290), approx_diad2_pos_3peaks=(1379, 1395, 1362)), path=None, filename=None, diad_array=None, filetype='Witec_ASCII', plot_figure=True)[source]
This function fits a spline to the spectral data. It then uses Scipy find peaks to identify the diad, HB and C13 peaks. It outputs parameters such as peak positions, peak prominences, and peak height ratios that can be very useful when splitting data into groups. It uses information stored in diad_id_config for most parameters
- Parameters:
config (dataclass) –
This is diad_id_config, which you should edit before you feed into the function to adjust peak finding parameters. E.g. to change prominence of peaks found by scipy, do config=pf.diad_id_config(prominence=12), then input into this function
- Things that can be adjusted in this config file
height: 1 (threshold height for Scipy find peaks) width: 0.5 (threshold width for Scipy find peaks) prominence: 50 (threshold prominence for scipy find peaks) exclude_range1: None (default ) or list - excludes a certain region of the spectra (e.g, [1200, 1250]) exclude_range2: None (default ) or list - excludes a certain region of the spectra (e.g, [1300, 1310]) Diad_window_width: 30 - Window to the left of Diad1 and to the right of Diad2 to look fo rdiad + hotbands. approx_diad2_pos: (1379, 1395) - Region for code to look for diad2 (+-Diad_window_width). approx_diad1_pos: (1275, 1290) - Region for code to look for diad1 (+-Diad_window_width). LH_bck_diad1, LH_bck_diad2, RH_bck_diad1, RH_bck_diad2: Positions either side of peaks to calculate prominences (e.g. LH_bck_diad1=[1180, 1220])
- path: str
Path where spectra is stored
- filename: str
Specific name of file being plotted
- filetype: str
choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
- plot_figure: bool
if True, plots a figure of peaks, and the identified peaks that Scipy found
- Returns:
pd.DataFrame – Dataframe of approximate peak positions, prominences, etc.
np.array of x and y coordinates of spectra
figure if plot_figure is True
- DiadFit.diads.loop_approx_diad_fits(*, spectra_path, config, Diad_Files, filetype, plot_figure)[source]
Uses the identify_diad_peaks function to loop over all files. :param config: See documentation for identify_diad_peaks :type config: dataclass :param spectra_path: Path where spectra is stored :type spectra_path: str :param Diad_Files: List of file names to loop over :type Diad_Files: list :param filetype: choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’,
‘HORIBA_txt’, ‘Renishaw_txt’
- Parameters:
plot_figure (bool) – if True, plots a figure of peaks, and the identified peaks that Scipy found
- Returns:
pd.DataFrame of the approximate peak parameters for each file.
np.array of y coordinates of every spectra.
- DiadFit.diads.loop_diad_skewness(*, Diad_files, path=None, filetype=None, skewness='abs', sort=False, int_cut_off=0.15, config_diad1: diad1_fit_config = diad1_fit_config(model_name='PseudoVoigtModel', fit_peaks=2, N_poly_bck_diad1=1, lower_bck_diad1=(1180, 1220), upper_bck_diad1=(1300, 1350), fit_gauss=False, gauss_amp=1000, diad_sigma=0.2, diad_sigma_min_allowance=0.2, diad_sigma_max_allowance=5, diad_prom=100, HB_prom=20, x_range_baseline=75, y_range_baseline=100, dpi=200, x_range_residual=20, return_other_params=False), config_diad2: diad2_fit_config = diad2_fit_config(model_name='PseudoVoigtModel', fit_peaks=3, N_poly_bck_diad2=1, lower_bck_diad2=(1300, 1360), upper_bck_diad2=(1440, 1470), fit_gauss=False, gauss_amp=1000, diad_sigma=0.2, diad_sigma_min_allowance=0.2, diad_sigma_max_allowance=5, diad_prom=100, HB_prom=20, C13_prom=10, x_range_baseline=75, y_range_baseline=100, plot_figure=True, dpi=200, x_range_residual=20, return_other_params=False), prominence=10, width=1, height=1)[source]
Loops over all supplied files to calculate skewness for multiple spectra using the functions assess_diad1_skewness() and assess_diad2_skewness()
- Parameters:
Diad_Files (list) – List of filenames to loop through
path (str) – Folder user wishes to read data from
filetype (str) – Identifies type of file choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
skewness ('abs' or 'dir') – If ‘abs’, gives absolute skewness (e.g., largest possible number regardless of direction) if ‘dir’ does skewness as one side over the other
- sort: bool
If true, sorts files alphabetically
- file_ext: str
File format. Default txt, could also enter csv etc.
- exclude_str: str
Excludes files with this string in their name. E.g. if exclude_str=’Ne’ it will exclude Ne lines
- exclude_type: str
Excludes files of this type, e.g. exclude_type=’png’ gets rid of image files.
- lower_baseline_diad1, lower_baseline_diad2: list
Region of spectra to use a baseline (LHS) for diad 1 and diad2
- upper_baseline_diad1, upper_baseline_diad2: list
Region of spectra to use a baseline (RHS) for diad 1 and diad 2
- N_poly_bck_diad1, N_poly_bck_diad2: int
Degree of polynomial for background fitting
- skewness: ‘abs’ or ‘dir’
If ‘abs’, gives absolute skewness (e.g., largest possible number regardless of direction) if ‘dir’ does skewness as one side over the other
- int_cut_off: float
Value of intensity at which to calculate the skewness (e.g. 0.15X the peak height).
- Returns:
pd.DataFrame wtih filename, skewness of Diad 1,
the x position of the LH and RH tie point at intensity=int_int_cut_off
- DiadFit.diads.peak_prominence(x, y, peak_position)[source]
This function calculates the approximate prominence of a peak, by finding the local minimum either side of the peak
- Parameters:
x (np.array) – x coordinates (wavenumber)
y (np.array) – y coordinates (intensity)
peak (position: int, float) – Peak position to calculate prominence for
- Returns:
float
- Return type:
prominence of peak in y coordinates.
- DiadFit.diads.plot_diad(*, path=None, filename=None, filetype='Witec_ASCII', Spectra_x=None, Spectra_y=None)[source]
This function makes a plot of the spectra for a specific file to allow visual inspectoin
- Parameters:
Either
path (str) – Folder where spectra is stored
filename (str) – Filename of specific spectra file of interest
filetype (str) – choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
Or – Spectra_x: np.array of x coordinates Spectra_y: np.array of y coordinates
Returns – plots figure with common peaks overlain
- DiadFit.diads.plot_diad_groups(*, x_cord, Weak_np=None, Medium_np=None, Strong_np=None, y_fig_scale=0.5)[source]
This makes a figure showing spectra classified into Weak, Medium and Strong groups on a 3 panel plot
- Parameters:
x_cord (np.array) – x_coordinate for all spectra (must be the same)
Weak_np (np.array) – y coordinates of spectra classified as weak
Medium_np (np.array) – y coordinates of spectra classified as medium
Strong_np (np.array) – y coordinates of spectra classified as strong
y_fig_scale (float) – The total height of the figure is 2+y_fig_scale*total number of spectra entered.
- Return type:
figure showing 3 diad groups.
- DiadFit.diads.plot_peak_params(fit_params, x_param='Diad1_pos', y1_param='approx_split', y2_param='Mean_Valley_prom', y3_param='C13_abs_prom', y4_param='HB2_abs_prom', fill_na=0)[source]
Makes a 2X2 figure of fit parameters, with the same x axis, and 4 different y axis
- Parameters:
fit_params (Pandas DataFrame) – dataframe of approximate fit params from function loop_approx_diad_fits
x_param (str) – parameter you want on the x axis of all plots
y1_param (str) – parameter on y axis of 1st subplot (top left)
y2_param (str) – parameter on y axis of 2nd subplot (top right)
y3_param (str) – parameter on y axis of 3rd subplot (bottom left)
y4_param (str) – parameter on y axis of 4th subplot (bottom right)
fill_na (int) – The integer value to fill Nans with to show up on plots.
- Return type:
4 part figure
- DiadFit.diads.plot_secondary_peaks(*, Diad_Files, path, filetype, xlim_plot=[1040, 1200], xlim_peaks=[1060, 1100], height=100, threshold=0.1, distance=10, prominence=5, width=6, prominence_filter=False, sigma=3, sigma_window=30, find_peaks_filter=False, just_plot=False, yscale=0.2, gaussian_smooth=False, smoothing_window=5, smooth_std=3)[source]
This function plots spectra stacked vertically ontop of each other, in a specific region. It also has two options to identify peak positions.
- Parameters:
Diad_Files (list) – List of filenames to plot
path (str) – Path where spectra are stored
filetype (str) – Choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
xlim_plot (Tuple[float, float]) – Range of x coordinates to make plot over
xlim_peaks (Tuple[float, float]) – Range of x coordinates to search for peaks in
peaks (Optional parameters for Smoothing spectra prior to finding) –
- gaussian_smooth:bool - If True, applies a gaussian smoothing filter using np.convolve defined by a
smoothing_window and a smooth_std.
If you want to actually ID peaks rather than just plot them, choose from:
- prominence_filter: bool
if True, finds max y coordinate in the window xlim_peaks. It then calculates the prominence by comparing this peak height to the average of the median of the 3 datapoints at the left and right of the window. If the peak is more than ‘prominence’ units above the average of these two end point medians, it is classified as a peak. Also need to specify prominence: float
find_peaks_filter: bool
if True, uses scipy find peaks to find peaks. Can tweak distance, prominence, width, threshold, and height.
- just_plot: bool
if True, just plots spectra, doesnt try to find peaks
- y_scale: float
Figure size is 2+y_scale * number of spectra
- Returns:
df_peaks, x_data, y_data, fig
df_peaks (pd.DataFrame for each file of peak position, height, prominence and file name. Filled with Nans if didnt find a peak)
x_data (np.array of x coordinates (Assumes the same for all files))
y_data (np.array of y coordinates of each spectra)
fig (figure it has plotted)
- DiadFit.diads.remove_diad_baseline(*, path=None, filename=None, Diad_files=None, filetype='Witec_ASCII', diad_array=None, exclude_range1=None, exclude_range2=None, N_poly=1, x_range_baseline=10, lower_bck=[1200, 1250], upper_bck=[1320, 1330], sigma=4, plot_figure=True)[source]
This function fits an Nth degree polynomial fitted through spectra points between two background positions. It returns baseline-subtracted data
- Parameters:
path (str) – Folder where spectra is stored
filename (str) – Specific file being read
OR
Diad_files (str) – Filename if in same folder as script.
filetype (str) – Identifies type of file. choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
exclude_range1 (None or list length 2) – Excludes a region, e.g. a cosmic ray
exclude_range2 (None or list length 2) – Excludes a region, e.g. a cosmic ray
sigma (int) – Filters out points that are more than this number of sigma from the median of the baseline.
N_poly (int) – Degree of polynomial to fit to the backgroun
lower_bck (list len 2) – wavenumbers of baseline region to the left of the peak of interest
upper_bck (list len 2) – wavenumbers to the right of the peak of interest
- plot_figure: bool
if True, plots figure, if False, doesn’t.
- Returns:
y_corr: Background subtracted y values trimmed in baseline range x: initial x values trimmed in baseline range Diad_short: Initial data (x and y) trimmed in baseline range Py_base: Fitted baseline for trimmed x coordinates Pf_baseline: polynomial model fitted to baseline Baseline_ysub: Baseline fitted for x coordinates in baseline Baseline_x: x co-ordinates in baseline. Baseline: Filtered to remove points outside sigma*std dev of baseline span: lower and upper x coordinate of region that isnt the baseline
- Return type:
y_corr, Py_base, x, Diad_short, Py_base, Pf_baseline, Baseline_ysub, Baseline_x, Baseline, span
Functions associated with fitting Ne lines
- class DiadFit.ne_lines.Ne_peak_config(model_name: str = 'PseudoVoigtModel', N_poly_pk1_baseline: float = 1, N_poly_pk2_baseline: float = 1, lower_bck_pk1: Tuple[float, float] = (-50, -25), upper_bck1_pk1: Tuple[float, float] = (8, 15), upper_bck2_pk1: Tuple[float, float] = (30, 50), lower_bck_pk2: Tuple[float, float] = (-44.2, -22), upper_bck1_pk2: Tuple[float, float] = (15, 50), upper_bck2_pk2: Tuple[float, float] = (50, 51), peaks_1: float = 2, peaks_2: float = 1, DeltaNe_ideal: float = 330.477634, x_range_baseline_pk1: float = 20, y_range_baseline_pk1: float = 200, x_range_baseline_pk2: float = 20, y_range_baseline_pk2: float = 200, pk1_sigma: float = 0.4, pk2_sigma: float = 0.4, x_range_peak: float = 15, x_range_residual: float = 7, LH_offset_mini: Tuple[float, float] = (1.5, 3), LH_offset_mini2: Tuple[float, float] = None, x_span_pk1: Tuple[float, float] | None = None, x_span_pk2: Tuple[float, float] | None = None)[source]
- class DiadFit.ne_lines.Neon_id_config(exclude_range_1: Tuple[float, float] | None = None, exclude_range_2: Tuple[float, float] | None = None, height: tuple = 10, distance: float = 1, prominence: float = 10, width: float = 1, threshold: float = 0.6, peak1_cent: float = 1117, peak2_cent: float = 1447, n_peaks: float = 6)[source]
- DiadFit.ne_lines.calculate_Ne_line_positions(wavelength=532.05, cut_off_intensity=2000)[source]
Calculates Raman shift for a given laser wavelength of Ne lines, using the datatable from NIST of the observed Ne line emissoin in air and the intensity of each line. Data from https://physics.nist.gov/PhysRefData/ASD/lines_form.html Ticked
- DiadFit.ne_lines.calculate_Ne_splitting(wavelength=532.05, line1_shift=1117, line2_shift=1447, cut_off_intensity=2000)[source]
Calculates ideal splitting in air between lines closest to user-specified line shift. E.g. if user enters 1117 and 1447 line, it looks for the nearest theoretical Ne line position based on your wavelength, and calculates the ideal splitting between these theoretical line positions. This is used to calculate the ideal Ne line splitting for doing the Ne correction routine of Lamadrid et al. (2017)
- Parameters:
- Return type:
df of theoretical splitting, line positions used for this, and entered Ne line positions
- DiadFit.ne_lines.filter_Ne_Line_neighbours(*, df_combo=None, Corr_factor=None, number_av=6, offset=5e-05, file_name_filt=None)[source]
This function discards Ne lines with a correction factor more than offset away from the median value of the N points (number_av) either side of it.
- DiadFit.ne_lines.find_closest(df, line1_shift)[source]
This function finds the closest Raman shift value in the inputted dataframe to the inputted line position
- df: pd.DataFrame
Dataframe of Ne line positions based on laser wavelength from the function calculate_Ne_line_positions, or calculate_Ar_line_positions
- line1_shift: int, float
input line position
- error_pk2
Closest theoretical line position
- DiadFit.ne_lines.find_max_row(df, target_shift, tol=2)[source]
This function is used to find the highest amplitude within a predefined ampl range for finding the right Ne line
- DiadFit.ne_lines.fit_Ne_lines(*, config: Ne_peak_config = Ne_peak_config(model_name='PseudoVoigtModel', N_poly_pk1_baseline=1, N_poly_pk2_baseline=1, lower_bck_pk1=(-50, -25), upper_bck1_pk1=(8, 15), upper_bck2_pk1=(30, 50), lower_bck_pk2=(-44.2, -22), upper_bck1_pk2=(15, 50), upper_bck2_pk2=(50, 51), peaks_1=2, peaks_2=1, DeltaNe_ideal=330.477634, x_range_baseline_pk1=20, y_range_baseline_pk1=200, x_range_baseline_pk2=20, y_range_baseline_pk2=200, pk1_sigma=0.4, pk2_sigma=0.4, x_range_peak=15, x_range_residual=7, LH_offset_mini=(1.5, 3), LH_offset_mini2=None, x_span_pk1=None, x_span_pk2=None), Ne_center_1=1117.1, Ne_center_2=1147, Ne_prom_1=100, Ne_prom_2=200, Ne=None, filename=None, path=None, prefix=True, plot_figure=True, loop=True, save_clipboard=False, close_figure=False, const_params=True, minimise='least_squares')[source]
This function reads in a user file, fits the Ne lines, and if required, saves an image into a new sub folder
- Parameters:
Ne (np.array) – x coordinate (wavenumber) and y coordinate (intensity)
path (filename and) – used to save filename in datatable, and to make a new folder.
prefix (bool) – Whether or not the filename has a prefix
Ne_center_1 (float) – Approximate peak position of the 1st and 2nd Ne line
Ne_center_2 (float) – Approximate peak position of the 1st and 2nd Ne line
Ne_prom_1 (float) – Approximate prominance of the 1st and 2nd Ne line
Ne_prom_2 (float) – Approximate prominance of the 1st and 2nd Ne line
plot_figure (bool) – Plots a figure, nice to inspect fits, makes it slower
loop (bool)
save_clipboard (bool) – Saves results to clipboard if true
close_figure (bool) – Closes figure if True (useful in some editors)
const_params (bool) – If true, means amplitude and peak sigma have to be closer to the guessed parameters (e.g. used after initial fit). E.g. forced within +-20% of estimated peak parameters if True, if false, +-100%.
minimise (str) – ‘weighted_least_squares’ or ‘least_squares’
config parameters from Ne_peak_config():
- model_name: str
allowed_models = “VoigtModel”, “PseudoVoigtModel”, “Pearson4Model”, “SkewedVoigtModel”
- N_poly_pk1_baseline, N_poly_pk1_baseline: int (default 1)
Degree of polynomial to fit to baseline around pk
lower_bck_pk1, upper_bck1_pk1, upper_bck2_pk1: tuple lower_bck_pk2, upper_bck1_pk2, upper_bck2_pk2: tuple
Background positions relative to estimated peak center.
- peaks_1, peaks_2: int
Number of peaks to fit to peak 1. If you need 2 overlapping peaks, put this Ne line as pk1 if not 1, you also need LH_offset_mini=(1.5, 3). Means second peak put between 1.5 and 3 units left of the 1st peak.
- DeltaNe_ideal: float
Ideal distance between two lines calculated for your laser wavelength.
- x_range_baseline_pk1, x_range_baseline_pk2: int, float
test
- y_range_baseline_pk1, y_range_baseline_pk2: int, float
test
- pk1_sigma, pk2_sigma: float (Default 0.4)
Estimated sigma of each peak
- x_range_peak: flt, int, or None
How much to either side of the peak to show on the final peak fitting plot
- x_range_residual: flt
How much either side of the peak is used for calculating residual
- x_span_pk1, x_span_pk2: float, default None
By default, function fits up to background, but can shrink that range here.
- Returns:
if loop is False – return df, Ne_pk1_reg_x_plot, Ne_pk1_reg_y_plot
if loop is True – return df
Also returns figure.
- DiadFit.ne_lines.fit_Ne_pk(x, y_corr, x_span=[-10, 8], Ne_center=1117.1, amplitude=98.1, pk1_sigma=0.28, LH_offset_mini=[1.5, 3], peaks_pk1=2, model_name='PseudoVoigtModel', block_print=True, const_params=True, spec_res=0.4, py_baseline=None, minimise='least_squares')[source]
This function fits the 1117 Ne line as 1 or two voigt peaks
- Parameters:
x (np.array) – x coordinate (wavenumber)
y (np.array) – Background corrected intensiy
x_span (list length 2. Default [-10, 8]) – Span either side of peak center used for fitting, e.g. by default, fits to 10 wavenumbers below peak, 8 above.
Ne_center (float (default=1117.1)) – Center position for Ne line being fitted
amplitude (integer (default = 98)) – peak amplitude
sigma (float (default =0.28)) – sigma of the voigt peak
peaks_pk1 (integer) – number of peaks to fit, e.g. 1, single voigt, 2 to get shoulder peak
LH_offset_mini (list length 2) – Forces second peak to be within -1.5 to -3 from the main peak position.
print_report (bool) – if True, prints fit report.
py_baseline – Output from the remove_ne_baseline function - needed for weighted least squares
minimise (str) – ‘weighted_least_squares’ or ‘least_squares’
- DiadFit.ne_lines.generate_Ne_corr_model(*, time, Ne_corr, N_poly=3, CI=0.67, bootstrap=False, std_error=True, N_bootstrap=500, save_fig=False, pkl_name='polyfit_data.pkl')[source]
This function takes time stamp and Ne correctoin data to make a predictive polynomial, which it then saves as a pkl file
- Parameters:
time (pd.Series) – Time through the run
Ne_corr (df or pd.Series) – if dataframe, has to have the column ‘Ne_corr’ Else, just Ne correction factor as a pd.Series itself
CI (float. Default 0.67) – Confidence interval to save as uncertainty on model.
Either
std_error (bool (True)) – calculates uncertainty on model using CI
Or
boostrap (bool (False)) – Used for testing of method, keep as False
save_fig (bool) – Saves figure.
pkl_name (str) – Name of model that is saved.
- DiadFit.ne_lines.identify_Ne_lines(*, config: Neon_id_config = Neon_id_config(exclude_range_1=None, exclude_range_2=None, height=10, distance=1, prominence=10, width=1, threshold=0.6, peak1_cent=1117, peak2_cent=1447, n_peaks=6), path=None, filename=None, filetype=None, plot_figure=True, print_df=False, Ne_array=None)[source]
Loads Ne line, uses scipy find peaks to identify peaks, overlays these, and returns approximate peak positions, prominences etc to feed into fitting algorithms
- Parameters:
config (from Neon_id_config) –
This is used to identify peaks using Scipy find peaks. Parameters that can be tweaked exclude_range_1: None, or Tuple[float, float]
Range to exclude (e.g, cosmic ray, instrument noise)
- exclude_range_2: None, or Tuple[float, float]
Range to exclude (e.g, cosmic ray, instrument noise)
- height, distance, prominence, width, threshold: float
Scipy find peak parameters you can tweak
- peak1_cent: float
Estimate of location of Ne line 1
- peak2_cent: float
Estimate of location of Ne line 2
- n_peaks: float
Looks through the largest N peaks of scipy in the entire spectra and identifies them on the plot
path (str) – Folder user wishes to read data from
filename (str) – Specific file being read
filetype (str) – choose from ‘Witec_ASCII’, ‘headless_txt’, ‘headless_csv’, ‘head_csv’, ‘Witec_ASCII’, ‘HORIBA_txt’, ‘Renishaw_txt’
plot_figure (bool) – If True, plots a figure highlighting the identified peaks
print_df (bool) – if True, prints the positions of the N biggest peaks it found
Ne_array (np.array) – Can also enter data as a numpy array, rather than as a filename, filepath and filetype
- Returns:
Ne, df_fit_params
Ne (np.array of spectral data (with ranges excluded))
df_fit_params (DataFrame of approximate peak positions, prominences etc.)
- DiadFit.ne_lines.reg_Ne_lines_time(df, fit='poly', N_poly=None, spline_fit=None)[source]
- Parameters:
df (pd.DataFrame) – dataframe of stitched Ne fits and metadata information from WITEC, must have columns ‘sec since midnight’ and ‘Ne_Corr’
fit (float 'poly', or 'spline') –
- If ‘poly’:
N_poly: int, degree of polynomial to fit (1 if linear)
- if ‘spline’:
spline_fit: The string has to be one of:
‘linear’, ‘nearest’, ‘nearest-up’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, ‘previous’. Look up documentation for interpld
- Returns:
figure of fit and data used to make it
Pf (fit model, can be used to evaluate unknown data (only within x range of df for spline fits).)
- DiadFit.ne_lines.remove_Ne_baseline_pk(Ne, N_poly_pk1_baseline=None, Ne_center_1=None, lower_bck=None, upper_bck1=None, upper_bck2=None, sigma_baseline=None)[source]
This function uses a defined range of values to fit a baseline of Nth degree polynomial to the baseline around a specified peak
- Parameters:
Ne (np.array) – np.array of x and y coordinates from the spectra
N_poly_pk1_baseline (int) – Degree of polynomial used to fit the background
Ne_center_1 (float) – Center position for Ne line being fitted
lower_bck (list (length 2). default [-50, -20]) – position used for lower background relative to peak, so =[-50, -20] takes a background -50 and -20 from the peak center
upper_bck1 (list (length 2). default [8, 15]) – position used for 1st upper background relative to peak, so =[8, 15] takes a background +8 and +15 from the peak center
upper_bck2 (list (length 2). default [30, 50]) – position used for 2nd upper background relative to peak, so =[30, 50] takes a background +30 and +50 from the peak center
- Returns:
y_corr, Py_base, x, Ne_short, Py_base, Baseline_y, Baseline_x
y_corr (numpy.ndarray) (The corrected y-values after subtracting the fitted polynomial baseline from the original data.)
Py_base (numpy.ndarray) (The y-values of the fitted polynomial baseline.)
x (numpy.ndarray) (The x-values of the trimmed data within the specified range.)
Ne_short (numpy.ndarray) (The trimmed data within the specified range.)
Baseline_y (numpy.ndarray) (The y-values of the baseline data points.)
Baseline_x (numpy.ndarray) (The x-values of the baseline data points)
Functions associated with fitting Ar lines
- DiadFit.argon_lines.calculate_Ar_line_positions(wavelength=532.05, cut_off_intensity=2000)[source]
Calculates Raman shift for a given laser wavelength of Ne lines, using the datatable from NIST of Ne line emissoin in air and the intensity of each line.
- DiadFit.argon_lines.calculate_Ar_splitting(wavelength=532.05, line1_shift=1117, line2_shift=1447, cut_off_intensity=2000)[source]
Calculates ideal splitting in air between lines closest to user-specified line shift. E.g. if user enters 1117 and 1447 line, it looks for the nearest theoretical Ne line position based on your wavelength, and calculates the ideal splitting between these theoretical line positions. This is used to calculate the ideal Ne line splitting for doing the Ne correction routine of Lamadrid et al. (2017)
- Parameters:
- Return type:
df of theoretical splitting, line positions used for this, and entered Ne line positions
- DiadFit.argon_lines.filter_Ar_Line_neighbours(*, df_combo=None, Corr_factor=None, number_av=6, offset=5e-05, file_name_filt=None)[source]
Filters Ar correction factors that deviate by more than offset from the local median of neighboring values (defined by number_av). Optionally excludes specified filenames via file_name_filt.
Functions for different densimeters (to convert splitting to density)
- DiadFit.densimeters.Francis_pureCO2(FDS, FDS_std, uncer_FDS, uncer_FDS_std=0)[source]
Returns density using a densimeter made from a single CO2 standard
- DiadFit.densimeters.apply_and_save_vertical_shift_to_model(*, pickle_in_path, new_x, new_y, pickle_out_path=None)[source]
Applies a vertical shift to a saved piecewise model based on new_x and new_y, then saves the shifted model to a new .pkl file.
- Parameters:
pickle_in_path (str) – Path to the original .pkl file (output from build_piecewise_poly_by_density).
new_x (array-like) – Corrected splitting values (x).
new_y (array-like) – Measured density values (y).
pickle_out_path (str, optional) – Where to save the new model. If None, appends ‘_shifted.pkl’ to the input path.
- Returns:
shift – Vertical shift applied to the model.
- Return type:
- DiadFit.densimeters.apply_and_save_vertical_shift_to_ucb_densimeter(new_x, new_y)[source]
Applies a vertical shift to a saved piecewise model based on new_x and new_y, then saves the shifted model to a new .pkl file in the same directory.
- Parameters:
filename (str) – Name of the original .pkl file (e.g., “smoothed_polyfit_June25_UCB.pkl”).
new_x (array-like) – Corrected splitting values (x).
new_y (array-like) – Measured density values (y).
pickle_out_name (str, optional) – Filename to save the new shifted model. If None, appends ‘_shifted.pkl’ to the input name.
- Returns:
shift – Vertical shift applied to the model.
- Return type:
- DiadFit.densimeters.build_piecewise_poly_by_density(x, y, y_bounds=(0.17, 0.65), degrees=(1, 3, 2), blend_width=0.05, save_path=None)[source]
Fits and optionally saves a smoothed piecewise polynomial model.
- Returns:
callable model_data : dict (can be pickled)
- Return type:
f_base
- DiadFit.densimeters.calculate_Densimeter_std_err_values(*, pickle_str, corrected_split, corrected_split_err, CI_dens=0.67, CI_split=0.67, str_d='LowD')[source]
This function propagates uncertainty from the densimeter polynomial fit and the overall error on the corrected splitting (from both the peak fitting and the Ne line correction).
- Parameters:
pickle_str (str) – Name of Pickle with regression model for a specific part of the densimeter. Need to be in the same path as the notebook you are calling this in.
corrected_split (pd.Series) – panda series of corrected splitting (cm-1)
corrected_split_err (pd. Series) – panda series of error on corrected splitting (contributions from both peak fitting and the Ne correction model if relevant)
str_d (str) – string of what density equation it came from, appended onto column headings.
CI_split (float) – confidence interval for splitting propagation. Should be set the same as CI_
CI_split
- Returns:
Dataframe of preferred density, and error from each source of uncertainty (Density_σ_dens=error from densimeter, Density_σ_split = error from spliting).
- Return type:
pd.DataFrame
- DiadFit.densimeters.calculate_Densimeter_std_err_values_smooth(*, model_data, corrected_split, corrected_split_err, CI_dens=0.67, CI_split=0.67, str_d='Smoothed', x=None, y=None)[source]
Calculates propagated uncertainty for a smoothed polynomial model.
- Parameters:
model_data (dict) – Dictionary from build_piecewise_poly_by_density including coeffs, blend_width, etc.
corrected_split (pd.Series or np.ndarray) – Corrected splitting values
corrected_split_err (float or pd.Series) – Uncertainty on splitting
CI_dens (float) – Confidence interval for uncertainty in the fit
CI_split (float) – Confidence interval for splitting uncertainty
str_d (str) – Prefix for column names
x (array-like (optional)) – Original data used to fit the model, if not included in model_data
y (array-like (optional)) – Original data used to fit the model, if not included in model_data
- Returns:
DataFrame of predicted values and propagated uncertainties
- Return type:
pd.DataFrame
- DiadFit.densimeters.calculate_dens_error(temp, Split, split_err)[source]
This function is used by the function calculate_density_cornell_old to calculate a max and min density corresponding to splitting values + - 1 sigma
- Parameters:
- Returns:
Dataframe with column for max calculated density and min calculated density.
- Return type:
pd.DataFrame
- DiadFit.densimeters.calculate_density_cornell(*, lab='CMASS', df_combo=None, temp='SupCrit', CI_split=0.67, CI_neon=0.67, Ne_pickle_str=None, pref_Ne=None, Ne_err=None, corrected_split=None, split_err=None)[source]
This function converts Diad Splitting into CO$_2$ density using the Cornell densimeters. Use lab=’CCMR’ for CCMR and lab=’CMASS’ for Esteban Gazels lab.
- Parameters:
lab (str. 'CMASS' or 'CCMR') – Name of the lab where the analy
Either
df_combo (pandas DataFrame) – data frame of peak fitting information
Or
corrected_split (pd.Series) – Corrected splitting (cm-1)
temp (str) – ‘SupCrit’ if measurements done at 37C ‘RoomT’ if measurements done at 24C - Not supported yet but could be added if needed.
CI_neon (float) – Default 0.67. Confidence interval to use, e.g. 0.67 returns 1 sigma uncertainties. If you use another number, note the column headings will still say sigma.
CI_split (float) – Default 0.67. Confidence interval to use, e.g. 0.67 returns 1 sigma uncertainties. If you use another number, note the column headings will still say sigma.
Either
- Ne_pickle_str: str
Name of Ne correction model
OR
- pref_Ne, Ne_err: float, int
For quick and dirty fitting can pass a preferred value for your instrument before you have a chance to regress the Ne lines (useful when first analysing new samples. )
- Returns:
Prefered Density (based on different equatoins being merged), and intermediate calculations
- Return type:
pd.DataFrame
- DiadFit.densimeters.calculate_density_cornell_old(*, temp='SupCrit', Split, split_err=None)[source]
This function converts Diad Splitting into CO$_2$ density using the densimeters of DeVitre et al. (2021) This should only be used for the Cornell Raman, not other Ramans at present. This is an older version of the function calculate_density_cornell - it does not have the same error propagation capabilities, but we keep it here for backwards compatability.
- Parameters:
- Returns:
Prefered Density (based on different equations being merged), and intermediate calculations
- Return type:
pd.DataFrame
- DiadFit.densimeters.calculate_density_ucb(*, Ne_line_combo='1117_1447', df_combo=None, temp='SupCrit', CI_split=0.67, CI_neon=0.67, Ne_pickle_str=None, Ar_pickle_str=None, pref_Ne=None, Ne_err=None, corrected_split=None, split_err=None)[source]
This function converts Diad Splitting into CO$_2$ density using the UC Berkeley calibration line developed by DeVitre and Wieser in 2023.
- Parameters:
Ne_line_combo (str, '1117_1447', '1117_1400', '1220_1447', '1220_1400', '1220_1567') – Combination of Ne lines used for drift correction
Either
df_combo (pandas DataFrame) – data frame of peak fitting information
Or
corrected_split (pd.Series) – Corrected splitting (cm-1)
temp (str) – ‘SupCrit’ if measurements done at 37C ‘RoomT’ if measurements done at 24C - Not supported yet but could be added if needed.
CI_neon (float) – Default 0.67. Confidence interval to use, e.g. 0.67 returns 1 sigma uncertainties. If you use another number, note the column headings will still say sigma.
CI_split (float) – Default 0.67. Confidence interval to use, e.g. 0.67 returns 1 sigma uncertainties. If you use another number, note the column headings will still say sigma.
Either
- Ne_pickle_str: str
Name of Ne correction model
OR
- pref_Ne, Ne_err: float, int
For quick and dirty fitting can pass a preferred value for your instrument before you have a chance to regress the Ne lines (useful when first analysing new samples. )
- Returns:
Prefered Density (based on different equatoins being merged), and intermediate calculations
- Return type:
pd.DataFrame
- DiadFit.densimeters.calculate_density_ucb_new(*, df_combo=None, temp='SupCrit', CI_split=0.67, CI_neon=0.67, Ne_pickle_str=None, pref_Ne=None, Ne_err=None, corrected_split=None, split_err=None, shift=0)[source]
This function converts Diad Splitting into CO$_2$ density using the UC Berkeley calibration line developed by DeVitre and Wieser in 2023.
- Parameters:
Ne_line_combo (str, '1117_1447', '1117_1400', '1220_1447', '1220_1400', '1220_1567') – Combination of Ne lines used for drift correction
Either
df_combo (pandas DataFrame) – data frame of peak fitting information
Or
corrected_split (pd.Series) – Corrected splitting (cm-1)
temp (str) – ‘SupCrit’ if measurements done at 37C ‘RoomT’ if measurements done at 24C - Not supported yet but could be added if needed.
CI_neon (float) – Default 0.67. Confidence interval to use, e.g. 0.67 returns 1 sigma uncertainties. If you use another number, note the column headings will still say sigma.
CI_split (float) – Default 0.67. Confidence interval to use, e.g. 0.67 returns 1 sigma uncertainties. If you use another number, note the column headings will still say sigma.
Either
- Ne_pickle_str: str
Name of Ne correction model
OR
- pref_Ne, Ne_err: float, int
For quick and dirty fitting can pass a preferred value for your instrument before you have a chance to regress the Ne lines (useful when first analysing new samples. )
- Returns:
Prefered Density (based on different equatoins being merged), and intermediate calculations
- Return type:
pd.DataFrame
- DiadFit.densimeters.calculate_errors_no_densimeter(*, df_combo, Ne_pickle_str='polyfit_data.pkl', temp='SupCrit', split_err=0, CI_split=0.67, CI_neon=0.67)[source]
This function calculates the error from just the Ne line correction method. It is largely redundant, still in use for some older projects using an old workflow
- DiadFit.densimeters.merge_fit_files(path)[source]
This function merges the files Discarded_df.xlsx, Weak_Diads.xlsx, Medium_Diads.xlsx, Strong_Diads.xlsx if they exist into one combined dataframe.
- Parameters:
path (str) – Path on computer where these files are stored
- Returns:
pandas dataframe where 3 sets of peak fits are merged.
- Return type:
pd.DataFrame
- DiadFit.densimeters.merge_in_carb_SO2(df_combo, file1_name='Carb_Peak_fits.xlsx', file2_name='SO2_Peak_fits.xlsx', prefix=False, str_prefix=' ', file_ext='.txt')[source]
This function checks for .xlsx files with secondary phases in the path with names file1_name and file2_name, if they are there it will merge them together into a dataframe.
It then merges this with the dataframe ‘df_combo’ of your other fits
- Parameters:
df_combo (pd.DataFrame) – Dataframe of peak fitting parameters for diads
file1_name (str) – Name of first excel spreadsheet of secondary phases to merge in
file2_name (str) – Name of second excel spreadsheet of secondary phases to merge in
prefix (bool) – If True, removes prefix followed by str_prefix from file name, e.g. if your spectra are called 01 SO2_acquisition, the file name would become SO2_acquisition for easier merging.
file_ext (str) – Removes from file name for merging dataframes (as above).
- Returns:
Dataframe of combined peak parameters from secondary phases and diads.
- Return type:
pd.DataFrame
- DiadFit.densimeters.propagate_error_split_argon_peakfit(*, df_fits, Ar_corr=None, Ar_err=None, pref_Ar=None)[source]
This function propagates errors in your Ar correction model and peak fits by quadrature.
- Parameters:
df_fits (pd.DataFrame) – Dataframe of peak fitting parameters. Must contain columns for ‘Diad1_cent_err’, ‘Diad2_cent_err’, ‘Splitting’
either (Choose)
Ar_corr (pd.DataFrame (Optional)) – Dataframe with columns for ‘upper_values’ and ‘lower values’, e.g. the upper and lower bounds of the error on the Ar correction model
Or
Ar_err (pref_Ar and) – A preferred value of the Ar correction factor and the error (e.g. pref_Ar=0.998, Ar_err=0.001). Used for rapid peak fitting before developing Ar lines.
- Returns:
two pd.Series
- Return type:
the error on the splitting, and the combined error from the splitting and the Ar correction model.
- DiadFit.densimeters.propagate_error_split_neon_peakfit(*, df_fits, Ne_corr=None, Ne_err=None, pref_Ne=None)[source]
This function propagates errors in your Ne correction model and peak fits by quadrature.
- Parameters:
df_fits (pd.DataFrame) – Dataframe of peak fitting parameters. Must contain columns for ‘Diad1_cent_err’, ‘Diad2_cent_err’, ‘Splitting’
either (Choose)
Ne_corr (pd.DataFrame (Optional)) – Dataframe with columns for ‘upper_values’ and ‘lower values’, e.g. the upper and lower bounds of the error on the Ne correction model
Or
Ne_err (pref_Ne and) – A preferred value of the Ne correction factor and the error (e.g. pref_Ne=0.998, Ne_err=0.001). Used for rapid peak fitting before developing Ne lines.
- Return type:
two pd.Series, the error on the splitting, and the combined error from the splitting and the Ne correction model.
Functions for propagating uncertainty
- DiadFit.error_propagation.calculate_temperature_density_MC(sample_i=0, N_dup=1000, CO2_dens_gcm3=None, error_CO2_dens=0, error_type_CO2_dens='Abs', error_dist_CO2_dens='normal', T_K=None, error_T_K=0, error_type_T_K='Abs', error_dist_T_K='normal', crust_dens_kgm3=None, error_crust_dens=0, error_type_crust_dens='Abs', error_dist_crust_dens='normal', XH2O=None, error_XH2O=None, error_type_XH2O='Abs', error_dist_XH2O='normal', model=None, neg_values=True)[source]
This function generates the range of T_K, CO2 densities, XH2O and crustal densities for 1 sample for performing Monte Carlo simulations using the function propagate_FI_uncertainty (e.g. this function makes the range of input parameters for each sample, but doesnt do the EOS calculations).
- Parameters:
sample_i (int) – The index of the sample
N_dup (int) – Number of synthetic inputs to do for each sample.
densities (Required input information about CO2)
CO2_dens_gcm3 (pd.Series, integer or float) – CO2 densities in g/cm3 to perform calculations with. Can be a column from your dataframe (df[‘density_g_cm3’]), or a single value (e.g.., 0.2)
error_CO2_dens (float) – Error in CO2 fluid density
error_type_CO2_dens (str) – Type of CO2 fluid density error. Can be ‘Abs’ or ‘Perc’
error_dist_CO2_dens (str) – Distribution of CO2 fluid density error. Can be ‘normal’ or ‘uniform’.
temperature (Required input information about fluid)
T_K (pd.Series, integer, float) – Temperature in Kelvin at which you think your fluid inclusion was trapped. Can be a column from your dataframe (df[‘T_K’]), or a single value (e.g.., 1500)
error_T_K (float) – Error in temperature.
error_type_T_K (str) – Type of temperature error. Can be ‘Abs’ or ‘Perc’.
error_dist_T_K (str) – Distribution of temperature error. Can be ‘normal’ or ‘uniform’.
either (Required input information for converting pressure to depth. Choose)
density (A fixed crustal) –
- crust_dens_kgm3: float
Density of the crust in kg/m^3.
- error_crust_dens: float, optional
Error in crust density.
- error_type_crust_dens: str
Type of crust density error. Can be ‘Abs’ or ‘Perc’.
- error_dist_crust_dens: str
Distribution of crust density error. Can be ‘normal’ or ‘uniform’.
OR a crustal density model
- model: str
see documentation for the function ‘convert_pressure_to_depth’ to see more detail about model options. If you select a model, it will just use this to calculate a depth, but it wont add any uncertainty to this model (as this is more likely systematic than random uncertainty that can be simulated with MC methods) For this function, you dont need any more info like d1, rho1, that comes in the propagate_FI_uncertainty function.
- Returns:
df_out – DataFrame containing information on temperature, CO2 density, crust density, and error.
- Return type:
pandas.DataFrame
- DiadFit.error_propagation.convert_co2_dens_press_depth(EOS='SW96', T_K=None, CO2_dens_gcm3=None, crust_dens_kgm3=None, output='kbar', g=9.81, model=None, XH2O=None, Hloss=True, d1=None, d2=None, d3=None, rho1=None, rho2=None, rho3=None, rho4=None, T_K_ambient=310.15)[source]
This function calculates pressure and depth based on input CO2 densities, temperatures, and crustal density information from the user
- Parameters:
EOS (str) – ‘SP94’ or ‘SW96’ - CO2 equation of state choosen
T_K (float, pd.Series) – Temperature in Kelvin at which fluid inclusion was trapped.
CO2_dens_gcm3 (float, pd.Series) – CO2 density of FI in g/cm3
choose (For Pressure to depth conversion)
crust_dens_kgm3 (float) – Crustal density in kg/m3
OR
model (str) –
- see documentation for the function ‘convert_pressure_to_depth’ to see more detail about model options.
If you select a model, it will just use this to calculate a depth, but it wont add any uncertainty to this model (as this is more likely systematic than random uncertainty that can be simulated with MC methods)
if model is two-step: If two step, must also define:
d1: Depth to first transition in km rho1: Density between surface and 1st transition d2: Depth to second transition in km (from surface) rho2: Density between 1st and 2nd transition
if model is three-step: If three step, must also define:
rho1: Density between surface and 1st transition d1: Depth to first transition in km rho2: Density between 1st and 2nd transition d2: Depth to second transition in km (from surface) rho3: Density between 2nd and 3rd transition depth. d3: Depth to third transition in km (from surface) rho4: Density between below d3
- Returns:
dataframe of pressure, depth, and input parameterss
- Return type:
pd.DataFrame
- DiadFit.error_propagation.loop_all_FI_MC(sample_ID, CO2_dens_gcm3, T_K, N_dup=1000, crust_dens_kgm3=None, d1=None, d2=None, rho1=None, rho2=None, rho3=None, error_crust_dens=0, error_type_crust_dens='Abs', error_dist_crust_dens='uniform', error_T_K=0, error_type_T_K='Abs', error_dist_T_K='normal', error_CO2_dens=0, error_type_CO2_dens='Abs', error_dist_CO2_dens='normal', plot_figure=False, fig_i=0)[source]
This is a redundant function, kept around for a while for backwards compatability
- DiadFit.error_propagation.propagate_FI_uncertainty(sample_ID, CO2_dens_gcm3, T_K, multiprocess='default', cores='default', EOS='SW96', N_dup=1000, plot_figure=False, fig_i=0, error_CO2_dens=0, error_type_CO2_dens='Abs', error_dist_CO2_dens='normal', crust_dens_kgm3=None, error_crust_dens=0, error_type_crust_dens='Abs', error_dist_crust_dens='uniform', model=None, d1=None, d2=None, rho1=None, rho2=None, rho3=None, error_T_K=0, error_type_T_K='Abs', error_dist_T_K='normal', XH2O=None, error_XH2O=0, error_type_XH2O='Abs', error_dist_XH2O='normal', Hloss=True, neg_values=True)[source]
This function performs Monte Carlo simulations of uncertainty in CO2 density, input temperature,
crustal density and XH2O. If XH2O is specified as not None, it will use Duan and Zhang 2006 EOS
It uses the function ‘calculate_temperature_density_MC’ to make the simulated variables, and then uses this to calculate a resulting pressure using the equation of state of choice
- multiprocess: ‘default’ or bool
Default uses multiprocessing for Duan and Zhang (2006) but not for Span and Wanger or Sterner and Pitzer. This is because these EOS are so fast, making the multiprocess has a time penalty. You can override this defualt by specifying True or False.
- cores: ‘default’
By default, if multiprocess, uses default number of cores selected by multiprocess, ca noverride.
- sample_ID: pd.Series
Panda series of sample names. E.g., select a column from your dataframe (df[‘sample_name’])
- N_dup: int
Number of Monte Carlo simulations to perform for each sample
- EOS: str
‘SW96’ or ‘SP94’ for the pure CO2 EOS
- plot_figure: bool
if True, plots a figure for one sample showing the range of input parameters. If this is True, also select: fig_i: int Which determins which sample is plotted (e.g. 0 is the 1st sample name, etc. )
- CO2_dens_gcm3: pd.Series, integer or float
CO2 densities in g/cm3 to perform calculations with. Can be a column from your dataframe (df[‘density_g_cm3’]), or a single value (e.g.., 0.2)
- error_CO2_dens: float
Error in CO2 fluid density
- error_type_CO2_dens: str
Type of CO2 fluid density error. Can be ‘Abs’ or ‘Perc’
- error_dist_CO2_dens: str
Distribution of CO2 fluid density error. Can be ‘normal’ or ‘uniform’.
- T_K: pd.Series, integer, float
Temperature in Kelvin at which you think your fluid inclusion was trapped. Can be a column from your dataframe (df[‘T_K’]), or a single value (e.g.., 1500)
- error_T_K: float
Error in temperature.
- error_type_T_K: str
Type of temperature error. Can be ‘Abs’ or ‘Perc’.
- error_dist_T_K: str
Distribution of temperature error. Can be ‘normal’ or ‘uniform’.
- XH2O: Default Nan, else pd.Series, integer, float
mol proportion of H2O in the fluid phase
- error_XH2O: float
Error in XH2O
- error_type_XH2O: str
Type of XH2O error. Can be ‘Abs’ or ‘Perc’.
- error_dist_XH2O: str
Distribution of XH2O error. Can be ‘normal’ or ‘uniform’.
For converting pressure to depth in the crust, choose either
A fixed crustal density
- crust_dens_kgm3: float
Density of the crust in kg/m^3.
- error_crust_dens: float, optional
Error in crust density.
- error_type_crust_dens: str
Type of crust density error. Can be ‘Abs’ or ‘Perc’.
- error_dist_crust_dens: str
Distribution of crust density error. Can be ‘normal’ or ‘uniform’.
OR a crustal density model
- model: str
see documentation for the function ‘convert_pressure_to_depth’ to see more detail about model options. If you select a model, it will just use this to calculate a depth, but it wont add any uncertainty to this model (as this is more likely systematic than random uncertainty that can be simulated with MC methods)
if model is two-step: If two step, must also define:
d1: Depth to first transition in km rho1: Density between surface and 1st transition d2: Depth to second transition in km (from surface) rho2: Density between 1st and 2nd transition
if model is three-step: If three step, must also define:
d1: Depth to first transition in km rho1: Density between surface and 1st transition d2: Depth to second transition in km (from surface) rho2: Density between 1st and 2nd transition d3: Depth to third transition in km (from surface) rho3: Density between 2nd and 3rd transition depth.
- neg_values: bool (default True)
if True, negative values of input parameters allowed, if False, makes neg values zero.
df_step, All_outputs,fig if plot_figure is true
- df_step: pd.DataFrame
has 1 row for each entered sample, with mean, median, and standard deviation for each parameter
- All_outputs: pd.DataFrame
has N_dup rows for each input, e.g. full simulation data for each sample. What is being shown in the figure
- fig: figure
Figure of simulated input and output parameters for 1 sample selected by the user.
- DiadFit.error_propagation.propagate_microthermometry_uncertainty(T_h_C, Sample_ID=None, error_T_h_C=0.3, N_dup=1000, error_dist_T_h_C='uniform', error_type_T_h_C='Abs', EOS='SW96', homog_to=None, set_to_critical=False)[source]
- This function propagates the uncertainty in measured temperature values to calculate the density of gas and
liquid CO2 using an equation of state (EOS).
It loops over more than 1 sample, using the function make_error_dist_microthermometry_1sam for each sample to generate the variable input parameters It generates a dataset of temperature measurements with random noise added to them based on the specified distribution and error type, calculates the CO2 density for each temperature value, and returns the mean and standard deviation of the density values.
- Parameters:
T_h_C (numeric or list of numeric values) – The measured temperature(s) of the sample(s) in degrees Celsius.
Sample_ID (str or pandas.Series) – The ID or a pandas Series of IDs of the sample(s). If not provided, the function uses an index number as the ID. Default is None.
error_T_h_C (float or pandas.Series) – The amount of error to add to the temperature measurement. If a pandas.Series is provided, the function takes the right one for each loop. Default value is 0.3.
N_dup (int, optional) – The number of duplicated samples to generate with random noise. Default value is 1000.
error_dist_T_h_C (str, optional) – The distribution of the random noise to be added to the temperature measurement. Can be either ‘normal’ or ‘uniform’. Default value is ‘uniform’.
error_type_T_h_C (str, optional) – The type of error to add to the temperature measurement. Can be either ‘Abs’ or ‘Perc’. Default value is ‘Abs’.
EOS (str, optional) – The equation of state to use for the calculation. Can be either ‘SW96’ or ‘SP94’. Default value is ‘SW96’.
- homog_tostr, optional
The phase to which the CO2 density is homogenized. Can be either ‘Gas’ or ‘Liq’. Default is None.
- set_to_critical: bool
Default False. If true, if you enter T_h_C which exceeds 30.9782 (the critical point of CO2) it replaces your entered Temp with that temp.
- Returns:
A tuple of two pandas DataFrames. The first DataFrame contains the mean and standard deviation of the gas and liquid CO2 density values. The second DataFrame contains the CO2 density values for each temperature value in the input. The output Std_density_Gas_gcm3_from_percentiles is a more representative output of the error than the standard deviation if the output distribution isn’t Gaussian.
- Return type:
pandas.DataFrame, pandas.DataFrame
Functions for performing CO2-H2O EOS calculations
- DiadFit.CO2_EOS.calc_prop_knownP_EOS_DZ2006(*, P_kbar=1, T_K=1200, XH2O=1)[source]
This function calculates the properties (molar volume, density, compressability factor, fugacity, and activity) for a mixed H2O-CO2 fluid using the EOS of Duan and Zhang (2006). It assumes you know P, T, and XH2O.
- DiadFit.CO2_EOS.calculate_CO2_density_homog_T(T_h_C, EOS='SW96', Sample_ID=None, homog_to=None, set_to_critical=False)[source]
Calculates CO2 density for a specified homogenization temperature in Celcius using eq 3.14 and 3.15 from the Span and Wanger (1996) equation of state.
- Parameters:
EOS (str, 'SW96') – Here for consistency with other functions, only supported for SW96
Optional
Sample_ID (int, pd.series) – Sample ID, will append to the final dataframe
homog_to (str ('L', 'V'), pd.series with strings. Optional) – If specified, returns an additional column ‘Bulk Density’ to choose between the liquid and gas.
set_to_critical (bool) – Default False. If true, if you enter T_h_C which exceeds 30.9782 (the critical point of CO2) it replaces your entered Temp with that temp.
- Returns:
colums for Liq_gcm3 and gas_gcm3, bulk density if homog_to specified
- Return type:
pd.DataFrame
- DiadFit.CO2_EOS.calculate_CO2_density_homog_T_SW96_NIST(T_h_C, Sample_ID=None, homog_to=None)[source]
Calculates CO2 density for a specified homogenization temperature in Celcius using the Span and Wanger (1996) equation of state. Parameterizes based on NIST webook. Very similar to above using the root equations (rounding error issues in webbook?)
- Parameters:
- Returns:
colums for Liq_gcm3 and gas_gcm3, bulk density if homog_to specified
- Return type:
pd.DataFrame
- DiadFit.CO2_EOS.calculate_Density_Sterner_Pitzer_1994(T_K, target_pressure_MPa)[source]
This function uses the objective function ‘objective_function_CO2_dens’ above to solve for CO2 density for a given pressure and Temp
- target_pressure_MPa: int, float
Pressure of CO2 fluid in MPa
- Return type:
CO2 density
- DiadFit.CO2_EOS.calculate_P_SP1994(*, T_K=1400, T_h_C=None, phase=None, CO2_dens_gcm3=None, return_array=False)[source]
This function calculates Pressure using Sterner and Pitzer (1994) using either a homogenization temp, or a CO2 density. You must also input a temperature of your system (e.g. 1400 K for a basalt)
- Parameters:
T_K (int, float, pd.Series) – Temperature in Kelvin to find P at (e.g. temp fluid was trapped at)
Either:
- T_h_C: int, float, pd.Series (optional)
homogenization temp during microthermometry
- phase: str
‘Gas’ or ‘Liquid’, the phase your inclusion homogenizes too
Or:
- CO2_dens_gcm3: int, float, pd.Series
Density of your inclusion in g/cm3, e.g. from Raman spectroscopy
- return_array: bool
if True, returns a pd.array not a df.
- Returns:
Has columns for T_K, T_h_C, Liq_density, Gas_density, P_MPa. Non relevant variables filled with NaN
- Return type:
Pandas.DataFrame
- DiadFit.CO2_EOS.calculate_P_for_rho_T(CO2_dens_gcm3, T_K, EOS='SW96', Sample_ID=None)[source]
This function calculates P in kbar for a specified CO2 density in g/cm3, a known T (in K), and a specified EOS
- Parameters:
CO2_dens_gcm3 (int, float, pd.Series, np.array) – CO2 density in g/cm3
T_K (int, float, pd.Series, np.array) – Temperature in Kelvin
EOS (str) – ‘SW96’ for Span and Wanger (1996), ‘SP94’ for Sterner and Pitzer (1994), or ‘DZ06’ for Pure CO2 from Duan and Zhang (2006)
Sample_ID (str, pd.Series) – Sample ID to be appended onto output dataframe
- Returns:
Pressure in kbar, MPa and input parameters
- Return type:
pd.DataFrame
- DiadFit.CO2_EOS.calculate_P_for_rho_T_DZ06(CO2_dens_gcm3, T_K)[source]
This function calculates P in kbar for a specified CO2 density in g/cm3 and a known T (in K) for the pure CO2 Duan and Zhang (2006) EOS (e.g. XH2O=0)
- DiadFit.CO2_EOS.calculate_P_for_rho_T_SP94(T_K, CO2_dens_gcm3, scalar_return=False)[source]
This function calculates P in kbar for a specified CO2 density in g/cm3 and a known T (in K) for the Sterner and Pitzer (1994) EOS
- DiadFit.CO2_EOS.calculate_P_for_rho_T_SW96(CO2_dens_gcm3, T_K)[source]
This function calculates P in kbar for a specified CO2 density in g/cm3 and a known T (in K) for the Span and Wanger (1996) EOS
- DiadFit.CO2_EOS.calculate_Pressure_DZ2006(*, mol_vol=None, density=None, T_K, XH2O)[source]
Used to calculate pressure in a loop for multiple inputs. Den - bulk density.
- Parameters:
mol_vol (molar volume in g/cm3)
density (density in g/cm3)
T_K (entrapment temperature in kelvin)
XH2O (molar fraction of H2O in the fluid)
Returns –
Pressure in bars
- DiadFit.CO2_EOS.calculate_Pressure_ind_DZ2006(*, mol_vol, T_K, XH2O, Pguess=None)[source]
This function calculates pressure (in bars) for a known molar volume, T in K and XH2O (mol frac) for a single value. It uses a look up table to get pressure, then a newton and raphson method (implemented in the function mixpressure) to find the best fit pressure. There are some densities, T_K and XH2O values where the volume is negative.
- DiadFit.CO2_EOS.calculate_SP19942(T_K, target_pressure_MPa)[source]
This function Solves for CO2 density for a given temp and pressure using the objective and minimise functions above.
- DiadFit.CO2_EOS.calculate_SP1994_Temp(CO2_dens_gcm3, target_pressure_MPa)[source]
This function Solves for Temp for a given CO2 density and pressure using the objective and minimise functions above.
- DiadFit.CO2_EOS.calculate_T_for_rho_P(CO2_dens_gcm3, P_kbar, EOS='SW96', Sample_ID=None)[source]
This function calculates P in kbar for a specified CO2 density in g/cm3, a known T (in K), and a specified EOS
- Parameters:
- Returns:
Temp in Kelvin and input parameters
- Return type:
pd.DataFrame
- DiadFit.CO2_EOS.calculate_T_for_rho_P_SP94(P_kbar, CO2_dens_gcm3)[source]
- This function calculates Temp for a known Pressure (in kbar) and a known CO2 density in g/cm3
using the Sterner and Pitzer EOS. it references the objective functions to solve for density.
- P_kbar: int, float, pd.Series, np.array
Pressure in kbar
- CO2_dens_gcm3: int, float, pd.Series, np.array
CO2 density in g/cm3
- pd.Series
Temp in K
- DiadFit.CO2_EOS.calculate_T_for_rho_P_SW96(CO2_dens_gcm3, P_kbar)[source]
This function calculates Temperature for a known CO2 density in g/cm3 and a known Pressure (in kbar) for the Span and Wagner (1996) EOS
- DiadFit.CO2_EOS.calculate_Temp_Sterner_Pitzer_1994(CO2_dens_gcm3, target_pressure_MPa)[source]
This function uses the objective function above to solve for Temp for a given pressure and CO2 density
- DiadFit.CO2_EOS.calculate_entrapment_P_XH2O(*, XH2O, CO2_dens_gcm3, T_K, T_K_ambient=310.15, fast_calcs=False, Hloss=True)[source]
“ This function calculates pressure for a measured CO$_2$ density, temperature and estimate of initial XH2O. It first corrects the density to obtain a bulk density for a CO2-H2O mix, assuming that H2O was lost from the inclusion. correcting for XH2O. It assumes that H2O has been lost from the inclusion (see Hansteen and Klugel, 2008 for method). It also calculates using other pure CO2 equation of states for comparison
- Parameters:
XH2O (float, pd.Series.) – The molar fraction of H2O in the fluid. Should be between 0 and 1. Can get an estimate from say VESical.
CO2_dens_gcm3 (float, pd.Series) – Measured CO2 density in g/cm3
T_K (float, pd.Series) – Temperature in Kelvin fluid was trapped at
T_K_ambient (pd.Series) – Temperature in Kelvin Raman measurement was made at.
fast_calcs (bool (default False)) – If True, only performs one EOS calc for DZ06, not 4 (with water, without water, SP94 and SW96). also specify H2Oloss=True or False
- Returns:
if fast_calcs is False
pd.DataFrame – Columns showing: P_kbar_pureCO2_SW96: Pressure calculated for the measured CO$_2$ density using the pure CO2 EOS from Span and Wanger (1996) P_kbar_pureCO2_SP94: Pressure calculated for the measured CO$_2$ density using the pure CO2 EOS from Sterner and Pitzer (1994) P_kbar_pureCO2_DZ06: Pressure calculated from the measured CO$_2$ density using the pure CO2 EOs from Duan and Zhang (2006) P_kbar_mixCO2_DZ06_Hloss: Pressure calculated from the reconstructed mixed fluid density using the mixed EOS from Duan and Zhang (2006) assuming H loss P_kbar_mixCO2_DZ06_noHloss: Pressure calculated from the reconstructed mixed fluid density using the mixed EOS from Duan and Zhang (2006) assuming H loss P Mix_Hloss/P Pure DZ06: Correction factor - e.g. how much deeper the pressure is from the mixed EOS with H loss P Mix_noHloss/P Pure DZ06: Correction factor - e.g. how much deeper the pressure is from the mixed EOS (assuming no H loss) rho_mix_calc_noHloss: Bulk density calculated (C+H) rho_mix_calc_Hloss: Bulk density calculated (C+H) after h loss CO2_dens_gcm3: Input CO2 density T_K: input temperature XH2O: input molar fraction of H2O
if fast_calcs is True – P_kbar_mixCO2_DZ06: Pressure calculated from the reconstructed mixed fluid density using the mixed EOS from Duan and Zhang (2006)
- DiadFit.CO2_EOS.calculate_molar_volume_DZ2006(*, P_kbar, T_K, XH2O)[source]
Used to calculate molar volume (cm3/mol) in a loop for multiple inputs
- DiadFit.CO2_EOS.calculate_molar_volume_ind_DZ2006(*, P_kbar, T_K, XH2O)[source]
This function calculates molar volume (cm3/mol) for a known pressure (kbar), T in K and XH2O (mol frac) for a single value
- DiadFit.CO2_EOS.calculate_rho_for_P_T(P_kbar, T_K, EOS='SW96')[source]
This function calculates CO2 density in g/cm3 for a known Pressure (in kbar), a known T (in K), and a specified EOS
- DiadFit.CO2_EOS.calculate_rho_for_P_T_H2O(P_kbar, T_K)[source]
This function calculates H2O density in g/cm3 for a known Pressure (in kbar), a known T (in K) using the Wanger and Pru (2002) EOS from CoolProp doi:10.1063/1.1461829.
- DiadFit.CO2_EOS.calculate_rho_for_P_T_SP94(P_kbar, T_K)[source]
This function calculates CO2 density in g/cm3 for a known Pressure (in kbar), a known T (in K) for the Sterner and Pitzer EOS it references the objective functions to solve for density.
- DiadFit.CO2_EOS.calculate_rho_for_P_T_SW96(P_kbar, T_K)[source]
This function calculates CO2 density in g/cm3 for a known Pressure (in kbar), a known T (in K) for the Span and Wagner (1996) EOS
- DiadFit.CO2_EOS.cbrt_calc(x)[source]
This function calculates the cubic root that can deal with negative numbers.
- DiadFit.CO2_EOS.density_to_mol_vol(*, density, XH2O)[source]
Converts density in g/cm3 to molar volume (mol/cm3) for a given XH2O
- DiadFit.CO2_EOS.get_EOS_params(P, TK)[source]
This function returns the EOS ‘constants’ if you know the pressure (in bars) and temperature (in Kelvin)
- DiadFit.CO2_EOS.ind_density_homog_T_h_CO2_SW96_loaded_phase_boundary_1sam(T_h_C, print=False, homog_to=None)[source]
Calculates CO2 density for a specified homogenization temperature in Celcius using the Span and Wanger (1996) equation of state. Only works for one sample, use loop_density_homog_T if you have more than one temperature
- Parameters:
- Returns:
colums for Liq_gcm3 and gas_gcm3
- Return type:
pd.DataFrame
- DiadFit.CO2_EOS.mixEOS(V, P, BVc, CVc2, DVc4, EVc5, FVc2, bmix, gVc2, TK)[source]
This function is like the one for the pureEOS. It calculates the compressability factor, and then calculates the compressability factor based on the P-V-T you entered. It returns the residual between those two values.
- DiadFit.CO2_EOS.mixEOS_CF(V, P, BVc, CVc2, DVc4, EVc5, FVc2, bmix, gVc2, TK)[source]
This function is like the one for the pureEOS. It calculates the compressability factor, and then calculates that based on the P-V-T you entered. It does not return the residual
- DiadFit.CO2_EOS.mix_fugacity(*, P_kbar, T_K, XH2O, Vmix)[source]
Used to calculate fugacity, compressability and activities for a panda series
- DiadFit.CO2_EOS.mix_fugacity_ind(*, P_kbar, T_K, XH2O, Vmix)[source]
This function calculates fugacity for a single sample. It returns the activity of each component (fugacity/fugacity in pure component)
- DiadFit.CO2_EOS.mix_lnphi(i, Zmix, BVc_prm, CVc2_prm, DVc4_prm, EVc5_prm, FVc2_prm, FVc2, bmix, b_prm, gVc2, gVc2_prm, Vmix)[source]
This function calculates lnphi values
- DiadFit.CO2_EOS.mixing_rules(B, C, D, E, F, Vc, Y, b, g, k1_temperature, k2_temperature, k3_temperature)[source]
This function applies the DZ06 mixing rules
- DiadFit.CO2_EOS.mixpressure(P, V, TK, Y)[source]
This function iterates in pressure space to get the best match in bars to the entered volume in cm3/mol using the mixEOS function above.
- DiadFit.CO2_EOS.mixvolume(V, P, BVc, CVc2, DVc4, EVc5, FVc2, bmix, gVc2, TK)[source]
This function iterates in volume space to get the best match volume (cm3/mol) to the entered pressure using the mixEOS function above.
- DiadFit.CO2_EOS.mol_vol_to_density(*, mol_vol, XH2O)[source]
Converts molar mass g/mol to densit g/cm3 for a given XH2O
- DiadFit.CO2_EOS.objective_function_CO2_dens(CO2_dens_gcm3, T_K, target_pressure_MPa)[source]
This function minimises the offset between the calculated and target pressure. It finds the temp and CO2 density that correspond to the target pressure
- DiadFit.CO2_EOS.objective_function_Temp(T_K, CO2_dens_gcm3, target_pressure_MPa)[source]
This function minimises the offset between the calculated and target pressure
- DiadFit.CO2_EOS.pureEOS(i, V, P, B, C, D, E, F, Vc, TK, b, g)[source]
This function calculates the compressability factor for a pure EOS using the modified Lee-Kesler equation.
i=0 for H2O, i=1 for CO2.
You input a volume, and it returns the difference between the compresability factor, and that calculated at the input P, V and T_K. E.g. gives the residual so that you can iterate.
- DiadFit.CO2_EOS.pureEOS_CF(i, V, P, B, C, D, E, F, Vc, TK, b, g)[source]
This function calculates the compressability factor for a pure EOS using the modified Lee-Kesler equation.
i=0 for H2O, i=1 for CO2.
You input a volume, and it returns the difference between the compresability factor, and that calculated at the input P, V and T_K. E.g. gives the residual so that you can iterate.
- DiadFit.CO2_EOS.pure_lnphi(i, Z, B, Vc, V, C, D, E, F, g, b)[source]
This function calculates the fugacity coefficient (kbar) from the equation of state for a pure fluid
Functions for calculating CO2 in vapour bubbles
- DiadFit.CO2_in_bubble_error.add_noise_to_variable(original_value, error, error_type, error_dist, N_dup, neg_values, neg_threshold)[source]
This function adds noise to each variable for the monte-carloing
- Parameters:
original_value (int, float) – Preferred value (e.g. center of distribution)
error_type (str) – ‘Abs’ if absolute error, ‘Perc’ if percent
error_dist (str) – ‘normal’ if normally distributed, ‘uniform’ if uniformly distributed.
N_dup (int) – number of duplicates
neg_values (bool) – whether negative values are replaced with zeros
- DiadFit.CO2_in_bubble_error.get_value(variable, i)[source]
This function returns the value if its not a series, otherwise returns the right row in the series
- DiadFit.CO2_in_bubble_error.propagate_CO2_in_bubble(*, N_dup=1000, sample_ID, vol_perc_bub=None, error_vol_perc_bub=None, error_type_vol_perc_bub='Abs', MI_x=None, MI_y=None, MI_z=None, VB_x=None, VB_y=None, VB_z=None, error_MI_x=None, error_MI_y=None, error_MI_z=None, error_VB_x=None, error_VB_y=None, error_VB_z=None, error_type_dimension='Abs', error_dist_dimension='normal', melt_dens_kgm3, CO2_bub_dens_gcm3, error_dist_vol_perc_bub='normal', error_CO2_bub_dens_gcm3=0, error_type_CO2_bub_dens_gcm3='Abs', error_dist_CO2_bub_dens_gcm3='normal', error_melt_dens_kgm3=0, error_type_melt_dens_kgm3='Abs', error_dist_melt_dens_kgm3='normal', plot_figure=True, fig_i=0, neg_values=True)[source]
This function propagates uncertainty in reconstruction of melt inclusion CO2 contents by feeding each row into propagate_CO2_in_bubble_ind. The returned standard deviation uses the 84th-16th percentile rather than the true standard deviation, as this is better for skewed distributions.
- Parameters:
SampleID (str, pd.series) – Sample_ID (e.g. sample names) which is returned on dataframe
N_dup (int) – Number of duplicates when generating errors for Monte Carlo simulations
Now
error (either you know the vol% of the bubble and the associated)
arguements (in which case enter these two)
vol_perc_bub (int, float, pd.series) – Volume proportion of sulfide in melt inclusion
error_vol_perc_bub (int, float, pd.Series) – Error for each variable, can be absolute or %
CO2_bub_dens_gcm3 (int, float, pd.Series) – Error for each variable, can be absolute or %
error_melt_dens_kgm3 (int, float, pd.Series) – Error for each variable, can be absolute or %
error_type_vol_perc_bub ('Abs' or 'Perc') – whether given error is perc or absolute
error_type_bub_dens_gcm3 ('Abs' or 'Perc') – whether given error is perc or absolute
error_type_melt_dens_kgm3 ('Abs' or 'Perc') – whether given error is perc or absolute
Or
polishing). (you have measured the x-y-z dimensions of the MI and VB (e.g. by perpendicular)
MI_x (int, float, pd.Series) – x, y, z dimensions of vapour bubble and melt inclusion
MI_y (int, float, pd.Series) – x, y, z dimensions of vapour bubble and melt inclusion
MI_Z (int, float, pd.Series) – x, y, z dimensions of vapour bubble and melt inclusion
VB_x (int, float, pd.Series) – x, y, z dimensions of vapour bubble and melt inclusion
VB_y (int, float, pd.Series) – x, y, z dimensions of vapour bubble and melt inclusion
VB_Z (int, float, pd.Series) – x, y, z dimensions of vapour bubble and melt inclusion
error_MI_x (int, float, pd.Series) – Error on x, y, z dimension meausurement
error_MI_y (int, float, pd.Series) – Error on x, y, z dimension meausurement
error_MI_Z (int, float, pd.Series) – Error on x, y, z dimension meausurement
error_VB_x (int, float, pd.Series) – Error on x, y, z dimension meausurement
error_VB_y (int, float, pd.Series) – Error on x, y, z dimension meausurement
error_VB_Z (int, float, pd.Series) – Error on x, y, z dimension meausurement
error_type_dimension ('Abs' or 'Perc') – Whether errors on all x-y-z are perc or abs
error_dist_dimension ('normal' or 'uniform') – Whether errors on all x-y-z are normally or uniformly distributed.
- melt_dens_kgm3:int, float, pd.series
Density of the silicate melt in kg/m3, e.g. from DensityX
- CO2_bub_dens_gcm3: int, float, pd.Series
Density of the vapour bubble in g/cm3
- error_dist_vol_perc_bub, error_dist_bub_dens_gcm3, error_dist_melt_dens_kgm3, error_dist_dimension: ‘normal’ or ‘uniform’
Distribution of simulated error
- plot_figure: bool
Default true - plots a figure of the row indicated by fig_i (default 1st row, fig_i=0)
- neg_values: bool
Default True - whether negative values are removed from MC simulations or not. False, replace all negative values with zeros.
- Returns:
pd.DataFrame – All outputs has calculations for every simulaiton df_step has average and standard deviation for each sample
- Return type:
df_step, All_outputs
- DiadFit.CO2_in_bubble_error.propagate_CO2_in_bubble_ind(sample_i=0, N_dup=1000, vol_perc_bub=None, CO2_bub_dens_gcm3=None, melt_dens_kgm3=None, MI_x=None, MI_y=None, MI_z=None, VB_x=None, VB_y=None, VB_z=None, error_MI_x=None, error_MI_y=None, error_MI_z=None, error_VB_x=None, error_VB_y=None, error_VB_z=None, error_type_dimension='Abs', error_dist_dimension='normal', error_vol_perc_bub=None, error_type_vol_perc_bub='Abs', error_dist_vol_perc_bub='normal', error_CO2_bub_dens_gcm3=0, error_type_CO2_bub_dens_gcm3='Abs', error_dist_CO2_bub_dens_gcm3='normal', error_melt_dens_kgm3=0, error_type_melt_dens_kgm3='Abs', error_dist_melt_dens_kgm3='normal', len_loop=1, neg_values=True)[source]
This function propagates uncertainty in reconstruction of melt inclusion bubble equivalent CO2 contents. and returns a dataframe. it does one sample at a time.
- Parameters:
- error_vol_perc_bub, CO2_bub_dens_gcm3, error_melt_dens_kgm3: int, float, pd.Series
Error for each variable, can be absolute or %
- error_type_vol_perc_bub, error_type_bub_dens_gcm3, error_type_melt_dens_kgm3: ‘Abs’ or ‘Perc’
whether given error is perc or absolute
- error_dist_vol_perc_bub, error_dist_bub_dens_gcm3, error_dist_melt_dens_kgm3: ‘normal’ or ‘uniform’
Distribution of simulated error
OR Enter melt inclusion and vapour bubble dimensions (diameter, not radii), and their errors MI_x, MI_y, MI_z, VB_x, VB_y, VB_z: int, float, series:
Diameter of melt inclusions.
- error_MI_x, error_MI_y, error_MI_z, error_VB_x, error_VB_y, error_VB_z:
Error on diameter of melt inclusions
- error_type_dimension=’Abs’ or ‘Perc’:
Specify whether errors on these dimensions are absolute or percentage
- error_dist_dimension=’normal’ or ‘uniform’:
Specify error distribution
- SampleID: str, pd.series
Sample_ID (e.g. sample names) which is returned on dataframe
- N_dup: int
Number of duplicates when generating errors for Monte Carlo simulations
- melt_dens_kgm3:int, float, pd.series
Density of the silicate melt in kg/m3, e.g. from DensityX
- CO2_bub_dens_gcm3: int, float, pd.Series
Density of the vapour bubble in g/cm3
- plot_figure: bool
Default true - plots a figure of the row indicated by fig_i (default 1st row, fig_i=0)
- neg_values: bool
Default True - whether negative values are removed from MC simulations or not. False, replace all negative values with zeros.
- Returns:
Input variable duplicated N_dup times with noise added.
- Return type:
pd.DataFrame
Functions for calculating molar gas proportions from Raman peak areas
- DiadFit.molar_gas_proportions.calculate_SO2_CO2_mol_prop_wave_indep(SO2_wavelength_ind, CO2_diad1_wavelength_ind, CO2_diad2_wavelength_ind, wavelength_nm, T_C, A_SO2, A_CO2_Tot)[source]
Takes wavelength independnet cross sections and CO2 and SO2 peak areas and converts them into SO2 mol proportions :param SO2_wavelength_ind: :type SO2_wavelength_ind: Wavelength independent cross section for SO2 :param CO2_diad1_wavelength_ind: :type CO2_diad1_wavelength_ind: Wavelength independent cross section for diad 1 (at 1285) :param CO2_diad2_wavelength_ind: :type CO2_diad2_wavelength_ind: Wavelength independent cross section for diad 2 (at 1388) :param wavelength_nm: :type wavelength_nm: Laser wavelenth of system in nm :param T_C: :type T_C: Temperature of analysis in C.
- Return type:
SO2 mol proportion
- DiadFit.molar_gas_proportions.calculate_SO2_CO2_ratio(SO2_area, diad1_area, diad2_area, SO2_cross_sec=5.3, diad1_cross_sec=0.89, diad2_cross_sec=1.4)[source]
Calculates SO2:CO2 ratio using the parameters from Marie-Camille Caumons lab
- DiadFit.molar_gas_proportions.calculate_wavelength_dependent_cross_section(wavelength_nm, T_C, Raman_shift_cm, wavelength_independent_cross_section)[source]
This function calculates the wavelength dependent cross section (lower case sigma) from the wavelength independent Raman scattering efficiency (Upper case sigma)
- Parameters
- wavelength_nm:
laser wavelength used in nm to calculate the wavelength dependent cross section
- wavelength_independent_cross_section:
Wavelength independent cross section
- T_K:
absolute temperature in Kelvin
- Raman_shift_cm:
Raman shift in cm-1 of peak of interest (e.g. 1151 for SO$_2$)
Wavelength dependent cross section
- DiadFit.molar_gas_proportions.calculate_wavelength_independent_cross_section(wavelength_nm, T_C, Raman_shift_cm, wavelength_dependent_cross_section)[source]
This function calculates the wavelength independent cross section (capital Sigma) from the wavelength dependent Raman scattering efficiency (lower case sigma)
- Parameters
- wavelength_nm:
laser wavelength used in nm to calculate the wavelength dependent cross section
- wavelength_dependent_cross_section:
Wavelength dependent cross section
- T_K:
absolute temperature in Kelvin
- Raman_shift_cm:
Raman shift in cm-1 of peak of interest (e.g. 1151 for SO$_2$)
Wavelength independent cross section
- DiadFit.molar_gas_proportions.convert_cross_section_wavelength1_wavelength2(wavelength_nm_1, wavelength_nm_2, Raman_shift_cm, wavelength_dependent_cross_section_wavelength1, T_C)[source]
This function calculates the wavelength dependent cross section (lower case sigma) for laser wavelength 2 from the wavelength-dependent cross section for laser wavelength 1.
- Parameters
- wavelength_nm_1:
laser wavelength used in nm to calculate the wavelength dependent cross section
- wavelength_nm_2:
laser wavelength of the system of interest you are trying to calculate the wavelength dependent cross section for
- wavelength_dependent_cross_section_wavelength1:
Wavelength dependent cross section
- T_K:
absolute temperature in Kelvin
- Raman_shift_cm:
Raman shift in cm-1 of peak of interest (e.g. 1151 for SO$_2$)
Wavelength dependent cross section for wavelength2
Functions for converting pressure to depth using different crustal profiles
- DiadFit.density_depth_crustal_profiles.convert_co2_dens_press_depth_old(T_K=None, CO2_dens_gcm3=None, crust_dens_kgm3=None, output='kbar', g=9.81, model=None, d1=None, d2=None, rho1=None, rho2=None, rho3=None, EOS='SW96')[source]
This is a now old function that isn’t used, kept for backwards functionality. Dont use unless you have built code relying on it!
- DiadFit.density_depth_crustal_profiles.convert_pressure_depth_2step(P_kbar=None, d1=None, rho1=None, rho2=None, g=9.81)[source]
Converts Pressure to depth using a 2 step profile for int or float
- DiadFit.density_depth_crustal_profiles.convert_pressure_depth_3step(P_kbar=None, d1=5, d2=14, rho1=2700, rho2=3000, rho3=3100, g=9.81)[source]
Converts Pressure to depth using a 3 step profile for int or float
- Parameters:
- rho3: int or float
Density (kg/m3) below 2nd step transition
- Returns:
Depth in km
- Return type:
- DiadFit.density_depth_crustal_profiles.convert_pressure_depth_4step(*, P_kbar=None, d1=None, d2=None, d3=None, rho1=None, rho2=None, rho3=None, rho4=None, g=9.81)[source]
Converts Pressure to depth using a 3 step profile for int or float
- Parameters:
- rho1: int or float
Density (kg/m3) down to first step transition
- rho2: int or float
Density (kg/m3) between first and second step transition
- rho3: int or float
Density (kg/m3) below 2nd step transition
- rho4: int or float
Density (kg/m3) below 3rd step transition
- Returns:
Depth in km
- Return type:
- DiadFit.density_depth_crustal_profiles.convert_pressure_to_depth(*, P_kbar=None, crust_dens_kgm3=None, g=9.81, d1=None, d2=None, d3=None, rho1=None, rho2=None, rho3=None, rho4=None, model=None)[source]
Converts pressure in kbar to depth in km using a variety of crustal density profiles
- Parameters:
P_kbar (int, float, pd.Series, np.ndarray) – Pressure in kbar
g (float) – gravitational constant, in m/s2
crust_dens_kgm3 –
If float: Crustal density in kg/m3
If model, choose from:
- ryan_lerner:
Parameterization of Ryan 1987, actual equation from Lerner et al. 2021 After 16.88 km (455 MPa), assume density is 2.746, as density turns around again. This profile is tweaked for Hawaii
- denlinger_lerner:
Parameterization of Denlinger and Flinders (2022) with the addition of 200 kg/m3 at shallow depths, from Lerner et al. (2024). After 15 km returns NaN following Lerner
- mavko_debari:
Parameterization of Mavko and Thompson (1983) and DeBari and Greene (2011) as given in Putirka (2017) Down the Crater Elements supplement.
OR
Else, just enter a crustal density in kg/m3, e.g., model=2700
- Return type:
Depth in km as a panda series
- DiadFit.density_depth_crustal_profiles.denlinger_lerner(P_kbar)[source]
Calculates depth for a given pressure using the Parameterization of Denlinger and Flinders (2022) with the addition of 200 kg/m3 at shallow depths, from Lerner et al. (2024). After 15 km returns NaN following Lerner
- DiadFit.density_depth_crustal_profiles.hill_zucca(P_kbar)[source]
Parameterization of Hill and Zucca (1987), as given in Putirka (2017) Down the Crater Elements supplement
- DiadFit.density_depth_crustal_profiles.loop_pressure_depth_2step(P_kbar=None, d1=14, rho1=2800, rho2=3100, g=9.81)[source]
Converts Pressure to depth using a 2 step profile for a pandas.Series of presssures
- DiadFit.density_depth_crustal_profiles.loop_pressure_depth_3step(P_kbar=None, d1=5, d2=14, rho1=2700, rho2=3000, rho3=3100, g=9.81)[source]
Converts Pressure to depth using a 3 step profile for pd.Series
- Parameters:
P_kbar (pd.Series) – Pressure in kbar
d1 (int or float) – depth in km of 1st step transition in density
d2 (int or float) – depth in km of 2nd step transition in density
rho1 (int or float) – Density (kg/m3) down to first step transition
rho2 (int or float) – Density (kg/m3) between first and second step transition
- rho3: int or float
Density (kg/m3) below 2nd step transition
- Returns:
Depth in km
- Return type:
pd.Series
- DiadFit.density_depth_crustal_profiles.loop_pressure_depth_4step(*, P_kbar=None, d1=5, d2=14, d3=20, rho1=2700, rho2=3000, rho3=3100, rho4=4000, g=9.81)[source]
Converts Pressure to depth using a 3 step profile for pd.Series
- Parameters:
P_kbar (pd.Series) – Pressure in kbar
d1 (int or float) – depth in km of 1st step transition in density
d2 (int or float) – depth in km of 2nd step transition in density
rho1 (int or float) – Density (kg/m3) down to first step transition
rho2 (int or float) – Density (kg/m3) between first and second step transition
- rho3: int or float
Density (kg/m3) below 2nd step transition
- Returns:
Depth in km
- Return type:
pd.Series
- DiadFit.density_depth_crustal_profiles.mavko_debari(P_kbar)[source]
Parameterization of Mavko and Thompson (1983) and DeBari and Greene (2011) as given in Putirka (2017) Down the Crater Elements supplement, used for Cascades
- DiadFit.density_depth_crustal_profiles.prezzi(P_kbar)[source]
Parameterization of Prezzi et al. (2009), as given in Putirka (2017) Down the Crater Elements supplement. Used for Andes.
Functions for modelling fluid inclusion re-equilibration
- DiadFit.relaxifi.calculate_DPdt(ascent_rate_ms, crustal_model_config=<DiadFit.relaxifi.config_crustalmodel object>, D_initial_km=None, D_final_km=None, D_step=100, initial_P_guess_kbar=0, tolerance=0.001)[source]
Calculate the decompression rate (dP/dt) during ascent.
Parameters: - ascent_rate_ms (float): Ascent rate in meters per second. - D_initial_km (float, optional): Initial depth in kilometers. Default is 30 km. - D_final_km (float, optional): Final depth in kilometers. Default is 0 km. - D_step (int, optional): Number of depth steps for calculation. Default is 100. - initial_P_guess_kbar (float, optional): Initial guess for pressure in kilobars (kbar). Default is 0. - tolerance (float, optional): Tolerance for pressure estimation. Default is 0.001.
Returns: - D (pd.Series): Depth values in kilometers. - Pexternal_steps_MPa (list): Lithostatic pressure values in megapascals (MPa) at each depth step. - dt (float): Time step for the integration.
- DiadFit.relaxifi.calculate_dR_dt(*, R_m, b_m, T_K, Pinternal_MPa, Pexternal_MPa)[source]
Calculate the rate of change of inclusion radius (dR/dt) based on power law creep.
Parameters: - R_m (float): Inclusion radius in meters. - b_m (float): Distance to the crystal defect structures. Wanamaker and Evans (1989) use R/b=1/1000. - T_K (float): Temperature in Kelvin. - Pinternal_MPa (float): Internal pressure in MPa. - Pexternal_MPa (float): External pressure in MPa.
Returns: - dR_dt (float): Rate of change of inclusion radius in meters per second.
- DiadFit.relaxifi.calculate_exponential_time_steps(totaltime_s, steps)[source]
Generate exponentially spaced time steps with smaller steps at the beginning and larger steps towards the end.
Parameters: - totaltime_s (float): The total time for the simulation. - steps (int): The number of steps in the simulation.
Returns: - dt_s_array (np.ndarray): Array of time intervals (dt_s) with smaller steps at the beginning.
- class DiadFit.relaxifi.config_crustalmodel(crust_dens_kgm3=None, d1=None, d2=None, rho1=None, rho2=None, rho3=None, model=None)[source]
A configuration class for specifying parameters of the crustal model.
Attributes: - crust_dens_kgm3 (float): The density of the crust in kilograms per cubic meter (kg/m^3). - d1 (float): The depth boundary for the first layer in kilometers (km). - d2 (float): The depth boundary for the second layer in kilometers (km). - rho1 (float): The density of the first layer in kilograms per cubic meter (kg/m^3). - rho2 (float): The density of the second layer in kilograms per cubic meter (kg/m^3). - rho3 (float): The density of the third layer in kilograms per cubic meter (kg/m^3). - model (str): The name of the model used for crustal calculations.
- DiadFit.relaxifi.find_P_for_kmdepth(target_depth_km, crustal_model_config=<DiadFit.relaxifi.config_crustalmodel object>, initial_P_guess_kbar=0, tolerance=0.1)[source]
Approximate the pressure (P_kbar) based on the target depth using the Newton-Raphson method.
Parameters: - target_depth_km (float, pd.Series, list): The desired depth(s) in kilometers (km). - initial_P_guess_kbar (float, optional): Initial guess for the pressure in kilobars (kbar). Default is 0. - crustal_model_config (object, optional): Configuration object containing crustal model parameters.
crust_dens_kgm3 (float, optional): The density of the crust in kilograms per cubic meter (kg/m^3). Default is None.
d1, d2 (float, optional): Depth boundaries for different layers (km). Default is None.
rho1, rho2, rho3 (float, optional): Densities for different layers (kg/m^3). Default is None.
model (str, optional): The name of the model used for depth calculation. Default is None.
tolerance (float, optional): Tolerance for the Newton-Raphson method. The pressure estimate should be within this tolerance of the true value. Default is 0.1.
Returns: - float or pd.Series or list: The estimated pressure(s) (P_kbar) that correspond to the target depth(s).
Notes: - If the target_depth_km is a single value, a float is returned. - If the target_depth_km is a Pandas Series, a Pandas Series is returned. - If the target_depth_km is a list, a list of floats is returned.
If crustal parameters are not provided in the crustal_model_config object, a single step model with a crustal density = 2750 kg/cm3 will be used.
- DiadFit.relaxifi.get_CO2dens_P(R_m, T_K, CO2_mass, EOS='SW96', return_volume=False)[source]
Calculate the density and pressure of CO2 inside a fluid inclusion (FI).
Parameters: - R_m (float): The radius of the fluid inclusion (FI), in meters. - T_K (float): The temperature, in Kelvin. - CO2_mass (float): The mass of CO2 within the FI, in grams (g). - EOS (str, optional): The equation of state (EOS) to use for density and pressure calculations.
Can be one of: ‘ideal’ (ideal gas), ‘SW96’ (Span and Wagner 1996), or ‘SP94’ (Sterner and Pitzer 1994). Defaults to ‘SW96’.
return_volume (bool, optional): Whether to return the volume of the FI along with density and pressure. Defaults to False.
Returns: - tuple or float: If return_volume is True, returns a tuple containing (V, CO2_dens, P), where:
V (float): The volume of the fluid inclusion (FI), in cubic meters (m³).
CO2_dens (float): The density of CO2 within the FI, in grams per cubic centimeter (g/cm³).
P (float): The pressure of CO2 within the FI, in MegaPascals (MPa).
If return_volume is False, returns a tuple containing (CO2_dens, P).
- DiadFit.relaxifi.get_R(R_m, b_m, T_K, Pinternal_MPa, Pexternal_MPa, dt_s, method='RK1')[source]
Find the radius R of an FI over a time step using the Runge-Kutta numerical method. Options are order 1 to 4 RK methods, such as RK1 (Euler), RK2 (Heun), RK3, or RK4.
Parameters: - R_m (float): Initial FI Radius in meters. - b_m (float): Distance to defect structures in meters. - T_K (float): Temperature in Kelvin. - Pinternal_MPa (float): Internal pressure in MPa. - Pexternal_MPa (float): External pressure in MPa. - dt_s (float): The time step for integration in seconds. - method (str, optional): The numerical integration method to use. Default is ‘RK1’.
Returns: - tuple: A tuple containing the updated value of R_m and the derivative dR_dt (rate of change of R).
- DiadFit.relaxifi.get_initial_CO2(R_m, T_K, P_MPa, EOS='SW96', return_volume=False)[source]
Calculate the initial density and mass of CO2 inside a fluid inclusion (FI).
Parameters: - R_m (float): The radius of the fluid inclusion (FI), in meters. - T_K (float): The temperature, in Kelvin. - P_MPa (float): The pressure, in MegaPascals (MPa). - EOS (str, optional): The equation of state (EOS) to use for density calculations.
Can be one of: ‘ideal’ (ideal gas), ‘SW96’ (Span and Wagner EOS 1996), or ‘SP94’ (Sterner and Pitzer EOS 1994). Defaults to ‘SW96’.
return_volume (bool, optional): Whether to return the volume of the FI along with density and mass. Defaults to False.
Returns: - tuple or float: If return_volume is True, returns a tuple containing (V, CO2_dens_initial, CO2_mass_initial), where:
V (float): The volume of the fluid inclusion (FI), in cubic meters (m³).
CO2_dens_initial (float): The initial density of CO2 within the FI, in grams per cubic centimeter (g/cm³).
CO2_mass_initial (float): The initial mass of CO2 within the FI, in grams (g).
If return_volume is False, returns a tuple containing (CO2_dens_initial, CO2_mass_initial).
- DiadFit.relaxifi.loop_R_b_constant_Pext(*, R_m_values, b_m_values, T_K, EOS, Pinternal_MPa, Pexternal_MPa, totaltime_s, steps, T4endcalc_PD, method='RK4', plotfig=False, crustal_model_config=<DiadFit.relaxifi.config_crustalmodel object>, step_type='linear', max_perc_change=5)[source]
Perform multiple simulations under constant external pressure with various R and b values.
Parameters: - R_m_values (list): A list of initial radius values for fluid inclusions (FI), in meters. - b_m_values (list): A list of initial distance values to a crystal defect (rim, crack, etc) from the FI center, in meters. - T_K (float): The temperature, in Kelvin. - EOS (str): The equation of state (EOS) to use for density calculations. Can be one of: ‘ideal’ (ideal gas),
‘SW96’ (Span and Wagner 1996), or ‘SP94’ (Sterner and Pitzer 1994).
Pinternal_MPa (float): The initial internal pressure of the FI, in MegaPascals (MPa).
Pexternal_MPa (float): The constant external pressure applied to the FI, in MegaPascals (MPa).
totaltime_s (float): The total simulation time, in seconds.
steps (int): The number of simulation steps.
T4endcalc_PD (float): The temperature at which to calculate the depths (Kelvin) at the end of the simulations.
method (str, optional): The numerical integration method to use for change in FI radius. Can be ‘RK1’ (also ‘Euler’), ‘RK2’ (also ‘Heun’), ‘RK3’, or ‘RK4’. Defaults to ‘RK4’.
plotfig (bool, optional): Whether to plot figures showing the changes in time, radius, and CO2 density. Defaults to False.
- crustal_model_config (dict, optional): Configuration parameters for the crustal model. Defaults to a predefined
configuration with a crustal density of 2750 kg/m³.
Returns: - dict: A dictionary containing simulation results for varying R and b values. The keys are in the format ‘R{index}_b{index}’,
where ‘index’ corresponds to the index of R_values and b_values, respectively. The values are DataFrames with simulation results.
- DiadFit.relaxifi.objective_function_depth(P_kbar, target_depth_km, crust_dens_kgm3, d1, d2, rho1, rho2, rho3, model)[source]
Calculate the difference between the current depth and the target depth given pressure (P_kbar) and other parameters.
Parameters: - P_kbar (float): The pressure in kilobars (kbar) to be used in the depth calculation. - target_depth_km (float): The desired depth in kilometers (km). - crust_dens_kgm3 (float): The density of the crust in kilograms per cubic meter (kg/m^3). - d1, d2 (float): Depth boundaries for different layers (km). - rho1, rho2, rho3 (float): Densities for different layers (kg/m^3). - model (str): The name of the model used for the depth calculation.
Returns: - float: The difference between the current depth and the target depth.
- class DiadFit.relaxifi.power_creep_law_constants[source]
Olivine power-law creep constants used in the stretching model (Wanamaker and Evans, 1989).
Attributes: - A (float): Creep law constant A (default: 3.9e3). - n (float): Creep law constant n (default: 3.6). - Q (float): Activation energy for dislocation motions in J/mol (default: 523000). - IgasR (float): Gas constant in J/(mol*K) (default: 8.314).
- DiadFit.relaxifi.stretch_at_constant_Pext(*, R_m, b_m, T_K, EOS='SW96', Pinternal_MPa, Pexternal_MPa, totaltime_s, steps, method='RK4', report_results='fullpath', plotfig=False, update_b=False, step_type='linear', max_perc_change=5)[source]
Simulate the stretching of a CO2 fluid inclusion (FI) under constant external pressure (e.g., quenching or storage).
Parameters: - R_m (float): The initial radius of the fluid inclusion (FI), in meters. - b_m (float): The initial distance to a crystal defect (rim, crack, etc) from the FI center, in meters. - T_K (float): The temperature, in Kelvin. - Pinternal_MPa (float): The initial internal pressure of the FI, in MegaPascals (MPa). - Pexternal_MPa (float): The constant external pressure applied to the FI, in MegaPascals (MPa). - totaltime_s (float): The total simulation time, in seconds. - steps (int): The number of simulation steps. - EOS (str): The equation of state (EOS) to use for density calculations. Can be one of: ‘ideal’ (ideal gas),
‘SW96’ (Span and Wagner 1996), or ‘SP94’ (Sterner and Pitzer 1994).
method (str, optional): The numerical integration method to use for change in FI radius. Can be ‘RK1’ (also ‘Euler’), ‘RK2’ (also ‘Heun’), ‘RK3’, or ‘RK4’. Defaults to ‘RK4’.
report_results (str, optional): The type of results to report. Can be ‘fullpath’ (all steps reported), ‘startendonly’ (only initial and end steps), or ‘endonly’ (only last step). Defaults to ‘fullpath’.
plotfig (bool, optional): Whether to plot figures showing the changes in time, radius, and CO2 density. Defaults to False.
update_b (bool, optional): Whether to update ‘b’ during the simulation. Defaults to False.
Returns: - pandas.DataFrame: A DataFrame containing the simulation results, including time, pressure, radius changes, and CO2 density.
Raises: - ValueError: If an unsupported EOS is specified.
- DiadFit.relaxifi.stretch_in_ascent(*, R_m, b_m, T_K, ascent_rate_ms, depth_path_ini_fin_step=[100, 0, 100], crustal_model_config=<DiadFit.relaxifi.config_crustalmodel object>, EOS, plotfig=True, report_results='fullpath', initial_P_guess_kbar=0, tolerance=0.001, method='RK4', update_b=False)[source]
Simulate the stretching of a CO2-dominated fluid inclusion (FI) during ascent.
Parameters: - R_m (float): The initial radius of the fluid inclusion (FI), in meters. - b_m (float): The initial distance to a crystal defect (rim, crack, etc) from the FI center, in meters. - T_K (float): The temperature, in Kelvin. - ascent_rate_ms (float): The ascent rate, in meters per second (m/s). - depth_path_ini_fin_step (list, optional): A list containing [initial_depth_km, final_depth_km, depth_step].
Defaults to [100, 0, 100], representing the depth path from initial to final depth in a number of steps.
crustal_model_config (dict, optional): Configuration parameters for the crustal model. Defaults to a predefined configuration with a crustal density of 2750 kg/m³.
EOS (str): The equation of state (EOS) to use for density calculations. Can be one of: ‘ideal’ (ideal gas), ‘SW96’ (Span and Wagner 1996), or ‘SP94’ (Sterner and Pitzer 1994).
plotfig (bool, optional): Whether to plot figures showing the changes in depth and CO2 density. Defaults to True.
report_results (str, optional): The type of results to report. Can be ‘fullpath’, ‘startendonly’, or ‘endonly’. Defaults to ‘fullpath’.
initial_P_guess_kbar (float, optional): Initial guess for internal pressure (Pinternal_MPa) in MPa. Defaults to 0.
tolerance (float, optional): Tolerance for pressure calculations. Defaults to 0.001.
method (str, optional): The numerical integration method to use for change in FI radius. Can be ‘RK1’ (also ‘Euler’), ‘RK2’ (also ‘Heun’), ‘RK3’, or ‘RK4’. Defaults to ‘RK4’.
update_b (bool, optional): Whether to update ‘b’ during the ascent. Defaults to False.
Returns: - pandas.DataFrame: A DataFrame containing the simulation results, including time, depth, pressure, radius changes,
and CO2 density.
Functions for fitting H2O peaks in Raman spectra
- DiadFit.H2O_fitting.check_if_spectra_negative(*, path, filename, Spectra=None, peak_pos_Host=None, tie_x_cord=2000, override=False, flip=False, plot_figure=True, dpi=200)[source]
This function checks if the unmixed specta is negative, based on two tie points. The first tie point is the mean y coordinate of the peak position of host +5 wavenumbers, and the second tie point (tie_x_cord) is an optional input. If the specta is inverted, this function inverts it.
- Parameters:
Spectra (np.array) – Spectra from the function make_evaluate_mixed_spectra
peak_pos_Host (list) – Host peak positions from the function find_host_peak_trough_pos
tie_x_cord (int or float) – X cooordinate to use as a tie point to ask whether the host peak’s y coordinate is higher or lower than this.
override (bool) – if False, function flips the spectra if its upsideown, if True, you can use the input ‘flip’ to manually flip the spectra
flip (bool) – If override is true, flip=False leaves spectra how it is, True flips the y axis.
- Returns:
Spectra – Flipped or unflipped spectra
- Return type:
np.array
- DiadFit.H2O_fitting.extract_xstal_MI_name(*, files, char_xstal, pos_xstal, char_MI, pos_MI, prefix=True, str_prefix=' ', file_ext='.txt')[source]
Extracts the names of the crystal and MI samples from a list of filenames
- Parameters:
(list) (files)
(str) (char_MI) – The character or string used to split the filenames into parts. E.g. ‘_’ if the filename is of the form ‘FM_7_MI’.
(str) – The character or string used to split the filenames into parts. E.g. ‘_’ if the filename is of the form ‘FM_7_MI’.
(int) (pos_MI)
(int)
(bool (prefix) – If True, removes prefix. E.g. WITEC instruments where 01 is appended onto the first file.
optional) (The file extension of the filenames. Default is ".txt".) – If True, removes prefix. E.g. WITEC instruments where 01 is appended onto the first file.
(str (file_type)
optional)
(str
optional)
Returns – df_out (pandas DataFrame): A dataframe with columns “filename”, “crystal_name”, and “MI_name”, containing the input filenames, the extracted crystal sample names, and the extracted MI sample names, respectively.
- DiadFit.H2O_fitting.find_host_peak_trough_pos(smoothed_host_y, x_new, height=1)[source]
“ This function identifies the peaks and troughs in the host mineral spectra
- Parameters:
path (smoothed_host_y) – Host spectra y values after applying a cubic spline, and trimming to the spectra region around the peaks (from function smooth_and_trim_around_host)
x_new (X values corresponding to y values in smoothed_ol_y)
height (int) – Height used for scipy find peaks function. May need tweaking
Returns
-----------
(x) (Peak positions)
(y) (peak heights)
position (trough y)
position
- DiadFit.H2O_fitting.fit_area_for_silicate_region(*, path, filename, Spectra=None, config1: sil_bck_pos_Schiavi_basalt(lower_range_sil=300, 340, mid_range1_sil=630, 640, mid_range2_sil=800, 830, upper_range_sil=1200, 1250, LW=400, 600, HW=800, 1200, N_poly_sil=3, sigma_sil=5), sigma_sil=5, exclude_range1_sil=None, exclude_range2_sil=None, plot_figure=True, save_fig=True, fit_sil='poly', dpi=200)[source]
Calculates the area of silicate peaks in a spectrum and fits a polynomial or spline curve to the baseline of the spectrum.
- Parameters:
path (str) – File path
filename (str) – File name
Spectra (numpy.ndarray, optional) – 2D array representing the spectrum data
config1 (object) – Configuration object for silicate peak and background positions. Default values stored in the dataclasses ‘sil_bck_pos_Schiavi_basalt’, _andesite, etc. Can tweak these as well.
exclude_range1_silicate (list of float, optional) – List of two numbers representing the start and end of a range of wavenumbers to be excluded from the spectrum
exclude_range2_silicate (list of float, optional) – List of two numbers representing the start and end of a second range of wavenumbers to be excluded from the spectrum
plot_figure (bool, optional) – Indicates whether or not to plot the fit (default is True)
dpi (int, optional) – Resolution of the plot in dots per inch (default is 200)
- Returns:
DataFrame with columns for ‘Silicate_Trapezoid_Area’, ‘Silicate_Simpson_Area’, as well as parameters for the selected background positions
- Return type:
pd.DataFrame
- DiadFit.H2O_fitting.fit_area_for_water_region(*, path, filename, Spectra=None, config1: water_bck_pos(fit_water='poly', N_poly_water=3, lower_bck_water=2750, 3100, upper_bck_water=3750, 4100), exclude_range1_water=None, exclude_range2_water=None, plot_figure=True, dpi=200)[source]
Calculates the area of water peaks in a spectrum and fits a polynomial or spline curve to the baseline of the spectrum.
- Parameters:
path (str) – File path
filename (str) – File name
Spectra (numpy.ndarray, optional) – 2D array representing the spectrum data
config1 (object) –
Configuration object for water peak and background positions. Default parameters stored in water_bck_pos, user can tweak. Parameters that need tweaking:
fit_water: str ‘poly’, N_poly_water: str, degree of polynomial to fit to background lower_bck_water: [float, float], background position to left of water peak upper_bck_water: [float, float], background position to right of water peak
exclude_range1_water (list of float, optional) – List of two numbers representing the start and end of a range of wavenumbers to be excluded from the spectrum
exclude_range2_water (list of float, optional) – List of two numbers representing the start and end of a second range of wavenumbers to be excluded from the spectrum
plot_figure (bool, optional) – Indicates whether or not to plot the fit (default is True)
dpi (int, optional) – Resolution of the plot in dots per inch (default is 200)
- Returns:
DataFrame with columns for different areas, and parameters for choosen background positions.
- Return type:
pd.DataFrame
- DiadFit.H2O_fitting.make_evaluate_mixed_spectra(*, path, filename, smoothed_host_y, smoothed_MI_y, Host_spectra, MI_spectra, x_new, peak_pos_Host, trough_x, trough_y, N_steps=20, av_width=2, X_min=0, X_max=1, plot_figure=True, dpi=200)[source]
This function unmixes glass and host spectra, and fits the best fit proportion where the host peak and trough disapears. Specifically, it calculates the mixed spectra by taking the measured MI spectra and subtracting X*Ol spectra, where X is the mixing proportions
- Parameters:
smoothed_host_y (np.array) – y coordinates of host around peak region (from the function smooth_and_trim_around_host)
smoothed_MI_y (np.array) – y coordinates of melt inclusion around peak region (from the function smooth_and_trim_around_host)
x_new (np.array) – x coordinates from smoothed Ol and MI curves (from the function smooth_and_trim_around_host)
Host_Spectra (np.array) – Full host spectra, not trimmed or smoothed (from the function get_data)
MI_Spectra (np.array) – Full MI spectra, not trimmed or smoothed (from the function get_data)
peak_pos_Host (list) – Peak positions (x) of Olivine peaks (from the function smooth_and_trim_around_host)
trough_x (float, int) – Peak position (x) of Olivine trough (from the function smooth_and_trim_around_host)
N_steps (int) – Number of mixing steps to use between X_Max and X_Max. E.g. Precisoin of mixed value.
av_width (int) – averages +- 1 width either side of the peak and troughs when doing assesment and regression
Returns:
- MI_Mix_Best: np.array
Spectra of best-fit unmixed spectra (e.g. where host peak and trough the smallest)
- ideal_mix: float
Best fit mixing proportion (i.e. X)
- Dist: float
Vertical distance between the host peak and trough (in intensity units)
- MI_Mix: np.array
Ubmixed spectra for each of the N_steps
- X: np.array
X coordinates of unmixed spectra (along with MI_Mix and X allows plots of unmixing)
if plot_figure is True, also returns a plot showing the unmixing process
- class DiadFit.H2O_fitting.sil_bck_pos_Schiavi_andesite(lower_range_sil: [Tuple[float, float]] = (230, 250), mid_range1_sil: [Tuple[float, float]] = (645, 700), mid_range2_sil: [Tuple[float, float]] = (825, 840), upper_range_sil: [Tuple[float, float]] = (1230, 1250), N_poly_sil: float = 3, sigma_sil: float = 5)[source]
Configuration object for Silicate background positions from Schiavi et al. (2018) for andesites
- Parameters:
here) (spectral range taken as background positions (from left to right as listed)
fit_silicate (str) – Type of fit for the baseline of the silicate region (‘poly’ or ‘spline’)
N_poly_silicate (int) – Degree of polynomial fit for the baseline of the silicate region
sigma_water (int) – Allow points on background within +-sigma_silicate * std dev of other background points
- class DiadFit.H2O_fitting.sil_bck_pos_Schiavi_basalt(lower_range_sil: [Tuple[float, float]] = (300, 340), mid_range1_sil: [Tuple[float, float]] = (630, 640), mid_range2_sil: [Tuple[float, float]] = (800, 830), upper_range_sil: [Tuple[float, float]] = (1200, 1250), LW: [Tuple[float, float]] = (400, 600), HW: [Tuple[float, float]] = (800, 1200), N_poly_sil: float = 3, sigma_sil: float = 5)[source]
Configuration object for Silicate background positions from Schiavi et al. (2018) for basalts
- Parameters:
here) (spectral range taken as background positions (from left to right as listed)
fit_silicate (str) – Type of fit for the baseline of the silicate region (‘poly’ or ‘spline’)
N_poly_silicate (int) – Degree of polynomial fit for the baseline of the silicate region
sigma_water (int) – Allow points on background within +-sigma_silicate * std dev of other background points
- class DiadFit.H2O_fitting.sil_bck_pos_Schiavi_basanite(lower_range_sil: [Tuple[float, float]] = (340, 360), mid_range1_sil: [Tuple[float, float]] = (630, 640), mid_range2_sil: [Tuple[float, float]] = (nan, nan), upper_range_sil: [Tuple[float, float]] = (1190, 1200), N_poly_sil: float = 3, sigma_sil: float = 5)[source]
Configuration object for Silicate background positions from Schiavi et al. (2018) for basanites
- Parameters:
here) (spectral range taken as background positions (from left to right as listed)
fit_silicate (str) – Type of fit for the baseline of the silicate region (‘poly’ or ‘spline’)
N_poly_silicate (int) – Degree of polynomial fit for the baseline of the silicate region
sigma_water (int) – Allow points on background within +-sigma_silicate * std dev of other background points
- class DiadFit.H2O_fitting.sil_bck_pos_Schiavi_rhyolite(lower_range_sil: [Tuple[float, float]] = (190, 210), mid_range1_sil: [Tuple[float, float]] = (670, 710), mid_range2_sil: [Tuple[float, float]] = (840, 845), upper_range_sil: [Tuple[float, float]] = (1250, 1260), N_poly_sil: float = 3, sigma_sil: float = 5)[source]
Configuration object for Silicate background positions from Schiavi et al. (2018) for rhyolites
- Parameters:
lower_range_sil (Tuple[float, float]) – spectral range taken as background positions (from left to right as listed here)
mid_range1_sil (Tuple[float, float]) – spectral range taken as background positions (from left to right as listed here)
mid_range2_sil (Tuple[float, float]) – spectral range taken as background positions (from left to right as listed here)
upper_range_sil (Tuple[float, float]) – spectral range taken as background positions (from left to right as listed here)
fit_silicate (str) – Type of fit for the baseline of the silicate region (‘poly’ or ‘spline’)
N_poly_silicate (int) – Degree of polynomial fit for the baseline of the silicate region
sigma_water (int) – Allow points on background within +-sigma_silicate * std dev of other background points
- DiadFit.H2O_fitting.smooth_and_trim_around_host(filename=None, x_range=[800, 900], x_max=900, Host_spectra=None, MI_spectra=None, plot_figure=True)[source]
Takes melt inclusion and host spectra, and trims into the region around the host peaks, and fits a cubic spline (used for unmixing spectra)
- Parameters:
x_range (list) – range of x coordinates to smooth between (e.g. [800, 900] by default
Host_spectra (nd.array) – numpy array of host spectra (x is wavenumber, y is intensity)
MI_spectra (nd.array) – numpy array of melt inclusion spectra (x is wavenumber, y is intensity)
filename (str) – name of file for saving figure
Returns:
x_new: x coordinates of smoothed curves y_cub_MI: smoothed y coordinates using a cubic spline for MI y_cub_Host: smoothed y coordinates using a cubic spline for Ol
peak_pos_Host: x coordinates of 2 host peaks peak_height_Host: y coordinates of 2 host peaks trough_x: x coordinate of minimum point between peaks trough_y: y coordinate of minimum point between peaks
- DiadFit.H2O_fitting.stitch_dataframes_together(df_sil, df_water, MI_file, Host_file=None, save_csv=False, path=False)[source]
This function stitches together results from the fit_area function for silicate and water peaks and returns a DataFrame with the combined results. The DataFrame includes peak areas and background positions for both silicate and water peaks, and adds columns for the ratios of water to silicate areas.
- Parameters:
- Returns:
DataFrame with columns for MI filename, Silicate_Trapezoid_Area, Silicate_Simpson_Area = Total Silicate Area using trapezoid or Simpson method LW_Trapezoid_Area, LW_Simpson_Area, MW_Trapezoid_Area, MW_Simpson_Area, HW_Trapezoid_Area, HW_Simpson_Area = Areas of LW, MW and HW Silicate areas (following Shiavi) using Trapezoid or Simpson integration methodWater_Trapezoid_area, Water_Simpson_area = Total area under water peak using Simpson or Trapezoid Integration method. Water_to_HW_ratio_Simpson, Water_to_HW_ratio_Trapezoid: Ratio of Total water area divided by the HW silicate area Water_to_Total_Silicate_ratio_Simpson, Water_to_Total_Silicate_ratio_Trapezoid: Ratio of Total water area divided by the HW silicate area
- Return type:
pd.DataFrame
- DiadFit.H2O_fitting.trough_or_peak_higher(spectra_x, spectra_y, peak_pos_x, trough_pos_x, trough_pos_y, av_width=1, plot=False, print_result=False)[source]
This function assesses whether the line between the 2 peaks is above or below the trough position Called by a loop to select the optimum unmixing ratio for host and melt
- Parameters:
spectra_x (x coordinates of spectra to test)
spectra_y (y coordinates of spectra to test)
peak_pos_x (x positions of 2 host peaks (from find_host_peak_trough_pos))
trough_pos_x (x position of trough (from find_host_peak_trough_pos))
trough_pos_y (y position of trough (from find_host_peak_trough_pos))
av_width (averages +- 1 width either side of the peak and troughs when doing assesment and regression)
plot (bool (default False))
print_result (bool (default False)) – prints whether it found the trough above or below the peak
plot – if True, Draws a plot showing the regression
Returns
-----------
Dist (Vertical height between the linear line (Y coordinate) between the 2 peaks at the x position of the trough,)
trough. (and the y coordinate of the)
- class DiadFit.H2O_fitting.water_bck_pos(fit_water: str = 'poly', N_poly_water: float = 3, lower_bck_water: [Tuple[float, float]] = (2750, 3100), upper_bck_water: [Tuple[float, float]] = (3750, 4100))[source]
Configuration object for water peak and background positions.
- Parameters
- lower_bck_water, upper_bck_water: Tuple[float, float]
2 coordinates for background to left of water peak, and to right.
- fit_waterstr
Type of fit for the baseline of the water region (‘poly’ or ‘spline’)
- N_poly_waterint
Degree of polynomial fit for the baseline of the water region
- sigma_water: int
Allow points on background within +-sigma*water * std dev of other background points
Functions for reading data from the P sensor
- DiadFit.Psensor.add_datetime_and_duration_cols(*, df, raman_cpu_offset='none', offset_hms=[0, 0, 0])[source]
Takes a DataFrame and adds columns for “Date and Time”, “unix_timestamp” and “duration_s”. The input frame should be either the complete DataFrame with spectra metadata and fits output by DiadFit or just the spectral metadata. “Date and Time” contains datetime objects; ‘unix_timestamp’ contains the numeric (float) timestamp for the date and time in standard UNIX time or seconds since epoch time (Jan 1st 1970 00:00:00 UTC), this is plottable; and “duration_s” is the duration of the analysis in seconds.
- Parameters:
df (pd.DataFrame) – The input DataFrame
- Returns:
df – The input DataFrame with additional columns for date and time, duration, and unix timestamp.
- Return type:
pd.DataFrame
- DiadFit.Psensor.get_files(*, path, filetype)[source]
Returns a list of files with specific file type(s) in the specified directory :param path: Path of the directory where the files are located. :type path: str :param filetype: Filetype(s) of the files to be included in the output list. :type filetype: str or list of str
- Returns:
file_ls – A list of files with the specified file type(s) in the directory
- Return type:
- DiadFit.Psensor.get_p_medians(*, pdata, sdata, export_all=False)[source]
Takes two DataFrames and returns a new DataFrame containing the median and median absolute deviation of the pressure values for each Raman analysis. It finds the closest matching rows in the two DataFrames based on timestamp and filters the pressure data between the matched timestamps. It then calculates the median and mean absolute deviation of the filtered pressure data. If export_all==True, it also includes the start time, end time, duration, and filename in the output DataFrame.
- Parameters:
pdata (pd.DataFrame) – The pressure DataFrame (output by read_pfiles)
sdata (pd.DataFrame) – The spectral DataFrame (loaded in by user)
export_all (bool (Optional)) – Indicates whether to include additional information in the output DataFrame.
- DiadFit.Psensor.read_pfiles(*, path, file, start_time, sn_name='0132212')[source]
Reads a csv or xlsx file of pressure data exported from ESI-TEC software and returns a dataframe with two extra columns “Date and Time” (datetime object) and “unix_timestamp” (timestamp expressed as UNIX time, or time in seconds since the epoch time Jan 1st, 1970 00:00:00 UTC) based on the start_time of the pressure recording and time since start in the file. It also renames the time column to Time_sincestart. :param path: Path of the directory where the file is located :type path: str :param file: The name of the file to be read :type file: str :param start_time: The starting time of the recording in the format ‘yyyy-mm-dd hh:mm:ss’, this can be obtained from the docx report by using the report_info function :type start_time: str :param sn_name: The serial number of the sensor. Default is ‘0132212’, this can be obtained from the docx report by using the report_info function. :type sn_name: str
- Returns:
data – DataFrame containing the data from the file along with “Date and Time” and “unix_timestamp” columns.
- Return type:
pd.DataFrame
- DiadFit.Psensor.report_info(*, path, report)[source]
Reads a word document report (exported from ESI-TEC software), extracts and returns the start date and time of the pressure recording and the serial number of the sensor. :param path: Path of the directory where the word document is located :type path: str :param report: The name of the word document :type report: str
- Returns:
start_time (datetime object) – Start time of the analysis in the report
sn_str (str) – Serial number of the sensor
Functions for fitting gas cell data to calibrate a densimeter
- DiadFit.densimeter_fitting.calculate_generic_std_err_values(*, pickle_str, new_x, CI=0.67)[source]
This function loads a model from a pickle file and calculates standard error values based on the model residuals, considering a confidence interval.
- Parameters:
(str) (pickle_str) – The path to the pickle file containing the model.
(np.array) (new_x) – x values for calculation.
(float (CI) – The confidence level for prediction intervals (default: 0.67).
optional) – The confidence level for prediction intervals (default: 0.67).
Returns
pandas.DataFrame (-)
- DiadFit.densimeter_fitting.plot_and_save_CO2cali_pickle(*, cali_data, CO2_dens_col='rho', Split_col='Split', split_error='split_err', CO2_dens_error='dens_err', density_range, N_poly=3, CI=0.67, std_error=True, save_fig=False, eq_division='ccmr', save_suffix='')[source]
This function calculates and saves three polynomial regression models as pickles for CO2 calibration data according to DeVitre et al. 2021 (low density, mid density and high density). The models are saved as pickles and the plot is saved as well if desired.
- Parameters:
(pandas.DataFrame) (cali_data) – A DataFrame containing calibration data (should at least contain a density and a fermi splitting column)
(str (save_suffix) – The column name corresponding to CO2 density in the cali_data DataFrame (default: ‘rho’).
optional) – The column name corresponding to CO2 density in the cali_data DataFrame (default: ‘rho’).
(str – The column name corresponding to Fermi splitting values in the cali_data DataFrame (default: ‘Split’).
optional) – The column name corresponding to Fermi splitting values in the cali_data DataFrame (default: ‘Split’).
(str – The column name corresponding to error in fermi splitting OR an array of splitting errors OR a float (default: ‘split_err’).
float – The column name corresponding to error in fermi splitting OR an array of splitting errors OR a float (default: ‘split_err’).
array-like (float or) – The column name corresponding to error in fermi splitting OR an array of splitting errors OR a float (default: ‘split_err’).
optional) – The column name corresponding to error in fermi splitting OR an array of splitting errors OR a float (default: ‘split_err’).
(str – The column name corresponding to error in CO2 density OR an array of density errors OR a float(default: ‘dens_err’).
array-like – The column name corresponding to error in CO2 density OR an array of density errors OR a float(default: ‘dens_err’).
optional) – The column name corresponding to error in CO2 density OR an array of density errors OR a float(default: ‘dens_err’).
(str) (density_range) – The density range to be fit (‘Low’, ‘Medium’, or ‘High’).
(int (N_poly) – The degree of the polynomial fit (default: 3).
optional) – The degree of the polynomial fit (default: 3).
(float (CI) – The confidence level for prediction intervals (default: 0.67).
optional) – The confidence level for prediction intervals (default: 0.67).
(bool (save_fig) – Whether to calculate and plot standard error (default: True).
optional) – Whether to calculate and plot standard error (default: True).
(bool – Whether to save the plot as an image (default: False).
optional) – Whether to save the plot as an image (default: False).
(str – Method for dividing the data based on density (‘ccmr’, ‘cmass’, or ‘cmass_24C’, default: ‘ccmr’). CCMR corresponds to the limits for each section as shown in DeVitre et al., 2021 (Chem. Geo), cmass is for those in DeVitre et al., 2023 (J. Volcanica).
optional) – Method for dividing the data based on density (‘ccmr’, ‘cmass’, or ‘cmass_24C’, default: ‘ccmr’). CCMR corresponds to the limits for each section as shown in DeVitre et al., 2021 (Chem. Geo), cmass is for those in DeVitre et al., 2023 (J. Volcanica).
(str – Suffix to be added to the saved file names (default: ‘’).
optional) – Suffix to be added to the saved file names (default: ‘’).
Returns – Saves a .pkl file, doesnt return anything.