Field-dependent magnetic SANS of misaligned magnetic moments in bulk ferromagnets - micromagnetic_misalignment_3D.c

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


//Core-shell form factor for anisotropy field (Hkx, Hky and Hkz), Nuc and lonitudinal magnetization Mz

static double
form_volume(double radius, double thickness)
{
 return M_4PI_3 * cube(radius + thickness);
}

static double
radius_effective(int mode, double radius, double thickness)
{
 switch (mode) {
   default:
   case 1: // outer radius
   return radius + thickness;
   case 2: // core radius
   return radius;
 }
}

static double fq(double q, double radius,
 double thickness, double core_sld, double shell_sld, double solvent_sld)
{
  const double form = core_shell_fq(q,
    radius,
    thickness,
    core_sld,
    shell_sld,
    solvent_sld);
  return form;
}


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); 

 }

 static double DMI_length(double Ms, double D, double qval)
 {

   return 2.0 * D * 4.0 * M_PI / Ms / Ms * qval ; //q in 10e10 m-1, A in 10e-3 J/m^2, mu0 in 4 M_PI 1e-7 
 }

//Mz is defined as the longitudinal magnetisation component along the magnetic field. In the approach to saturation this component is (almost) constant with magnetic field and simplfy reflects the nanoscale variations in the saturation magnetisation value in the sample.
// The misalignment of the magnetisation due to perturbing magnetic anisotropy or dipolar magnetic fields in the sample enter Mx and My, the two transversal magnetisation components, reacting to a magnetic field.
//The micromagnetic solution for the magnetisation are from Michels et al. PRB 94, 054424 (2016).

 static double fqMxreal( double x,double y, double z, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
 {
  const double qsq = x * x + y * y + z * z; 
  const double q=sqrt(qsq);
  const double f = reduced_field(q, Ms, Hi, A)*(Hkx*(1.0+reduced_field(q, Ms, Hi, A)*y*y/qsq)-Ms*Mz*x*z/qsq*(1.0+reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D, q)*DMI_length(Ms, D, q))-Hky*reduced_field(q, Ms, Hi, A)*x*y/qsq)/(1.0+reduced_field(q, Ms, Hi, A)*(x*x+y*y)/qsq -square(reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z)));
  return f;
}

static double fqMximag( double x, double y, double z, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
{
  const double qsq=x*x+y*y+z*z; 
  const double q=sqrt(qsq);    
  const double f = -reduced_field(q, Ms, Hi, A)*(Ms*Mz*(1.0+reduced_field(q, Ms, Hi, A))*DMI_length(Ms, D,y)+Hky*reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z))/(1.0+reduced_field(q, Ms, Hi, A)*(x*x+y*y)/qsq -square(reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z)));
  return f;
}

static double fqMyreal( double x, double y, double z, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
{
  const double qsq=x*x+y*y+z*z; 
  const double q=sqrt(qsq);   
  const double f = reduced_field(q, Ms, Hi, A)*(Hky*(1.0+reduced_field(q, Ms, Hi, A)*x*x/qsq)-Ms*Mz*y*z/qsq*(1.0+reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,q)*DMI_length(Ms, D,q))-Hkx*reduced_field(q, Ms, Hi, A)*x*y/qsq)/(1.0+reduced_field(q, Ms, Hi, A)*(x*x+y*y)/qsq -square(reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z)));
  return f;
}

static double fqMyimag( double x, double y, double z, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
{
  const double qsq=x*x+y*y+z*z; 
  const double q=sqrt(qsq);    
  const double f = reduced_field(q, Ms, Hi, A)*(Ms*Mz*(1.0+reduced_field(q, Ms, Hi, A))*DMI_length(Ms, D,x)-Hkx*reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z))/(1.0+reduced_field(q, Ms, Hi, A)*(x*x+y*y)/qsq -square(reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z)));
  return f;
}


//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];
}

static double SCALAR_VEC( double *vec1, double *vec2)
{
 return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2];
}

static double MAG_VEC( double *vec)
{
 return sqrt(SCALAR_VEC(vec,vec));
}

