xgeodes.m

GPS software toolbox for GPS receiver development

%                               xgeodes.m
%  Scope:   This MATLAB program determines geodesic when the departure and
%           destination points are specified [1]; WGS-84 earth model is used.
%  Usage:   xgeodes
%  Inputs:  - latitude and longitude data for the departure and destination 
%             points; default data are given. The data are specified in the 
%             following order (each row represents a departure-destination
%             pair):
%             - lat1, departure point latitude, in degrees
%             - lon1, departure point longitude, in degrees
%             - lat2, departure point latitude, in degrees
%             - lon2, departure point longitude, in degrees
%           - name of the input datafile (optional)
%           - departure-destination data (optional) 
%           - name of the output file (optional)
%           - selection of the longitude-latitude plot of the selected 
%             waypoints (optional)
%  Outputs: - input/output data stored on the selected output file or
%             displayed on screen. For each pair of departure-destination
%             point the range, departure bearing and destination bearing
%             are determined
%           - plot of the longitude-latitude of the specified waypoints (with
%             geodesic distance) - optional
%  Remark:  The program can be easily modified to input latitude/longitude data 
%           in radians or degrees/minutes/seconds.
%  Reference:
%           [1] Anon., Minimum operation performance standards for airborne
%               supplemental navigation equipment using global positioning
%               system (GPS). RTCA/DO-208, July 1991, Appendix B.
%  External Matlab macros used:  convcon, geodes, wgs84con
%  Last update:  11/07/00
%  Copyright (C) 1998-00 by LL Consulting. All Rights Reserved.

clear 
close all
yes = 'y';
icase = 0;

% global flatness b_smaxis eprime2
wgs84con;
% deg_rad
convcon;

%  Select the input/output media

disp('  ');
answer1 = input('Do you want to use the default data? (y/n)[n] --> ','s');
disp('  ');
if (strcmp(answer1,yes) == 1)
   disp('Enter data from:  1) data from input data file  xgeodes1.dat');
   disp('                  2) LAX to JFK airports'); 
   answer2 = input('Make selection (1 or 2) --> ');
   disp('  ');
   if  (answer2 == 1)
      finp ='xgeodes1.dat';
   else
      tt(1,1) =   33.9422222222;      %  LAX latitude, in degrees
      tt(1,2) = -118.4072222222;      %  LAX longitude, in degrees
      tt(1,3) =   40.64;              %  JFK latitude, in degrees
      tt(1,4) =  -73.7780555555;      %  JFK longitude, in degrees
      icase = 1;   
   end   
else
   answer3 = input('Do you want to enter the data from keyboard? (y/n)[n] --> ','s');
   disp('  ');
	if (strcmp(answer3,yes) == 1)
      tt(1,1) = input('Enter departure latitude, in degreees --> ');
      tt(1,2) = input('Enter departure longitude, in degreees --> ');
      tt(1,3) = input('Enter destination latitude, in degreees --> ');
      tt(1,4) = input('Enter destination longitude, in degreees --> ');
      icase = 1;
   else
      finp = input('Specify the input data filename, e.g. xgeodes1.dat -->   ','s');
   end
   disp('  ');
end

if (icase == 0)
	tt = load(finp);
	[nrow,ncol] = size(tt);
   if  ncol ~= 4
      disp('XGEODES - Error: check the input data (number of columns)');
      break
	end
end

answer4 = input('Do you want to save the generated data? (y/n)[y] --> ','s');
disp('  ');
if  isempty(answer4)
   answer4 = 'y';
end
if strcmp(answer4,yes) == 1
   fout = input('Specify the output filename --> ','s');
   disp('  ');
else
   fout = 1;     %   output to the screen
end
[nrow,ncol] = size(tt);

%  Determine the geodetic and save the data for each set of data 

range_total = 0.;
x = zeros(2*nrow,1);
y = zeros(2*nrow,1);
kk = 1;
for k = 1:nrow
   
   lat1 = tt(k,1);
   lon1 = tt(k,2);
   lat2 = tt(k,3);
   lon2 = tt(k,4);
   x(kk) = lon1;
   x(kk+1) = lon2;
   y(kk) = lat1;
   y(kk+1) = lat2;
   kk = kk + 2;
   
%  Determine the geodetic data
   
   [range,bearing1,bearing2] = geodes(lat1,lon1,lat2,lon2);

%  Save the generated data into a specified file or displayed on screen

   fprintf(fout,'************************************************************\n\n');
   fprintf(fout,'*****Input data %2.0f\n',k);
   fprintf(fout,'Departure latitude    = %22.14e degrees\n',tt(k,1));
   fprintf(fout,'Departure longitude   = %22.14e degrees\n',tt(k,2));
   fprintf(fout,'Destination latitude  = %22.14e degrees\n',tt(k,3));
   fprintf(fout,'Destination longitude = %22.14e degrees\n',tt(k,4));
   fprintf(fout,'\n\n');
   fprintf(fout,'*****Output data %2.0f\n',k);
   fprintf(fout,'Range                 = %22.14e meters\n',range);
   fprintf(fout,'Departure bearing     = %22.14e degrees\n',bearing1);
   fprintf(fout,'Destination bearing   = %22.14e degrees\n\n',bearing2);
   fprintf(fout,'************************************************************\n\n');
   
   range_total = range_total + range;
   
end

answer5 = input('Do you want to plot the waypoints? (y/n)[y] --> ','s');
disp('  ');
if  isempty(answer5)
   answer5 = 'y';
end
if strcmp(answer5,yes) == 1
   
%  Execute the plot containing the waypoints; the total range is displayed  

   temp = ['Total range = ' num2str(range_total) ' meters' ];
   disp('Place the text with total range on the graph; ');
   tempp = input('      Press <Enter> or <Return> key to start ... ');
   figure(1)
   plot(x,y), grid,...
   title('Longitude - Latitude plot of the selected waypoints (with geodesic distance)');
   xlabel('Longitude, in degrees');
   ylabel('Latitude, in degrees');
   hold on
   plot(x,y,'*');
   gtext(temp)           %  mouse placement of the text on a graph
   disp('   ');
   temp = input('Press <Enter> or <Return> key to continue ... ');
   disp('   ');

end

disp('End of the program  XGEODES.M');
disp('  ');