Two dimensional NMR simulations are also a snap to perform using GAMMA. For example, consider the above E-COSY pulse sequence which is used to ascertain scalar coupling values of "passive" spins. The proper phase cycle is essential to this sequence.

GAMMA provides pulses, delays, and acquisition functions which the programmer may readily assemble in a program. Coupled with GAMMA's interface to plotting routines and access to a full programming language (e.g. looping over the phase cycle) it is relatively easy to construct programs for e-cosy and a wide variety of other pulse sequences. Such programs can be "black boxes" which will simulate the spectrum when given a spin system and set of acquisition parameters. Below are output spectra from a couple of input systems fed into the same program.

Shown are results from some simple 3-spin proton systems, but the program itself is not limited protons nor to 3 spins. It is general for any input system and the user can specify all acquisition parameters. It can even produce the above figures (and sub-section plots) directly.

Due to the inherent requirement for robust phase cycling, e-cosy simulation programs tend to be a bit much for the beginning GAMMA user. If this is too long to follow just have a look at some simpler COSY like programs and peruse the example documentation which should have many such examples.

#include <gamma.h>

// 			   Define Constants

  const P_cycl  = 12;				// Phase cycle length
  const double P_mix[12] = 			// Pulse phase cycle
	{0,60,0,60,0,60,0,300,120,300,240,300};

//-------------------------------------------------------------------
// 			   Begin Program
//-------------------------------------------------------------------