void ORTH_VEC(double *result_vec, double *vec1, double *vec2)
{
 result_vec[0] = vec1[0] - SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2) * vec2[0];
 result_vec[1] = vec1[1] - SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2) * vec2[1];
 result_vec[2] = vec1[2] - SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2) * vec2[2];
}


//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[8]) //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.real
   weight[1] = weight[0]; // dd.imag
   weight[2] = in_spin * out_spin / norm;             // uu.real
   weight[3] = weight[2];             // uu.imag
   weight[4] = (1.0-in_spin) * out_spin / norm;       // du.real
   weight[5] = weight[4];       // du.imag 
   weight[6] = in_spin * (1.0-out_spin) / norm;       // ud.real
   weight[7] = weight[6];       // ud.imag
 }

//transforms scattering vector q in polarisation/magnetisation coordinate system
 static void set_scatvec(double *qmag, double q,
  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

  qmag[0] = q*(cos_alpha * cos_theta);
  qmag[1] = q*(cos_theta * sin_alpha*sin_beta + 
  cos_beta * sin_theta);
  qmag[2] = q*(-cos_beta * cos_theta* sin_alpha + 
  sin_beta * sin_theta);

}


//Evaluating the magnetic scattering vector (Halpern Johnson vector) for general orientation of q and collecting terms for the spin-resolved (POLARIS) cross sections. Mz is along the applied magnetic field direction, which is also the polarisation direction.
static void mag_sld(
   // 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 x, double y, double z,
  double mxreal, double myreal, double mximag, double myimag, double mz, double nuc, double sld[8])
{
  double vector[3];
  //The (transversal) magnetisation and hence the magnetic scattering sector is here a complex quantity. The spin-flip (magnetic) scattering amplitude is given with
  //  MperpPperpQ \pm i MperpP  (Moon-Riste-Koehler Phys Rev 181, 920, 1969) with Mperp and MperpPperpQ the
  //magnetisation scattering vector components perpendicular to the polarisation/field direction. Collecting terms in SF that are real (MperpPperpQreal + SCALAR_VEC(MperpPimag,qvector) )  and imaginary (MperpPperpQimag \pm SCALAR_VEC(MperpPreal,qvector) )
  double Mvectorreal[3];
  double Mvectorimag[3];
  double Pvector[3];
  double Mperpreal[3];
  double MperpPreal[3];
  double MperpPperpQreal[3];
  double Mperpimag[3];
  double MperpPimag[3];
  double MperpPperpQimag[3];

  const double q = sqrt(x*x + y*y + z*z);
  SET_VEC(vector, x/q, y/q, z/q);	
   //Moon-Riste-Koehler choose z as pointing along field/polarisation axis  
  SET_VEC(Mvectorreal, mxreal, myreal, mz);
  SET_VEC(Mvectorimag, mximag, myimag, 0);
  SET_VEC(Pvector, 0, 0, 1);
   //Magnetic scattering vector Mperp could be simplified like in Moon-Riste-Koehler
   //leave the generic computation just to check
  ORTH_VEC(Mperpreal, Mvectorreal, vector);
  ORTH_VEC(MperpPreal, Mperpreal, Pvector);
  ORTH_VEC(MperpPperpQreal, MperpPreal, vector);
  ORTH_VEC(Mperpimag, Mvectorimag, vector);
  ORTH_VEC(MperpPimag, Mperpimag, Pvector);
  ORTH_VEC(MperpPperpQimag, MperpPimag, vector);
  
  
   sld[0] = nuc - SCALAR_VEC(Pvector,Mperpreal); // dd.real => sld - D Pvector \cdot Mperp 
   sld[1] = 0; //dd.imag  = nuc_img - SCALAR_VEC(Pvector,Mperpimg); nuc_img only exist for noncentrosymmetric nuclear structures; Mperp along P vector is per definition ~ Mz assumed to be real here
   sld[2] = nuc + SCALAR_VEC(Pvector,Mperpreal);              // uu => sld + D Pvector \cdot Mperp
   sld[3] = 0; //uu.imag

   sld[4] = MAG_VEC(MperpPperpQreal)+SCALAR_VEC(MperpPimag,vector);       // du.real => length of vector MperpPperpQ:  | MperpP- (MperpP\cdot qvector)  qvector | and Mperp perpendicular to P and along q
   sld[5] = MAG_VEC(MperpPperpQimag)-SCALAR_VEC(MperpPreal,vector);       // du.imag => | MperpPimag- (MperpPimag\cdot qvector)  qvector | - i MperpPreal \cdot qvector
   sld[6] = MAG_VEC(MperpPperpQreal)-SCALAR_VEC(MperpPimag,vector);      // ud.real =>  length of vector MperpPperpQ
   sld[7] = MAG_VEC(MperpPperpQimag)+SCALAR_VEC(MperpPreal,vector);         // ud.imag => | MperpPimag- (MperpPimag\cdot qvector)  qvector | + i MperpPreal \cdot qvector

 }


