Aggregate and chain structure-factor models

sasmodels/models/compact_polydisperse_cluster.c

@@ -0,0 +1,144 @@
+/* SCHULZ ------------------------------------------------------------ */
+double SCHULZ(double R, double RA, double Z)
+{
+    double dum =
+        (Z+1.0)*log(fabs((Z+1.0)/RA))
+      + Z*log(fabs(R))
+      - (Z+1.0)*R/RA
+      - sas_gammaln(Z+1.0);
+
+    return exp(dum);
+}
+
+/* S_HS -------------------------------------------------------------- */
+double S_HS(double Q, double RHS, double ETA)
+{
+    double ALN = pow(1.0 - ETA, 4.0);
+    double AL = pow(1.0 + 2.0*ETA, 2.0) / ALN;
+    double BE = -6.0 * ETA * pow(1.0 + 0.5*ETA, 2.0) / ALN;
+    double GA = 0.5 * ETA * AL;
+
+    double AR = 2.0 * RHS * Q;
+    double GG;
+
+    if (AR < 0.4)
+    {
+        GG =  AL*(1.0/3.0 - AR*AR/30.0)
+            + BE*(1.0/4.0 - AR*AR/36.0)
+            + GA*(1.0/6.0 - AR*AR/48.0);
+    }
+    else
+    {
+        double SA = sin(AR);
+        double CA = cos(AR);
+
+        GG =  AL*(SA - AR*CA)/pow(AR,3.0)
+            + BE*(2.0*AR*SA + (2.0 - AR*AR)*CA - 2.0)/pow(AR,4.0)
+            + GA*(-pow(AR,4.0)*CA
+                 + 4.0*((3.0*AR*AR - 6.0)*CA
+                      + (pow(AR,3.0) - 6.0*AR)*SA + 6.0)) / pow(AR,6.0);
+    }
+
+    return 1.0 / (1.0 + 24.0*ETA*GG);
+}
+/* PSCHULZ ----------------------------------------------------------- */
+double PSCHULZ(double Q, double RAV, double SIGMA)
+{
+    double Z = 1.0/(SIGMA*SIGMA) - 1.0;
+    double pi = acos(-1.0);
+    double FIPI3 = 4.0*pi/3.0;
+
+    double ALP  = (Z+1.0)/(Q*RAV);
+    double ALP1 = (Q*RAV)/(Z+1.0);
+    double ATTA = atan(2.0*ALP1);
+
+    double PSI1 = (Z+1.0)*ATTA;
+    double PSI2 = (Z+2.0)*ATTA;
+    double PSI3 = (Z+3.0)*ATTA;
+
+    double ALP2   = ALP*ALP;
+    double ALP2_4 = 1.0/(4.0 + ALP2);
+
+    double C1 = pow(ALP2*ALP2_4,0.5*Z+0.5)*pow(ALP/(Z+1.0),6.0);
+    double T1 = pow(ALP/(Z+1.0),6.0) - C1*cos(PSI1);
+
+    double C2 = pow(ALP2*ALP2_4,0.5*Z+1.5)*pow(ALP/(Z+1.0),4.0);
+    double T2 = (Z+2.0)/(Z+1.0)*( pow(ALP/(Z+1.0),4.0) + C2*cos(PSI3) );
+
+    double C3 = pow(ALP2*ALP2_4,0.5*Z+1.0)*pow(ALP/(Z+1.0),5.0);
+    double T3 = -2.0*C3*sin(PSI2);
+
+    double V2 = pow(FIPI3,2.0) *(Z+6.0)*(Z+5.0)*(Z+4.0)*
+                (Z+3.0)*(Z+2.0)/pow(Z+1.0,5.0);
+
+    return (8.0*pi*pi*(T1+T2+T3))/V2;
+}
+
+/* P_POLY ------------------------------------------------------------ */
+double P_POLY(double Q, double RL, double SIGL, double d, double n_aggreg)
+{
+    double PI43 = 4.1887892;
+    double RS = d;
+    double RDIF = RL - RS;
+
+    double Z = 1.0/(SIGL*SIGL) - 1.0;
+
+    double SVL  = PI43*pow(RDIF,3.0)*(Z+3.0)*(Z+2.0)/((Z+1.0)*(Z+1.0));
+    double SVL2 = PI43*PI43 * pow(RDIF,6.0)
+                  * (Z+6.0)*(Z+5.0)*(Z+4.0)*(Z+3.0)*(Z+2.0)
+                  / pow(Z+1.0,5.0);
+
+    double VS   = PI43*pow(RS,3.0);
+    double ETA = n_aggreg * (SVL * VS)/SVL2;
+    double ALPP = ETA/VS;
+
+    int NPOI = 100;
+    double RMAX = RS + 6.0*RL*SIGL;
+    double DELR = RMAX/NPOI;
+
+    double SSQ = 0.0;
+    double SSQN= 0.0;
+
+    for(int i=1;i<=NPOI;i++)
+    {
+        double R  = RS + (i-0.5)*DELR;
+        double RD = R - RS;
+        double D  = SCHULZ(RD, RDIF, Z);
+
+        double RH = R - 2.0*RS;
+        if(RH < 0) RH = 0.0;
+
+        double FR = (pow(RD,3.0) - pow(RH,3.0))/pow(RD,3.0);
+    if(FR < 0.0) FR = 0.0;
+    if(FR > 1.0) FR = 1.0;
+
+    double ETAEFF = ETA*(FR*0.5 + (1.0-FR));
+
+        SSQ  += D * pow(RD,3.0) * S_HS(Q,RS,ETAEFF);
+        SSQN += D * pow(RD,3.0) + 1.0e-10;
+    }
+
+    SSQ = SSQ/SSQN;
+
+    double PS   = square(VS * sas_3j1x_x(Q*RS));
+    double PLAV = PSCHULZ(Q, RDIF, SIGL);
+
+    double RN  = SVL  * ALPP;
+    double RNX = sqrt(SVL2)*ALPP;
+
+    return (RN*PS*SSQ + PS*RNX*(RNX-1.0)*PLAV)/(RN*PS);
+}
+
+/* ----------- SASVIEW ENTRY POINT ----------- */
+double SQ_COMPACT(double Q, double RL, double SIGL, double d, double n_aggreg)
+{
+    return P_POLY(Q, RL, SIGL, d, n_aggreg);
+}
+
+/* SasView required function */
+double Iq(double Q, double radius_effective, double volfraction,
+          double radius_cluster, double radius_cluster_sigma_relative, double n_aggreg)
+{
+    (void)volfraction;
+    return SQ_COMPACT(Q, radius_cluster, radius_cluster_sigma_relative, radius_effective, n_aggreg);
+}

