%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function reads the surface uniformity scan data                %
% and plots the results the results of the analysis.                  %
%                                                                     %
% Dependencies: optim.m                                               %
%                                                                     %
% Authors: Cornelia Bengea,                                           %
%          Mihai Cara, mcara@purdue.edu                               %
% Last Updated by Mihai Cara on 06/29/2004                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [v, nunif, figh, UserID, TestDate, TubeManufacturer,...
        TubeModel, SerialNumber, TubeGain, HighVoltage, LSPower, atten]...
    = pamunif(fd, file_ver, delta_r, delta_sig, az_radius, MaxIter, Weighted, attenr)
%Reading file for the uniformity test
fgetl(fd); % NL
fgetl(fd); % ** Info: **
UserID = fgetl(fd);%user identification
TestDate = fgetl(fd);%test date
TubeManufacturer = fgetl(fd);%the manufacturer of the tube
TubeModel = fgetl(fd);%the model of the tube
SerialNumber = fgetl(fd);%the serial number
TubeGain = fgetl(fd);%the gain of the tube
Comments = fgetl(fd);%Comments section related to the experiment
HighVoltage = fgetl(fd);%the high voltage setting
LSPower = fgetl(fd);%the Light Source Power setting

switch file_ver
case 1
    nl_skip = 62;
case 2
    nl_skip = 64;
case 3
    nl_skip = 66;
case 4
    nl_skip = 67;
otherwise
    disp('Unknown File Version!');
    return;
end

for i=1:nl_skip
    s = fgetl(fd); % reading of the rest of the header
end
temp = fscanf(fd, 'Data Section. Rows: %d, Columns: %d');
N = temp(1); % number of rows = size(Y)
M = temp(2); % number of columns = size(X)
s = fgetl(fd); % skip line
s = fgetl(fd); % skip line

s = fscanf(fd, 'X coordinates:');
X = fscanf(fd, '%g', M);%X coordinates
s = fgetl(fd); % skip line
s = fscanf(fd, 'Y coordinates:');
Y = fscanf(fd, '%g', N);%Y coordinates
s = fgetl(fd); % skip line
data = reshape(fscanf(fd, '%g', N*M), M, N)'; %the data aranged in a N*M matrix
e = reshape(fscanf(fd, '%g', N*M), M, N)';%the errors arranged in a N*M matrix
fclose(fd);

[v, bk] = optim([mean(mean(data)) 1.0 11.0 0.0 0.0],...
    MaxIter, N, M, X, Y, data, e, Weighted);
data = data - bk;
idata = ideal_surface(X, Y, v);
az_radius = v(3); %the effective radius used in the azimuthal profile
scrnsz = get(0,'ScreenSize');
figh = figure('Position', [5 42 scrnsz(3)-10 scrnsz(4)-100]);

%the surface graph of the data collected for the uniformity test
subplot(2,2,2);

surf(X,Y,data);
% dt=data;
% hold on;
% for i=1:N
%     for j= 1:M
%         dt(i,j)=ideal_surface(X(i),Y(j),v);
%     end
% end
% surf(X,Y,dt);
% hold off;
view(0,90);

xlabel('X-axis')
ylabel('Y-axis')
title('2-D View')
axis square;
axis tight;
colorbar;

%the subplot representing the radial profile for the collected data
rmax=sqrt(max(abs(X-v(4)))^2 + max(abs(Y-v(5)))^2); %the maximum radius
[x, y] = meshgrid((X-v(4)).^2, (Y-v(5)).^2);
rds = sqrt(x + y);
xbin = 0:delta_r:(rmax+2*delta_r);
binmax = size(xbin,2);
ybin = zeros(1,binmax);
count = zeros(1,binmax);
for i = 1:N
    for j = 1:M
        Bin = floor(rds(i,j)/delta_r)+1; %the number of bins
        ybin(Bin) = ybin(Bin)+data(i,j); %the current's intensity corresponding to each point in the matrix with N*M elements
        count(Bin) = count(Bin)+1; %the count
    end
end

x=[]; y=[]; k=1;
for i = 1:binmax
    if (count(i) >= 1)
        x(k) = xbin(i);
        y(k) = ybin(i)/count(i);
        k = k + 1;
    end
end

% Compute QE attenuation factors:
if (size(attenr,1)>size(attenr,2))
    attenr = attenr';
end
sizex = size(x,2);
sizeattenr = size(attenr, 2);
% get the average in the middle circle:
av0 = avrg(sizex, xbin, ybin, count, attenr(1));
atten = [];
for i = 2:sizeattenr
    av = avrg(sizex, xbin, ybin, count, attenr(i));
    atten(i-1,1) = attenr(i);
    atten(i-1,2) = av/av0;
end

