function [plon,plat,rnlon] = rotparm (clon,clat,rotang)
%
%=======================================================================
%                                                                    ===
% function [plon,plat,rnlon] = rotparm (clon,clat,rotang);           ===
%                                                                    ===
% This routine computes the parameters for a spherical rotation      ===
% such that the new equator passes through the center point,         ===
% making the prescribed angle with respect to East.                  ===
%                                                                    ===
% On Input:                                                          ===
%                                                                    ===
%     CLON, CLAT   The longitude and latitude of point center of the ===
%                  grid.  (radians)                                  ===
%     ROTANG       The angle of rotation with respect to East.       ===
%                  (radians)                                         ===
%                                                                    ===
% On Output:                                                         ===
%                                                                    ===
%     PLON         The longitude of the rotated pole in the original ===
%                  system.  (radians)                                ===
%     PLAT         The latitude of the rotated pole in the original  ===
%                  system.  (radians)                                ===
%     RNLON        The longitude of the original pole in the rotated ===
%                  system.  (radians)                                ===
%                                                                    ===
%=======================================================================

%=======================================================================
%                                                                    ===
% Calls:  get_ewpt                                                   ===
%                                                                    ===
%=======================================================================

qpi = pi/4;

%-----------------------------------------------------------------------
% Compute the "eastern" and "western" points along the rotated equator.
%-----------------------------------------------------------------------

[elon,elat,wlon,wlat] = get_ewpt (clon,clat,rotang);

%-----------------------------------------------------------------------
% Compute the trigonometric functions of east and west latitudes and
% of the difference of their longitudes.
%-----------------------------------------------------------------------

sine = sin(elat);
sinw = sin(wlat);
cose = cos(elat);
cosw = cos(wlat);

cosemw = cos(elon-wlon);
sinemw = sin(elon-wlon);

%-----------------------------------------------------------------------
% Compute the trigonometric functions of the interior angles.
%-----------------------------------------------------------------------

sin_pen = sinw*cose - sine*cosw*cosemw;
cos_pen = cosw*sinemw;
sin_pwn = sine*cosw - sinw*cose*cosemw;
cos_pwn = cose*sinemw;

%-----------------------------------------------------------------------
% Compute latitude of rotated pole.
%-----------------------------------------------------------------------

plat = asin( cose*cosw*sinemw );

%-----------------------------------------------------------------------
% Compute longitude of rotated pole.  Keep all longitudes in the
% range (-pi, pi].
%-----------------------------------------------------------------------

numer = sin_pen;
denom = - cos_pen*sine;
lon1 = elon + atan2(numer,denom);
lon1=sign(lon1)*(rem(abs(lon1),2*pi) - 2*pi*fix(rem(abs(lon1),2*pi)/pi));

numer = sin_pwn;
denom = - cos_pwn*sinw;
lon2 = wlon - atan2(numer,denom);
lon2=sign(lon2)*(rem(abs(lon2),2*pi) - 2*pi*fix(rem(abs(lon2),2*pi)/pi));

plon = ( lon1 + lon2 )/2;
plon=sign(plon)*(rem(abs(plon),2*pi) - 2*pi*fix(rem(abs(plon),2*pi)/pi));

%-----------------------------------------------------------------------
% Compute longitude of north pole in rotated coordinates.  Keep all
% longitudes in the range (-pi, pi].
%-----------------------------------------------------------------------

numer = sin_pen*cose;
denom = sine;
lon1 = qpi - atan2(numer,denom);
lon1=sign(lon1)*(rem(abs(lon1),2*pi) - 2*pi*fix(rem(abs(lon1),2*pi)/pi));

numer = sin_pwn*cosw;
denom = sinw;
lon2 = atan2(numer,denom) - qpi;
lon2=sign(lon2)*(rem(abs(lon2),2*pi) - 2*pi*fix(rem(abs(lon2),2*pi)/pi));

rnlon = ( lon1 + lon2 )/2;
rnlon=sign(rnlon)*(rem(abs(rnlon),2*pi)-2*pi*fix(rem(abs(rnlon),2*pi)/pi));
