Skip to content

Finalize cone scripts #270

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
KarenDeLannoye merged 1 commit into from
Feb 24, 2026
Merged

Finalize cone scripts #270

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 22 additions & 47 deletions Scripts/MaCFP-4/Cone_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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]')
Expand Down Expand Up @@ -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)
Expand All @@ -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('_')
Expand All @@ -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))

Expand Down Expand Up @@ -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'}
Expand Down Expand Up @@ -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
)

Expand Down
56 changes: 41 additions & 15 deletions Scripts/MaCFP-4/MCC_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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()

Expand Down Expand Up @@ -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))
Expand Down
5 changes: 5 additions & 0 deletions Scripts/MaCFP-4/Utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down