function varargout = daspectm(zunits,vfac,lat,lon,az,angleunits,radius)

%DASPECTM figure DataAspectRatio property for a map  
%
%     DASPECTM('zunits') sets the figure DataAspectRatio property so
%     that the z axis is in proportion to the x and y projected coordinates.
%     This permits elevation data to be displayed without vertical distortion.
%     The string 'zunits' specifies the units of the elevation data, and can
%     be any string recognized by DISTDIM.
%
%     DASPECTM('zunits',vfac) sets the DataAspectRatio property so
%     that the z axis is vertically exaggerated by the factor vfac. If
%     omitted, the default is no vertical exaggeration.
%     
%     DASPECTM('zunits',vfac,lat,long) sets the aspect ratio based 
%     on the local map scale at the specified geographic location. 
%     If omitted, the default is the center of the map limits.
%     
%
%     DASPECTM('zunits',vfac,lat,long,az) also specifies the direction 
%     along which the scale is computed. If omitted, 90 degrees (west) 
%     is assumed.
%
%     DASPECTM('zunits',vfac,lat,long,az,'gunits') also specifies the 
%     units in which the geographic position and direction are given. 
%     If omitted, 'degrees' is assumed.
%
%     DASPECTM('zunits',vfac,lat,long,az,'gunits',radius) 
%     uses the last input to determine the radius of the sphere.  If radius is 
%     a string, then it is evaluated as an ALMANAC body to determine the 
%     spherical radius.  If numerical, it is the radius of the desired sphere 
%     in zunits. If omitted, the default radius of the Earth is used.
%
%     See also  DASPECT, PAPERSCALE

%  Copyright 1996-2000 Systems Planning and Analysis, Inc. and The MathWorks, Inc.
%  Written by: W. Stumpf, A. Kim, T. Debole

if nargin < 1 | nargin > 7 | nargin == 3
	error('Incorrect number of arguments')
end

% fill in values of optional arguments

if nargin < 7; radius = 'earth'; end
if nargin < 6; angleunits = 'degrees';end
if nargin < 5; az = 90;end
if nargin < 4;     
	lat = mean(getm(gca,'mapLatLim'));
    lon = mean(getm(gca,'mapLonLim'));
    angleunits = getm(gca,'angleunits');
    lat = angledim(lat,angleunits,'deg');
    lon = angledim(lon,angleunits,'deg');
end
if nargin < 2; vfac = 1;end

% check for globe

if strmatch(getm(gca,'mapprojection'),'globe')
	error('Aspect ratio changes not appropriate for Globe')
end	


% Cartesian coordinates of starting point

[xo,yo] = mfwdtran(lat,lon); 

% Geographical and paper coordinates of a point downrange

sdistdeg = distdim(1,zunits,'deg',radius);

[nlat,nlon] = reckon(lat,lon,sdistdeg,az);
[xn,yn] = mfwdtran(nlat,nlon); % cartesian coordinates

% Distance between two points in cartesian coordinats

cdist = sqrt((xo-xn)^2+(yo-yn)^2);

% Set the DataAspectRatio

zratio = 1/(cdist*vfac);
dataspectratios = [1 1 zratio];
set(gca,'DataAspectRatio',dataspectratios)	% set dataaspectratio