% compute ki^2 within effective radius:
ki2 = 0.0;
count = 0;
for i = 1:N
    for j = 1:M
        if (rds(i,j) < v(3))
            ki2 = ki2 + ((data(i,j)-idata(i, j)))^2;
            count = count + 1;
        end
    end
end
nunif = 100*sqrt(ki2 / (count-5))/abs(v(1));

subplot(2,2,3);
plot(x,y);
hold on;
x = 0:0.25:max(x);
y1=[];
for i=1:size(x,2);
    y1(i) = ideal_surface1(x(i),v);
end
plot(x,y1,'r:');
clear y1;
hold off;
plot_axis = axis;
plot_axis(2) = ceil(max(x)+delta_r);
plot_axis(3) = 1.5*min([0 min(y)]);
plot_axis(4) = 1.05*max([v(1) max(y)]);
axis(plot_axis);

line([v(3) v(3)],[0 v(1)],'Color','r') %the vertical line on the graph 
line([0 v(3)],[v(1) v(1)],'Color','r') %the horizontal line on the graph 
xlabel('Radius (mm)')   %representing the  mean of our collected data
ylabel('Average Current')
title('Radial profile')
% text(0.6*v(3)/2, 0.9*v(1), ['Av = ' num2str(v(1))])%,'FontSize',18)
text(v(3), v(1), [' \leftarrow Av = ' num2str(v(1))])%,'FontSize',18)
text(v(3),v(1)/2,[' \leftarrow R_{eff} = ' num2str(v(3))])%,'FontSize',18)

%The subplot representing the azimuthal profile for the collected data
subplot(2,2,4);
bmax = floor(360/delta_sig); %maximum number of bins 
counta = zeros(bmax,1); %vector assigning the maximum length of the count
S = zeros(bmax,1);%vector assigning the maximum length of the current's intensity
sigmas = (delta_sig/2):delta_sig:(bmax*delta_sig); %vector of the angles in degrees

for n = 1:N
    for m = 1:M
        r = sqrt(X(m)^2+Y(n)^2); %the radia corresponding to each point in the matrix with N*M elements
        if (r > az_radius)     %disregards all the points greater than the azimuthal radius
            continue;
        end
        sig = atan2(Y(n), X(m));%calculates the four-quadrant inverse tangent
        if (sig < 0.0)
            sig = 2*pi + sig;
        end

        sig = 180*sig/pi;  %transformation of values of angles from radians to degrees 
        Bina = floor(sig/delta_sig)+1; %the number of bins
        S(Bina) = S(Bina)+data(n,m);   %the current's intensity corresponding to each point in the matrix with N*M elements
        counta(Bina) = counta(Bina)+1; %the count
    end
end
x=[]; y=[]; k=1;
for i = 1:bmax
    if (counta(i) >= 1)
        x(k) = sigmas(i);
        y(k) = S(i)/counta(i);
        k = k + 1;
    end
end
% plot(sigmas,S./counta)
plot(x,y);
plot_axis = axis;
plot_axis(2) = 360;
axis(plot_axis);
% axis tight;
% hold on;
% warning off MATLAB:polyfit:RepeatedPointsOrRescale;
% p=polyfit(x,y, floor(k/2)); y=polyval(p, x); plot(x,y, 'r');
% hold off;
xlabel('Angle(degrees)')
ylabel('Average Current')
title('Azimuthal profile')

%the subplot of the information about the test tube being tested, our settings 
%and the date at which the experiment was performed

subplot(2,2,1);
text(0, 0.5, sprintf(['******* PMT info: *******\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n'...
                      '******* Analysis: *******\n'...
                      'Non-uniformity, (%%): \\epsilon = %.3g\n'...
                      'Background, (A) = %.3e\n'...
                      'Eff. avrg., (A) = %g\nSlope, (mm^{-1}) = %g\n'...
                      'Eff. radius, (mm) = %g\n'...
                      'x_0, (mm) = %g\ny_0, (mm) = %g'],...
    TestDate,...
    TubeManufacturer, ...
    TubeModel,...
    SerialNumber,...
    HighVoltage,...
    LSPower,...
    TubeGain,...
    nunif,...
    bk,...
    v(1), v(2), v(3), v(4), v(5)));
axis off;

function is = ideal_surface1(r,v)
is = (0.5*v(1)).*(1-tanh(v(2).*(r-v(3))));

function is = ideal_surface2(r,v)
is = r*(0.5*v(1)).*(1-tanh(v(2).*(r-v(3))));

function is = ideal_surface(x,y,v)
[x1,y1]=meshgrid((x-v(4)).^2, (y-v(5)).^2);
is = (0.5*v(1)).*(1-tanh(v(2).*(sqrt(x1+y1)-v(3))));

function av = avrg(sz, xbin, ybin, count, r)
s = 0; k = 0;
for i = 1:sz
    if (xbin(i) <= r)
        s = s + ybin(i);
        k = k + count(i);
    end
end
av = s/k;
