function [MG,FB,X,Xref,S] = mgaruncmp(cmd,shot,varargin) 
% function [MG,FB,X,Xref,S] = mgaruncmp(cmd,shot,fort_equiv,tol,verbose) 
% MGA.F vs MGA.M comparisons

 %% Argument Parsing
 p = inputParser;
 p.addRequired('cmd');
 p.addRequired('shot');
 p.addOptional('fort_equiv',false);
 p.addOptional('tol'       ,  0.1);
 p.addOptional('verb'      ,false);
 p.addOptional('saveshot'  ,mgu.mytestshot());
 p.addParameter('fbtem'     ,false);
 
 p.parse(cmd,shot,varargin{:});
 
 cmd        = p.Results.cmd;
 shot       = p.Results.shot;
 fort_equiv = p.Results.fort_equiv;
 tol        = p.Results.tol;
 verb       = p.Results.verb;
 saveshot   = p.Results.saveshot;
 fbtem      = p.Results.fbtem;
 
 %% Function body
 shotdesign_dir = fileparts(mfilename('fullpath')); %Path to the current shotdesign version
 
 % load MG parameters
 fprintv(verb,'loading parameters: mga(%d)\n',shot);
 try
  MG = mgp(shot);
  if ~MG.check, MG.fix;end
 catch ME
  if startsWith(ME.message,'MatrixGenerator:mdso')
    error('mgaruncmd:noShotFiles','MDS files for shot %d could not be open',shot);
  else
    rethrow(ME);
  end
 end
 
 % Flag cases where comparison is difficult
 assert(~fbtem || (MG.nvvel == 19 && MG.nfast <= 0)   ,'mgaruncmd:fbtmeasMissing'    ,'nvvel~=19 or nfast > 0 are incompatible with mgaruncmd running fbte.m');
 assert(abs(MG.ohcorr) < 1e4                          ,'mgaruncmd:largeOhcorr'       ,'Large ohcorr, skipping');
 assert(MG.iscramb >= 0                               ,'mgaruncmd:voltageControl'    ,'Voltage control shot (ISCRAMB<0), skipping');
 assert(~ismember(mod(MG.mvloop,10),[5,6])            ,'mgaruncmd:unsupportedMvloop' ,'unsupported mvloop=%d',MG.mvloop);
 assert(numel(MG.toft) == 14 && numel(MG.diohdt) == 14,'mgaruncmd:unsupportedToft'   ,'unsupported toft/diohdt input');
 assert(ismember(MG.iscramb,[-2:0 2:12 61:63])        ,'mgaruncmd:unsupportedIscramb','iscramb=%d',MG.iscramb)
 assert(MG.mvloop ~= 8                                ,'mgaruncmd:doubletShot'       ,'Doublet shot (MIDPLAN=8), skipping');
 
 if shot >= 45000
  Xsav = mgaget(shot);
  if std(diff(Xsav.TOH)) < 1e-6
   mgaver = 'm';
  else
   mgaver = 'f';
  end
  fprintv(verb,'#%d was prepared with mga.%1s\n',shot,mgaver)
 else
  fprintv(verb,'MDS storage of MGA results is unavailable for shot #%d\n',shot);
 end
 
 fprintv(verb,'save to personal shot %d\n',saveshot);
 MG.save(saveshot)
 
 if fbtem
  [~,~,LY] = fbte(dest_shot,'itamax',100);
  assert(numel(LY.t) == MG.nruns && all(LY.isconverged),'mgaruncmd:noconvFBT','FBT did not converge for at least one run, cannot continue');
  
  % fbtmeas.dat needed for MGA.F
  testsdir = fullfile(fileparts(mfilename('fullpath')),'tests');
  copyfile(fullfile(testsdir,'fbtmeas.dat'),'.','f');
 else
  % run FBTE (.F)
  fbtcmd = sprintf('%s/fbt %d > fbt.log 2>&1',shotdesign_dir,saveshot);
  cleanupobj = onCleanup(@() delete('fbt.log'));
  fprintv(verb,'running %s\n',fbtcmd);
  assert(system(fbtcmd) == 0,'mgaruncmd:errorFBTFortran','Error when running fbt');
  if system('grep "no convergence" fbt.log -q') ~= 1
   fprintv(verb,'FBT did not converge for at least one run, will try once more with increased number of iterations\n');
   MG.itamax(:) = 2*max(MG.itamax);
   MG.save(saveshot);
   assert(system(fbtcmd) == 0,'mgaruncmd:errorFBTFortran','Error when running fbt');
   assert(system('grep "no convergence" fbt.log -q') == 1,'mgaruncmd:noConvFBT','FBT did not converge for at least one run, cannot continue');
  end
 end
 assert(logical(exist('fbtmeas.dat','file')),'fbtmeas.dat missing, cannot continue');
 
 % run MGA.F
 mgacmd = sprintf('%s/mga %d > mga.log 2>&1',shotdesign_dir,saveshot);
 cleanupobj2 = onCleanup(@() delete('mga.log'));
 fprintv(verb,'running %s\n',mgacmd);
 assert(system(mgacmd) == 0,'mgaruncmd:errorMGAFortran','Error when running mga');
 Xref = mgaget(saveshot);

 % run MGA.M
 fprintv(verb,'running matlab mga(%d)\n',saveshot);
 [~,FB,X] = mga(saveshot,fort_equiv);
 
 assert(X.MIDPLAN  ~= 8,'mgaruncmd:doubletShot','Doublet shot (MIDPLAN=8), skipping');
 
 % Compare
 fprintv(verb,'analyzing...\n\n');
 [~,S] = mgacmp(X,Xref,cmd,tol,tol,verb);
 
end

function fprintv(verb,varargin)
% verbose printing
if verb
  fprintf(varargin{:});
end
end
