%MEQG  MEQ geometry and electromagnetic parameters
% G = MEQG(G,P,varargin) defines missing fields in the structure G using parameters
% in structure P.
% 
% Inputs: 
% G structure with some of the following fields:
% .rx,zx            Computational grid
% .rl,zl            Limiter contour
% .rw,zw,ww,hw,Twa  Active coil filament position, connection matrix:    Iw=Twa*Ia
% .rv,zv,wv,hv,Tvu  Vessel filament position, general vessel description matrix: Iv=Tvu*Iu
% .rf,zf            Flux loop position
% .rm,zm,am         Magnetic probe position and angle [ctr-clockwise in R,Z plane]
% .rW,zW,oW,aW      Wall gap origin, orientation angle, length
%  Notes:
%   * instead of .wv,.hv, also .dv can be specfied for circular filaments
%   * .Tivs can be passed as the matrix such that Is = Tivs*Iv for vessel segments
%
% P structure with meq parameters
%
% Optional additional Green's functions may be requested via varargin
%   e.g. G = MEQG(G,P,'Mxv')
%
% Output: G structure with added Green's functions such as
%         .Mby,Mbc,Mxa,Maa,Mxu,Mfx,Mfa,Mfu,Bmx,Bma,Bmu,dzMfx,dzBmx,Brxa,Brxu
%         depending on selection
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function G = meqg(G,P,varargin)

assert(isstruct(G)&&isstruct(P),'first two arguments must be G,P structures');

% check genlib presence
if isempty(which('greenem'))
  error('greenem.m not found, please add genlib to the matlab path')
end

% extra arguments
xargs = varargin;
if P.gsxe == 2
  xargs = [xargs,{'Mboa','Mbou','Mxia','Mxiu'}]; % for GS solution method with exact external currents
end 
if P.izgrid % extra fields for additional z grid
  xargs = [xargs,{'Mza','Mzu','Mzy'}];
end
% default fields to get
defaultfields = {'Mbc','Mxa','Mxu','Mfx','Mfu','Mfa','Bmx','Bmu','Bma'};

% final fields, sort for efficient retrieval
gfields = sortxargs(unique([defaultfields,xargs]));

% Windings & Vessel consistency: either rectangular (w,h) or circular (d)
for obj={'v','w'}
  dd = ['d',obj{:}]; ww=['w',obj{:}]; hh=['h',obj{:}];
  if ~isfield(G,dd), G.(dd)=nan(size(G.(ww))); end   % default d if not specified
  assert(~any(any(isnan([G.(ww),G.(hh)]))) || ...  % all w,h specified, OR
    (all(xor(all(isnan([G.(ww),G.(hh)]),2),isnan(G.(dd))))),... d specified for nan [w,h] entries
    'all windings must specify either [%s,%s] or %s, and the other must be nan',ww,hh,dd);
end

%% Mesh
if all(~isfield(G,{'rx','zx'}))
  [G.rx,G.zx] = meqgx(P,G); % subfunction
elseif any(~isfield(G,{'rx','zx'}))
  error('G provided to meqg contains only one of rx and zx fields.');
end
[G.rrx,G.zzx] = meshgrid(G.rx,G.zx);

if P.izgrid % optional extended grid
  if ~all(isfield(P,{'rz','zz'}))
    % default grid wrapping all conducting structures, with margin
    G.rz = linspace(min([G.rw;G.rv])*0.9,max([G.rw;G.rv])*1.1,31)';
    G.zz = linspace(min([G.zw;G.zv])*1.1,max([G.zw;G.zv])*1.1,41)';
  else
    G.rz = reshape(P.rz,[],1); G.zz = reshape(P.zz,[],1); % reshape to column vector
    assert(G.rz(1)>0,'all R grid values must be greater than 0')
  end
  [G.rrz,G.zzz] = meshgrid(G.rz,G.zz);
end

