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

% Calculates latitude and longitude on standard grid from input arrays of 
% latitude and longitude on equatorial latitude-longitude (eq) grid used
% in regional models. Both input and output latitudes and longitudes are in degrees.
%
% CEB Jun 2010
% Converted from John's eqtoll.pro
% Written by Pete Clark, translating A. Dickinson FORTRAN to WAVE
%
% Documentation: The transformation formulae are described in
% unified model on-line documentation paper S1.
%
% Inputs:
% phi_eq - latitude in equatorial lat-lon coords  model
% lambda_eq - longitude in equatorial lat-lon coords 
% phi_pole - latitude of equatorial lat-lon pole
% lambda_pole longitude of equatorial lat-lon pole
%
% Outputs:
% phi - latitude, normal
% lambda - longitude (0 =< lon < 360)


small=1.0e-7; % ??

pi_over_180 = pi/180; % conversion factor degrees to radians
recip_pi_over_180 = 1./(pi/180); % conversion factor radians to degrees

% initialise local constants
% latitude of zeroth meridian
if (lambda_pole > 90) | (lambda_pole < -90)
  lambda_zero=lambda_pole+180;   
else    
  lambda_zero=lambda_pole;
end

% 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 equatorial to standard latitude-longitude
% scale eq longitude to range -180 to +180 degs

e_lambda=lambda_eq;
i=find(e_lambda > 180.0);
e_lambda(i)=e_lambda(i)-360.0;
i=find(e_lambda < -180.0);
e_lambda(i)=e_lambda(i)+360.0;

% convert eq latitude & longitude to radians

e_lambda=pi_over_180.*e_lambda;
e_phi=pi_over_180.*phi_eq;

% compute latitude using equation (4.7)
arg=cos_phi_pole.*cos(e_lambda).*cos(e_phi) + sin(e_phi).*sin_phi_pole;
i=find(arg > 1.0);
arg(i)=1.0;
i=find(arg < -1.0);
arg(i)=-1.0;
a_phi=asin(arg);
phi=recip_pi_over_180.*a_phi;

% compute longitude using equation (4.8)

term1 =(cos(e_phi).*cos(e_lambda).*sin_phi_pole - sin(e_phi).*cos_phi_pole);
term2=cos(a_phi);

i1=find(term2<small);
i2=find(term2>small);

a_lambda=term1; % dummy to set up space
a_lambda(i1)=0.0;

arg=term1/term2;
i=find(arg>(1.0-small));
arg(i) = 1.0 -small;

a_lambda(i2)=recip_pi_over_180.*acos(arg(i2));
a_lambda(i)=0;
a_lambda(i2) = apply_sgn(a_lambda(i2),e_lambda(i2));
a_lambda(i2) = a_lambda(i2)+lambda_zero(i2);

% scale longitude to range 0 to 360 degs

i1=find(a_lambda > 360.0);
a_lambda(i1)=a_lambda(i1)-360.0;
i1=find(a_lambda < 0.0);
a_lambda(i1)=a_lambda(i1)+360.0;
lambda=a_lambda;