Author: bugman
Date: Wed Jul 31 16:56:37 2013
New Revision: 20490
URL: http://svn.gna.org/viewcvs/relax?rev=20490&view=rev
Log:
Converted the lib.dispersion.tp02 module from Matlab code to Python.
The code has also been made fail-safe and repetitive calculations have been
shifted outside of the
loop to speed things up.
Modified:
branches/relax_disp/lib/dispersion/tp02.py
Modified: branches/relax_disp/lib/dispersion/tp02.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/lib/dispersion/tp02.py?rev=20490&r1=20489&r2=20490&view=diff
==============================================================================
--- branches/relax_disp/lib/dispersion/tp02.py (original)
+++ branches/relax_disp/lib/dispersion/tp02.py Wed Jul 31 16:56:37 2013
@@ -32,52 +32,74 @@
- Martin Tollinger,
http://article.gmane.org/gmane.science.nmr.relax.devel/4276
"""
-function residual = funTrottPalmer(optpar)
-
-% TrottPalmer's equ. acc. to Korzhnev (JBiomolNMR, 26, 39-48, 2003)
-
-global nu_0 x y Rcalc rms nfields offset w1 R1_51 R1_81
-%keyboard
-Rcalc=zeros(nfields,size(w1,2));
-tau_ex=optpar(1);
-pb=optpar(2);
-pa=1-pb;
-kex=1/tau_ex;
-
-for k=1:nfields
- % we assume that A resonates at 0 [s^-1], without loss of generality
- dw=nu_0(k)*optpar(3)*2*pi; % [s^-1]
- Wa=0*2*pi; % Larmor frequ. [s^-1]
- Wb=dw; % Larmor frequ. [s^-1]
- Wsl=offset*2*pi; % Larmor frequ. of spin lock
[s^-1]
- W=pa*Wa+pb*Wb; % pop-averaged Larmor frequ. [s^-1]
- da=Wa-Wsl; % offset of sl from A
- db=Wb-Wsl; % offset of sl from B
- d =W -Wsl; % offset of sl from
pop-average
- waeff2=w1.^2+da.^2; % effective field at A
- wbeff2=w1.^2+db.^2; % effective field at B
- weff2=w1.^2+d.^2; % effective field at pop-average
- theta=atan(w1/d);
- theta_degrees=atan(w1/d)*360/(2*pi); % not used, just to check
-
- Rinf=optpar(3+k);
- xx=pa*pb*dw^2*kex;
- yytrottpalmer=waeff2.*wbeff2./weff2+kex^2;
-
yytrottpalmer_extended=waeff2.*wbeff2./weff2+kex^2-2*(sin(theta).^2)*pa*pb*dw^2;
-
- %keyboard
- if k==1;
-
Rcalc(k,:)=(cos(theta).^2)*R1_51+(sin(theta).^2)*(Rinf)+(sin(theta).^2).*xx./yytrottpalmer_extended;
- end
- if k==2;
-
Rcalc(k,:)=(cos(theta).^2)*R1_81+(sin(theta).^2)*(Rinf)+(sin(theta).^2).*xx./yytrottpalmer_extended;
- end
+# Python module imports.
+from math import atan, cos, pi, sin
-end
+def r1rho_TP02(r1rho_prime=None, pA=None, pB=None, dw=None, kex=None,
theta=pi/2, R1=0.0, spin_lock_fields=None, back_calc=None, num_points=None):
+ """Calculate the R1rho' values for the TP02 model.
-if (size(Rcalc)==size(y))
- residual=sum(sum((y-Rcalc).^2));
- rms=sqrt(residual/(size(y,1)*size(y,2)));
-end
+ See the module docstring for details. This is the Trott and Palmer
(2002) equation according to Korzhnev (J. Biomol. NMR (2003), 26, 39-48).
+
+ @keyword r1rho_prime: The R1rho_prime parameter value (R1rho with
no exchange).
+ @type r1rho_prime: float
+ @keyword pA: The population of state A.
+ @type pA: float
+ @keyword pB: The population of state B.
+ @type pB: float
+ @keyword dw: The chemical exchange difference between
states A and B in rad/s.
+ @type dw: float
+ @keyword kex: The kex parameter value (the exchange rate
in rad/s).
+ @type kex: float
+ @keyword theta: The rotating frame tilt angle.
+ @type theta: float
+ @keyword R1: The R1 relaxation rate.
+ @type R1: float
+ @keyword spin_lock_fields: The CPMG nu1 frequencies.
+ @type spin_lock_fields: numpy rank-1 float array
+ @keyword back_calc: The array for holding the back calculated
R1rho values. Each element corresponds to one of the spin-lock fields.
+ @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 spin_lock_fields and back_calc arguments.
+ @type num_points: int
+ """
+
+ # Repetitive calculations (to speed up calculations).
+ Wa = 0.0 # Larmor frequency [s^-1].
+ Wb = dw # Larmor frequency [s^-1].
+ sin_theta2 = sin(theta)**2
+ R1_cos_theta2 = R1 * cos(theta)**2
+ R1rho_prime_sin_theta2 = r1rho_prime * sin_theta2
+
+ # The numerator.
+ numer = pA * pB * dw**2 * kex
+
+ # Loop over the dispersion points, back calculating the R1rho values.
+ for i in range(num_points):
+ # Catch zeros (to avoid pointless mathematical operations).
+ if numer == 0.0:
+ back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2
+ continue
+
+ # We assume that A resonates at 0 [s^-1], without loss of generality.
+ Wsl = spin_lock_fields[i] * 2.0 * pi # Larmor frequency of spin
lock [s^-1].
+ W = pA*Wa + pB*Wb # Pop-averaged Larmor frequency [s^-1].
+ da = Wa - Wsl # Offset of spin-lock from A.
+ db = Wb - Wsl # Offset of spin-lock from B.
+ d = W - Wsl # Offset of spin-lock from pop-average.
+ waeff2 = w1**2 + da**2 # Effective field at A.
+ wbeff2 = w1**2 + db**2 # Effective field at B.
+ weff2 = w1**2 + d**2 # Effective field at pop-average.
+ theta = atan(w1/d)
+
+ # Denominator.
+ denom = waeff2 * wbeff2 / weff2 + kex**2
+ #denom_extended = waeff2*wbeff2/weff2+kex**2-2*sin_theta2*pA*pB*dw**2
+
+ # Avoid divide by zero.
+ if denom == 0.0:
+ back_calc[i] = 1e100
+ continue
+
+ # R1rho calculation.
+ back_calc[i] = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 *
numer / denom
r20490 - /branches/relax_disp/lib/dispersion/tp02.py