function mgams(mode,arg1,arg2,arg3,arg4)
 % MGAMS launches GUI
 
 global H             % GUI handles                     [struct]
 global MGclip        % Clipboard for equil. copy/paste [struct]
 global MGceq         % Current equilibrium index       [double]
 global MGgas         % Current gas++ traces            [struct]
 global MGgrp         % Clipboard for mgams groups      [struct]
 global CBgas         % Clipboard for gas               [struct]
 global MP            % Current mgp instance            [mgp]
 global MGerr         % object for MGAMS GUI errors     [MException]
 
 persistent CBgas_bd  % Clipboard for BD gas            [struct]
 global FBTE_eq
 
 NUVAR = 10; % number of variables in user selected list
 NGVAR = 20; % number of variables in one group
 
 %% Initialisation
 if nargin == 0 % init
  
  % Parameter initialisation
  MGclip.neq_ = 0;
  
  % GUI handles
  H.tplt = zeros(1,5); % time plot axes
  H.bndl = zeros(1,8); % x-section lines
  H.tlin = zeros(1,10); % time plot lines
  H.tlab = []; % time plot labels
  [H.varl,H.vare] = deal(zeros(1,NUVAR+NGVAR)); % variable name labels and fields
  [H.dlgl,H.dlge] = deal(zeros(1,4)); % modal dialog box labels and fields
  [H.msge,... % message text
   H.enum,... % EQ number
   H.shot,... % shot number
   H.time,... % EQ time
   H.varm,... % variable menu
   H.xplt,... % x-section axes
   H.dlgb,... % modal dialog box button
   H.cpos] = deal(0); % x-section cursor position text
  
  % Variables
  MP = mgp;
  
  % Error
  MGerr = MException('','');
  
  % groups
  MGgrp(1).name = 'GPS';
  MGgrp(1).var = {'gpr','gpz','gpb','gpfa','gpfb','gpfe','gpbr','gpbz','gpba','gpbe','gpcr','gpcz','gpca','gpce','gpia','gpie','gpfd','gpbd','gpcd','gpid'};
  MGgrp(2).name = 'Feedback';
  MGgrp(2).var = {'gainz','gainvz','aipgain','gain','ggain','hgain','gainr','fastm','nfast'};
  MGgrp(3).name = 'Breakdown';
  MGgrp(3).var = {'aipdel','aipipz','delipz','diohdt','dpszero','brmzero','brpzero','bzmzero','bzpzero','hpla','ierat','if36fb','iohfb','ohcorr','psifac','toft','wpla'};
  MGgrp(4).name = 'Vertical';
  MGgrp(4).var = {'iscramb','zup','zlp','nzaxel','wscal','vscal','uc2','ipr16','nzaxre','nfast','ikriz','inova','delipz','aipipz','zeecorr','ohcorr'};
  MGgrp(5).name = 'Elongation';
  MGgrp(5).var = {'ifour','iscramb','gainext','nelz','nzaxre','zlp','zup'};
  MGgrp(6).name = 'SCD ne control';
  MGgrp(6).var = {'dneon','dnekp','dneki','dnekd','dneka','dnefil','dnemode','dnescal1','dnescal2'};
  MGgrp(7).name = 'Debug nodes';
  MGgrp(7).var = strcat('debug_p',cellstr(int2str((1:9).')).'); % debug_p1 ... debug_p9
  
  H.fig = figure('position',[10 10 935 865],'name','MGAMS - Matrix Generator','tag','MGAMS','numbertitle','off','resize','off','menubar','none','windowbuttondownfcn','mgams mouse');
  set(H.fig,'color','k',...
    'defaultaxescolororder','default','defaultaxescolor','k','defaultaxesxcolor','w','defaultaxesycolor','w','defaultaxesfontsize',9,...
    'defaultuicontrolhorizontalalignment','left','defaultuicontrolfontsize',9,...
    'defaulttextfontsize',9)
  
  % menus
  mbcol = [.7 .7 1];
  list = {
   -1		   'Model           '
   100000	'Standard (Ip<0 B<0 z=+0.23)'
   100001	'Ip<0 B<0 z= 0.00'
   100002	'Ip<0 B<0 z=-0.23'
   100003	'Ip>0 B<0 z=+0.23'
   100004	'Ip>0 B<0 z= 0.00'
   };
  h(1) = uimenu('label','Shot');
  h(2) = uimenu(h(1),'label','Load');
  uimenu(h(2),'label','Shot #','callback','mgams load','accelerator','L');
  for k = 1:size(list,1)
   uimenu(h(2),'label',list{k,2},'callback',['mgams load ' int2str(list{k,1})])
  end
  h(2) = [];
  uimenu(h(1),'label','Save in model','callback','mgams save -1','accelerator','M','foregroundcolor','r')
  uimenu(h(1),'label','Save in shot #','callback','mgams save','accelerator','S','foregroundcolor','r')
  uimenu(h(1),'label','Prepare shot','callback','mgams prepare','accelerator','P','foregroundcolor','r')
  uimenu(h(1),'label','Blank recipe','callback','mgams blank_recipe','accelerator','B','foregroundcolor','r')
  h(2) = uimenu(h(1),'label','Configure preparation');
  uimenu(h(2),'label','FBTE-Fortran','checked','off','userdata','s = unix(''fbt # > /tmp/$USER/fbt-#.log && ! grep "no convergence" /tmp/$USER/fbt-#.log '');')
  uimenu(h(2),'label','PSITBX'      ,'checked','off','userdata','s = 0; psitbxput(1,#)')
  uimenu(h(2),'label','FBTE-Matlab' ,'checked','on' ,'userdata','s = 0; s = run_fbte(#);')
  %
  uimenu(h(2),'label','MGA-Fortran'                  ,'checked','off','userdata','s = 0; s = run_mga(#,''Fortran'');','Separator','on')
  uimenu(h(2),'label','MGA-Matlab compatibility mode','checked','on' ,'userdata','s = 0; s = run_mga(#,''Matlab compatibility mode'');')
  uimenu(h(2),'label','MGA-Matlab'                   ,'checked','off','userdata','s = 0; s = run_mga(#,''Matlab'');')
  uimenu(h(2),'label','MGA-Matlab checking vs .F'    ,'checked','off','userdata','s = 0; s = run_mga(#,''Matlab checked'');')
  %
  uimenu(h(2),'label','FBTE18-Matlab'            ,'checked','on' ,'userdata','s = 0; fbte18(#)','Separator','on')
  uimenu(h(2),'label','FBTE18-Matlab, rtdtacq'   ,'checked','off','userdata','s = 0; fbte18(#,108,''rtdtacq'')')
  uimenu(h(2),'label','FBTE18-Matlab, rtdtacqden','checked','off','userdata','s = 0; fbte18(#,108,''rtdtacqden'')')
  uimenu(h(2),'label','FBTE18-Matlab, rtdtacqfps','checked','off','userdata','s = 0; fbte18(#,108,''rtdtacqfps'')')
  %
  uimenu(h(2),'label','Blank-Matlab'       ,'checked','on' ,'userdata','s = tcv_recipe(''fbte'',#);','Separator','on')
  %
  set(get(h(2),'children'),'tag','prepconfig','callback','set(gcbo,''checked'',ifc(strcmp(get(gcbo,''checked''),''on''),''off'',''on''))')
  uimenu(h(2),'label','Restore default selection','callback','mgams restore_prep_config','Separator','on');
  for k = 1:mgu.size
   l = uimenu(h(1),'label',sprintf('#%6d %s',mgu.shot1(k),mgu.user(k)));
   if k == 3, set(l,'separator','on'), end
   if strcmp(mgu.user(k),getenv('USER')), set(l,'foregroundcolor',[0 .5 0]), end
  end
  uimenu(h(1),'label','Show FBTE equilibria','checked','off','accelerator','H','separator','on',...
    'tag','show_fbte','callback','set(gcbo,''checked'',ifc(strcmp(get(gcbo,''checked''),''on''),''off'',''on''));mgams equil select;')
  uimenu(h(1),'label','Type','callback','mgams type','separator','on')
  uimenu(h(1),'label','Print','callback','mgams print')
  uimenu(h(1),'label','Quit','callback','mgams quit','accelerator','Q','separator','on')
  h(1) = uimenu('label','Equilibrium');
  uimenu(h(1),'label','Copy...','callback','mgams copy');
  uimenu(h(1),'label','Paste...','callback','mgams paste');
  uimenu(h(1),'label','Duplicate...','callback','mgams duplicate');
  uimenu(h(1),'label','Delete...','callback','mgams delete');
  uimenu(h(1),'label','Interpolate...','callback','mgams interpolate');
  uimenu(h(1),'label','Ramp down...','callback','mgams rampdown');
  uimenu(h(1),'label','Shift...','callback','mgams shift','separator','on');
  uimenu(h(1),'label','Scale...','callback','mgams scale');
  uimenu(h(1),'label','Mirror...','callback','mgams mirror');
  uimenu(h(1),'label','Discretise...','callback','mgams discretise');
  uimenu(h(1),'label','Half analytic...','callback','mgams halfanalytic');
  uimenu(h(1),'label','Analytic fit...','callback','mgams analyticfit');
  uimenu(h(1),'label','Set point...','callback','mgams setpoint');
  uimenu(h(1),'label','Plasma current ramp...','callback','mgams currentramp','separator','on');
  uimenu(h(1),'label','Set absolute time...','callback','mgams setabsolutetime');
  uimenu(h(1),'label','Set flat top','callback','mgams setflattop','accelerator','F','separator','on')
  h(1) = uimenu('label','Gas/Density');
  uimenu(h(1),'label','Copy','callback','mgams gascopy');
  uimenu(h(1),'label','Paste','callback','mgams gaspaste');
  uimenu(h(1),'label','Proportional to Plasma Current','callback','mgams gasip');
  uimenu(h(1),'label','Copy Gas Breakdown','callback','mgams gascopy_bd','separator','on');
  uimenu(h(1),'label','Paste Gas Breakdown','callback','mgams gaspaste_bd');
  h(1) = uimenu('label','Macro');
  uimenu(h(1),'label','Use FPS...','callback','mgams usefps');
  uimenu(h(1),'label','Do not use FPS...','callback','mgams nofps');
  uimenu(h(1),'label','Copy ramp up equilibria...','callback','mgams copyrampup');
  uimenu(h(1),'label','Paste ramp up equilibria...','callback','mgams pasterampup');
  uimenu(h(1),'label','Activate SCD ne controller...','callback','mgams activate_dne');
  H.varm = uimenu('label','Variables');
  for k = 1:NUVAR
   h(1) = uimenu(H.varm,'label',['Variable ' int2str(k)]);
   col = 'kb';
   for ka = 'a':'z'
    h(2) = uimenu(h(1),'label',char(ka));
    for l = find(strncmp(char(ka),MP.name,1))'
     uimenu(h(2),'label',MP.name(l),'foregroundcolor',col(MP.eqdp(l)+1),'callback',sprintf('mgams(''var'',%d,%d)',k,l))
    end
   end
  end
  h(1) = uimenu('label','Groups');
  for k = 1:numel(MGgrp)
   uimenu(h(1),'label',MGgrp(k).name,'callback',sprintf('mgams(''grp'',%d)',k));
  end
  H.fbtplot = findobj(H.fig,'tag','show_fbte');
  
  % axes
  pos = [50 235 200 505];
  H.xplt = axes('units','pixel','position',pos,'Parent',H.fig,...
                'DefaultLineHitTest','off','DefaultPatchHitTest','off');
  % Get static data
  mdsopen('static',-1);
  r_v_in = mdsdata('static("\\r_v:in")');
  z_v_in = mdsdata('static("\\z_v:in")');
  r_v_out = mdsdata('static("\\r_v:out")');
  z_v_out = mdsdata('static("\\z_v:out")');
  r_t = mdsdata('static("\\r_t")');
  z_t = mdsdata('static("\\z_t")');
  i_t = find(r_t > 0.88 & z_t < 0,1,'first');
  r_t = r_t([i_t:end,1:i_t-1]);
  z_t = z_t([i_t:end,1:i_t-1]);
  patch([r_v_in;flipud(r_v_out)],[z_v_in;flipud(z_v_out)],...
   [.5 .5 .5],'edgecolor',[.5 .5 .5],'Parent',H.xplt)
  patch([r_t;r_v_in],[z_t;z_v_in],...
   [.7 .7 .7],'edgecolor',[.7 .7 .7],'Parent',H.xplt)
  mdsclose;
  [x,y] = meshgrid(0.5:.1:1.5,-1:.1:1);
  line(x(:),y(:),'marker','.','markersize',1,'linestyle','none','color','w','Parent',H.xplt)
  line(repmat(.624,1,8),[-.6055 -.4325 -.2595 -.0865 .0865 .2595 .4325 .6055],'linestyle','none','marker','>','color','r','Parent',H.xplt)
  line(repmat(.624,1,9),[-.661 -.519 -.346 -.173 0 .173 .346 .519 .661],'linestyle','none','marker','<','color','m','Parent',H.xplt)
  for k = 1:numel(H.bndl), H.bndl(k) = line(NaN,NaN,'Parent',H.xplt); end
  set(H.xplt,'dataaspectratio',[1,1,1],'plotboxaspectratio',[.65,1.64,1],'xlim',[.55 1.2],'ylim',[-.82 .82],'xtick',0.5:.1:1.5,'ytick',-1:.1:1,'box','on')
  for k = 1:numel(H.tplt)
   H.tplt(k) = axes('units','pixel','position',[pos(1)+pos(3)+50 pos(2)+pos(4)-k*pos(4)/numel(H.tplt) 415 pos(4)/numel(H.tplt)],'xticklabel','','yticklabel','','box','on','Parent',H.fig,...
                    'DefaultLineHitTest','off','DefaultPatchHitTest','off');
   H.tlin(k)   = line('XData',NaN,'YData',NaN,'Parent',H.tplt(k),'Linestyle','-' ,'marker','o'   ,'Color','b');
   H.tlin(k+5) = line('XData',NaN,'YData',NaN,'Parent',H.tplt(k),'Linestyle','--','marker','none','Color','w','YLimInclude','off');
  end
  
  % shot and equilibrium selection
  pos = [50 750 200 110]; bcol = [.7 .7 .7];
  pp = uipanel('units','pixels','position',pos,'backgroundcolor',bcol,'title',['TCV MGAMS as ' getenv('USER')]);
  H.msge = uicontrol('Parent',pp,'style','text','position',[5 45 190 40],'max',2,'backgroundcolor',bcol);
  uicontrol('Parent',pp,'style','text','position',[5 25 50 20],'string','Shot #','backgroundcolor',bcol);
  H.shot = uicontrol('Parent',pp,'style','edit','position',[55 25 50 20],'string','-1','backgroundcolor',bcol);
  H.time = uicontrol('Parent',pp,'style','text','position',[105 25 50 20],'string','ms','backgroundcolor',bcol);
  uicontrol('Parent',pp,'style','text','position',[5 5 50 20],'string','Equil #','backgroundcolor',bcol);
  uicontrol('Parent',pp,'style','pushbutton','position',[55 5 25 20],'string','<<','callback','mgams equil previous10','horizontalalignment','center','backgroundcolor',bcol);
  uicontrol('Parent',pp,'style','pushbutton','position',[80 5 25 20],'string','<','callback','mgams equil previous','horizontalalignment','center','backgroundcolor',bcol);
  H.enum = uicontrol('Parent',pp,'style','edit','position',[105 5 40 20],'callback','mgams equil select','backgroundcolor',bcol);
  uicontrol('Parent',pp,'style','pushbutton','position',[145 5 25 20],'string','>','callback','mgams equil next','horizontalalignment','center','backgroundcolor',bcol);
  uicontrol('Parent',pp,'style','pushbutton','position',[170 5 25 20],'string','>>','callback','mgams equil next10','horizontalalignment','center','backgroundcolor',bcol);
  
  % user selected variables
  pos = [300 750 415 110]; col = [.7 .7 0];
  pp = uipanel('units','pixels','position',pos,'backgroundcolor',col);
  for k = 1:NUVAR
   mgamsitem('edit',[5+fix((k-1)/5)*200 85-rem(k-1,5)*20 100],'',col,pp)
  end
  
  % group variables
  pos = [725 700-NGVAR*20 210 NGVAR*20+40]; bcol = [.3 1 .3];
  pp=uipanel('units','pixels','position',pos,'backgroundcolor',bcol,'tag','MGgrp');
  uicontrol('Parent',pp,'style','pushbutton','position',[105 NGVAR*20+5 50 20],'string','Copy','callback','mgams copyall','backgroundcolor',bcol)
  uicontrol('Parent',pp,'style','pushbutton','position',[155 NGVAR*20+5 50 20],'string','Paste','callback','mgams pastegrp','backgroundcolor',bcol)
  for k = 1:NGVAR
   mgamsitem('text',[  5 (NGVAR-k)*20+5 100],'',bcol,pp)
  end
  
  % modal
  pos = [725 750 210 110];
  pp=uipanel('units','pixels','position',pos,'backgroundcolor',mbcol);
  H.dlgb = uicontrol('parent',pp,'style','pushbutton','position',[5 85 200 20],'string','Parametrised Action','backgroundcolor',mbcol,'Tooltip','No action selected');
  for k = 1:numel(H.dlgl)
   H.dlgl(k) = uicontrol('parent',pp,'style','text','position',[5 85-k*20 150 20],'backgroundcolor',mbcol,'Visible','off');
   H.dlge(k) = uicontrol('parent',pp,'style','edit','position',[155 85-k*20 50 20],'backgroundcolor','y','Tooltip','','Visible','off');
  end
  
  % shape
  pos = [5 5 930 185]; bcol = [1 .3 .3];
  pp=uipanel('units','pixels','position',pos,'backgroundcolor',bcol,'Title','Shape','backgroundcolor',bcol);
  H.cpos = uicontrol('parent',pp,'style','text','position',[pos(1:2)+[5 105] 100 40],'backgroundcolor',bcol);
  mgamsitem('text',[105 145  50],'timeeq',bcol,pp);
  mgamsitem('text',[255 145  50],'iansha',bcol,pp);
  mgamsitem('text',[405 145  50],'ilia',bcol,pp);
  mgamsitem('text',[105 125  50],'placu1',bcol,pp);
  mgamsitem('text',[105 105  50],'vsec',bcol,pp);
  mgamsitem('text',[255 125  50],'rmajo1',bcol,pp);
  mgamsitem('text',[405 125  50],'zmajo1',bcol,pp);
  mgamsitem('text',[555 125  50],'rmino1',bcol,pp);
  mgamsitem('text',[255 105  50],'cappa1',bcol,pp);
  mgamsitem('text',[405 105  50],'delta1',bcol,pp);
  mgamsitem('text',[555 105  50],'hlamd1',bcol,pp);
  mgamsitem('text',[  5  85 300],'rbro',bcol,pp);
  mgamsitem('text',[405  85 300],'zbro',bcol,pp);
  mgamsitem('text',[  5  65 300],'rbzo',bcol,pp);
  mgamsitem('text',[405  65 300],'zbzo',bcol,pp);
  mgamsitem('text',[  5  45 300],'rlim1',bcol,pp);
  mgamsitem('text',[405  45 300],'zlim1',bcol,pp);
  mgamsitem('text',[  5  25 820],'rlia1',bcol,pp);
  mgamsitem('text',[  5   5 820],'zlia1',bcol,pp);
  
  mgams var
  mgams('grp',1)
  status = copyfile('/home/tcv_oper/mgams/base_files/fbtmeas.dat','./');
  mgams load
  set(H.fig,'handlevisibility','callback')
  mgamsdisplay;
  return
 end
 
 %% callbacks
 try
  if ~check(MP), fix(MP); end
  assert(check(MP),'Could not fix incoherent dimensions in equilibrium dependent parameters')
  switch mode
   
   %% mouse
   case 'mouse'
    % Check this was a standard left button mouse click
    if ~strcmp(get(H.fig,'SelectionType'),'normal'), return; end
    % Check we clicked an axes object
    obj = get(H.fig,'CurrentObject');
    if ~strcmp(obj.Type,'axes'), return; end
    % Discriminate between time plots and cross-section plot
    if any(obj == H.tplt)
     x = get(obj,'currentpoint'); x = x(1);
     mgams('equil',iround(mgamstimebase,x(1)));
    elseif obj == H.xplt
     x = get(H.xplt,'currentpoint'); x = x(1,[1 2]);
     tmp = sprintf('%+6.3f,%+6.3f',x(1),x(2));
     if any(~isnan(MP.rlia1(:,MGceq)))
      tmp = [tmp ',#' int2str(imin(MP.rlia1(:,MGceq)-x(1)+1i*MP.zlia1(:,MGceq)-1i*x(2)))];
     end
     set(H.cpos,'string',tmp)
    end
    
    %% load [shot]
   case 'load'
    if nargin == 2, set(H.shot,'string',arg1), end
    mgamsmsge(sprintf('Loading shot #%s',get(H.shot,'string')),'r')
    shot = str2double(get(H.shot,'string'));
    assert(isfinite(shot),'MatrixGenerator:GUI','Invalid shot #')
    % Get list of current gas traces from MODEL
    mgc.mdso(-1,'r');
    tmp = mdsdata('GETNCI(["\\PCS::DRAW_REFS:REF_021","\\PCS::DRAW_FEEDFOR_GAS:*"],"PATH")');
    % Load all mgams parameters
    MP.load(shot)
    % Now load gas traces
    MGgas = repmat(struct([]),numel(tmp),1);
    for k = 1:numel(tmp)
     MGgas(k).node = tmp{k};
     MGgas(k).x    = mdsdata(['IF_ERROR(' tmp{k} ',[0,0])']);
     MGgas(k).t    = mdsdata(['IF_ERROR(DIM_OF(' tmp{k} '),[0,1])']);
    end
    FBTE_eq = mgamsloadfbte(shot);
    
    % convert Ip feedback to IOH feedback during breakdown (optionally)
    % use reference breakdown and prefill (optionally)
%    options.Interpreter='tex';options.Default='OK';
%    shn=get(H.shot,'string');
%    shl=numel(shn);
%    % if shl > 3 & shl < 6 & MP.ierat~=4,
%    if shl > 3 && (shl<6 || str2double(shn)>100010),
%     %  selection = questdlg('Convert to IOH feedback for breakdown?', ...
%     %    '','OK','No, thanks',options);
%     selection = questdlg('Use reference breakdown and prefill for this scenario?', ...
%      '','OK','No, thanks',options);
%     if strcmp(selection,'OK'),
%      %    ioh_bd_control_replicate(str2double(get(H.shot,'string')));
%      bd_status=blessed_breakdown;
%      if ~bd_status,
%       disp('No reference breakdown available for this scenario');
%      else
%       set(H.shot,'string','100100');
%       mgams save
%       set(H.shot,'string',shn);
%      end
%     end
%    end
    mgams('equil',MP.numeq);
    
    wgame pan
    % Display actual shot loaded
    mgamsmsge(sprintf('Shot #%d loaded',shot));
    if shot>0&&shot<100000&&numel(MP.zeecorr)>1&&~isnan(MP.zeecorr(2))&&MP.zeecorr(2), ...
	questdlg(['Consider setting zeecorr(2) = ',int2str(shot), ...
		  ' for optimal positioning'], ...
	        'Vertical optimization','Got it','Got it');
    end
    
    %% save [shot]
   case 'save'
    if nargin == 2, set(H.shot,'string',arg1), end
    mgamsmsge(sprintf('Saving in shot #%s',get(H.shot,'string')),'r')
    shot = str2double(get(H.shot,'string'));
    assert(isfinite(shot),'MatrixGenerator:GUI','Invalid shot #')
    MP.save(shot)
    for k = 1:numel(MGgas)
     if mdsdata('WAVE_PUT($1,F_FLOAT($2),F_FLOAT($3))',MGgas(k).node,MGgas(k).x,MGgas(k).t)
      mgamsmsge(['Cannot put wave ' MGgas(k).node],'f'), return
     end
    end
    wgame pan
    mgamsmsge(sprintf('Shot #%s saved',get(H.shot,'string')))
    
    %% prepare
   case 'prepare'
    % Check for discarded FBTE runs (numeq < nruns)
    if MP.numeq < MP.nruns
     answer = questdlg(sprintf('Equilibria %d to %d are currently deselected\nand will not be part of the computed references.\n\n - Choose "Proceed" if this was planned\n - Choose "Fix" to select up to the last eq. (this is the default option)\n - Choose "Cancel" to perform a different selection using the "Set flat top" macro',MP.numeq+1,MP.nruns),'Warning: MGAMS equilibrium selection','Proceed','Fix','Cancel','Fix');
     switch answer
      case 'Cancel', return;
      case 'Fix', MP.numeq = MP.nruns; mgams('save');
      case 'Proceed'
      otherwise, error('Something awful happenened');
     end
    end
    % Protect MODEL shot
    if strcmp(get(H.shot,'string'),'-1') && ~strcmp(getenv('USER'),'tcv_oper')
     mgamsmsge('Only tcv_oper can prepare shot -1','f');return
    end
    mgamsmsge('Preparing shot','c')
    % Restrict search to current MGAMS figure
    h = findobj(H.fig,'checked','on','-and','tag','prepconfig');
    % [AM 04.jul.2017] It seems the creation order for handles has changed
    % We should not rely on this.
% [SC 30.11.18] In Matlab 9.2 at least the order of the object vector itself
%   simply appears to be backwards (whereas the order of the handles is
%   jumbled)
%    if ~verLessThan('matlab','8.4.0')
%      direction = 'descend';
%    else
%      direction = 'ascend';
%    end
    % Sort handles
%    [~,order] = sort(double(h),direction);
%    h = h(order);
    cmd = cellstr(strrep(get(h,'userdata'),'#',get(H.shot,'string')));
    lab = cellstr(get(h,'label'));
    % Create user folder for log files
    [status,~,~] = mkdir(fullfile('/tmp',getenv('USER')));
    assert(status, ['Could not create /tmp/',getenv('USER'),' folder. Aborting ...']);
    for k = numel(cmd):-1:1
     if ~strncmp(lab{k},'Blank',5) || (strcmp(get(H.shot,'string'),'-1') && strcmp(questdlg('Do you want to run the blank recipe?','','Yes','No','Yes'),'Yes'))
      mgamsmsge(['Doing ' lab{k}],'c')
      try
       eval(cmd{k})
       assert(s==0)
       mgamsmsge([lab{k} ' done'])
      catch ME
       % Display error message for MATLAB functions.
       if s==0,disp(getReport(ME));end
       mgamsmsge([lab{k} ' failed'],'f')
       return
      end
     end
    end
    if numel(MP.zeecorr)>1&&~isnan(MP.zeecorr(2))&&MP.zeecorr(2)
     mgamsmsge('Doing vertical optimization','c')
     try
      pos_corr(MP.zeecorr(2),[],[],[],[],1,str2num(get(H.shot,'string')));
     catch ME
      mgamsmsge(ME.message,'f')
      disp(ME.getReport)
      disp('Did Liuqe run for your reference shot?');
      disp('The a-la-carte vertical reference option can be canceled by setting zeecorr(2) to zero');
     end
    end
    mgamsmsge('Shot prepared')
    
    % Re-load data for FBTE plots
    if any(strncmp(lab,'FBTE-',5))
     shot = str2double(get(H.shot,'string'));
     FBTE_eq = mgamsloadfbte(shot);
    end
    
    % Re-open shot
    if any(strcmp(lab,'FBTE18-Matlab'))
     shot = str2double(get(H.shot,'string'));
     MP.mdso(shot,'r');
    end
    % In case of the model tree, fbte18.m opens the results tree to set the
    % default LIUQE parameters. As a consequence, further WGAME calls would
    % fail if shot was not re-loaded first.
    
    %% blank recipe standalone
   case 'blank_recipe'
    if strcmp(get(H.shot,'string'),'-1')
      mgamsmsge('Doing Blank-Matlab','c')
      recstat=tcv_recipe('fbte');
      if recstat
       mgamsmsge('Blank recipe failed','f')
      else
       mgamsmsge('Blank recipe done')
      end
    end
    
    %% restore default prepare configuration
   case 'restore_prep_config'
    Hfig = get(H.xplt,'Parent');
    % Deselect all
    h = findobj(Hfig,'tag','prepconfig');
    set(h,'Checked','off');
    % Select FBTE-Matlab
    h = findobj(Hfig,'tag','prepconfig','-and','Text','FBTE-Matlab');
    set(h,'Checked','on');
    % Select MGA-Matlab compatibility mode
    h = findobj(Hfig,'tag','prepconfig','-and','Text','MGA-Matlab compatibility mode');
    set(h,'Checked','on');
    % Select FBTE18-Matlab
    h = findobj(Hfig,'tag','prepconfig','-and','Text','FBTE18-Matlab');
    set(h,'Checked','on');
    % Select Blank-Matlab
    h = findobj(Hfig,'tag','prepconfig','-and','Text','Blank-Matlab');
    set(h,'Checked','on');

    %% type print
   case {'type','print'}
    mgamsmsge('printing','c')
    shot = str2double(get(H.shot,'string'));
    fid = max(fopen('mgams.lis','w'),1);
    fprintf(fid,'TCV MGAMS #%d %s\n*** Constants\n',shot,date);
    l = 0;
    for k = find(MP.size == 0)
     x = MP.(MP.name{k});
     if numel(x) == 1
      if l == 3; fprintf(fid,'\n'); l = 0; end
      fprintf(fid,'%-s %10.4g  ',MP.name{k},x);
      l = l + 1;
     end
    end
    for k = find(MP.size == 0)
     x = MP.(MP.name{k});
     if numel(x) > 1
      fprintf(fid,'\n%-s',MP.name{k});
      fprintf(fid,' %.4g',x);
     end
    end
    for keq = 1:numel(MP.timeeq)
     fprintf(fid,'\n*** Equilibrium #%d\n',keq);
     l = 0;
     for k = find(MP.size == 1)
      x = MP.(MP.name{k}); x = x(keq);
      if l == 3; fprintf(fid,'\n'); l = 0; end
      fprintf(fid,'%-s %10.4g  ',MP.name{k},x);
      l = l + 1;
     end
     for k = find(MP.size == 2)
      x = MP.(MP.name{k}); x = x(:,keq);
      fprintf(fid,'\n%-s',MP.name{k});
      fprintf(fid,' %.4g',x);
     end
    end
    fclose(fid);
    if strcmp(mode,'print')
     system('lpr mgams.lis');
    else
     type mgams.lis
    end
    mgamsmsge('Print completed')
    
    %% quit
   case 'quit'
    MGceq = [];
    varlist = MP.name(cell2mat(get(H.vare(1:NUVAR),'userdata')));
    save ~/.mgamspref.mat varlist
    close(H.fig)
    wgame quit
    return
    
    %% mgams copy <#eq>
   case 'copy'
    if nargin == 1
     mgamsmodal('Copy','copy','Copy equilibria',...
                'Number of EQ','1','Maximum number of eq. to copy.')
     return
    end
    mgamsmsge('Copying equilibrium','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    if isnan(arg1) || arg1 <= 0, mgamsmsge('Wrong number of EQ','f'), return, end
    k = MGceq:min(MGceq+arg1-1,numel(MP.timeeq));
    for l = find(MP.eqdp)
     MGclip.(MP.name{l}) = MP.(MP.name{l})(:,k);
    end
    MGclip.neq_ = numel(k);
    mgamsmsge(sprintf('%d EQ copied',numel(k)))
    
    %% mgams paste <dt> <paste time> <offset>
   case 'paste'
    if nargin == 1
     mgamsmodal('Paste','paste','Paste equilibria',...
                'Time separation','NaN',['Time difference between successive eq. if "Paste time" is 0\n',...
                                         'If NaN eq. will be equally spaced between the current eq. and the next\n',...
                                         'or 50ms will be used if there is no next eq.\n',...
                                         'An error will be triggered if there is no room to accomodate\n',...
                                         'the number of copied eq. and time separation'],...
                 'Paste time',0,['If true and time vectors are compatible\n',...
                                 'the time values from the copied eq. will be preserved.'],...
                 'Offset','NaN','Set to -1 to paste eq. before the first eq.')
     return
    end
    mgamsmsge('Pasting equilibrium','r')
    if strcmp(arg1,'eval'), [arg1,arg2,arg3] = mgamsmodaleval; end
    if isnan(arg3), arg3 = 0; end
    neq = MGclip.neq_;
    if neq == 0; mgamsmsge('No EQ copied','f'), return, end
    if arg2 && ...
      (MGceq+arg3 > 0                 && MGclip.timeeq(1)   <= MP.timeeq(MGceq+arg3) || ...
      MGceq+arg3 < numel(MP.timeeq) && MGclip.timeeq(neq) >= MP.timeeq(MGceq+arg3+1))
     mgamsmsge('Cannot paste time','f')
     return
    end
    mgams('duplicate',neq,arg1)
    for l = find(MP.eqdp)
     if ~arg2 && strcmp(MP.name{l},'timeeq'), continue; end
     x = MGclip.(MP.name{l}); if size(x,2) ~= neq, x = x'; end
     tmp = MP.(MP.name{l});
     X = [tmp; NaN(size(x,1)-size(tmp,1),size(tmp,2))];
     X(:,MGceq+arg3:MGceq+arg3+neq-1) = [x; NaN(size(X,1)-size(x,1),neq)];
     MP.(MP.name{l}) = X;
    end
    
    mgamsmsge(sprintf('%d EQ pasted',neq))
    
    %% mgams copyrampup
   case 'copyrampup'
    if nargin == 1
     mgamsmodal('Copy ramp up','copyrampup','Copy ramp up equilibria',...
                'Plasma current','160','Copy eq. up to this value of Ip')
     return
    end
    mgamsmsge('Copying ramp up','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    k = min([find(MP.placu1 > arg1) numel(MP.placu1)+1])-1;
    mgams('copyall')
    mgams('equil',1)
    mgams('copy',k)
    mgamsmsge(sprintf('%d ramp up EQ copied',k))
    
    %% mgams pasterampup <dt>
   case 'pasterampup'
    if nargin == 1
     mgamsmodal('Paste ramp up','pasterampup','Replace ramp up equilibria',...
                'Plasma current','160','Replace eq. up to this value of Ip with copied eq.')
     return
    end
    mgamsmsge('Pasting ramp up','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    k = min([find(MP.placu1 > arg1) numel(MP.placu1)+1])-1;
    mgams('pastegrp',find(strcmp('Breakdown',{MGgrp.name})));
    mgams('equil',1)
    mgams('delete',k)
    mgams('paste',0,1,-1)
    mgamsmsge(sprintf('%d/%d ramp up EQ deleted/pasted',k,MGclip.neq_))
    
    %% mgams duplicate <#eq> <dt>
   case 'duplicate'
    if nargin == 1
     mgamsmodal('Duplicate','duplicate','Duplicate current equilibrium',...
                'Number of EQ','1','Number of eq. to add',...
                'Time separation','NaN',['Time difference between successive eq.\n',...
                                         'If NaN eq. will be equally spaced between the current eq. and the next\n',...
                                         'or 50ms will be used if there is no next eq.\n',...
                                         'An error will be triggered if there is no room to accomodate\n',...
                                         'the desired number of eq. and time separation'])
     return
    end
    mgamsmsge('Duplicating equilibrium','r')
    if strcmp(arg1,'eval'), [arg1,arg2] = mgamsmodaleval; end
    if isnan(arg1) || arg1 <= 0, mgamsmsge('Wrong number of EQ','f'), return, end
    neq = numel(MP.timeeq);
    if arg2 < 0 || MGceq < neq && (arg1+1)*arg2 >= MP.timeeq(MGceq+1)-MP.timeeq(MGceq)
     mgamsmsge('Wrong time separation','f'), return
    end
    if isnan(arg2)
     if MGceq == neq, arg2 = 0.05;
     else,            arg2 = (MP.timeeq(MGceq+1)-MP.timeeq(MGceq))/(arg1+1); end
    end
    
    for k = find(MP.eqdp)
     MP.(MP.name{k}) = MP.(MP.name{k})(:,[1:MGceq zeros(1,arg1)+MGceq MGceq+1:neq]);
    end
    if MP.numeq > MGceq, MP.numeq = MP.numeq + arg1; end
    MP.timeeq(MGceq+(1:arg1)) = MP.timeeq(MGceq) + (1:arg1)*arg2;
    MP.nruns = numel(MP.timeeq);
    mgams('equil',MGceq+1);
    
    mgamsmsge(sprintf('EQ duplicated %d times',arg1))
    
    %% rampdown
   case 'rampdown'
    if nargin == 1
     mgamsmodal('Ramp down','rampdown','Convert program to explicit ramp-down',...
                'Flat top','NaN',['Length of flat top for explicit ramp-down.\n',...'
                                  'Uses variable "flattop" if NaN.'],...
                'Time factor','NaN',['Ratio of ramp down time to ramp up time.\n',...
                                     'Uses variable "amagic" if NaN.'])
     return
    end
    mgamsmsge('Ramp down','r')
    if strcmp(arg1,'eval'), [arg1,arg2] = mgamsmodaleval; end
    if isnan(arg1), arg1 = MP.flattop; end
    if isnan(arg2), arg2 = MP.amagic; end
    neq = numel(MP.timeeq);
    if neq < MP.numeq, mgamsmsge('Would erase EQs','f'), return, end
    if arg1 <= 0, mgamsmsge('Wrong flat top','f'), return, end
    if arg2 <= 0, mgamsmsge('Wrong time factor','f'), return, end
    for k = find(MP.eqdp)
     MP.(MP.name{k}) = MP.(MP.name{k})(:,[1:neq neq:-1:1]);
    end
    MP.numeq = 2 * neq;
    MP.nruns = 2 * neq;
    MP.amagic = 0;
    MP.flattop = NaN;
    t = mgamstimeeq2t(MP.timeeq);
    t(neq+1:2*neq) = t(neq) + arg1 + arg2*(t(neq) - t(neq:-1:1));
    MP.timeeq = mgamst2timeeq(t);
    
    mgamsmsge(sprintf('Ramp down %5.3fs %5.3fs',t(neq+1),t(2*neq)))
    
    %% mgams delete <#eq>
   case 'delete'
    if nargin == 1
     mgamsmodal('Delete','delete','Delete equilibria',...
                'Number EQ','1','Number of eq. to delete')
     return
    end
    mgamsmsge('Deleting EQ','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    if isnan(arg1) || arg1 <= 0, mgamsmsge('Wrong number of EQ','f'), return, end
    neq = numel(MP.timeeq);
    if neq == 1, mgamsmsge('Only one EQ left','f'), return, end
    for k = find(MP.eqdp)
     MP.(MP.name{k}) = MP.(MP.name{k})(:,[1:MGceq-1 MGceq+arg1:neq]);
    end
    if MP.numeq > MGceq+arg1-1
     MP.numeq = MP.numeq - arg1;
    elseif MP.numeq >= MGceq && MP.numeq <= MGceq+arg1-1
     MP.numeq = max(MGceq-1,1);
    end
    MP.nruns = numel(MP.timeeq);
    mgams('equil',max(MGceq-1,1));
    
    mgamsmsge(sprintf('%d EQ deleted',arg1))
    
    %% interpolate
   case 'interpolate'
    if nargin == 1
     mgamsmodal('Interpolate','interpolate','Insert new equilibria based on interpolation\nof the current and next one',...
                'Number of EQ','1','Number of eq. to add by interpolation')
     return
    end
    mgamsmsge('Interpolating EQ','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    if isnan(arg1) || arg1 <= 0, mgamsmsge('Wrong number of EQ','f'), return, end
    neq = numel(MP.timeeq);
    if MGceq == neq, mgamsmsge('Impossible for last EQ','f'), return, end
    if ~MP.iansha(MGceq) && ~MP.iansha(MGceq+1) && MP.ilia(MGceq) ~= MP.ilia(MGceq+1)
     mgamsmsge('Incompatible EQ types','f')
     return
    end
    for k = find(MP.eqdp)
     MP.(MP.name{k}) = MP.(MP.name{k})(:,[1:MGceq zeros(1,arg1)+MGceq MGceq+1:neq]);
    end
    k = MGceq:MGceq+arg1+1;
    for l = {'botlim','bzero','capaj','dissi','eftlim','ell1','emm1','omega',...
      'placex','placu1','ppal','ppfac','psirat','qzero','rbro','rbzo',...
      'ritlim','rlim1','strki','timeeq','toplim','ttfac1','vsec','weitam',...
      'weitex','zbro','zbzo','zlim1'}
     MP.(l{1})(k) = linspace(MP.(l{1})(MGceq),MP.(l{1})(k(end)),arg1+2);
    end
    if MP.iansha(MGceq) && MP.iansha(k(end)) % analytic to analytic
     for l = {'rmajo1','zmajo1','rmino1','cappa1','delta1','hlamd1'}
      MP.(l{1})(k) = linspace(MP.(l{1})(MGceq),MP.(l{1})(k(end)),arg1+2);
     end
    else
     if MP.iansha(MGceq+arg1+1)
      nx = MP.ilia(MGceq);
     else
      nx = MP.ilia(MGceq+arg1+1);
     end
     x = 2*pi/nx*(0:nx-1)';
     MP.iansha(MGceq+1:MGceq+arg1) = zeros(1,arg1);
     MP.rlia1(:,MGceq+1:MGceq+arg1) = zeros(size(MP.rlia1,1),arg1) + NaN;
     MP.zlia1(:,MGceq+1:MGceq+arg1) = zeros(size(MP.zlia1,1),arg1) + NaN;
     if MP.iansha(MGceq)
      rz = mgamsanalytic(x,MGceq);
     else
      rz = MP.rlia1(1:nx,MGceq) + 1i*MP.zlia1(1:nx,MGceq);
     end
     if MP.iansha(MGceq+arg1+1)
      rz = [rz, mgamsanalytic(x,MGceq+arg1+1)];
     else
      rz = [rz MP.rlia1(1:nx,MGceq+arg1+1)+1i*MP.zlia1(1:nx,MGceq+arg1+1)];
     end
     x = unwrap(angle(rz - ones(nx,1)*mean(rz))); % sort out points
     [~,k] = sort(x(:,1)); rz(:,1) = rz(k,1);
     [~,k] = sort(x(:,2)); rz(:,2) = rz(k,2);
     x = zeros(1,nx); for k = 1:nx, x(k) = norm(rz(:,1)-rz([k:nx 1:k-1],2)); end
     [~,k] = min(x); rz(:,2) = rz([k:nx 1:k-1],2);
     rz = rz(:,1)*ones(1,arg1) + (rz(:,2)-rz(:,1))*(1:arg1)/(arg1+1);
     MP.rlia1(1:nx,MGceq+1:MGceq+arg1) = real(rz);
     MP.zlia1(1:nx,MGceq+1:MGceq+arg1) = imag(rz);
    end
    if MP.numeq > MGceq, MP.numeq = MP.numeq + arg1; end
    MP.nruns = numel(MP.timeeq);
    
    mgamsmsge(sprintf('%d EQ interpolated',arg1))
    
    %% discretise
   case 'discretise'
    if nargin == 1
     mgamsmodal('Discretise','discretise','Convert analytic boundary to discretised form',...
                'Number of points','12','Number of (rlia1,zlia1) points to use for discretised boundary')
     return
    end
    mgamsmsge('Discretising EQ','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    if isnan(arg1) || arg1 <= 0, mgamsmsge('Wrong number of points','f'), return, end
    x = mgamsanalytic(2*pi/arg1*(0:arg1-1)',MGceq);
    MP.rlia1 = [MP.rlia1; zeros(numel(x)-size(MP.rlia1,1),size(MP.rlia1,2))+NaN];
    MP.zlia1 = [MP.zlia1; zeros(numel(x)-size(MP.zlia1,1),size(MP.zlia1,2))+NaN];
    x = [x; zeros(size(MP.rlia1,1)-numel(x),1)+NaN];
% temporarily replaced following line to fix a problem with methods
%   MP.rlia1(:,MGceq) = real(x); MP.zlia1(:,MGceq) = imag(x);
    MP.rlia1(1:numel(x),MGceq) = real(x); MP.zlia1(1:numel(x),MGceq) = imag(x);
    MP.iansha(MGceq) = 0; MP.ilia(MGceq) = arg1;
    
    mgamsmsge(sprintf('EQ discretised in %d points',arg1))
    
    %% shift
   case 'shift'
    if nargin == 1
     mgamsmodal('Shift','shift','Shift equilibria in the poloidal plane',...
                'Number of EQ','1','Number of eq. to modify',...
                'Delta-R','0','Distance to shift in the R direction',...
                'Delta-Z','0','Distance to shift in the Z direction')
     return
    end
    mgamsmsge('Shifting EQ','r')
    if strcmp(arg1,'eval'), [arg1,arg2,arg3] = mgamsmodaleval; end
    if isnan(arg1) || isnan(arg2) || isnan(arg3) || arg1 < 0, mgamsmsge('Wrong arguments','f'), return, end
    k = MGceq:min(MGceq+arg1-1,numel(MP.timeeq));
    MP.rmajo1(k)  = MP.rmajo1(k)  + arg2; MP.zmajo1(k)  = MP.zmajo1(k)  + arg3;
    MP.rbro(:,k)  = MP.rbro(:,k)  + arg2; MP.zbro(:,k)  = MP.zbro(:,k)  + arg3;
    MP.rbzo(:,k)  = MP.rbzo(:,k)  + arg2; MP.zbzo(:,k)  = MP.zbzo(:,k)  + arg3;
    MP.rlim1(:,k) = MP.rlim1(:,k) + arg2; MP.zlim1(:,k) = MP.zlim1(:,k) + arg3;
    MP.rlia1(:,k) = MP.rlia1(:,k) + arg2; MP.zlia1(:,k) = MP.zlia1(:,k) + arg3;
    
    mgamsmsge(sprintf('%d EQ shifted by (%.4g,%.4g)',arg1,arg2,arg3))
    
    %% scale
   case 'scale'
    if nargin == 1
     mgamsmodal('Scale','scale','Scale equilibria in the poloidal plane',...
                'Number EQ','1','Number of eq. to modify',...
                'R-axis','NaN',['R coordinate of the center of transformation.\n',...
                                'If NaN R coordinates will be left untouched.'],...
                'Z-axis','NaN',['Z coordinate of the center of transformation.\n',...
                                'If NaN Z coordinates will be left untouched.'],...
                'Factor','1','Scaling factor for the transformation.')
     return
    end
    mgamsmsge('Scaling EQ','r')
    if strcmp(arg1,'eval'), [arg1,arg2,arg3,arg4] = mgamsmodaleval; end
    if isnan(arg1) || arg1 < 0, mgamsmsge('Wrong number of EQ','f'), return, end
    k = MGceq:min(MGceq+arg1-1,numel(MP.timeeq));
    if ~isnan(arg2)
     MP.rbro(:,k)  = arg2 + (MP.rbro(:,k)  - arg2) * arg4;
     MP.rbzo(:,k)  = arg2 + (MP.rbzo(:,k)  - arg2) * arg4;
     MP.rlim1(:,k) = arg2 + (MP.rlim1(:,k) - arg2) * arg4;
     MP.rlia1(:,k) = arg2 + (MP.rlia1(:,k) - arg2) * arg4;
    end
    if ~isnan(arg3)
     MP.zbro(:,k)  = arg3 + (MP.zbro(:,k)  - arg3) * arg4;
     MP.zbzo(:,k)  = arg3 + (MP.zbzo(:,k)  - arg3) * arg4;
     MP.zlim1(:,k) = arg3 + (MP.zlim1(:,k) - arg3) * arg4;
     MP.zlia1(:,k) = arg3 + (MP.zlia1(:,k) - arg3) * arg4;
    end
    if ~isnan(arg2) && ~isnan(arg3)
     MP.rmajo1(k) = arg2 + (MP.rmajo1(k) - arg2) * arg4;
     MP.zmajo1(k) = arg3 + (MP.zmajo1(k) - arg3) * arg4;
     MP.rmino1(k) = MP.rmino1(k) * abs(arg4);
     MP.delta1(k) = MP.delta1(k) * sign(arg4);
    elseif ~isnan(arg2)
     MP.rmajo1(k) = arg2 + (MP.rmajo1(k) - arg2) * arg4;
     MP.rmino1(k) = MP.rmino1(k) * abs(arg4);
     MP.cappa1(k) = MP.cappa1(k) / abs(arg4);
     MP.delta1(k) = MP.delta1(k) * sign(arg4);
    elseif ~isnan(arg3)
     MP.zmajo1(k) = arg3 + (MP.zmajo1(k) - arg3) * arg4;
     MP.cappa1(k) = MP.cappa1(k) * abs(arg4);
    end
    
    mgamsmsge(sprintf('%d EQ scaled around (%.4g,%.4g) by %.4g',arg1,arg2,arg3,arg4))
    
    %% mirror
   case 'mirror'
    if nargin == 1
     mgamsmodal('Mirror','mirror','Mirror equilibria in the poloidal plane',...
                'Number EQ','1','Number of eq. to modify',...
                'R-axis','NaN',['R coordinate of the center of transformation.\n',...
                                'If NaN R coordinates will be left untouched.'],...
                'Z-axis','NaN',['Z coordinate of the center of transformation.\n',...
                                'If NaN Z coordinates will be left untouched.'])
     return
    end
    mgamsmsge('Mirroring EQ','r')
    if strcmp(arg1,'eval'), [arg1,arg2,arg3] = mgamsmodaleval; end
    mgams('scale',arg1,arg2,arg3,-1)
    mgamsmsge(sprintf('%d EQ mirrored around (%.4g,%.4g)',arg1,arg2,arg3))
    
    %% halfanalytic
   case 'halfanalytic'
    if nargin == 1
     mgamsmodal('Half Analytic','halfanalytic','Set control points using analytic formula above or below a certain Z value',...
                'Above Z','NaN',['Only points above this value will be modified.\n',...
                                 'Leave NaN to use "Below Z" option.'],...
                'Below Z','NaN','Only points below this value will be modified.')
     return
    end
    mgamsmsge('Computing half analytic EQ','r')
    if strcmp(arg1,'eval'), [arg1,arg2] = mgamsmodaleval; end
    if MP.iansha(MGceq)
     mgamsmsge('EQ is analytic','f')
     return
    end
    if ~isnan(arg1)
     k = find(MP.zlia1(:,MGceq) >= arg1 & ~isnan(MP.zlia1(:,MGceq)));
     h = 'top';
    elseif ~isnan(arg2)
     k = find(MP.zlia1(:,MGceq) <= arg2 & ~isnan(MP.zlia1(:,MGceq)));
     h = 'bottom';
    else
     k = [];
     h = 'no';
    end
    if ~isempty(k)
     dR = MP.rlia1(k,MGceq)-MP.rmajo1(MGceq);
     dZ = MP.zlia1(k,MGceq)-MP.zmajo1(MGceq);
     % Try to conserve poloidal distribution of points
     thetag = unwrap(atan2(dZ,dR));
     % First get geometrical angle as function of x
     nt = 65;
     switch h, case 'top',    t = linspace(0  ,pi,nt);
               case 'bottom', t = linspace(-pi,0 ,nt);
     end
     y = mgamsanalytic(t,MGceq);
     dR_t = real(y)-MP.rmajo1(MGceq);
     dZ_t = imag(y)-MP.zmajo1(MGceq);
     thetag_t = unwrap(atan2(dZ_t,dR_t));
     % Now interpolate to get x
     t = interp1(thetag_t,t,thetag,'spline');
     x = mgamsanalytic(t,MGceq);
    else, x = [];
    end
    MP.rlia1(k,MGceq) = real(x); MP.zlia1(k,MGceq) = imag(x);
    
    mgamsmsge(sprintf('%s-half analytic EQ',h))
    
    %% analyticfit
   case 'analyticfit'
    if nargin == 1
     mgamsmodal('Analytic fit','analyticfit','Fit control points above or below a certain Z value using analytic formula',...
                'Above Z','NaN',['Only points above this value will be fitted.\n',...
                                 'Leave NaN to use "Below Z" option.'],...
                'Below Z','NaN','Only points below this value will be fitted.')
     return
    end
    mgamsmsge('Computing analytic fit for EQ','r')
    if strcmp(arg1,'eval'), [arg1,arg2] = mgamsmodaleval; end
    if MP.iansha(MGceq)
     mgamsmsge('EQ is analytic','f')
     return
    end
    if ~isnan(arg1)
     k = find(MP.zlia1(:,MGceq) >= arg1 & ~isnan(MP.zlia1(:,MGceq)));
     h = 'top-half';
    elseif ~isnan(arg2)
     k = find(MP.zlia1(:,MGceq) <= arg2 & ~isnan(MP.zlia1(:,MGceq)));
     h = 'bottom-half';
    else
     k = ~isnan(MP.zlia1(:,MGceq));
     h = 'full';
    end
    R = MP.rlia1(k,MGceq);
    Z = MP.zlia1(k,MGceq);
    
    x0 = [MP.rmajo1(MGceq);MP.zmajo1(MGceq);MP.rmino1(MGceq);MP.cappa1(MGceq);MP.delta1(MGceq);MP.hlamd1(MGceq);...
          unwrap(atan2(Z-MP.zmajo1(MGceq),R-MP.rmajo1(MGceq)))];
    
    lsqOpts = optimoptions('lsqnonlin', 'Display', 'off', 'FunctionTolerance', 0, 'SpecifyObjectiveGradient', true);
    [x,~,~,exitflag] = lsqnonlin(@(x)analyticfit(R,Z,x),x0,[],[],lsqOpts);
    
    if exitflag <= 0
     mgamsmsge('analytic fit failed','f')
     return
    end
    
    MP.rmajo1(MGceq) = x(1);
    MP.zmajo1(MGceq) = x(2);
    MP.rmino1(MGceq) = x(3);
    MP.cappa1(MGceq) = x(4);
    MP.delta1(MGceq) = x(5);
    MP.hlamd1(MGceq) = x(6);
    
    mgamsmsge(sprintf('%s analytic fit',h))
    
    %% setpoint
   case 'setpoint'
    if nargin == 1
     mgamsmodal('Get/set point','setpoint','Get or set control point coordinate',...
                'Point','NaN','Index of (rlia1,zlia1) point to get/set',...
                'R','NaN',['New value of R coordinate for selected point.\n',...
                           'NaN to get current coordinates'],...
                'Z','NaN',['New value of Z coordinate for selected point.\n',...
                           'NaN to get current coordinates'])
     return
    end
    mgamsmsge('Getting/setting point','r')
    if strcmp(arg1,'eval'), [arg1,arg2,arg3] = mgamsmodaleval; end
    if MP.iansha(MGceq)
     mgamsmsge('EQ is analytic','f')
     return
    end
    if ~isnan(arg1) && arg1 > 0 && arg1 <= size(MP.rlia1,1)
     x = [MP.rlia1(arg1,MGceq) MP.zlia1(arg1,MGceq)];
     if isnan(arg2) || isnan(arg3)
      set(H.dlge(2),'string',mgamsformat(x(1)));
      set(H.dlge(3),'string',mgamsformat(x(2)));
      return
     end
     MP.rlia1(arg1,MGceq) = arg2;
     MP.zlia1(arg1,MGceq) = arg3;
     
     mgamsmsge(sprintf('Point moved to (%s,%s)',mgamsformat(arg2),mgamsformat(arg3)))
    else
     mgamsmsge('Wrong point number','f')
    end
    
    %% currentramp
   case 'currentramp'
    if nargin == 1
     mgamsmodal('Plasma Current Ramp','currentramp','Introduce plasma current ramp between desired equilibria',...
                '1st EQ','1','Index of initial eq. for the ramp',...
                'Last EQ','Inf','Index of the final eq. for the ramp. Inf uses the last eq..',...
                '1st current','NaN','Initial current value. NaN uses the value at the initial eq.',...
                'Last current','NaN','Final current value. NaN uses the value at the final eq.')
     return
    end
    mgamsmsge('Computing plasma current ramp','r')
    if strcmp(arg1,'eval'), [arg1,arg2,arg3,arg4] = mgamsmodaleval; end
    if isnan(arg1) || isnan(arg2), mgamsmsge('Wrong EQ number','f'), return, end
    arg1 = max(arg1,1); arg2 = min(arg2,numel(MP.timeeq));
    if isnan(arg3), arg3 = MP.placu1(arg1); end
    if isnan(arg4), arg4 = MP.placu1(arg2); end
    t = mgamstimebase;
    MP.placu1(arg1:arg2) = arg3 + (t(arg1:arg2)-t(arg1))*(arg4-arg3)/(t(arg2)-t(arg1));
    
    mgamsmsge(sprintf('Plasma current ramp %d/%d',arg3,arg4))
    
    %% setabsolutetime
   case 'setabsolutetime'
    if nargin == 1
     mgamsmodal('Set absolute time','setabsolutetime','Change equilibrium absolute time',...
                'Time',sprintf('%.3f',mgamstimeeq2t(MP.timeeq(MGceq))),'Absolute time value')
     return
    end
    mgamsmsge('Setting absolute time','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    t = mgamst2timeeq(arg1);
    if isnan(t) || (MGceq > 1 && t <= MP.timeeq(MGceq-1)) ...
      || (MGceq < numel(MP.timeeq) && t >= MP.timeeq(MGceq+1))
     mgamsmsge('Wrong time','f'), return
    end
    MP.timeeq(MGceq) = t;
    
    mgamsmsge(sprintf('Absolute time set to %.3f',arg1))
    
    %% setflattop
   case 'setflattop'
    MP.numeq = MGceq;
    mgams('equil',MGceq);
    
    mgamsmsge(sprintf('End of evolution set at #%d',MP.numeq))
    
    %% gascopy
   case 'gascopy'
    mgamsmsge('Copying gas/density','r')
    CBgas = MGgas;
    mgamsmsge('Gas/density copied')
    
    %% gaspaste
   case 'gaspaste'
    mgamsmsge('Pasting gas/density','r')
    if isempty(CBgas), mgamsmsge('Gas/density clipboard empty','f'), return, end
    MGgas = CBgas;
    mgamsmsge('Gas/density pasted')
    
    %% gascopy_bd
   case 'gascopy_bd'
    % Copying all valves (assuming 1st trace is density)
    mgamsmsge('Copying gas BD','r')
    CBgas_bd = MGgas;
    for ii = 2:numel(MGgas)
      mask = MGgas(ii).t <= 0.01;
      CBgas_bd(ii).t = MGgas(ii).t(mask);
      CBgas_bd(ii).x = MGgas(ii).x(mask);
      if ~isempty(CBgas_bd(ii).t) && CBgas_bd(ii).t(end) < 0
        CBgas_bd(ii).t(end+1,1) = 0;
        CBgas_bd(ii).x(end+1,1) = interp1(MGgas(ii).t,MGgas(ii).x,0,'linear',0); % 0 if extrapolate (should never happen for the main valve)
      end
    end
    mgamsmsge('Gas BD copied')
    
    %% gaspaste_bd
   case 'gaspaste_bd'
    % Pasting all valves (assuming 1st trace is density)
    mgamsmsge('Pasting gas BD','r')
    if isempty(CBgas_bd), mgamsmsge('Gas BD clipboard empty','f'), return, end
    for ii = 2:numel(MGgas)
      mask = MGgas(ii).t > 0.01;
      MGgas(ii).t = [CBgas_bd(ii).t;MGgas(ii).t(mask)];
      MGgas(ii).x = [CBgas_bd(ii).x;MGgas(ii).x(mask)];
    end
    mgamsmsge('Gas BD pasted')
    % NB: wgame will not be updated until the trace is saved in the model.
    
    %% gasip
   case 'gasip'
    if nargin == 1
     mgamsmodal('Density Proportional to Plasma Current','gasip','Set density proportional to plasma current',...
                'Flat top density','1','Density value in fringes for the last point of the program.')
     return
    end
    mgamsmsge('Setting ne ~ Ip','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    x = [0 MP.placu1(1:MP.numeq)/MP.placu1(MP.numeq)*arg1]; x = [x fliplr(x)];
    t = mgamsallt; t = [0 t t(end)+0.01];
    MGgas(1).x = x;
    MGgas(1).t = t;
    mgamsmsge(sprintf('ne ~ Ip set to %4.2f',arg1))
    
    %% usefps
   case 'usefps'
    if nargin == 1
     mgamsmodal('Use FPS','usefps','Enable FPS use',...
                'fastm','0.4','Fast-coil gain in M-matrix')
     return
    end
    mgamsmsge('Setting FPS use','r')
    if strcmp(arg1,'eval'), arg1 = mgamsmodaleval; end
    MP.nfast = -3; MP.rampt = 0.5; MP.ikriz = 2; MP.fastm = arg1;
    
    mgamsmsge(sprintf('Using FPS with fastm = %.3f',arg1))
    
    %% nofps
   case 'nofps'
    if nargin == 1
     mgamsmodal('Do not use FPS','nofps','Disable FPS use')
     return
    end
    mgamsmsge('Setting FPS use','r')
    MP.nfast = 0; MP.rampt = NaN; MP.ikriz = 0; MP.fastm = 0;
    
    mgamsmsge('FPS off')
    
    %% activate_dne
   case 'activate_dne'
    if nargin == 1
     if (isempty(MP.dneon)     || isnan(MP.dneon)),     dneon  = '0.07'; else, dneon  = num2str(MP.dneon);     end
     if (isempty(MP.dnekp)     || isnan(MP.dnekp)),     dnekp  = '5';    else, dnekp  = num2str(MP.dnekp);     end
     if (isempty(MP.dneki)     || isnan(MP.dneki)),     dneki  = '60';   else, dneki  = num2str(MP.dneki);     end
     if (isempty(MP.dnefil(1)) || isnan(MP.dnefil(1))), dnefil = '0.2';  else, dnefil = num2str(MP.dnefil(1)); end
     mgamsmodal('Activate SCD ne control','activate_dne','Set Wouter''s SCD density controller settings',...
                'Activation time',dneon,'Absolute time for activation of SCD gas controller',...
                'Proportional gain',dnekp,'',...
                'Integral gain',dneki,'',...
                'Filter activation time',dnefil,'Absolute time for reset of integral errors')
     return
    end
    mgamsmsge('Setting SCD ne control parameters','r')
    if strcmp(arg1,'eval'), [arg1,arg2,arg3,arg4] = mgamsmodaleval; end
    MP.dneon = arg1; MP.dnekp = arg2; MP.dneki = arg3; MP.dnefil(1) = arg4;
    % Sanity checking on parameters:
    scd_param.dneon.default    = 0.07; scd_param.dneon.limits    = [0 3];   scd_param.dneon.discrete = 0;
    scd_param.dnekp.default    = 5;    scd_param.dnekp.limits    = [0 100]; scd_param.dnekp.discrete = 0;
    scd_param.dneki.default    = 60;   scd_param.dneki.limits    = [0 300]; scd_param.dneki.discrete = 0;
    scd_param.dnekd.default    = 0;    scd_param.dnekd.limits    = [0 100]; scd_param.dnekd.discrete = 0;
    scd_param.dneka.default    = 10;   scd_param.dneka.limits    = [0 20];  scd_param.dneka.discrete = 0;
    scd_param.dnemode.default  = 0;    scd_param.dnemode.limits  = [0 1];   scd_param.dnemode.discrete = 1;
    scd_param.dnescal1.default = 1;    scd_param.dnescal1.limits = [0 3];   scd_param.dnescal1.discrete = 0;
    scd_param.dnescal2.default = 1;    scd_param.dnescal2.limits = [0 3];   scd_param.dnescal2.discrete = 0;
    % This node has been deleted from the MGAMS section of the PCS tree
    % scd_param.dne1mm.default   = 0;    scd_param.dne1mm.limits   = [0 2];   scd_param.dne1mm.discrete = 1;
    scd_param.dnefil.default   = [0.2 25 3]; scd_param.dnefil.limits = [0 3; 10 100; 1 4]; scd_param.dnefil.discrete = [0 0 1];
    params = fieldnames(scd_param);
    for k=1:numel(params)
     for n=1:numel(scd_param.(params{k}).default)
      % Check if parameter is set correctly:
      if (n > numel(MP.(params{k})) || isempty(MP.(params{k})(n,1)) || isnan(MP.(params{k})(n,1)))
       MP.(params{k})(n,1) = scd_param.(params{k}).default(n); % Set to default value
      else
       % If so, check if set values are sane:
       if (scd_param.(params{k}).discrete(n) == 1)
        if(~any(MP.(params{k})(n) == scd_param.(params{k}).limits(n,1):scd_param.(params{k}).limits(n,2)))
         mgamsmsge(['Parameter ' params{k} ' set to invalid value'],'r')
        end
       else
        % Clamp to specified limits:
        if ( MP.(params{k})(n) < scd_param.(params{k}).limits(n,1))
         mgamsmsge(['WARNING: clamped ' params{k} ' to limit value'])
         MP.(params{k})(n,1) = scd_param.(params{k}).limits(n,1);
        elseif (MP.(params{k})(n,1) > scd_param.(params{k}).limits(n,2))
         mgamsmsge(['WARNING: clamped ' params{k} ' to limit value'])
         MP.(params{k})(n) = scd_param.(params{k}).limits(n,2);
        end
       end
      end
     end
    end
    mgamsmsge('SCD ne control set, run proper recipe')
    
    %% equil
   case 'equil' % mgams equil 'select'|'previous10'|'previous'|'next'|'next10'|'eq#'|eq#
    switch arg1
     case 'select',     MGceq = str2double(get(H.enum,'string'));
     case 'previous10', MGceq = MGceq - 10;
     case 'previous',   MGceq = MGceq - 1;
     case 'next',       MGceq = MGceq + 1;
     case 'next10',     MGceq = MGceq + 10;
     otherwise,         if ischar(arg1), MGceq = str2double(arg1); else, MGceq = arg1; end
    end
    if isnan(MGceq), MGceq = 1; end
    MGceq = max(1,min(MGceq,numel(MP.timeeq)));
    set(H.enum,'string',int2str(MGceq))
    set(H.time,'string',sprintf('%6.3fs',mgamstimeeq2t(MP.timeeq(MGceq))))
    
    %% modify
   case 'modify'
    mgamsmsge('Modifying variable','r')
    % Get callback object
    ind = get(gcbo,'userdata');
    par = MP.name(ind);
    x = get(gcbo,'string');
    if ~strncmp(par,'debug_p',7) % Not debug nodes
     try
      x = eval(['[' x ']']);
     catch
      error('MatrixGenerator:GUI','Wrong syntax for "%s"',par)
     end
    end
    MP.set(x,ind,MGceq);
    mgamsmsge(sprintf('"%s" is now %s',par,mgamsformat(x)))
    
    %% mgams('var',pos#[,var#])
   case 'var'
    if nargin == 1
     if exist('~/.mgamspref.mat','file')
      load ~/.mgamspref.mat varlist
     else
      varlist = {'placu1','zmajo1','bzero','qzero','ppfac',...
                 'delta1','testa','iscramb','zeecorr','dexpcode'};
     end
     varlist(NUVAR+1:end) = [];
     mgams('var',1:numel(varlist),MP.ind(varlist))
     mgamsmsge('Pref variable list')
     return
    elseif nargin == 2
     arg2 = get(gcbo,'string');
     arg2 = MP.ind(arg2);
     if isempty(arg2)
      mgamsmsge(['Unknown variable "' get(gcbo,'string') '"'],'f')
      set(gcbo,'string','','userdata',[])
     end
    end
    col = 'kb';
    for k = 1:numel(arg2)
     set([H.varl(arg1(k)) H.vare(arg1(k))],'userdata',arg2(k),'string',MP.name(arg2(k)),'foregroundcolor',col(MP.eqdp(arg2(k))+1),'visible','on')
    end
    
    %% mgams grp <grp#>
   case 'grp'
    mgamsmsge('Changing group','r')
    set(findobj('tag','MGgrp'),'Title',MGgrp(arg1).name,'userdata',arg1)
    set([H.varl(NUVAR+(1:NGVAR)) H.vare(NUVAR+(1:NGVAR))],'visible','off')
    mgams('var',NUVAR+(1:NGVAR),MP.ind(MGgrp(arg1).var))
    mgamsmsge('Group changed')
    
    %% mgams copyall
   case 'copyall'
    mgamsmsge('Copying all MP','r')
    for k = 1:MP.size
     MGclip.(MP.name{k}) = MP.(MP.name{k});
    end
    mgamsmsge('All MP copied')
    
    %% mgams pastegrp <grp#>
   case 'pastegrp'
    if nargin == 1, arg1 = get(findobj('tag','MGgrp'),'userdata'); end
    mgamsmsge(['Pasting group "' MGgrp(arg1).name '"'],'r')
    for k = 1:numel(MGgrp(arg1).var)
     MP.(MGgrp(arg1).var{k}) = MGclip.(MGgrp(arg1).var{k});
    end
    mgamsmsge(['Group "' MGgrp(arg1).name '" pasted'])
    
   otherwise, error('Invalid direct mgams call')
  end
 catch MGerr
  mgamsmsge(MGerr.message,'f')
  disp(MGerr.getReport)
 end
 fix(MP)
 % Update plots only for the top-level call
 stack = dbstack(1);
 recurs = ~isempty(stack) && any(strcmp({stack.name},'mgams'));
 if ~recurs, mgamsdisplay;end
 
end

%% update figure:
function mgamsdisplay
 
 global H MP MGceq FBTE_eq
 
 if isempty(MGceq), return, end

 % Shape panel
 x = linspace(0,2*pi,51);
 if MP.iansha(MGceq)
  rz = mgamsanalytic(x,MGceq);
  set(H.bndl(5)  ,'xdata',real(rz),'ydata',imag(rz),'linestyle','-','marker','none','color','y')
  set(H.bndl(6:8),'xdata',NaN     ,'ydata',NaN)
 else
  set(H.bndl(5),'xdata',MP.rlim1(:,MGceq),'ydata',MP.zlim1(:,MGceq),'marker','+','color','y','markersize', 8,'linestyle','none')
  set(H.bndl(6),'xdata',MP.rlia1(:,MGceq),'ydata',MP.zlia1(:,MGceq),'marker','o','color','y','markersize', 8,'linestyle','none')
  set(H.bndl(7),'xdata',MP.rbro (:,MGceq),'ydata',MP.zbro (:,MGceq),'marker','x','color','y','markersize',12,'linestyle','none')
  set(H.bndl(8),'xdata',MP.rbzo (:,MGceq),'ydata',MP.zbzo (:,MGceq),'marker','+','color','y','markersize',12,'linestyle','none')
 end
 h = H.fbtplot;
 if strcmp(h.Checked,'off')
  % Display next eq. control points
  if MGceq == numel(MP.timeeq)
   % There is no next eq.
   set(H.bndl(1:4),'xdata',NaN,'ydata',NaN)
  else
   if MP.iansha(MGceq+1)
    rz = mgamsanalytic(x,MGceq+1);
    set(H.bndl(1),'xdata',real(rz),'ydata',imag(rz),'linestyle',':','color','g','marker','none')
    set(H.bndl(2:4),'xdata',NaN,'ydata',NaN)
   else
    set(H.bndl(1),'xdata',MP.rlim1(:,MGceq+1),'ydata',MP.zlim1(:,MGceq+1),'marker','+','color','g','markersize',4,'linestyle','none','Linewidth',1)
    set(H.bndl(2),'xdata',MP.rlia1(:,MGceq+1),'ydata',MP.zlia1(:,MGceq+1),'marker','o','color','g','markersize',4,'linestyle','none','Linewidth',1)
    set(H.bndl(3),'xdata',MP.rbro (:,MGceq+1),'ydata',MP.zbro (:,MGceq+1),'marker','x','color','g','markersize',6,'linestyle','none','Linewidth',1)
    set(H.bndl(4),'xdata',MP.rbzo (:,MGceq+1),'ydata',MP.zbzo (:,MGceq+1),'marker','+','color','g','markersize',6,'linestyle','none','Linewidth',1)
   end
  end
 else
  % Display psi-contours from FBTE results
  if ~isempty(FBTE_eq)
   neq = numel(FBTE_eq.r_axis);
   if (MGceq <= neq)
    set(H.bndl(1),'XData', FBTE_eq.r_axis(:,MGceq) ,'YData', FBTE_eq.z_axis(:,MGceq) ,'marker','+'   ,'color',[0.4,0.4,1],'markersize',6,'linestyle','none','Linewidth',1);
    set(H.bndl(2),'XData', FBTE_eq.r_xpts(:,MGceq) ,'YData', FBTE_eq.z_xpts(:,MGceq) ,'marker','x'   ,'color',[0.4,0.4,1],'markersize',6,'linestyle','none','Linewidth',1);
    set(H.bndl(3),'XData', FBTE_eq.Cr{2,MGceq}     ,'YData', FBTE_eq.Cz{2,MGceq}     ,'marker','none','color',[0.4,0.4,1],'markersize',6,'linestyle','-'   ,'Linewidth',2);
    set(H.bndl(4),'XData',[FBTE_eq.Cr{1,MGceq},NaN,FBTE_eq.Cr{3,MGceq}],...
                  'YData',[FBTE_eq.Cz{1,MGceq},NaN,FBTE_eq.Cz{3,MGceq}]              ,'marker','none','color',[0.4,0.4,1],'markersize',6,'linestyle','--'  ,'Linewidth',1);
   else
    % Do not display anything if there's not enough FBTE eq. (e.g. when new eq. have been inserted)
    set(H.bndl(1:4),'XData',NaN,'YData',NaN);
   end
  end
 end

 % Time panel
 t = mgamstimebase;
 for k = 1:numel(H.tplt)
  kpar = get(H.varl(k),'userdata');
  % For used panels, change data and YTickLabelMode
  if isempty(kpar) || kpar == 0 || ~MP.eqdp(kpar)
   set(H.tlin(k)  ,'XData',NaN,'YData',NaN);
   set(H.tlin(k+5),'XData',NaN,'YData',NaN);
   set(H.tplt(k)  ,'YTickLabel','');
  else
   y = MP.get(kpar);
   set(H.tlin(k)  ,'XData',t               ,'YData',y);
   yscale = 10.^ceil(log10(max(abs(y))));
   set(H.tlin(k+5),'XData',t([MGceq MGceq]),'YData',yscale*[-1,1]);
   set(H.tplt(k)  ,'YTickLabelMode','Auto');
  end
 end
 set(H.tplt,'xlim',[min(t)-.01 max(t)+.01],'xtick',unique(sort(t)),'xticklabel','','box','on');
 ax = H.tplt(end);
 
 % Check if labels changed before updating
 labels = arrayfun(@(x) sprintf('%4.0f ',x),t*1000,'UniformOutput',false);
 if ~isequal(labels.',get(H.tlab,'String'))
   delete(H.tlab);
   % Behavior of the text function seems to have changed for versions > 8.3.0
   H.tlab = text(t,repmat(min(get(ax,'ylim')),1,numel(t)),labels,...
     'rotation',90,'horizontalalignment','right','fontsize',10,'color','c','Parent',ax);
 end
 % Set colors for text labels
 set(H.tlab( logical(MP.iansha))     ,'color','m'); % analytic eq.
 set(H.tlab(~logical(MP.iansha))     ,'color','c'); % manual eq.
 set(H.tlab(MP.numeq+1:numel(H.tlab)),'color','w'); % inactive eq.
   
 % Replace all strings in textboxes in one call
 h = H.vare(H.vare ~= 0);
 k = get(h,'userdata');
 k(cellfun(@isempty,k)) = {0};
 set(h,{'string'},(MP.format([k{:}],MGceq)).');
 
end

%% format variable: str = mgamsformat(nvar,neq) or str = mgamsformat(x)
function str = mgamsformat(nvar,neq)
 global MP
 str = '';
 h = '%.4g';
 if nargin == 1
  if ischar(nvar), str=['{',nvar,'}']; return, end % Debug nodes
  x = nvar(:);
 else
  if isempty(nvar) || nvar == 0, return, end
  x = MP.(MP.name{nvar}); if MP.eqdp(nvar), x = x(:,neq); end
  if any(strcmp(MP.name{nvar},{'rmajo1' 'zmajo1' 'rmino1' 'rlim1' 'zlim1' 'rlia1' 'zlia1' 'rbro' 'zbro' 'rbzo' 'zbzo'}))
   h = '%+6.3f';
  end
 end
 x = x(~isnan(x));
 if numel(x) > 1
  str = sprintf([h ','],x); str = ['[' str(1:end-1) ']'];
 elseif numel(x) == 1
  str = sprintf(h,x);
 end
end

% open shot file: status = mgamsopenshot(shot,'r'|'w')
% function stat = mgamsopenshot(shot,mode)
%
% global H MGusers
%
% stat = 0;
% if ischar(shot), shot = str2double(shot); end
% if isnan(shot)
%  mgamsmsge('Wrong shot number','f')
%  return
% end
% for k = 1:numel(MGusers)
%  if mode == 'w' && shot >= MGusers(k).shot1 && shot <= MGusers(k).shot2 ...
%     && ~strcmp(getenv('USER'),MGusers(k).user)
%   mgamsmsge(sprintf('Shot #%d protected',shot),'f')
%   return
%  end
% end
% if shot >= 100000 && mode == 'w' && (isempty(mdsopen('PCS',shot)) || ...
%  uint64(mdsdata('MAXVAL(USING(GETNCI(GETNCI("\\MGAMS.DATA:*","PATH"),"TIME_INSERTED"),,-1,"PCS"))')) > ...
%  uint64(mdsdata('MAXVAL(      GETNCI(GETNCI("\\MGAMS.DATA:*","PATH"),"TIME_INSERTED")           )')))
%  mgamsmsge(sprintf('Creating new shot #%d',shot),'r')
%  mdsdata('TREECLOSE("PCS",$)',shot)
%  mdslogical('PCS$DATA','MGAMS_DEVEL')
%  if mdstcl('SET TREE PCS'), end % to mask wrong errors
%  if mdstcl(sprintf('DELETE PULSE %d/ALL',shot)), end
%  if mdstcl(sprintf('CREATE PULSE %d',shot)), end
%  mdslogical('PCS$DATA')
% end
% if isempty(mdsopen('PCS',shot))
%  mgamsmsge(sprintf('Shot #%d cannot be open(created)',shot),'f')
%  return
% end
% if ~rem(mdssetdef('\MGAMS.DATA'),2)
%  mgamsmsge(sprintf('Shot #%d has no MGAMS branch',shot),'f')
%  return
% end
% set(H.shot,'string',int2str(mdsdata('$SHOT')))
% stat = 1;
% end

% item: mgamsitem(mode,pos,param,backgroundcolor)
function mgamsitem(mode,pos,par,col,parent)
 
 global H MP
 
 if nargin<5, parent=gcf;end
 par = MP.ind(par);
 k = find(H.varl,1,'last');
 if isempty(k), k = 1; else, k = k+1; end
 H.varl(k) = uicontrol('parent',parent,'style', mode ,'position',[pos(1) pos(2) 95 20],'string',MP.name(par),'callback',sprintf('mgams(''var'',%d)',k),'userdata',par,'backgroundcolor',col,'horizontalalignment','right');
 H.vare(k) = uicontrol('parent',parent,'style','edit','position',[pos(1)+100 pos(2) pos(3) 20],'callback','mgams modify','userdata',par,'background','y');
end

% analytic shape: r+j*z = mgamsanalytic(theta,neq)
function rz = mgamsanalytic(th,neq)
 global MP
 t = th(:);dum = zeros(numel(t),1);
 x0 = [MP.rmajo1(neq);MP.zmajo1(neq);MP.rmino1(neq);MP.cappa1(neq);MP.delta1(neq);MP.hlamd1(neq);t];
 F = analyticfit(dum,dum,x0);
 rz = reshape(F(1:numel(t))+1i*F(numel(t)+1:end),size(th));
end

function [F,J] = analyticfit(R,Z,x)
 t = x(7:end);
 Rf = x(1)+x(3)*cos(t+x(5)*sin(t)-x(6)*sin(2*t));
 Zf = x(2)+x(3)*x(4)*sin(t);
 F = [Rf-R;Zf-Z];
 if nargout > 1
  J = [ones(numel(R),1),...                                                                               % dRf/drmajo
       zeros(numel(R),1),...                                                                              % dRf/dzmajo
       cos(t+x(5)*sin(t)-x(6)*sin(2*t)),...                                                               % dRf/drmino
       zeros(numel(R),1),...                                                                              % dRf/dkappa
       -x(3)*sin(t+x(5)*sin(t)-x(6)*sin(2*t)).*sin(t),...                                                 % dRf/ddelta
       +x(3)*sin(t+x(5)*sin(t)-x(6)*sin(2*t)).*sin(2*t),...                                               % dRf/dhlamd
       diag(-x(3)*sin(t+x(5)*sin(t)-x(6)*sin(2*t)).*(1+x(5)*cos(t)-2*x(6)*cos(2*t)));... % dRf/dt
       zeros(numel(Z),1),...                                                                              % dZf/drmajo
       ones(numel(Z),1),...                                                                               % dZf/dzmajo
       x(4)*sin(t),...                                                                                    % dZf/drmino
       x(3)*sin(t),...                                                                                    % dZf/dkappa
       zeros(numel(Z),1),...                                                                              % dZf/ddelta
       zeros(numel(Z),1),...                                                                              % dZf/dhlamd
       diag(x(3)*x(4)*cos(t)                                                       ),...  % dZf/dt
       ];
  end
end
   

% time base: t = mgamstimebase
function t = mgamstimebase
 global MP
 t = mgamstimeeq2t(MP.timeeq);
end

% time base: t = mgamstimeeq2t(timeeq)
function t = mgamstimeeq2t(timeeq)
 global MP
 t = 0.01 + MP.timefac*(timeeq-MP.timeeq(1));
end

% time base: timeeq = mgamst2timeeq(t)
function timeeq = mgamst2timeeq(t)
 global MP
 timeeq = (t - 0.01) / MP.timefac + MP.timeeq(1);
end

% time base: t = mgamsallt
function t = mgamsallt
 global MP
 t1 = MP.timeeq(1:MP.numeq);
 t2 = fliplr(t1);
 t1 = 0.01 + MP.timefac * (t1 - t1(1));
 t2 = t1(end) + MP.flattop + MP.amagic * (t2(1) - t2);
 t = [t1 t2];
end

% message display: mgamsmsge(msge,level)
function mgamsmsge(msge,level)
 global H
 x = rem(round(clock),100);
 x = sprintf('%02d.%02d.%02d %02d:%02d:%02d',x([3 2 1 4 5 6]));
 if nargin == 1
  set(H.msge,'string',{x,msge},'foreground',[0 .5 0])
 elseif level == 'f'
  set(H.msge,'string',{x,msge},'foreground','r'), fprintf('\7')
 else
  set(H.msge,'string',{x,msge},'foreground','y')
 end
 drawnow
 fprintf('\r%s\t%s\n',x,msge);
end

% set up modal dialog box: mgamsmodal(button_string,callback_mode,help_str,arg1_name,arg1_default,arg1_help...)
function mgamsmodal(name,callback,help,varargin)
 global H
 if isempty(help), help=name;end
 set(H.dlgb,'string',name,'callback',['mgams ' callback ' eval'],'Tooltip',sprintf(help));
 set([H.dlgl,H.dlge],'string','','Visible','off');
 for k = 1:(nargin-3)/3
  tooltip=sprintf(varargin{3+(k-1)*3});
  if isempty(tooltip), tooltip=varargin{1+(k-1)*3};end
  set(H.dlgl(k),'Visible','on','string',varargin{1+(k-1)*3})
  set(H.dlge(k),'Visible','on','string',varargin{2+(k-1)*3},'Tooltip',tooltip);
 end
end

% evaluate arguments of modal dialog box
function varargout = mgamsmodaleval(varargin)
 % [AM 30.06.2017] Modifed to allow expressions in modal boxes
 % The result must evaluate to a scalar numeric value.
 global H MP MGceq %#ok<NUSED>
 % MP and MGceq can be used in expressions
 if nargout > numel(H.dlge)
   error('MatrixGenerator:modaleval','Too many output arguments');
 end
 varargout = cell(1,nargout);
 for k = 1:nargout
  name = get(H.dlgl(k),'string');
  try
   x = eval(['[' get(H.dlge(k),'string') ']']);
  catch ME
   x = [];
  end
  if isscalar(x) && isnumeric(x)
    varargout{k} = x;
    % Update string with evaluated value?
    if ~isempty(name), mgamsmsge(sprintf('[modal] "%s" is now %s',name,mgamsformat(x)));end
  else
    exception = MException('MatrixGenerator:modaleval','Wrong syntax for "%s"',name);
    if exist('ME','var'), exception = addCause(exception,ME);end
    throw(exception);
  end
 end
end

function O = mgamsloadfbte(shot)

 mdsopen('pcs',shot);
 psi    = tdi('tcv_eq($1,$2)','psi'     ,'runs');
 p_axis = tdi('tcv_eq($1,$2)','psi_axis','runs');
 r_axis = tdi('tcv_eq($1,$2)','r_axis'  ,'runs');
 z_axis = tdi('tcv_eq($1,$2)','z_axis'  ,'runs');
 r_xpts = tdi('tcv_eq($1,$2)','r_xpts'  ,'runs');
 z_xpts = tdi('tcv_eq($1,$2)','z_xpts'  ,'runs');
 
 has_fbte = isodd(psi.status);
 if has_fbte
  neq    = numel(psi.dim{3});
 else
  neq    = 0;
 end
 single_axis = (numel(p_axis.data) == neq);

 if has_fbte && size(psi.data,3) == neq && single_axis % Exclude doublets for now
  psin   = 1 - psi.data./reshape(p_axis.data,[1,1,neq]);
   
  v = [0.9,1.0,1.1];
  ncont = numel(v);
   
  Cr = cell(ncont,neq);
  Cz = cell(ncont,neq);
  for ieq = 1:neq
   C = cell(ncont,1);
   if ncont == 1
    cont = v*[1 1];
   else
    cont = v;
   end
   Cmat = contourc(psi.dim{1},psi.dim{2},psin(:,:,ieq).',cont);
   ik = 1;
   while ik < size(Cmat,2)
    kk = find(v == Cmat(1,ik));
    C{kk} = [C{kk},[NaN;NaN],Cmat(:,ik+1:ik+Cmat(2,ik))];
    ik = ik+Cmat(2,ik)+1;
   end
   for kk=1:ncont
    if isempty(C{kk}), C{kk} = [NaN;NaN];end
    Cr{kk,ieq} = C{kk}(1,:);
    Cz{kk,ieq} = C{kk}(2,:);
   end
  end
   
  O.r_axis = reshape(r_axis.data,[],neq);
  O.z_axis = reshape(z_axis.data,[],neq);
  O.r_xpts = reshape(r_xpts.data,[],neq);
  O.z_xpts = reshape(z_xpts.data,[],neq);
  O.Cr     = Cr;
  O.Cz     = Cz;
  
 else
   
  O.r_axis = NaN;
  O.z_axis = NaN;
  O.r_xpts = NaN;
  O.z_xpts = NaN;
  O.Cr     = {NaN;NaN;NaN};
  O.Cz     = {NaN;NaN;NaN};
   
 end
   
 
end

function s = run_mga(shot,version)
 % RUN_MGA run the MGA part of the shot preparation

 global MP

 switch version
  case 'Fortran'
   s = unix(sprintf('mga %d > /tmp/$USER/mga-%d.log',shot,shot));
  case 'Matlab'
   s = 0; mga(shot);
  case 'Matlab compatibility mode'
   s = 0; mga(shot,true);
  case 'Matlab checked'
   % Run comparison tool - will overwrite shot XXX999
   fort_equiv = true;
   [~,~,~,~,S] = mgaruncmp('all',shot,fort_equiv,Inf);
   % Check match within allowed tolerance
   tol  = 0.002;
   tolV = 0.010;
   ok = true;
   for F=fieldnames(S)'
     switch F{:}
       case 'EFVOLT', ok = ok && all(S.(F{:})(:) < tolV);
       otherwise,     ok = ok && all(S.(F{:})(:) < tol );
     end
   end
   % Run on actual shot
   if ok
    s = 0;mga(shot,fort_equiv);
   else
    % Display info message
    mgamsmsge('Comparison failed, saving offending parameters and falling back to MGA-Fortran','w');
    % First save offending parameter structure
    MG = mgp(shot);
    filename = sprintf('/home/tcv_oper/mgacmp-fail/mg_%s_%s.mat',getenv('USER'),datestr(now,'YYYYmmDDHHMMSS'));
    save(filename,'MG');
    %
    s = unix(sprintf('mga %d > /tmp/$USER/mga-%d.log',shot,shot));
   end
 end

 if MP.iscramb>1000
  if isempty(which('optimize_zctrl_mga'))
   oldpath = path;
   restore_path = onCleanup(@() path(oldpath)); % restore path even upon errors
   % tmp folders for Pesamosca vertical control (while not deployed yet)
   newpath = {'/home/pesamosc/VC_tmp/meq',...
              '/home/pesamosc/VC_tmp/rzip',...
              '/home/pesamosc/VC_tmp/optimized_vertical_control',...
              '/home/pesamosc/VC_tmp/optimized_vertical_control/operations'};
   addpath(newpath{:}); % Single addpath call
   add_rzip_paths;
  end
  % run optimization 
  optimize_zctrl_mga(shot);
 end
end

function s = run_fbte(shot)
 % RUN_FBTE run the FBTE part of the shot preparation

 global MP
 
 list = {'r','z','b','fa','fb','fe','br','bz','ba','be','cr','cz','ca','ce'};
 nvar = numel(list);
 
 rows0 = size(MP.(['gp',list{1}]),1);
 for ii = 2:nvar
   assert(size(MP.(['gp',list{ii}]),1) == rows0, 'MGAMS:Inconsistent number of rows between gp%s and gp%s',list{1},list{ii});
 end
 
 [~,~,LY] = fbte(shot);
 
 s = ~all(LY.isconverged);
 if s
   mask = find(~LY.isconverged);
   list = strjoin(cellstr(num2str([mask;LY.t(mask)].','%d[t=%5.3fs]')),' ');
   fprintf('FBTE did not converge for slices %s\n',list);
 end
 
end
