xgeoidh.m

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

%                             xgeoidh.m
%  Scope:   This MATLAB program computes the WGS-84 geoid height correction for
%           the current location specified by latitude and longitude. Optional,
%           the worldwide WGS-84 geoid height correction contour map is plotted.
%  Usage:   xgeoidh
%  Inputs:  - name of the output file containing the results
%           - default - WGS-84 geoid heights data file (tgeoid84.dat)
%           - latitude in radians,  - pi/2 to pi/2
%           - longitude in radians,  0 to 2*pi (automatically adjustable 
%             -2*pi to 4*pi)
%           - selection of the contour map (optional)
%  Outputs: - input/output data stored on the selected output file or
%             displayed on screen; geoid height correction in meters 
%           - plot of the worldwide WGS-84 geoid height correction contour map
%  Remark:  The program can be easily modified to enter a generic geoid 
%           height correction earth model instead of the WGS-84 earth model
%           (default).
%  Reference:
%           [1] Anonym., Department of Defense World Geodetic System 1984.
%               DMA Technical Report 8350.2, Second Edition, September 1, 1991.
%  External Matlab macros used:  geoidh 
%  Last update:  11/02/00
%  Copyright (C) 1997-00 by LL Consulting. All Rights Reserved. 

clear
close all

yes = 'y';

%  Save the generated data into a specified file if desired

disp('  ');
answer2 = input('Do you want to save the generated data? (y/n)[y] ','s');
if  isempty(answer2)
   answer2 = 'y';
end
if strcmp(answer2,yes) == 1
   f2 = input('Specify the output file name  --> ','s');
else
   f2 = 1;        %  output to the screen
end
disp('  ');

%  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.

fname = 'tgeoid84.dat'; % it can be changed to another geoid heights data file
tgeoid = load(fname);

answer = 'y';
while (strcmp(answer,yes) == 1)
 
%  Enter latitude and longitude

   flag = 0;
   while (flag == 0)
      lat = input('Enter latitude in radians ( -pi/2 to pi/2 ) -->  ');
      if  ( lat >= -pi/2 ) & ( lat <= pi/2 )
          flag = 1;
      end
   end 
   flag = 0;
   while (flag == 0)
      lon = input('Enter longitude in radians  ( 0 to 2*pi )   -->  ');
      if  ( lon >= 0. ) & ( lon <= 2*pi )
          flag = 1;
      end
   end 
   disp('  ');
   
%  Compute geoid height correction

   xoutput  = geoidh(lat,lon,tgeoid);

%  Save input data and results into an external file or display on screen

   fprintf(f2,'\n********************************************************\n\n');
   fprintf(f2,'***  COMPUTATION OF WGS-84 GEOID HEIGHT CORRECTION  ****');
   fprintf(f2,'\n\n*****   Input data   ****** \n\n');
   fprintf(f2,'Latitude in radians  --> %16.7f\n',lat);
   fprintf(f2,'Longitude in radians --> %16.7f\n',lon);
   fprintf(f2,'\n*****   Output data   ***** \n\n');
   fprintf(f2,'Geoid height correction in meters = %12.5f\n',xoutput);  
   fprintf(f2,'\n********************************************************\n\n');

%  Select another computation, if desired

   answer = input('Do you want another computation? (y/n)[n] ','s');

end

%  Plot the worldwide contour map for the WGS-84 geoid height correction

disp('  ');
answer1 = input('Do you want to map the WGS-84 geoid height correction? (y/n)[y] ','s');
if  isempty(answer1)
   answer1 = 'y';
end
if strcmp(answer1,yes) == 1
   xlongrid = 0:10:360;
   ylatgrid = -90:10:90;
   v = -60:10:60;    %  10 degrees - contour lines
   tgeoid = [tgeoid  tgeoid(:,1)];
   figure(1)
   cc = contour(xlongrid,ylatgrid,tgeoid,v);
   xlabel('Longitude, in degrees');
   ylabel('Latitude, in degrees');
   title('WGS-84 Geoid Height Correction (in meters)');
   clabel(cc);
end

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