import sandy
import pandas as pd
from os.path import join
import numpy as np
import random, sys
za = 92235
tape = sandy.get_endf6_file("jeff_33", "nfpy", za * 10)
nfpy = sandy.Fy.from_endf6(tape)
energies = nfpy.data.E.unique()
nsmp = 10 # sample size
smp = {}
for e in energies:
idx = nfpy.data.query(f"E=={e} & MT==454").index
ify = nfpy.data.loc[idx]
cov = sandy.CategoryCov(np.diag(ify.DFY**2), index=ify.ZAP, columns=ify.ZAP) # Diagonal covariance matrix
seed = random.randrange(2**32 - 1) # create a seed
print(f"sampling IFY for energy {e:.3e} eV...")
smp[e] = cov.sampling(nsmp, seed=seed) # need to change the seed for the different energies
WARNING: Large condition number of covariance matrix: 1.36e+11 WARNING: Large condition number of covariance matrix: 1.38e+11
sampling IFY for energy 2.530e-02 eV... sampling IFY for energy 4.000e+05 eV... sampling IFY for energy 1.400e+07 eV...
WARNING: Large condition number of covariance matrix: 1.28e+11
with pd.ExcelWriter(f'PERT_{za}_MF8_MT454.xlsx') as writer:
for e, s in smp.items():
s.data.to_excel(writer, sheet_name=f'{e:.3e}')
Skip the part above if you already have the file of perturbations.
smp = pd.read_excel(f'PERT_{za}_MF8_MT454.xlsx', sheet_name=None, index_col=0)
za = 92235
tape = sandy.get_endf6_file("jeff_33", "nfpy", za * 10)
nfpy = sandy.Fy.from_endf6(tape)
### run only if you want consistent CFYs
# tape_rdd = sandy.get_endf6_file("jeff_33", "decay", "all")
# rdd = sandy.DecayData.from_endf6(tape_rdd) # this can take a while
smp_min = 0 # write ENDF-6 file only in the sample range [smp_min, smp_max]
smp_max = 9
file_template = "u235_fy_{}.jeff33"
for ismp in range(smp_min, smp_max+1):
file = file_template.format(ismp)
f = sandy.Fy(nfpy.data.copy())
for e, s in smp.items():
idx_ify = nfpy.data.query(f"E=={float(e)} & MT==454").index
idx_cfy = nfpy.data.query(f"E=={float(e)} & MT==459").index
f.data.loc[idx_ify, "DFY"] = f.data.loc[idx_ify, "FY"] # just for me, i copy the original IFYs where uncertainties should be, so i can compare them to the perturbed ones (anyways I don't use uncertainties)
f.data.loc[idx_cfy, "DFY"] = f.data.loc[idx_cfy, "FY"] # same but for CFYs
f.data.loc[idx_ify, "FY"] *= s[ismp].values # IMPORTANT, this does not update the CFYs, which in random ENDF-6 file are inconsistent with the perturbed IFYs
#f = f.apply_qmatrix(922350, e, rdd, keep_fy_index=True) # Run this if you want to update the CFYs (slower), or else comment it out
print(f"writing file '{file}'...")
f.to_endf6(tape).to_file(file)
writing file 'u235_fy_0.jeff33'... writing file 'u235_fy_1.jeff33'... writing file 'u235_fy_2.jeff33'... writing file 'u235_fy_3.jeff33'... writing file 'u235_fy_4.jeff33'... writing file 'u235_fy_5.jeff33'... writing file 'u235_fy_6.jeff33'... writing file 'u235_fy_7.jeff33'... writing file 'u235_fy_8.jeff33'... writing file 'u235_fy_9.jeff33'...