Skip to content

Commit

Permalink
refactor model
Browse files Browse the repository at this point in the history
  • Loading branch information
@dehoni
dehoni committed Sep 10, 2024
1 parent f9c7c5f commit 5414c55
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 78 deletions.
114 changes: 46 additions & 68 deletions sasmodels/models/micromagnetic_FF_3D.c
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
//Core-shell form factor for anisotropy field (Hkx, Hky and Hkz), Nuc and lonitudinal magnetization Mz
static double
form_volume(double nuc_radius, double nuc_thickness, double mag_radius, double mag_thickness, double hk_radius)
form_volume(double radius, double thickness)
{
if( nuc_radius+nuc_thickness>mag_radius+ mag_thickness)
return M_4PI_3 * cube(nuc_radius + nuc_thickness);
if( radius+thickness>radius+ thickness)
return M_4PI_3 * cube(radius + thickness);
else
return M_4PI_3 * cube(mag_radius + mag_thickness);
return M_4PI_3 * cube(radius + thickness);

}

Expand All @@ -27,10 +27,9 @@ static double fq(double q, double radius,
static double reduced_field(double q, double Ms, double Hi,
double A)
{
if( Hi > 1.0e-6 )
return Ms / (Hi + 2.0 * A * 4.0 * M_PI / Ms * q * q * 10.0);//q in 10e10 m-1, A in 10e-12 J/m, mu0 in 1e-7
else
return Ms / (1.0e-6 + 2.0 * A * 4.0 * M_PI / Ms * q * q * 10.0);
// q in 10e10 m-1, A in 10e-12 J/m, mu0 in 1e-7
return Ms / (fmax(Hi, 1.0e-6) + 2.0 * A * 4.0 * M_PI / Ms * q * q * 10.0);

}

static double DMI_length(double Ms, double D, double qval)
Expand All @@ -51,7 +50,9 @@ static double fqMxreal( double qx, double qy, double qz, double Mz, double Hkx,
const double qsq = qx*qx + qy*qy + qz*qz;
const double q = sqrt(qsq);
const double Hr = reduced_field(q, Ms, Hi, A);
const double f = Hr*(Hkx*(qsq+Hr*qy*qy)-Ms*Mz*qx*qz*(1.0+Hr*DMI_length(Ms, D,q)*DMI_length(Ms, D,q))-Hky*Hr*qx*qy)/(qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI_length(Ms, D,qz)*q));
const double DMI = DMI_length(Ms, D,q);
const double denominator = (qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI*q))/Hr;
const double f = (Hkx*(qsq+Hr*qy*qy)-Ms*Mz*qx*qz*(1.0+Hr*square(DMI))-Hky*Hr*qx*qy)/denominator;
return f;
}

Expand All @@ -60,7 +61,9 @@ static double fqMximag(double qx, double qy, double qz, double Mz, double Hkx, d
const double qsq = qx*qx + qy*qy + qz*qz;
const double q = sqrt(qsq);
double Hr = reduced_field(q, Ms, Hi, A);
const double f = -Hr*qsq*(Ms*Mz*(1.0+Hr)*DMI_length(Ms, D,qy)+Hky*Hr*DMI_length(Ms, D,qz))/(qsq+Hr*(qx*qx+qy*qy) -square(Hr*DMI_length(Ms, D,qz)*q));
const double DMI = DMI_length(Ms, D,q);
const double denominator = (qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI*q))/Hr;
const double f = -qsq*(Ms*Mz*(1.0+Hr)*DMI+Hky*Hr*DMI)/denominator;
return f;
}

Expand All @@ -69,7 +72,9 @@ static double fqMyreal( double qx, double qy, double qz, double Mz, double Hkx,
const double qsq = qx*qx + qy*qy + qz*qz;
const double q = sqrt(qsq);
const double Hr = reduced_field(q, Ms, Hi, A);
const double f = Hr*(Hky*(qsq+Hr*qx*qx)-Ms*Mz*qy*qz*(1.0+Hr*DMI_length(Ms, D,q)*DMI_length(Ms, D,q))-Hkx*Hr*qx*qy)/(qsq+Hr*(qx*qx+qy*qy) -square(Hr*DMI_length(Ms, D,qz)*q));
const double DMI = DMI_length(Ms, D,q);
const double denominator = (qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI*q))/Hr;
const double f = (Hky*(qsq+Hr*qx*qx)-Ms*Mz*qy*qz*(1.0+Hr*square(DMI))-Hkx*Hr*qx*qy)/denominator;
return f;
}

