Superparamagnetic Core-Shell Spheres with 3D field orientation - magnetic_langevin_core_shell_3D.c

```    ```  static double clipp(double value, double low, double high) //from kernel_iq.c
{
return (value < low ? low : (value > high ? high : value));
}

static double langevin(
double x) {
// cotanh(x)-1\x

if (x < 0.00001) {
// avoid dividing by zero
return 1.0/3.0*x;
} else {
return 1.0/tanh(x)-1/x;
}
}

static double langevinoverx(
double x) {
// cotanh(x)-1\x

if (x < 0.00001) {
// avoid dividing by zero
return 1.0/3.0;
} else {
return langevin(x)/x;
}
}

//from coremultishell.c

static double
{
int n = (int)(fp_n+0.5);
for (int i=0; i < n; i++) {
r += thickness[i];
}
return r;
}

static double
form_volume(double radius, double fp_n, double thickness[])
{
}

{
switch (mode) {
default:
}
}

static double fqnucsq(double q, double sld_core, double radius,
double sld_solvent, double fp_n, double sld[], double thickness[])
{
const int n = (int)(fp_n+0.5);
double f, r, last_sld;
last_sld = sld_core;
f = 0.;
for (int i=0; i<n; i++) {
f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sas_3j1x_x(q*r);
last_sld = sld[i];
r += thickness[i];
}
f += M_4PI_3 * cube(r) * (sld_solvent - last_sld) * sas_3j1x_x(q*r);
return 1e-4 *f*f;
}

static double fqMtranssq(double q, double sld_core, double radius,
double sld_solvent, double fp_n, double sld[], double thickness[],double eta_core,double eta_solvent,double eta[], double delta_solvent,double delta[])
{
const int n = (int)(fp_n+0.5);
double f,fsq_perp, r, last_sld;

last_sld = sld_core*sqrt(langevinoverx(eta_core));
f = 0.;
fsq_perp=0;
for (int i=0; i<n; i++) {
f += M_4PI_3 * cube(r) *(sld[i]*sqrt(langevinoverx(eta[i])*delta[i]) - last_sld) * sas_3j1x_x(q*r);//transversal magnetisation component coaligned with core transversal magnetisation
fsq_perp += square(M_4PI_3 * cube(r) *sld[i]*sqrt(langevinoverx(eta[i])*(1-delta[i])) * sas_3j1x_x(q*r));//transversal magnetisation component perpendicular with core transversal magnetisation, adding incoherently (only of squared amplitude)
last_sld = sld[i]*sqrt(langevinoverx(eta[i])*delta[i]);
r += thickness[i];
}
f += M_4PI_3 * cube(r) * (sld_solvent*sqrt(langevinoverx(eta_solvent)*delta_solvent) - last_sld) * sas_3j1x_x(q*r);//transversal magnetisation component coaligned with core magnetisation
fsq_perp += square(M_4PI_3 * cube(r) * sld_solvent*sqrt(langevinoverx(eta_solvent)*(1-delta_solvent))  * sas_3j1x_x(q*r));
return 1e-4 *(f*f +fsq_perp);
}

static double fqMzsq(double q, double sld_core, double radius,
double sld_solvent, double fp_n, double sld[], double thickness[],double eta_core,double eta_solvent,double eta[])
{
const int n = (int)(fp_n+0.5);
double f, r, last_sld;
last_sld = sld_core*sqrt(1-2*langevinoverx(eta_core));// sqrt() needed for correct scale to intensity later
f = 0.;
for (int i=0; i<n; i++) {
f += M_4PI_3 * cube(r) * (sld[i]*sqrt(1-2*langevinoverx(eta[i])) - last_sld) * sas_3j1x_x(q*r);
last_sld = sld[i]*sqrt(1-2*langevinoverx(eta[i]));
r += thickness[i];
}
f += M_4PI_3 * cube(r) * (sld_solvent*sqrt(1-2*langevinoverx(eta_solvent)) - last_sld) * sas_3j1x_x(q*r);
return 1e-4 *f*f;
}

static double fqnuc(double q, double sld_core, double radius,
double sld_solvent, double fp_n, double sld[], double thickness[])
{
const int n = (int)(fp_n+0.5);
double f, r, last_sld;
last_sld = sld_core;
f = 0.;
for (int i=0; i<n; i++) {
f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sas_3j1x_x(q*r);
last_sld = sld[i];
r += thickness[i];
}
f += M_4PI_3 * cube(r) * (sld_solvent - last_sld) * sas_3j1x_x(q*r);
return 1e-2 *f;
}

