# define parameters
	%Machine
	omegarf=-424288800;
	omega1=0.005*abs(omegarf);

	%Acquisition parameters
	pipulse=0.7*10^(-6);
	pw=1*pipulse;
	delay=0.025;
	nopts=1*1024;
	sw=1/(delay/nopts);

	%Processing parameters
	freqaxis=linspace(-sw/2,sw/2,nopts);
	
	%species a
	omega0a=omegarf+2*pi*10000;
	T2a= 10^(-3); T1a=10^(-3);
	mza0=5; 

	%species b
	omega0b=omegarf-2*pi*5000;
	T2b= 10^(-3); T1b=10^(-3);
	mzb0=10; 

	%kinetic parameters
	k=1*10^(1); pab=mza0/mzb0;
	
	# bake parameters into arrays for easy handling
	omega=[omegarf,omega1,omega0a,omega0b];
	T2=[T2a,T2b];
	T1=[T1a,T1b];
	mz=[mza0,mzb0];
##################

%pulse
	t=linspace(0,pw,nopts);
	x0=[0,0,mz(1),0,0,mz(2)];
	y=getpulse(x0,omega,T1,T2,mz,k,t);
	pulse=[t' y];

%observe
	t=linspace(0,delay,nopts);
	x0=[y(nopts,1),y(nopts,2),y(nopts,3),y(nopts,4),y(nopts,5),y(nopts,6)];
	#turn off rf power
	omega(2)=0;
	#simulate acquisition
	y=getpulse(x0,omega,T1,T2,mz,k,t);
	observe=[t' y];

#Process spectrum
	spectrum=imag(fftshift(fft(observe(:,2)+i*observe(:,3)+observe(:,5)+i*observe(:,6))));
	
#Plot spectrum
subplot(3,1,1)
	plot(observe(:,1),pulse(:,2:7))
subplot(3,1,3)
	plot(freqaxis,spectrum)
subplot(3,2,3)
	plot(observe(:,1),observe(:,2)+observe(:,5),observe(:,1),observe(:,3)+observe(:,6))
subplot(3,2,4)
	plot(observe(:,1),observe(:,2)+observe(:,5),observe(:,1),observe(:,3)+observe(:,6))

#abd=[(freqaxis./1000)' spectrum];

