function [x,y] = ll2xy (lon,lat,coord,nx,ny,gridx,gridy,rlngd,rlatd, ...
                        delx,dely,thetad)
%
%=======================================================================
%                                                                    ===
% function [x,y] = ll2xy (lon,lat,coord,nx,ny,gridx,gridy,rlngd, ... ===
%                                        rlatd,delx,dely,thetad);    ===
%                                                                    ===
% This routine computes the grid coordinates from the given          ===
% longitude and latitude.                                            ===
%                                                                    ===
% ------                                                             ===
% Input:                                                             ===
% ------                                                             ===
%                                                                    ===
%     LON, LAT  the longitude and latitude of the point    (deg)     ===
%     COORD     flag for coordinate type                             ===
%                  [0]  Cartesian (tangent plane) grid.              ===
%                  [1]  unrotated spherical grid.                    ===
%                  [2]  rotated spherical grid.                      ===
%     DELX,                                                          ===
%     DELY,     Offset between transform & grid centers.   (m | deg) ===
%     GRIDX     grid spacing in x-direction                (m | deg) ===
%     GRIDY     grid spacing in y-direction                (m | deg) ===
%     NX, NY    number of grid points in x & y directions            ===
%     RLNGD     longitude at center of transformation      (deg)     ===
%     RLATD     latitude at center of transformation       (deg)     ===
%     THETAD    angle grid is rotates with respect to east (deg)     ===
%                                                                    ===
% -------                                                            ===
% Output:                                                            ===
% -------                                                            ===
%                                                                    ===
%     X, Y      the corresponding grid coordinates                   ===
%                                                                    ===
%=======================================================================

%=======================================================================
%                                                                    ===
% Calls:  rotparm                                                    ===
%                                                                    ===
%=======================================================================

deg2rad = pi/180;
xc      = (nx+1)/2 - (delx/gridx);
yc      = (ny+1)/2 - (dely/gridy);

%-----------------------------------------------------------------------
%  Cartesian (tangent plane) system.
%-----------------------------------------------------------------------

if (coord==0),
   re_d2r = 6371315.0*deg2rad;
   f1x    = re_d2r/gridx;
   f1y    = re_d2r/gridy;
   f2x    = f1x*cos(rlatd*deg2rad);
   f2y    = f1y*cos(rlatd*deg2rad);
   sterm  = sin(thetad*deg2rad);
   cterm  = cos(thetad*deg2rad);
   dlat   = lat-rlatd;
   dlon   = lon-rlngd;
   x      = xc+f1x.*dlat.*sterm+f2x.*dlon.*cterm;
   y      = yc+f1y.*dlat.*cterm-f2y.*dlon.*sterm;

%-----------------------------------------------------------------------
%  Unrotated Spherical system.
%-----------------------------------------------------------------------

  elseif (coord==1),
   x = xc+lon./gridx;
   y = yc+lat./gridy;

%-----------------------------------------------------------------------
%  Rotated spherical system.
%-----------------------------------------------------------------------

  else,
   rad2deg = 180/pi;
   [plon,plat,rnlon] = rotparm (rlngd*deg2rad,rlatd*deg2rad,thetad*deg2rad);
   cosp = cos(plat);
   sinp = sin(plat);
   cosl = cos(lat.*deg2rad);
   sinl = sin(lat.*deg2rad);
   rlat = asin(sinl.*sinp+cosl.*cosp.*cos(lon.*deg2rad-plon));
   sin1 = cosl.*sin(lon.*deg2rad-plon);
   cos1 = sinl.*cosp-sinp.*cosl.*cos(lon.*deg2rad-plon);
   rlon = rnlon-atan2(sin1,cos1);
   ind  = find(abs(rlon)>pi);
   if (~isempty(ind)),
      rlon(ind) = sign(rlon(ind)).*(rem(abs(rlon(ind)),2*pi) - ...
                             2.*pi.*fix(rem(abs(rlon(ind)),2*pi)./pi));
   end;
   rlon = rlon.*rad2deg;
   rlat = rlat.*rad2deg;
   x    = xc+rlon./gridx;
   y    = yc+rlat./gridy;

end;
