Code covered by the BSD License

# Geodetic Toolbox

### Mike Craymer (view profile)

13 Jun 2007 (Updated )

Toolbox for angle, coordinate and date conversions and transformations. Version 2.97.

DirInv.m
```%echo on
%--------------------------------------------------
% DirInv
%   Script to perform direct or inverse geodetic
%   problem.
% 01 May 96
%
%                    refell
%--------------------------------------------------
%clear

%----- Get type of problem

'Direct',...
'Inverse');
if (iprob==1)
elseif (iprob==2)
else
error('Undefined geodetic problem');
end

%----- Get reference ellipsoid

'WGS84 (GPS)',...
if (iref==1)
[a,b,e2,finv]=refell('GRS80');
elseif (iref==2)
[a,b,e2,finv]=refell('WGS84');
elseif (iref==3)
[a,b,e2,finv]=refell('CLK86');
end

%----- Get lat,lon,ht of 1st (from) points

disp(' ');
disp('Enter matrix of lats (D M S), longs (D M S), hts (m) in []:');
data=input('');
h1=data(7);
n1=length(h1);

%------------------------------------------------------
% For direct problem
%------------------------------------------------------

if (iprob==1)

%----- Get az,va,dist to 2nd points
disp(' ');
disp('Enter matrix of az (D M S), va (D M S), dist (m) in []:');
data=input('');
d=data(7);
n2=length(d);
if (n1~=n2)
error('Different numbers of from points and obs');
end

%----- Compute coords of 2nd points
[lat2,lon2,h2]=direct(lat1,lon1,h1,az,va,d,a,e2);

%----- List results
for i=1:n1
fprintf(1,'\nTo Point   %4.0f\n',i);
fprintf(1,'Latitude   %4.0f %2.0f %9.6f\n',latdms(i,1),latdms(i,2),latdms(i,3));
fprintf(1,'Longitude  %4.0f %2.0f %9.6f\n',londms(i,1),londms(i,2),londms(i,3));
fprintf(1,'Height     %9.4f\n',h2(i));
end

%------------------------------------------------------
% For inverse problem
%------------------------------------------------------

elseif (iprob==2)

%----- Get coordinates of 2nd points
disp(' ');
disp('Enter matrix of lats (D M S), longs (D M S), hts (m) in []:');
data=input('');
h2=data(7);
n2=length(h2);
if (n1~=n2)
error('Different number of from and to points');
end

%----- Compute inverse
[az,va,d]=inverse(lat1,lon1,h1,lat2,lon2,h2,a,e2);

%----- List results
for i=1:n1
fprintf(1,'\nLine      %3.0f\n',i);
fprintf(1,'Azimuth   %3.0f %2.0f %9.6f\n',azdms(i,1),azdms(i,2),azdms(i,3));