Expand All @@ -78,27 +83,24 @@ static double fqMyimag( double qx, double qy, double qz, double Mz, double Hkx,
const double qsq = qx*qx + qy*qy + qz*qz;
const double q = sqrt(qsq);
const double Hr = reduced_field(q, Ms, Hi, A);
const double f = Hr*qsq*(Ms*Mz*(1.0+Hr)*DMI_length(Ms, D,qx)-Hkx*Hr*DMI_length(Ms, D,qz))/(qsq+Hr*(qx*qx+qy*qy) -square(Hr*DMI_length(Ms, D,qz)*q));
const double DMI = DMI_length(Ms, D,q);
const double denominator = (qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI*q))/Hr;
const double f = qsq*(Ms*Mz*(1.0+Hr)*DMI-Hkx*Hr*DMI)/denominator;
return f;
}

static double
Iqxy(double qx, double qy, double nuc_radius, double nuc_thickness, double mag_radius, double mag_thickness, double hk_radius, double nuc_sld_core, double nuc_sld_shell, double nuc_sld_solvent, double mag_sld_core, double mag_sld_shell, double mag_sld_solvent, double hk_sld_core, double Hi, double Ms, double A, double D, double up_i, double up_f, double alpha, double beta)
Calculate_Scattering(double q, double cos_theta, double sin_theta, double radius, double thickness, double nuc_sld_core, double nuc_sld_shell, double nuc_sld_solvent, double mag_sld_core, double mag_sld_shell, double mag_sld_solvent, double hk_sld_core, double Hi, double Ms, double A, double D, double up_i, double up_f, double alpha, double beta)
{
const double q = sqrt(qx*qx + qy*qy);
if (q > 1.0e-16 ) {
const double cos_theta=qx/q;
const double sin_theta=qy/q;

double qrot[3];
set_scatvec(qrot,q,cos_theta, sin_theta, alpha, beta);
// 0=dd.real, 1=dd.imag, 2=uu.real, 3=uu.imag, 4=du.real, 6=du.imag, 7=ud.real, 5=ud.imag
double weights[8];
set_weights(up_i, up_f, weights);

double mz=fq(q, mag_radius, mag_thickness, mag_sld_core, mag_sld_shell, mag_sld_solvent);
double nuc=fq(q, nuc_radius, nuc_thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent);

double mz=fq(q, radius, thickness, mag_sld_core, mag_sld_shell, mag_sld_solvent);
double nuc=fq(q, radius, thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent);
double hk=fq(q, radius, 0, hk_sld_core, 0, 0);

double cos_gamma, sin_gamma;
double sld[8];
Expand All @@ -111,9 +113,9 @@ Iqxy(double qx, double qy, double nuc_radius, double nuc_thickness, double mag_r
//Only the core of the defect/particle in the matrix has an effective
//anisotropy (for simplicity), for the effect of different, more complex
// spatial profile of the anisotropy see Michels PRB 82, 024433 (2010)
double Hkx= fq(q, hk_radius, 0, hk_sld_core, 0, 0)*sin_gamma;
double Hky= fq(q, hk_radius, 0, hk_sld_core, 0, 0)*cos_gamma;
double Hkx= hk*sin_gamma;
double Hky= hk*cos_gamma;

double mxreal=fqMxreal(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
double mximag=fqMximag(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
double myreal=fqMyreal(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
Expand All @@ -131,14 +133,32 @@ Iqxy(double qx, double qy, double nuc_radius, double nuc_thickness, double mag_r
}
}
total_F2 += GAUSS_W[i] * form ;
return total_F2;
}
}


static double
Iqxy(double qx, double qy, double radius, double thickness, double nuc_sld_core, double nuc_sld_shell, double nuc_sld_solvent, double mag_sld_core, double mag_sld_shell, double mag_sld_solvent, double hk_sld_core, double Hi, double Ms, double A, double D, double up_i, double up_f, double alpha, double beta)
{
const double q = sqrt(qx*qx + qy*qy);
if (q > 1.0e-16 ) {
const double cos_theta=qx/q;
const double sin_theta=qy/q;
} else {
const double cos_theta=0.0;
const double sin_theta=0.0;

double total_F2 = Calculate_Scattering(q, cos_theta, sin_theta, radius, thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent, mag_sld_core, mag_sld_shell, mag_sld_solvent, hk_sld_core, Hi, Ms, A, D, up_i, up_f, alpha, beta);

return 0.5*1.0e-4*total_F2;
}
}



static double
Iq(double q, double nuc_radius, double nuc_thickness, double mag_radius, double mag_thickness, double hk_radius,double nuc_sld_core, double nuc_sld_shell, double nuc_sld_solvent, double mag_sld_core, double mag_sld_shell, double mag_sld_solvent, double hk_sld_core, double Hi, double Ms, double A, double D, double up_i, double up_f, double alpha, double beta)
Iq(double q, double radius, double thickness, double nuc_sld_core, double nuc_sld_shell, double nuc_sld_solvent, double mag_sld_core, double mag_sld_shell, double mag_sld_solvent, double hk_sld_core, double Hi, double Ms, double A, double D, double up_i, double up_f, double alpha, double beta)
{
// slots to hold sincos function output of the orientation on the detector plane
double sin_theta, cos_theta;
Expand All @@ -148,49 +168,7 @@ Iq(double q, double nuc_radius, double nuc_thickness, double mag_radius, double
const double theta = M_PI * (GAUSS_Z[j] + 1.0); // 0 .. 2 pi
SINCOS(theta, sin_theta, cos_theta);

double qrot[3];
set_scatvec(qrot,q,cos_theta, sin_theta, alpha, beta);
// 0=dd.real, 1=dd.imag, 2=uu.real, 3=uu.imag, 4=du.real, 5=du.imag, 6=ud.real, 7=ud.imag
double weights[8];
set_weights(up_i, up_f, weights);

double mz=fq(q, mag_radius, mag_thickness, mag_sld_core, mag_sld_shell, mag_sld_solvent);
double nuc=fq(q, nuc_radius, nuc_thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent);


double cos_gamma, sin_gamma;
double sld[8];

//loop over random anisotropy axis with isotropic orientation gamma for Hkx and Hky
//To be modified for textured material see also Weissmueller et al. PRB 63, 214414 (2001)
double total_F2 = 0.0;
for (int i=0; i<GAUSS_N ;i++) {
const double gamma = M_PI * (GAUSS_Z[i] + 1.0); // 0 .. 2 pi
SINCOS(gamma, sin_gamma, cos_gamma);
//Only the core of the defect/particle in the matrix has an effective
//anisotropy (for simplicity), for the effect of different, more complex
//spatial profile of the anisotropy see Michels PRB 82, 024433 (2010).
double Hkx= fq(q, hk_radius, 0, hk_sld_core, 0, 0)*sin_gamma;
double Hky= fq(q, hk_radius, 0, hk_sld_core, 0, 0)*cos_gamma;

double mxreal=fqMxreal(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
double mximag=fqMximag(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
double myreal=fqMyreal(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
double myimag=fqMyimag(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);

mag_sld(qrot[0], qrot[1], qrot[2], mxreal, mximag, myreal, myimag, mz, 0, nuc, sld);
double form = 0.0;
for (unsigned int xs=0; xs<8; xs++) {
if (weights[xs] > 1.0e-8 ) {
// Since the cross section weight is significant, set the slds
// to the effective slds for this cross section, call the
// kernel, and add according to weight.
// loop over uu, ud real, du real, dd, ud imag, du imag
form += weights[xs]*sld[xs]*sld[xs];
}
}
total_F2 += GAUSS_W[i] * form ;
}
double total_F2 = Calculate_Scattering(q, cos_theta, sin_theta, radius, thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent, mag_sld_core, mag_sld_shell, mag_sld_solvent, hk_sld_core, Hi, Ms, A, D, up_i, up_f, alpha, beta);
total_F1D += GAUSS_W[j] * total_F2 ;
}
//convert from [1e-12 A-1] to [cm-1]
Expand Down
14 changes: 4 additions & 10 deletions sasmodels/models/micromagnetic_FF_3D.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,11 +136,8 @@

# pylint: disable=bad-whitespace, line-too-long
# ["name", "units", default, [lower, upper], "type","description"],
parameters = [["nuc_radius", "Ang", 50., [0, inf], "volume", "Structural radius of the core"],
["nuc_thickness", "Ang", 40., [0, inf], "volume", "Structural thickness of shell"],
["mag_radius", "Ang", 50., [0, inf], "volume", "Magnetic radius of the core"],
["mag_thickness", "Ang", 40., [0, inf], "volume", "Magnetic thickness of shell"],
["hk_radius", "Ang", 50., [0, inf], "volume", "Anisotropy radius of the core"],
parameters = [["radius", "Ang", 50., [0, inf], "volume", "Structural radius of the core"],
["thickness", "Ang", 40., [0, inf], "volume", "Structural thickness of shell"],
["nuc_sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", "Core scattering length density"],
["nuc_sld_shell", "1e-6/Ang^2", 1.7, [-inf, inf], "", "Scattering length density of shell"],
["nuc_sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", "Solvent scattering length density"],
Expand Down Expand Up @@ -177,11 +174,8 @@ def random():
radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
thickness = outer_radius - radius
pars = dict(
nuc_radius=radius,
nuc_thickness=thickness,
mag_radius=radius,
mag_thickness=thickness,
hk_radius=radius
radius=radius,
thickness=thickness,

)
return pars
Expand Down

0 comments on commit 5414c55

Please sign in to comment.