%% Raman Scanning Analysis Program by PAUL KLIMOV, Summer 2009.
close all; clear all; clc

%% THESE FEW LINES REQUIRE EDITING EVERY TIME A NEW FILE IS USED

load 2D_23X65_25mW_5s.txt; % title of the main file
load wavenumber.txt %title of the file that contains wavenumber info
Map=cell(46,110);
WN=wavenumber 
WNsmooth=linspace(min(WN),max(WN),1000)
DATA=X2D_23X65_25mW_5s';
clear data1 data2 wavenumber
xdim=46 % the number of spectra in x direction
ydim=110 % the number of spectra in y direction

%while i load the data, i will find noise and subtract it out
NoiseInd=find(WN>1800&WN<2300); %this selects a region of noise
for i = 1:ydim;
    for k = 1:xdim;
        Map{k,i}=DATA(:,k+(i-1)*xdim);
        NoiseVals=Map{k,i}(NoiseInd);
        MeanNoise=mean(NoiseVals);
        Map{k,i}=Map{k,i}-MeanNoise;
        Map2{k,i}=spline(WN,Map{k,i},WNsmooth); %interpolation
        
        %Normalize over 2D mode - this will be used later
        NormFactor=max(Map2{k,i}(find(WNsmooth>2500&WNsmooth<3000)));
        Normal2DMap{k,i}=Map2{k,i}/NormFactor;
        
        %Normalize over G mode - this will be used later
        NormFactorG=max(Map2{k,i}(find(WNsmooth>1000& WNsmooth<2100)));
        NormalGMap{k,i}=Map2{k,i}/NormFactorG;

    end
end
%% Fitting Lorentzian to the G mode

GInd=find(WNsmooth>1000& WNsmooth<2100);
PeakG=cell(xdim,ydim);
gamma=cell(xdim,ydim);
WidthG=cell(xdim,ydim);
CenterG=cell(xdim,ydim);
LorG=cell(xdim,ydim);
xG=WNsmooth(GInd)

for i = 1:ydim
    for k=1:xdim
    y=Map2{k,i}(GInd);

        [sigma,mu,AG]=mygaussfit(xG,y,.4); %optimized for .4
        Gauss=AG*exp( -(xG-mu).^2 / (2*sigma^2) ); % gaussian fit
        CenterG=xG(find(Gauss==max(Gauss))); % center of gaussian

        %Lorentzian fitting - finding gamma and the amplitude
        gammavec=linspace(0,50,200); 
        peakvec=linspace(AG+5, AG-5,10);
        wavenumbervec=linspace(CenterG-3, CenterG+3,10);
        
        % my fitting algorithm. It will scan through gammas, peak heights,
        % and frequencies to find the best fit. This part takes a long time
        % to run
        for u =1:10
            for l = 1:10
                for j = 1:200
                    testPeakG=peakvec(l); 
                    testgamma=gammavec(j);
                    testwn=wavenumbervec(u);
                    testLor=testPeakG*testgamma^2./((xG-testwn).^2+testgamma^2);
                    residual(j,u,l)=sum(abs(y-testLor)); % residuals in a 3 dim matrix
                end
            end
        end
        %next line will find the minima of the 3 dimensional residual
        %matrix and then from the indicies, my program will determine the
        %corresponding gamma, frequency, and peak height.
        [minj,minu,minl]=ind2sub(size(residual),find(residual==min(min(min(residual)))));
        
        gammaG{k,i}=gammavec(minj); %HWHM of Lorentzian
        PeakG{k,i}=peakvec(minl); %Peak intensity of Lorentzian
        wnG{k,i}=wavenumbervec(minu); %frequency of Lorentzian
        
        % The following lines look for errors and get rid of them without
        % crashing the program. They look for empty matricies, or matricies
        % with more than one element, which can result from spectra with
        % only noise.
        if sum(size(gammaG{k,i}))>2
            gammaG{k,i}=0;
            PeakG{k,i}=0;
            wnG{k,i}=0;
            errormatG{k,i}=1; %errormatG is a structure with 1's where there were problems
        end
        if sum(size(gammaG{k,i}))==0
            gammaG{k,i}=0;
            PeakG{k,i}=0;
            wnG{k,i}=0;
            errormatG{k,i}=1;
        end

        FWHMGaussG(k,i)=2.3548*sigma; %The gaussian FWHM. This is not used  anymore
        LorG{k,i}=PeakG{k,i}*gammaG{k,i}^2./((xG-wnG{k,i}).^2+gammaG{k,i}^2); %the lorentzian fit
        WidthG{k,i}=gammaG{k,i}; %HWHM of lorentzian
    i
    end