sasmodels/models/compact_polydisperse_cluster.py

@@ -0,0 +1,89 @@
+r"""
+Definition
+----------
+Structure factor for a polydisperse spherical cluster with internal correlations.
+
+WARNING: Should not be used in combination with very anisotropic particle shapes.
+
+Polydispersity described by Schulz distribution
+and expressions taken from Pedersen, J. S., Møller, T. L., Raak, N., & Corredig,
+M. (2022). A model on an absolute scale for the small-angle X-ray scattering
+from bovine casein micelles. Soft Matter, 18(45), 8613-8625.
+
+The implemented expressions are described in
+Pedersen, J. S., Møller, T. L., & Corredig, M. (2026). Scattering from 'Babinet 'particles (or not…):
+spherical particles made up of spheres and spherical particles with s
+pherical voids. Applied Crystallography, 59(1).
+
+Modifications have been done for adapting to monodisperse distance between
+points in cluster and for allowing weight-average aggregation number to be a fit parameter.
+
+See also
+Larsen, A. H., Pedersen, J. S., & Arleth, L. (2020). Assessment of
+structure factors for analysis of small-angle scattering data from
+desired or undesired aggregates. Applied Crystallography, 53(4), 991-1005.
+
+
+Parameters
+----------
+radius_effective : effective scatterer radius (Å), half the minimum center-to-center distance; first parameter for :math:`P@S` wiring (see sasmodels structure-factor conventions)
+volfraction : unused in this :math:`S(q)`; required as second parameter for :math:`P@S` products
+radius_cluster : mean radius of large clusters
+radius_cluster_sigma_relative : relative polydispersity of radius_cluster
+n_aggreg : weight-average aggregation number
+
+References
+----------
+
+# Pedersen, J. S., Møller, T. L., & Corredig, M. (2026). Scattering from 'Babinet 'particles (or not…):
+spherical particles made up of spheres and spherical particles with s
+pherical voids. Applied Crystallography, 59(1).
+
+Authorship and Verification
+----------------------------
+
+* **Author:** Jan Skov Pedersen
+* **Last Modified by:** Jan Skov Pedersen, April 12, 2026
+* **Last Reviewed by:**  **Date:**
+"""
+
+
+import numpy as np
+
+name = "compact_polydisperse_cluster"
+title = "Compact Polydisperse Cluster Structure Factor"
+description = """
+Structure factor for a compact polydisperse cluster
+with internal correlations.
+"""
+
+category = "structure-factor"
+structure_factor = True
+single = False
+
+# Must match C: Iq(Q, radius_effective, volfraction, radius_cluster, radius_cluster_sigma_relative, n_aggreg)
+parameters = [
+    ["radius_effective", "Ang", 10.0, [0.0, np.inf], "",
+     "effective scatterer radius (half center-to-center distance)"],
+    ["volfraction", "", 0.2, [0.0, 1.0], "",
+     "unused in S(q); required for P@S products"],
+    ["radius_cluster", "Ang", 40.0, [0.0, np.inf], "", "Average cluster radius"],
+    ["radius_cluster_sigma_relative", "", 0.4, [0.0, 1.0], "", "Relative size polydispersity"],
+    ["n_aggreg", "", 50.0, [10.0, 100.0], "", "Weight-average aggregation number"],
+]
+
+source = ["lib/sas_gammainc.c", "lib/sas_3j1x_x.c", "compact_polydisperse_cluster.c"]
+
+def random():
+    import random
+
+    return {
+        "radius_effective": random.uniform(2.5, 15.0),
+        "volfraction": random.uniform(0.01, 0.3),
+        "radius_cluster": random.uniform(50, 200),
+        "radius_cluster_sigma_relative": random.uniform(0.01, 0.3),
+        "n_aggreg": random.uniform(10, 100),
+    }
+
+def test():
+    print("compact_polydisperse_cluster model loaded successfully.")

