From eabe6f33d8be10e726cf6499d7915c28a348051c Mon Sep 17 00:00:00 2001
From: astellhorn <93910032+astellhorn@users.noreply.github.com>
Date: Fri, 19 Jan 2024 17:07:31 +0100
Subject: [PATCH] parallelepiped chains based on MagneticOrientedChains

have inserted the option for chains (taken form MagneticOrientedCHains) into the parallelepiped script from sasview
---
 sasmodels/models/ParallelepipedChain.c  | 218 ++++++++++++++++++++++++
 sasmodels/models/ParallelepipedChain.py |  92 ++++++++++
 2 files changed, 310 insertions(+)
 create mode 100644 sasmodels/models/ParallelepipedChain.c
 create mode 100644 sasmodels/models/ParallelepipedChain.py

diff --git a/sasmodels/models/ParallelepipedChain.c b/sasmodels/models/ParallelepipedChain.c
new file mode 100644
index 00000000..4676a024
--- /dev/null
+++ b/sasmodels/models/ParallelepipedChain.c
@@ -0,0 +1,218 @@
+static double
+form_volume(double length_a, double length_b, double length_c)
+{
+    return length_a * length_b * length_c;
+}
+
+static double
+radius_from_excluded_volume(double length_a, double length_b, double length_c)
+{
+    double r_equiv, length;
+    double lengths[3] = {length_a, length_b, length_c};
+    double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2]));
+    double length_1 = lengthmax;
+    double length_2 = lengthmax;
+    double length_3 = lengthmax;
+
+    for(int ilen=0; ilen<3; ilen++) {
+        if (lengths[ilen] < length_1) {
+            length_2 = length_1;
+            length_1 = lengths[ilen];
+            } else {
+                if (lengths[ilen] < length_2) {
+                        length_2 = lengths[ilen];
+                }
+            }
+    }
+    if(length_2-length_1 > length_3-length_2) {
+        r_equiv = sqrt(length_2*length_3/M_PI);
+        length  = length_1;
+    } else  {
+        r_equiv = sqrt(length_1*length_2/M_PI);
+        length  = length_3;
+    }
+
+    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length)));
+}
+
+static double
+radius_effective(int mode, double length_a, double length_b, double length_c)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent cylinder excluded volume
+        return radius_from_excluded_volume(length_a,length_b,length_c);
+    case 2: // equivalent volume sphere
+        return cbrt(length_a*length_b*length_c/M_4PI_3);
+    case 3: // half length_a
+        return 0.5 * length_a;
+    case 4: // half length_b
+        return 0.5 * length_b;
+    case 5: // half length_c
+        return 0.5 * length_c;
+    case 6: // equivalent circular cross-section
+        return sqrt(length_a*length_b/M_PI);
+    case 7: // half ab diagonal
+        return 0.5*sqrt(length_a*length_a + length_b*length_b);
+    case 8: // half diagonal
+        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c);
+    }
+}
+
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
+    double sld,
+    double solvent_sld,
+    double length_a,
+    double length_b,
+    double length_c)
+{
+    const double mu = 0.5 * q * length_b;
+
+    // Scale sides by B
+    const double a_scaled = length_a / length_b;
+    const double c_scaled = length_c / length_b;
+
+    // outer integral (with gauss points), integration limits = 0, 1
+    double outer_total_F1 = 0.0; //initialize integral
+    double outer_total_F2 = 0.0; //initialize integral
+    for( int i=0; i<GAUSS_N; i++) {
+        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
+        const double mu_proj = mu * sqrt(1.0-sigma*sigma);
+
+        // inner integral (with gauss points), integration limits = 0, 1
+        // corresponding to angles from 0 to pi/2.
+        double inner_total_F1 = 0.0;
+        double inner_total_F2 = 0.0;
+        for(int j=0; j<GAUSS_N; j++) {
+            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
+            double sin_uu, cos_uu;
+            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
+            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
+            const double si2 = sas_sinx_x(mu_proj * cos_uu);
+            const double fq = si1 * si2;
+            inner_total_F1 += GAUSS_W[j] * fq;
+            inner_total_F2 += GAUSS_W[j] * fq * fq;
+        }
+        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
+        inner_total_F1 *= 0.5;
+        inner_total_F2 *= 0.5;
+
+        const double si = sas_sinx_x(mu * c_scaled * sigma);
+        outer_total_F1 += GAUSS_W[i] * inner_total_F1 * si;
+        outer_total_F2 += GAUSS_W[i] * inner_total_F2 * si * si;
+    }
+    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
+    outer_total_F1 *= 0.5;
+    outer_total_F2 *= 0.5;
+
+    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
+    const double V = form_volume(length_a, length_b, length_c);
+    const double contrast = (sld-solvent_sld);
+    const double s = contrast * V;
+    *F1 = 1.0e-2 * s * outer_total_F1;
+    *F2 = 1.0e-4 * s * s * outer_total_F2;
+}
+
+static double
+Iqabc(double qa, double qb, double qc,
+    double sld,
+    double solvent_sld,
+    double mag_sld,
+    double length_a,
+    double length_b,
+    double length_c,
+    double Length, double Singlets, double Doubles, double Trimers, double Quadramers, double Pentamers)
+{
+    const double siA = sas_sinx_x(0.5*length_a*qa);
+    const double siB = sas_sinx_x(0.5*length_b*qb);
+    const double siC = sas_sinx_x(0.5*length_c*qc);
+    const double V = form_volume(length_a, length_b, length_c);
+    const double drho = (sld - solvent_sld);
+    const double mrho = mag_sld
+    const double Amp = V * drho * siA * siB * siC;
+    const double MAmp = V * mrho * siA * siB * siC;
+
+    // adding several particles in the chain:
+    double VolumeFraction = 1.0;
+    double SingletFraction = Singlets;
+    double DimerFraction = Doubles;
+    double TrimerFraction = Trimers;
+    double QuadramerFraction = Quadramers;
+    double PentamerFraction = Pentamers;
+
+    // leave out polydispersity and angular orientation
+
+    // Just start parameters:
+    double SingletIntensity = 0;
+    double MSingletIntensity = 0;
+    double DimerIntensity = 0;
+    double MDimerIntensity = 0;
+    double TrimerIntensity = 0;
+    double MTrimerIntensity = 0;
+    double QuadramerIntensity = 0;
+    double MQuadramerIntensity = 0;
+    double PentamerIntensity = 0;
+    double MPentamerIntensity = 0;
+
+    // we take out: (i) polydispersity of chains (called "anglewt/sigma" before), (ii) orientation of chain (called "Viewing angle")
+    // Now: chain only oriented along Qx
+    // What does it mean for the magnetic orientation of the particles?: at the moment not oriented / random
+
+    double Q_X = q*cos(0);
+    double Q_Y = q*sin(0);
+
+
+    double ChainProjX = cos(0);
+    double ChainProjY = cos(phi);
+    double ChainProjZ = sin(phi);
+
+
+    double real_phase = 1.0;
+    double img_phase = 0.0;
+    double mreal_phase = 1.0;
+    double mimg_phase = 0.0;
+
+    for(int k=1; k<5; k++){
+        real_phase += cos(k*Length*(Q_X*ChainProjX + Q_Y*ChainProjY));
+        img_phase += sin(k*Length*(Q_X*ChainProjX + Q_Y*ChainProjY));
+
+    if(k==1){
+    DimerIntensity  += pow(Amp*img_phase,2))/(2.0*Vol);
+    MDimerIntensity  += pow(MAmp*mimg_phase,2))/(2.0*Vol);
+    }
+    if(k==2){
+    TrimerIntensity  += pow(Amp*img_phase,2))/(3.0*Vol);
+    MTrimerIntensity  += pow(MAmp*mimg_phase,2))/(3.0*Vol);
+    }
+    if(k==3){
+    QuadramerIntensity  += pow(Amp*img_phase,2))/(4.0*Vol);
+    MQuadramerIntensity  += pow(MAmp*mimg_phase,2))/(4.0*Vol);
+    }
+    if(k==4){
+    PentamerIntensity  += pow(Amp*img_phase,2))/(5.0*Vol);
+    MPentamerIntensity  += pow(MAmp*mimg_phase,2))/(5.0*Vol);
+    }
+    }
+    //end k loop for dimers
+
+    double FractionScale = SingletFraction + DimerFraction + TrimerFraction + QuadramerFraction + PentamerFraction;
+    if(FractionScale == 0){FractionScale = 1.0;}
+
+    double SIntensity = SingletFraction*SingletIntensity + DimerFraction*DimerIntensity + TrimerFraction*TrimerIntensity + QuadramerFraction*QuadramerIntensity + PentamerFraction*PentamerIntensity;
+    double MIntensity = 0.0;
+    if(MVar <= 1){
+        MIntensity = MSingletIntensity*(SingletFraction + DimerFraction + TrimerFraction + QuadramerFraction + PentamerFraction);
+    }
+    else{
+        MIntensity = SingletFraction*MSingletIntensity + DimerFraction*MDimerIntensity + TrimerFraction*MTrimerIntensity + QuadramerFraction*MQuadramerIntensity + PentamerFraction*MPentamerIntensity;
+    }
+
+    double Intensity = (SIntensity+MIntensity)*(1E4)/FractionScale;
+
+    return Intensity;
+
+}
+
diff --git a/sasmodels/models/ParallelepipedChain.py b/sasmodels/models/ParallelepipedChain.py
new file mode 100644
index 00000000..48c07b47
--- /dev/null
+++ b/sasmodels/models/ParallelepipedChain.py
@@ -0,0 +1,92 @@
+r"""
+
+Definition
+----------
+This plug-in model calculates chains of oriented parallelepipeds, with the option of 
+adding a magnetic SLD. The chain scattering is the incoherent 
+sum of a user-defined combination of singletons, dimers, trimers, 
+quadramers, and pentamers. Note that no matter the numerical values 
+selected for the amount of each chain type, the fraction of each will be 
+normalized such that the sum of the chain type fractions is unity. 
+
+The chains are oriented along the x-direction. 
+
+The magnetism of the chains is given by a uniform magnetic sld within each parallelepiped particle.
+
+The scattering amplitude form factor is calculated in same way as the parallelepiped (s. https://www.sasview.org/docs/user/models/parallelepiped.html), and it is then multiplied by a complex structure factor that depends on chain length:
+
+.. math::
+
+    I(q)= scale*V*((nuc_sld - sld_solvent)+mag_sld)^2*P(q,alpha)*(  ext{real  phase}^2 + 	ext{img  phase}^2)  +background
+        P(q,alpha) = integral from 0 to 1 of ...
+           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
+        with
+            phi(mu,a) = integral from 0 to 1 of ..
+            (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du
+            S(x) = sin(x)/x
+            mu = q*B
+        V: Volume of the rectangular parallelepiped
+        alpha: angle between the long axis of the
+            parallelepiped and the q-vector for 1D
+
+and
+.. math::
+    	ext{real  phase} =  1.0 + sum_{k=0}^{4} sum_{n=0}^{k} cos(k*Length*(Q_X*ChainProjX))
+and
+.. math::
+    	ext{img  phase} =  0.0 + sum_{k=0}^{4} sum_{n=0}^{k} sin(k*Length*(Q_X*ChainProjX ))
+
+
+Here $V$ is the volume of one parallelepiped, $nuc_sld$ is the univorm nuclear sld of one parallelepiped, $sld_solvent$ the surrounding solvent sld, $sld_mag$ the uniform magnetic sld, $k$ defines the number of parallelepipeds in the chain which is oriented along Q_X. This model is a combination of a single parallelepiped (s. https://www.sasview.org/docs/user/models/parallelepiped.html) together with the chain model described in https://marketplace.sasview.org/models/143/, but with fixed oriented direction along x, without polydisperse orientation of the chain, and with a uniform magnetic sld.
+"""
+import numpy as np
+from sasmodels.special import *
+from numpy import inf
+
+name = "ParallelepipedChain"
+title = "Base-script: Rectangular parallelepiped with uniform scattering length density. Add-on: Chain of parallelepipeds along x-direction, with uniform magnetic scattering length density (no applied magnetic field)"
+description = """User model for chains of parallelepipeds oriented along X-axis"""
+
+category = "shape:parallelepiped"
+
+parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld",
+               "Parallelepiped scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld",
+               "Solvent scattering length density"],
+              ["mag_sld", "1e-6/Ang^2", 1, [-inf, inf], "sld",
+               "magnetic scattering length density"],
+              ["length_a", "Ang", 35, [0, inf], "length",
+               "Shorter side of the parallelepiped"],
+              ["length_b", "Ang", 75, [0, inf], "length",
+               "Second side of the parallelepiped"],
+              ["length_c", "Ang", 400, [0, inf], "length",
+               "Larger side of the parallelepiped"],
+              ["Length", "Ang", 400, [0, inf], "length",
+               "Length of the distance between parallelepipeds"],
+              ['singlets', 'Fraction of singlets', 1.0, [0, 100], '', ''],
+              ['doublets', 'Fraction of doubles', 1.0, [0, 100], '', ''],
+              ['trimers', 'Fraction of trimers', 1.0, [0, 100], '', ''],
+              ['quadramers', 'Fraction of quadramers', 1.0, [0, 100], '', ''],
+              ['pentamers', 'Fraction of pentamers', 1.0, [0, 100], '', ''],
+             ]
+
+source = ["lib/gauss76.c", "ParallelepipedChain.c"]
+have_Fq = True
+radius_effective_modes = [
+    "equivalent cylinder excluded volume", "equivalent volume sphere",
+    "half length_a", "half length_b", "half length_c",
+    "equivalent circular cross-section", "half ab diagonal", "half diagonal",
+    ]
+
+def random():
+    """Return a random parameter set for the model."""
+    length = 10**np.random.uniform(1, 4.7, size=3)
+    pars = dict(
+        length_a=length[0],
+        length_b=length[1],
+        length_c=length[2],
+    )
+    return pars
+
+
+