main (int argc, char* argv[])

  {
  cout << "\n\tSimulation of an E-COSY spectrum in TPPI mode\n";
  cout << "\t (Hard Pulses, No Relaxation, MATLAB Output)\n";

//          	   Read in Spectral Parameters

  int qn = 1;					// Query number
  spin_system sys;				// A spin system
  sys.ask_read(argc, argv, qn++);		// Request/Read system 
  int t1pts, t2pts;
  query_parameter(argc, argv, qn++,		// Get number FID of points
              "\n\tAcquisition Size? ", t2pts);
  t1pts = 2*t2pts;				// Set t1 size for TPPI
  string Isotype = sys.symbol(0);		// Get isotope type of first spin
  int isoset = query_isotope(sys, Isotype);     // Choose an isotope type
  query_offset(sys, isoset, 1); 		// Ask for an offset frequency

// 		Set Up the Isotropic Hamiltonian

  string J;
  query_parameter(argc, argv, qn++,		// Weak or strong coupling
      "\n\tWeak or strong coupling (w/s)? ", J);
  gen_op H;					// Set Hamilitonian for
  if (J == "w")	H = Hcs(sys) + HJw(sys); 	// strong or weak coupling
  else          H = Hcs(sys) + HJ(sys);

//       Set Up the Relaxation Superoperator, Liouvillian

  super_op L = complexi*Hsuper(H);		// L = -i*[Ho, ] (rad/sec)
  double lwhh = 0.5;				// Expected dipolar linewidth
  cout << sys;					// See system before Nyquist
  double NyqF = query_Nyquist(sys,isoset,lwhh); // Choose a Nyquist frequency
  double t2dt = 1.0/(2.0*NyqF);			// Dwell time, quadrature
  double t1dt = t2dt/2.0;			// t1 time increment, TPPI

//       Set Up Other Necessary Operators, Superoperators

  super_op eLt1 = exp(L, -t1dt);		// Delay t1 relaxation superoperator
  gen_op sigma0 = sigma_eq(sys);		// Density matrix at equilibrium
  set_trace(sigma0, 1.0);
  super_op G_d1 = R_prop(eLt1, sigma0);		// Delay t1 relaxation propagator
  gen_op RTPPI = Rz(sys,90);			// Rotation for TPPI phasing
  RTPPI.Op_base(H);				// Put operator in Ho eigenbasis
  G_d1 = U_transform(RTPPI)*G_d1;		// Delay t1 + TPPI propagator
  gen_op D = Fm(sys);				// Detector to F-
  gen_op Upx = Ixypuls_U(sys, 0, +90);		// Propagator for x pulse
  gen_op Upmx = Ixypuls_U(sys, 0, -90);		// Propagator for -x pulse
  acquire ACQ(D, L, t2dt);			// Prepare for acquisitions
  gen_op U_mix;					// Temporary mixing propagator
  super_op G_mix;				// Phase cycle superoperator
  gen_op sigma3;

//		Construct Mixing, Phase Cycle Superoperator
//
  double conv = acos(-1)/180;
  double P_det = 0;
  cout << "\n\tConstructing Phase Cycle Superoperator\n";
  cout << "\n\t# scan	beta	reference\n";
  for(int i=0; i<P_cycl; i++ )
    {
    cout <<"\t" << i << "	"
         << P_mix[i] << "	"
         << P_det    <<"\n";
    U_mix  = Rz(sys,-P_det);			// Detector phase cycle
    U_mix *= Upmx;				// 3rd 90 Pulse (-x)
    U_mix *= gen_op(Rz(sys,+P_mix[i]));		// Phase  pulses 1 & 2
    U_mix *= Upx;				// 2nd 90 Pulse
    U_mix.Op_base(H);				// Put in eigenbasis of Ho
    G_mix += U_transform(U_mix);		// Add to U transform superoperator
    P_det = acos(-cos(P_det*conv))/conv;	// Adjust detector phase
    }

//			Apply Pulse Sequence

  row_vector blk(t2pts);			// Set 1D block for output
  matrix data(t1pts,t2pts);			// Array for output
  gen_op sigma1 = evolve(sigma0, Upx);		// Apply first (PI/2)x pulse
  gen_op sigma2 = sigma1;			// Initial sigma2 (t1 = 0)
  for(int t1=0; t1w2)
    data.put_block(t1,0,blk);			// 	Store FID in data matrix
    evolve_ip(sigma2, G_d1);			//	Evolution next t1
    }


  int t1ptsA = t1pts/2; 			// We'll only get 1/2 points in t1-w1
  row_vector shtblk(t1ptsA);			// A vector for those points
  complex z, zph;				// Some scratch complex numbers
  zph = exp(complexi*0.5*PI);			// For 90 zero order phase correction
  int j = t1pts-1;				// Required index
  for(int t2=0; t2<t2pts; t2++)		// Loop over all t2 increments
    {
    blk=transpose(data.get_block(0,t2,t1pts,1));//	Clip t1pts column at point t2
    exponential_multiply(blk, -8);		//	Apodize interferrogram in t1
    for(int i=0; iw1)
    for(i=0, j=t1pts-1; i<t1ptsA; i++, j--)	//	Only 1/2 the points are good
      { 					//	Reform complex point following
      z = complex(blk.getRe(i),blk.getRe(j));	//	the real FFT.  The restore it
      shtblk.put(z*zph,i);			//	in short block, with phase applied
      }
    data.put_block(0, t2, transpose(shtblk));	//	Replace the column
    }
  data = data.get_block(0,0,t1ptsA,t2pts);	// Only 1/2 of data is good
  MATLAB("ecosy.mat", "f1f2mx", data);		// Output array to MATLAB
  }

// These are the MATLAB processing commands that are of consequence
// (I assume npts on t2 was 512, replace 512 below with your value)
// Note that this processing demands that the spectrum point order
// is adjusted (split, flipped, rejoined) te get the proper frequency
// ordering due to the nature of MATLAB's FFT algorithm.

//	load ecosy.mat			Load the output file 	(1024, 512)
//	t1f2mx = real(t1f2mx);		Use only the reals	(1024, 512)
//	f1f2mx = fft(t1f2mx);		FFT columns (t1->f1) 	(1024, 512)
//      f1f2 = f1f2mx(1:512,1:512);	Clip out the 2nd half   (512, 512)
//	f1f2 = imag(f1f2);		Use only imaginaries	(512, 512)
//	f1f2 = flipud(f1f2);		Reverse the row order	(512, 512)
//      contour(f1f2, 5);		Plot contours

// There should be an ecosyML0.m file around (that came with this program)
// which contains the MATLAB commands for generating a contour plot.  That
// can be run directly.  In MATLAB, assuming ecosy.mat and ecosyML0.m are
// found by MATLAB, the single command ecosyML0 should produce the plot.