clear

close all

yes = 'y';

% Initialize the iono constants alpha and beta

ionocon

deg_to_rad = pi / 180.; % degrees to radians transformation

% Read the input file containing satellites data (almanac)

disp(' ');

disp('Enter almanac data - the default is data file wk749.dat');

answer1 = input('Do you want to use the default data file? (y/n)[y] --> ','s');

if isempty(answer1)

answer1 = yes;

end

disp(' ');

if (strcmp(answer1,yes) == 1)

ff = 'wk749.dat';

else

ff = input('Specify the input filename (with extension) --> ','s');

disp(' ');

end

almanac = load(ff);

[npt,nmax] = size(almanac);

nr_sat = npt; % number of satellites in the almanac

ind_sv = almanac(:,1); % id of the svs in the almanac

% Read the geographic location for the test

disp('Enter geographic location data - the default is the data file

locat1.dat');

answer2 = input('Do you want to use the default data file? (y/n)[y] ','s');

if isempty(answer2)

answer2 = yes;

end

if (strcmp(answer2,yes) == 1)

fff = 'locat1.dat';

hh = load(fff);

[npt,nmax] = size(hh);

lat = hh(1); % only the first point is selected

lon = hh(2);

alt = hh(3);

else

lat = input('Enter latitude in degrees --> ');

lon = input('Enter longitude in degrees --> ');

alt = input('Enter altitude in meters --> ');

end

lat_rad = lat * deg_to_rad; % in radians

lon_rad = lon * deg_to_rad; % in radians

% Compute the geographic locations in ECEF

loc_xyz = tgdecef(lat_rad,lon_rad,alt);

% Specify the elevation angle limit

disp(' ');

disp('Enter elevation angle limit - the default value is 5 degrees');

answer3 = input('Do you want to use the default value? (y/n)[y] ','s');

if isempty(answer3)

answer3 = yes;

end

if (strcmp(answer3,yes) == 1)

ele_limd = 5.;

else

ele_limd = input('Specify the elevation angle limit (in degrees) --> ');

end

ele_lim = ele_limd * deg_to_rad; % elevation angle limit in radians

% Select number of time samples and time step

disp(' ');

disp('Enter number of time samples - the default value is 10 ');

answer5 = input('Do you want to use the default value? (y/n)[y] ','s');

if isempty(answer5)

answer5 = yes;

end

if (strcmp(answer5,yes) == 1)

t_sample = 10; % default - number of time samples

else

t_sample = input('Enter number of time samples --> ' );

end

t_incr = 0.;

if t_sample > 1

disp(' ');

disp('Enter time step - the default value is 300. ');

answer6 = input('Do you want to use the default value? (y/n)[y] ','s');

if isempty(answer6)

answer6 = yes;

end

if (strcmp(answer6,yes) == 1)

t_incr = 300.; % default - time step

else

t_incr = input('Enter time step, in seconds --> ' );

end

end

disp(' ');

% Start main computation loop - time loop

nr_vis = zeros(1,t_sample);

t_iono = zeros(nr_sat,t_sample);

vis_sat = zeros(t_sample,nr_sat);

for kkk = 1:t_sample % cycle for all time samples

fprintf('***** Computation in progress for time sample = %4.0f\n',kkk);

% Determine satellite position

for k = 1:nr_sat, % cycle for all satellites

temp = almanac(k,2:9);

psat_temp = (svpalm(tsim,temp))';

psat(1:3,k) = psat_temp;

end

% Compute elevation and azimuth angles

for k = 1:nr_sat, % cycle for all satellites

temp2 = psat(:,k);

[temp3,temp1,ulos,range] = elevaz(loc_xyz,temp2);

eangle(k) = temp3;

aangle(k) = temp1;

end

% Determine the number of visible satellites and the iono correction

% for each visible satellite; if the satellite is not visible the iono

% correction is set to 0.

for k = 1:nr_sat,

if eangle(k) >= ele_lim

nr_vis(kkk) = nr_vis(kkk) + 1;

vis_sat(k,kkk) = k;

el = eangle(k);

az = aangle(k);

ionocorr = ionoc(lat_rad,lon_rad,el,az,tsim,alpha,beta);

t_iono(k,kkk) = ionocorr;

else

vis_sat(k,kkk) = 0;

t_iono(k,kkk) = 0.;

end

end

tsim = tsim + t_incr;

end % end of the time loop (index kkk)

% Save the iono correction data (optional)

disp(' ');

answer7 = input('Do you want to save the iono correction data? (y/n)[n]

','s');

if isempty(answer7)

answer7 = 'n';

end

if (strcmp(answer7,yes) == 1)

disp(' ');

answer8 = input('Do you want to use the default file, xionoc1.out? (y/n)[y]

','s');

if isempty(answer8)

answer8 = yes;