function [Lonr,Latr,pm,pn,ang] = ...
           oct_easy_grid(nx,ny,dl,size_x,size_y,tra_lon,tra_lat,rot);
%
%==================================================================== 
%   Easy grid make an rectangular, orthogonal grid with minimal 
%   gridsize variation. It uses a Mercator projection around the 
%   equator and then rotates the sphere around its 3 axis to position 
%   the grid where it is desired.
%   
%   Inputs:
%   nx:      Number of grid point in the x direction
%   ny:      Number of grid point in the y direction
%   size_x:  Domain size in x-direction
%   size_y:  Domain size in y-direction
%   tra_lon: Desired longitude of grid center
%   tra_lat: Desired latitude of grid center
%   rot:     Rotation of grid direction (0: x direction is west-east)
%
%  (c) 2008, Jeroen Molemaker, UCLA
%      2019, P. Marchesiello, modif. for croco_tools
%==================================================================== 
%
R_earth=6367442.76;   % Earth radius
deg2rad=pi/180;
rad2deg=180/pi;

if rot==0

%
% Get the longitude and latitude as in croco_tools
%
dlon=size_x/(R_earth*cos(tra_lat*deg2rad))*rad2deg;
dlat=size_y/R_earth*rad2deg;
lonmin=tra_lon-0.5*dlon;
lonmax=tra_lon+0.5*dlon;
latmin=tra_lat-0.5*dlat;
latmax=tra_lat+0.5*dlat;
lonr=(lonmin:dl:lonmax);
i=1;
latr(i)=latmin;
while latr(i)<=latmax
  i=i+1;
  latr(i)=latr(i-1)+dl*cos(latr(i-1)*pi/180);
end
[Lonr,Latr]=meshgrid(lonr,latr);
[Lonu,Lonv,Lonp]=oct_rho2uvp(Lonr); 
[Latu,Latv,Latp]=oct_rho2uvp(Latr);

else

%
% Perform Mercator projection around the equator
% before rotation and translation to center lat,lon
%
if (size_y>size_x) 
  length = size_y; nl = ny;
  width  = size_x; nw = nx;
else
  length = size_x; nl = nx;
  width  = size_y; nw = ny;
end
   
dlon = length/R_earth;
lonr = dlon*[-0.5:1:nl+0.5]/nl - dlon/2;  

dlat = width/R_earth;
y1 = log(tan(pi/4-dlat/4));
y2 = log(tan(pi/4+dlat/4));
y = (y2-y1)*[-0.5:1:nw+0.5]/nw + y1;  
latr = 2*atan(exp(y)) - pi/2; 
dlat_cen = 0.5*(latr(round(nw/2)+1)-latr(round(nw/2)-1));
dlon_cen = dlon/nl;
mul = dlat_cen/dlon_cen*length/width*nw/nl;
latr = latr/mul;

[Lonr,Latr] = meshgrid(lonr,latr);
Lonu = 0.5*(Lonr(:,1:end-1)+Lonr(:,2:end));
Latu = 0.5*(Latr(:,1:end-1)+Latr(:,2:end));
Lonv = 0.5*(Lonr(1:end-1,:)+Lonr(2:end,:));
Latv = 0.5*(Latr(1:end-1,:)+Latr(2:end,:));

if (size_y>size_x) 
 [Lonr,Latr] = rot_sphere(Lonr,Latr,90);
 [Lonu,Latu] = rot_sphere(Lonu,Latu,90);
 [Lonv,Latv] = rot_sphere(Lonv,Latv,90);
 Lonr = flipdim(Lonr,1)';
 Latr = flipdim(Latr,1)';
 Lonu_tmp = flipdim(Lonv,1)';
 Latu_tmp = flipdim(Latv,1)';
 Lonv = flipdim(Lonu,1)';
 Latv = flipdim(Latu,1)';
 Lonu = Lonu_tmp;
 Latu = Latu_tmp;
end

[Lonr,Latr] = rot_sphere(Lonr,Latr,rot);
[Lonu,Latu] = rot_sphere(Lonu,Latu,rot);
[Lonv,Latv] = rot_sphere(Lonv,Latv,rot);
  
[Lonr,Latr] = tra_sphere(Lonr,Latr,tra_lat);
[Lonu,Latu] = tra_sphere(Lonu,Latu,tra_lat);
[Lonv,Latv] = tra_sphere(Lonv,Latv,tra_lat);

Lonr = Lonr*rad2deg + tra_lon;
Lonu = Lonu*rad2deg + tra_lon;
Lonv = Lonv*rad2deg + tra_lon;
Latr = Latr*rad2deg;
Latu = Latu*rad2deg;
Latv = Latv*rad2deg;

Lonr(Lonr<-180) = Lonr(Lonr<-180) + 360;
Lonu(Lonu<-180) = Lonu(Lonu<-180) + 360;
Lonv(Lonv<-180) = Lonv(Lonv<-180) + 360;

end % rot==0