//calculate 2D from _fq
 static double
 Iqxy(double qx, double qy, double radius, double thickness,double core_nuc, double shell_nuc, double solvent_nuc, double core_Ms, double shell_Ms, double solvent_Ms, double core_hk,  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 qmag[3];
     set_scatvec(qmag,q,cos_theta, sin_theta, alpha, beta);

   double weights[8];  // 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 mz=fq(q, radius, thickness, core_Ms, shell_Ms, solvent_Ms);
   double nuc=fq(q, radius, thickness, core_nuc, shell_nuc, solvent_nuc);

   double cos_gamma, sin_gamma;
   double sld[8];
   double sld_res[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;
   double total_F2_res = 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);	
      
      double Hkx= fq(q, radius, thickness, core_hk, 0, 0)*sin_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 Hky= fq(q, radius, thickness, core_hk, 0, 0)*cos_gamma;

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




      
      mag_sld(qmag[0], qmag[1], qmag[2], mxreal, myreal, mximag, myimag, mz, nuc, sld);
      mag_sld(qmag[0], qmag[1], qmag[2], 0, 0, 0, 0, mz, nuc, sld_res);//Magnetic SLD at saturation
      double form = 0.0;
      double form_res = 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];
         form_res += weights[xs]*sld_res[xs]*sld_res[xs];
       }
     }
     total_F2 += GAUSS_W[i] * form ;
     total_F2_res += GAUSS_W[i] * form_res ;     
   }

   return 0.5*1.0e-4*(total_F2-total_F2_res);}
 }



//calculate 1D by averaging over theta
//TODO: choose orientation and sector width for averaging
// 2D to 1D
 static double
 Iq(double q, double radius, double thickness,double core_nuc, double shell_nuc, double solvent_nuc, double core_Ms, double shell_Ms, double solvent_Ms, double core_hk, double Hi, double Ms, double A, double D,  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;
  double total_F1D_res = 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 qmag[3];
      set_scatvec(qmag,q,cos_theta, sin_theta, alpha, beta);

   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);
   
  double mz=fq(q, radius, thickness, core_Ms, shell_Ms, solvent_Ms);
 double nuc=fq(q, radius, thickness, core_nuc, shell_nuc, solvent_nuc);

   double cos_gamma, sin_gamma;
   double sld[8];
   double sld_res[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;
   double total_F2_res = 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);	
      
      double Hkx= fq(q, radius, thickness, core_hk, 0, 0)*sin_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 Hky= fq(q, radius, thickness, core_hk, 0, 0)*cos_gamma;

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




      
      mag_sld(qmag[0], qmag[1], qmag[2], mxreal, myreal, mximag, myimag, mz, nuc, sld);
      mag_sld(qmag[0], qmag[1], qmag[2], 0, 0, 0, 0, mz, nuc, sld_res);//Magnetic SLD at saturation
      double form = 0.0;
      double form_res = 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];
        form_res += weights[xs]*sld_res[xs]*sld_res[xs];        
      }
    }
    total_F2 += GAUSS_W[i] * form ;
    total_F2_res += GAUSS_W[i] * form_res ;    
  }

  total_F1D += GAUSS_W[j] * total_F2 ;
  total_F1D_res += GAUSS_W[j] * total_F2_res ; 
}

   //convert from [1e-12 A-1] to [cm-1]
return 0.25*1.0e-4*(total_F1D-total_F1D_res);
}







Back to Model Download