diff --git a/Scripts/MaCFP-4/Cone_analysis.py b/Scripts/MaCFP-4/Cone_analysis.py index 64cb448..5aba53f 100644 --- a/Scripts/MaCFP-4/Cone_analysis.py +++ b/Scripts/MaCFP-4/Cone_analysis.py @@ -12,7 +12,7 @@ from scipy.signal import savgol_filter from Utils import device_data, get_series_names, make_institution_table, device_subset, label_def, format_latex -from Utils import format_with_uncertainty, format_temperature, format_regular, extract_heating_rate, extract_atmosphere, get_condition_key +from Utils import format_ignition, format_regular, extract_heating_rate, extract_atmosphere, get_condition_key from Utils import DATA_DIR @@ -255,7 +255,7 @@ def Calculate_dm_dt(df:pd.DataFrame): zorder =1 else: zorder =5 - ax2.plot(df['Time (s)'], df['HRR (kW/m2)'], '.', label = label, color=color, zorder=zorder) + ax2.plot(df['Time (s)'], df['HRR (kW/m2)'], '-', label = label, color=color, zorder=zorder) ax1.set_ylim(bottom=0) ax1.set_xlabel('Time [s]') @@ -382,13 +382,13 @@ def Calculate_dm_dt(df:pd.DataFrame): for i, path in enumerate(paths): df = pd.read_csv(path) df = calculate_int_HRR(df) - ax1.plot(df['Time (s)'], df['HRR (kW/m2)'], '.', color = color[flux], alpha=0.7, markersize = 0.1, zorder=4) + ax1.plot(df['Time (s)'], df['HRR (kW/m2)'], '.', color = color[flux], alpha=0.5, markersize = 0.1, zorder=5) df_average = average_cone_series(series) - ax1.plot(df_average['Time (s)'], df_average['HRR (kW/m2)'], label = flux + '/m$^2$', color = color[flux], zorder = 3) + ax1.plot(df_average['Time (s)'], df_average['HRR (kW/m2)'], label = flux + '/m$^2$', color = color[flux], zorder = 2) ax1.fill_between(df_average['Time (s)'], df_average['HRR (kW/m2)']-2*df_average['unc HRR (kW/m2)'], df_average['HRR (kW/m2)']+2*df_average['unc HRR (kW/m2)'], - color=color[flux], alpha = 0.3, zorder=2) + color=color[flux], alpha = 0.3, zorder=3) ax1.set_ylim(bottom=0,top=250) ax1.set_xlim(right=2500) @@ -402,7 +402,7 @@ def Calculate_dm_dt(df:pd.DataFrame): # Back side temperature plots for all unique atmospheres and heating rates (when available) -linestyle = ['-','--',':'] +linestyle = ['--',':','-'] for series in unique_conditions_cone_material: fig1, ax1 = plt.subplots(figsize=(6, 4)) parts = series.split('_') @@ -416,13 +416,27 @@ def Calculate_dm_dt(df:pd.DataFrame): temp_col = f'TC back {i} (K)' if temp_col in df.columns: ax1.plot(df['Time (s)'], df[temp_col], label=label, color=color, linestyle = linestyle[i-1]) + if 'TC Top (K)' in df.columns: + ax1.plot(df['Time (s)'], df['TC Top (K)'], label=label, color=color, linestyle = '-') ax1.set_ylim(bottom=250) ax1.set_xlabel('Time [s]') ax1.set_ylabel('Temperature [K]') fig1.tight_layout() - ax1.legend() + + # Get unique (label, style) combinations + handles, labels = ax1.get_legend_handles_labels() + unique = {} + for handle, label in zip(handles, labels): + # Create a key based on label and visual properties + key = (label, handle.get_linestyle(), handle.get_color()) + if key not in unique: + unique[key] = (handle, label) + + unique_handles = [v[0] for v in unique.values()] + unique_labels = [v[1] for v in unique.values()] + ax1.legend(unique_handles, unique_labels) fig1.savefig(str(base_dir) + '/Cone/Cone_{}_{}_{}_BackT.{}'.format(material, flux,orient,ex)) @@ -506,45 +520,6 @@ def Calculate_dm_dt(df:pd.DataFrame): -# Indivdual Mass and mass loss rate plots for all unique institutions atmospheres and heating rates (gasification) -for set in Gas_sets: - fig1, ax1 = plt.subplots(figsize=(6, 4)) - fig2, ax2 = plt.subplots(figsize=(6, 4)) - paths_GAS_set = list(DATA_DIR.glob(f"*/{set}_[rR]*.csv")) - for path in paths_GAS_set: - df_raw = pd.read_csv(path) - df=Calculate_dm_dt(df_raw) - label, color = label_def(path.stem.split('_')[0]) - if institute == 'TIFP+UCT': - ax1.plot(df['Time (s)'],savgol_filter(df['dm/dt']/0.01,41,3),'-', label = label, color=color) - elif institute == 'FSRI': - ax1.plot(df['Time (s)'],savgol_filter(df['dm/dt']/0.00385,41,3),'-', label = label, color=color) - ax2.plot(df['Time (s)'], df['Mass (g)'], '.', label = label, color=color) - - ax1.set_ylim(bottom=0) - ax1.set_xlabel('Time [s]') - ax1.set_ylabel('Mass loss rate [g s$^{-1}$ m$^{-2}$]') - fig1.tight_layout() - handles1, labels1 = ax1.get_legend_handles_labels() - by_label1 = dict(zip(labels1, handles1)) - ax1.legend(by_label1.values(), by_label1.keys()) - - ax2.set_ylim(bottom=0) - ax2.set_xlabel('Time [s]') - ax2.set_ylabel('Mass [g]') - fig2.tight_layout() - handles2, labels2 = ax2.get_legend_handles_labels() - by_label2 = dict(zip(labels2, handles2)) - ax2.legend(by_label2.values(), by_label2.keys()) - - fig1.savefig(str(base_dir) + '/Cone/Gasification_{}_{}_MLR.{}'.format(material, flux,ex)) - fig2.savefig(str(base_dir) + '/Cone/Gasification_{}_{}_Mass.{}'.format(material, flux,ex)) - - - plt.close(fig1) - plt.close(fig2) - - # Plots comparing parallel verus perpendicular color = {'perpendicular':'black', 'parallel':'red'} @@ -716,7 +691,7 @@ def Calculate_dm_dt(df:pd.DataFrame): # Format ignition time final_table_sorted['ignitiont_formatted'] = final_table_sorted.apply( - lambda row: format_regular(row['ignition time'], row['std ignition time']), + lambda row: format_ignition(row['ignition time'], row['std ignition time']), axis=1 ) diff --git a/Scripts/MaCFP-4/MCC_analysis.py b/Scripts/MaCFP-4/MCC_analysis.py index 7861088..e97a42a 100644 --- a/Scripts/MaCFP-4/MCC_analysis.py +++ b/Scripts/MaCFP-4/MCC_analysis.py @@ -95,7 +95,7 @@ def calculate_int_HRR(df:pd.DataFrame) -> pd.DataFrame: total_hrr = np.zeros(len(df)) for i in range(1, len(df)): total_hrr[i] = total_hrr[i-1] + 0.5 * (df['HRR (W/g)'].iloc[i-1] + df['HRR (W/g)'].iloc[i]) * (df['Time (s)'].iloc[i] - df['Time (s)'].iloc[i-1]) - df['Int HRR'] = total_hrr + df['Int HRR'] = total_hrr/1000 return df @@ -151,15 +151,41 @@ def average_MCC_series(series_name: str, exclude:Optional[Union[str, List[str]]] dTdt_cols = merged_df.filter(regex=r'^dTdt').columns intHRR_cols = merged_df.filter(regex=r'^Int HRR').columns - df_average = pd.DataFrame({ - 'Temperature (K)': merged_df['Temperature (K)'], - 'HRR (W/g)': merged_df[hrr_cols].mean(axis=1), - 'HRR_std': merged_df[hrr_cols].std(axis=1, skipna=True, ddof=0), - 'dTdt (K/min)': merged_df[dTdt_cols].mean(axis=1), - 'dTdt_std': merged_df[dTdt_cols].std(axis=1, skipna=True, ddof=0), - 'int HRR': merged_df[intHRR_cols].mean(axis=1), - 'int HRR_std': merged_df[intHRR_cols].std(axis=1, skipna=True, ddof=0), - }).dropna(subset=['HRR (W/g)'], how='all') + df_average = pd.DataFrame({'Temperature (K)': merged_df['Temperature (K)']}) + n=2 + sum = merged_df[hrr_cols].rolling(2*n+1, min_periods=1,center=True).sum().sum(axis=1) + cnt = merged_df[hrr_cols].rolling(2*n+1, min_periods=1,center=True).count().sum(axis=1) + df_average['HRR (W/g)'] = sum / cnt # Series: mean of all non-NaN values in rows i-2..i+2 across all columns + + diff = merged_df[hrr_cols].sub(df_average['HRR (W/g)'], axis=0)**2 + sum_diff = diff.rolling(2*n+1, min_periods=1,center=True).sum().sum(axis=1) + df_average['HRR_std'] = np.sqrt(sum_diff/(cnt*(cnt-1))) + + sum = merged_df[dTdt_cols].rolling(2*n+1, min_periods=1,center=True).sum().sum(axis=1) + cnt = merged_df[dTdt_cols].rolling(2*n+1, min_periods=1,center=True).count().sum(axis=1) + df_average['dTdt (K/min)'] = sum / cnt # Series: mean of all non-NaN values in rows i-2..i+2 across all columns + + diff = merged_df[dTdt_cols].sub(df_average['dTdt (K/min)'], axis=0)**2 + sum_diff = diff.rolling(2*n+1, min_periods=1,center=True).sum().sum(axis=1) + df_average['dTdt_std'] = np.sqrt(sum_diff/(cnt*(cnt-1))) + + sum = merged_df[intHRR_cols].rolling(2*n+1, min_periods=1,center=True).sum().sum(axis=1) + cnt = merged_df[intHRR_cols].rolling(2*n+1, min_periods=1,center=True).count().sum(axis=1) + df_average['int HRR'] = sum / cnt # Series: mean of all non-NaN values in rows i-2..i+2 across all columns + + diff = merged_df[intHRR_cols].sub(df_average['int HRR'], axis=0)**2 + sum_diff = diff.rolling(2*n+1, min_periods=1,center=True).sum().sum(axis=1) + df_average['int HRR_std'] = np.sqrt(sum_diff/(cnt*(cnt-1))) + + # df_average = pd.DataFrame({ + # 'Temperature (K)': merged_df['Temperature (K)'], + # 'HRR (W/g)': merged_df[hrr_cols].mean(axis=1), + # 'HRR_std': merged_df[hrr_cols].std(axis=1, skipna=True, ddof=0), + # 'dTdt (K/min)': merged_df[dTdt_cols].mean(axis=1), + # 'dTdt_std': merged_df[dTdt_cols].std(axis=1, skipna=True, ddof=0), + # 'int HRR': merged_df[intHRR_cols].mean(axis=1), + # 'int HRR_std': merged_df[intHRR_cols].std(axis=1, skipna=True, ddof=0), + # }).dropna(subset=['HRR (W/g)'], how='all') return df_average @@ -217,7 +243,7 @@ def average_MCC_series(series_name: str, exclude:Optional[Union[str, List[str]]] ax2.set_ylim(bottom=0) ax2.set_xlim(right=900) ax2.set_xlabel('Temperature (K)') - ax2.set_ylabel('Integral HRR [J g$^{-1}$]') + ax2.set_ylabel('Integral HRR [kJ g$^{-1}$]') fig2.tight_layout() handles2, labels2 = ax2.get_legend_handles_labels() by_label2 = dict(zip(labels2, handles2)) @@ -355,7 +381,7 @@ def average_MCC_series(series_name: str, exclude:Optional[Union[str, List[str]]] T_onset_list.append(T_onset) T_onset10_list.append(T_onset10) FGC_list.append(FGC_v) - HR_total_list.append(HR_total/1000) #kJ/g + HR_total_list.append(HR_total) #kJ/g HR_capacity_list.append(HR_Capacity) ax_HRR.plot(df['Temperature (K)'], df['HRR (W/g)'], '.',color ='black',markersize=0.00000000000002) @@ -382,7 +408,7 @@ def average_MCC_series(series_name: str, exclude:Optional[Union[str, List[str]]] # Axes labels ax_HRR.set_xlabel('Temperature (K)') ax_HRR.set_ylabel('HRR [W g$^{-1}$]') - ax_intHRR.set_ylabel('Integral HRR [J g$^{-1}$]') + ax_intHRR.set_ylabel('Integral HRR [kJ g$^{-1}$]') # Figure title plt.title(Duck+"\n"+Conditions) @@ -497,7 +523,7 @@ def plot_hrr_and_onset_vs_peak_temp(df): ax2.set_ylim(bottom=0) ax2.set_xlim(350,1000) ax2.set_xlabel('Temperature (K)') -ax2.set_ylabel('Integral HRR [J/g]') +ax2.set_ylabel('Integral HRR [kJ/g]') fig2.tight_layout() ax2.legend() @@ -611,7 +637,7 @@ def plot_hrr_and_onset_vs_peak_temp(df): ax2.set_ylim(bottom=0) ax2.set_xlim(350, 1000) ax2.set_xlabel('Temperature (K)') -ax2.set_ylabel('Integral HRR [J/g]') +ax2.set_ylabel('Integral HRR [kJ/g]') # Remove duplicate legend entries handles2, labels2 = ax2.get_legend_handles_labels() by_label2 = dict(zip(labels2, handles2)) diff --git a/Scripts/MaCFP-4/Utils.py b/Scripts/MaCFP-4/Utils.py index 03f6415..1ec970e 100644 --- a/Scripts/MaCFP-4/Utils.py +++ b/Scripts/MaCFP-4/Utils.py @@ -274,6 +274,11 @@ def format_regular(mean, std): return f"${mean:.2f}$" else: return f"${mean:.2f} \\pm {std:.2f}$" + + +def format_ignition(mean, std): + """Format regular numbers (no scientific notation)""" + return f"${mean:.2f} \\pm {std:.2f}$" # Extract sorting keys from conditions def extract_heating_rate(conditions):