ionoc.m

GPS TOOLBOX包含以下内容： 1、GPS相关常量和转换因子； 2、角度变换； 3、坐标系转换： &#61656 点变换； &#61656 矩阵变换； &#61656 向量变换

%                           ionoc.m
%  Scope:   This MATLAB macro computes L1 iono correction for a specified user 
%           by using the Klobuchar model; WGS-84 constants are used.
%  Usage:   ionocorr = ionoc(lat,lon,el,az,tgps,alpha,beta)
%  Description of parameters:
%           lat      - input, user's geodetic latitude, in radians
%           lon      - input, user's geodetic longitude, in radians
%           el       - input, user's elevation angle, in radians
%           az       - input, user's azimuth angle, in radians
%           tgps     - input, GPS system time, in seconds
%           alpha    - input, coefficients of cubic fit to the amplitude of 
%                      vertical delay (array with 4 elements)
%           beta     - input, coefficients of cubic fit to the period of model
%                      (array with 4 elements)
%           ionocorr - output, L1 iono correction, in meters
%  Remark:  L2 iono correction can be determined as follows
%               L2 Iono correction = (L1 iono correction) * ((77/60)**2)
%  References:
%           [1] Parkinson, B. W., Spilker, J. J., Editors, Global Positioning 
%               System: Theory and Applications. AIAA, 1996 (Chapter 12, 
%               Ionospheric effects on GPS by J. A. Klobuchar).
%           [2] Farrell, J., Barth, M., The Global Positioning System and 
%               Inertial Navigation. McGraw Hill, 1998, Appendix E.
%  External Matlab macros used:  wgs84con
%  Last update:  07/22/00
%  Copyright (C) 1999-00 by LL Consulting. All Rights Reserved.

function   ionocorr = ionoc(lat,lon,el,az,tgps,alpha,beta)

wgs84con;
% global constants used: c_speed

%  Calculate angles in semicircles

elsm = el / pi;
azsm = az / pi;
latsm = lat / pi;
lonsm = lon / pi;

%  Compute the earth-centered angle

psi = 0.0137 / (elsm + 0.11) - 0.022;

%  Compute the subionospheric latitude

iono_lat = latsm + psi * cos(az);
if (iono_lat > 0.416)
   iono_lat = 0.416;
elseif (iono_lat < -0.416) 
   iono_lat = -0.416;
end

%  Compute the subionospheric longitude

iono_lon = lonsm + (psi * sin(az)) / (cos(iono_lat * pi));

%  Find the local time at the subionospheric point

localt = 43200. * iono_lon + tgps;
kk = 0;
while ( (localt >= 86400.) & (kk < 10) )
  localt = localt - 86400.;
  kk = kk + 1;
end
while ( (localt < 0.) & (kk < 10) )
  localt = localt + 86400.;
  kk = kk + 1;
end
if (kk == 10) 
   error('IONOC.m - error in local time computation');
end

%  Calculate the geomagnetic latitude of the earth projection of the 
%  ionospheric intersection point 

latm = iono_lat + 0.064 * cos((iono_lon - 1.617) * pi);

%  Calculate the period by using the beta terms from the GPS message

per = ((beta(4)*latm + beta(3))*latm + beta(2))*latm + beta(1); 
if (per < 72000.) 
   per = 72000.;
end

%  Calculate the argument of the cosine term

x = (pi + pi) * (localt - 50400.) / per;

%  Calculate the slant factor

slantf = 0.53 - elsm;
slantf = 1. + 16. * slantf * slantf * slantf;

%  Calculate the amplitude by using the alpha terms from the GPS message

amp = ((alpha(4)*latm + alpha(3))*latm + alpha(2))*latm + alpha(1); 
if (amp < 0.) 
   amp = 0.;
end

%  Calculate the L1 ionospheric correction (in meters)

if (abs(x) < 1.57)
   xx = x*x;
   ionocorr = slantf * (5.e-9 + amp * ( 1. - 0.5*xx + (xx*xx/24) ) ) * c_speed;
else
   ionocorr = slantf * 5.e-9 * c_speed;
end