
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%                                                        &
%                        3-D SIFT                        &
%                                                        &
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% last update, mar1709, y.matsuki
% last update, may2812, y.matsuki
% last update, jun2612, m.t.eddy
% last update, dec2814, ym
% last update, jan3016, ym

%% General instructions: 
% This run script guides you through the SIFT processing. In each section,
% edit arguments for your data, and execute the commands in the section.
% The "section" is deliniated by the symbol "%%" and horizontal lines in
% the script.

% There are number of ways to execute the commands: 
% 1. Place the cursor in a "section", and press the "Run Section" button in
% the EDITOR tab at the top of this script viewer. An equivalent shortcut
% is "Apple+return".

% 2. "Evaluate Selection" command is also very useful for executing only
% the highlighted lines. You can set up a shortcut for the "Evaluate
% Selection" at Preferences/Keyboard/Shortcus.

% 3. Above is same to simply copy and paste the relevant lines into the
% "Command Window".


% Usually, fold all the if-blocks in this run script by "Apple+'='". ;use
% "Apple+shift+'='" to unfold all the if-blocks to make internal changes.

% If if-blocks don't fold, go to Preferences/Editor/Code Folding, and make
% if-blocks foldable.


%%
%\\\\\\\\\\\\\\\
% 0. Preparation
%\\\\\\\\\\\\\\\

% 0-1. 
% Convert your data to nmrPipe format, and determine phase correction value
% for the direct acquisition dimension.

% Refer to the section 0. of the instruction manual for more detail.





%%
%\\\\\\\\\\\\\\\
% 1. Set Parameters
%\\\\\\\\\\\\\\\

%% 1-1.
% Set the path to the SIFT script,
    cd ~/SIFT3DinMATLAB/SIFTsuite3D
    
%% 1-2.
% Read NUS table as
    load ~/SIFT3DinMATLAB/example/sched
    NUStab=sched;        
    
%% 1-3.
% Set other acquisition and processing parameters

    nt1     =12;        % # of full grids after SIFT reconstruction
    nt2     =12;        
    nt3     =1024;      

    nzf1    =nt1*8;     % after zero filling (power of two)
    nzf2    =nt2*8;     
    nzf3    =nt3;      

    cent1   =177;       % center(ppm) for F1 dimension
    cent2   =118;       % center(ppm) for F2 dimension
    cent3   =100;       % center(ppm) for F3 dimension

    obs1    =175.937;   % carrier(MHz) in F1 dimension
    obs2    =70.899;    % carrier(MHz) in F2 dimension
    obs3    =175.937;   % carrier(MHz) in F3 dimension

    sw1     =5278.716;  % SW(Hz) for F1 dimension
    sw2     =2836.075;  % SW(Hz) for F2 dimension
    sw3     =88028.172; % SW(Hz) for F3 dimension
        
    ph1     =[140.4,281];    % standard phase or
    ph2     =[-107,209];     % that found during the nmrPipe processing
    ph3     =[57.4,202];     % found in xyphase, all peaks inverted
    
    wnd1    =1;         % qsine window function for t1, on=1, off=0
    wnd2    =1;         % qsine window function for t2
    wnd3    =1;         % qsine window function for t3

    fctr1   =1.0;       % qsine window function ends at pi*fct1
    fctr2   =1.0;       
    fctr3   =1.0;       

    lab1    ='13Co';     % axis label for F1 dimension
    lab2    ='15N';     % axis label for F2 dimension
    lab3    ='13Cx';     % axis label for F3 dimension

    lowedge1 =(sw1/obs1)/2 + cent1; 
    lowedge2 =(sw2/obs2)/2 + cent2; 
    lowedge3 =(sw3/obs3)/2 + cent3;
disp('======================')
disp('section 2 is executed!')
disp('======================')





%% 
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
% 2. Import pre-SIFT time domain data
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

%% 2-1.

%--arguments--
    inname= '~/SIFT3DinMATLAB/example/pre.fid'; 
%--commands--
[read1,read2,read3,read4]=read3dnl(nt1,nt2,nt3,2048,inname,NUStab);
    pre_cc=read1;
    pre_sc=read2;
    pre_cs=read3;
    pre_ss=read4;

%% 2-2.
% Check if the data is correctly imported and reshuffled.

% Set the data just read in, and execute the section.
%--arguments--
    in1=pre_cc;
    in2=pre_sc;
    in3=pre_cs;
    in4=pre_ss;