end
'done' %tells you when the program has finished running. 
%It can take over 30 numbers depending on the resolution


%% run this to clear some unimportant variables and free up memory
clear residual testLor testwn testgamma testPeakG gammavec peakvec wavenumbervec 
clear NoiseVals MeanNoise DATA Map
%% run this if you want to save out data to the current directory. You
%% Might want to add to these names to be more specific!
save WidthG2 WidthG
save LorG2 LorG
save PeakG2 PeakG
save gammaG2 gammaG
save wnG2 wnG

%% Run this if you want to load data from the current directory
%% If you changed names earlier, make sure to change these variables
load WidthG WidthG
load LorG LorG
load PeakG PeakG
load gammaG gammaG
load wnG wnG
%% Fitting a single Lorentzian to the 2D mode. Same algorithm as for G mode.
TwoDInd=find(WNsmooth>2550&WNsmooth<2850);
%fit the 2D mode with a lorentzian
PeakTwoD=cell(xdim,ydim);
gamma=cell(xdim,ydim);
WidthTwoD=cell(xdim,ydim);
FWHMGaussTwoD=cell(xdim,ydim);
CenterTwoD=cell(xdim,ydim);
Lor2D=cell(xdim,ydim)
x2D=WNsmooth(TwoDInd)

for i = 1:ydim
    for k=1:xdim
    y=Map2{k,i}(TwoDInd);

    [sigma,mu,A2D]=mygaussfit(x2D,y,.4); %optimized for .4
        Gauss=A2D*exp( -(x2D-mu).^2 / (2*sigma^2) );
        Center2D=x2D(find(Gauss==max(Gauss)));

        %Lorentzian fitting - finding gamma and the amplitude
        
        gammavec=linspace(0,50,200);
        peakvec=linspace(A2D+5, A2D-5,10);
        wavenumbervec=linspace(Center2D-3, Center2D+3,10);
        
        %residual=zeros(1,1000)
        for u =1:10
            for l = 1:10
                for j = 1:200
                    testPeak2D=peakvec(l);
                    testgamma=gammavec(j);
                    testwn=wavenumbervec(u);
                    testLor=testPeak2D*testgamma^2./((x2D-testwn).^2+testgamma^2);
                    residual(j,u,l)=sum(abs(y-testLor)); 
                end
            end
        end
        [minj,minu,minl]=ind2sub(size(residual),find(residual==min(min(min(residual)))));
        
        gamma2D{k,i}=gammavec(minj);
        Peak2D{k,i}=peakvec(minl);
        wn2D{k,i}=wavenumbervec(minu);
        
        if sum(size(gamma2D{k,i}))>2
            gamma2D{k,i}=0;
            Peak2D{k,i}=0;
            wn2D{k,i}=0;
            errormat2D{k,i}=1;
        end
        if sum(size(gamma2D{k,i}))==0
            gamma2D{k,i}=0;
            Peak2D{k,i}=0;
            wn2D{k,i}=0;
            errormat2D{k,i}=1;
        end
        FWHMGaussG(k,i)=2.3548*sigma;
        Lor2D{k,i}=Peak2D{k,i}*gamma2D{k,i}^2./((x2D-wn2D{k,i}).^2+gamma2D{k,i}^2);
        Width2D{k,i}=gamma2D{k,i};
    i
    end
