Magnetic Whirl - meron.c
//cylinder cross section for nuclear scattering
static double
form_volume(double length, double radius)
{
return M_PI*radius*radius*length;
}
static double
radius_from_excluded_volume( double length,double radius)
{
return 0.5*cbrt(0.75*radius*(2.0*radius*length
+ (radius + length)*(M_PI*radius + length)));
}
static double
radius_from_volume(double length, double radius)
{
return cbrt(M_PI*radius*radius*length/M_4PI_3);
}
static double
radius_from_diagonal(double length, double radius)
{
return sqrt(radius*radius + 0.25*length*length);
}
static double
radius_effective(int mode, double length, double radius)
{
switch (mode) {
default:
case 1:
return radius_from_excluded_volume(length, radius );
case 2:
return radius_from_volume(length, radius);
case 3:
return radius;
case 4:
return 0.5*length;
case 5:
return (radius < 0.5*length ? radius : 0.5*length);
case 6:
return (radius > 0.5*length ? radius : 0.5*length);
case 7:
return radius_from_diagonal(length, radius);
}
}
static double
fq(double x, double y, double z, double length, double radius, double sld, double solvent_sld)
{
const double q=MAG_VEC(x,y,0);
const double s = (sld - solvent_sld);
return s *sas_2J1x_x(q*radius) * sas_sinx_x(z*0.5*length);
}
//Definition of special integral function on Bessel functions see
//Supplemental Material of Metlov and Michels Sci. Rep. 6, 25055 (2016)
static double total_F2(double k)
{
return 0.5 * k * M_PI * (sas_J1(k)* struve(0,k)- sas_J0(k)*struve(1,k));
}
static double muperp0( double k)
{
return total_F2(k)/k/k;
}
static double muz0(double cos_alpha,double sin_alpha, double k)
{
return muperp0(k)*sin_alpha; //complex quantity does not mix with first order terms
}
static double muy0(double cos_alpha,double sin_alpha, double k)
{
return -muperp0(k)*cos_alpha; //complex quantity does not mix with first order terms
}
// Mx is component along the length of the cylinder, i.e. perpendicular to the vortex
//structure in the disc.
static double fqMxreal( double x, double y, double z, double L, double R, double Ms)
{
const double q=MAG_VEC(x, y,0);
const double k=R*q;
return Ms*sas_2J1x_x(k)* sas_sinx_x(z*0.5*L);
}
static double fqMyimag( double x, double y, double z, double L, double R, double Ms)
{
const double q=MAG_VEC(x, y,0);
const double cos_alpha=x/q;
const double sin_alpha=y/q;
const double k=R*q;
const double fx=Ms*sas_sinx_x(0.5*L*z);
return M_PI*fx*muy0(cos_alpha, sin_alpha, k);
}
static double fqMzimag( double x, double y, double z, double L, double R, double Ms)
{
const double q=MAG_VEC(x, y,0);
const double cos_alpha=x/q;
const double sin_alpha=y/q;
const double k=R*q;
const double fx=Ms*sas_sinx_x(0.5*L*z);
return M_PI*fx*muz0(cos_alpha, sin_alpha, k);
}
//calculate 2D from _fq
static double
Iqxy(double qx, double qy, double length, double radius, double core_nuc, double solvent_nuc, double magnetic_sld_disc, double magnetic_sld_length,double up_i, double up_f, double alpha, double beta)
{
const double q = MAG_VEC(qx,qy,0 );
if (q > 1.0e-16 ) {
const double cos_theta=qx/q;
const double sin_theta=qy/q;
double qrot[3];
double sld[8];
double weights[8];
set_scatvec(qrot,q,cos_theta, sin_theta, alpha, beta);
// double nuc=fq(qrot[0], qrot[1], qrot[2],length, radius, core_nuc, solvent_nuc);
double nuc=fqMxreal(qrot[0], qrot[1], qrot[2],length, radius, core_nuc);
// 0=dd.real, 1=dd.imag, 2=uu.real, 3=uu.imag, 4=du.real, 6=du.imag, 7=ud.real, 5=ud.imag
set_weights(up_i, up_f, weights);
double mxreal=fqMxreal(qrot[0], qrot[1], qrot[2], length, radius, magnetic_sld_length);
double myimag=fqMyimag(qrot[0], qrot[1], qrot[2], length, radius, magnetic_sld_disc);
double mzimag=fqMzimag(qrot[0], qrot[1], qrot[2], length, radius, magnetic_sld_disc);
mag_sld(qrot[0], qrot[1], qrot[2], 0, mzimag, 0, myimag, mxreal, 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];
}
}
return 1.0e-4*form* square(form_volume(length, radius));
}
}
static double
Iq(double q,double length, double radius, double core_nuc, double solvent_nuc, double magnetic_sld_disc, double magnetic_sld_length, double up_i, double up_f, double alpha, double beta)
{
double sin_theta, cos_theta; // slots to hold sincos function output of the orientation on the detector plane
double total_F1D = 0.0;
for (int j=0; j<GAUSS_N ;j++) {
const double theta = M_PI * (GAUSS_Z[j] + 1.0); // 0 .. 2 pi
SINCOS(theta, sin_theta, cos_theta);
double sld[8];
double qrot[3];
double weights[8];
// 0=dd.real, 1=dd.imag, 2=uu.real, 3=uu.imag, 4=du.real, 5=du.imag, 6=ud.real, 7=ud.imag
set_weights(up_i, up_f, weights);
set_scatvec(qrot,q,cos_theta, sin_theta, alpha, beta);
// double nuc=fq(qrot[0], qrot[1], qrot[2], length, radius, core_nuc, solvent_nuc);
double nuc=fqMxreal(qrot[0], qrot[1], qrot[2],length, radius, core_nuc);
double mxreal=fqMxreal(qrot[0], qrot[1], qrot[2], length, radius, magnetic_sld_length);
double myimag=fqMyimag(qrot[0], qrot[1], qrot[2], length, radius, magnetic_sld_disc);
double mzimag=fqMzimag(qrot[0], qrot[1], qrot[2], length, radius, magnetic_sld_disc);
mag_sld(qrot[0], qrot[1], qrot[2], 0, mzimag, 0, myimag, mxreal, 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_F1D += GAUSS_W[j] * form ;
}
//convert from [1e-12 A-1] to [cm-1]
return 0.5*1.0e-4*total_F1D* square(form_volume(length, radius ));
}
}
Back to Model
Download