static double fqMz(double q, double sld_core, double radius,
double sld_solvent, double fp_n, double sld[], double thickness[],double eta_core,double eta_solvent,double eta[])
{
const int n = (int)(fp_n+0.5);
double f, r, last_sld;
last_sld = sld_core*langevin(eta_core);
f = 0.;
for (int i=0; i<n; i++) {
f += M_4PI_3 * cube(r) * (sld[i]*langevin(eta[i]) - last_sld) * sas_3j1x_x(q*r);
last_sld = sld[i]*langevin(eta[i]);
r += thickness[i];
}
f += M_4PI_3 * cube(r) * (sld_solvent*langevin(eta_solvent) - last_sld) * sas_3j1x_x(q*r);
return 1e-2 *f;
}

//Magnetization component wrt to field are fixed and defined by Langevin function.
//M_trans=M_x=M_y=0
//M_z= langevin(xi)
//z axis oriented always along field/polarisation axis.

//some basic vector algebra

void SET_VEC(double *vector, double v0, double v1, double v2) {
vector[0] = v0;
vector[1] = v1;
vector[2] = v2;
}

void SCALE_VEC(double *vector, double a) {
vector[0] = a*vector[0];
vector[1] = a*vector[1];
vector[2] = a*vector[2];
}

void ADD_VEC(double *result_vec, double *vec1, double *vec2) {
result_vec[0] = vec1[0] + vec2[0];
result_vec[1] = vec1[1] + vec2[1];
result_vec[2] = vec1[2] + vec2[2];
}

//Rotation of Detector around the magnetic reference frame
//Magnetic coordinates and scattering coordinates coincides when z/field axis along the neutron beam.
//then azimuthal angle theta=0Deg, rotation beta=arbitrary

static double square_mag_scat(
double cos_theta, double sin_theta,
double alpha, double beta) {
double cos_alpha, sin_alpha;
double cos_beta, sin_beta;

SINCOS(alpha*M_PI/180, sin_alpha, cos_alpha);
SINCOS(beta*M_PI/180, sin_beta, cos_beta);
//field is defined along (0,0,1), orientation of detector
//is precessing in a cone around B with an inclination of theta
double scat_vector[3];
double perp_scat_vector[3];

SET_VEC(scat_vector, cos_alpha* cos_theta,
cos_theta* sin_alpha*sin_beta +
cos_beta* sin_theta, -cos_beta* cos_theta* sin_alpha +
sin_beta* sin_theta);

return square(scat_vector[2]);

}

//weight of the magnetisation vector for the scattering, done by evaluating the magnetic scattering vector (Halpern Johnson vector) for general orientation of q. Inserting result in the spin-resolved scattering cross section
//and collecting terms

//Calculation of the 4 spin-resolved scattering cross sections for superparamagnetic spheres, field perpendicular to beam

//spin-resolved (POLARIS) cross sections
//NSF++ (F_N-Mz*(1-rotated_scat_vector[3]^2))^2+ Mx^2*(1-rotated_scat_vector[3]^2)*rotated_scat_vector[3]^2
static double Idd(double q, double sld_core,double magnetic_sld_core,double eta_core,double radius,
double sld_solvent,double magnetic_sld_solvent,double eta_solvent,double delta_solvent,
double fp_n, double sld[],double magnetic_sld[],double eta[], double delta[], double thickness[], double cos_theta, double sin_theta, double alpha, double beta) {

return fqnucsq(q, sld_core, radius, sld_solvent, fp_n, sld, thickness)-2*fqnuc(q, sld_core, radius, sld_solvent, fp_n, sld, thickness)*fqMz(q, magnetic_sld_core, radius, magnetic_sld_solvent, fp_n, magnetic_sld, thickness,eta_core,eta_solvent,eta)*(1-square_mag_scat(cos_theta, sin_theta,alpha, beta))+fqMzsq(q, magnetic_sld_core, radius, magnetic_sld_solvent, fp_n, magnetic_sld, thickness,eta_core,eta_solvent,eta)*square(1-square_mag_scat(cos_theta, sin_theta,alpha, beta))+fqMtranssq(q, magnetic_sld_core, radius, magnetic_sld_solvent, fp_n, magnetic_sld, thickness,eta_core,eta_solvent,eta, delta_solvent,delta)*(1-square_mag_scat(cos_theta, sin_theta,alpha, beta))*square_mag_scat(cos_theta, sin_theta,alpha, beta);
}

