xtropoc1.m

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

%                                xtropoc1.m
%  Scope:  This Matlab program determines tropospheric delay contours for a 
%          specified location (latitude/longitude), elevation range, and 
%          altitude range by using a simplified model [1]; WGS-84 constants are 
%          used.
%  Usage:  xtropoc1
%  Inputs: - WGS-84 geoid heights data file, by default tgeoid84.dat
%          - location (latitude/longitude), in radians
%          - other parameters are initialized by default
%  Outputs: - plot of tropo corrections (contour grapg) for the specified location
%  Reference:
%          [1] Farrell, J., Barth, M., The Global Positioning System and 
%              Inertial Navigation. McGraw Hill, 1998.
%  External Matlab macros used: geoidh, tropoc1
%  Last update:  11/15/00
%  Copyright (C) 1999-00 by LL Consulting. All Rights Reserved.

clear 
close all
yes = 'y';
answer1 = yes;

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

ang_incr = pi / 18.;    %  angle increment, in radians    (default)
alt_incr = 1000.;       %  altitude increment, in meters  (default)

%  Load the WGS-84 geoid heights based on WGS-84 geoid correction earth model.       
%  The table tgeoid84.dat contains WGS-84 geoid heights (in meters) for a 19 
%  by 36 Latitude/Longitude grid; the latitudes go northward from -90 to 90 
%  degrees, and longitudes go eastward from 0 to 350 degrees.

f1 = 'tgeoid84.dat';    %  default file
tgeoid = load(f1);

disp('  ');
kfig = 1;
while strcmp(answer1,yes)

%  Read the geographic location for the test

   lat = input('Enter latitude, in degrees  -->  ');
   lon = input('Enter longitude, in degrees -->  ');
   disp('  ');  
   lat_rad   = lat * deg_to_rad;         %  in radians
   lon_rad   = lon * deg_to_rad;         %  in radians

%  Start main computation loop 

   xx = zeros(10,11);

   for kk = 0:1:9              %  cycle for all elevation steps
 	  elev = ang_incr*kk;       %  elevation angle, in radians
     for k = 0:1:10            %  cycle for all altitude steps
         alt = alt_incr*k;     %  altitude above ellipsoid, in meters
         temp = tropoc1(lat_rad,lon_rad,alt,elev,tgeoid);
         xx(kk+1,k+1) = temp;
      end
   end

   figure(kfig)
   [c,cs] = contour(xx,[10 4 3 2 1.5 1 0.8]);
   clabel(c); grid
   set(cs,'MarkerSize',.2,'LineWidth',1);    
   h1 = gca;
   set(h1,'XLim',[1 11],'XTick',[1:2:10],...
          'XTickLabel',['    0.';' 2000.';' 4000.';' 6000.';' 8000.'])
   set(h1,'YLim',[1 10],'YTick',[1:1:10],...
          'YTickLabel',[' 0';'10';'20';'30';'40';'50';'60';'70';'80';'90'])
   set(get(gca,'xlabel'),'String',(['Altitude, in meters']),...       
          'FontName','Helvetica','FontSize',12)
   set(get(gca,'ylabel'),'String',('Elevation Angle, in degrees'),....    
          'FontName','Helvetica','FontSize',12)
   temp1 = [num2str(lat) ' deg.'];
   temp2 = [num2str(lon) ' deg.'];
   name = ['Tropo Delay Contours (meters), for lat=' temp1 ',lon=' temp2];
   set(get(h1,'Title'),'String',name,'FontName','Helvetica','FontSize',12);
      
   answer1 = input('Do you want to execute another computation/plot? (y/n)[n] --> ','s');
   if isempty(answer1)
      answer1 = 'n';
   end
   kfig = kfig + 1;
   disp('  ');
end

disp('End of the program  XTROPOC1');
disp('  ');