电离层误差模型的matlab代码

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