function [elon,elat,wlon,wlat] = get_ewpt (clon,clat,rotang)
%
%=======================================================================
%                                                                    ===
% function [elon,elat,wlon,wlat] = get_ewpt (clon,clat,rotang);      ===
%                                                                    ===
% This routine computes two points to determine a spherical          ===
% rotation.  In rotated system, the points would lie on the          ===
% "equator" at "longitudes" 45E and 45W.                             ===
%                                                                    ===
% On Input:                                                          ===
%                                                                    ===
%     CLAT, CLON   The center of the grid.  (radians)                ===
%     ROTANG       The angle the new "equator" make with geographic  ===
%                  east at the center of the grid.  (radians)        ===
%                                                                    ===
% On Output:                                                         ===
%                                                                    ===
%     ELAT, ELON   The geographic coordinates of the "east point".   ===
%                  New latitude=0, New longitude=45.  (radians)      ===
%     WLAT, WLON   The geographic coordinates of the "west point".   ===
%                  New latitude=0, New longitude=-45. (radians)      ===
%                                                                    ===
%=======================================================================

%=======================================================================
%                                                                    ===
% Calls:  none                                                       ===
%                                                                    ===
%=======================================================================

rrt2 = 1/sqrt(2);
hpi = pi/2;

%-----------------------------------------------------------------------
% Find "eastern" point.
%-----------------------------------------------------------------------

colat = acos( rrt2*(sin(clat) + cos(clat)*sin(rotang)) );
elat = hpi - colat;

if ( (colat~=0) & (colat~=pi) ),

   num = rrt2*(cos(clat) - sin(clat)*sin(rotang));
   if (abs(rotang)<=pi/2),
      elon = clon + acos( num/sin(colat) );
     else,
      elon = clon - acos( num/sin(colat) );
   end;

  else

   elon = clon;

end;

%-----------------------------------------------------------------------
% Find "western" point.
%-----------------------------------------------------------------------

colat = acos( rrt2*(sin(clat) - cos(clat)*sin(rotang)) );
wlat = hpi - colat;

if ( (colat~=0) & (colat~=pi) ),

   num = rrt2*(cos(clat) + sin(clat)*sin(rotang));
   if (abs(rotang)<=pi/2),
      wlon = clon - acos( num/sin(colat) );
     else,
      wlon = clon + acos( num/sin(colat) );
   end;

  else,

   wlon = clon;

end;

%-----------------------------------------------------------------------
% Make sure latitudes are in the range [-pi/2, pi/2].
%    Step 1:  ensure that latitudes are in the range [-pi, pi].
%    Step 2:  further correct to [-pi/2, pi/2] (this involves
%             the longitudes).
% Make sure longitudes are in the range [-pi, pi].
%-----------------------------------------------------------------------

elat=sign(elat)*(rem(abs(elat),2*pi) - 2*pi*fix(rem(abs(elat),2*pi)/pi));
wlat=sign(wlat)*(rem(abs(wlat),2*pi) - 2*pi*fix(rem(abs(wlat),2*pi)/pi));

if (elat>pi/2),
   elat = pi - elat;
   elon = elon + pi;
  elseif (elat<-pi/2),
   elat = mpi - elat;
   elon = elon + pi;
end;

if (wlat>pi/2),
   wlat = pi - wlat;
   wlon = wlon + pi;
  elseif (wlat<-pi/2),
   wlat = mpi - wlat;
   wlon = wlon + pi;
end;

elon=sign(elon)*(rem(abs(elon),2*pi) - 2*pi*fix(rem(abs(elon),2*pi)/pi));
wlon=sign(wlon)*(rem(abs(wlon),2*pi) - 2*pi*fix(rem(abs(wlon),2*pi)/pi));
