clear all;close all; clc;

% Load Lidar observation data
load LidarD.mat                 % Load Data
R_site      = 0.418;            % [km]
dR          = 0.03;             % [km] Altitude step
n_pow       = 141;
pow         = 0;
Rmin        = 2.3;
Rmax        = 5.9;
R           = R_site:dR:15;      % Alt [km]
eps         = 0.01;              % Error
Gr_alt      = 0;                 % Overlap Altitude [km]
figN        = 1;                 % Figure Count
n           = 6;                 % Polynomial degree
indm = round((Rmin - R_site)/dR + 1);
indx = round((Rmax - R_site)/dR + 1);
Pow_vec = {raw_mat1 raw_mat2 raw_mat3 raw_mat4};

clear raw_mat1  raw_mat2 raw_mat3 raw_mat4
labels = {'Parallel - Analog';'Parallel - Digital';...
    'Perpendicular - Analog';'Perpendicular - Digital'};

for f_ind = 1:1

    raw_mat = Pow_vec{1,f_ind};
    %Calculating average power profile
    pow = 0;
    for i=1:n_pow
        pow = pow + raw_mat(1:indx,i);
    end
    pow = pow'/n_pow;                 % Average power profile
    plog = log(pow.*R(1:indx).^2);    % Average log power profile
    %plog = log(R(1:indx).^2);
   % plog = log(R(1:indx).^2)
    figure(figN);
    plot(pow,R(1:indx),'.');
    hold on; grid on;
    ylabel('Altitude [km] ');
    xlabel('Power');
    title(['Average Power Profile(N = 141)' labels(f_ind,1)]);
    hold off;
    figN = figN + 1;
%     % Average log power profile
%     figure(figN);
%     plot(plog,R(1:indx),'.g');
%     hold on; grid on;
%     ylabel('Altitude [km] ');
%     xlabel('log(P*R^2)');
%     title(['Average Log Power Profile(N = 141) log(P*R^2)' labels(f_ind,1)]);
%     hold off;
%     figN = figN + 1;

Polynomial fit to the power profile

    A = polyfit(R(indm:indx),pow(indm:indx),n);
    Pint = 0;
    for i =1:n+1
         Pint  = Pint + A(i)*R(1:indx).^(n+1-i);
    end
    figure(figN);hold on; grid on;
    plot(pow(indm:indx),R(indm:indx),'*b');
    plot(Pint,R(1:indx),'.r');
    ylabel('Altitude [km] ');
    xlabel('Power');
    title(['Polynomial fit to average power profile' labels(f_ind,1)]);
    hold off;
    legend('Actual', 'Fit')
    figN = figN + 1;

Geometric factor derived from average power profile

    S   =  log(R(1:indx).^2.*pow);
    Sint = log(R(1:indx).^2.*Pint);
    GR_int = exp(S - Sint);         %Geometric factor,Overlap factor

    figure(figN);hold on; grid on;
    plot(R(1:indx),GR_int,'.r');
    xlabel('Altitude [km] ');
    ylabel('Geometric factor G(R)');
    title(['Geometric factor derived from average power profile'...
       labels(f_ind,1)]);
    grid on;
    hold off;
    figN = figN + 1;

Calculating geometric form factor for each power profile

     for j = 1:n_pow
        pow = raw_mat(1:indx,j);
        pow = pow';
        A = polyfit(R(indm:indx),pow(indm:indx),n);
        Pint = 0;
        for i = 1:n+1
            Pint  = Pint + A(i)*R(1:indx).^(n+1-i);
        end
        S   =  log(R(1:indx).^2.*pow);
        Sint = log(R(1:indx).^2.*Pint);
        GR_int = exp(S - Sint);     %Geometric factor
        for k = 1:indx
            if (abs(GR_int(k) - 1.0) < eps )
                Gr_alt(j) = R(k);
                break;
            end
        end
     end
    figure(figN);
    hist(Gr_alt,20);
    xlabel('Altitude [km]');
    ylabel('Number of Power profiles');
    title(['G(R) altitude variation for various power profiles',...
        labels(f_ind,1)]);
    figN = figN + 1;
%     %% Mean altitude overlapping occurs, Sea Level
%     mean_Gr_alt = mean(Gr_alt);
%     std_Gr_alt =  std(Gr_alt);
%     display(labels(f_ind,1));
%     fprintf('Overlap Alt = %4.2f +- %4.2f [km] \n\n',mean_Gr_alt,std_Gr_alt);

end