%--commands--
    if isempty(in1)==0 &&isempty(in2)==0 &&isempty(in3)==0 &&isempty(in4)==0
    figure()
    subplot(2,1,1)
    plot(real(in1(:,1,1)))
    subplot(2,1,2)
	temp2=sineb3dv2(in1,90,1,1);
    NFFT=2^nextpow2(size(temp2,1));
    temp2=fft(temp2(:,1,1),NFFT);
    temp2=mphase3dv3(temp2,sw3,ph3,1);
    temp2=fftshift(real(temp2));
    plot(temp2)
    end
    if isempty(in1)==0 &&isempty(in2)==0 &&isempty(in3)==0 &&isempty(in4)==0
    clear temp
    [temp,t1F3t2,t2F1F3]=proc3dv2(in1,in2,in3,in4,nt1,nt2,nt3,sw1,sw2,sw3,ph1,ph2,ph3,wnd1,wnd2,wnd3,fctr1,fctr2,fctr3);
    [a,b]=max(real(temp));
    [c,d]=max(squeeze(a));
    [e,f]=max(squeeze(c));
    g=squeeze(b);
	pick=g(d(f),f);
    figure()
    subplot(2,1,1)
    plot(real(squeeze(t1F3t2(pick,:,1))),'.-')
    xlim([1 max(NUStab(:,1))])
    legend('t1')
    subplot(2,1,2)
    plot(real(squeeze(t2F1F3(pick,1,:))),'.-')
    xlim([1 max(NUStab(:,2))])
    legend('t2')
    figure()
    scatter(NUStab(:,1),NUStab(:,2),'.')
    grid on
    axis([1 max(NUStab(:,1)) 1 max(NUStab(:,2))])
    xlabel('t1')
    ylabel('t2')
    end

%% 2-3. 
% Axis swap or inversion

% Change flags (0=False, 1=True) according to the acquisition mode you
% used in the experiment, and execute this section.
%--arguments--
    in1=pre_cc;
    in2=pre_sc;
    in3=pre_cs;
    in4=pre_ss;
                        % in t1 -------
    alt1=0;             % States-TPPI/ Echo-Antiecho (axis left/right swap)
    neg1=0;             % Echo-Antiecho (axis inversion)
                        % in t2 -------
    alt2=0;             % States-TPPI/ Echo-Antiecho
    neg2=0;             % Echo-Antiecho
%--commands--
    [pre_cc,pre_sc,pre_cs,pre_ss]=acqmode(pre_cc,pre_sc,pre_cs,pre_ss,alt1,neg1,alt2,neg2);

    
    
    
    
    
    
    
%%    
%\\\\\\\\\\\\\\\\\\\\\\\\\\\
% 3. Scout for bright region
%\\\\\\\\\\\\\\\\\\\\\\\\\\\

% SIFT needs you to speficy bright spectral regions, or signal-containing
% regions, so that SIFT can use the dark spectral points to restore the NUS
% data.

%% 3-1.
% Process the pre-SIFT data

%--arguments--
    in1=pre_cc;
    in2=pre_sc;
    in3=pre_cs;
    in4=pre_ss;
    min_cont=0.1; % minimum contour level, ratio to the tallest peak
    max_cont=1.0; % maximum
    n_cont=10;
%--commands--
    if isempty(in1)==0 && isempty(in2)==0 && isempty(in3)==0 && isempty(in4)==0 
    [projf3f1_]=proj3dF3F1v3(in1,in2,[1 nzf3],nt1,nzf3,sw1,sw3,ph1,ph3,wnd1,wnd3,fctr1,fctr3);
    [projf3f2_]=proj3dF3F2v3(in1,in2,in3,in4,[1 nzf3],nt2,nzf3,sw2,sw3,ph2,ph3,wnd2,wnd3,fctr2,fctr3);
    clear in1 in2 in3 in4
    spec1=real(projf3f1_);
    spec2=real(projf3f2_);
    end
    if isempty(spec1)==0 && isempty(spec2)==0 
    cont=logspace(log10(min_cont),log10(max_cont),n_cont).*max(max(spec1));
    figure()
    subplot(1,2,1)
    mplp2d(squeeze(real(spec1)),sw3,sw1,cent3,cent1,obs3,obs1,cont);
    legend('projf3f1')
    xlabel('F1 (ppm)','fontsize',16)
    ylabel('F3 (ppm)','fontsize',16)
    subplot(1,2,2)
    mplp2d(squeeze(real(spec2)),sw3,sw2,cent3,cent2,obs3,obs2,cont);
    legend('projf3f2')
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F3 (ppm)','fontsize',16)
    end
    
