电离层误差模型的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
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;
end
if (strcmp(answer8,yes) == 1)
fout = 'xionoc1.out'; % default file name
else
disp(' ');
fout = input('Enter the name of the output file (with extension) --> ','s');
end
for k = 1:t_sample
for k = 1:nr_sat
fprintf(fout,'%10.5f ',t_iono(k,kkk));
end
fprintf(fout,'\n');
end
end
% Execute the plots
time = 1:t_sample;
ind_sv_max = max(ind_sv);
t_iono_sv = zeros(ind_sv_max,t_sample);
ind_sv_sv = zeros(1,ind_sv_max);
for k = 1:nr_sat
t_iono_sv(ind_sv(k),:) = t_iono(k,:);
ind_sv_sv(ind_sv(k)) = ind_sv(k);
end
% Number of visible satellites versus time sample
if (t_sample > 1)
figure(1)
plot(time,nr_vis,'*')
grid, axis([1 t_sample 3 12]);
aa =['Time sample number (time step = ' num2str(t_incr) ' seconds, ']; aa = [aa 'start time tow = ' num2str(tweek) ' seconds)'];
xlabel(aa);
ylabel('Number of visible satellites');
cc = 'Number of visible satellites versus Time sample number (';
cc = [cc 'almanac ' ff ')'];
title(cc);
else
fprintf('\nNumber of visible satellites = %3.0f\n',nr_vis(1));
end
% Iono correction versus time sample number and satellite id
if (t_sample > 1)
figure(2)
surf(t_iono_sv);
axis([0 t_sample 0 ind_sv_max 0 16.]);
xlabel('Time sample number');
ylabel('Satellite number');
zlabel('Iono correction, in meters');
title('Iono correction versus Time sample and Satellite number');
else
for k = 1:nr_sat
if (t_iono(k,1) ~= 0.0)
fprintf('\nIono correction for satellite %3.0f is %11.4f
meters\n',ind_sv(k),t_iono(k,1));
end
end
end
% Iono correction versus satellite id for a selected time sample
fprintf('\nNumber of time sample available = %3.0f\n',t_sample);
sel_sample = input('Select the time sample number to be analyzed --> '); figure(3)
plot(ind_sv,t_iono(:,sel_sample),'*');
grid, axis([1 ind_sv_max -1 16]);
xlabel(['Satellite number (almanac ' ff ')']);
ylabel('Iono correction, in meters (0 if satellite is not visible)'); title(['Iono correction versus Satellite number, for Time sample # '
num2str(sel_sample)]);
% Iono correction versus time sample for a selected satellite
fprintf('\nThe following satellites are in the almanac:\n');
for k = 1:nr_sat
fprintf('%3.0f',ind_sv(k));
end
disp(' ');
sel_sat = input('Select the satellite id --> ');
if (t_sample > 1)
figure(4)
plot(time,t_iono_sv(sel_sat,:));
hold on
plot(time,t_iono_sv(sel_sat,:),'*');
grid
xlabel(aa);
ylabel('Iono correction, in meters');
cc = ['Iono correction versus Time sample number, for satellte # ' num2str(sel_sat)];
cc = [cc ' (almanac ' ff ')'];
title(cc);
else
fprintf('\nIono correction for satellite %3.0f is %11.4f
meters\n',sel_sat,t_iono_sv(sel_sat,1));
end
disp(' ');
disp('End of the program XIONOC');
disp(' ');。