-
Notifications
You must be signed in to change notification settings - Fork 27
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
parallelepiped chains based on MagneticOrientedChains
have inserted the option for chains (taken form MagneticOrientedCHains) into the parallelepiped script from sasview
- Loading branch information
1 parent
59f5bae
commit eabe6f3
Showing
2 changed files
with
310 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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; | ||
|
||
} | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 | ||
|
||
|
||
|