%% 3-2
% Set the bright regions

% Look up the spectral limits from the spectra plotted in the section 3-1
%--arguments--
    brightLimits_F1=[168 180]; % [high-field limit, low-field limit]
    brightLimits_F2=[108 123]; 
    brightLimits_F3=[0   200]; % values found during preparation using nmrPipe 
%--commands--
if isempty(brightLimits_F1)==0 && isempty(brightLimits_F2)==0 && isempty(brightLimits_F3)==0 
    sigf1=[];
    sigf2=[];
    sigf3=[];
    nr=length(brightLimits_F1(:,1));
    for n=1:nr
        [limit1]=ppm2pts(brightLimits_F1(n,:),sw1,nzf1,cent1,obs1,[]);
        [limit2]=ppm2pts(brightLimits_F2(n,:),sw2,nzf2,cent2,obs2,[]);
        [limit3]=ppm2pts(brightLimits_F3(n,:),sw3,nzf3,cent3,obs3,[]);
        
        if cent1>brightLimits_F1(n,1) && cent1<brightLimits_F1(n,2)...
                && cent2>brightLimits_F2(n,1) && cent2<brightLimits_F2(n,2)...
                && cent3>brightLimits_F3(n,1) && cent3<brightLimits_F3(n,2)
            sigf1=[sigf1;
                1            limit1(1,2);
                1            limit1(1,2);
                1            limit1(1,2);
                1            limit1(1,2);
                limit1(1,1)  nzf1;
                limit1(1,1)  nzf1;
                limit1(1,1)  nzf1;
                limit1(1,1)  nzf1
                ];
            sigf2=[sigf2;
                1            limit2(1,2);
                1            limit2(1,2);
                limit2(1,1)  nzf2;
                limit2(1,1)  nzf2;
                1            limit2(1,2);
                1            limit2(1,2);
                limit2(1,1)  nzf2;
                limit2(1,1)  nzf2
                ];
            sigf3=[sigf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3
                ];
        elseif cent1>brightLimits_F1(n,1) && cent1<brightLimits_F1(n,2)...
                && cent2>brightLimits_F2(n,1) && cent2<brightLimits_F2(n,2)
            sigf1=[sigf1;
                1            limit1(1,2);
                1            limit1(1,2);
                limit1(1,1)  nzf1;
                limit1(1,1)  nzf1
                ];
            sigf2=[sigf2;
                1            limit2(1,2);
                limit2(1,1)  nzf2;
                1            limit2(1,2);
                limit2(1,1)  nzf2
                ];
            sigf3=[sigf3;
                limit3(1,1)  limit3(1,2);
                limit3(1,1)  limit3(1,2);
                limit3(1,1)  limit3(1,2);
                limit3(1,1)  limit3(1,2)
                ];
        elseif cent1>brightLimits_F1(n,1) && cent1<brightLimits_F1(n,2)...
                && cent3>brightLimits_F3(n,1) && cent3<brightLimits_F3(n,2)
            sigf1=[sigf1;
                1            limit1(1,2);
                1            limit1(1,2);
                limit1(1,1)  nzf1;
                limit1(1,1)  nzf1
                ];
            sigf2=[sigf2;
                limit2(1,1)  limit2(1,2);
                limit2(1,1)  limit2(1,2);
                limit2(1,1)  limit2(1,2);
                limit2(1,1)  limit2(1,2)
                ];
            sigf3=[sigf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3
                ];
        elseif cent2>brightLimits_F2(n,1) && cent2<brightLimits_F2(n,2)...
                && cent3>brightLimits_F3(n,1) && cent3<brightLimits_F3(n,2)
            sigf1=[sigf1;
                limit1(1,1)  limit1(1,2);
                limit1(1,1)  limit1(1,2);
                limit1(1,1)  limit1(1,2);
                limit1(1,1)  limit1(1,2)
                ];
            sigf2=[sigf2;
                1            limit2(1,2);
                1            limit2(1,2);
                limit2(1,1)  nzf2;
                limit2(1,1)  nzf2
                ];
            sigf3=[sigf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3
                ];
        elseif cent1>brightLimits_F1(n,1) && cent1<brightLimits_F1(n,2)
            sigf1=[sigf1;
                1            limit1(1,2);
                limit1(1,1)  nzf1
                ];
            sigf2=[sigf2;
                limit2(1,1)  limit2(1,2);
                limit2(1,1)  limit2(1,2)
                ];
            sigf3=[sigf3;
                limit3(1,1)  limit3(1,2);
                limit3(1,1)  limit3(1,2)
                ];
        elseif cent2>brightLimits_F2(n,1) && cent2<brightLimits_F2(n,2)
            sigf1=[sigf1;
                limit1(1,1)  limit1(1,2);
                limit1(1,1)  limit1(1,2)
                ];
            sigf2=[sigf2;
                1            limit2(1,2);
                limit2(1,1)  nzf2
                ];
            sigf3=[sigf3;
                limit3(1,1)  limit3(1,2);
                limit3(1,1)  limit3(1,2)
                ];
        elseif cent3>brightLimits_F3(n,1) && cent3<brightLimits_F3(n,2)
            sigf1=[sigf1;
                limit1(1,1)  limit1(1,2);
                limit1(1,1)  limit1(1,2)
                ];
            sigf2=[sigf2;
                limit2(1,1)  limit2(1,2);
                limit2(1,1)  limit2(1,2)
                ];
            sigf3=[sigf3;
                1            limit3(1,2);
                limit3(1,1)  nzf3
                ];
        else
            sigf1=[sigf1;
                limit1(1,1)  limit1(1,2)
                ];
            sigf2=[sigf2;
                limit2(1,1)  limit2(1,2)
                ];
            sigf3=[sigf3;
                limit3(1,1)  limit3(1,2)
                ];
        end
    end
