# Implementation of Allard et al.'s homogeneous master equation 
# extended for arbritary numbers of spins
# runs on GNU/Octave - not meant for matlab

clear 

#System

#####################################################################################	
# Edit the following:
#####################################################################################	
	pw=1*10^(-5) #duration of RF pulse - a pi/2 pulse is ca 1*10^(-5) (10 us) for omegax=1.5*10^5
	d=[100*10^(-1)] #duration of acquisition 
	si=8*1024; #number of points in spectrum
	pl=0; #power level attenuation

	C0=[10 10 10 10]; #relative concentrations
	T2=[0.16 0.16 0.16 0.16]; #T2 
	T1=[0.24 0.24 0.24 0.24]; # T1
	Omega=[250*2*pi 150*2*pi 0*2*pi -150*2*pi]; #center frequencies in Hz times 2*pi
	Omegax=[1.5*10^5 1.5*10^5 1.5*10^(5) 1.5*10^5]./10^(pl); # omega1 x component i.e. 'power'
	Omegay=[0 0 0 0]; # omega1, y component. Only use to change phase of pulse
#	k format for 4 spins in equilibrium:
#	[0 k_ab k_ac k_ad; 0 0 k_bc k_bd; 0 0 0 k_cd; 0 0 0 0]		
	k=[0 0.1 10.1 250; 0 0 100 0.1; 0 0 0 0.1; 0 0 0 0] # example
#	k=[0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0]  # no exchange
	
#####################################################################################	
	
	for n=1:columns(C0)
		k(:,n)=k(n,:)'.*C0(n)./C0';
	endfor
	k
	k=kron(k,eye(3,3));
	
	for n=1:3*columns(C0)
		k(n,n)=-sum(k(:,n));
	endfor
	k
	params=[C0;T2;T1;-Omega;Omegax;Omegay]
	spins=columns(C0);
#pulse

	HME=hmegen(params,k)
	testfunc= @(X,t) HME*X;

	t=linspace(0,pw,1*1024);

	x0=[1/2]
	for n=1:spins
		x0=[x0 0 0 C0(n) ];
	endfor
	x0=x0';

	HMESOL1=lsode(testfunc,x0,t);

#observe
	t=linspace(0,d(1),si);
	x0=HMESOL1(rows(HMESOL1),:)'

	newparams=[]

	for n=1:spins
			newparams=[newparams; params(1:4,n)' zeros(1,2)]
	endfor
	params=newparams';
	size(params),size(newparams')	
	disp('here')
#	params=[params(1:4,1)' zeros(1,columns(C0)); params(1:4,2)' zeros(1,columns(C0)); params(1:4,3)' zeros(1,columns(C0))]'

	HME2=hmegen(params,k)
	testfunc= @(X,t) HME2*X;
	HMESOL=lsode(testfunc,x0,t);

	a=0;
for n=1:spins
	a=a+HMESOL(:,(n-1)*3+2)+i*HMESOL(:,(n-1)*3+3);
endfor

spec=[real(fftshift( fft(a) )) imag(fftshift(fft(a)))];
freq=linspace(-(si+1)/(2*d(1)),(si-1)/(2*d(1)),si);
plot(freq,spec(:,2))
