SasView / sasmodels

Package for calculation of small angle scattering models using OpenCL.
BSD 3-Clause "New" or "Revised" License
16 stars 28 forks source link

profile plots for polydisperse models #584

Closed pkienzle closed 1 year ago

pkienzle commented 1 year ago

With polydispersity there is a distribution of profiles. How can this be represented on a plot?

A simple approach is to generate many different profiles and over-plot them. Here is an example for a core-shell sphere:

import numpy as np
from sasmodels.weights import get_weights
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

def cs_sphere_profile(radius, thickness, sld_core, sld_shell, sld_solvent):
    z = [0, radius, radius, radius+thickness, radius+thickness, (radius+thickness)*1.25]
    rho = [sld_core, sld_core, sld_shell, sld_shell, sld_solvent, sld_solvent]
    return np.asarray(z), np.asarray(rho)

radius, radius_pd = 40, 0.2
thickness, thickness_pd = 70, 0.1
sld_core, sld_shell, sld_solvent = 2.07, 4.5, -0.5

nz, nr, nt, nsigmas = 1000, 150, 25, 3
r, Pr = get_weights(
    disperser='gaussian',
    n=nr, width=radius_pd, nsigmas=nsigmas,
    value=radius, limits=(0, np.inf), relative=True)
t, Pt = get_weights(
    disperser='gaussian',
    n=nt, width=thickness_pd, nsigmas=nsigmas,
    value=thickness, limits=(0, np.inf), relative=True)

max_z = r[-1]+t[-1]
z = np.linspace(0, max_z, nz)

nlines = nr*nt
rhos = np.empty((nlines, nz))
weights = np.empty(nlines)
index = 0
for rj, Prj in zip(r, Pr):
    for tk, Ptk in zip(t, Pt):
        zjk, rhojk = cs_sphere_profile(rj, tk, sld_core, sld_shell, sld_solvent)
        rhos[index] = np.interp(z, zjk, rhojk) # default is rhojk[-1] for z > zjk[i]
        weights[index] = Prj*Ptk
        index += 1
rho_avg = np.average(rhos, axis=0, weights=weights)

# plot individual lines with transparency proportional to weights
colors = np.zeros((nlines, 4))
colors[:,2] = 1 # blue
colors[:,3] = weights*15
lines = np.empty((nlines, nz, 2))
lines[:,:,0] = z[None, :]
lines[:,:,1] = rhos
h = LineCollection(lines, colors=colors)
plt.plot(z, rho_avg)
plt.gca().add_collection(h)

The result is the following plot:

image
pkienzle commented 1 year ago

Closing this ticket because I don't see this as useful in general.

Here's a core-multishell model with large PD on the core. The details on the shells are smoothed away in the profile even though the multilayer shell clearly shows up as a peak in the scattering.

Shifting the profile to align on the core-shell boundary or on the center of the third shell would introduce other artifacts.

import numpy as np
from scipy.stats import truncnorm
from sasmodels.weights import get_weights
from sasmodels.models.core_multi_shell import profile
import matplotlib.pyplot as plt
%matplotlib inline

sld_core, sld_solvent = 2.2, 2.2
radius, radius_pd = 564, 0.3
thickness = np.array([14.8, 14.9, 5.7, 14.9, 14.8])
thickness_pd = 0.1
sld = np.array([2.0, 2.27, 0.15, 2.27, 2.0])
n = len(thickness)

# Cheating because all thickness have the same pd. Otherwise we need
# to scale each thickness by 1 + 3*pd to get the max possible value
z_max = radius*(1 + 3*radius_pd) + sum(thickness)*(1 + 3*thickness_pd)

nlines, nz = 50000, 200
z = np.linspace(0, z_max, nz)
rhos = np.empty((nlines, nz))
pars = np.empty((nlines, n+1))
for k in range(nlines):
    # Truncated normal defines bounds using number of sigmas from the mean.
    # That is, the lower limit with mean + sigma*a = 0 requires a = -mean/sigma.
    # For polydisperse parameters, sigma = pd*mean, therefore a = -1/pd
    rk = truncnorm.rvs(-1/radius_pd, np.inf, radius, radius*radius_pd)
    tk = truncnorm.rvs(-1/thickness_pd, np.inf, thickness, thickness*thickness_pd)
    pars[k, 0] = rk
    pars[k, 1:] = tk
    zk, rhok = profile(sld_core, rk, sld_solvent, n, sld, tk)
    rhos[k] = np.interp(z, zk, rhok)
rho_avg = np.mean(rhos, axis=0)

# Average profile, dominated by pd in radius
nplot = min(400, nlines)
plt.plot(z, rho_avg, color='b')

# Fuzziness plot
plt.figure()
plt.plot(z, rhos[:nplot].T, color='b', alpha=0.1)

# Histograms of the individual parameters, to check that they are truncated normal
plt.figure()
for k in range(n+1):
    plt.subplot(2,3,k+1)
    plt.hist(pars[:, k], bins=50)
None
image image image