clearvars; clc;

%--------------------------------------------------------
%--------------------------------------------------------
alpha = [1.8491D-08  2.4361D-08 -4.0024D-07  6.1181D-07];
beta =  [1.3108D+05  1.5212D+05 -5.5311D+05  3.8104D+06];

%alpha = [7.7941D-09 -7.9182D-09 -2.1120D-07 -2.8892D-07];
%beta =  [9.0339D+04 -6.6516D+04 -9.8950D+05 -4.3387D+06];

% Time (UTC)
year = 2014;
month = 6;
day = 29;
hour = 0;
MIN = 0;
sec = 0;  
elev = 90;                        % Elevation (deg)
azimuth = 0;                      % Azimuth   (deg)
%--------------------------------------------------------
%--------------------------------------------------------  

lat_rng = -87.5 : 87.5;               % Latitude  (deg)
lon_rng = -180 : 180;             % Latitude  (deg)
              
i = 1;
for lat = lat_rng
    j = 1;
    for lon = lon_rng

        % Compute Time-of-Week in seconds (UTC)
        [week, tow] = cal2gpstime(year, month, day, hour, MIN, sec);

        % Ionospheric L1 correction
        dIon1(i,j) = klobuchar(lat, lon, elev, azimuth, tow, alpha, beta);
        
        j = j + 1;
    end
    i = i + 1;
end

[M, c] = contourf(lon_rng, lat_rng, dIon1, 100, 'LineColor','none');
set(gcf, 'Position',[100 100 1200 600])
axis equal
colormap('jet')
cb = colorbar;
caxis([0, 10]);

xlabel('longitude [°]', 'FontSize',12)
ylabel('latitude [°]', 'FontSize',12)
ylabel(cb, 'Ionospheric L1 correction [m]', 'FontSize',12)
title(sprintf('Global Ionospheric Map - Klobuchar @ %d-%02d-%02d %02d:%02d UTC', year, month, day, hour, MIN), 'FontSize',14)


