function G = meqgana(P)
% Geometry for anamak
% See also MEQPANA
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

rhoCu = 1.72e-8; % resistivity of Copper
rhoSS = 6.9e-7; % resistivity of Stainless Steel

if P.fastcoil, nc = P.nc-1;
else,          nc = P.nc;
end

% theta grids
thv = linspace(0,2*pi,P.nv+1)'; thv = thv(1:end-1);
thl = linspace(0,2*pi,P.nl+1)'; thl = thl(1:end-1);
thf = linspace(0,2*pi,P.nf+1)'; thf = thf(1:end-1);
thm = linspace(0,2*pi,P.nm+1)'; thm = thm(1:end-1);
thc = linspace(0,2*pi,  nc+1)'; thc = thc(1:end-1);
thW = linspace(0,2*pi,P.nW+1)'; thW = thW(1:end-1);

% vessel
G.rv = P.r0 +          P.av*cos(thv);
G.zv =    0 + P.cappav*P.av*sin(thv);
G.wv = mean(sqrt(diff(G.rv).^2 + diff(G.zv).^2)) * ones(P.nv,1);
G.hv = G.wv; % Square conductors
G.Rv = rhoSS * 2*pi*G.rv./G.wv.^2; % Copper vessel

% limiter
G.rl = P.r0 +          P.al*sin(thl);
G.zl =    0 + P.cappav*P.al*cos(thl);

% windings/coils
G.rw = P.r0 +          P.ac*cos(thc);
G.zw =    0 + P.cappav*P.ac*sin(thc);
G.dw = nan(size(G.rw));
G.ww = 0.01*ones(size(G.rw)); % winding width (one winding per coil for now)
G.hw = G.ww; % Square conductors
G.Rw = rhoCu * 2*pi*G.rw./G.ww.^2; % Copper winding resistance
G.Twa = eye(numel(G.rw)); % interconnection matrix windings-coils

if P.fastcoil
  G.rw = [G.rw; P.r0 + [1; 1]*sqrt(2)/2*P.av*0.98];
  G.zw = [G.zw;    0 + [1;-1]*sqrt(2)/2*P.av*0.98];
  G.dw = [G.dw;G.ww(end)*[1;1]];
  G.ww = [G.ww;nan ;nan ];
  G.hw = [G.hw;nan ;nan ];
  G.Rw = [G.Rw;rhoCu*2*pi*G.rw(end-1:end)./(pi*G.dw(end-1:end).^2)];
  antiseries = zeros(numel(G.rw),1); antiseries(end-1:end)=[1;-1];
  G.Twa = [[G.Twa;zeros(2,size(G.Twa,2))],antiseries];
end

G.rc = G.rw; G.zc = G.zw;
G.hc = P.wc*ones(size(G.rc)); 
G.wc = P.wc*ones(size(G.rc));

G.Tca = G.Twa; 

% flux loops
G.rf = P.r0 + P.af * cos(thf);
G.zf =    0 + P.af * sin(thf) * P.cappav;
G.rm = P.r0 + P.am * cos(thm);
G.zm =    0 + P.am * sin(thm) * P.cappav;
G.am =      thm+pi/2; % direction

% Wall gaps
G.rW  = P.r0 +          P.al*cos(thW);
G.zW  =    0 + P.cappav*P.al*sin(thW);
G.oW  = -thW;
G.aW  = repmat(2*P.av,size(thW));

% dimension labels
if ~P.fastcoil
  G.dima = cellstr(num2str((1:size(G.Twa,2))','PF_%03d'));
else
  G.dima = [cellstr(num2str((1:size(G.Twa,2)-1)','PF_%03d'));{'VS_001'}];
end
G.dimc = cellstr(num2str((1:size(G.Tca,1))','Coil_%03d'));
G.dimf = cellstr(num2str((1:numel(G.rf))','Ff_%03d'));
G.dimm = cellstr(num2str((1:numel(G.rm))','Bm_%03d'));
G.dimv = cellstr(num2str((1:numel(G.rv))','vv_%03d'));
G.dimW = cellstr(num2str((1:numel(G.rW))','gap_%03d'));

% vessel specification needed in liu
if isfield(P,'selu')
  switch P.selu
    case 's'
      if ~isfield(P,'nu'), P.nu = P.nv; end
      G.Tvs = eye(P.nu); G.dims=cellstr(num2str((1:numel(P.nu))','vs_%03d'));
  end
end

% Power supplies
G.Valim   = inf(numel(G.dima),1);
G.Ialim   = inf(numel(G.dima),1);
G.Talim   = inf(numel(G.dima),1);
end
