Changeset 318b5bbb in sasview for sansmodels/src/c_models/libfunc.c

Timestamp:

Dec 18, 2012 10:55:24 AM (12 years ago)

Author:

Jae Cho <jhjcho@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

6550b64

Parents:

0203ade

Message:

Added polarization and magnetic stuffs

File:

• sansmodels/src/c_models/libfunc.c (modified) (4 diffs)

sansmodels/src/c_models/libfunc.c

@@ -103 +103 @@
     return -tmp+log(2.5066282746310005*ser/x);
 }
+
+// calculate magnetic sld and return total sld
+// bn : contrast (not just sld of the layer)
+// m0: max mag of M; mtheta: angle from x-z plane;
+// mphi: angle (anti-clock-wise)of x-z projection(M) from x axis
+// spinfraci: the fraction of UP among UP+Down (before sample)
+// spinfracf: the fraction of UP among UP+Down (after sample and before detector)
+// spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis
+// Note: all angles are in degrees.
+polar_sld cal_msld(int isangle, double qx, double qy, double bn,
+                                double m01, double mtheta1, double mphi1,
+                                double spinfraci, double spinfracf, double spintheta)
+{
+        //locals
+        double q_x = qx;
+        double q_y = qy;
+        double sld = bn;
+        int is_angle = isangle;
+        double pi = 4.0*atan(1.0);
+        double s_theta = spintheta * pi/180.0;
+        double m_max = m01;
+        double m_phi = mphi1;
+        double m_theta = mtheta1;
+        double in_spin = spinfraci;
+        double out_spin = spinfracf;
+
+        double m_perp = 0.0;
+        double m_perp_z = 0.0;
+        double m_perp_y = 0.0;
+        double m_perp_x = 0.0;
+        double m_sigma_x = 0.0;
+        double m_sigma_z = 0.0;
+        double m_sigma_y = 0.0;
+        double b_m = 0.0;
+        double q_angle = 0.0;
+        double mx = 0.0;
+        double my = 0.0;
+        double mz = 0.0;
+        polar_sld p_sld;
+        p_sld.uu = sld;
+        p_sld.dd = sld;
+        p_sld.re_ud = 0.0;
+        p_sld.im_ud = 0.0;
+        p_sld.re_du = 0.0;
+        p_sld.im_du = 0.0;
+
+        //No mag means no further calculation
+        if (isangle>0){
+                if (m_max < 1.0e-32){
+                        p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+                        p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+                        return p_sld;
+                }
+        }
+        else{
+                if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){
+                        p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+                        p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+                        return p_sld;
+                }
+        }
+
+        //These are needed because of the precision of inputs
+        if (in_spin < 0.0) in_spin = 0.0;
+        if (in_spin > 1.0) in_spin = 1.0;
+        if (out_spin < 0.0) out_spin = 0.0;
+        if (out_spin > 1.0) out_spin = 1.0;
+
+        if (q_x == 0.0) q_angle = pi / 2.0;
+        else q_angle = atan(q_y/q_x);
+        if (q_y < 0.0 && q_x < 0.0) q_angle -= pi;
+        else if (q_y > 0.0 && q_x < 0.0) q_angle += pi;
+
+        q_angle = pi/2.0 - q_angle;
+        if (q_angle > pi) q_angle -= 2.0 * pi;
+        else if (q_angle < -pi) q_angle += 2.0 * pi;
+
+        if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){
+                m_perp = 0.0;
+                }
+        else {
+                m_perp = m_max;
+                }
+        if (is_angle > 0){
+                m_phi *= pi/180.0;
+                m_theta *= pi/180.0;
+                mx = m_perp * cos(m_theta) * cos(m_phi);
+                my = m_perp * sin(m_theta);
+                mz = -(m_perp * cos(m_theta) * sin(m_phi));
+        }
+        else{
+                mx = m_perp;
+                my = m_phi;
+                mz = m_theta;
+        }
+        //ToDo: simplify these steps
+        // m_perp1 -m_perp2
+        m_perp_x = (mx) *  cos(q_angle);
+        m_perp_x -= (my) * sin(q_angle);
+        m_perp_y = m_perp_x;
+        m_perp_x *= cos(-q_angle);
+        m_perp_y *= sin(-q_angle);
+        m_perp_z = mz;
+
+        m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta));
+        m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta));
+        m_sigma_z = (m_perp_z);
+
+        //Find b
+        p_sld.uu -= m_sigma_x;
+        p_sld.dd += m_sigma_x;
+        p_sld.re_ud = m_sigma_y;
+        p_sld.re_du = m_sigma_y;
+        p_sld.im_ud = m_sigma_z;
+        p_sld.im_du = -m_sigma_z;
+
+        p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+        p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+
+        p_sld.re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.re_ud;
+        p_sld.im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.im_ud;
+        p_sld.re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.re_du;
+        p_sld.im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.im_du;
+
+        return p_sld;
+}
+
 
 /** Modifications below by kieranrcampbell@gmail.com
@@ -131 +258 @@
     ap = a;
     del = sum = 1.0/a;
-    
+
     for(n=1;n<=ITMAX;n++) {
       ++ap;
@@ -150 +277 @@
 
 /**
-   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction 
-   representation 
+   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction
+   representation
 **/
 
@@ -196 +323 @@
 }
 
-/** 
+/**
     Implementation of the error function, erf(x)
 **/