sasmodels/models/fractal_aggregate.c

@@ -0,0 +1,59 @@
+/*
+ * Fractal structure factor S(q)
+ * Normalization: S(0) = N
+ *
+ * Parameters:
+ *   q  : scattering vector
+ *   radius_effective : effective scatterer radius (half center-to-center distance)
+ *   volfraction      : unused (required for P@S parameter order)
+ *   fractal_dim          : fractal dimension
+ *   n_aggreg            : number of particles in cluster
+ */
+
+static double fractal_sq_N(double q, double r, double D, double N)
+{
+    const double eps = 1e-12;
+
+    /* Special cases */
+    if (N <= 1.0) return 1.0;     /* no clustering */
+    if (D <= 0.0) return 1.0;
+
+    /* q ? 0 limit */
+    if (fabs(q) < eps) {
+        return N;
+    }
+
+    const double Dm1 = D - 1.0;
+
+    /* Avoid singularity at D = 1 */
+    if (fabs(Dm1) < 1e-8) {
+        /* Limit D ? 1:
+           S(q) = 1 + (N-1)/N * atan(q*r*N)/(q*r)
+        */
+        double x = q * r;
+        double y = x * N;  /* since N^(1/D) ? N when D?1 */
+        return 1.0 + (N-1.0)/N * atan(y)/(x);
+    }
+
+    /* General case */
+    double x = q * r;
+    double Nd = pow(N, 1.0/D);
+    double y = x * Nd;   /* = q*r*N^(1/D) */
+
+    double term =
+        sin(Dm1 * atan(y))
+        * pow(x, -D)
+        * pow(1.0 + 1.0/(y*y), -0.5*Dm1);
+
+    double prefactor = (N - 1.0) / (N * Dm1);
+
+    return 1.0 + prefactor * term;
+}
+
+/* ========== Sasmodels interface ========== */
+double Iq(double q, double radius_effective, double volfraction,
+          double fractal_dim, double n_aggreg)
+{
+    (void)volfraction;
+    return fractal_sq_N(q, radius_effective, fractal_dim, n_aggreg);
+}

