r21223 - in /branches/relax_disp: lib/dispersion/mq_ns_cpmg_2site.py target_functions/relax_disp.py -- October 24, 2013 - 17:14

Author: bugman
Date: Thu Oct 24 17:14:07 2013
New Revision: 21223

URL: http://svn.gna.org/viewcvs/relax?rev=21223&view=rev
Log:
Bug fix for the 'MQ NS CPMG 2-site' model.

This was found with the aid of private feedback from Dmitry Korzhnev and him 
emailing his cpmg_fitd9
program.  The problem is that he defines the 'n' parameter as half of a CPMG 
block.  The code was
however assuming that 'n' is a full CPMG block!


Modified:
    branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py
    branches/relax_disp/target_functions/relax_disp.py

Modified: branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py?rev=21223&r1=21222&r2=21223&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py (original)
+++ branches/relax_disp/lib/dispersion/mq_ns_cpmg_2site.py Thu Oct 24 
17:14:07 2013
@@ -71,7 +71,7 @@
     matrix[1, 1] = -k_BA - 1.j*(dwH + dw) - r20
 
 
-def r2eff_mq_ns_cpmg_2site(M0=None, F_vector=array([1, 0], float64), 
m1=None, m2=None, r20=None, pA=None, pB=None, dw=None, dwH=None, k_AB=None, 
k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, 
power=None):
+def r2eff_mq_ns_cpmg_2site(M0=None, F_vector=array([1, 0], float64), 
m1=None, m2=None, r20=None, pA=None, pB=None, dw=None, dwH=None, k_AB=None, 
k_BA=None, inv_tcpmg=None, tcp=None, back_calc=None, num_points=None, n=None):
     """The 2-site numerical solution to the Bloch-McConnell equation.
 
     This function calculates and stores the R2eff values.
@@ -107,8 +107,8 @@
     @type back_calc:        numpy rank-1 float array
     @keyword num_points:    The number of points on the dispersion curve, 
equal to the length of the tcp and back_calc arguments.
     @type num_points:       int
-    @keyword power:         The matrix exponential power array.
-    @type power:            numpy int16, rank-1 array
+    @keyword n:             The n value whereby one CPMG block is defined at 
2n.
+    @type n:                numpy int16, rank-1 array
     """
 
     # Populate the m1 and m2 matrices (only once per function call for 
speed).
@@ -136,9 +136,9 @@
         M2_M1_M1_M2_star = dot(M2_M1_star, M1_M2_star)
 
         # Matrices for even n.
-        if power[i] % 2 == 0:
+        if n[i] % 2 == 0:
             # The power factor (only calculate once).
-            fact = int(power[i] / 2)
+            fact = int(n[i] / 2)
 
             # (M1.M2.M2.M1)^(n/2)
             A = square_matrix_power(M1_M2_M2_M1, fact)
@@ -155,7 +155,7 @@
         # Matrices for odd n.
         else:
             # The power factor (only calculate once).
-            fact = int((power[i] - 1) / 2)
+            fact = int((n[i] - 1) / 2)
 
             # (M1.M2.M2.M1)^((n-1)/2).M1.M2
             A = square_matrix_power(M1_M2_M2_M1, fact)

Modified: branches/relax_disp/target_functions/relax_disp.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/relax_disp/target_functions/relax_disp.py?rev=21223&r1=21222&r2=21223&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Thu Oct 24 17:14:07 
2013
@@ -217,6 +217,10 @@
             self.power = zeros(self.num_disp_points, int16)
             for i in range(self.num_disp_points):
                 self.power[i] = int(round(self.cpmg_frqs[i] * 
self.relax_time))
+
+        # The strange n definition of Korzhnev.
+        if model == MODEL_MQ_NS_CPMG_2SITE:
+            self.n = 2 * self.power
 
         # Convert the spin-lock data to rad.s^-1.
         if spin_lock_nu1 != None:
@@ -896,7 +900,7 @@
                 dwH_frq = dwH[spin_index] * self.frqs[spin_index, frq_index]
 
                 # Back calculate the R2eff values.
-                r2eff_mq_ns_cpmg_2site(M0=self.M0, m1=self.m1, m2=self.m2, 
r20=R20[r20_index], pA=pA, pB=pB, dw=dw_frq, dwH=dwH_frq, k_AB=k_AB, 
k_BA=k_BA, inv_tcpmg=self.inv_relax_time, tcp=self.tau_cpmg, 
back_calc=self.back_calc[spin_index, frq_index], 
num_points=self.num_disp_points, power=self.power)
+                r2eff_mq_ns_cpmg_2site(M0=self.M0, m1=self.m1, m2=self.m2, 
r20=R20[r20_index], pA=pA, pB=pB, dw=dw_frq, dwH=dwH_frq, k_AB=k_AB, 
k_BA=k_BA, inv_tcpmg=self.inv_relax_time, tcp=self.tau_cpmg, 
back_calc=self.back_calc[spin_index, frq_index], 
num_points=self.num_disp_points, n=self.n)
 
                 # For all missing data points, set the back-calculated value 
to the measured values so that it has no effect on the chi-squared value.
                 for point_index in range(self.num_disp_points):