%
% Compute pn and pm (as in croco_tools)
%
[Mp,L]=size(Latu);
[M,Lp]=size(Latv);
Lm=L-1;
Mm=M-1;
dx=zeros(Mp,Lp);
dy=zeros(Mp,Lp);
dx(:,2:L)=oct_spheric_dist(Latu(:,1:Lm),Latu(:,2:L),...
                       Lonu(:,1:Lm),Lonu(:,2:L));
dx(:,1)=dx(:,2);
dx(:,Lp)=dx(:,L);

dy(2:M,:)=oct_spheric_dist(Latv(1:Mm,:),Latv(2:M,:),...
                       Lonv(1:Mm,:),Lonv(2:M,:));
dy(1,:)=dy(2,:);
dy(Mp,:)=dy(M,:);

pm=1./dx;
pn=1./dy;
%
% Compute angles of local grid positive x-axis relative to east
% (as in croco_tools)
%
ang=oct_get_angle(Latu,Lonu)*rad2deg;

return

%
%==================================================================== 
%  Rotate sphere around its y-axis
%  Part of Easy Grid
%  (c) 2008, Jeroen Molemaker, UCLA
%==================================================================== 
%
function [lon2,lat2] = rot_sphere(lon1,lat1,rot)


[n,m] = size(lon1);
rot = rot*pi/180;
eps = 1.e-20;
%
% Translate into x,y,z
% conventions:  (lon,lat) = (0,0)  corresponds to (x,y,z) = ( 0,-r, 0)
%               (lon,lat) = (0,90) corresponds to (x,y,z) = ( 0, 0, r)
%
x1 = sin(lon1).*cos(lat1);
y1 = cos(lon1).*cos(lat1);
z1 = sin(lat1);
%
% We will rotate these points around the small circle defined by 
% the intersection of the sphere and the plane that
% is orthogonal to the line through (lon,lat) (0,0) and (180,0)
%
% The rotation is in that plane around its intersection with
% aforementioned line.
%
% Since the plane is orthogonal to the y-axis (in my definition at least),
% Rotations in the plane of the small circle maintain constant y and are around
% (x,y,z) = (0,y1,0)
%
rp1 = sqrt(x1.^2+z1.^2);

ap1 = pi/2*ones(n,m);
ap1(abs(x1)>eps) = atan(abs(z1(abs(x1)>eps)./x1(abs(x1)>eps)));
ap1(x1<0) = pi - ap1(x1<0);
ap1(z1<0) = -ap1(z1<0);

ap2 = ap1 + rot;
x2 = rp1.*cos(ap2);
y2 = y1;
z2 = rp1.*sin(ap2);

%
% transform back from (x,y,z) to (lat,lon)
%
lon2 = pi/2*ones(n,m);
lon2(abs(y2)>eps) = atan(abs(x2(abs(y2)>eps)./y2(abs(y2)>eps)));
lon2(y2<0) = pi - lon2(y2<0);
lon2(x2<0) = -lon2(x2<0);

pr2 = sqrt(x2.^2+y2.^2);
lat2 = pi/2*ones(n,m);
lat2(abs(pr2)>eps) = atan(abs(z2(abs(pr2)>eps)./pr2(abs(pr2)>eps)));
lat2(z2<0) = -lat2(z2<0);

return

%
%==================================================================== 
%  Translate rotated grid
%  Part of oct_easy grid
%  (c) 2008, Jeroen Molemaker, UCLA
%==================================================================== 
%
function [lon2,lat2] = tra_sphere(lon1,lat1,tra)


[n,m] = size(lon1);
tra = tra*pi/180;  %% translation in latitude direction
%
% translate into x,y,z
% conventions:  (lon,lat) = (0,0)  corresponds to (x,y,z) = ( 0,-r, 0)
%               (lon,lat) = (0,90) corresponds to (x,y,z) = ( 0, 0, r)
%
x1 = sin(lon1).*cos(lat1);
y1 = cos(lon1).*cos(lat1);
z1 = sin(lat1);

rp1 = sqrt(y1.^2+z1.^2);

ap1 = pi/2*ones(n,m);
ap1(abs(y1)>1e-5) = atan(abs(z1(abs(y1)>1e-5)./y1(abs(y1)>1e-5)));
ap1(y1<0) = pi - ap1(y1<0);
ap1(z1<0) = -ap1(z1<0);

ap2 = ap1 + tra;
x2 = x1;
y2 = rp1.*cos(ap2);
z2 = rp1.*sin(ap2);
%
% transform back from (x,y,z) to (lat,lon)
%
lon2 = pi/2*ones(n,m);
lon2(abs(y2)>1e-5) = atan(abs(x2(abs(y2)>1e-5)./y2(abs(y2)>1e-5)));
lon2(y2<0) = pi - lon2(y2<0);
lon2(x2<0) = -lon2(x2<0);

pr2 = sqrt(x2.^2+y2.^2);
lat2 = pi/2*ones(n,m);
lat2(abs(pr2)>1e-5) = atan(abs(z2(abs(pr2)>1e-5)./pr2(abs(pr2)>1e-5)));
lat2(z2<0) = -lat2(z2<0);

return