(* Burkhard Luy, Oliver Schedletzkyi, and Steffen J. Glaser, Analytical Polarization Transfer Functions for Four Coupled Spins 1/2 under Isotropic Mixing Conditions Journal of Magnetic Resonance 138, 19-27 (1999). A Hamiltonian of the form H=2 pi J (I1xI2x + I1yI2y + I1zI2z) is assumed (expanded to 4 spins, of course). The scalar couplings between spins 1 and 2 is called j12 etc.. The coupling constants have to be defined right after this comment. The polarization transfer I1z -> I2z is stored in tr12 and the transfer I1z -> I2z is stored in tr11. Please cite the above mentioned reference whenever you use this program for scientific work. *) (*Kopplungskonstanten *) J12=10 J13=-20 J14=0 J23=-20 J24=0 J34=0 (* J12=4.1 J13=9.4 J14=6.8 J23=5.3 J24=8.2 J34=-4.6 *) (* Vorstufen *) sigma=(J12+J34)/4 delta=(J12-J34)/4 A=(J13+J14+J23+J24)/4 B=(J13-J14+J23-J24)/4 H=(-J13-J14+J23+J24)/4 K=(J13-J14-J23+J24)/4 R=Sqrt[(A-2*sigma)^2+6*B^2+12*delta^2+6*H^2+3*K^2]/3 W=Sqrt[A^2+3*K^2+4*sigma(sigma-A)] Q=(A^3 + 9*A*B^2 - 54*B^2*delta - 36*A*delta^2 + 9*A*H^2 + 54*delta*H^2 - 54*B*H*K - 9*A*K^2 - 6*A^2*sigma - 18*B^2*sigma + 72*delta^2*sigma - 18*H^2*sigma + 18*K^2*sigma + 12*A*sigma^2 - 8*sigma^3)/27 (* Winkel *) phi=ArcCos[-Q/R^3] (* If[(2*sigma-A-W)==0,If[K>0,psi=Pi/2,psi=-Pi/2], psi=ArcTan[(2*sigma-A+W)/(Sqrt[3]*K)] ] *) If[(2*sigma-A+W)==0,If[K>0,psi=Pi/2,psi=-Pi/2], psi=ArcTan[(Sqrt[3]*K)/(2*sigma-A+W)] ] (* lambda,mue,nues *) lambda=2*Pi*(sigma+A) mue1=-2*Pi*((sigma+A)/3-2*R*Cos[phi/3]) mue2=-2*Pi*((sigma+A)/3+2*R*Cos[(phi-Pi)/3]) mue3=-2*Pi*((sigma+A)/3+2*R*Cos[(phi+Pi)/3]) nu1=-2*Pi*(sigma+A+W) nu2=-2*Pi*(sigma+A-W) (* Zur Eigenvektorbestimmung *) mu1=-((sigma+A)/3-2*R*Cos[phi/3]) mu2=-((sigma+A)/3+2*R*Cos[(phi-Pi)/3]) mu3=-((sigma+A)/3+2*R*Cos[(phi+Pi)/3]) (* Vorstufen Yu Vorfaktoren *) v1=Cos[psi] u2=Cos[psi] u1=Sin[psi] v2=-Sin[psi] alpha1=Sqrt[2]*(-2*delta-sigma)*B-Sqrt[2]*H*K-Sqrt[2]*B*mu1 alpha2=Sqrt[2]*(-2*delta-sigma)*B-Sqrt[2]*H*K-Sqrt[2]*B*mu2 alpha3=Sqrt[2]*(-2*delta-sigma)*B-Sqrt[2]*H*K-Sqrt[2]*B*mu3 beta1=-2*B*H+(sigma-A)*K-K*mu1 beta2=-2*B*H+(sigma-A)*K-K*mu2 beta3=-2*B*H+(sigma-A)*K-K*mu3 gamma1=2*H^2-(-2*delta-sigma)*(sigma-A)+(-2*delta-A)*mu1-mu1^2 gamma2=2*H^2-(-2*delta-sigma)*(sigma-A)+(-2*delta-A)*mu2-mu2^2 gamma3=2*H^2-(-2*delta-sigma)*(sigma-A)+(-2*delta-A)*mu3-mu3^2 alp1=K*(K-Sqrt[2]*H-Sqrt[2]*B)+4*delta^2-sigma^2+Sqrt[8]*delta*(H- B)-Sqrt[2]*sigma*(H+B)-(2*sigma+Sqrt[2]*H+Sqrt[2]*B)*mu1-mu1^2 alp2=K*(K-Sqrt[2]*H-Sqrt[2]*B)+4*delta^2-sigma^2+Sqrt[8]*delta*(H- B)-Sqrt[2]*sigma*(H+B)-(2*sigma+Sqrt[2]*H+Sqrt[2]*B)*mu2-mu2^2 alp3=K*(K-Sqrt[2]*H-Sqrt[2]*B)+4*delta^2-sigma^2+Sqrt[8]*delta*(H- B)-Sqrt[2]*sigma*(H+B)-(2*sigma+Sqrt[2]*H+Sqrt[2]*B)*mu3-mu3^2 bet1=Sqrt[2]*B*(Sqrt[2]*B-Sqrt[2]*H-K)+sigma^2-2*delta*sigma+ 2*delta*A-sigma*A+Sqrt[2]*H*(2*delta-sigma)+K*(sigma-A)+ (2*delta-A-K-Sqrt[2]*H)*mu1-mu1^2 bet2=Sqrt[2]*B*(Sqrt[2]*B-Sqrt[2]*H-K)+sigma^2-2*delta*sigma+ 2*delta*A-sigma*A+Sqrt[2]*H*(2*delta-sigma)+K*(sigma-A)+ (2*delta-A-K-Sqrt[2]*H)*mu2-mu2^2 bet3=Sqrt[2]*B*(Sqrt[2]*B-Sqrt[2]*H-K)+sigma^2-2*delta*sigma+ 2*delta*A-sigma*A+Sqrt[2]*H*(2*delta-sigma)+K*(sigma-A)+ (2*delta-A-K-Sqrt[2]*H)*mu3-mu3^2 gam1=Sqrt[2]H(Sqrt[2]*H-Sqrt[2]*B-K)+(sigma-A)*(2*delta+sigma+K)- Sqrt[2]*B*(2*delta+sigma)-(2*delta+A+K+Sqrt[2]*B)*mu1-mu1^2 gam2=Sqrt[2]H(Sqrt[2]*H-Sqrt[2]*B-K)+(sigma-A)*(2*delta+sigma+K)- Sqrt[2]*B*(2*delta+sigma)-(2*delta+A+K+Sqrt[2]*B)*mu2-mu2^2 gam3=Sqrt[2]H(Sqrt[2]*H-Sqrt[2]*B-K)+(sigma-A)*(2*delta+sigma+K)- Sqrt[2]*B*(2*delta+sigma)-(2*delta+A+K+Sqrt[2]*B)*mu3-mu3^2 mag1=Sqrt[alp1^2+bet1^2+gam1^2] mag2=Sqrt[alp2^2+bet2^2+gam2^2] mag3=Sqrt[alp3^2+bet3^2+gam3^2] al1=alp1/mag1 al2=alp2/mag2 al3=alp3/mag3 be1=bet1/mag1 be2=bet2/mag2 be3=bet3/mag3 ga1=gam1/mag1 ga2=gam2/mag2 ga3=gam3/mag3 magn1=Sqrt[alpha1^2+beta1^2+gamma1^2] magn2=Sqrt[alpha2^2+beta2^2+gamma2^2] magn3=Sqrt[alpha3^2+beta3^2+gamma3^2] tx1=alpha1/magn1 tx2=alpha2/magn2 tx3=alpha3/magn3 tY1=beta1/magn1 tY2=beta2/magn2 tY3=beta3/magn3 tZ1=gamma1/magn1 tZ2=gamma2/magn2 tZ3=gamma3/magn3 x1=al1 x2=al2 x3=al3 Y1=be1 Y2=be2 Y3=be3 Z1=ga1 Z2=ga2 Z3=ga3 (* die Vorfaktoren *) a1=-5*(x1^2-2*Y1^2)/12 a2=-5*(x2^2-2*Y2^2)/12 a3=-5*(x3^2-2*Y3^2)/12 b12=-(x1*x2(Z1*Z2-Y1*Y2)-(x1^2*Y2^2+x2^2*Y1^2)/2+(x1^2*x2^2)/4+Z1^2*Z2^2) b13=-(x1*x3(Z1*Z3-Y1*Y3)-(x1^2*Y3^2+x3^2*Y1^2)/2+(x1^2*x3^2)/4+Z1^2*Z3^2 ) b23=-(x2*x3(Z2*Z3-Y2*Y3)-(x2^2*Y3^2+x3^2*Y2^2)/2+(x2^2*x3^2)/4+Z2^2*Z3^2 ) c11=-((u1^2*(2*x1^2-Y1^2))/6+Z1*v1*((Y1*u1)/Sqrt[3]-(Z1*v1)/2)) c12=-((u2^2*(2*x1^2-Y1^2))/6+Z1*v2*((Y1*u2)/Sqrt[3]-(Z1*v2)/2)) c21=-((u1^2*(2*x2^2-Y2^2))/6+Z2*v1*((Y2*u1)/Sqrt[3]-(Z2*v1)/2)) c22=-((u2^2*(2*x2^2-Y2^2))/6+Z2*v2*((Y2*u2)/Sqrt[3]-(Z2*v2)/2)) c31=-((u1^2*(2*x3^2-Y3^2))/6+Z3*v1*((Y3*u1)/Sqrt[3]-(Z3*v1)/2)) c32=-((u2^2*(2*x3^2-Y3^2))/6+Z3*v2*((Y3*u2)/Sqrt[3]-(Z3*v2)/2)) (* aa1 = (5*x1^2)/12 + (5*x1*z1)/(3*2^(1/2)) + (5*z1^2)/6; aa2 = (5*x2^2)/12 + (5*x2*z2)/(3*2^(1/2)) + (5*z2^2)/6; aa3 = (5*x3^2)/12 + (5*x3*z3)/(3*2^(1/2)) + (5*z3^2)/6; bb12 = (x1^2*x2^2)/4 + x1*x2*Z1*y2 + y1^2*y2^2 - (x1*x2^2*z1)/2^(1/2) - 2^(1/2)*x2*y1*y2*z1 + (x2^2*z1^2)/2 - (x1^2*x2*z2)/2^(1/2) - 2^(1/2)*x1*y1*y2*z2 + x1*x2*z1*z2 + (x1^2*z2^2)/2; bb13 = (x1^2*x3^2)/4 + x1*x3*y1*y3 + y1^2*y3^2 - (x1*x3^2*z1)/2^(1/2) - 2^(1/2)*x3*y1*y3*z1 + (x3^2*z1^2)/2 - (x1^2*x3*z3)/2^(1/2) - 2^(1/2)*x1*y1*y3*z3 + x1*x3*z1*z3 + (x1^2*z3^2)/2; bb23 = (x2^2*x3^2)/4 + x2*x3*y2*y3 + y2^2*y3^2 - (x2*x3^2*z2)/2^(1/2) - 2^(1/2)*x3*y2*y3*z2 + (x3^2*z2^2)/2 - (x2^2*x3*z3)/2^(1/2) - 2^(1/2)*x2*y2*y3*z3 + x2*x3*z2*z3 + (x2^2*z3^2)/2; xx1 = u1 xx2 = u2 yy1 = v1 yy2 = v2 cc11 = (xx1^2*x1^2)/3 + (2/3)^(1/2)*xx1*x1*yy1*y1 + (yy1^2*y1^2)/2 - (2^(1/2)*xx1^2*x1*Y1)/3 - (xx1*yy1*y1*Y1)/3^(1/2) + (xx1^2*Y1^2)/6; cc12 = (xx2^2*x1^2)/3 + (2/3)^(1/2)*xx2*x1*yy2*y1 + (yy2^2*y1^2)/2 - (2^(1/2)*xx2^2*x1*Y1)/3 - (xx2*yy2*y1*Y1)/3^(1/2) + (xx2^2*Y1^2)/6; cc21 = (xx1^2*x2^2)/3 + (2/3)^(1/2)*xx1*x2*yy1*y2 + (yy1^2*y2^2)/2 - (2^(1/2)*xx1^2*x2*Y2)/3 - (xx1*yy1*y2*Y2)/3^(1/2) + (xx1^2*Y2^2)/6; cc22 = (xx2^2*x2^2)/3 + (2/3)^(1/2)*xx2*x2*yy2*y2 + (yy2^2*y2^2)/2 - (2^(1/2)*xx2^2*x2*Y2)/3 - (xx2*yy2*y2*Y2)/3^(1/2) + (xx2^2*Y2^2)/6; cc31 = (xx1^2*x3^2)/3 + (2/3)^(1/2)*xx1*x3*yy1*y3 + (yy1^2*y3^2)/2 - (2^(1/2)*xx1^2*x3*Y3)/3 - (xx1*yy1*y3*Y3)/3^(1/2) + (xx1^2*Y3^2)/6; cc32 = (xx2^2*x3^2)/3 + (2/3)^(1/2)*xx2*x3*yy2*y3 + (yy2^2*y3^2)/2 - (2^(1/2)*xx2^2*x3*Y3)/3 - (xx2*yy2*y3*Y3)/3^(1/2) + (xx2^2*Y3^2)/6; *) aa1=5*((x1+Sqrt[2]*Y1)^2)/48 aa2=5*((x2+Sqrt[2]*Y2)^2)/48 aa3=5*((x3+Sqrt[2]*Y3)^2)/48 bb12=0.25*((x1*x2)/2 + (Z1*Z2) - (x1*Y2)/Sqrt[2] - (x2*Y1)/Sqrt[2])^2 bb13=0.25*((x1*x3)/2 + (Z1*Z3) - (x1*Y3)/Sqrt[2] - (x3*Y1)/Sqrt[2])^2 bb23=0.25*((x2*x3)/2 + (Z2*Z3) - (x2*Y3)/Sqrt[2] - (x3*Y2)/Sqrt[2])^2 cc11=0.25*((x1*u1)/Sqrt[3] + (Z1*v1)/Sqrt[2] - (Y1*u1)/Sqrt[6])^2 cc21=0.25*((x2*u1)/Sqrt[3] + (Z2*v1)/Sqrt[2] - (Y2*u1)/Sqrt[6])^2 cc31=0.25*((x3*u1)/Sqrt[3] + (Z3*v1)/Sqrt[2] - (Y3*u1)/Sqrt[6])^2 cc12=0.25*((x1*u2)/Sqrt[3] + (Z1*v2)/Sqrt[2] - (Y1*u2)/Sqrt[6])^2 cc22=0.25*((x2*u2)/Sqrt[3] + (Z2*v2)/Sqrt[2] - (Y2*u2)/Sqrt[6])^2 cc32=0.25*((x3*u2)/Sqrt[3] + (Z3*v2)/Sqrt[2] - (Y3*u2)/Sqrt[6])^2 (* Die Transferfunktion *) sum1=a1(1-Cos[(lambda-mue1)*tau])+a2(1-Cos[(lambda-mue2)*tau])+a3(1-Cos[(lambda-mue3)*tau]) sum2=b12(1-Cos[(mue1-mue2)*tau])+b13(1-Cos[(mue1-mue3)*tau])+b23(1-Cos[(mue2-mue3)*tau]) sum3=c11(1-Cos[(mue1-nu1)*tau])+c12(1-Cos[(mue1-nu2)*tau])+ c21(1-Cos[(mue2-nu1)*tau])+c22(1-Cos[(mue2-nu2)*tau])+ c31(1-Cos[(mue3-nu1)*tau])+c32(1-Cos[(mue3-nu2)*tau]) ssum1=aa1(1-Cos[(lambda-mue1)*tau])+aa2(1-Cos[(lambda-mue2)*tau])+aa3(1-Cos[(lambda-mue3)*tau]) ssum2=bb12(1-Cos[(mue1-mue2)*tau])+bb13(1-Cos[(mue1-mue3)*tau])+bb23(1-Cos[(mue2-mue3)*tau]) ssum3=cc11(1-Cos[(mue1-nu1)*tau])+cc12(1-Cos[(mue1-nu2)*tau])+ cc21(1-Cos[(mue2-nu1)*tau])+cc22(1-Cos[(mue2-nu2)*tau])+ cc31(1-Cos[(mue3-nu1)*tau])+cc32(1-Cos[(mue3-nu2)*tau]) transfer11=1-ssum1-ssum2-ssum3 transfer12=sum1+sum2+sum3 tr11=Simplify[N[transfer11]] tr12=Simplify[N[transfer12]]/4 Plot[tr11, {tau,0,0.400}] Plot[tr12, {tau,0,0.400}] (* Plot[a1(1-Cos[(lambda-mue1)*tau]), {tau,0,0.400}] Plot[a2(1-Cos[(lambda-mue1)*tau]), {tau,0,0.400}] Plot[a3(1-Cos[(lambda-mue1)*tau]), {tau,0,0.400}] Plot[b12(1-Cos[(mue2-mue1)*tau]), {tau,0,0.400}] Plot[b13(1-Cos[(mue3-mue1)*tau]), {tau,0,0.400}] Plot[b23(1-Cos[(mue3-mue2)*tau]), {tau,0,0.400}] Plot[c11(1-Cos[(nu1-mue1)*tau]), {tau,0,0.400}] Plot[c21(1-Cos[(nu1-mue2)*tau]), {tau,0,0.400}] Plot[c31(1-Cos[(nu1-mue3)*tau]), {tau,0,0.400}] Plot[c12(1-Cos[(nu2-mue1)*tau]), {tau,0,0.400}] Plot[c22(1-Cos[(nu2-mue2)*tau]), {tau,0,0.400}] Plot[c32(1-Cos[(nu2-mue3)*tau]), {tau,0,0.400}] *) (* Plot[N[ssum1], {tau,0,0.756}] Plot[N[ssum2], {tau,0,0.756}] Plot[N[ssum3], {tau,0,0.756}] *)