end

disp('======================')
disp('section 3-2 is executed!')
disp('======================')







%%
%\\\\\\\\\
% 4. SIFT
%\\\\\\\\\

%% 4-1.
% Exploratory processing
%--arguments--
    ncyc=20;
    in1=pre_cc;
    in2=pre_sc;
    in3=pre_cs;
    in4=pre_ss;
%--commands--
    if isempty(in1)==0 && isempty(in2)==0 && isempty(in3)==0 && isempty(in4)==0
    ofst=0; % normally, set to zero.
    [junk1,junk2,junk3,junk4,sn]=sift3dv10(in1,in2,in3,in4,ncyc,NUStab,sigf1,sigf2,sigf3,nzf1,nzf2,nzf3,ofst,[],sw1,sw2);
    clear junk1 junk2 junk3 junk4
    end
    if isempty(sn)==0
    figure()
    subplot(2,1,1)
    plot(sn(:,1))
    legend('max')
    xlabel('Number of SIFT cycle','fontsize',16)
    subplot(2,1,2)
    plot(sn(:,3),'Color',[1 0 0])
    legend('s/n')
    xlabel('Number of SIFT cycle','fontsize',16)
    end
    
%% 4-2
% Rerun with the optimal number of cycles

%--arguments--
    ncyc=10;
    in1=pre_cc;
    in2=pre_sc;
    in3=pre_cs;
    in4=pre_ss;
%--commands--
    if isempty(in1)==0 && isempty(in2)==0 && isempty(in3)==0 && isempty(in4)==0
    [out1,out2,out3,out4]=sift3dv10(in1,in2,in3,in4,ncyc,NUStab,sigf1,sigf2,sigf3,nzf1,nzf2,nzf3,ofst,pick,sw1,sw2);
    end
    post_cc=out1;
    post_sc=out2;
    post_cs=out3;
    post_ss=out4;
    clear in1 in2 in3 in4 out1 out2 out3 out4

%% 
% Save all data in the Matlab working memory for later use.
    save ~/SIFT3DinMATLAB/example/workspace.mat
    
%% 
% To load saved set of data:
    load ~/SIFT3DinMATLAB/example/workspace.mat

    
    
    
    
    
    
    
%% 
%\\\\\\\\\\\\\
% 5. Analysis
%\\\\\\\\\\\\\
% You can output the time domain data to process them in nmrPipe, or
% process directly in Matlab for a quick overview.

%% 5-1 
% Output time domain data to bruker format, for processing data with
% nmrPipe:

% pre-SIFT data:

%--arguments--
    out1=pre_cc;
    out2=pre_sc;
    out3=pre_cs;
    out4=pre_ss;
	dest='~/SIFT3DinMATLAB/example/pre_reshuffled.bruker';
%--commands--
    if isempty(out1)==0 && isempty(out2)==0 && isempty(out3)==0 && isempty(out4)==0
    scale=1e5/max(max(max(real(out1))));
    in1=out1.*scale;
    in2=out2.*scale;
    in3=out3.*scale;
    in4=out4.*scale;
    brukout3dtimev1(in1,in2,in3,in4,nt1,nt2,nt3,dest)
    end

%%
% post-SIFT data

%--arguments--
    out1=post_cc;
    out2=post_sc;
    out3=post_cs;
    out4=post_ss;
	dest='~/SIFT3DinMATLAB/example/post.bruker';
%--commands--
    if isempty(out1)==0 && isempty(out2)==0 && isempty(out3)==0 && isempty(out4)==0
    scale=1e5/max(max(max(real(out1))));
    in1=out1.*scale;
    in2=out2.*scale;
    in3=out3.*scale;
    in4=out4.*scale;
    brukout3dtimev1(in1,in2,in3,in4,nt1,nt2,nt3,dest)
    end
    
    
    
%% 5-2
% Process data in Matlab.

% pre-SIFT data:
    in1=pre_cc;
    in2=pre_sc;
    in3=pre_cs;
    in4=pre_ss;
    [pre_ft3d]=proc3dv2(in1,in2,in3,in4,nzf1,nzf2,nzf3,sw1,sw2,sw3,ph1,ph2,ph3,wnd1,wnd2,wnd3,fctr1,fctr2,fctr3);
    clear in1 in2 in3 in4
%%
% post-SIFT data
    in1=post_cc;
    in2=post_sc;
    in3=post_cs;
    in4=post_ss;
    [post_ft3d]=proc3dv2(in1,in2,in3,in4,nzf1,nzf2,nzf3,sw1,sw2,sw3,ph1,ph2,ph3,wnd1,wnd2,wnd3,fctr1,fctr2,fctr3);
    clear in1 in2 in3 in4

%% 5-3
% Output spectra to Sparky format

    home='/Applications/Sparky.app/Contents/Resources/bin/matrix2ucsf'; % Sparky's home
    in1=pre_ft3d;
    ucsfout1= '~/SIFT3DinMATLAB/example/pre_ft3d';
    in2=post_ft3d;
    ucsfout2= '~/SIFT3DinMATLAB/example/post_ft3d';
    if isempty(in1)==0 && isempty(in2)==0
    binout3dfreq(in1,nzf2,nzf1,nzf3,obs2,obs1,obs3,sw2,sw1,sw3,lowedge2,lowedge1,lowedge3,lab2,lab1,lab3,ucsfout1,home);
    binout3dfreq(in2,nzf2,nzf1,nzf3,obs2,obs1,obs3,sw2,sw1,sw3,lowedge2,lowedge1,lowedge3,lab2,lab1,lab3,ucsfout2,home);
    end

%% 5-4
% Plot spectra in Matlab

% F2/F3 and F1/F3 projections
    in1=pre_cc;
    in2=pre_sc;
    in3=pre_cs;
    in4=pre_ss;
%     in1=post_cc;
%     in2=post_sc;
%     in3=post_cs;
%     in4=post_ss;
    label='pre-SIFT';
    min_cont=0.025;
    max_cont=1.0;
    n_cont=10;
    if isempty(in1)==0 && isempty(in2)==0 && isempty(in3)==0 && isempty(in4)==0 
    [projf3f1_]=proj3dF3F1v3(in1,in2,[1 nzf3],nzf1,nzf3,sw1,sw3,ph1,ph3,wnd1,wnd3,fctr1,fctr3);
    [projf3f2_]=proj3dF3F2v3(in1,in2,in3,in4,[1 nzf3],nzf2,nzf3,sw2,sw3,ph2,ph3,wnd2,wnd3,fctr2,fctr3);
    clear in1 in2 in3 in4
    spec1=real(projf3f1_);
    spec2=real(projf3f2_);
    end
    if isempty(spec1)==0 && isempty(spec2)==0 
    cont=logspace(log10(min_cont),log10(max_cont),n_cont).*max(max(spec1));
    figure()
    subplot(1,2,1)
    mplp2d(squeeze(real(spec1)),sw3,sw1,cent3,cent1,obs3,obs1,cont);
    legend('projf3f1')
    xlabel('F1 (ppm)','fontsize',16)
    ylabel('F3 (ppm)','fontsize',16)
    subplot(1,2,2)
    mplp2d(squeeze(real(spec2)),sw3,sw2,cent3,cent2,obs3,obs2,cont);
    legend('projf3f2')
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F3 (ppm)','fontsize',16)
    end

