Skip to content

fixing unit issue in Shkl. Units are now in degrees and correctly scaled to be additive with U and Y respectively. #850

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
merged 1 commit into from
Oct 30, 2025
Merged

fixing unit issue in Shkl. Units are now in degrees and correctly scaled to be additive with U and Y respectively. #850

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
40 changes: 22 additions & 18 deletions hexrd/powder/wppf/peakfunctions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,9 @@ def _gaussian_fwhm(uvw, P, gamma_ani_sqr, eta_mixing, tth, dsp):
tanth = np.tan(th)
cth2 = np.cos(th) ** 2.0
sig2_ani = gamma_ani_sqr * (1.0 - eta_mixing) ** 2 * dsp**4
# sig2_ani is in radians and does not have the 1E4 factor
# built in to U. we will add it here
sig2_ani = np.degrees(sig2_ani) * 1e4
sigsqr = (U + sig2_ani) * tanth**2 + V * tanth + W + P / cth2
if sigsqr <= 0.0:
sigsqr = 1.0e-12
Expand Down Expand Up @@ -99,6 +102,9 @@ def _lorentzian_fwhm(xy, xy_sf, gamma_ani_sqr, eta_mixing, tth, dsp):
tanth = np.tan(th)
cth = np.cos(th)
sig_ani = np.sqrt(gamma_ani_sqr) * eta_mixing * dsp**2
# sig2_ani is in radians and does not have the 1E2 factor
# built in to Y. we will add it here
sig_ani = np.degrees(sig_ani) * 1e2
gamma = (X + xy_sf) / cth + (Y + sig_ani) * tanth
return gamma * 1e-2

Expand Down Expand Up @@ -212,16 +218,17 @@ def _unit_lorentzian(p, x):
f = gamma / ((x - x0) ** 2 + gamma**2)
return f


@njit(cache=True, nogil=True)
def _heaviside(x, x0):
y = np.empty_like(x)
for ii in np.arange(x.size):
if x[ii] < 0.:
y[ii] = 0.
elif x[ii] == 0.:
if x[ii] < 0.0:
y[ii] = 0.0
elif x[ii] == 0.0:
y[ii] = x0
else:
y[ii] = 1.
y[ii] = 1.0

return y

Expand All @@ -245,24 +252,24 @@ def _split_unit_gaussian(p, x):
y : numpy.ndarray
intensity of split gaussian function
'''
x0 = p[0]
x0 = p[0]
fwhm_l = p[1]
fwhm_r = p[2]

sigma_l = fwhm_l / gauss_width_fact
sigma_r = fwhm_r / gauss_width_fact

heav_l = _heaviside(x0-x, 1.0)
heav_r = _heaviside(x-x0, 1.0)
heav_l = _heaviside(x0 - x, 1.0)
heav_r = _heaviside(x - x0, 1.0)

p_l = np.array([x0, fwhm_l])
p_r = np.array([x0, fwhm_r])

gauss_l = _unit_gaussian(p_l, x)
gauss_r = _unit_gaussian(p_r, x)

return (gauss_l*heav_l +
gauss_r*heav_r)
return gauss_l * heav_l + gauss_r * heav_r


# =========================================================
# 1-D split pseudo-voight functions
Expand Down Expand Up @@ -299,8 +306,8 @@ def _split_unit_pv(p, x):
fwhm_g_r = p[3]
fwhm_l_r = p[4]

heav_l = _heaviside(x0-x, 1.0)
heav_r = _heaviside(x-x0, 1.0)
heav_l = _heaviside(x0 - x, 1.0)
heav_r = _heaviside(x - x0, 1.0)

eta_l, fwhm_l = _mixing_factor_pv(fwhm_g_l, fwhm_l_l)
eta_r, fwhm_r = _mixing_factor_pv(fwhm_g_r, fwhm_l_r)
Expand All @@ -313,16 +320,17 @@ def _split_unit_pv(p, x):
gamma_r = fwhm_r / lorentz_width_fact

g_l = _unit_gaussian(np.array([x0, fwhm_l]), x)
l_l = _unit_lorentzian(np.array([x0, fwhm_l]), x)*gamma_l
l_l = _unit_lorentzian(np.array([x0, fwhm_l]), x) * gamma_l

g_r = _unit_gaussian(np.array([x0, fwhm_r]), x)
l_r = _unit_lorentzian(np.array([x0, fwhm_r]), x)*gamma_r
l_r = _unit_lorentzian(np.array([x0, fwhm_r]), x) * gamma_r

pv_l = eta_l * l_l + (1.0 - eta_l) * g_l

pv_r = eta_r * l_r + (1.0 - eta_r) * g_r

return pv_l*heav_l + pv_r*heav_r
return pv_l * heav_l + pv_r * heav_r


@njit(cache=True, nogil=True)
def _mixing_factor_pv(fwhm_g, fwhm_l):
Expand Down Expand Up @@ -390,7 +398,6 @@ def _func_h(tau, tth_r):

@njit(cache=True, nogil=True)
def _func_W(HoL, SoL, tau, tau_min, tau_infl, tth):

if tth < np.pi / 2.0:
if tau >= 0.0 and tau <= tau_infl:
res = 2.0 * min(HoL, SoL)
Expand Down Expand Up @@ -615,7 +622,6 @@ def computespectrum_pvfcj(
np.array([Iobs.shape[0], tth.shape[0], dsp.shape[0], hkl.shape[0]])
)
for ii in prange(nref):

II = Iobs[ii]
t = tth[ii]
d = dsp[ii]
Expand Down Expand Up @@ -648,7 +654,6 @@ def computespectrum_pvtch(
np.array([Iobs.shape[0], tth.shape[0], dsp.shape[0], hkl.shape[0]])
)
for ii in prange(nref):

II = Iobs[ii]
t = tth[ii]
d = dsp[ii]
Expand Down Expand Up @@ -691,7 +696,6 @@ def computespectrum_pvpink(
np.array([Iobs.shape[0], tth.shape[0], dsp.shape[0], hkl.shape[0]])
)
for ii in prange(nref):

II = Iobs[ii]
t = tth[ii]
d = dsp[ii]
Expand Down
Loading