function x = ionoKlobuchar(Latu, Lonu, Az, El, Ttr, klob_coeff)
%Computation of the ionospheric correction due to Klobuchar broadcast 
%ionospheric coefficients for observation on GPS L1 frequency

%INPUT:
%Latu: geodetic latitude [dgree]  
%Lonu: geodetic longitude [dgree]
%Az: Azimuth [dgree]
%El: Elevation [dgree]
%Ttr: Time of observation [s], seconds of week
%klob_coeff: Matrix (2x4), alpha (1st row) and beta (2nd row) Klobuchar 
%               coefficients from navigation message
 
%OUTPUT:
%x: ionospheric correction [m]

%Copyright (c) 2023, M.F. Glaner
%Adapted by GCC Group
%--------------------------------------------------------------------------

C=299792458;

alpha = klob_coeff(1,:);
beta  = klob_coeff(2,:);

% Convert angles from dgree into semicircles
Latu = Latu/180;
Lonu = Lonu/180;
Az = Az/180;
El = El/180;

% Broadcast ionospheric Klobuchar coefficients
a0 = alpha(1,1);
a1 = alpha(1,2);
a2 = alpha(1,3);
a3 = alpha(1,4);
b0 = beta(1,1);
b1 = beta(1,2);
b2 = beta(1,3);
b3 = beta(1,4);

psi = 0.0137/(El + 0.11) - 0.022;       % Earth-centered angle [sc]
Lati = Latu + psi*cos(Az*pi);           % Subionospheric latitude [sc]
if Lati > 0.416
    Lati = 0.416;
else
    if Lati < -0.416
        Lati = -0.416;
    end
end
Loni = Lonu + psi*sin(Az*pi)/cos(Lati*pi);  % Subionospheric longitude [sc]
Latm = Lati + 0.064*cos((Loni -1.617)*pi);  % Geomagnetic latitude of the ionosphere intersection point [sc]
Ttr = mod(Ttr, 86400);                      % seconds of day
T = 43200*Loni + Ttr;                       % Local time at subionospheric point [sec]
if T >= 86400
    T = T - 86400;
else
    if T < 0
        T = T +86400;
    end
end
F = 1 + 16*(0.53 - El)^3;                       % Slant Factor
per = b0 + b1*Latm + b2*Latm^2 + b3*Latm^3;     % Period of model
if per < 72000
    per = 72000;
end
amp = a0 + a1*Latm + a2*Latm^2 + a3*Latm^3;     % Amplitud of the model
if amp < 0
    amp = 0;
end
x = 2*pi*(T - 50400)/per;       % Phase of the model (Max at 14.00 = 50400 sec local time)
if abs(x) > 1.57                % ionospheric correction
    diono = F*5*10^-9*C;
else
    diono = F*(5*10^-9 + amp*(1 - 0.5*(x^2) + x^4/24))*C;
end

x = diono;          % ionospheric delay in [m]
end

