%% Fitting routie for beam profiling measurments
%
% This model takes a set of data values for scaned beam widths as a
% function of distances an computes the waist and its position using the
% GaussProFit routine.  

clear all % clear the decks
close all

addpath('~/Box Sync/TCNlab_shared/standardMatlablibrary/') % Adds path to standard Matlab functions used for the TCN labwork

lambda = 1064e-9; % [m] Set wavelength of laser light


% The data measured with a WinCamD profiling camera in steps of distance
% away from the laser head

% Here the data is for the North laser (measured on June 17 2016) operating
% at 701.1 mW (Adj = 0).  The laser was allowed time to warm up and a pick
% off of the beam was taken by first reflecting off the front surface of a
% W1-PW1-1025-UV-1064-45UND UV grade fused silica window and then a
% W2-PW1-1037-UV-45P UV grade fused silica window with AR coating on back
% and front. The resulting light was ~200 uW.

z= [0 25 50 75 100 125 150 175 200 225 250 275 300 325 350 375 400 425]*1e-3 + (72e-3+25e-3+62e-3); % Distance from the laser head. The fixed added on value at the end is the distance to the first measurment point
W_horz = [768 829 866.2 915.7 965.3 1018.0 1072.3 1135.7 1180.0 1231.7 1285.1 1346.9 1402.5 1462.6 1517.0 1568.1 1648.3 1700.1]*1e-6/2; % Horizonal beam radius
W_vert = [760.4 762.8 870.1 971.8 983.2 1119.2 1223.5 1231.2 1353.5 1428.1 1522.5 1558.4 1619.4 1729.9 1813.9 1856.0 1888.6 1969.1]*1e-6/2; % Vertical beam radius

[z0h,w0h] = GaussProFit(z,W_horz);
[z0v,w0v] = GaussProFit(z,W_vert);




display(['Horz. beam waist = ' num2str(z0h*10^6) ' um'])
display(['Horz. beam waist position = ' num2str(w0h*10^3) ' mm'])

display(['Vert. beam waist = ' num2str(z0v*10^6) ' um'])
display(['Vert. beam waist position = ' num2str(w0v*10^3) ' mm'])

% And plot some values
w  = @(z,w0,z0) w0.*sqrt(1+(lambda.*(z-z0)./pi./w0.^2).^2); %beam radius as a function of distance from the waist z located at z0 and the beam radius at the waist z0
zPlot = linspace(-200e-3,max(z),1000); %Generate vector for plotting

figure(1)
plot(z*1e3, W_horz*1e6, 'ob',z*1e3, W_vert*1e6, 'or', zPlot*1e3, w(zPlot,z0h,w0h)*1e6,'-b',zPlot*1e3, w(zPlot,z0v,w0v)*1e6,'-r')
title('Beam profile data points for Horz. and vert. axis with fits as a function of z')
xlabel('z [mm]')
ylabel('beam radius [\mu m]')
legend('Vert. 1/e^2 radius data points','Horz. 1/e^2 radius data points','Vert. profile fit','Horz. profile fit')
set(gca,'fontsize', 18)
hold off