Merge branch 'pyramidff' into develop

Core/FormFactors/inc/FormFactorAnisoPyramid.h

@@ -59,6 +59,11 @@ protected:
     virtual void init_parameters();
 
 private:
+    complex_t fullAnisoPyramidPrimitive(complex_t a, complex_t b, complex_t c,
+                                        double d, double z) const;
+    complex_t g(complex_t x, complex_t c, complex_t bd, double z) const;  // helper function
+    complex_t h(complex_t x, complex_t bd, double z) const;               // helper function
+    complex_t k(complex_t x, double d, double z) const;                   // helper function
     double m_length;
     double m_width;
     double m_height;

Core/FormFactors/inc/FormFactorPyramid.h

@@ -55,6 +55,8 @@ protected:
 
 private:
     complex_t fullPyramidPrimitive(complex_t a, complex_t b, complex_t c, double z) const;
+    complex_t g(complex_t x, complex_t c, double z) const;  // helper function
+    complex_t h(complex_t x, double z) const;  // helper function
     double m_length;
     double m_height;
     double m_alpha;

Core/FormFactors/src/FormFactorAnisoPyramid.cpp

@@ -31,6 +31,26 @@ FormFactorAnisoPyramid::FormFactorAnisoPyramid(
 bool FormFactorAnisoPyramid::check_initialization() const
 {
     bool result(true);
+    if(m_alpha<0.0 || m_alpha>Units::PID2) {
+        std::ostringstream ostr;
+        ostr << "FormFactorAnisoPyramid() -> Error in class initialization with parameters ";
+        ostr << " length:" << m_length;
+        ostr << " width:" << m_width;
+        ostr << " height:" << m_height;
+        ostr << " alpha:" << m_alpha << "\n\n";
+        ostr << "Check for '0 <= alpha <= pi/2' failed.";
+        throw Exceptions::ClassInitializationException(ostr.str());
+    }
+    if(m_length<0.0 || m_width<0.0 || m_height<0.0) {
+        std::ostringstream ostr;
+        ostr << "FormFactorAnisoPyramid() -> Error in class initialization with parameters ";
+        ostr << " length:" << m_length;
+        ostr << " width:" << m_width;
+        ostr << " height:" << m_height;
+        ostr << " alpha:" << m_alpha << "\n\n";
+        ostr << "Check for '0 <= length, 0<= width and 0 <= height' failed.";
+        throw Exceptions::ClassInitializationException(ostr.str());
+    }
     if(2.*m_height > std::min(m_length, m_width)*std::tan(m_alpha)) {
         std::ostringstream ostr;
         ostr << "FormFactorAnisoPyramid() -> Error in class initialization ";
@@ -72,104 +92,85 @@ complex_t FormFactorAnisoPyramid::evaluate_for_q(const cvector_t& q) const
     complex_t qy = q.y();
     complex_t qz = q.z();
 
-    complex_t F;
     const complex_t im(0,1);
+    if (L < W) {
+        complex_t full = fullAnisoPyramidPrimitive(qx/tga, qy/tga, qz, (W-L)*tga/2.0, -L*tga/2.0);
+        complex_t top = fullAnisoPyramidPrimitive(qx/tga, qy/tga, qz, (W-L)*tga/2.0, H-L*tga/2.0);
+        return std::exp(im*qz*L*tga/2.0) * (full-top) / (tga*tga);
+    } else {
+        complex_t full = fullAnisoPyramidPrimitive(qy/tga, qx/tga, qz, (L-W)*tga/2.0, -W*tga/2.0);
+        complex_t top = fullAnisoPyramidPrimitive(qy/tga, qx/tga, qz, (L-W)*tga/2.0, H-W*tga/2.0);
+        return std::exp(im*qz*W*tga/2.0) * (full-top) / (tga*tga);
+    }
+}
 
-    if( std::abs(qx) > Numeric::double_epsilon
-            && std::abs(qy) > Numeric::double_epsilon ) {
-        // General case
-        complex_t q1, q2, q3, q4;
-        q1=(H/2.)*((qx-qy)/tga + qz);
-        q2=(H/2.)*((qx-qy)/tga - qz);
-        q3=(H/2.)*((qx+qy)/tga + qz);
-        q4=(H/2.)*((qx+qy)/tga - qz);
-        complex_t K1,K2,K3,K4;
-        K1 = MathFunctions::Sinc(q1)*std::exp(im*q1)
-                + MathFunctions::Sinc(q2)*std::exp(-im*q2);
-        K2 = -MathFunctions::Sinc(q1)*std::exp(im*q1)*im
-                + MathFunctions::Sinc(q2)*std::exp(-im*q2)*im;
-        K3 = MathFunctions::Sinc(q3)*std::exp(im*q3)
-                + MathFunctions::Sinc(q4)*std::exp(-im*q4);
-        K4 = -MathFunctions::Sinc(q3)*std::exp(im*q3)*im
-                + MathFunctions::Sinc(q4)*std::exp(-im*q4)*im;
-        F = K1*std::cos( (qx*L-qy*W)/2. ) + K2*std::sin( (qx*L-qy*W)/2.  )
-                - K3*std::cos( (qx*L+qy*W)/2. ) - K4*std::sin( (qx*L+qy*W)/2. );
-        F = F*H/(qx*qy);
-    } else if(std::abs(qx) <= Numeric::double_epsilon
-              && std::abs(qy) <= Numeric::double_epsilon) {
-
-        if (std::abs(qz) <= Numeric::double_epsilon) {
-        //Volume qx=qy=qz=0
-        F = L*W*H - (L + W)*H*H/tga + 4.0/3.0*H*H*H/(tga*tga);
+complex_t FormFactorAnisoPyramid::fullAnisoPyramidPrimitive(complex_t a, complex_t b, complex_t c,
+                                                            double d, double z) const
+{
+    const complex_t im(0, 1);
+    if (std::norm(a*z) > Numeric::double_epsilon && std::norm(b*z) > Numeric::double_epsilon) {
+        if (std::abs((a-b)*(a-b)-c*c)*z*z > Numeric::double_epsilon &&
+            std::abs((a+b)*(a+b)-c*c)*z*z > Numeric::double_epsilon) {
+            complex_t phase = std::exp(im * c * z);
+            complex_t numerator = std::sin(a*z) * (b*(a*a - b*b + c*c)*std::cos(b*(d-z))
+                                                   - im*c*(a*a + b*b - c*c)*std::sin(b*(d-z)))
+                              + a*std::cos(a*z) * ((a*a - b*b - c*c)*std::sin(b*(d-z))
+                                                   + 2.0*im*b*c*std::cos(b*(d-z)));
+            complex_t denominator = a*b*(a - b - c)*(a + b - c)*(a - b + c)*(a + b + c);
+            return -4.0 * phase * numerator / denominator;
         } else {
-            //qx=0 qy=0 qz!=0
-             F=im*(
-                - 8.0/std::pow(tga,2) - 2.0*im*qz*(L+W)/tga + std::pow(qz,2)*L*W
-                - std::exp(im*H*qz)*(L*W*qz*qz - 2.0*(L+W)*qz/tga*(H*qz+im)
-                + 4.0*(std::pow(qz*H,2)+2.0*im*qz*H-2.0)/std::pow(tga,2))
-                   )/std::pow(qz,3);
+            return 2.0 * (g(a-b, c, b*d, z) - g(-a-b, c, b*d, z)) / (a*b);
         }
-
+    } else if (std::norm(a*z) <= Numeric::double_epsilon
+               && std::norm(b*z) <= Numeric::double_epsilon) {
+        if (std::norm(c*z) <= Numeric::double_epsilon) {
+            return 2.0*(d - 2.0*z/3.0)*z*z;
+        } else
+            return 4.0*(2.0*im -c*d -im*std::exp(im*c*z)*(c*c*(d-z)*z + im*c*(d-2.0*z) + 2.0))
+                    / std::pow(c, 3);
     } else {
-
-        complex_t qxy, q5, q6;
-         double R0, Rxy;
-        if(std::abs(qy) <= Numeric::double_epsilon
-                && std::abs(qx) > Numeric::double_epsilon) {
-            // qx!=0 qy=0
-            qxy=qx;
-            Rxy = L/2.0;
-            R0 = W/2.0;
+        if (std::norm(a*z) <= Numeric::double_epsilon) {
+            return 2.0 * (h(c+b, b*d, z) - h(c-b, -b*d, z)) / b;
         } else {
-            // qy!=0 qx=0
-            qxy=qy;
-            R0 = L/2.0;
-            Rxy = W/2.0;
+            return 2.0 * (k(c+a, d, z) - k(c-a, d, z)) / a;
         }
+    }
+}
 
-        q5 = (qz - qxy/tga)/2.;
-        q6 = (qz + qxy/tga)/2.;
-
-        if (std::abs(q5) <= Numeric::double_epsilon) {
-            //q5 = 0, q6!=0
-        F= -2.0*H*im/qxy*(
-                    (R0 - H/(2.0*tga))*std::exp(im*qxy*Rxy)
-                    - std::exp(im*q6*H-im*qxy*Rxy)*
-                    ( R0*MathFunctions::Sinc(q6*H)
-                    + im*( std::exp(im*q6*H)
-                           - MathFunctions::Sinc(q6*H) )/(2.0*q6*tga))
-                    );
-
+complex_t FormFactorAnisoPyramid::g(complex_t x, complex_t c, complex_t bd, double z) const
+{
+    const complex_t im(0, 1);
+    if (std::abs((x*x-c*c)*z*z) > Numeric::double_epsilon) {
+        return (im*c*std::cos(bd) + x*std::sin(bd)
+                - std::exp(im*c*z)*(im*c*std::cos(bd+x*z) + x*std::sin(bd+x*z))) / (x*x-c*c);
+    } else {
+        if (std::norm(c*z) > Numeric::double_epsilon) {
+            return (im * (std::exp(2.0*im*c*z) + 2.0*im*c*z - 1.0) * std::cos(bd)
+                       - (std::exp(2.0*im*c*z) - 2.0*im*c*z - 1.0) * std::sin(bd) ) / (4.0*c);
         } else {
-
-            if (std::abs(q6) <= Numeric::double_epsilon) {
-
-                //q5!= 0, q6=0
-                F=-2.0*H*im/qxy*(
-                            std::exp(im*q5*H+im*qxy*Rxy)*
-                            ( R0*MathFunctions::Sinc(q5*H)
-                            + im*( std::exp(im*q5*H)
-                                   - MathFunctions::Sinc(q5*H) )/(2.0*q5*tga))
-                            +(H/(2.0*tga)-R0)*std::exp(-im*qxy*Rxy)
-                            );
-
-            } else {
-
-                //q5!= 0, q6!=0
-                F=-2.0*H*im/qxy*(
-                            std::exp(im*q5*H+im*qxy*Rxy)*
-                            (R0*MathFunctions::Sinc(q5*H)
-                            + im*( std::exp(im*q5*H)
-                                   - MathFunctions::Sinc(q5*H) )/(2.0*q5*tga))
-                            - std::exp(im*q6*H-im*qxy*Rxy)
-                            *(R0*MathFunctions::Sinc(q6*H)
-                            + im*( std::exp(im*q6*H)
-                                   - MathFunctions::Sinc(q6*H))/(2.0*q6*tga))
-                            );
-            }
+            return -z * std::cos(bd);
+        }
+    }
 }
+
+complex_t FormFactorAnisoPyramid::h(complex_t x, complex_t bd, double z) const
+{
+    const complex_t im(0, 1);
+    if (std::norm(x*z) > Numeric::double_epsilon) {
+        return std::exp(-im*bd)*(std::exp(im*x*z)*(im+x*z) - im) / (x*x);
+    } else {
+        return im*std::exp(-im*bd)*z*z/2.0;
+    }
 }
-    return F;
+
+complex_t FormFactorAnisoPyramid::k(complex_t x, double d, double z) const
+{
+    const complex_t im(0, 1);
+    if (std::norm(x*z) > Numeric::double_epsilon) {
+        return (x*d - im + std::exp(im*x*z)*(im + x*(z-d))) / (x*x);
+    } else {
+        return im*(d - z/2.0)*z;
+    }
 }
 
 

Core/FormFactors/src/FormFactorPyramid.cpp

@@ -30,6 +30,24 @@ FormFactorPyramid::FormFactorPyramid(
 bool FormFactorPyramid::check_initialization() const
 {
     bool result(true);
+    if(m_alpha<0.0 || m_alpha>Units::PID2) {
+        std::ostringstream ostr;
+        ostr << "FormFactorPyramid() -> Error in class initialization with parameters ";
+        ostr << " length:" << m_length;
+        ostr << " height:" << m_height;
+        ostr << " alpha:" << m_alpha << "\n\n";
+        ostr << "Check for '0 <= alpha <= pi/2' failed.";
+        throw Exceptions::ClassInitializationException(ostr.str());
+    }
+    if(m_length<0.0 || m_height<0.0) {
+        std::ostringstream ostr;
+        ostr << "FormFactorPyramid() -> Error in class initialization with parameters ";
+        ostr << " length:" << m_length;
+        ostr << " height:" << m_height;
+        ostr << " alpha:" << m_alpha << "\n\n";
+        ostr << "Check for '0 <= length and 0 <= height' failed.";
+        throw Exceptions::ClassInitializationException(ostr.str());
+    }
     if(2.*m_height > m_length*std::tan(m_alpha)) {
         std::ostringstream ostr;
         ostr << "FormFactorPyramid() -> Error in class initialization with parameters ";
@@ -64,64 +82,73 @@ complex_t FormFactorPyramid::evaluate_for_q(const cvector_t& q) const
     double H = m_height;
     double R = m_length/2.;
     double tga = std::tan(m_alpha);
-    //TODO: check for tga==0 or tga<0 or tga->+-infinity
+    if (m_height == 0.0) return 0.0;
+    if (m_alpha == 0.0) return 0.0;
+    if (m_length == 0.0) return 0.0;
 
     complex_t qx = q.x();
     complex_t qy = q.y();
     complex_t qz = q.z();
 
-    complex_t F;
     const complex_t im(0, 1);
-    //TODO: the following checks came from the previous implementation and will only catch
-    //part of the numerical instabilities
-    if (std::norm(qx) > Numeric::double_epsilon && std::norm(qy) > Numeric::double_epsilon) {
-        complex_t full = fullPyramidPrimitive(qx/tga, qy/tga, qz, -R*tga);
-        complex_t top = fullPyramidPrimitive(qx/tga, qy/tga, qz, H-R*tga);
-        F = std::exp(im*qz*R*tga)*(full-top)/(tga*tga);
-    } else if (std::norm(qx) <= Numeric::double_epsilon
-               && std::norm(qy) <= Numeric::double_epsilon) {
-        if (std::norm(qz) <= Numeric::double_epsilon)
-            F = 4. / 3. * tga * R * R * R
-                * (1. - (1. - H / R / tga) * (1. - H / R / tga) * (1. - H / R / tga));
-        else
-            F = 4. * im * (-2. / tga / tga - 2. * im * qz * R / tga + qz * qz * R * R
-                           - std::exp(im * H * qz) * ((-1. + im + H * qz) / tga - qz * R)
-                             * ((1. + im + H * qz) / tga - qz * R)) / std::pow(qz, 3);
+    complex_t full = fullPyramidPrimitive(qx/tga, qy/tga, qz, -R*tga);
+    complex_t top = fullPyramidPrimitive(qx/tga, qy/tga, qz, H-R*tga);
+    return std::exp(im*qz*R*tga) * (full-top) / (tga*tga);
+}
+
+complex_t FormFactorPyramid::fullPyramidPrimitive(complex_t a, complex_t b, complex_t c,
+                                                          double z) const
+{
+    const complex_t im(0, 1);
+    if (std::norm(a*z) > Numeric::double_epsilon && std::norm(b*z) > Numeric::double_epsilon) {
+        if (std::abs((a-b)*(a-b)-c*c)*z*z > Numeric::double_epsilon &&
+            std::abs((a+b)*(a+b)-c*c)*z*z > Numeric::double_epsilon) {
+            complex_t phase = std::exp(im * c * z);
+            complex_t numerator = std::sin(a*z) * (b*(a*a - b*b + c*c)*std::cos(b*z)
+                                                   + im*c*(a*a + b*b - c*c)*std::sin(b*z))
+                              + a*std::cos(a*z) * ((-a*a + b*b + c*c)*std::sin(b*z)
+                                                   + 2.0*im*b*c*std::cos(b*z));
+            complex_t denominator = a * b * (a - b - c) * (a + b - c) * (a - b + c) * (a + b + c);
+            return -4.0 * phase * numerator / denominator;
+        } else {
+            return 2.0 * (g(a-b, c, z) - g(a+b, c, z)) / (a*b);
+        }
+    } else if (std::norm(a*z) <= Numeric::double_epsilon
+               && std::norm(b*z) <= Numeric::double_epsilon) {
+        if (std::norm(c*z) <= Numeric::double_epsilon) {
+            return -4.0*std::pow(z, 3)/3.0;
+        } else
+            return 4.0*im * (2.0 + std::exp(im*c*z)*(c*c*z*z + 2.0*im*c*z - 2.0)) / std::pow(c, 3);
     } else {
-        complex_t qxy;
-        if (std::norm(qy) <= Numeric::double_epsilon && std::norm(qx) > Numeric::double_epsilon) {
-            qxy = qx;
+        complex_t abmax;
+        if (std::norm(b*z) <= Numeric::double_epsilon) {
+            abmax = a;
         } else {
-            qxy = qy;
+            abmax = b;
         }
-        F = (4. * (qxy * tga * (-(qxy * qxy * R) + qz * tga * (complex_t(0.0, -2.0) + qz * R * tga))
-                   * std::cos(qxy * R)
-                   - std::exp(im * H * qz) * qxy
-                     * (H * std::pow(qxy, 2) - qxy * qxy * R * tga
-                        - qz * (complex_t(0.0, 2.0) + H * qz) * std::pow(tga, 2)
-                        + std::pow(qz, 2) * R * std::pow(tga, 3)) * std::cos(qxy * (R - H / tga))
-                   + tga * (std::pow(qxy, 2) * (1. - complex_t(0.0, 1.0) * qz * R * tga)
-                            + std::pow(qz, 2) * std::pow(tga, 2)
-                              * (1. + complex_t(0.0, 1.0) * qz * R * tga)) * std::sin(qxy * R)
-                   + complex_t(0.0, 1.0) * std::exp(im * H * qz) * tga
-                     * (std::pow(qz, 2) * std::pow(tga, 2)
-                        * (complex_t(0.0, 1.0) + H * qz - qz * R * tga)
-                        + std::pow(qxy, 2) * (complex_t(0.0, 1.0) - H * qz + qz * R * tga))
-                     * std::sin(qxy * (R - H / tga))))
-            / (qxy * std::pow(qxy - qz * tga, 2) * std::pow(qxy + qz * tga, 2));
+        return 2.0 * (h(c - abmax, z) - h(c + abmax, z)) / abmax;
     }
-    return F;
 }
 
-complex_t FormFactorPyramid::fullPyramidPrimitive(complex_t a, complex_t b, complex_t c,
-                                                  double z) const
+complex_t FormFactorPyramid::g(complex_t x, complex_t c, double z) const
 {
     const complex_t im(0, 1);
-    complex_t phase = std::exp(im * c * z);
-    complex_t numerator = std::sin(a * z) * (b * (a * a - b * b + c * c) * std::cos(b * z)
-                                             + im * c * (a * a + b * b - c * c) * std::sin(b * z))
-                          + a * std::cos(a * z) * ((-a * a + b * b + c * c) * std::sin(b * z)
-                                                   + 2.0 * im * b * c * std::cos(b * z));
-    complex_t denominator = a * b * (a - b - c) * (a + b - c) * (a - b + c) * (a + b + c);
-    return -4.0 * phase * numerator / denominator;
+    if (std::abs((x*x-c*c)*z*z) > Numeric::double_epsilon) {
+        return (im*c - std::exp(im*c*z)*(x*std::sin(x*z) + im*c*std::cos(x*z))) / (x*x-c*c);
+    } else {
+        if (std::norm(c*z) > Numeric::double_epsilon) {
+            return im * (std::exp(2.0*im*c*z) + 2.0*im*c*z - 1.0) / (4.0*c);
+        } else {
+            return -z;
+        }
+    }
+}
+
+complex_t FormFactorPyramid::h(complex_t x, double z) const
+{
+    const complex_t im(0, 1);
+    if (std::norm(x*z) > Numeric::double_epsilon) {
+        return (im - (im+x*z)*std::exp(im*x*z))/(x*x);
+    }
+    return -im*z*z/2.0 + z*z*z*x/3.0;
 }

Doc/UserManual/FormFactors.tex

@@ -1082,34 +1082,21 @@ They must fulfill
 \paragraph{Form factor etc}\strut\\
 Notation:
 \begin{displaymath}
-  \ell\coloneqq L/2,\quad
-  h\coloneqq H/2,\quad
-  f_\pm(z)\coloneqq \exp(\pm i z)\sinc(z).
+  a\coloneqq q_x/\tan\alpha,\quad
+  b\coloneqq q_y/\tan\alpha,\quad
+  c\coloneqq q_z.
 \end{displaymath}
 Results:
 \begin{equation*}
-\begin{array}{@{}l@{}l@{}}
-\DS F=
-\frac{H}{q_xq_y} \Big\{
-   &+\DS  f_+\left(\left(\frac{q_x-q_y}{\tan\alpha} +q_z\right)h\right)
-        \exp(-i(q_x-q_y)\ell)
-\\[3.6ex]
-   &\DS+ f_-\left(\left(\frac{q_x-q_y}{\tan\alpha} -q_z\right)h\right)
-        \exp(+i(q_x-q_y)\ell)
-\\[3.6ex]
-   &\DS- f_+\left(\left(\frac{q_x+q_y}{\tan\alpha} +q_z\right)h\right)
-        \exp(-i(q_x+q_y)\ell)
-\\[3.6ex]
-   &\DS- f_-\left(\left(\frac{q_x+q_y}{\tan\alpha} -q_z\right)h\right)
-        \exp(+i(q_x+q_y)\ell)
-\Big\},
-\end{array}
+  \DS F= \frac{\exp(iq_zL\tan\alpha/2)}{\tan^2\alpha} \Big\{ I_p(H-L\tan\alpha/2) - I_p(-L\tan\alpha/2) \Big\} 
 \end{equation*}
+Where $I_p(z)$ is an antiderivative for the final z--integration in the Fourier transform of the non--truncated pyramid shape function:
 \begin{equation*}
-  V= H \Big[L^2 - \frac{2LH}{\tan\alpha} + \dfrac{4}{3} \dfrac{H^2}{\tan^2\alpha}\Big].
-\end{equation*}
-\begin{equation*}
-  S=L^2.
+\begin{array}{@{}l@{}l@{}}
+  \DS I_p(z) = & 4\exp(icz)\Big\{ \sin(az)\big[b(a^2-b^2+c^2)\cos(bz) + ic(a^2+b^2-c^2)\sin(bz)\big] \\
+  & +a\cos(az)\big[(-a^2+b^2+c^2)\sin(bz) + 2ibc\cos(bz) \big] \Big\} \\
+  & / \Big\{ab\big[(a+b)^2-c^2\big]\big[(a-b)^2-c^2\big]\Big\}
+\end{array}
 \end{equation*}
 \begin{equation*}
   V = \dfrac{1}{6}  L^3 \tan\alpha\left[ 1

Examples/python/utils/plot_intensity_data_diff.py

+#!/usr/bin/env python
+# Plots intensity data difference stored in BornAgain "*.int" or "*.int.gz" format
+# Usage: python plot_intensity_data.py intensity_reference.int.gz intensity_other.int.gz
+
+import numpy
+import matplotlib
+import pylab
+from bornagain import *
+
+
+def plot_intensity_data(ref, data):
+    phi_min = rad2deg(ref.getAxis(0).getMin())
+    phi_max = rad2deg(ref.getAxis(0).getMax())
+    alpha_min = rad2deg(ref.getAxis(1).getMin())
+    alpha_max = rad2deg(ref.getAxis(1).getMax())
+    im = pylab.imshow(numpy.rot90(data, 1), norm=matplotlib.colors.LogNorm(),
+                      extent=[phi_min, phi_max, alpha_min, alpha_max])
+    cb = pylab.colorbar(im)
+    cb.set_label(r'Intensity (arb. u.)', size=16)
+    pylab.xlabel(r'$\phi_f (^{\circ})$', fontsize=16)
+    pylab.ylabel(r'$\alpha_f (^{\circ})$', fontsize=16)
+    pylab.show()
+
+
+if __name__ == '__main__':
+    if len(sys.argv)!=3:
+        exit("Usage: python plot_intensity_data_diff.py intensity_reference.int.gz intensity_other.int.gz")
+    intensity_ref = IntensityDataIOFactory.readIntensityData(sys.argv[1])
+    intensity_other = IntensityDataIOFactory.readIntensityData(sys.argv[2])
+    data = numpy.abs((intensity_ref.getArray() - intensity_other.getArray())/intensity_ref.getArray())
+    plot_intensity_data(intensity_ref, data)