% propagation under a static Hamiltonian d=750; % dipolar coupling sw=1600; % spectral width sys=[1/2 1/2]; % two spin-1/2 spin system tsteps=64; % points in FID % Hamiltonian H=d*(2*S(sys,1,'z')*S(sys,2,'z')-S(sys,1,'x')*S(sys,2,'x')-S(sys,1,'y')*S(sys,2,'y')); sigma0=S(sys,1,'x'); % initial density matrix detect=F(sys,'x'); % detection operator dwell=1/sw; sigma=sigma0; % initial density matrix U=expm(-i*2*pi*dwell*H); % dwell time propagator FID=zeros(1,tsteps); % data array for j=1:tsteps % propagation loop FID(j)=trace(detect*sigma); % detection sigma=U*sigma*U'; % propagation end FID1=FID; % EXAMPLE 1 [V,evals]=eig(H); % eigenvectors and eigenvalues phases=exp(-2*pi*i*diag(evals)*dwell); % convert eigenvalues to a vector in angular units detectH=V'*detect*V; % detection operator in H eigenbasis sigma0H=V'*sigma0*V; % initial density matrix in H eigenbasis A= sigma0H.* (detectH.'); % amplitudes matrix FID=zeros(1,tsteps); n=length(evals); % dimensionality of Hilbert space FID=FID+trace(A); % add constant terms for j=1:n for k=j+1:n % loop over \$k>j\$ transitions v=2*A(j,k); % signal amplitude if abs(v)>1e-10 % ignore transitions of negligible amplitude p=phases(j)*conj(phases(k)); % phase factor for t=1:tsteps % loop over dwell times FID(t)=FID(t)+real(v); % accummulate FID v=v*p; % time evolution of signal end end end end FID2=FID; tscale=[0:tsteps-1]*dwell*1000; plot(tscale,real(FID1),'r-',tscale,real(FID2),'b*'); xlabel('time / ms');