end
'done'
%% Run this to save out 2D mode information
save Width2D_2 Width2D  %Format: save filename variable
save Lor2D_2 Lor2D
save Peak2D_2 Peak2D
save gamma2D_2 gamma2D
save wn2D_2 wn2D
%% Run this to load  2D mode information
load Width2D_2 Width2D
load Lor2D_2 Lor2D
load Peak2D_2 Peak2D
load gamma2D_2 gamma2D
load wn2D_2 wn2D
%% Run this to clear out some memory
clear A2D DATA Gauss MeanNoise NoiseInd NoiseVals NormFactor NormFactorG
clear testLor sigma residual peakvec testPeak2D testgamma testwn wavenumbervec
%% Ratio of integrated intensities of 2D and G modes
Gintegrationind=find(WNsmooth>1570&WNsmooth<1610); %finds region of integration
TwoDintegrationind=find(WNsmooth>2620&WNsmooth<2780); %finds region of integration

for i = 1:ydim
    for k=1:xdim
    IntegratedG(k,i)=sum(abs((Map2{k,i}(Gintegrationind)))); 
    Integrated2D(k,i)=sum(abs(Map2{k,i}(TwoDintegrationind)));
    end
end
Ratio2DtoG=Integrated2D./IntegratedG; % The integrated intensities ratio
imagesc(Ratio2DtoG,[0,5]);colorbar
%% Save Integrated Intensity data
save IntegratedG IntegratedG
save Integrated2D Integrated2D
save Ratio2DtoG Ratio2DtoG
%% Load Integrated Intensity data
load IntegratedG IntegratedG
load Integrated2D Integrated2D
load Ratio2DtoG Ratio2DtoG
%% Central Frequency calculations
WN2DcmInd=find(WNsmooth>2670&WNsmooth<2737)

for i = 1:ydim
    for k = 1:xdim
        Intensitysum(k,i)=sum(Normal2DMap{k,i}(WN2DcmInd));
        %first moment
        WNcm(k,i)=sum((WNsmooth(WN2DcmInd)).*Normal2DMap{k,i}(WN2DcmInd))/Intensitysum(k,i);
    end
end