sasmodels/models/fractal_aggregate.py

@@ -0,0 +1,100 @@
+# fractal_aggregate.py
+# Sasmodels plugin wrapper for the fractal structure factor
+#
+# Matches the format and syntax of free_rotating_chain.py
+r"""
+Definition
+----------
+This model calculates the structure factor of a fractal-like aggregates.
+
+WARNING: Should not be used in combination with very anisotropic particle shapes.
+
+It uses the following equation:
+
+.. math::
+
+    S(q) = 1 + \frac{D_f\, \Gamma(D_f-1)}{\left[1+1/(q\xi)^2\right]^{(D_f-1)/2}}
+    \frac{\sin\left[(D_f-1)\tan^{-1}(q\xi)\right]}{(q R_0)^{D_f}}
+
+where $\xi$ is the correlation length representing the cluster size and $D_f$
+is the fractal dimension, representing the self similarity of the structure.
+
+The expression has been reformulated so that the aggregation number $N_{agg}$
+is a fit parameter instead of $\xi$.
+
+When combined with a polydisperse form factor, simple $P(q)\,S(q)$ with
+$S(q)$ evaluated at an effective radius differs from averaging
+$P(q,r)\,S(q,r)$ over $r$. For polydisperse primary particles, the
+$\beta$-/decoupling approximation (*structure_factor_mode = 1*) is usually
+more appropriate.
+
+Parameters
+----------
+radius_effective : effective scatterer radius (Å), half the center-to-center distance; use unconstrained *radius_effective_mode* to fit independently of the form-factor radius
+volfraction : required for $P@S$ products; **not used** in this $S(q)$ (see ``sasmodels.product`` until optional structure-factor volfraction exists)
+fractal_dim : fractal dimension
+n_aggreg : number of particles in the fractal cluster
+
+See also
+Larsen, A. H., Pedersen, J. S., & Arleth, L. (2020). Assessment of
+structure factors for analysis of small-angle scattering data from
+desired or undesired aggregates. Applied Crystallography, 53(4), 991-1005.
+
+
+
+Validation
+----------
+
+Translated FORTRAN code
+
+References
+----------
+
+# Teixeira, J. (1988). Small-angle scattering by fractal systems. Applied Crystallography, 21(6), 781-785.
+
+Authorship and Verification
+----------------------------
+
+* **Author:** Jan Skov Pedersen
+* **Last Modified by:** Jan Skov Pedersen, April 12, 2026
+* **Last Reviewed by:** Reviewer Name Here **Date:**
+"""
+
+
+
+import numpy as np
+
+name = "fractal_aggregate"
+title = "Fractal aggregate structure factor"
+description = """
+Fractal structure factor S(q)
+"""
+
+category = "structure-factor"
+structure_factor = True
+
+# Must match C kernel: Iq(q, radius_effective, volfraction, fractal_dim, n_aggreg)
+parameters = [
+    ["radius_effective", "Ang", 10.0, [0.0, np.inf], "",
+     "effective scatterer radius (half center-to-center distance)"],
+    ["volfraction", "", 0.2, [0.0, 1.0], "",
+     "unused in S(q); required as second parameter for P@S products"],
+    ["fractal_dim", "", 2.0, [1.0, 3.0], "", "Fractal dimension"],
+    ["n_aggreg", "", 50.0, [1.0, np.inf], "", "Number of particles"],
+]
+
+# Kernel source
+source = ["fractal_aggregate.c"]
+
+def random():
+    import random
+
+    return {
+        "radius_effective": random.uniform(2.5, 25.0),
+        "volfraction": random.uniform(0.01, 0.3),
+        "fractal_dim": random.uniform(1.1, 2.9),
+        "n_aggreg": random.uniform(5.0, 300.0),
+    }
+
+def test():
+    print("fractal_aggregate plugin loaded correctly.")