function G = meqgload(G,fname,verb)
%MEQGLOAD Load geometry data from file to avoid recomputing it
% G = meqgload(G,fname,verbosity)
%
% INPUTS:
%    G: Geometry structure. The geometry data in the file is compared with
%    this structure to determine compatibility.
%    For example Mvv will be loaded only if rv,zv,hv,wv are compatible 
%    between the stored file and G.
%    fname: File name to be loaded
%    v: verbosity. 0: silent 1:basic 2:list loaded variables
% OUTPUTS:
%    G: Updated geometry structure with additional loaded fields
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin<=2; verb=0; end % verbosity

Gs = load(fname);

% find which quantities are equal between loaded and target values
a = isfield(G,'Twa') && isfield(Gs,'Twa') && isequal(Gs.Twa,G.Twa); % circuits equal
u = isfield(G,'Tvu') && isfield(Gs,'Tvu') && isequal(Gs.Tvu,G.Tvu); % vessel mode equal
w = paramsequal(Gs,G,'w'); % windings equal
v = paramsequal(Gs,G,'v'); % vessel equal
x = paramsequal(Gs,G,'x'); % computational grid equal
z = paramsequal(Gs,G,'z'); % extended grid equal
m = paramsequal(Gs,G,'m'); % magnetic probes equal
f = paramsequal(Gs,G,'f'); % flux loops equal

if verb
  name = ['a','u','w','v','x','z','m','f'];
  mask = [ a , u , w , v , x , z , m , f ];

  comp = name(mask); ncomp=name(~mask);
  
  if ~isempty( comp), fprintf('meqgload:     loading compatible              %s\n', comp); end
  if ~isempty(ncomp), fprintf('meqgload: not loading incompatible or missing %s\n',ncomp); end
end
      
% Optional loading 
if v, G=doload('Mvv',G,Gs,verb);  end
if w, G=doload('Mww',G,Gs,verb);  end
if a, G=doload('Maa',G,Gs,verb);  end

if w && v, G=doload({'Mwv'},G,Gs,verb); end
if w && u, G=doload({'Mwu'},G,Gs,verb); end
if a && v, G=doload({'Mav'},G,Gs,verb); end
if a && u, G=doload({'Mau'},G,Gs,verb); end

if m && x, G=doload('Bmx',G,Gs,verb); end
if f && x, G=doload('Mfx',G,Gs,verb); end
if m && a, G=doload('Bma',G,Gs,verb); end
if f && a, G=doload('Mfa',G,Gs,verb); end
if m && v, G=doload('Bmv',G,Gs,verb); end
if f && v, G=doload('Mfv',G,Gs,verb); end
if m && u, G=doload('Bmu',G,Gs,verb); end
if f && u, G=doload('Mfu',G,Gs,verb); end

if x && v , G=doload({'Mxv','Brxv','Bzxv','dzMvx','Mbv'},G,Gs,verb); end
if x && w , G=doload({'Mxw','Brxw','Bzxw','dzMwx','Mbw'},G,Gs,verb); end
if x && a , G=doload({'Mxa','Brxa','Bzxa','dzMax'      },G,Gs,verb); end
if x && u , G=doload({'Mxu','Brxu','Bzxu','dzMux'      },G,Gs,verb); end

if z && v , G=doload({'Mzv','Brzv','Bzzv'},G,Gs,verb); end
if z && w , G=doload({'Mzw','Brzw','Bzzw'},G,Gs,verb); end
if z && a , G=doload({'Mza','Brza','Bzza'},G,Gs,verb); end
if z && u , G=doload({'Mzu','Brzu','Bzzu'},G,Gs,verb); end
if z && x , G=doload({'Mzy'              },G,Gs,verb); end

end

function eq = paramsequal(G1,G2,element)
  eq = true;
  for val={'r','z','w','h','d'}
    myfield = [val{:},element];
    eq = eq && ...
      (...
      (~isfield(G1,myfield) && ~isfield(G2,myfield)) || ...
      (isfield(G1,myfield) && isfield(G2,myfield) && isequaln(G1.(myfield),G2.(myfield)))...
      );
  end
end

function Gtarget = doload(fields,Gtarget,Gsrc,verb)
% Load required fields from Gsrc to Gtarget
if ~iscell(fields), fields={fields}; end
for ifield=1:numel(fields)
  field = fields{ifield};
  if isfield(Gsrc,field)
    Gtarget.(field) = Gsrc.(field);
  end
  if (verb-1)>0, fprintf('loaded %s\n',field); end
end
end