%% 
% F1/F2 plane cut at the tallest peak
    spec1=real(post_ft3d);
    spec2=real(pre_ft3d);
    min_cont=0.1;
    max_cont=1.0;
    n_cont=15;    
    label1='post-SIFT';
    label2='pre-SIFT';
if isempty(spec1)==0
    [a,b]=max(real(spec1),[],3);
    [c,d]=max(a,[],2);
    [e,f]=max(c);
    pick=f
    del3=(sw3/obs3)/(nzf3-1);
if pick>nzf3/2
    pickppm=cent3-del3*(nzf3-pick+1);
else
    pickppm=cent3+del3*(pick-1);
end

    figure()
    subplot(1,2,1)
    spec11=real(squeeze(spec1(pick,:,:)));
    con=logspace(log10(min_cont),log10(max_cont),n_cont).*max(max(max(real(spec1))));
    mplp2d(squeeze(real(spec11)),sw1,sw2,cent1,cent2,obs1,obs2,con);
    legend(label1)
    title(['F3=',num2str(pickppm)])
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F1 (ppm)','fontsize',16)

    subplot(1,2,2)
    spec22=real(squeeze(spec2(pick,:,:)));
    con=logspace(log10(min_cont),log10(max_cont),n_cont).*max(max(max(real(spec2))));
    mplp2d(squeeze(real(spec22)),sw1,sw2,cent1,cent2,obs1,obs2,con);
    legend(label2)
    title(['F3=',num2str(pickppm)])
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F1 (ppm)','fontsize',16)
end

%%
% 1D slices in the selected F1/F2 plane.
        F1=177.5;     % ppm
        F2=113;     % ppm
        if isempty(F1)==0 && isempty(F2)==0
            %-- ppm to pts
            del1=(sw1/obs1)/(nzf1-1);
            if F1>=cent1
                F1pts=floor((F1-cent1)/del1)+1;
            else
                F1pts=nzf1-floor((cent1-F1)/del1);
            end
            F1pts
            %----
            del2=(sw2/obs2)/(nzf2-1);
            if F2>=cent2
                F2pts=floor((F2-cent2)/del2)+1;
            else
                F2pts=nzf2-floor((cent2-F2)/del2);
            end
            F2pts

            figure()
            subplot(2,2,1)
            mplp2dslcv2(squeeze(spec11),sw1,sw2,cent1,cent2,obs1,obs2,1,F1pts);
            title(label1)
            legend(['F1=',num2str(F1)])
            xlabel('F2 (ppm)','fontsize',16)
            subplot(2,2,3)
            mplp2dslcv2(squeeze(spec11),sw1,sw2,cent1,cent2,obs1,obs2,2,F2pts);
            legend(['F2=',num2str(F2)])
            xlabel('F1 (ppm)','fontsize',16)

            subplot(2,2,2)
            mplp2dslcv2(squeeze(spec22),sw1,sw2,cent1,cent2,obs1,obs2,1,F1pts);
            title(label2)
            legend(['F1=',num2str(F1)])
            xlabel('F2 (ppm)','fontsize',16)
            subplot(2,2,4)
            mplp2dslcv2(squeeze(spec22),sw1,sw2,cent1,cent2,obs1,obs2,2,F2pts);
            legend(['F2=',num2str(F2)])
            xlabel('F1 (ppm)','fontsize',16)
        end

        % Change spectral limits in 1D slices
        subplot(2,2,1);
        xlim([lowedge2-sw2/obs2 lowedge2])
%         ylim([-2e4 14e4])

        subplot(2,2,3);
        xlim([lowedge1-sw1/obs1 lowedge1])
%         ylim([-2e4 14e4])

        subplot(2,2,2);
        xlim([lowedge2-sw2/obs2 lowedge2])
%         ylim([-2e4 14e4])

        subplot(2,2,4);
        xlim([lowedge1-sw1/obs1 lowedge1])
%         ylim([-2e4 14e4])
    
