clc;clear;close all;
tic;
%% Plot coeffcients
FontName = 'Times New Roman';
FontSize = 9;

%limit setting
PlotMinW = 1;
PlotMinH = 1;
PlotMaxW = 3.4;
PlotMaxH = 1.5;

%% Fig.position
Fig1 = figure;
Fig1.Units = 'inches';
Fig1.InnerPosition = [PlotMinW, PlotMinH, PlotMaxW, PlotMaxH];

ax1 = axes;
ax1.InnerPosition = [0.14 0.25 0.985-0.14 0.97-0.25];%position with labels

%% load experimental results
%SMF
load('power_500_SMF.mat');

Power_dB(isnan(Power_dB))=[];
BinCount = 11;
h1=histogram(ax1,Power_dB,BinCount,'Normalization','pdf',...
            'DisplayName', 'Measured');%hold on;

%% fit curve parameters
StepMult = 10; %interpolation. should be even
IsStep = (max(h1.BinLimits)-min(h1.BinLimits))/BinCount/StepMult;
IsMin  = min(h1.BinLimits);%IsStep/2;
IsMax  = max(h1.BinLimits);
Is_dB = IsMin : IsStep : IsMax;
IsCircleIndex = StepMult/2+1 : StepMult : length(Is_dB);
IsAbs = 10.^(Is_dB/10);

%% generate GG curve.
sigmaR2 = 9.01;
I0Abs = 0.3165;

% a and b
beta02 = 0.4*sigmaR2; %spherical wave Rytov variance.
sigmaR = sqrt(sigmaR2);
beta0  = sqrt(beta02);

d = 1e-3; %diameter of receiver aperture

%sperical wave. 
a = (exp(0.49*beta0^2 / (1 + 0.18*d^2 + 0.56*beta0^(12/5))^(7/6)) - 1)^(-1);
b = (exp(0.51*beta0^2 * (1 + 0.69*beta0^(12/5))^(-5/6) ...
                 / (1+0.9*d^2+0.62*(d^2)*beta0^(12/5))) - 1)^(-1);

% Gamma-Gamma analytic.
fI = I0Abs^(-1) * 2/(gamma(a)*gamma(b)) * (a*b)^((a+b)/2) ...
.*(IsAbs/I0Abs).^((a+b)/2-1) .* besselk(a-b,2*sqrt(a*b*IsAbs/I0Abs));

fI = fI.*IsAbs.*log(10); 
  %f_dB*[10*log10(x+dx)-10*log10(x)]=f_abs*dx
  %f_dB = f_abs/10*log10(x)' = f_abs*I*ln(10)
  
fI = fI ./ sum(fI(IsCircleIndex)) ./ ((max(h1.BinLimits)-min(h1.BinLimits)) / BinCount);
  %normalization
sigmaI2 = 1/a + 1/b + 1/(a*b);

%% plot
plot(Is_dB(IsCircleIndex), h1.BinCounts/sum(h1.BinCounts)...
                          /((max(h1.BinLimits)-min(h1.BinLimits))/BinCount),...
     'rs', 'DisplayName', 'SMF experiment');hold on;
plot(ax1,Is_dB,fI,'r-','DisplayName', 'SMF curve fitting');hold on;
delete(h1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% load experimental results. FMF
load('power_500_FMF.mat');

Power_dB(isnan(Power_dB))=[];
BinCount = 11;
h2=histogram(ax1,Power_dB,BinCount,'Normalization','pdf',...
            'DisplayName', 'Measured');%hold on;

%% fit curve parameters
StepMult = 10; %interpolation. should be even
IsStep = (max(h2.BinLimits)-min(h2.BinLimits))/BinCount/StepMult;
IsMin  = min(h2.BinLimits);%IsStep/2;
IsMax  = max(h2.BinLimits);
Is_dB = IsMin : IsStep : IsMax;
IsCircleIndex = StepMult/2+1 : StepMult : length(Is_dB);
IsAbs = 10.^(Is_dB/10);

%% generate GG curve.
sigmaR2 = 1.45;
I0Abs = 0.9950;

% a and b
beta02 = 0.4*sigmaR2; %spherical wave Rytov variance.
sigmaR = sqrt(sigmaR2);
beta0  = sqrt(beta02);

d = 1e-3;%diameter of receiver aperture

%sperical wave.
a = (exp(0.49*beta0^2 / (1 + 0.18*d^2 + 0.56*beta0^(12/5))^(7/6)) - 1)^(-1);
b = (exp(0.51*beta0^2 * (1 + 0.69*beta0^(12/5))^(-5/6) ...
                 / (1+0.9*d^2+0.62*(d^2)*beta0^(12/5))) - 1)^(-1);

% Gamma-Gamma analytic.
fI = I0Abs^(-1) * 2/(gamma(a)*gamma(b)) * (a*b)^((a+b)/2) ...
.*(IsAbs/I0Abs).^((a+b)/2-1) .* besselk(a-b,2*sqrt(a*b*IsAbs/I0Abs));

fI = fI.*IsAbs.*log(10); 
  %f_dB*[10*log10(x+dx)-10*log10(x)]=f_abs*dx
  %f_dB = f_abs/10*log10(x)' = f_abs*I*ln(10)

fI = fI ./ sum(fI(IsCircleIndex)) ./ ((max(h2.BinLimits)-min(h2.BinLimits)) / BinCount);
  %normalization
sigmaI2 = 1/a + 1/b + 1/(a*b);

%% plot
plot(Is_dB(IsCircleIndex), h2.BinCounts/sum(h2.BinCounts)...
                          /((max(h2.BinLimits)-min(h2.BinLimits))/BinCount),...
     'bo', 'DisplayName', 'FMF experiment');hold on;
plot(ax1,Is_dB,fI,'b--', 'DisplayName', 'FMF curve fitting');hold on;
delete(h2);
xlabel('Received power [dB]');
ylabel('PDF');
legend('Location','northwest');
axis([-30 5 0 0.15]);
ax1.FontName = FontName;
ax1.FontSize = FontSize;

toc;