function simulAttenuation_dlr
% function simulAttenuation_dlr
%
% The program simulates total signal attenuation for DLR-Oberpfaffenhofen for various
% time periods (spring, summer, fall, winter).
% The choice is made by typing the parameters as requested by the program.
% The programm call functions which contain the polynomial fittings for the specific channel model
% parameters (mean values and standard deviations). The name of these functions is created by this
% program in line 27 and 28

% input parameters
% -------------------------------------------------------------------------
groundStation = input('Ground station (dlr) : ', 's');

period = input(['Season (spring, summer, fall or winter) '], 's');

nbDaysMax = 0;
while nbDaysMax<10
    nbDaysMax = input('Number of simulation days (>10): ');
end


% creates string for the function pMn and pSig
% -------------------------------------------------------------------------
pMn = ['mu_' groundStation '_' period];
pSig = ['sigma_' groundStation '_' period];


% call the simulator
% -------------------------------------------------------------------------
try
    simulator(nbDaysMax, pMn, pSig);
catch
    disp(['Error in the parameters definition: ' pMn ' and ' pSig ' are not found, error: ' lasterr])
end



% **************************************************************************
function simulator(nbDaysMax, pMn, pSig)
% function simulator(path, nbDaysMax, pMn, pSig)
% nbDaysMax : number of days for the simulation
% pMn       : function for mean
% pSig      : function for sigma
%
% Returns three graphics
%       1. the ten first days of simulation
%       2. the probability of exceeding attenuation
%       3. the fade duration probability for xdB = 5, 10, 15, 20, 25 and 30 dB
%


% initialization
% -------------------------------------------------------------------------
delta = 1;      % to define the segment type
dBSize = 0.1;   % resolution of the attenuation scale (for the cumulative attenutation)    

% determines 2 values to start and initialize the data vector
d2_val = rand(1) + 1;
d1_val = rand(1) + 1;
signalData = [d2_val d1_val];

% initializes variables
Number_data = 2;
nbDays = 0;

% set up constants
maxSimul = 50;
minSimul = 0.01; % trick : not 0 for the indexing of "attDistrib"

% number of samples for one day
modelStep = 64; % in secondes
nbSamplesDay = (24 * 3600)/modelStep;
nbPlotDays = 10;
nbTotSamples = nbDaysMax * nbSamplesDay;

% set up scales
timeScale = modelStep * [0 : 1 : (nbPlotDays*nbSamplesDay-1)];

% initialize attDistrib
attDistrib = zeros(1,ceil(maxSimul/dBSize));

% initilize fadingDistrib
maxFading = 100000;
stepFading_Distrib = modelStep;
fadingLimits = [5:5:30]';     % fading limits
fadingDur  = zeros(length(fadingLimits),1);       % duration of fading
fadingDistrib = zeros(size(fadingLimits,1), ceil(maxFading/stepFading_Distrib));    %count the time for each limit


% simulator
% -------------------------------------------------------------------------
while nbDays < nbDaysMax
    
    Number_data = Number_data + 1;
    if ( mod(Number_data, nbSamplesDay) == 0)
        nbDays = nbDays + 1
    end

    % Calculate the new value
    nextValue = callNextValue(d1_val, d2_val, delta, pMn, pSig);
    
    % Limit the value if necessairy
    if nextValue > maxSimul
        nextValue = maxSimul;
    elseif nextValue < minSimul
        nextValue = minSimul;
    end
    
    % add the value if during the "nbPlotDays" first days
    if Number_data <= (nbPlotDays*nbSamplesDay)
        signalData = [signalData nextValue];
    end
    
    % For the attenuation distribution
    try
    attDistrib(ceil(nextValue/dBSize)) = attDistrib(ceil(nextValue/dBSize)) + 1;
catch
    d=3;
end
    % For the fading duration
    for j=1:(length(fadingLimits))
        if ((d1_val < fadingLimits(j)) & (nextValue >= fadingLimits(j)))
            fadingDur(j) = 1;
        elseif ((d1_val >= fadingLimits(j)) & (nextValue < fadingLimits(j)))
            fadingDistrib(j, fadingDur(j)) = fadingDistrib(j, fadingDur(j)) +1;
            fadingDur(j) = 0;
        elseif fadingDur(j) ~= 0
            fadingDur(j) = fadingDur(j) + 1;
        end
    end
    
    % memory for the next run
    d2_val = d1_val;
    d1_val = nextValue;
end


% ******************************************************************************
% plot the results
% ******************************************************************************

% distrib   : distribution of attenuation
% Step_att  : step of attenuation
% Max_att   : maximum of attenuation
% refCul_att: reference distribution of attenuation
% refCul_prob : 
% Axes_handle

% for results of the first ten days 
figure;
plot(timeScale/(3600*24), -signalData)

