function [xgrid,ygrid] = ll2grid(rlat,rlon)
%   ll2grid converts a latitude and longitude to the GLOBEC grid
%   using a conversion to Universal Transverse Mercator projection
%   
%   j. klinck, Feb 2001 (klinck@ccpo.odu.edu)

global SECS COSROT SINROT ONCE
global UTMXOFF UTMYOFF XOFF YOFF ROTANGLE IZONE ISPHERE
global AXIS BXIS AK0 RADSEC

%      set up transformation variables
%if exist('ONCE') == 0
if isempty(ONCE)
%         need to do only once
  ONCE=1;
  setuputm;
end
[m,n]=size(rlat);
if n==1
   rlat=rlat';
   rlon=rlon';
end

%         protect input and convert to seconds
%         change sign of lon for conversion
lat=rlat*SECS; lon=-rlon*SECS;

if IZONE == 5010
   i = find (lon < 0);
   lon(i)=lon(i)+1.296d6;
end
[xutm,yutm]=ll2utm(lat,lon);

xutm=xutm-UTMXOFF;
yutm=yutm-UTMYOFF;

xgrid=  COSROT*xutm + SINROT*yutm;
ygrid= -SINROT*xutm + COSROT*yutm;

%           convert to km and offset
%xgrid=round(XOFF + .001*xgrid);
%ygrid=round(YOFF + .001*ygrid);
xgrid=fix(XOFF + .001*xgrid);
ygrid=fix(YOFF + .001*ygrid);

function [xx,yy]=ll2utm(slat,slon);
%   convert from lat, lon to a universal transverse mercator
%   based on the parameters given in SETUPUTM.
%     universal transverse mercator conversion
%
%     rewritten 6-16-83 by j.f. waananen using formulation
%     by j.p. snyder in bulletin 1532, pages 68-69
%
%      rewrite utm.f for the case idir=0, which convert lat lon to utm
%       2/5/2001 j. klinck
%
global SECS COSROT SINROT ONCE
global UTMXOFF UTMYOFF XOFF YOFF ROTANGLE IZONE ISPHERE
global AXIS BXIS AK0 RADSEC

a = AXIS(ISPHERE);
b = BXIS(ISPHERE);
es= (a*a-b*b)/(a*a);

iizone = IZONE;

%     compute utm zone(izone) and central meridian in seconds for
%     geodetic to utm conversion where zone is not input.
izz=zeros(size(slon));
iutz=zeros(size(slon));
cm=zeros(size(slon));
iizone=zeros(size(slon));

if(IZONE ==  0) 
  izz=round(slon/2.16d4);
  iizone = 30 - izz;
  i=find(slon < 0);
     iizone(i)=31 - izz(i);
  i=find(abs(iizone) > 30);
    iutz(i) = abs(iizone(i))-30;
    cm(i)=((iutz(i)*6.0d0)-3.0d0)*(-3600.0d0);
  j=find(abs(iizone) <= 30);
    iutz(j) = 30-abs(iizone(j));
    cm(j)=((iutz(j)*6.0d0)+3.0d0)*3600.0d0;
else
  if(abs(IZONE) > 30)
      iutz = abs(IZONE)-30;
      cm=((iutz*6.0d0)-3.0d0)*(-3600.0d0);
  else
      iutz = 30-abs(IZONE);
      cm=((iutz*6.0d0)+3.0d0)*3600.0d0;
  end
  iizone=IZONE;
end

i=find(slat < 0);
  iizone(i)=-abs(IZONE);

i=find(iizone < 0);
   slat(i) = -abs(slat(i));

if(IZONE == 0) 
  iizone=IZONE;
end
i=find (slat <0 & iizone > 0);
  iizone(i)=-iizone(i);

%                 begin computation
phi = slat/RADSEC;
dlam = -(slon-cm)/RADSEC;
epri = es./(1.0d0-es);
en = a./sqrt(1.0d0-es*sin(phi).^2);
t = tan(phi).^2;
c = epri*cos(phi).^2;
aa = dlam.*cos(phi);
s2 = sin(2.0d0*phi);
s4 = sin(4.0d0*phi);
s6 = sin(6.0d0*phi);
f1 = (1.d0-(es/4.d0)-(3.d0*es.*es/6.4d1)-(5.d0*es.*es.*es/2.56d2));
f2 = ((3.d0*es/8.d0)+(3.d0*es.*es/3.2d1)+(4.5d1*es.*es.*es./1.024d3));
f3 = ((1.5d1*es.*es./2.56d2)+(4.5d1*es.*es.*es./1.024d3));
f4 = (3.5d1*es.*es.*es/3.072d3);
em = a.*(f1*phi - f2*s2 + f3*s4 - f4*s6);
xx = AK0*en.*(aa + (1.d0-t+c).*aa.^3/6.d0 + ...
      (5.d0-1.8d1*t+t.*t+7.2d1*c-5.8d1*epri).*aa.^5/1.2d2);
xx = xx + 5.0d5 ;
yy = AK0.*(em + en.*tan(phi).*((aa.*aa./2.d0) + ...
       (5.d0-t+9.d0*c+4.d0*c.*c).*aa.^4/2.4d1 + ...
       (6.1d1-5.8d1*t+t.*t+6.d2*c-3.3d2*epri).*aa.^6/7.2d2));
i=find(iizone < 0 | slat < 0.0d0);
   yy(i) = yy(i) + 1.0d7;

%  ---- latitude within 84 degrees
i=find(abs(slat) > 302400);
if (~isempty(i))
   disp(' *************   error exit from utm. ');
   disp('   latitude value is poleward of 84 degrees');
   disp('   utm is not valid.');
   disp('   calculation has continued but values ');
   disp('   may not be valid.');
   slat(i)/3600.d0
   slon(i)/3600.d0
end

%  ---- delta longitude within 0.16 radians of central meridian

%    print disabled to avoid extra output. Some points on the south west
%    side of the GLOBEC grid are just outside of the permissible area.

%i=find(abs(dlam) > 1.6d-1);
%if (~isempty(i))
%   disp(' *************   error exit from utm. ');
%   disp('  d lon not within 0.16 radians of central merid.');
%   disp('   calculation has continued but values ');
%   disp('   may not be valid.');
%   slat(i)/3600.d0
%   slon(i)/3600.d0
%end

%    clean up dummy variables
clear iiz iizone iutz cm phi dlam epri en t 
clear c aa s2 s4 s6 f1 f2 f3 f4 em

% end of ll2grid.m 