%%
% F1/F2 plane sliced at a given F3 (ppm): 
    spec1=post_ft3d;
    spec2=pre_ft3d;
    F3=178; %ppm
    con_min=0.1;
    con_max=1.00;
    ncon=15;
    label1='post-SIFT';
    label2='pre-SIFT';    
if isempty(spec1)==0

    del3=(sw3/obs3)/(nzf3-1);
if F3>=cent3
    F3pts=floor((F3-cent3)/del3)+1;
else
    F3pts=nzf3-floor((cent3-F3)/del3);
end
F3pts

    figure()
    subplot(1,2,1)
    spec1=real(squeeze(spec1(F3pts,:,:)));
    con=logspace(log10(con_min),log10(con_max),ncon).*max(max(max(real(spec1))));
    mplp2d(squeeze(real(spec1)),sw1,sw2,cent1,cent2,obs1,obs2,con);
    legend(label1)
    title(['F3=',num2str(F3)])
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F1 (ppm)','fontsize',16)

    subplot(1,2,2)
    spec2=real(squeeze(spec2(F3pts,:,:)));
    con=logspace(log10(con_min),log10(con_max),ncon).*max(max(max(real(spec2))));
    mplp2d(squeeze(real(spec2)),sw1,sw2,cent1,cent2,obs1,obs2,con);
    legend(label2)
    title(['F3=',num2str(F3)])
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F1 (ppm)','fontsize',16)
end
%%
% 1D slices in the selected F1/F2 plane.
        F1=179.5;     % ppm
        F2=113;     % ppm
        if isempty(F1)==0 && isempty(F2)==0
            %-- ppm to pts
            del1=(sw1/obs1)/(nzf1-1);
            if F1>=cent1
                F1pts=floor((F1-cent1)/del1)+1;
            else
                F1pts=nzf1-floor((cent1-F1)/del1);
            end
            F1pts
            %----
            del2=(sw2/obs2)/(nzf2-1);
            if F2>=cent2
                F2pts=floor((F2-cent2)/del2)+1;
            else
                F2pts=nzf2-floor((cent2-F2)/del2);
            end
            F2pts

            figure()
            subplot(2,2,1)
            mplp2dslcv2(squeeze(spec1),sw1,sw2,cent1,cent2,obs1,obs2,1,F1pts);
            title(label1)
            legend(['F1=',num2str(F1)])
            xlabel('F2 (ppm)','fontsize',16)
            subplot(2,2,3)
            mplp2dslcv2(squeeze(spec1),sw1,sw2,cent1,cent2,obs1,obs2,2,F2pts);
            legend(['F2=',num2str(F2)])
            xlabel('F1 (ppm)','fontsize',16)

            subplot(2,2,2)
            mplp2dslcv2(squeeze(spec2),sw1,sw2,cent1,cent2,obs1,obs2,1,F1pts);
            title(label2)
            legend(['F1=',num2str(F1)])
            xlabel('F2 (ppm)','fontsize',16)
            subplot(2,2,4)
            mplp2dslcv2(squeeze(spec2),sw1,sw2,cent1,cent2,obs1,obs2,2,F2pts);
            legend(['F2=',num2str(F2)])
            xlabel('F1 (ppm)','fontsize',16)
        end

        % Change spectral limits in 1D slices
        subplot(2,2,1);
        xlim([lowedge2-sw2/obs2 lowedge2])
%         ylim([-4e3 4e3])

        subplot(2,2,3);
        xlim([lowedge1-sw1/obs1 lowedge1])
%         ylim([-4e3 4e3])

        subplot(2,2,2);
        xlim([lowedge2-sw2/obs2 lowedge2])
%         ylim([-4e3 4e3])

        subplot(2,2,4);
        xlim([lowedge1-sw1/obs1 lowedge1])
%         ylim([-4e3 4e3])

%%
% F3/F1 plane sliced at a given F2 (ppm): 
    spec=post_ft3d;
    F2=113; %ppm
    con_min=0.03;
    con_max=1.00;
    ncon=15;
    label='post-SIFT';
if isempty(spec)==0

    del2=(sw2/obs2)/(nzf2-1);
if F2>=cent2
    F2pts=floor((F2-cent2)/del2)+1;
else
    F2pts=nzf2-floor((cent2-F2)/del2);