% Set up the display
xlabel('Time (days)');
ylabel('Simulated received power (dB)');
title('Simulated received signal','Fontsize',12);
set(gcf, 'numbertitle', 'off', 'name', 'Simulated received signal');
axis([0 10 -40 0]);
grid('on');


% for results of the cumulative distribution
% ---------------------------------------------------------------------------
cumDistrib = fliplr(cumsum(fliplr(attDistrib)));
cumDistribProb = cumDistrib/sum(attDistrib);

Max_N = ceil(maxSimul/dBSize);
Att_scale = dBSize * [0:1:Max_N-1];

% Plot the current statistics
figure
semilogy(Att_scale, cumDistribProb, 'k.-');

% Set up the display
xlabel('Attenuation (dB)');
ylabel('Probability of exceeding attenuation');
title('P(x)','Fontsize',12);
set(gcf, 'numbertitle', 'off', 'name', 'Cumulative distribution');
axis([0 maxSimul 0.000001 1]);
grid('on');


% plot the duration cumulative function
% ----------------------------------------------------------------------------
durTimeScale = modelStep * [0:size(fadingDistrib, 2)-1];
durTimeScale(1) = 2;

temp=[];
PDeltaT = [];
PDeltaTEvent = [];
for j=1:size(fadingDistrib, 1)
    %calculate the probability to be above x dB of attenuation
    probaXdB = (fadingDistrib(j,:) *  (modelStep *[1:length(fadingDistrib(j,:))]')) ...
                ./ (modelStep * (nbTotSamples-1));
    % calculate the total time of fading fo each fading increment
    temp = [(fadingDistrib(j,:) .* [1:length(fadingDistrib)]) ];
    if sum(temp)~=0 % to avoid error if dividing by 0
        PDeltaT = [PDeltaT ; probaXdB * fliplr(cumsum(fliplr(temp))) ...
                                    ./ ( sum(temp) .* ones(size(temp)) )];
    else
        PDeltaT = [PDeltaT ; zeros(size(temp))];
    end
end

%display the current statistics
figure
for j=1:size(PDeltaT, 1)
    loglog(durTimeScale, PDeltaT(j,:), '+-', 'color', 'k');
    hold on
end

% Set up the display
xlabel('\Deltat (seconds)');
ylabel('P(x, \Deltat)');
title('P(x, \Deltat)', 'Fontsize', 12);
axis([0 maxFading 0.00001 1]);
set(gcf, 'numbertitle', 'off', 'name', 'P(x, Delta t)');
grid('on');
        


% *************************************************************************
function nextValue = callNextValue(val_d1, val_d2, DeltadB, pMn, pSig)
% function Next_value5 = cal_next_value5(val_d1, val_d2, DeltadB, pMn, pSig)
% val_d1    : value at the distance d1 from the point to calculate
% val_d2    : value at the distance d2 from the point to calculate
% DeltadB   : delta for segment type determination
% pMn       : polynom for mean determination
% pSig      : polynom for sigma determination
%
% Determines the next value from functions or polynomials pMn and pSig without
% any quantization
%

% Determination of the segment type
% -------------------------------------------------------------------------
diff = val_d1-val_d2;
if abs(diff) <= DeltadB
    segStyle = 2; % CONSTANT
elseif diff > DeltadB
    segStyle = 1; % DOWN
else
    segStyle = 3; % UP
end

% case of only one segment type (first order model)
if ~isa(pMn, 'char') & (size(pMn,2)==1) & (size(pSig,2)==1)
    segStyle=1;
end

% Determination of mean and sigma values for the random lottery
% -------------------------------------------------------------------------
m = feval(pMn, val_d1, segStyle);
s = feval(pSig, val_d1, segStyle);

% Calculate the next value
% -------------------------------------------------------------------------
nextValue = normrnd(m, s);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SUB FUNCTIONS TO DEFINE THE PARAMTERS SIGMA AND MU
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%***************************************************************************
% DLR 
%***************************************************************************

% spring
%---------------------------------------------------------------------------

function mu = mu_dlr_spring(val, segType)
% function mu = mu_dlr_spring(val, segType)
% mu        : mean value
% val       : value (vector or single value) to evaluate the mean
% segType   : segment type (1, 2 or 3 for D, C, or U)
%
% Definition of the analytical function for mean at the dlr, 
% season "spring".
%
% This is defined as a real Matlab function like cos or exp since the 
% input "val2" can be a vector or a single value.
% So the function can be plotted. For example for a Down segment:
%       plot([0:0.01:50], mu_dlr_spring([0:0.01:50], 2))
%
% The polynom is defined by:
%   mu = x0 + x1.val + x2.val^2
% for each segment type D, C and U.
% 
% pMn mapping:
% pMn = [   D2    C2    U2
%           D1    C1    U1
%           D0    C0    U0  ]
%

mu = zeros(size(val));

pMu = [ -0.0024 -0.0024 -0.0024;
         1.0700  0.9970  0.9200;
         0.1450  0.0160  0.0500];
     
mu = polyval(pMu(:,segType), val);


function sig = sigma_dlr_spring(val, segType)
% function sig = sigma_dlr_spring(val, segType)
% sig       : standard deviation value
% val       : value (vector or single value) to evaluate the standard deviation
% segType   : segment type (1, 2 or 3 for D, C, or U)
%
% Definition of the analytical function for standard deviation at the
% DLR, season "spring".
%
% This is defined as a real matlab function like cos or exp since the 
% input "val2" can be a vector or a single value.
% So the function can be plotted. For example for a Down segment:
%       plot([0:0.01:50], sigma_dlr_spring([0:0.01:50], 2))
%

sig = zeros(size(val));

if segType == 1
    % val < 25
    ind1 = (val<=25);
    sig(ind1) = polyval([0.000 0.240 0.100], val(ind1));
    % val >25
    ind2 = (val>25);
    sig(ind2) = 6.100;
elseif segType==2
    % val < 2
    ind1 = (val<=2);
    sig(ind1) = 0.120;
    % 2 < val <= 25
    ind2 = (val>2 & val<=25);
    sig(ind2) = polyval([0.000 0.140 -0.160], val(ind2));
    % val > 25
    ind3 = (val>25);
    sig(ind3) = 3.340;
elseif segType==3
    % val <= 25
    ind1 = (val<=25);
    sig(ind1) = polyval([0.000 0.160 0.020], val(ind1));
    % val > 25
    ind2 = (val>25);
    sig(ind2) = 4.020;
end


% summer
%---------------------------------------------------------------------------

function mu = mu_dlr_summer(val, segType)
% function mu = mu_dlr_summer(val, segType)
% mu        : mean value
% val       : value (vector or single value) to evaluate the mean
% segType   : segment type (1, 2 or 3 for D, C, or U)
%
% Definition of the analytical function for mean at the DLR, 
% season "summer".
%
% This is defined as a real Matlab function like cos or exp since the 
% input "val" can be a vector or a single value.
% So the function can be plotted. For example for a Down segment:
%       plot([0:0.01:50], mu_dlr_summer([0:0.01:50], 2))
%
% The polynom is defined by:
%   mu = x0 + x1.val + x2.val^2
% for each segment type D, C and U.
% 
% pMn mapping:
% pMn = [   D2    C2    U2
%           D1    C1    U1
%           D0    C0    U0  ]
%

pMu = [ -0.0015 -0.0007 -0.0005;
         1.0800  0.9930  0.8990;
         0.1900  0.0220  0.1000];
     
mu = polyval(pMu(:,segType), val);


function sig = sigma_dlr_summer(val, segType)
% function sig = sigma_summer_dlr(val, segType)
% sig       : standard deviation value
% val       : value (vector or single value) to evaluate the standard deviation
% segType   : segment type (1, 2 or 3 for D, C, or U)
%
% Definition of the analytical function for standard deviation at the
% DLR, season "summer".
%
% This is defined as a real Matlab function like cos or exp since the 
% input "val" can be a vector or a single value.
% So the function can be plotted. For example for a Down segment:
%       plot([0:0.01:50], sigma_dlr_summer([0:0.01:50], 2))
%

sig = zeros(size(val));

if segType == 1
    % val <= 25
    ind1 = (val<=25);
    sig(ind1) = polyval([0.000 0.190 0.500], val(ind1));
    % val > 25
    ind2 = (val>25);
    sig(ind2) = 5.250;
elseif segType==2
    % val <= 2
    ind1 = (val<=2);
    sig(ind1) = 0.120;
    % 2 < val <= 25
    ind2 = (val>2 & val<=25);
    sig(ind2) = polyval([0.000 0.140 -0.160], val(ind2));
    % val > 25
    ind3 = (val>25);
    sig(ind3) = 3.34;
elseif segType==3
    % val <= 25
    ind1 = (val<=25);
    sig(ind1) = polyval([0.000 0.155 0.170 ], val(ind1));
    % val > 25
    ind2 = (val>25);
    sig(ind2) = 4.045;
end


% fall
%---------------------------------------------------------------------------

function mu = mu_dlr_fall(val, segType)
% function mu = mu_dlr_fall(val, segType)
% mu        : mean value
% val       : value (vector or single value) to evaluate the mean
% segType   : segment type (1, 2 or 3 for D, C, or U)
%
% Definition of the analytical function for mean at the dlr, 
% season "fall".
%
% This is defined as a real Matlab function like cos or exp since the 
% input "val2" can be a vector or a single value.
% So the function can be plotted. For example for a Down segment:
%       plot([0:0.01:50], mu_dlr_fall([0:0.01:50], 2))
%
% The polynom is defined by:
%   mu = x0 + x1.val + x2.val^2
% for each segment type D, C and U.
% 
% pMn mapping:
% pMn = [   D2    C2    U2
%           D1    C1    U1
%           D0    C0    U0  ]
%

mu = zeros(size(val));

pMu = [ -0.0010 -0.0025 0.0000;
         1.0400  0.9990 0.9000;
         0.1000  0.0155 0.0500];
     
mu = polyval(pMu(:,segType), val);


function sig = sigma_dlr_fall(val, segType)
% function sig = sigma_dlr_fall(val, segType)
% sig       : standard deviation value
% val       : value (vector or single value) to evaluate the standard deviation
% segType   : segment type (1, 2 or 3 for D, C, or U)
%
% Definition of the analytical function for standard deviation at the
% DLR, season "fall".
%
% This is defined as a real Matlab function like cos or exp since the 
% input "val" can be a vector or a single value.
% So the function can be plotted. For example for a Down segment:
%       plot([0:0.01:50], sigma_dlr_fall([0:0.01:50], 2))
%

sig = zeros(size(val));

if segType == 1         % DOWN
    % val <= 25
    ind1 = (val<=25);
    sig(ind1) = polyval([0.000 0.140 0.500], val(ind1));
    % val > 25
    ind2 = (val>25);
    sig(ind2) = 4.00;
elseif segType==2       % CONSTANT
    % val <= 2
    ind1 = val<=2;
    sig(ind1) = 0.090;
    % 2 < val <=25
    ind2 = (val>2 & val<=25);
    sig(ind2) = polyval([0.000 0.135 -0.180], val(ind2));
    % val > 25
    ind3 = (val>25);
    sig(ind3) = 3.195;
elseif segType==3       % UP
    % val <= 25
    ind1 = (val<=25);
    sig(ind1) = polyval([0.000 0.140 0.100], val(ind1));
    % val > 25
    ind2 = (val>25);
    sig(ind2) = 3.600;
end



% winter
%---------------------------------------------------------------------------

function mu = mu_dlr_winter(val, segType)
% function mu = mu_dlr_winter(val, segType)
% mu        : mean value
% val       : value (vector or single value) to evaluate the mean
% segType   : segment type (1, 2 or 3 for D, C, or U)
%
% Definition of the analytical function for mean at the DLR, 
% season "winter".
%
% This is defined as a real Matlab function like cos or exp since the 
% input "val" can be a vector or a single value.
% So the function can be plotted. For example for a Down segment:
%       plot([0:0.01:50], mu_dlr_winter([0:0.01:50], 2))
%
% The polynom is defined by:
%   mu = x0 + x1.val + x2.val^2
% for each segment type D, C and U
% 
% pMn mapping:
% pMn = [   D2    C2    U2
%           D1    C1    U1
%           D0    C0    U0  ]
%

pMu = [ -0.01300   -0.00070    0.00000;
         1.15000    0.99150    0.83000;
        -0.50000    0.01800    0.35000];
     
mu = polyval(pMu(:,segType), val);



function sig = sigma_dlr_winter(val, segType)
% function sig = sigma_dlr_winter(val, segType)
% sig       : standard deviation value
% val       : value (vector or single value) to evaluate the standard deviation
% segType   : segment type (1, 2 or 3 for D, C, or U)
%
% Definition of the analytical function for standard deviation at the
% DLR, season "winter".
%
% This is defined as a real Matlab function like cos or exp since the 
% input "val" can be a vector or a single value.
% So the function can be plotted. For example for a Down segment:
%       plot([0:0.01:50], sigma_dlr_winter([0:0.01:50], 2))
%

sig = zeros(size(val));

if segType == 1     % DOWN
    % val <= 25
    ind1 = (val<=25);
    sig(ind1) = polyval([0.000 0.150 0.100], val(ind1));
    % val > 25
    ind2 = (val>25);
    sig(ind2) = 2.650;
elseif segType==2       % CONSTANT
    % val <= 2
    ind1 = (val<=2);
    sig(ind1) = 0.090;
    % 2 < val <= 25
    ind2 = (val>2 & val<=25);
    sig(ind2) = polyval([0.000 0.120 -0.150], val(ind2));
    % val > 25
    ind3 = (val>25);
    sig(ind3) = 2.85;
elseif segType==3           % UP
    % val <= 25
    ind1 = (val<=25);
    sig(ind1) = polyval([0.000 0.105 0.150], val(ind1));
    % val > 25
    ind2 = (val>25);
    sig(ind2) = 2.775;
end