function D=beamdata(M,guess)
%Imports array of beam characterization measurements. Structure of M is 
% [z, x, I, a] where z is the displacement of the beam blockage along
% the optical axis, x is the coordinate of razor edge, I is the measured
% output of the photodetector and a is the ambient light level
%and guess is an estimate of the parameters [Amplitude x0 width] for the
%first measurement
%Output Structure [z A x0 w residual_2norm]
thisfile=mfilename('fullpath');
thisdir=strrep(thisfile,mfilename(),'');

zdat=M(:,1);
xdat=M(:,2)-M(:,4);
Idat=M(:,3);
measnum=size(xdat);
measnum=measnum(:,1);

D=[];
n=1;
N=[];
loopend=0;

while (n<(measnum)*loopend==0)==1
    i=1;
    if n+18<measnum
    while zdat(n+i-1)==zdat(n+i)
        i=i+1;
    end
    else
        loopend=1;
    end
    guess(1,2)=xdat(n+9,1);
    [x resnorm]=lsqcurvefit(@gsbeam,guess,xdat(n:n+i-1),Idat(n:n+i-1));
    D=[D;M(n,1) x resnorm];
    if n+18<measnum
        N=[N;n n+i-1];
    else
        N=[N;n measnum];
    end
    n=n+i;
    guess=x;
end
%plot fitted curves
[zpoints Dc]=size(D);
Dt=D;
Dt(:,3)=zeros(zpoints,1);
n=1;
xp=(round(min(xdat)-10):0.1:round(max(xdat)));
figure
grid on
hold
while n<(zpoints+1)
    yp=gsbeam(Dt(n,2:4),xp);
    plot(xp,yp,'Color',[Dt(1,1)/Dt(n,1) abs(0.5-Dt(1,1)/Dt(n,1)) 1-Dt(1,1)/Dt(n,1)]);
    n=n+1;
end
title('Gaussian Fit to Beam Data: All fit Curves (legend is z in mm)')
xlabel('x Distance of Razor from Beam Center (mm)')
ylabel('Photodetector Output (V)')
text(0.0,-0.1,[thisfile,'.m'],'FontSize',8,'Interpreter','none','Units','normalized')
legend(num2str(D(:,1)),'Location','Southeast')
%print to file
orient landscape
print('-dpdf',[thisfile '.pdf'])
hold

%plot with scatterplot
n=1;
while n<(zpoints+1)
    xp=round(min(M(N(n,1):N(n,2),2)))-2:0.1:round(max(M(N(n,1):N(n,2),2)))+2;
    figure
    hold
    grid on
    xlabel('Stage Translation (mm)')
    ylabel('Photodetector Output (V)')
    text(0.8,0.2,['A = ' num2str(D(n,2))],'FontSize',12,'Interpreter','none','Units','normalized');
    text(0.8,0.15,['x0 = ' num2str(D(n,3))],'FontSize',12,'Interpreter','none','Units','normalized');
    text(0.8,0.1,['w = ' num2str(D(n,3))],'FontSize',12,'Interpreter','none','Units','normalized');
    text(0.0,-0.13,[thisfile,'.m'],'FontSize',8,'Interpreter','none','Units','normalized')
    title(['Gaussian Fit to Beam Data: Fit Curves and Data: z=' num2str(D(n,1)) 'mm'])
    yp=gsbeam(D(n,2:4),xp);
    plot(xp,yp,'Color',[1 0 0]);
    scatter(M(N(n,1):N(n,2),2),M(N(n,1):N(n,2),3))
    orient landscape
    print('-dpdf',[thisfile 'z' num2str(n) '.pdf'])
    n=n+1;
    hold
    close all
end

%plot widths
ruler=D(1:25,1);
widths=D(1:25,4);
linreg=polyfit(ruler,widths,1);
intercept=linreg(1,2);
slope=linreg(1,1);

figure
hold
scatter(ruler,widths,'Marker','o')
xp=(min(ruler)-10):10:(max(ruler)+10);
yp=xp*slope+intercept;

plot(xp,yp)
title('Beam Width vs. Ruler Location');
xlabel('Location of Stage on Ruler (mm)');
ylabel('Calculated Beam Width (mm)');
text(0.8,0.1,['Slope=' num2str(slope)],'FontSize',12,'Interpreter','none','Units','normalized')
text(0.8,0.15,['Intercept=' num2str(intercept)],'FontSize',12,'Interpreter','none','Units','normalized')
text(0.8,0.05,['(last point not included in calculation)'],'FontSize',8,'Interpreter','none','Units','normalized')
text(0.0,-0.13,[thisfile,'.m'],'FontSize',8,'Interpreter','none','Units','normalized')
orient landscape
print('-dpdf',[thisfile 'widths.pdf'])

end