% get the beam size from the HWS ETM source D1800125-v5_sn01

[out,r] = system('tar -xf HWS*.tar');

% load the files
dist = [1,10,29,51,84,105,140,180,240,295,351,435,490,565]; % beam propagation distance
files = dir('*.raw');
close all

% get indices
y = 1:1024;
y2 = y(:)*ones(1,1024);
x2 = y2';
p = ((1:1024) - 512)*12E-3;

wAll = [];
for ii = 1:numel(files)
    fid = fopen(files(ii).name, 'r');
    arr = fread(fid, [1024,1024], 'uint16')';
    fclose(fid);
    
    % get centroid
    arr2 = arr > 400;
    xc = round(sum(sum(arr2.^2.*x2))/sum(sum(arr2.^2)));
    yc = round(sum(sum(arr2.^2.*y2))/sum(sum(arr2.^2)));

    % get beam size
    arr2 = arr-50;
    wx = sqrt(sum(sum(arr2.^2.*(x2-xc).^2))/sum(sum(arr2.^2)));
    wy = sqrt(sum(sum(arr2.^2.*(y2-yc).^2))/sum(sum(arr2.^2)));
    
    dN = 420;
    imagesc(p, p, arr(yc-dN:yc+dN, xc-dN:xc+dN));
    axis square
    grid on
    W = 12E-3*mean([wx wy]);
    title(['Propagation distance, d = ', num2str(dist(ii)), 'mm, w0 = ', num2str(round(W*1e2)/1E2), 'mm'])
    wAll = [wAll; W];
    xlabel('Coordinates (mm)')
    ylabel('Coordinates (mm)')

    caxis([50 2500])
    pause(0.1);
    
end

figure('DefaultAxesFontSize', 14)
a = polyfit(dist(:), wAll(:), 1);
x = 0:1500;
plot(dist(:)', wAll(:), 'o')
hold all
plot(x, a(1)*x + a(2))
xlabel('Propagation distance (mm)')
ylabel('Second moment (mm)')
title(['Second moment vs distance from fiber [12,3",4,125] config - ', mfilename()], 'Interpreter', 'none')
grid on
orient landscape
print('-dpdf', [mfilename(), '.pdf'])