      function [ t1 p1 ] = givloc(lat1,lon1,del,az)
% Given a starting location on the Earth, find the point on the Earth 
% at the specified distance and azimuth from that starting point.
%
% input:
%   lat1,lon1 = starting latitude and longitude (degrees)
%   del      = distance to target, in degrees
%   az       = azimuth to target, in degrees
% output:
%   t1 = target latitude   ( + = N, - = S)
%   p1 = station longitude ( + = E, - = W)
%
      rad  = 57.29578;
      delr = del/rad;
      azr  = az/rad;
      
      colat = 90.0 - lat1;
      colon = lon1;
      if (colon < 0.0)
          colon = 360.0 + colon;
      end;

% convert to geocentric
      t0=geoc4(colat/rad);
      ctheta=sin(delr)*sin(t0)*cos(azr) + cos(t0)*cos(delr);
      t1=acos(ctheta);
      if t0 == 0.0  % source at north pole
        p1=az;      % set station longitude same as azimuth
      elseif t1 == 0.0   % station at north pole
        p1=0.0;          % set station longitude to 0
      else       % compute station longitude
        sphi=sin(delr)*sin(azr)/sin(t1);
        cphi=(cos(delr) - cos(t0)*ctheta)/(sin(t0)*sin(t1));
        p1=colon + atan2(sphi,cphi)*rad;
      end;
% convert colatitude to geograf. latitude
      t1 = 90.0-geog4(t1)*rad;
% assume p1 never > 720
      if p1 > 360.0
          p1 = p1 - 360.0;
      end;
% convert colongitude to longitude
      if p1 > 180.0
          p1 = p1 - 360.0;
      end;