//NSF-- (F_N+Mz*(1-rotated_scat_vector[3]^2))^2+ Mx^2*(1-rotated_scat_vector[3]^2)*rotated_scat_vector[3]^2
static double Iuu(double q, double sld_core,double magnetic_sld_core,double eta_core,double radius,
double sld_solvent,double magnetic_sld_solvent,double eta_solvent,double delta_solvent,
double fp_n, double sld[],double magnetic_sld[],double eta[], double delta[], double thickness[], double cos_theta, double sin_theta, double alpha, double beta) {
return fqnucsq(q, sld_core, radius, sld_solvent, fp_n, sld, thickness)+2*fqnuc(q, sld_core, radius, sld_solvent, fp_n, sld, thickness)*fqMz(q, magnetic_sld_core, radius, magnetic_sld_solvent, fp_n, magnetic_sld, thickness,eta_core,eta_solvent,eta)*(1-square_mag_scat(cos_theta, sin_theta,alpha, beta))+fqMzsq(q, magnetic_sld_core, radius, magnetic_sld_solvent, fp_n, magnetic_sld, thickness,eta_core,eta_solvent,eta)*square(1-square_mag_scat(cos_theta, sin_theta,alpha, beta))+fqMtranssq(q, magnetic_sld_core, radius, magnetic_sld_solvent, fp_n, magnetic_sld, thickness,eta_core,eta_solvent,eta, delta_solvent,delta)*(1-square_mag_scat(cos_theta, sin_theta,alpha, beta))*square_mag_scat(cos_theta, sin_theta,alpha, beta);
}

//SF Mz^2*(1-rotated_scat_vector[3]^2)*rotated_scat_vector[3]^2 +Mx^2*(1+rotated_scat_vector[3]^4)
static double Idu(double q, double sld_core,double magnetic_sld_core,double eta_core,double radius,
double sld_solvent,double magnetic_sld_solvent,double eta_solvent,double delta_solvent,
double fp_n, double sld[],double magnetic_sld[],double eta[], double delta[], double thickness[], double cos_theta, double sin_theta, double alpha, double beta) {
return fqMtranssq(q, magnetic_sld_core, radius, magnetic_sld_solvent, fp_n, magnetic_sld,thickness,eta_core,eta_solvent,eta, delta_solvent,delta)*(1+square(square_mag_scat(cos_theta, sin_theta,alpha, beta)))+fqMzsq(q, magnetic_sld_core, radius, magnetic_sld_solvent, fp_n,magnetic_sld, thickness,eta_core,eta_solvent,eta)* (1-square_mag_scat(cos_theta, sin_theta,alpha, beta))*square_mag_scat(cos_theta, sin_theta,alpha, beta);
}

static double Iud(double q, double sld_core,double magnetic_sld_core,double eta_core,double radius,
double sld_solvent,double magnetic_sld_solvent,double eta_solvent,double delta_solvent,
double fp_n, double sld[],double magnetic_sld[],double eta[], double delta[], double thickness[], double cos_theta, double sin_theta, double alpha, double beta) {
return Idu(q, sld_core,magnetic_sld_core, eta_core, radius, sld_solvent,magnetic_sld_solvent, eta_solvent,delta_solvent,
fp_n, sld,magnetic_sld,eta,delta,  thickness, cos_theta, sin_theta, alpha, beta);
}

//weighting of spin resolved cross sections to reconstruct partially polarised beam with imperfect optics using up_i/up_f.
static void set_weights(double in_spin, double out_spin, double weight[4]) //from kernel_iq.c
{
double norm=out_spin;

in_spin = clipp(sqrt(square(in_spin)), 0.0, 1.0);//opencl has ambiguities for abs()
out_spin = clipp(sqrt(square(out_spin)), 0.0, 1.0);

if (out_spin < 0.5){norm=1-out_spin;}

// The norm is needed to make sure that the scattering cross sections are
//correctly weighted, such that the sum of spin-resolved measurements adds up to
// the unpolarised or half-polarised scattering cross section. No intensity weighting
// needed on the incoming polariser side assuming that a user has normalised
// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively.

weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd
weight[1] = (1.0-in_spin) * out_spin / norm;       // du
weight[2] = in_spin * (1.0-out_spin) / norm;       // ud
weight[3] = in_spin * out_spin / norm;             // uu
}

