% File: NMRHourteenPulse.m
%
% Date:   05-Mar-02
% Author: M. Steffen <msteffen@snowmass.stanford.edu
% 
% usage: NMRCalib(pw90,expflag)
% 
% where: pw90  -  array containing the proton and carbon pulse widths in
%                 microseconds (eg. [14.25 9])
%        expflag - 0 does experiment plus analysis (DEFAULT VALUE)
%		   1 does everything except for experiment
%		   (useful if analysis screws up)       
 

function spect = NMRCalib(pw90,expflag) 

     expname=100;  % SAVES DATA IN EXP100 !!!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check for correct array size and inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(nargin~=2)
     expflag=0;
end
if(nargin<1)
     error('NMRCalib requires at least 1 input!');
end
if(length(pw90)~=2)
     error('pw90 needs to be an array of 2 values!');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET PULSE SEQUENCE AND DEFINE SOME PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

seqfil='''juniorcalib''';
d1=100;
rcvrs='''yynn''';
sw=400;
mrfb=1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vjexp(6);
vsetparam('seqfil',seqfil,'%s');

vsetparam('d1',d1,'%d');
vsetparam('rcvrs',rcvrs,'%s');
vsetparam('sw',sw,'%d');
vsetparam('mrfb',mrfb,'%d');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN EXPERIMENT AND GET DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(expflag==0)
fprintf(' %%%% starting experiment now. \n');
unix('sleep 1');
vgawait(6);
s=sprintf('vgetdatnoplot(6,%d)',expname);
eval(s);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ IN DATA AND FFT USING V2SDAT2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=sprintf('sdp%d=v2sdat2j(%d,1,0)',expname,expname); eval(s);
s=sprintf('sdc%d=v2sdat2j(%d,2,0)',expname,expname); eval(s);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate peaks and phase thermal spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=sprintf('peakfrq_p%d = peaksj(sdp%d,0);',expname,expname); eval(s);
s=sprintf('peakfrq_c%d = peaksj(sdc%d,0);',expname,expname); eval(s);

s=sprintf('peaks_p%d=pintj(sdp%d,peakfrq_p%d,10);',expname,expname,expname); eval(s);
s=sprintf('phase_p%d=sum(phase(peaks_p%d))/2',expname,expname); eval(s);
s=sprintf('peaks_c%d=pintj(sdc%d,peakfrq_c%d,10);',expname,expname,expname); eval(s);
s=sprintf('phase_c%d=sum(phase(peaks_c%d))/2',expname,expname); eval(s);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ASSIGN SPECT STRUCTURE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=sprintf('spect.hfreq =  sdp%d.fdat;',expname); eval(s);
s=sprintf('spect.cfreq =  sdc%d.fdat;',expname); eval(s);
s=sprintf('spect.hspect = sdp%d.spect*exp(-i*phase_p%d);',expname, expname); eval(s);
s=sprintf('spect.cspect = sdc%d.spect*exp(-i*phase_c%d);',expname, expname); eval(s);
s=sprintf('spect.hsfo =   sdp%d.sfo;',expname); eval(s);
s=sprintf('spect.csfo =   sdc%d.dfo;',expname); eval(s);
s=sprintf('spect.dt =     sdc%d.dt;',expname); eval(s);
s=sprintf('spect.hpeaks = peaks_p%d*exp(-i*phase_p%d);',expname,expname); eval(s);
s=sprintf('spect.cpeaks = peaks_c%d*exp(-i*phase_c%d);',expname,expname); eval(s);
s=sprintf('spect.hphase = phase_p%d;',expname); eval(s);
s=sprintf('spect.cphase = phase_c%d;',expname); eval(s);
s=sprintf('spect.pp =     sdc%d.pp;',expname); eval(s);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT SPECTRA TO MAKE SURE PHASING WORKED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
plot(spect.hfreq,real(spect.hspect));
title('proton spectrum');
ylabel('signal [arb. units]');
xlabel('frequency from hsfo');

figure(2)
plot(spect.cfreq,real(spect.cspect));
title('carbon spectrum');
ylabel('signal [arb. units]');
xlabel('frequency from csfo');