end
F2pts

    figure()
    spec=real(squeeze(spec(:,:,F2pts)));
    con=logspace(log10(con_min),log10(con_max),ncon).*max(max(max(real(spec))));
    mplp2d(squeeze(real(spec)),sw3,sw1,cent3,cent1,obs3,obs1,con);
    legend(label)
    title(['F2=',num2str(F2)])
    xlabel('F1 (ppm)','fontsize',16)
    ylabel('F3 (ppm)','fontsize',16)
end

%% Plot 1D slices in the F3/F1 plane.
        F3=171.8;     % ppm
        F1=172;     % ppm
        if isempty(F3)==0 && isempty(F1)==0
            %-- ppm to pts
            del1=(sw1/obs1)/(nzf1-1);
            if F1>=cent1
                F1pts=floor((F1-cent1)/del1)+1;
            else
                F1pts=nzf1-floor((cent1-F1)/del1);
            end
            F1pts
            %----
            del3=(sw3/obs3)/(nzf3-1);
            if F3>=cent3
                F3pts=floor((F3-cent3)/del3)+1;
            else
                F3pts=nzf3-floor((cent3-F3)/del3);
            end
            F3pts


            figure()
            subplot(2,1,1)
            mplp2dslcv2(squeeze(spec),sw3,sw1,cent3,cent1,obs3,obs1,1,F3pts);
            title(label)
            legend(['F3=',num2str(F3)])
            xlabel('F1 (ppm)','fontsize',16)
            subplot(2,1,2)
            mplp2dslcv2(squeeze(spec),sw3,sw1,cent3,cent1,obs3,obs1,2,F1pts);
            legend(['F1=',num2str(F1)])
            xlabel('F3 (ppm)','fontsize',16)
        end

        % Change spectral limits in 1D slices
        subplot(2,1,1);
        xlim([lowedge1-sw1/obs1 lowedge1])
        %     ylim([-1e8 4e8])

        subplot(2,1,2);
        xlim([lowedge3-sw3/obs3 lowedge3])
        %     ylim([-1e8 4e8])
    
%%
% F3/F2 plane sliced at a given F1 (ppm): 
    spec=post_ft3d;
    F1=177.5; %ppm
    con_min=0.025;
    con_max=1.00;
    ncon=15;
    label='post-SIFT';
if isempty(spec)==0

    del1=(sw1/obs1)/(nzf1-1);
if F1>=cent1
    F1pts=floor((F1-cent1)/del1)+1;
else
    F1pts=nzf1-floor((cent1-F1)/del1);
end
F1pts

    figure()
    spec=real(squeeze(spec(:,F1pts,:)));
    con=logspace(log10(con_min),log10(con_max),ncon).*max(max(max(real(spec))));
    mplp2d(squeeze(real(spec)),sw3,sw2,cent3,cent2,obs3,obs2,con);
    legend(label)
    title(['F1=',num2str(F1)])
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F3 (ppm)','fontsize',16)
end

%% Plot 1D slices in the F3/F2 plane.
        F3=177.5;     % ppm
        F2=113;     % ppm
        if isempty(F3)==0 && isempty(F2)==0
            %-- ppm to pts
            del3=(sw3/obs3)/(nzf3-1);
            if F3>=cent3
                F3pts=floor((F3-cent3)/del3)+1;
            else
                F3pts=nzf3-floor((cent3-F3)/del3);
            end
            F3pts
            %----
            del2=(sw2/obs2)/(nzf2-1);
            if F2>=cent2
                F2pts=floor((F2-cent2)/del2)+1;
            else
                F2pts=nzf2-floor((cent2-F2)/del2);
            end
            F2pts


            figure()
            subplot(2,1,1)
            mplp2dslcv2(squeeze(spec),sw3,sw2,cent3,cent2,obs3,obs2,1,F3pts);
            title(label)
            legend(['F3=',num2str(F3)])
            xlabel('F2 (ppm)','fontsize',16)
            subplot(2,1,2)
            mplp2dslcv2(squeeze(spec),sw3,sw2,cent3,cent2,obs3,obs2,2,F2pts);
            legend(['F2=',num2str(F2)])
            xlabel('F3 (ppm)','fontsize',16)
        end

        % Change spectral limits in 1D slices
        subplot(2,1,1);
        xlim([lowedge2-sw2/obs2 lowedge2])
        %     ylim([-1e8 4e8])
        subplot(2,1,2);
        xlim([lowedge3-sw3/obs3 lowedge3])
        %     ylim([-1e8 4e8])


% end of run-script.
    
    