//calculate 2D from _fq
static double
Iqxy(double qx, double qy, double sld_core,double magnetic_sld_core,double eta_core,double radius,
double sld_solvent,double magnetic_sld_solvent,double eta_solvent,double delta_solvent,
double fp_n, double sld[],double magnetic_sld[],double eta[], double delta[], double thickness[], double up_i, double up_f, double alpha, double beta)
{
const double q = sqrt(qx*qx + qy*qy);
double cos_theta=qx/q;
double sin_theta=qy/q;

double weights[8];  // uu, ud, du, dd, fill,fill, fill, fill (make memory alloc happy)
set_weights(up_i, up_f, weights);

fp_n, sld,magnetic_sld,eta,delta,thickness,cos_theta ,sin_theta, alpha, beta) + weights[1]*Idu(q, sld_core,magnetic_sld_core,eta_core,radius,sld_solvent,magnetic_sld_solvent,eta_solvent,delta_solvent,
fp_n, sld,magnetic_sld,eta,delta,thickness,cos_theta ,sin_theta, alpha, beta) + weights[2]*Iud(q, sld_core,magnetic_sld_core,eta_core,radius,sld_solvent,magnetic_sld_solvent,eta_solvent,delta_solvent,
fp_n, sld,magnetic_sld,eta,delta,thickness,cos_theta ,sin_theta, alpha, beta) + weights[3]*Idd(q, sld_core,magnetic_sld_core,eta_core,radius,sld_solvent,magnetic_sld_solvent,eta_solvent,delta_solvent,
fp_n, sld,magnetic_sld,eta,delta,thickness,cos_theta ,sin_theta, alpha, beta);

return form;
}

//calculate 1D by averaging over theta
//TODO: choose orientation and sector width for averaging
// 2D to 1D
static double
Iq(double q, double sld_core,double magnetic_sld_core,double eta_core,double radius,
double sld_solvent,double magnetic_sld_solvent,double eta_solvent,double delta_solvent,
double fp_n, double sld[],double magnetic_sld[],double eta[], double delta[], double thickness[], double up_i, double up_f, double alpha, double beta)
{
double weights[8];  // uu, ud, du, dd, fill,fill, fill, fill (make memory alloc happy)
set_weights(up_i, up_f, weights);
double sin_theta, cos_theta; // slots to hold sincos function output of the orientation on the detector plane
double total_F2 = 0.0;
for (int i=0; i<GAUSS_N ;i++) {

const double theta = M_PI * (GAUSS_Z[i] + 1.0); // 0 .. 2 pi
SINCOS(theta, sin_theta, cos_theta);
const double form = weights[0]*Iuu(q, sld_core,magnetic_sld_core,eta_core,radius,sld_solvent,magnetic_sld_solvent,eta_solvent,delta_solvent,
fp_n, sld,magnetic_sld,eta,delta,thickness,cos_theta ,sin_theta, alpha, beta) + weights[1]*Idu(q, sld_core,magnetic_sld_core,eta_core,radius,sld_solvent,magnetic_sld_solvent,eta_solvent,delta_solvent,
fp_n, sld,magnetic_sld,eta,delta,thickness,cos_theta ,sin_theta, alpha, beta) + weights[2]*Iud(q, sld_core,magnetic_sld_core,eta_core,radius,sld_solvent,magnetic_sld_solvent,eta_solvent,delta_solvent,
fp_n, sld,magnetic_sld,eta,delta,thickness,cos_theta ,sin_theta, alpha, beta) + weights[3]*Idd(q, sld_core,magnetic_sld_core,eta_core,radius,sld_solvent,magnetic_sld_solvent,eta_solvent,delta_solvent,
fp_n, sld,magnetic_sld,eta,delta,thickness,cos_theta ,sin_theta, alpha, beta);

total_F2 += GAUSS_W[i] * form ;
}

//convert from [1e-12 A-1] to [cm-1]
return 0.5*total_F2;
}

```
```