imagesc(WNcm,[2700 2705])
colorbar; axis equal tight
title('first moment 2D')
set(gcf,'color',[1,1,1]); colormap hot
%% Save Moments
save WNcm WNcm
%% Plots
%%
%%
%%
%% Gmode Plots
subplot(1,2,1)
imagesc(cell2mat(WidthG)',[6.4 8])
colormap hot
colorbar
ylabel(texlabel('Gamma_G (cm^-^1)'))
set(gcf,'color',[1,1,1]); axis equal tight
set(gca,'YTickLabel',[]);
set(gca,'XTickLabel',[]);

%these next loops get rid of 'broken' data points, so we can plot the G
%frequencies without crashing the program
for i = 1:ydim
    for k = 1:xdim
        if length(find(LorG{k,i}==max(LorG{k,i})))>1
            shiftG(k,i)=0
            continue
        end
        shiftG(k,i)=wnG{find(LorG{k,i}==max(LorG{k,i}))};
    end
end

subplot(1,2,2)
imagesc(shiftG',[1590.2641 1590.2641])
colorbar;set(gcf,'color',[1,1,1])
colormap hot; axis equal tight
ylabel(texlabel('omega_G (cm^-^1)'))
set(gca,'YTickLabel',[]);
set(gca,'XTickLabel',[]);

figure
imagesc(cell2mat(PeakG),[1475 1760])
colorbar
title('G Mode Intensity')
set(gcf,'color',[1,1,1])
colormap hot

%2D to G integrated intensity
figure
subplot(3,1,1)
climratio=[1,1.76]*10^4
imagesc(IntegratedG,climratio)
colorbar;title('IG');colormap hot;set(gcf,'color',[1,1,1])

subplot(3,1,2)
climratio=[2,5]*10^4
imagesc(Integrated2D,climratio)
colorbar;title('I2D');colormap hot;set(gcf,'color',[1,1,1])

subplot(3,1,3)
climratio=[2,4]
imagesc(Integrated2D./IntegratedG,climratio)
colorbar;title('I2D/IG');colormap hot;set(gcf,'color',[1,1,1])

%% Widths at different heights for 2D mode and one G. This measures the
%% spectra directly, not the Lorentzians corresponding to the data
TwoDInd=find(WNsmooth>2550&WNsmooth<2850); %indecies of 2D mode
GInd2=find(WNsmooth>1570&WNsmooth<1620); %indicies of the G mode
for i = 1:ydim
    for k = 1:xdim
        %The following two lines will get indicies where the G and 2D modes
        %have greater than .75 of max intensity
        seventyfiveInd2D=find(Normal2DMap{k,i}(TwoDInd)>.75);
        seventyfiveIndG=find(NormalGMap{k,i}(GInd2)>..75);
        
        if isempty(seventyfiveInd2D)==1
            errorind(k,i)=1
            continue
        end

        if isempty(seventyfiveIndG)==1
            errorindG(k,i)=1
            continue
        end
        % the following two lines measure the width from the spectra
        % directly
        seventyfivewidth2D(k,i)=abs(WNsmooth(min(seventyfiveInd2D))-WNsmooth(max(seventyfiveInd2D)));
        seventyfivewidthG(k,i)=abs(WNsmooth(min(seventyfiveIndG))-WNsmooth(max(seventyfiveIndG)));
    end
end

%this sets bad data points to have very high widths so they take on
%'extreme'colors in the colormaps
seventyfivewidth(find(errorind==1))=100
seventyfivewidth2D(find(errorind==1))=100

figure
climsseventyfive=[34 49]
imagesc(seventyfivewidth,climsseventyfive)
title('seventy five width')
colorbar
axis equal tight; colormap hot
set(gcf,'color',[1,1,1])

climsseventyfiveG=[8 19]
imagesc(seventyfivewidthG,climsseventyfiveG)
title('fifty width G')
colorbar
axis equal tight
set(gcf,'color',[1,1,1])
colormap hot

%% The 2D mode averagedd in the regions. 
ABA2D=zeros(1,1000)
ABC2D=zeros(1,1000)

%the following code will add the spectra together in the denoted grid
for i = 1:12
    for j = 1:12
        tempspec=Normal2DMap{10+(i-1),15+j-1}; % must change indicies on different samples
        ABA2D= ABA2D + tempspec;
        tempspec=Normal2DMap{10+(i-1),53+j-1}; %must change these indicies as well
        ABC2D= ABC2D + tempspec;           
    end
end

%the following normalizes the resulting spectra
ABA2D=ABA2D/max(ABA2D);
ABC2D=ABC2D/max(ABC2D);

plot(WNsmooth,ABA2D,WNsmooth,ABC2D,'m','linewidth',2);axis([2600,2800,0,1]);
set(gca,'YTickLabel',[]);xlabel('cm^-^1'); legend('ABA','ABC')

%% histograms in the grid for the width
% to use the histogram functions, you have to have a one dim vector, so
% these next couple of lines reshape the vector for this purpose
AWidth=reshape(seventyfivewidth2D(10:22,15:27),[prod(size(seventyfivewidth2D(10:22,15:27))) 1])
BWidth=reshape(seventyfivewidth2D(10:22,53:65),[prod(size(seventyfivewidth2D(10:22,15:27))) 1])
XWidth=linspace(20,60,140)


[HAWidth,XAWidth]=hist(AWidth,11)%this gives you each x compnonent and the counts for each x
[HBWidth,XBWidth]=hist(BWidth,19)

XAWidth(1)-XAWidth(2) %this tells you the bin size
XBWidth(1)-XBWidth(2) %this tells you the bin size
% you would like both bin sizes to be very close to the same, so that they
% are essentially normalized. You have to play around with these lines

%You must adjust the amplitude of this curve so that it matches the histfit
%fit. The reason we have to do this is to extract the numerical data from
%the fit given by the histfit function. For some unknown reason, MATLAB
%does not allow one to do this...
YAWidth=88*sqrt(1/2/pi/std(AWidth)^2)*exp(-((XWidth-mean(AWidth)).^2)/(2*std(AWidth)^2))
histfit(AWidth,11); hold on
plot(XWidth,YAWidth,XAWidth,HAWidth,'o');
hold off

%repeat the same process here as the above.
YBWidth=88*sqrt(1/2/pi/std(BWidth)^2)*exp(-((XWidth-mean(BWidth)).^2)/(2*std(BWidth)^2))
histfit(BWidth,19); hold on
plot(XWidth,YBWidth,XBWidth,HBWidth,'o');
hold off

NormAWidth=YAWidth/sum(YAWidth)
NormBWidth=YBWidth/sum(YBWidth)

plot(XWidth,YAWidth,XWidth,YBWidth)
plot(XWidth,NormAWidth,XWidth,NormBWidth)

ABA.widthmean=mean(AWidth)
ABA.widthstd=std(AWidth)
ABA.widthintersect=41.3376
ABA.widthintegral=sum(NormAWidth(find(XWidth>ABA.widthintersect)))

ABC.widthmean=mean(BWidth)
ABC.widthstd=std(BWidth)
ABC.widthintegral=sum(NormBWidth(find(XWidth<ABA.widthintersect)))

%% histograms in the grid for first moment
% For an explanation of what to do here, please read the instructions given
% in the cell above this one.

AMom=reshape(WNcm(10:22,15:27),[prod(size(WNcm(10:22,15:27))) 1])
BMom=reshape(WNcm(10:22,53:65),[prod(size(WNcm(10:22,15:27))) 1])
XMom=linspace(2690,2710,400)

[HAMom,XAMom]=hist(AMom,13);
[HBMom,XBMom]=hist(BMom,10);

XAMom(1)-XAMom(2)
XBMom(1)-XBMom(2)

YAMom=37*sqrt(1/2/pi/std(AMom)^2)*exp(-((XMom-mean(AMom)).^2)/(2*std(AMom)^2))
histfit(AMom,13); hold on
plot(XMom,YAMom,XAMom,HAMom,'o');
hold off

YBMom=36*sqrt(1/2/pi/std(BMom)^2)*exp(-((XMom-mean(BMom)).^2)/(2*std(BMom)^2))
histfit(BMom,10); hold on
plot(XMom,YBMom,XBMom,HBMom,'o');
hold off

NormAMom=YAMom/sum(YAMom)
NormBMom=YBMom/sum(YBMom)

plot(XMom,YAMom,XMom,YBMom)
plot(XMom,NormAMom,XMom,NormBMom)

ABA.mommean=mean(AMom)
ABA.momstd=std(AMom)
ABA.momintersect=2701.0637 % this has to be entered manually
ABA.momintegral=sum(NormAMom(find(XMom>ABA.momintersect)))

ABC.mommean=mean(BMom)
ABC.momstd=std(BMom)
ABC.momintegral=sum(NormBMom(find(XMom<ABA.momintersect)))



%% Histograms for central frequency
AFreq=reshape(cell2mat(wn2D(10:22,15:27)),[prod(size(cell2mat(wn2D(10:22,15:27)))) 1])
BFreq=reshape(cell2mat(wn2D(10:22,53:65)),[prod(size(cell2mat(wn2D(10:22,15:27)))) 1])
XFreq=linspace(2690,2710,400)

[HAFreq,XAFreq]=hist(AFreq,15);
[HBFreq,XBFreq]=hist(BFreq,11);

XAFreq(1)-XAFreq(2)
XBFreq(1)-XBFreq(2)

YAFreq=89*sqrt(1/2/pi/std(AFreq)^2)*exp(-((XFreq-mean(AFreq)).^2)/(2*std(AFreq)^2))
histfit(AFreq,15); hold on
plot(XFreq,YAFreq,XAFreq,HAFreq,'o');
hold off

YBFreq=89*sqrt(1/2/pi/std(BFreq)^2)*exp(-((XFreq-mean(BFreq)).^2)/(2*std(BFreq)^2))
histfit(BFreq,11); hold on
plot(XFreq,YBFreq,XBFreq,HBFreq,'o');
hold off

NormAFreq=YAFreq/sum(YAFreq)
NormBFreq=YBFreq/sum(YBFreq)

plot(XFreq,YAFreq,XFreq,YBFreq)
plot(XFreq,NormAFreq,XFreq,NormBFreq)

ABA.freqmean=mean(AFreq)
ABA.freqstd=std(AFreq)
ABA.freqintersect=2696.1463 % this has to be entered manually
ABA.freqintegral=sum(NormAFreq(find(XFreq>ABA.freqintersect)))

ABC.freqmean=mean(BFreq)
ABC.freqstd=std(BFreq)
ABC.freqintegral=sum(NormBFreq(find(XFreq<ABA.freqintersect)))
%% Integrating two halves
% This data will integrate the left half of the 2D mode, then the right
% half, then take the ratio. 

for i = 1:ydim;
    for k=1:xdim
        
Data=Normal2DMap{k,i}(TwoDInd);
newindicies=find(Data>.15); %integrate from the point where the curves 
%reach 15% of maximum intensity

newdata=Data(newindicies);
centerindex=round(length(newindicies)/2); %this finds center index
LeftInt=sum(newdata(1:centerindex-3));
if sum(size(newdata))==1
    continue 
end
RightInt=sum(newdata(centerindex+3:end));
Ratio(k,i)=LeftInt/RightInt;
    end
end

imagesc(Ratio',[.5,1]); colorbar; axis equal tight; colormap hot
%% histograms in the grid for Ratio (asymmetry parameter)
ARat=reshape(Ratio(10:22,15:27),[prod(size(Ratio(10:22,15:27))) 1]);
BRat=reshape(Ratio(10:22,53:65),[prod(size(Ratio(10:22,15:27))) 1]);
XRat=linspace(0,2,400);

[HARat,XARat]=hist(ARat,10);
[HBRat,XBRat]=hist(BRat,11);

XARat(1)-XARat(2)
XBRat(1)-XBRat(2)

YARat=4.9*sqrt(1/2/pi/std(ARat)^2)*exp(-((XRat-mean(ARat)).^2)/(2*std(ARat)^2))
histfit(ARat,10); hold on
plot(XRat,YARat,XARat,HARat,'o');
hold off

YBRat=5.1*sqrt(1/2/pi/std(BRat)^2)*exp(-((XRat-mean(BRat)).^2)/(2*std(BRat)^2))
histfit(BRat,11); hold on
plot(XRat,YBRat,XBRat,HBRat,'o');
hold off

NormARat=YARat/sum(YARat)
NormBRat=YBRat/sum(YBRat)

plot(XRat,YARat,XRat,YBRat)
plot(XRat,NormARat,XRat,NormBRat)

ABA.ratmean=mean(ARat)
ABA.ratstd=std(ARat)
ABA.ratintersect=2.78 % this has to be entered manually
ABA.ratintegral=sum(NormARat(find(XRat>ABA.ratintersect)))

ABC.ratmean=mean(BRat)
ABC.ratstd=std(BRat)
ABC.ratintegral=sum(NormBRat(find(XRat<ABA.ratintersect)))

%% Residual parameter
smooth2Dmode=Normal2DMap{15,51}(TwoDInd) %This finds just the asymmetric ABC 2D Mode
for i = 1:ydim;
    for k=1:xdim
Data=Normal2DMap{k,i}(TwoDInd);
SUMDAT(k,i)=sum(abs(Data-smooth2Dmode)); %the main algorithm is given here
    end
end
clear Data %clears this variable
imagesc(SUMDAT,[0,8]);colorbar; colormap hot
%% histograms in the grid for SUM
%performs the same statistical procedure as given above
ASum=reshape(SUMDAT(10:22,15:27),[prod(size(SUMDAT(10:22,15:27))) 1]);
BSum=reshape(SUMDAT(10:22,53:65),[prod(size(SUMDAT(10:22,15:27))) 1]);
XSum=linspace(0,5,400);

[HASum,XASum]=hist(ASum,10);
[HBSum,XBSum]=hist(BSum,10);

XASum(1)-XASum(2)
XBSum(1)-XBSum(2)

YASum=54*sqrt(1/2/pi/std(ASum)^2)*exp(-((XSum-mean(ASum)).^2)/(2*std(ASum)^2))
histfit(ASum,10); hold on
plot(XSum,YASum,XASum,HASum,'o');
hold off

YBSum=52*sqrt(1/2/pi/std(BSum)^2)*exp(-((XSum-mean(BSum)).^2)/(2*std(BSum)^2))
histfit(BSum,10); hold on
plot(XSum,YBSum,XBSum,HBSum,'o');
hold off

NormASum=YASum/sum(YASum)
NormBSum=YBSum/sum(YBSum)

plot(XSum,YASum,XSum,YBSum)
plot(XSum,NormASum,XSum,NormBSum)

ABA.summean=mean(ASum)
ABA.sumstd=std(ASum)
ABA.sumintersect=2.78 % this has to be entered manually
ABA.sumintegral=sum(NormASum(find(XSum>ABA.sumintersect)))

ABC.summean=mean(BSum)
ABC.sumstd=std(BSum)
ABC.sumintegral=sum(NormBSum(find(XSum<ABA.sumintersect)))

%% Final Histogram  plots
subplot(1,4,4)
plot(XWidth,YAWidth,'b',XWidth,YBWidth,'m',XAWidth,HAWidth,'bo',...
    XBWidth,HBWidth,'mo','linewidth',3)
xlabel(texlabel('Gamma_2_D (cm^-^1)'))
subplot(1,4,1)
plot(XMom,YAMom,'b',XMom,YBMom,'m',XAMom,HAMom,'bo',...
    XBMom,HBMom,'mo','linewidth',3)
legend('ABA','ABC')
xlabel(texlabel('<omega_2_D> (cm^-^1)'));ylabel('counts')
subplot(1,4,2)
plot(XRat,YARat,'b',XRat,YBRat,'m',XARat,HARat,'bo',...
    XBRat,HBRat,'mo','linewidth',3)
xlabel(texlabel('I_R/I_L (a.u.)'))
set(gcf,'color',[1,1,1])
subplot(1,4,3)
plot(XSum,YASum,'b',XSum,YBSum,'m',XASum,HASum,'bo',...
    XBSum,HBSum,'mo','linewidth',3)
xlabel(texlabel('eta (a.u.)'))
set(gcf,'color',[1,1,1])

%% Multiple gaussian fits, using a downloaded algorithm. It works poorly,
%% however
for i = 1:ydim;
    for k = 1:xdim;
Data=Normal2DMap{k,i};
Data=Data(TwoDInd);
[u,sig,t,iter]=fit_mix_gaussian(Data',2);
U{k,i}=u;
SIG{k,i}=sig;
DIST(k,i)=abs(u(1)-u(2));
i
    end
end
imagesc(DIST',[.4,.52]);axis equal tight; colorbar

%% Plots for long strip 2

subplot(1,4,4)
imagesc(seventyfivewidth',[34 49])
colorbar; colormap hot; axis equal tight
set(gcf,'color',[1,1,1]);set(gca,'YTickLabel',[]);
set(gca,'XTickLabel',[]);
ylabel(texlabel('Gamma_2D (cm^-^1)'))
%set(gca,'YTickLabel',[]);xlabel('cm^-^1')

subplot(1,4,3)
imagesc(SUMDAT',[0,8]); axis equal tight
colorbar; colormap hot;
set(gcf,'color',[1,1,1]);set(gca,'YTickLabel',[]);
set(gca,'XTickLabel',[]);
ylabel(texlabel('eta (a.u.)'))

subplot(1,4,1)
imagesc(WNcm',[2700 2705])
colorbar; axis equal tight
ylabel(texlabel('<omega_2D> (cm^-^1)'))
set(gcf,'color',[1,1,1]); colormap hot;set(gca,'YTickLabel',[]);
set(gca,'XTickLabel',[]);

subplot(1,4,2)
imagesc(Ratio',[.5,1]); colorbar; axis equal tight; colormap hot
ylabel(texlabel('I_R/I_L (a.u.)'))
set(gcf,'color',[1,1,1]); colormap hot;set(gca,'YTickLabel',[]);
set(gca,'XTickLabel',[]);