% Optional load of cached .mat with geometry information
if ~strcmp(P.tok,'file') && ~isempty(P.gfile)
  G = meqgload(G,P.gfile,P.debug);
end

%% Limiter
if ~isfield(G,'rl') % center of edge grid cells
  rxh = 0.5*(G.rx(1:end-1)+G.rx(2:end));
  zxh = 0.5*(G.zx(1:end-1)+G.zx(2:end));
  nrx = numel(G.rx);
  nzx = numel(G.zx);
  G.rl = [rxh;rxh(end)*ones(nzx-3,1);rxh(end:-1:1);rxh(1)*ones(nzx-3,1)];
  G.zl = [zxh(1)*ones(nrx-2,1);zxh;zxh(end)*ones(nrx-3,1);zxh(end:-1:2)];
end

%% Points in between the limiter points
G.rlh = 0.5*(G.rl + G.rl([2:end,1]));
G.zlh = 0.5*(G.zl + G.zl([2:end,1]));

%% Coils


%% Generalized vessel description
if ~all(isfield(G,{'Tvu','dimu'}))
  
  if ~isequal(P.selu,'n')
    assert(all(isfield(G,{'rv','zv'})),'need to specify vessel description rv,zv')
  elseif ~all(isfield(G,{'rv','zv'}))
    G.rv = []; G.zv = []; % default empty vessel
  end
  G.nv = numel(G.rv);

  switch P.selu
    case {'s','sx'}
      if ~isfield(G,'Tvu')
        assert(isfield(G,'Tvs'),'must pass Tvu or Tvs for custom segment description'); 
        G.Tvu = G.Tvs; 
      end 
      % segment description
      if     isfield(G,'Rs'); G.Ru = G.Rs;
      elseif isfield(G,'Rv'), G.Ru = diag(G.Tvu'*diag(G.Rv)*G.Tvu);
      else,   G.Ru = [];
      end
      G.dimu = G.dims;
    case {'e','ex'}
      % calculate eigenmode description from v
      nu = P.nu;
      if  isempty(nu)
        % Use all eigenvalues by default
        nu = G.nv;
      else
        assert(nu>0 && nu<=G.nv,...
          'P.nu must be nonzero, numeric and smaller than G.nv = %d',G.nv)
      end
      % Equations from JMM TCV toolbox documentation
      G.Mvv = getg(G,'Mvv');
      [Tve] = vveig(G.Mvv,G.Rv); % use vveig from genlib
      G.Tvu = Tve(:,1:nu);
      if isfield(G,'Rv'),  G.Ru = diag(G.Tvu'*diag(G.Rv)*G.Tvu);
      else, G.Ru = [];
      end
      G.dimu = cellstr(num2str((1:nu)','e%03d'));
    case {'v','vx'}
      assert(isempty(P.nu)||P.nu == G.nv,'Must set P.nu=%d (or empty) for this vessel description %s.',G.nv,P.selu)
      G.Tvu = eye(G.nv);
      if isfield(G,'dimv')
        G.dimu = G.dimv;
      else
        G.dimu = cellstr(num2str((1:G.nv)','v%03d'));
      end
      if isfield(G,'Rv'),  G.Ru = G.Rv;
      else, G.Ru = [];
      end
    case 'n'
      % No vessel - empty structures
      G.nu = 0; 
      G.Tvu = zeros(G.nv,G.nu);
      G.dimu = {''};
    otherwise
      error('invalid vessel mode')
  end
end
G.nu = size(G.Tvu,2); % nu size

% If Tivs is specified, compute Tius such that Is=Tius*Iu
if isfield(G,'Tivs')
  G.Tius = G.Tivs*G.Tvu; % this is identity if selu='s'
  G.ns   = size(G.Tius,1);
else
  G.Tius = zeros(0,G.nu); % no Is
  G.ns   = 0;
end

%% Mappings & sizes

% Defaults
if ~isfield(G,'Twa'), G.Twa = eye(numel(G.rw)); end
if ~isfield(G,'Tff'), G.Tff = eye(numel(G.rf)); end
if ~isfield(G,'Tmm'), G.Tmm = eye(numel(G.rm)); end

% Sizes
G.na = size(G.Twa,2);
G.nl = numel(G.rl);
G.nf = size(G.Tff,1);
G.nm = size(G.Tmm,1);

%% Get Green's function fields
for fields = gfields
 G.(fields{:}) = getg(G,fields{:});
end

%% Wall gaps
for wfield = {'rW','zW','aW','oW'}
  if isfield(P,wfield{:})
    G.(wfield{:}) = P.(wfield{:});
  end
end
  
if isfield(G,'rW')
  G.nW = numel(G.rW);
  assert(all(G.nW==[numel(G.aW),numel(G.rW),numel(G.zW),numel(G.oW)]),'all wall entries G.?W must have same size')
else
  G.nW=0;
end

%% Dipole matrix
if ~isfield(G,'Dda')
  G.Dda = -eye(G.na-1,G.na); G.Dda(G.na:G.na:end) = 1; % Diagonal (1) and upper diagonal (-1)
end

%%
G = orderfields(G);

%% remove rrx,zzx
G = rmfield(G,{'rrx','zzx'});

end

function val = getg(G,field)
% generic G field getter
if isfield(G,field), val=G.(field); return; end

% get auxiliary quantities as required
if any(contains(field,{'b','y','c'})), [lxy,lxb,lxc] = getlxybc(G); 
end
if any(contains(field,{'o','i'})), [lwi,lvi] = geti(G);
end

% short-circuit cases where no vessel is used
if strcmp(field(end),'u') && isempty(G.Tvu)
  o1 = field(end-1); if contains(o1,{'i','o'}), o1=field(end-2); end % other object
  switch o1
    case 'b',   no = sum(lxb(:));
    case 'x',   no = numel(G.rrx);
    case 'z',   no = numel(G.rrz);
    otherwise,  no = G.(sprintf('n%s',o1)); % number of rows needed
  end
  val = zeros(no,0); % directly assign empty value of correct size
  return
end

switch field
  case 'Mxv',  val = greenem('m' ,G.rrx,G.zzx,G.rv,G.zv);
  case 'Mxw',  val = greenem('m' ,G.rrx,G.zzx,G.rw,G.zw);
  case 'Mzv',  val = greenem('m' ,G.rrz,G.zzz,G.rv,G.zv);
  case 'Mzw',  val = greenem('m' ,G.rrz,G.zzz,G.rw,G.zw);
  case 'Mwv',  val = greenem('m' ,G.rw ,G.zw, G.rv,G.zv);
  case 'Brxv', val = greenem('br',G.rrx,G.zzx,G.rv,G.zv);
  case 'Bzxv', val = greenem('bz',G.rrx,G.zzx,G.rv,G.zv);
  case 'Brxw', val = greenem('br',G.rrx,G.zzx,G.rw,G.zw);
  case 'Bzxw', val = greenem('bz',G.rrx,G.zzx,G.rw,G.zw);   
  case 'Brzv', val = greenem('br',G.rrz,G.zzz,G.rv,G.zv);
  case 'Bzzv', val = greenem('bz',G.rrz,G.zzz,G.rv,G.zv);
  case 'Brzw', val = greenem('br',G.rrz,G.zzz,G.rw,G.zw);
  case 'Bzzw', val = greenem('bz',G.rrz,G.zzz,G.rw,G.zw);   
  case 'Mvv',  val = meqSelf(G.rv,G.zv,G.wv,G.hv,G.dv);
  case 'Mww',  val = meqSelf(G.rw,G.zw,G.ww,G.hw,G.dw);
  case 'Maa',  val = G.Twa'*getg(G,'Mww'   ) * G.Twa;
  case 'Mau',  val =        getg(G,'Mav'   ) * G.Tvu;
  case 'Mav',  val = G.Twa'*getg(G,'Mwv'   )        ;
  case 'Mxa',  val =        getg(G,'Mxw'   ) * G.Twa;
  case 'Muu',  val = G.Tvu'*getg(G,'Mvv'   ) * G.Tvu;
  case 'Mxu',  val =        getg(G,'Mxv'   ) * G.Tvu;
  case 'Mwu',  val =        getg(G,'Mwv'   ) * G.Tvu;
  case 'Mza',  val =        getg(G,'Mzw'   ) * G.Twa;
  case 'Mzu',  val =        getg(G,'Mzv'   ) * G.Tvu;
  case 'Mfa',  val =        getg(G,'Mfw'   ) * G.Twa;
  case 'Mfu',  val =        getg(G,'Mfv'   ) * G.Tvu;
  case 'Bma',  val =        getg(G,'Bmw'   ) * G.Twa;
  case 'Bmu',  val =        getg(G,'Bmv'   ) * G.Tvu;
  case 'Brxa', val =        getg(G,'Brxw'  ) * G.Twa;
  case 'Bzxa', val =        getg(G,'Bzxw'  ) * G.Twa;
  case 'Brxu', val =        getg(G,'Brxv'  ) * G.Tvu;
  case 'Bzxu', val =        getg(G,'Bzxv'  ) * G.Tvu;
  case 'Brza', val =        getg(G,'Brzw'  ) * G.Twa;
  case 'Bzza', val =        getg(G,'Bzzw'  ) * G.Twa;
  case 'Brzu', val =        getg(G,'Brzv'  ) * G.Tvu;
  case 'Bzzu', val =        getg(G,'Bzzv'  ) * G.Tvu;
  case 'dzMxa',val =        getg(G,'dzMxw' ) * G.Twa;
  case 'dzMax',val = G.Twa'*getg(G,'dzMwx' )        ;
  case 'dzMux',val = G.Tvu'*getg(G,'dzMvx' )        ;
  case 'dzMxu',val =        getg(G,'dzMxv' ) * G.Tvu;
  case 'drMxa',val =        getg(G,'drMxw' ) * G.Twa;
  case 'drMxu',val =        getg(G,'drMxv' ) * G.Tvu;
  case 'Mxia', val =        getg(G,'Mxiw'  ) * G.Twa( lwi,:); % from windings inside grid to grid
  case 'Mxiu', val =        getg(G,'Mxiv'  ) * G.Tvu( lvi,:);  % From vessel inside grid to grid  
  case 'Mbou', val =        getg(G,'Mbov'  ) * G.Tvu(~lvi,:);
  case 'Mboa', val =        getg(G,'Mbow'  ) * G.Twa(~lwi,:);
  case 'Brboa',val =        getg(G,'Brbow' ) * G.Twa(~lwi,:); % From windings outside grid to boundary
  case 'Brbou',val =        getg(G,'Brbov' ) * G.Tvu(~lvi,:); % From vessel outside grid to boundary
  case 'drBrxa',val =       getg(G,'drBrxw') * G.Twa;
  case 'dzBrxa',val =       getg(G,'dzBrxw') * G.Twa;
  case 'drBzxa',val =       getg(G,'drBzxw') * G.Twa;
  case 'dzBzxa',val =       getg(G,'dzBzxw') * G.Twa;
  case 'drBrxu',val =       getg(G,'drBrxv') * G.Tvu;
  case 'dzBrxu',val =       getg(G,'dzBrxv') * G.Tvu;
  case 'drBzxu',val =       getg(G,'drBzxv') * G.Tvu;
  case 'dzBzxu',val =       getg(G,'dzBzxv') * G.Tvu;
  case 'drdrMxa',val =       getg(G,'drdrMxw') * G.Twa;
  case 'drdzMxa',val =       getg(G,'drdzMxw') * G.Twa;
  case 'dzdrMxa',val =       getg(G,'dzdrMxw') * G.Twa;
  case 'dzdzMxa',val =       getg(G,'dzdzMxw') * G.Twa;
  case 'drdrMxu',val =       getg(G,'drdrMxv') * G.Tvu;
  case 'drdzMxu',val =       getg(G,'drdzMxv') * G.Tvu;
  case 'dzdrMxu',val =       getg(G,'dzdrMxv') * G.Tvu;
  case 'dzdzMxu',val =       getg(G,'dzdzMxv') * G.Tvu;
  case 'Ra',     val =       diag(G.Twa' * diag(G.Rw) * G.Twa);
  case 'Mbw'
    if isfield(G,'Mxw'), val = G.Mxw(lxb,:);
    else,                val = greenem('m',G.rrx(lxb),G.zzx(lxb),G.rw,G.zw);
    end
  case 'Mbv'
    if isfield(G,'Mxv'), val = G.Mxv(~lxy,:);
    else,                val = greenem('m',G.rrx(lxb),G.zzx(lxb),G.rv,G.zv);
    end
  case 'Mby'
    if isfield(G,'Mxx'), val = G.Mxx(~lxy,lxy);
    else,                val = greenem('m',G.rrx(lxb),G.zzx(lxb),G.rrx(lxy),G.zzx(lxy));
    end
  case 'Mbc'
    if isfield(G,'Mxx'), val = G.Mxx(~lxy,lxc);
    else,                val = greenem('m',G.rrx(lxb),G.zzx(lxb),G.rrx(lxc),G.zzx(lxc));
    end
  case 'Mfv',  val = G.Tff*greenem('m',G.rf,G.zf,G.rv,G.zv);
  case 'Mfw',  val = G.Tff*greenem('m',G.rf,G.zf,G.rw,G.zw);
  case 'Mfx',  val = G.Tff*greenem('m',G.rf,G.zf,G.rrx,G.zzx);
  case 'Bmv'
               val = G.Tmm*(cos(G.am).*greenem('br'   ,G.rm,G.zm,G.rv,G.zv)  +  ...
                            sin(G.am).*greenem('bz'   ,G.rm,G.zm,G.rv,G.zv));
  case 'Bmw'
               val = G.Tmm*(cos(G.am).*greenem('br'   ,G.rm,G.zm,G.rw,G.zw)  +  ...
                            sin(G.am).*greenem('bz'   ,G.rm,G.zm,G.rw,G.zw));
  case 'Bmx'
               val = G.Tmm*(cos(G.am).*greenem('br'   ,G.rm,G.zm,G.rrx,G.zzx)  +  ...
                            sin(G.am).*greenem('bz'   ,G.rm,G.zm,G.rrx,G.zzx));
  case 'dzMvx',   val = greenem('dz1m',G.rrx,G.zzx,G.rv,G.zv)';
  case 'dzMwx',   val = greenem('dz1m',G.rrx,G.zzx,G.rw,G.zw)';                       
  case 'dzMxv',   val = greenem('dz1m',G.rrx,G.zzx,G.rv,G.zv);
  case 'dzMxw',   val = greenem('dz1m',G.rrx,G.zzx,G.rw,G.zw);
  case 'drMxv',   val = greenem('dr1m',G.rrx,G.zzx,G.rv,G.zv);
  case 'drMxw',   val = greenem('dr1m',G.rrx,G.zzx,G.rw,G.zw);
  case 'drBrxw',  val = greenem('dr1br',G.rrx,G.zzx,G.rw,G.zw);
  case 'dzBrxw',  val = greenem('dz1br',G.rrx,G.zzx,G.rw,G.zw);
  case 'drBzxw',  val = greenem('dr1bz',G.rrx,G.zzx,G.rw,G.zw);
  case 'dzBzxw',  val = greenem('dz1bz',G.rrx,G.zzx,G.rw,G.zw);
  case 'drBrxv',  val = greenem('dr1br',G.rrx,G.zzx,G.rv,G.zv);
  case 'dzBrxv',  val = greenem('dz1br',G.rrx,G.zzx,G.rv,G.zv);
  case 'drBzxv',  val = greenem('dr1bz',G.rrx,G.zzx,G.rv,G.zv);
  case 'dzBzxv',  val = greenem('dz1bz',G.rrx,G.zzx,G.rv,G.zv);
  case 'drdrMxw', val = greenem('dr1r1m',G.rrx,G.zzx,G.rw,G.zw);
  case 'drdzMxw', val = greenem('dr1z1m',G.rrx,G.zzx,G.rw,G.zw);
  case 'dzdrMxw', val = greenem('dz1r1m',G.rrx,G.zzx,G.rw,G.zw);
  case 'dzdzMxw', val = greenem('dz1z1m',G.rrx,G.zzx,G.rw,G.zw);
  case 'drdrMxv', val = greenem('dr1r1m',G.rrx,G.zzx,G.rv,G.zv);
  case 'drdzMxv', val = greenem('dr1z1m',G.rrx,G.zzx,G.rv,G.zv);
  case 'dzdrMxv', val = greenem('dz1r1m',G.rrx,G.zzx,G.rv,G.zv);
  case 'dzdzMxv', val = greenem('dz1z1m',G.rrx,G.zzx,G.rv,G.zv); 
  case 'dzMfx',   val = getg(G,'Tff') *   greenem('dz1m' ,G.rf,G.zf,G.rrx,G.zzx);
  case 'dzBmx',   val = getg(G,'Tmm') * (cos(G.am).*greenem('dbrdz',G.rm,G.zm,G.rrx,G.zzx) +  ...
                                        sin(G.am).*greenem('dbzdz',G.rm,G.zm,G.rrx,G.zzx));
  case 'Mbow'
    if isfield(G,'Mxw'),  val = G.Mxw(lxb,~lwi);
    else,                 val = greenem('m',G.rrx(lxb),G.zzx(lxb),G.rw(~lwi),G.zw(~lwi));
    end
  case 'Mbov' 
    if isfield(G,'Mxv'),  val = G.Mxv(lxb,~lvi);
    else,                 val = greenem('m',G.rrx(lxb),G.zzx(lxb),G.rv(~lvi),G.zv(~lvi));
    end
  case 'Brbow'
    if isfield(G,'Brxw'), val = G.Brxw(lxb,~lwi);
    else,                 val = greenem('br',G.rrx(lxb),G.zzx(lxb),G.rw(~lwi),G.zw(~lwi));
    end
  case 'Brbov' 
    if isfield(G,'Brxv'), val = G.Brxv(lxb,~lvi);
    else,                 val = greenem('br',G.rrx(lxb),G.zzx(lxb),G.rv(~lvi),G.zv(~lvi));
    end
  case 'Mxiv'
    if ~isfield(G,'Mxv'), val = greenem('mut',G.rrx,G.zzx,G.rv(lvi),G.zv(lvi));
    else,                 val = G.Mxv(:,lvi);
    end
  case 'Mxiw'
    if ~isfield(G,'Mxw'), val = greenem('mut',G.rrx,G.zzx,G.rw(lwi),G.zw(lwi));
    else,                 val = G.Mxw(:,lwi);
    end
  case 'Mzy',             val = greenem('mut',G.rrz,G.zzz,G.rrx(lxy),G.zzx(lxy));
  otherwise
    error('undefined variable request %s for %s',field,mfilename);
end
end

function M = meqSelf(r,z,w,h,d)
M = greenem('mut'  ,r ,z,[w,h]);
icirc = find(isnan(w) & isnan(h)); % not rectangular windings
if any(icirc), assert(nargin>4), else, return; end
isub = sub2ind(size(M),icirc,icirc);
M(isub) = greenem('self'  ,r(icirc) ,z(icirc),d(icirc)); % substitute self
end

function xargs = sortxargs(xargs)
% Order  optional fields so following fields can re-use previous results
% e.g. 'Mua' is fast to compute if 'Mvw' is known;
%
% First   those without 'a','u','c','b' (these can be computed from others
% Then    those with    'a','u','c','b'  and also 'v','w','x'
% Finally those with    'a','u','c','b'  and   no 'v','w','x'
xau = contains(xargs,{'a','u','c','b'}); % xargs with 'a', 'u', 'c' or 'b'
xvw = contains(xargs,{'v','w','x'}); % xargs with 'v' or 'w' or 'x'
x1 = ~xau; x2 = xvw & xau; x3 = ~xvw & xau;
xargs = [xargs(x1),xargs(x2),xargs(x3)];
end

function [lxy,lxb,lxc] = getlxybc(G)
lxy = false(size(G.rrx)); lxy(2:end-1,2:end-1) = true ;
[lxb,lxc] = deal(~lxy); lxc([1 end],[1 end]) = false;
end

function [lwi,lvi] = geti(G)
lwi = G.rw > G.rx(1) & G.rw < G.rx(end) & G.zw > G.zx(1) & G.zw < G.zx(end); % windings inside grid
lvi = G.rv > G.rx(1) & G.rv < G.rx(end) & G.zv > G.zx(1) & G.zv < G.zx(end); % vessel filaments inside grid
end

function [rx,zx] = meqgx(P,G)
% Sets computational grid from parameters in P and limiter contour in G if available

% Check for FE settings
hasFE = isfield(P,'rip');
if hasFE, [rip,rop,zlp,zup] = deal(P.rip,P.rop,P.zlp,P.zup);
else,     [rip,rop,zlp,zup] = deal(NaN);
end
% Check if limiter is defined
haslim = all(isfield(G,{'rl','zl'}));
%
nz = P.nz;
if haslim
  ril = min(G.rl);
  rol = max(G.rl);
  zll = min(G.zl);
  zul = max(G.zl);
  % Solve for [zl;zu] assuming half-grid offset with respect to limiter if unknown
  A = eye(2)+1/(4*nz)*[-1,1;1,-1]; b = [zll;zul];
  if ~isnan(P.zl), A(1,:) = [1,0]; b(1) = P.zl; end % zl is known
  if ~isnan(P.zu), A(2,:) = [0,1]; b(2) = P.zu; end % zu is known
  x = A\b; zl = x(1); zu = x(2);
  % Same for [ri;ro]
  if isnan(P.nr)
    dz = (zu-zl)/(2*nz);
    A = eye(2); b = [ril;rol]+dz/2*[-1;1];
  else
    A = eye(2)+1/(2*P.nr)*[-1,1;1,-1]; b = [ril;rol];
  end
  if ~isnan(P.ri), A(1,:) = [1,0]; b(1) = P.ri; end % ri is known
  if ~isnan(P.ro), A(2,:) = [0,1]; b(2) = P.ro; end % ro is known
  x = A\b; ri = x(1); ro = x(2);
else
  % Without a limiter, use FE limits
  if ~isnan(P.zl), zl = P.zl; else, zl = zlp; end
  if ~isnan(P.zu), zu = P.zu; else, zu = zup; end
  if ~isnan(P.ri), ri = P.ri; else, ri = rip; end
  if ~isnan(P.ro), ro = P.ro; else, ro = rop; end
end
assert(~any(isnan([zl,zu,ri,ro])),'Computed grid boundaries contain NaN values, check parameters');
  
% Compute nr if needed
if isnan(P.nr), nr = (ro - ri)/(zu - zl)*(2*nz);
else,           nr = P.nr;
end
  
% create rx,zx
rx = linspace(ri,ro,  nr+1)';
zx = linspace(zl,zu,2*nz+1)';
end
