function [phi_eq,lambda_eq] = lltoeq(phi,lambda,phi_pole,lambda_pole)

% Calculates latitude and longitude on equatorial latitude-longitude (eq) grid used
% in regional models from a standard grid.
%
% CEB Sept 2010
% Converted from John's lltoeq.pro
%  written by a. dickinson 11/8/89
%
%  reviewed by r. swinbank 19/01/90
%
%  purpose:
%            calculates latitude and longitude on equatorial
%            latitude-longitude (eq) grid used in regional
%            models from input arrays of latitude and
%            longitude on standard grid. both input and output
%            latitudes and longitudes are in degrees.
%
%  documentation:
%            the transformation formulae are described in
%            unified model on-line documentation paper s1.
%  24/09/2002 - Name of routine sgn changed to apply_sgn
%
%  ------------------------------------------------------------------
%  Arguments:--------------------------------------------------------
%
% phi               !in latitude
% lambda            !in longitude (0 =< lon < 360)
% lambda_eq         !out longitude in equatorial lat-lon coords
% phi_eq            !out latitude in equatorial lat-lon coords
% phi_pole          !in  latitude of equatorial lat-lon pole
% lambda_pole       !in  longitude of equatorial lat-lon pole

small=1.0e-6;

pi=3.14159265358979323846;    % pi
pi_over_180 =pi/180.0;        % conversion factor degrees to radians
recip_pi_over_180 = 180.0/pi; % conversion factor radians to degrees
% latitude of zeroth meridian
lambda_zero=lambda_pole+180;
% sine and cosine of latitude of eq pole
sin_phi_pole=sin(pi_over_180*phi_pole);
cos_phi_pole=cos(pi_over_180*phi_pole);

% transform from standard to equatorial latitude-longitude
% scale longitude to range -180 to +180 degs   
a_lambda=lambda-lambda_zero;
i=find(a_lambda > 180.0);
if length(i)>0
  a_lambda(i)=a_lambda(i)-360.0;      
end    
i=find(a_lambda < -180.0);
if length(i)>0
  a_lambda(i)=a_lambda(i)+360.0;      
end 

% convert latitude & longitude to radians
a_lambda=pi_over_180.*a_lambda;
a_phi=pi_over_180.*phi;

% compute eq latitude using equation (4.4)
arg=-1.*cos_phi_pole.*cos(a_lambda).*cos(a_phi) + sin(a_phi).*sin_phi_pole;

i=find(arg>1.0);
if length(i)>0
  arg(i)=1.0;  
end  
i=find(arg<-1.0);
if length(i)>0
  arg(i)=-1.0;  
end 
e_phi=asin(arg);
phi_eq=recip_pi_over_180.*e_phi;

% compute eq longitude using equation (4.6)
term1 =cos(a_phi).*cos(a_lambda).*sin_phi_pole + sin(a_phi).*cos_phi_pole;
term2 =cos(e_phi);
i1=find(term2 < small);
i2=find(term2 > small);
e_lambda=lambda;
if length(i1)>0
  e_lambda(i1)=0.0;  
end    
if length(i2)>0
  arg=term1(i2)./term2(i2);
  i=find(arg>1,1);
  %i=find(arg>1,0);
  if length(i)>1.0
    arg(i)=1.0;  
  end 
  i=find(arg<-1,1);
  %i=find(arg<-1,0);
  if length(i)>1.0
    arg(i)=-1.0;  
  end  
  e_lambda(i2)=recip_pi_over_180.*acos(arg);
  e_lambda(i2)=apply_sgn(e_lambda(i2),a_lambda(i2));  
end


% scale longitude to range 0 to 360 degs
i=find(e_lambda > 360.0);
if size(i)>0
  e_lambda(i)=e_lambda(i)-360.0;  
end    
i=find(e_lambda < 0.0);
if size(i)>0
  e_lambda(i)=e_lambda(i)+360.0;  
end   
lambda_eq=e_lambda;