%MEQMESH  Computational grid related quantities
% L = MEQMESH(L,R1,R2,NR,Z1,Z2,NZ,MBC) defines the fields of L:
%  Computational grid
%  .rx   = LINSPACE(R1,R2,NR)
%  .zx   = LINSPACE(Z1,Z2,NZ)
%  .nrx  = NR
%  .nzx  = NZ
%  .drx  = rx(2)-rx(1)
%  .dzx  = zx(2)-zx(1)
%  .idrx = 1/drx
%  .idzx = 1/dzx
%  .dsx  = drx*dzx
%  .idsx = 1/dsx,
%  .d2x  = drx^2+dzx^2
%  Inner mesh
%  .ry   = rx(2:end-1) 
%  .zy   = zx(2:end-1)
%  .nry  = NR-2
%  .nzy  = NZ-2
%  .iry  = 1/ry
%  Grid boundary (stored as [F(1:nz,1) F(1,2) F(nz,2) ... F(1:nz,nr)])
%  .rb,.zb boundary location
%  .lxb    for Fb=Fx(lxb)
%  Normal flux derivative on boundary,
%  dF(x)/dx=(4F(x+dx)-F(x+2dx)-3(F(x)=0))/2/dx for boundary condition
%  .kxc,.Tbc for Fb=Tbc*(Fx(kxc)*[4;-1])
%  L = MEQMESH(...,'P1',...}) optionally defines the fields
%  .rrx,.zzx = MESHGRID(RX,ZX)
%  .rry,.zzy = MESHGRID(RY,ZY)
%
% For details, see: [J-M. Moret et al. Fus.Eng.Des 2015], Sections 4.1,4.2
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function L = meqmesh(L,rx,zx,Mbc,varargin)
 
 %% Checks
 nr = numel(rx); nz = numel(zx);
 assert(~mod(log2(nz-1),1),'(numel(zx)-1) must be a power of 2')

 %% Mesh
 L.rx = rx;
 L.zx = zx;
 L.drx = rx(2)-rx(1); L.idrx = 1/L.drx; L.nrx = nr;
 L.dzx = zx(2)-zx(1); L.idzx = 1/L.dzx; L.nzx = nz;
 L.dsx = L.drx*L.dzx; L.idsx = 1/L.dsx; L.d2x = L.drx^2+L.dzx^2; L.nx = L.nrx*L.nzx;
 [rrx,zzx] = meshgrid(rx,zx);
 
 %% Boundary
 lxb = true(L.nzx,L.nrx); lxb(2:end-1,2:end-1) = false; nb = 2*(L.nzx+L.nrx)-4; % indices of x on boundary
 L.zb = zzx(lxb); L.rb = rrx(lxb);
 L.lxy = ~lxb; % points on x that are also on y
 
 %% Inner mesh
 L.ry  = L.rx(2:end-1); L.zy = L.zx(2:end-1);
 L.nry = L.nrx - 2; L.nzy = L.nzx - 2; L.ny = L.nry*L.nzy;
 L.iry = 1./L.ry;
 L.rry = rrx(2:end-1,2:end-1); L.zzy = zzx(2:end-1,2:end-1);

 %% For normal flux derivative on boundary (no corner)
 lxc = lxb; lxc([1 end],[1 end]) = false;
 rc = rrx(lxc);
 cc = [repmat(L.dzx/L.drx,L.nzx-2,1);repmat(L.drx/L.dzx,2*L.nrx-4,1);repmat(L.dzx/L.drx,L.nzx-2,1)] ./ rc / (4*pi*mu0);
 drc = [ zeros(      L.nzx-2,1);repmat(L.drx,2*L.nrx-4,1); zeros(      L.nzx-2,1)];
 dzc = [repmat(L.dzx,L.nzx-2,1); zeros(      2*L.nrx-4,1);repmat(L.dzx,L.nzx-2,1)];
 Lcc = mu0*rc.*((1+(1/32)*(dzc./rc).^2+(1/96)*(drc./rc).^2).*log(8*rc./(drc+dzc))+((1/128)*(dzc./rc).^2-0.5)); % FH used mu0*r*(ln(16r/ds)-1+(ds/r)^2*(ln(16r/ds)-2/3))
 Mbc(sub2ind([nb,nb-4],find(lxc(lxb))',1:nb-4)) = Lcc;
 L.Tbc = Mbc.*cc';

 for k = 1:length(varargin)
  switch varargin{k}
   case 'rrx', L.rrx = rrx;
   case 'zzx', L.zzx = zzx;
  end
 end
 
 %% for field computation using finite elements
 if L.P.ifield
    L.i2pirxdzx = 1./(4*pi*L.dzx*L.rx');
    L.i2pirxdrx = 1./(4*pi*L.drx*L.rx');
 end
 
 %% Optional additional grid for post-processing
 if L.P.izgrid
   L.nrz = numel(L.G.rz);  L.nzz = numel(L.G.zz);
   L.drz = L.G.rz(2)-L.G.rz(1); L.idrz = 1/L.drz; L.nrz = numel(L.G.rz);
   L.dzz = L.G.zz(2)-L.G.zz(1); L.idzz = 1/L.dzz; L.nzz = numel(L.G.zz);
   L.rz = L.G.rz; L.zz = L.G.zz;
   [L.rrz,L.zzz] = meshgrid(L.rz,L.zz);
   if L.P.ifield
     L.i2pirzdzz = 1./(4*pi*L.dzz*L.rz');
     L.i2pirzdrz = 1./(4*pi*L.drz*L.rz');
   end
 end
 
end