function [MG,FB,MGA] = mga(shot,fort_equiv)
%function [MG,FB,MGA] = mga(shot,fort_equiv)

%%% J.Sinha 15 Nov. 2017 %%%
% New vairable nbackoff introduced to include the new single-axis
% breakdown scenario
% nbackoff=16 to use the new scenario
% NVVEL=1038 does not work for the doublet - problem with fbt.f???

if nargin < 2 || isempty(fort_equiv), fort_equiv = false;
else,                                 fort_equiv = logical(fort_equiv);end

%% Read MGAMS parameters
MG=mgp(shot);
if fort_equiv
 % Convert to single precision, to avoid round-off differences
 MG = single(MG);
 MG = fix_int(MG);
end

%% Convert to struct with additional parameters from debug_p2
MG = struct_with_debug(MG,2);

 %% iscramb decoding
 % iscramb values larger than 1000 indicate that we are using the full 
 % F-coil set for Ipz control [Work by Federico Pesamosca]
 % pending proper parameter implementation.
 MG.iscramb = rem(MG.iscramb,1000);
 
 % Use iscramb also to set another temp variable, pending proper
 % implementation as mgams parameters

 % if abs(iscramb)>=100 : don't use use R part in M matrix - otherwise use it (legacy)
 inores     = logical(fix(MG.iscramb/100 ));   
 MG.iscramb = rem(MG.iscramb,100);
 
 %% Discretize analytical equilibria
 for k = find(MG.iansha)
  w = (0:MG.ilia(k)-1)*(2*pi/MG.ilia(k));
  MG.rlia1(1:MG.ilia(k),k) = MG.rmajo1(k) + MG.rmino1(k)*cos(w+MG.delta1(k)*sin(w)-MG.hlamd1(k)*sin(w*2));
  MG.zlia1(1:MG.ilia(k),k) = MG.zmajo1(k) + MG.rmino1(k)*MG.cappa1(k)*sin(w);
 end
%  if isnan(MG.ohcorr), MG.ohcorr = 0; end
%  MG.gainext(end+1:3) = 0;
%  MG.aipgain(end+1:2) = 0;
 
 %% Read FBT, limit to numeq equilibria
 mdso(shot)
 FB.rtp = mdsr('\FBT:RAX'        )';
 FB.ztp = mdsr('\FBT:ZAX'        )';
 FB.rbp = mdsr('\FBT:RAX2'       )';
 FB.zbp = mdsr('\FBT:ZAX2'       )';
 FB.Fa  = mdsr('\FBT:SHVSEC'     )';
 FB.Iea = mdsr('\FBT:GCURR'      ) ; FB.Iea = [FB.Iea(1:16,:);zeros(2,size(FB.Iea,2))];
 FB.Bm  = mdsr('\FBT:BEE'        ) ;
 FB.Ff  = mdsr('\FBT:FLUX'       ) ;
 FB.rze = mdsr('\FBT:RZE'        )';
 FB.Ip  = mdsr('\FBTE:I_P'       )';
 FB.li  = mdsr('\FBTE:L_I'       )';
 FB.Ap  = mdsr('\FBTE:AREA'      )';
 FB.Kp  = mdsr('\FBTE:KAPPA_EDGE')';
 for k = fieldnames(FB)', FB.(k{1}) = FB.(k{1})(:,1:MG.numeq); end
 FB.Ff  = FB.Ff-FB.Ff(1,:);
 FB.t   = MG.toft(end) + (MG.timeeq(1:MG.numeq) - MG.timeeq(1))*MG.timefac;
 FB     = orderfields(FB);
 
 %% Checks
 if (fort_equiv && MG.numeq  > MG.nruns), MG.numeq  = MG.nruns; end
 if (fort_equiv && MG.nzaxre > MG.numeq), MG.nzaxre = MG.numeq; end
 if (fort_equiv && MG.nzaxel > MG.numeq), MG.nzaxel = MG.numeq; end
 assert(ismember(MG.iscramb,[-2:0 2:12 61:63])                ,'TCVCS:MGA:InvalidInputData','iscramb'                             )
 assert(MG.zshift(1) == 0 || all(MG.zshift(2:3)) || fort_equiv,'TCVCS:MGA:InvalidInputData','zshift = [~0,0,0]'                   )
 assert(MG.nzaxre <= MG.numeq                                 ,'TCVCS:MGA:InvalidInputData','nzaxre > numeq'                      )
 assert(MG.nzaxel <= MG.numeq                                 ,'TCVCS:MGA:InvalidInputData','nzaxel > numeq'                      )
 assert(MG.gainvze == 0                                       ,'TCVCS:MGA:InvalidInputData','gainvze ~= 0'                        )
 assert(MG.gain >= 0                                          ,'TCVCS:MGA:InvalidInputData','gain < 0'                            )
 assert(ismember(MG.nfast,[-5 -3 -2 0])                       ,'TCVCS:MGA:InvalidInputData','nfast = -1,-4'                       )
 assert(~ismember(rem(MG.mvloop,10),[5 6])                    ,'TCVCS:MGA:InvalidInputData','mvloop = 5,6,15,16'                  )
 assert(MG.ierat~= 4||~ismember(MG.mvloop,[11 21])            ,'TCVCS:MGA:InvalidInputData','mvloop = 11,21, ierat=4'             )
 assert(MG.ierat~=-1||~ismember(MG.mvloop,[11 21])            ,'TCVCS:MGA:InvalidInputData','mvloop = 11,21, ierat=-1'            )
 assert(MG.iscramb>-2||~MG.nfast                              ,'TCVCS:MGA:InvalidInputData','iscramb < -1, nfast~=0'              )
 assert(numel(MG.toft) == numel(MG.diohdt)                    ,'TCVCS:MGA:InvalidInputData','toft and diohdt have different sizes')
 assert(~(fort_equiv && MG.iscramb<0)                         ,'TCVCS:MGA:InvalidInputData','iscramb <0, fort_equiv=true'         )
 assert(~(fort_equiv && MG.mvloop==8)                         ,'TCVCS:MGA:InvalidInputData','mvloop=8 (doublet), fort_equiv=true' )
 assert(all(FB.rtp)                                           ,'TCVCS:FBT:NoConvergence'   ,'FBT no convergence'                  )
 assert(all(FB.rze>MG.ri & FB.rze<MG.ro)                      ,'TCVCS:MGA:InvalidInputData','rze is out of the standard FBT grid' )
 
 if ~exist('nbackoff','var')
  nbackoff=0;
 end
 
 %% Negative plasma current
 if MG.iohfb == -1
  MG.dpszero   = -MG.dpszero;
  MG.rshift3   = -MG.rshift3;
  MG.delipz    = -MG.delipz;
  MG.aipipz    = -MG.aipipz;
  MG.zeecorr(1)= -MG.zeecorr(1);
  MG.relovo    = -MG.relovo;
  MG.brmzero   = -MG.brmzero;
  MG.brpzero   = -MG.brpzero;
  MG.brmzerb   = -MG.brmzerb;
  MG.brpzerb   = -MG.brpzerb;
  MG.bzmzero   = -MG.bzmzero;
  MG.bzpzero   = -MG.bzpzero;
  MG.bzmzerb   = -MG.bzmzerb;
  MG.bzpzerb   = -MG.bzpzerb;
  MG.aipdel    = -MG.aipdel;
  MG.diohdt    = -MG.diohdt;
  MG.zshift(1) = -MG.zshift(1);
  FB.Fa        = -FB.Fa;
  FB.Iea       = -FB.Iea;
  FB.Bm        = -FB.Bm;
  FB.Ff        = -FB.Ff;
 end
 
 %% Geometry
 if fort_equiv
  % Use translation from FORTRAN version
  G = fbtmeas(MG.nvvel);
 else
  % Dimension/system selection
  G.na = 18;
  dima = mdsdata('["E_"//III(1..8),"F_"//III(1..8),"OH_"//III(1..2)]');
  dimt = mdsdata('"T_"//III(1..3)');
  dimf = mdsdata('III(1..38)');
  dimm = mdsdata('III(1..38)');
  G.nv = MG.nvvel;
  switch G.nv
   case  19,  selv =  'R'; dimv = mdsdata('STATIC("DIM_R")');
   case  38
    if MG.mvloop==8 || nbackoff==16
              selv = 'SX'; dimv = mdsdata('STATIC("DIM_SX")');
    else
              selv =  'S'; dimv = mdsdata('STATIC("DIM_S")');
    end
   case 256,  selv =  'V'; dimv = mdsdata('STATIC("DIM_V")');
   case 1038, selv = 'SX'; dimv = mdsdata('STATIC("DIM_SX")'); G.nv=38;
  end
  if G.nv < 0,selv =  'E'; dimv = mdsdata('"E_"//III(1..$1)',-G.nv); G.nv = -G.nv; end
  
  % Magnetics
  G.nm = numel(dimm);
  G.nf = numel(dimf);
  G.rm = mdsdata('STATIC("R_M")[$1]',dimm);
  G.zm = mdsdata('STATIC("Z_M")[$1]',dimm);
  
  % Mesh
  G.nrx = mdsdata('STATIC("NA_XF")');
  G.nzx = mdsdata('STATIC("NB_XF")');
  G.nx  = G.nrx*G.nzx;
  G.rx  = mdsdata('STATIC("R_XF:AXIS")');
  G.zx  = mdsdata('STATIC("Z_XF:AXIS")');
  G.ax  = mdsdata('STATIC("A_XF")');
  G.bx  = mdsdata('STATIC("B_XF")');
  [G.rrx,G.zzx] = meshgrid(G.rx,G.zx);
  
  % Conductors
  G.rw = mdsdata('STATIC("R_W")');
  G.zw = mdsdata('STATIC("Z_W")');
  G.rt = mdsdata('STATIC("R_C")[$1]',dimt); % can be improved using w
  G.zt = mdsdata('STATIC("Z_C")[$1]',dimt);
  G.ru = mdsdata('STATIC("R_V")');
  G.zu = mdsdata('STATIC("Z_V")');
  
  
  % Green's functions
  % transformation matrices
  G.Twa = mdsdata('STATIC("T_W_C")')*mdsdata('STATIC("T_C_A")[*,$1]',dima);
  if G.nv == 256, G.Tuv = 1;
  else
   G.Tuv = mdsdata('STATIC("T_V_"//$1)[*,$2]',selv(1),dimv); 
  end
  G.Ttt = mdsdata('STATIC("T_C_A")[$1,"TOR_001"]',dimt);
  % resistance and mutual inductance
  G.Ra  = mdsdata('STATIC("RES_A")[$1]',dima);
  G.Rv  = mdsdata('STATIC("RES_"//$1)[$2]',selv,dimv);
  G.Maa = mdsdata('STATIC("MUT_A_A")[$1,$2]',dima,dima);
  G.Mvv = mdsdata('STATIC("MUT_"//$1//"_"//$1)[$2,$2]',selv,dimv);
  G.Mva = mdsdata('STATIC("MUT_"//$1//"_A")[$2,$3]',selv,dimv,dima);
  % Rv and Ra used in .FOR
  % tmp = [1.2711 1.2738 1.2738 1.2738 1.2738 1.2738 1.2895 1.3213 1.6773 2.0924 2.3368 2.7809 2.6548 2.1638 3.0275 3.5829 3.3532 3.2480 3.4017 ...
  %        3.7648 3.4017 3.2480 3.3153 3.4573 2.9508 2.0713 2.3129 2.5959 2.5119 2.1984 1.7093 1.3174 1.2824 1.2738 1.2738 1.2738 1.2738 1.2738]' / 1000;
  % switch nv
  %  case 19, tmp = 1./tmp; G.Rv = 1./(.5*tmp([38 2:2:36])+tmp(1:2:37)+0.5*tmp(2:2:38));
  %  case 38, G.Rv = tmp;
  % end
  % G.Ra(OH1) = G.Ra(OH1) - 0.00070;
  % probe field and loop flux from conductors
  G.Mfa = mdsdata('STATIC("MUT_F_A")[$1,$2]',dimf,dima);
  G.Mfv = mdsdata('STATIC("MUT_F_"//$1)[$2,$3]',selv,dimf,dimv);
  G.Mfx = mdsdata('STATIC("MUT_F_XF")[$1,*]',dimf);
  G.Mfa = G.Mfa-G.Mfa(1,:);
  G.Mfv = G.Mfv-G.Mfv(1,:);
  G.Mfx = G.Mfx-G.Mfx(1,:);
  G.Bma = mdsdata('STATIC("BT_M_A")[$1,$2]',dimm,dima);
  G.Bmv = mdsdata('STATIC("BT_M_"//$1)[$2,$3]',selv,dimm,dimv);
  G.Bmx = mdsdata('STATIC("BT_M_XF")[$1,*]',dimm);
  
  % Function handle to compute green functions
  G.greenem = @greenem;
 end
 % Coil indexing
 tmp = num2cell(1:G.na);
 [E1,E2,E3,E4,E5,E6,E7,E8,F1,F2,F3,F4,F5,F6,F7,F8,OH1,OH2] = deal(tmp{:});
 EF = E1:F8; OH = [OH1 OH2];

 %% Midplane selection
 if MG.mvloop == 8 || MG.mvloop == 18, midplan =  8; % doublet
  if fort_equiv
   %NOTE: not sure how this impacts J.Sinha doublet shots
   if     FB.rbp(1) < 0.01 && FB.ztp(1) < 0, iplatop = -1;
   elseif FB.rbp(1) < 0.01 && FB.ztp(1) > 0, iplatop =  1;
   else,                                     iplatop =  0;
   end
  else
   if     FB.rbp(1) < 0.01 && FB.ztp(1) < 0, iplatop = -1;
   elseif FB.rbp(1) > 0.01 && FB.ztp(1) < 0, iplatop =  1;
   else,                                     iplatop =  0;
   end
  end
  k = find(FB.rbp < 0.01); FB.rbp(k) =  FB.rtp(k); FB.zbp(k) = -FB.ztp(k);
  k = find(FB.ztp < 0   ); FB.ztp(k) = -FB.ztp(k); FB.zbp(k) = -FB.zbp(k);
 elseif fort_equiv
  % NOTE: Is midplan==[5,35] a .m novelty?
  if     FB.ztp(MG.nzaxre) < -0.11, midplan = 37;
  elseif FB.ztp(MG.nzaxre) >  0.11, midplan =  3;
  else,                             midplan =  1;
  end
 else
  if     FB.ztp(        1) < -0.35, midplan = 35;
  elseif FB.ztp(        1) < -0.11, midplan = 37;
  elseif FB.ztp(        1) >  0.35, midplan =  5;
  elseif FB.ztp(        1) >  0.11, midplan =  3;
  else,                             midplan =  1;
  end
 end

 %% Finite elements
 if midplan == 8, zy0 = 0;
 else,            zy0 = FB.ztp(MG.nzaxel); end
 if MG.iscramb ~= 4
  zutil = 0.75; tmp = MG.zup - MG.zlp;
  MG.zlp = zy0 - 0.5*tmp; MG.zup = zy0 + 0.5*tmp;
  if fort_equiv
   % This shrinks the FE extent to fit within the allowed space and keeps
   % it centered around zy0
   if MG.zlp < -zutil, MG.zlp = -zutil        ; MG.zup =  zutil + 2*zy0; end
   if MG.zup >  zutil, MG.zlp = -zutil + 2*zy0; MG.zup =  zutil        ; end
  else
   % This shifts the FE extent to fit within the allowed space and keeps
   % its original size set by MG.zup - MG.zlp
   if MG.zlp < -zutil, MG.zlp = -zutil      ; MG.zup = -zutil + tmp; end
   if MG.zup >  zutil, MG.zlp =  zutil - tmp; MG.zup =  zutil      ; end
  end
 end
 if midplan == 8
  MG.nelz  = 6;
  MG.xip =  0.65; MG.xop = 1.11;
  MG.zlp = -0.70; MG.zup = 0.70;
 end
 G.ay = (MG.xop - MG.xip)/3;
 G.by = (MG.zup - MG.zlp)/(MG.nelz+1);
 G.ry = MG.xip + (1:2)*G.ay;
 G.zy = MG.zlp + (1:MG.nelz)*G.by;
 % mesh to finite elements
 G.Txy = reshape(max(1-abs(reshape(G.zx,[G.nzx,    1,      1,      1]) - ...
                           reshape(G.zy,[    1,    1,MG.nelz,      1]))/G.by,0) .* ...
                 max(1-abs(reshape(G.rx,[    1,G.nrx,      1,      1]) - ...
                           reshape(G.ry,[    1,    1,      1,      2]))/G.ay,0),...
                 G.nx,2*MG.nelz);
 % probe field and loop flux from conductors
 G.Mfy = G.Mfx * G.Txy; % Relative to loop 1
 G.Bmy = G.Bmx * G.Txy;
 
 %% Ip observer
 Ipm = sqrt(diff(G.rm([end,1:end])).^2 + diff(G.zm([end,1:end])).^2).';
 Ipm = (Ipm+Ipm([2:end 1]))/(8e-7*pi);
 
 %% Radial observer
 if midplan == 8, kRad = [5 16;35 24];
 else,            kRad = mod((midplan+min(max(MG.iscale,-2),2))*[1 -1]-[1 18],38)+1;
 end
 
 Radf = zeros(size(kRad,1),G.nf); Radf(1,kRad(1,:)) =      [       1         -1]                  *MG.psifac;
 Radm = zeros(size(kRad,1),G.nm); Radm(1,kRad(1,:)) = 2*pi*[MG.gapin -MG.gapout].*G.rm(kRad(1,:))'*MG.psifac;
 Radf(end,kRad(end,:)) = Radf(1,kRad(1,:));
 Radm(end,kRad(end,:)) = Radm(1,kRad(1,:));
 
 %% Vertical observer
 if midplan == 8
  zp0  = FB.ztp(MG.nzaxre); % z-position of top plasma
  zp02 = FB.zbp(MG.nzaxre); % z-position of bottom plasma
 else
  zp0 = FB.ztp(MG.nzaxre);
 end
 Wf = repmat(MG.wscal,G.nf,1);
 Wm = repmat(MG.vscal,G.nm,1);
 if MG.ipr16   == 0, Wf([10:12 28:30]) = 0; end
 % Always deselect probes and loops close to G-coils
 Wf([13:16 24:27]) = 0;
 Wm([13:16 24:27]) = 0;
 Wd = Wm; Wd([1 7 9 10 12:16 20 24:28 30 31 33]) = 0.;
 if fort_equiv
  Dzy = repmat((G.zy - zp0)*G.ay/G.ax*G.by/G.bx,1,2);
  Ipy = G.ay/G.ax*G.by/G.bx*ones(size(Dzy));
 else
  Dzy = (G.zzx(:)' - zp0)*G.Txy;
  Ipy = sum(G.Txy);
 end
 [Ipzf ,Ipzm ,Ipza ,Ipzv ] = Iy(Wf,Wm,MG.uc2*[1 1],G.Mfa,G.Mfv,G.Mfy,G.Bma,G.Bmv,G.Bmy,Ipm,Ipy,Dzy);
 [~    ,Ipzmm,Ipzma,Ipzmv] = Iy(0 ,Wm,MG.uc2*[1 1],G.Mfa,G.Mfv,G.Mfy,G.Bma,G.Bmv,G.Bmy,Ipm,Ipy,Dzy);
 if fort_equiv
  [~    ,Ipzdm            ] = Iy(0 ,Wd,MG.uc2*[0 1],G.Mfa,G.Mfv,G.Mfy,G.Bma,G.Bmv,G.Bmy,Ipm,Ipy,Dzy);
 else
  [~    ,Ipzdm            ] = Iy(0 ,Wd,MG.uc2*[0 0],G.Mfa,G.Mfv,G.Mfy,G.Bma,G.Bmv,G.Bmy,Ipm,Ipy,Dzy);
 end
 if midplan == 8
  if MG.ikriz == 1 % choose b/w new vertical observer or F.Hofmann's observer for doublet
   % new vertical observer for doublet
   % remove the flux loops at the top and bottom
   Wf([10:12 28:30]) = 0;
   TT = reshape(G.Txy,[G.nzx G.nrx 2*MG.nelz]);
   T_t = reshape(TT(G.zx>0,:,:),[32*G.nrx 2*MG.nelz]);
   T_b = reshape(TT(G.zx<0,:,:),[32*G.nrx 2*MG.nelz]);
   Dzy1 = reshape(G.zzx(G.zx>0,:),[1 32*G.nrx])-zp0;
   Dzy2 = reshape(G.zzx(G.zx<0,:),[1 32*G.nrx])-zp02;
   Dzy_t = Dzy1*T_t;%dZ for top plasma
   Dzy_b = Dzy2*T_b;%dZ for bottom plasma
   % vertical observer for top plasma
   [Ipzf_t ,Ipzm_t ,Ipza_t ,~ ] = Iy(Wf,Wm,MG.uc2*[1 1],G.Mfa,G.Mfv,G.Mfy,G.Bma,G.Bmv,G.Bmy,Ipm,Ipy,Dzy_t);
   % vertical observer for bottom plasma
   [Ipzf_b ,Ipzm_b ,Ipza_b ,~ ] = Iy(Wf,Wm,MG.uc2*[1 1],G.Mfa,G.Mfv,G.Mfy,G.Bma,G.Bmv,G.Bmy,Ipm,Ipy,Dzy_b); 
  else % F.Hofmann's observer for doublet
   v_fc  = G.Mfa(2:end,[ 9 12]      )*[1;-1];
   %   v_fc  = sum(G.Mfa(2:end,[ 9 12]      ),2);
   v_fx1 = sum(G.Mfy(2:end,[ 1  7]      ),2);
   v_fx2 = sum(G.Mfy(2:end,[ 2  8]      ),2);
   v_fx3 = sum(G.Mfy(2:end,[ 6 12  5 11]),2);
   v_fx4 =     G.Mfy(2:end,[ 6 12  5 11])*[1;1;-1;-1];
   if MG.ipr16 == 0
    k = [9 10 11 27 28 29];
    v_fc (k,:) = 0;
    v_fx1(k,:) = 0;
    v_fx2(k,:) = 0;
    v_fx3(k,:) = 0;
    v_fx4(k,:) = 0;
   end
   A = [eye(37) v_fx1 v_fx2 v_fc v_fx3; [[v_fx1 v_fx2 v_fc v_fx3]' zeros(4)]];
   b = [v_fx4 ; zeros(4,1)];
   S = diag([(1e-4*ones(37,1));ones(4,1)]); % scaling for numerical conditioning
   v = (S*A)\(S*b);
   
   IpDz_xf_t = [0 2*G.ay*G.by/(G.ax*G.bx)*G.by/(v(1:37)'*v(1:37)).*v(1:37)'];
   IpDz_xf_b = [IpDz_xf_t(1) fliplr(IpDz_xf_t(2:38))];
   %   IpDz_xf_b = fliplr(IpDz_xf_t);
   v_mc  = G.Bma(:,[9 12]       )*[1;-1];

%   v_mc  = sum(G.Bma(:,[9 12]       ),2);
   v_mx1 = sum(G.Bmy(:,[ 1 7]       ),2);
   v_mx2 = sum(G.Bmy(:,[ 2 8]       ),2);
   v_mx3 = sum(G.Bmy(:,[ 6 12  5 11]),2);
   v_mx4 =     G.Bmy(:,[ 6 12  5 11])*[1;1;-1;-1];
   
   
   A = [eye(38) v_mx1 v_mx2 v_mc v_mx3; [[v_mx1 v_mx2 v_mc v_mx3]' zeros(4)]];
   b = [v_mx4 ; zeros(4,1)];
   S = diag([(1e-4*ones(38,1));ones(4,1)]); % scaling for numerical conditioning
   v = (S*A)\(S*b);
   IpDz_xm_t = 2*G.ay*G.by/(G.ax*G.bx)*G.by/(v(1:38)'*v(1:38)).*v(1:38)';
   IpDz_xm_b = [IpDz_xm_t(1) fliplr(IpDz_xm_t(2:38))];
%   IpDz_xm_b = fliplr(IpDz_xm_t);
  end
 end
 
 %% Elongation observer
 if fort_equiv
  if midplan == 8, Dky = repmat(G.zy.^2,1,2);
  else
   switch rem(MG.ifour,10)
    case 4, Dky = abs(G.zy - zp0);
            Dky = repmat(Dky*G.ay/G.ax*G.by/G.bx,1,2);
    case 5, Dky =    (G.zy - zp0).^2;
            Dky = repmat(Dky*G.ay/G.ax*G.by/G.bx,1,2);
    case 6
     rb = [MG.rlim1(1:MG.ilie(MG.nzaxre),MG.nzaxre) ; MG.rlia1(1:MG.ilia(MG.nzaxre),MG.nzaxre)];
     zb = [MG.zlim1(1:MG.ilie(MG.nzaxre),MG.nzaxre) ; MG.zlia1(1:MG.ilia(MG.nzaxre),MG.nzaxre)];
     k = [imax(zb) imin(zb) imax(rb) imin(rb)];
     Dky = 69000*[1 1 -1 -1]*G.greenem('mut',rb(k),zb(k),G.rrx,G.zzx)*G.Txy;
    otherwise, Dky = zeros(1,2*MG.nelz);
   end
  end
 else
  if midplan == 8, Dky = G.zzx(:)'.^2*G.Txy;
  else
   switch rem(MG.ifour,10)
    case 2, Dky = abs(G.zzx(:)' - zp0)   *G.Txy;
    case 3, Dky =    (G.zzx(:)' - zp0).^2*G.Txy;
    case 4, Dky = abs(repmat(G.zy,1,2) - zp0)   .*sum(G.Txy);
    case 5, Dky =    (repmat(G.zy,1,2) - zp0).^2.*sum(G.Txy);
    case 6
     rb = [MG.rlim1(1:MG.ilie,MG.nzaxre) ; MG.rlia1(1:MG.ilia,MG.nzaxre)];
     zb = [MG.zlim1(1:MG.ilie,MG.nzaxre) ; MG.zlia1(1:MG.ilia,MG.nzaxre)];
     k = [imax(zb) imin(zb) imax(rb) imin(rb)];
     Dky = 69000*[1 1 -1 -1]*G.greenem('mut',rb(k),zb(k),G.rrx,G.zzx)*G.Txy;
    otherwise, Dky = zeros(1,2*MG.nelz);
   end
  end
 end
 [Ipkf,Ipkm,Ipka] = Iy(Wf,Wm,MG.uc2*[1 1],G.Mfa,G.Mfv,G.Mfy,G.Bma,G.Bmv,G.Bmy,Ipm,Ipy,Dky);
 
 %% Backoff and coils for quadrupole setting (BACKOFF)
 % initial conditions were computed for OHCUR=10kA and VLOOP=7V corresponding to dioh=15504A/s
 [Baa,Bav] = deal(zeros(1,8));
 switch midplan
  case  1        
   if nbackoff == 16
    kq  = [  E4   E5   F4   F5 ];
    kb  =    [E3  E4  E5   E6   F3    F4    F5    F6];
    Baa =    [15  -4  -6   27   73    220   303   40];
    Bav =    [32  38  30   55  -374  -144   280  -641];
    Ba_tor = [5   -2  -1  -2    215  -193   82   -74];
   else
    kq  = [  E4   E5   F3   F6 ];
    kb  = [  F2   F4   F5   F7   E1];
    if MG.if36fb == 1
     Baa = [ 244  167  237  232  190];
     Bav = [-442 -172 -149 -450    0];
     switch MG.iscramb
      case 5
       if MG.bzero(1) < 1.3, Bav = [-292 -122  -99 -300    0]; end
      case 7,                Bav = [-532  -92 -124 -440    0];
      case {8 10},           Bav = [-476 -178 -116 -439    0];
     end
    else
     Baa = [ 112  249  240  102  -29]; % Caa = [ 342  174  273  342  135];
     Bav = [-477 -205 -188 -494 -119]; % Cav = [-586 -167 -195 -412    0];
    end
   end
  case 37
   if nbackoff==16
    kq  = [  E3   E4   F3   F4 ];
    kb  =    [E2     E3   E4   E5    F2   F3      F4    F5];
    Baa =    [145    28   21   27  -333   807    -333  473];
    Bav =    [64     17   22   7   -463  -11     -1   -478];
    Ba_tor = [-2    -2   -1   -3    149  -88     51   -53];   
   else
    kq  = [  E3   E4   F2   F5 ];
    kb  = [  F1   F3   F4   F6   E1];
    if MG.if36fb == 1
     Baa = [ -34  225  202  330  100];
     Bav = [-442  -90 -162 -405  100];
    else
     Baa = [-154  225  202  230   20]; % Caa = [ -54  244  188  226   48];
     Bav = [-442  -90 -162 -405   20]; % Cav = [-597  -68 -205 -511  102];
    end
   end
  case 35  %z=-0.4m
       kq  = [E2   E3   F2   F3];
       kb = [E1 E2 E3 E4 F1 F2 F3 F4];
       Baa=[280         175          65         178       -1585        1790        -864        1100];
       Bav=[336  122  117  161   -771  220   97   -599]; 
       Ba_tor=[7.4667   5.0008   5.2259   7.0172 122.7498   0.6806  23.4234   -55.4160];    
  case  3
      if nbackoff ==16
       kq  = [  E5   E6   F5   F6 ];
       kb  =    [E4    E5   E6   E7   F4    F5     F6    F7];
       Baa =    [27    21   28   145  473  -333    807  -333];
       Bav =    [10    25   20   67  -490   4     -18   -435];
       Ba_tor = [1    -5   -3   -6    89   -131    94    -53];    
      else   
       kq  = [  E5   E6   F4   F7 ];
       kb  = [  F3   F5   F6   F8   E8   F1];
       if MG.if36fb == 1
        Baa = [ 310  224  215  -51  100    0];
        Bav = [-561 -247 -221 -115  100    0];
       else               
        Baa = [ 230  202  225  -34    0  -80]; % Caa = [ 267  180  253  -46   23  -91];
        Bav = [-419 -172  -82 -329    0  -80]; % Cav = [-222 -290  -17 -563  -17 -529];
       end
      end
  case 5 %z=+0.4m
       kq = [E6 E7 F6 F7];
       kb=[ E5 E6 E7 E8 F5 F6 F7 F8];
       Baa=[178 65 175 280 1100 -864 1790 -1585];
       Bav=[180   125   132   354  -614    59   263  -775];
       Ba_tor=[-8.5723    -9.9462    -7.7926   -16.8144  50.1193   -72.5403  73.3451   -79.0939];    
  case 8
   if MG.ifour == 8
    kq  = [  F1   F2   F3   F4   F5   F6   F7   F8];
    kb  = [  E2   E3   E6   E7   E1   E4   E5   E8];
    Bav = [   0    0    0    0 -100 -200 -200 -100];
   else 
    kq=[E2 E3 E6 E7 F2 F3 F6 F7];
    kb=[E1 E2 E3 E4 E5 E6 E7 E8 F1 F2 F3 F4 F5 F6 F7 F8];
    %new back-off and lower gradient Ebias=0;
    Baa=[689   374   317   267   267   317   374   689,...
        -187   283   309    37    37   309   283  -187];
    Bav=[320    99    99    81    90    98   102   316,...
        -996   489  -259  -143   -2  -396   644  -1082]/1.1587;
    Ba_tor=[1.9517   2.0239    1.9070    2.0065    1.6709    1.9707    1.8689    2.7846,...
           98.4779   30.1930   -8.3251  -15.9585  -13.6802   -6.2601   -6.8285   -9.3872];
   end
 end
 Baa = MG.ohback/ 21267*Baa';
 Bav = MG.veback/168780*Bav';
 
 %% Control coil combinations
 om = 1 - MG.alpha;
 op = 1 + MG.alpha;
 if midplan == 8
%   kc  = [F1 F4 F5 F8];
%   c = [ 0.540 -0.580  0.264 -0.064; % non-orthogonal
%         0.226  0.908 -0.765  0.212;
%         0.212 -0.765  0.908  0.226;
%         0.064 -0.264  0.580 -0.540];
  kc  = [F2 F4 F5 F7];
  c   = [ 0.540 -0.580  0.264 -0.064; % non-orthogonal
          0.226  0.908 -0.765  0.212;
          0.212 -0.765  0.908  0.226;
          0.064 -0.264  0.580 -0.540];
  r = c;
 else
  switch MG.iscramb
   case {0,-1,-2}
    kc = [];
    c = [];
    r = [];
   case  2
    if nbackoff==16
     switch midplan
      case  1, kc = [F2 F3 F6 F7];
      case  3, kc = [F3 F4 F7 F8]; 
      case 37, kc = [F1 F2 F5 F6];
      case 35, kc = [F1 F2 F5 F6];
     end
    else,      kc = [kb(1) kq(3) kq(4) kb(4)];
    end
    c = [  1   1   0   0
           0   0   1   1
           1  -1   0   0
           0   0  -1   1];
    r = [  1   1   1   1
          -1  -1   1   1];
   case  3
    kc = [F2 F3 F4 , F6 F7 F8];
    if fort_equiv
     % This seems to be new in .F (2018)
     if FB.ztp(MG.nzaxre) < 0, kc = kc-1; end
    else
     if FB.ztp(1) < 0, kc = kc-1; end
    end
    c = [ op   0  om  om   0  op
          -1  -1   0   0   1   1
           1  -1   0   0   1  -1
          om   0 -op -op   0  om
           0   0  -1   1   0   0
           0   1   0   0   1   0];
    r = c(1:2,:);
    r(2,r(2,:)~=0) = r(2,r(2,:)~=0).*[0.0448 0.0448 0.0556 0.0556]/0.05;
   case  4
    kc = [F1 F2 F3 , F6 F7 F8];
    c = [  0   1   1   1   1   0
          -1  -1   0   0   1   1
           0   1  -1  -1   1   0
          -1   1   0   0  -1   1
           0   0  -1   1   0   0
           1   0   0   0   0   1];
    r = c(1:2,:);
   case  5
    kc = [F2 F3 F4 , F6 F7 F8];
    c = [ op   0  om  om   0  op
           0  -1  -1   1   1   0
           0  -1   1  -1   1   0
         -om   0  op  op   0 -om
          -1   0   0   0   0   1
           0   1   0   0   1   0];
    r = c(1:2,:);
   case  6
    kc = [F1 F3 F4 , F5 F6 F8];
    c = [ op  om   0   0  om  op
           0  -1  -1   1   1   0
         -om  op   0   0  op -om
           0   1  -1   1  -1   0
          -1   0   0   0   0   1
           0   0   1   1   0   0];
    r = c(1:2,:);
   case  7
    kc  = [F3 F4 , F6 F7];
    c = [  1   1   0   0
           0   0   1   1
           1  -1   0   0
           0   0  -1   1];
    r = [  1   1   1   1
          -1  -1   1   1];
   case  8
    kc = [F4 F5 , F6 F7];
    c = [ op  om  om  op , 1/4
          -1  -1   1   1 , 1/4
           1  -1   1  -1 , 1/4
          om -op -op  om , 1/4];
    r = c(1:2,1:end-1);
   case  9
    kc = [F3 F4 F5 , F6 F7 F8];
    c = [  0  op  om  om  op   0
          -1  -1   0   0   1   1
           0  om -op -op  om   0
          -1   1   0   0  -1   1
           0   0  -1   1   0   0
           1   0   0   0   0   1];
    r = c(1:2,:);
   case  10
    kc = [F5 F6 , F7 F8];
    c = [ op  om  om  op , 1/4
          -1  -1   1   1 , 1/4
           1  -1   1  -1 , 1/4
          om -op -op  om , 1/4];
    r = c(1:2,1:end-1);
   case  11
    kc = [F4 F5 , F7 F8];
    c = [ op  om  om  op , 1/4
          -1  -1   1   1 , 1/4
           1  -1   1  -1 , 1/4
          om -op -op  om , 1/4];
    r = c(1:2,1:end-1);
   case 12
    kc = [kb(1) kq(3) kb(2) kb(3) kq(4) kb(4)];
    c = [  1   1   1   1   1   1 ,  1/6
          -1  -1   0   0   1   1 ,  1/4
           1   0  -1  -1   0   1 ,  1/4
          -1   1   0   0  -1   1 ,  1/4
           0   0  -1   1   0   0 ,  1/2
          .5  -1  .5  .5  -1  .5 ,  1/3];
    r = c(1:3,1:end-1);
   case  61
    kc = [F5 F6 , F7 F8];
    c = [  1   1   0   0 , 1/2
           0   0   1   1 , 1/2
           1  -1   0   0 , 1/2
           0   0  -1   1 , 1/2];
    r = [  1   1   1  1
          -1  -1   1  1];
   case  62
    kc = [F3 F4 , F5 F6];
    c = [  1   1   0   0 , 1/2
           0   0   1   1 , 1/2
           1  -1   0   0 , 1/2
           0   0  -1   1 , 1/2];
    r = [  1   1   1  1
          -1  -1   1  1];
   case  63
    kc = [F2 F3 , F5 F6];
    c = [  1   1   0   0 , 1/2
           0   0   1   1 , 1/2
           1  -1   0   0 , 1/2
           0   0  -1   1 , 1/2];
    r = [  1   1   1  1
          -1  -1   1  1];
  end
 end
 nr = size(r,1);
 if 1==0 % future option for automatically determining the orthogonal combinations
  [~,~,tmp] = svd(c(1:nr,1:length(kc)));
  c(nr+1:end,1:length(kc)) = tmp(:,nr+1:end)';
 end
 if size(c,1) == size(c,2), f = 1./diag(c*c');
 else                     , f = c(:,end)     ; c(:,end) = [];
 end
 [Tca,Nca] = deal(eye(16));
 Tca(kc,kc) = c;
 Nca(kc,kc) = c.*f;
 Tra = zeros(nr,16); Tra(:,kc) = r;
 
 %% Quadrupole points (BRBZEQ)
 if midplan == 8
  rq = [repmat(FB.rtp,2,1);repmat(FB.rbp,2,1)]+MG.deltar*[-1;1;-1;1];
  zq = [repmat(FB.ztp,2,1);repmat(FB.zbp,2,1)];
  brq = [MG.brmzero;MG.brpzero;MG.brmzerb;MG.brpzerb];
  bzq = [MG.bzmzero;MG.bzpzero;MG.bzmzerb;MG.bzpzerb];
 elseif fort_equiv
  rq =        FB.rze + MG.deltar*[-1;1];
  zq = repmat(FB.ztp,2,1);
  brq = [MG.brmzero;MG.brpzero];
  bzq = [MG.bzmzero;MG.bzpzero];
 else
  rq =        FB.rtp + MG.deltar*[-1;1];
  zq = repmat(FB.ztp,2,1);
  brq = [MG.brmzero;MG.brpzero];
  bzq = [MG.bzmzero;MG.bzpzero];
 end
 Brq = sum(reshape(G.greenem('br',rq,zq,G.rw,G.zw)*G.Twa,[size(rq) G.na]).*permute(FB.Iea,[3 2 1]),3);
 Bzq = sum(reshape(G.greenem('bz',rq,zq,G.rw,G.zw)*G.Twa,[size(rq) G.na]).*permute(FB.Iea,[3 2 1]),3);
 %% Scenario (SCENARI)
 SC.kbd = numel(find(MG.toft <= 0));
 if fort_equiv
  SC.kru = SC.kbd+10;
 else
  SC.kru = numel(MG.toft);
 end
 SC.kft = SC.kru + MG.numeq-1;
 SC.nrd = (MG.amagic > 0) * MG.numeq;
 SC.kps = SC.kft + SC.nrd + 1;
 SC.nt  = SC.kps + 1;
 
 if fort_equiv
  % This only avoids roundoff differences with FORTRAN version
  tmp = MG.timeeq(1:MG.numeq)*MG.timefac;
  SC.t = MG.toft(1:SC.kru);
  SC.t(end+(1:MG.numeq-1)) = SC.t(end) + (tmp(2:end) - tmp(1));
  tmp = (tmp(end) - tmp(end-1:-1:1))*(MG.amagic/MG.timefac);
  if SC.nrd>0
   SC.t(end+1)            = SC.t(end) + MG.flattop;
   SC.t(end+(1:SC.nrd-1)) = SC.t(end) + tmp(1:SC.nrd-1);
  end
  SC.t(end+1)            = SC.t(end) + 0.015;
  SC.t(end+1)            = SC.t(end) + 0.25;
 else
  tmp  = [MG.flattop fliplr(diff(MG.timeeq(1:MG.numeq)))*MG.amagic];
  SC.t = [MG.toft MG.toft(end)+cumsum([diff(MG.timeeq(1:MG.numeq))*MG.timefac tmp(1:SC.nrd) 0.015 0.25])];
 end
 
 tmp = (MG.toft(SC.kbd+1:SC.kru)-MG.toft(SC.kbd))/(MG.toft(SC.kru)-MG.toft(SC.kbd));
 %   %pre            %breakdown            %ramp-up           %ramp-down     %post
 k = [ ones(1,SC.kbd) ones(1,SC.kru-SC.kbd) 2:MG.numeq         SC.nrd:-1:1     ones(1,2)];
 f = [zeros(1,SC.kbd) tmp                   ones(1,MG.numeq-1) ones(1,SC.nrd) zeros(1,2)];
 a = 1.63;
 % MGA.F uses sqrt(f) for Ip, f.^2 seems more appropriate - only important at BD
 if fort_equiv, ap = 0.5; fq = 0;
 else,          ap = 2.0; fq = 1;end
 SC.Ip  = f.^ap.*FB.Ip (  k);
 SC.Ia  = f.^2 .*FB.Iea(:,k);
 SC.Ff  = f.^2 .*FB.Ff (:,k);
 SC.Bm  = f.^2 .*FB.Bm (:,k);
 % BD quadrupole maintained to comply with .f "error" (replace fq with 1 to remove this)
 SC.Brq = f.^a.*Brq   (:,k) + (1-fq*f.^a).*brq; if ~fort_equiv, SC.Brq(:,[1 end-1 end]) = 0;end
 SC.Bzq = f.^a.*Bzq   (:,k) + (1-fq*f.^a).*bzq; if ~fort_equiv, SC.Bzq(:,[1 end-1 end]) = 0;end
 SC.Bt  = MG.bzero(k); SC.Bt([1 end]) = 0;
 if midplan == 8, SC.rp = [FB.rtp(k) ; FB.rbp(k)];
  SC.zp = [FB.ztp(k) ; FB.zbp(k)];
 elseif fort_equiv
  SC.rp =  FB.rze(k);
  SC.zp =  FB.ztp(k);
 else
  SC.rp =  FB.rtp(k);
  SC.zp =  FB.ztp(k);
 end
 SC = orderfields(SC);
 
 dioh = [MG.diohdt(1:SC.kru)*1000 zeros(1,SC.nt-SC.kru)];
 if fort_equiv
  vind = diff((G.greenem('self',FB.rze,FB.ztp,sqrt([FB.Ap./FB.Kp ; FB.Ap.*FB.Kp]),20)+2e-7*pi*FB.li.*FB.rze).*FB.Ip+FB.Fa)./diff(MG.timeeq(1:MG.numeq));
 else
  vind = diff((G.greenem('self',FB.rtp,FB.ztp,sqrt([FB.Ap./FB.Kp ; FB.Ap.*FB.Kp])   )+2e-7*pi*FB.li.*FB.rtp).*FB.Ip+FB.Fa)./diff(MG.timeeq(1:MG.numeq));
 end
 k = find(SC.t(SC.kru+1:SC.kft-1) > ifc(MG.ierat == 2,0.050,-Inf));
 dioh(k + SC.kru) = -15504*(vind(k+1)/MG.timefac+MG.relovo);
 dioh(    SC.kft) = -15504*MG.relovo;
 if SC.nrd ~= 0
  dioh(SC.kft+1:SC.kps-2) =  15504*(vind(end:-1:1)  /MG.amagic-MG.relovo);
 end
 dioh(SC.kps-1) = dioh(SC.kps-2);
 SC.Ia(OH,:) = repmat([0 cumsum(dioh(1:end-2).*diff(SC.t(1:end-1))) 0],2,1);
 if fort_equiv
  dioh(end-1) = -SC.Ia(OH(1),end-1)/diff(SC.t(end-1:end));
 end
 
 % E bias
 if midplan == 8, Ebias = 0;
 else,            Ebias = 100*MG.iohfb;  end
 k = setdiff(E1:E8,kq);
 if fort_equiv
  Ebias = [repmat(Ebias,1,SC.nt-2)   0 0]; % Should be only before BD?
 else
  Ebias = [0 repmat(Ebias,1,SC.nt-3) 0 0]; % Should be only before BD?
 end
 SC.Ia(k,:) = SC.Ia(k,:)+Ebias;
 % Bz ramp
 if MG.iohfb == 1, Bzramp = 0;
 else
  if midplan == 8, Bzramp = +0.015;    % Bz ramp modified by JS for doublets (used to be uc1 for iohfb == 1)
  else,            Bzramp = -0.01; end
 end
 if isfield(MG,'bzramp'), Bzramp = MG.bzramp;end
 SC.Bzq(:,2:3) = SC.Bzq(:,2:3) + Bzramp;
 
 %% Resampling (REFER)
 if fort_equiv
  TV.t = SC.t(1);
  while TV.t(end) < SC.t(2)-0.02       , TV.t(end+1) = TV.t(end)+0.002;end
  while TV.t(end) < SC.t(SC.kru+4)     , TV.t(end+1) = TV.t(end)+0.001;end
  while TV.t(end) < SC.t(SC.kft)+0.02  , TV.t(end+1) = TV.t(end)+0.002;end
  while TV.t(end) < SC.t(SC.kft+1)-0.02, TV.t(end+1) = TV.t(end)+0.010;end
  while TV.t(end) < SC.t(end)+0.002    , TV.t(end+1) = TV.t(end)+0.002;end
  TV.t(MG.ntmax+1:end) = [];
  interp = @(x,y,z) interp1(x,y,z,'linear','extrap');
 else
  TV.dt      = 0.002;
  TV.t       = [-fliplr(TV.dt:TV.dt:-min(SC.t)) 0 TV.dt:TV.dt:max(SC.t)];
  interp = @interp1q;
 end
 TV.dt      = diff(TV.t);
 TV.nt      = numel(TV.t);
 TV.Ip      = interp(SC.t',SC.Ip' ,TV.t')';
 TV.Bt      = interp(SC.t',SC.Bt' ,TV.t')';
 TV.rp      = interp(SC.t',SC.rp' ,TV.t')';
 TV.zp      = interp(SC.t',SC.zp' ,TV.t')';
 TV.Bzq     = interp(SC.t',SC.Bzq',TV.t')';
 TV.Brq     = interp(SC.t',SC.Brq',TV.t')';
 TV.Ff      = interp(SC.t',SC.Ff' ,TV.t')';
 TV.Bm      = interp(SC.t',SC.Bm' ,TV.t')';
 TV.Ia      = interp(SC.t',SC.Ia' ,TV.t')';
 TV.Iv      = zeros(G.nv,TV.nt);
 TV         = orderfields(TV);
 
 %% Vessel and compensation coil currents (CURRENT)
 % OH backoff
 TV.Ia(kb,:) = TV.Ia(kb,:) + Baa*TV.Ia(OH1,:);
 
 % Time varying Green's functions
 rq = [TV.rp-MG.deltar;TV.rp+MG.deltar]; zq = repmat(TV.zp,2,1);
 % for doublets
 if midplan ==8
  rq = [rq(1,:);rq(3,:);rq(2,:);rq(4,:)]; zq = [zq(1,:);zq(3,:);zq(2,:);zq(4,:)];
  Map1 = permute(reshape(G.greenem('m' ,TV.rp(1,:),TV.zp(1,:),G.rw,G.zw)*G.Twa,[size(TV.rp(1,:)),G.na]),[3 1 2]);
  Map2 = permute(reshape(G.greenem('m' ,TV.rp(2,:),TV.zp(2,:),G.rw,G.zw)*G.Twa,[size(TV.rp(2,:)),G.na]),[3 1 2]);
  Map  = Map1+Map2;
  Mvp1 = permute(reshape(G.greenem('m' ,TV.rp(1,:),TV.zp(1,:),G.ru,G.zu)*G.Tuv,[size(TV.rp(1,:)),G.nv]),[3 1 2]);
  Mvp2 = permute(reshape(G.greenem('m' ,TV.rp(2,:),TV.zp(2,:),G.ru,G.zu)*G.Tuv,[size(TV.rp(2,:)),G.nv]),[3 1 2]);
 else
  if fort_equiv
   [wpla,hpla] = deal(0.1*ones(size(TV.rp)));
   Map = permute(reshape(G.greenem('mut',TV.rp,TV.zp,wpla,hpla,G.rw,G.zw,G.aw,G.bw,8)*G.Twa,[size(TV.rp),G.na]),[3 1 2]);
   Mvp = permute(reshape(G.greenem('mut',TV.rp,TV.zp,wpla,hpla,G.ru,G.zu,G.au,G.bu,8)*G.Tuv,[size(TV.rp),G.nv]),[3 1 2]);
  else
   Map = permute(reshape(G.greenem('mut',TV.rp,TV.zp,G.rw,G.zw)*G.Twa,[size(TV.rp),G.na]),[3 1 2]);
   Mvp = permute(reshape(G.greenem('mut',TV.rp,TV.zp,G.ru,G.zu)*G.Tuv,[size(TV.rp),G.nv]),[3 1 2]);
  end
 end
 Brqa = permute(reshape(G.greenem('br',   rq,   zq,G.rw,G.zw)*G.Twa,[size(   rq),G.na]),[1 3 2]);
 Bzqa = permute(reshape(G.greenem('bz',   rq,   zq,G.rw,G.zw)*G.Twa,[size(   rq),G.na]),[1 3 2]);
 Brqv = permute(reshape(G.greenem('br',   rq,   zq,G.ru,G.zu)*G.Tuv,[size(   rq),G.nv]),[1 3 2]);
 Bzqv = permute(reshape(G.greenem('bz',   rq,   zq,G.ru,G.zu)*G.Tuv,[size(   rq),G.nv]),[1 3 2]);
 % fort_equiv not needed, these values aren't used.
 Brqt =         reshape(G.greenem('br',   rq,   zq,G.rt,G.zt)*G.Ttt,[size(   rq), 1]);
 Bzqt =         reshape(G.greenem('bz',   rq,   zq,G.rt,G.zt)*G.Ttt,[size(   rq), 1]);
 
 % TFC backoff
 switch midplan
  case    3, c = 1;                         cr = -420e-6; cz =  150e-6;
  case    1, c = 1;                         cr = -690e-6; cz = -210e-6;
  case   37, c = 1;                         cr = -910e-6; cz = -860e-6;
  case    8, c = ([1+iplatop 1-iplatop])/2; cr =    0   ; cz =    0   ;
  otherwise, c = 1;                         cr =    0   ; cz =    0   ;
 end
 cr = cr*MG.if36fb/1.43; cz = cz*MG.if36fb/1.43;
 if midplan == 8,   Brqt    =      Brqt*65e3;                Bzqt    =      Bzqt*65e3;
 elseif fort_equiv, Brqt(:) = cr*MG.bzero(1);                Bzqt(:) = cz*MG.bzero(1);
 else,              Brqt    = repmat(cr*TV.Bt,size(rq,1),1); Bzqt    = repmat(cz*TV.Bt,size(rq,1),1);
 end
 % for doublets
 if midplan == 8
  Uvp1(:,2:TV.nt) = c(1)*diff(Mvp1.*reshape(TV.Ip,[1,size(TV.Ip)]),1,3);
  Uvp2(:,2:TV.nt) = c(1)*diff(Mvp2.*reshape(TV.Ip,[1,size(TV.Ip)]),1,3);
  Uvp (:,2:TV.nt) = Uvp1(:,2:TV.nt)+Uvp2(:,2:TV.nt);
 else
  Uvp (:,2:TV.nt) = c*diff(Mvp.*reshape(TV.Ip,[1,size(TV.Ip)]),1,3);
 end
 
 if nbackoff==16 && MG.iohfb == -1
  % Empiric vertical field to shift nulls in IP<0 scenarios (JS)
  switch midplan
   case 37
    kim=[F6 F7 F8];
    Iim=[100 100 100];
   case 3
    kim=[F1 F2 F3];
    Iim=[100 100 100];
   case 1
    kim=[F1 F2 F7 F8];
    Iim=[100 100 100 100];
   case 5
    kim = [F1 F2 F3];
    Iim = [100 100 100];
   case 35
    kim = [F6 F7 F8];
    Iim = [100 100 100];
  end
  TV.Ia(kim,:)=TV.Ia(kim,:)+repmat(Iim.',[1,TV.nt]);
 end
 
 %TFC back-off not done using quadrupole coils
 if nbackoff ==16 || midplan==8
  Brqt = zeros(1,TV.nt);
  Bzqt = zeros(1,TV.nt);
 end
 
 % try to decrease the quadrupole field value from 0 to 6ms
 if nbackoff==16
  TV.Bzq(:,142:144) = TV.Bzq(:,141)*[0.75,0.5,0.25];
  TV.Brq(:,142:144) = TV.Brq(:,141)*[0.75,0.5,0.25];
 end

 if fort_equiv, kt_end = TV.nt;
 else,          kt_end = find(TV.t > SC.t(SC.kps),1);end
 if fort_equiv
  TV.Ia(EF,1) = 0;
  TV.Ia(kq,:) = 0; % FORTRAN uses full currents as variables
  for kt = 2:kt_end
   % FORTRAN uses previous total vessel current
   TV.Ia(kb,kt) = TV.Ia(kb,kt) + Bav*sum(TV.Iv(:,kt-1));
   tmpA = [diag(G.Rv)*TV.dt(kt-1)+G.Mvv    G.Mva(:,kq);
           Brqv(:,:,kt)                  Brqa(:,kq,kt);
           Bzqv(:,:,kt)                  Bzqa(:,kq,kt)];
   tmpb = [G.Mvv*TV.Iv(:,kt-1)-Uvp(:,kt)-G.Mva*(TV.Ia(:,kt)-TV.Ia(:,kt-1));
            TV.Brq(:,kt)-Brqt(:,kt)-Brqa(:,:,kt)*TV.Ia(:,kt);
            TV.Bzq(:,kt)-Bzqt(:,kt)-Bzqa(:,:,kt)*TV.Ia(:,kt)];
   tmp = single(double(tmpA)\double(tmpb));
   TV.Iv( :,kt) =                  tmp(     1:G.nv);
   TV.Ia(kq,kt) =                  tmp(G.nv+1:end );
  end
 else
  for kt = 2:kt_end
   tmp = [diag(G.Rv)*TV.dt(kt-1)+G.Mvv+G.Mva(:,kb)*repmat(Bav,1,G.nv)   G.Mva(:,kq);
          Brqv(:,:,kt)+Brqa(:,kb,kt)*repmat(Bav,1,G.nv)               Brqa(:,kq,kt);
          Bzqv(:,:,kt)+Bzqa(:,kb,kt)*repmat(Bav,1,G.nv)               Bzqa(:,kq,kt)] \ ...
         [G.Mvv*TV.Iv(:,kt-1)-Uvp(:,kt)-G.Mva*(TV.Ia(:,kt)-TV.Ia(:,kt-1));
          TV.Brq(:,kt)-Brqt(:,kt)-Brqa(:,:,kt)*TV.Ia(:,kt);
          TV.Bzq(:,kt)-Bzqt(:,kt)-Bzqa(:,:,kt)*TV.Ia(:,kt)];
   TV.Iv( :,kt) =                tmp(     1:G.nv);
   TV.Ia(kq,kt) = TV.Ia(kq,kt) + tmp(G.nv+1:end );
   TV.Ia(kb,kt) = TV.Ia(kb,kt) + Bav*sum(TV.Iv(:,kt));
  end
 end
 if ~fort_equiv
  TV.Ia(EF,kt_end:TV.nt) = TV.Ia(EF,kt)*linspace(1,0,TV.nt-kt_end+1);
 end
 if nbackoff==16 || midplan==8
  TV.Ia(kb,:)=TV.Ia(kb,:)+repmat(MG.if36fb*Ba_tor.',[1,TV.nt]); 
 end
 % Compute Iv from Ia evolution (assumes Ip does not affect Iv)
 for kt = kt_end+1:TV.nt
  TV.Iv(:,kt) = (diag(G.Rv)*TV.dt(kt-1)+G.Mvv) \ (G.Mvv*TV.Iv(:,kt-1)-G.Mva*(TV.Ia(:,kt)-TV.Ia(:,kt-1)));
 end
 
 % Inductive coil voltage given Ia, Iv, Ip evolution
 TV.Ua(:,2:TV.nt) = diff(G.Maa*TV.Ia + G.Mva'*TV.Iv + squeeze(c(1)*Map.*reshape(TV.Ip,[1,size(TV.Ip)])),1,2)./TV.dt;
 
 %% Shape references (WAVGEN1)
 f     = min(max(TV.t/MG.toft(SC.kru),0),1);
 TV.Ff = (1-f).*(G.Mfa*TV.Ia+G.Mfv*TV.Iv) + f.*TV.Ff;
 TV.Bm = (1-f).*(G.Bma*TV.Ia+G.Bmv*TV.Iv) + f.*TV.Bm;
 Ipon  = TV.t > 0;
 
 % Plasma current reference
 Ipref = TV.Ip*MG.curfac + MG.aipdel;
 % Radial reference
 f = (TV.t-SC.t(SC.kbd))/(SC.t(SC.kps-1)-SC.t(SC.kbd));
 Rcor = -MG.psifac*MG.dpszero + 1e-12*MG.psifac*MG.rshift3*TV.Ip.^2 + 1e-6*MG.psifac*TV.Ip.*(MG.rshift1+diff([MG.rshift1 MG.rshift2])*f);
 if fort_equiv, f = min(max(TV.t/MG.toft(SC.kru),0),1);end
 f =  (MG.dpsfac - 1)*f + 1;
 if midplan == 8
  Rref  = deal(Rcor - f.*(Radm(1,:)*TV.Bm + Radf(1,:)*TV.Ff));
  Rbref = deal(Rcor - f.*(Radm(2,:)*TV.Bm + Radf(2,:)*TV.Ff));
 else
  [Rref,Rbref] = deal(Rcor - f.*(Radm*TV.Bm + Radf*TV.Ff));
 end
%  % Radial reference
%  f = (TV.t-SC.t(SC.kbd))/(SC.t(SC.kps-1)-SC.t(SC.kbd));
%  Rcor = -MG.psifac*MG.dpszero + 1e-12*MG.psifac*MG.rshift3*TV.Ip.^2 + 1e-6*MG.psifac*TV.Ip.*(MG.rshift1+diff([MG.rshift1 MG.rshift2])*f);
%  f =  (MG.dpsfac - 1)*f + 1;
%  [Rref,Rbref] = deal(Rcor - f.*(Radm*TV.Bm + Radf*TV.Ff));

 % choosing the modififed radial ref. (MG.radsel=1) indep.of Ipref
 if MG.radsel == 1
  Rref  = Rref./TV.Ip;   Rref(isinf(Rref))= 0;
  Rbref = Rbref./TV.Ip;  Rbref(isinf(Rbref))= 0;
 end
 
 % Vertical reference
 [zbref,zref2,zref3,kref] = deal(zeros(1,TV.nt));
 zcor = [MG.aipipz MG.zeecorr(1)].*[5000 1e-12]*[Ipon; TV.Ip.^2];
 if midplan == 8
  if MG.ikriz ==1 %choose b/w the z-references constructed using new vertical observer or F.Hofmann's observer 
   % z references created from the new vertical observer for doublet
   zref1 = +2*MG.rampt*Ipzf_t*TV.Ff + 2*(1-MG.rampt)*Ipzm_t*TV.Bm -Ipza_t*TV.Ia-3600;% -3600;% +MG.delipz +zcor; % top plasma
   zbref = +2*MG.rampt*Ipzf_b*TV.Ff + 2*(1-MG.rampt)*Ipzm_b*TV.Bm -Ipza_b*TV.Ia+3600;% +3600;% -MG.delipz +zcor; % bottom plasma
  else
   % why not -2*MG.rampt and 2*(1-MG.rampt) as in observer
   zref1 = MG.rampt*IpDz_xf_t*TV.Ff + (1-MG.rampt)*IpDz_xm_t*TV.Bm + 1*MG.delipz + -1*zcor+10000;
   zbref = MG.rampt*IpDz_xf_b*TV.Ff + (1-MG.rampt)*IpDz_xm_b*TV.Bm -   MG.delipz + -1*zcor+5000;
   % zref1 = -MG.rampt*IpDz_xf_t*TV.Ff - (1-MG.rampt)*IpDz_xm_t*TV.Bm - MG.delipz + zcor;
   % zbref = -MG.rampt*IpDz_xf_b*TV.Ff - (1-MG.rampt)*IpDz_xm_b*TV.Bm + MG.delipz + zcor;
  end
 else
  % Vessel and coil correction
  zref = (TV.zp-zp0).*TV.Ip + MG.delipz + zcor;
  tmp = zeros(4,TV.nt);
  tmp(1,:) = zref                                + Ipza      (  OH).*(  MG.ohcorr*[-1 1])*TV.Ia(OH,:).*Ipon + Ipzv      *TV.Iv;
  tmp(2,:) = zref + Ipza      (  EF)*TV.Ia(EF,:) + Ipza      (  OH).*(1+MG.ohcorr*[-1 1])*TV.Ia(OH,:).*Ipon + Ipzv      *TV.Iv;
  tmp(3,:) = zref + Ipzma     (  EF)*TV.Ia(EF,:) + Ipzma     (  OH).*(1+MG.ohcorr*[-1 1])*TV.Ia(OH,:).*Ipon + Ipzmv     *TV.Iv;
  tmp(4,:) = zref + Ipzm*G.Bma(:,EF)*TV.Ia(EF,:) + Ipzm*G.Bma(:,OH).*(1+MG.ohcorr*[-1 1])*TV.Ia(OH,:).*Ipon + Ipzm*G.Bmv*TV.Iv;
  % zref1 is the Ipz reference for slow proportional control (Ref 20)
  if ismember(MG.inova,[3 4])          , zref1 = tmp(1,:);
  else                                 , zref1 = tmp(2,:);
  end
  % zref2 is the Ipz reference for slow and/or fast derivative control (Ref 23)
  if fort_equiv
   % NOTE: This could perhaps be simplified since zref2 is not always in use
   if (MG.ikriz == 2 || MG.nfast == 0) && MG.inova == 4, zref2 = tmp(2,:);
   elseif (MG.ikriz == 2 && MG.nfast == -2)            , zref2 = tmp(3,:);
   elseif (MG.ikriz == 2 && ismember(MG.nfast,[-3,-5])), zref2 = tmp(4,:);
   else                                                , zref2 = zref;
   end
  else
   if MG.inova == 4                     , zref2 = tmp(2,:);
   elseif ismember(MG.nfast,[-3 -5])    , zref2 = tmp(4,:);
   elseif MG.ikriz == 2 || MG.nfast == 0, zref2 = tmp(3,:);
   else                                 , zref2 = zref;
   end
  end
  % zref3 is the Ipz reference for fast derivative control (Ref 3) - must be inova = 4
  if fort_equiv
   % NOTE: This could perhaps be simplified since zref3 is not always in use
   if MG.ikriz == 2 && MG.inova == 4
    if MG.nfast == -2                   , zref3 = tmp(3,:);
    elseif ismember(MG.nfast,[-3,-5])   , zref3 = tmp(4,:);
    else                                , zref3 = zref;
    end
   else                                 , zref3 = zref;
   end
  else
   if ismember(MG.nfast,[-3 -5])        , zref3 = tmp(4,:);
   elseif MG.ikriz == 2 || MG.nfast == 0, zref3 = tmp(3,:);
   else                                 , zref3 = zref;
   end
  end
  
  % Wobbling
  if MG.zshift(1)
   if fort_equiv && ~MG.zshift(2), MG.zshift(2) = SC.t(SC.kft  ); end
   if fort_equiv && ~MG.zshift(3), MG.zshift(3) = SC.t(SC.kft+1); end
   k = TV.t > MG.zshift(2) & TV.t < MG.zshift(3);
   if fort_equiv && MG.ipripz > 0
    Zcor = sin(2*pi*abs(MG.ipripz)*(TV.t-SC.t(SC.kft)));
   else
    Zcor = sin(2*pi*abs(MG.ipripz)*(TV.t-MG.zshift(2)));
   end
   if MG.ipripz > 0, Zcor = sign(Zcor); end
   if fort_equiv
    zref1(k) = zref1(k) + MG.zshift(1)*Zcor(k).*TV.Ip(k).^2*1e-12;
    zref2(k) = zref2(k) + MG.zshift(1)*Zcor(k).*TV.Ip(k).^2*1e-12;
    zref3(k) = zref3(k) + MG.zshift(1)*Zcor(k).*TV.Ip(k).^2*1e-12;
   else
    zref1(k) = zref1(k) + MG.zshift(1)*Zcor(k);
   end
  end
 end
 
 % Kappa reference
 if MG.ifour < 20
  if fort_equiv
   TV.Ff = interp(SC.t',SC.Ff' ,TV.t')';
   TV.Bm = interp(SC.t',SC.Bm' ,TV.t')';
  end
  kref = Ipkf*TV.Ff + Ipkm*TV.Bm;
 else
  switch rem(MG.ifour,10)
   case 4, kref = abs(TV.zp-zp0)   .*TV.Ip;
   case 5, kref =    (TV.zp-zp0).^2.*TV.Ip;
  end
 end
 switch floor(MG.ifour/10)
  case 1, kref = kref - Ipka*TV.Ia;
  case 2, kref = kref + Ipka*TV.Ia;
 end
%   kref = kref + MG.gainext(2) + 1e-6*MG.gainext(3)*abs(TV.Ip);
  
 %% References (WAVGEN2)
 % filter
 % f = [0 0.03 0.12 0.24 0.37 0.5 0.63 0.76 0.88 0.97 ones(1,TV.nt-10)];
 % [B,A] = butter(2,150*2*TV.dt); RF.Ua = filtfilt(B,A,(TV.Ua.*f)')';
 % [B,A] = butter(1,450*2*TV.dt); RF.Ia = filtfilt(B,A,(TV.Ia.*f)')';
 %                                zref1 = filtfilt(B,A,zref1      );
 %                                zref2 = filtfilt(B,A,zref2      );
 %                                zref3 = filtfilt(B,A,zref3      );
 %                                Rref  = filtfilt(B,A,Rref       );
 %                                Rbref = filtfilt(B,A,Rbref      );
 %                                kref  = filtfilt(B,A,kref       );
 if fort_equiv
  % startup
  f = interp1(SC.t(1)+2e-3*(0:12),[0,0,0.03,0.12,0.24,0.37 0.5 0.63 0.76 0.88 0.97 1 1],TV.t,'linear','extrap');
  TV.Ua          = f.*TV.Ua;
  TV.Ia(E1:F8,:) = f.*TV.Ia(E1:F8,:);
  zref1          = f.*zref1;
  zref2          = f.*zref2;
  zref3          = f.*zref3;
  zbref          = f.*zbref;
  Rref           = f.*Rref;
  Rbref          = f.*Rbref;
  kref           = f.*kref;
  %
  % filter
  RF.t = NaN(1,TV.nt);
  RF.Ua = TV.Ua;
  RF.Ia = TV.Ia;
  % First phase (pre)
  ista = 4;
  iend = 2*find(TV.t(2:2:end-1) < SC.t(2)-0.02,1,'last');
  ksta = 2;
  kend = iend/2;
  RF.t (     1:kend) =  TV.t(   ista-2:2:iend);
  RF.Ua(:,ksta:kend) = (TV.Ua(:,ista-1:2:iend) + 2*TV.Ua(:,ista:2:iend) + TV.Ua(:,ista+1:2:iend+1))/4;
  RF.Ia(:,ksta:kend) = (TV.Ia(:,ista-1:2:iend) + 2*TV.Ia(:,ista:2:iend) + TV.Ia(:,ista+1:2:iend+1))/4;
  zref1(:,ksta:kend) = (zref1(:,ista-1:2:iend) + 2*zref1(:,ista:2:iend) + zref1(:,ista+1:2:iend+1))/4;
  zref2(:,ksta:kend) = (zref2(:,ista-1:2:iend) + 2*zref2(:,ista:2:iend) + zref2(:,ista+1:2:iend+1))/4;
  zref3(:,ksta:kend) = (zref3(:,ista-1:2:iend) + 2*zref3(:,ista:2:iend) + zref3(:,ista+1:2:iend+1))/4;
  zbref(:,ksta:kend) = (zbref(:,ista-1:2:iend) + 2*zbref(:,ista:2:iend) + zbref(:,ista+1:2:iend+1))/4;
  Rref (:,ksta:kend) = (Rref (:,ista-1:2:iend) + 2*Rref (:,ista:2:iend) + Rref (:,ista+1:2:iend+1))/4;
  Rbref(:,ksta:kend) = (Rbref(:,ista-1:2:iend) + 2*Rbref(:,ista:2:iend) + Rbref(:,ista+1:2:iend+1))/4;
  kref (:,ksta:kend) = (kref (:,ista-1:2:iend) + 2*kref (:,ista:2:iend) + kref (:,ista+1:2:iend+1))/4;
  Ipref(:,ksta:kend) = (Ipref(:,ista-1:2:iend) + 2*Ipref(:,ista:2:iend) + Ipref(:,ista+1:2:iend+1))/4;
  % Second phase (breakdown)
  ista = iend+2;
  iend = find(TV.t(1:end-2) < SC.t(SC.kru+4),1,'last');
  ksta = kend+1;
  kend = ksta+(iend - ista);
  RF.t (ksta:kend) = TV.t(ista:iend);
  RF.Ua(:,ksta:kend) = (TV.Ua(:,ista-2:iend-2) + 2*TV.Ua(:,ista-1:iend-1) + 2*TV.Ua(:,ista:iend) + 2*TV.Ua(:,ista+1:iend+1) + TV.Ua(:,ista+2:iend+2))/8;
  RF.Ia(:,ksta:kend) = (TV.Ia(:,ista-1:iend-1)                            + 2*TV.Ia(:,ista:iend)                            + TV.Ia(:,ista+1:iend+1))/4;
  zref1(:,ksta:kend) = (zref1(:,ista-1:iend-1)                            + 2*zref1(:,ista:iend)                            + zref1(:,ista+1:iend+1))/4;
  zref2(:,ksta:kend) = (zref2(:,ista-1:iend-1)                            + 2*zref2(:,ista:iend)                            + zref2(:,ista+1:iend+1))/4;
  zref3(:,ksta:kend) = (zref3(:,ista-1:iend-1)                            + 2*zref3(:,ista:iend)                            + zref3(:,ista+1:iend+1))/4;
  zbref(:,ksta:kend) = (zbref(:,ista-1:iend-1)                            + 2*zbref(:,ista:iend)                            + zbref(:,ista+1:iend+1))/4;
  Rref (:,ksta:kend) = (Rref (:,ista-1:iend-1)                            + 2*Rref (:,ista:iend)                            + Rref (:,ista+1:iend+1))/4;
  Rbref(:,ksta:kend) = (Rbref(:,ista-1:iend-1)                            + 2*Rbref(:,ista:iend)                            + Rbref(:,ista+1:iend+1))/4;
  kref (:,ksta:kend) = (kref (:,ista-1:iend-1)                            + 2*kref (:,ista:iend)                            + kref (:,ista+1:iend+1))/4;
  Ipref(:,ksta:kend) = (Ipref(:,ista-1:iend-1)                            + 2*Ipref(:,ista:iend)                            + Ipref(:,ista+1:iend+1))/4;
  % Third phase (to the end)
  ista = iend+1;
  iend = TV.nt-2;
  ksta = kend+1;
  kend = ksta-1+numel(ista:3:iend);
  RF.t (ksta:kend) = TV.t(ista:3:iend);
  RF.Ua(:,ksta:kend) = (TV.Ua(:,ista-2:3:iend-2) + 2*TV.Ua(:,ista-1:3:iend-1) + 3*TV.Ua(:,ista:3:iend) + 2*TV.Ua(:,ista+1:3:iend+1) + TV.Ua(:,ista+2:3:iend+2))/9;
  RF.Ia(:,ksta:kend) = (TV.Ia(:,ista-2:3:iend-2) + 2*TV.Ia(:,ista-1:3:iend-1) + 3*TV.Ia(:,ista:3:iend) + 2*TV.Ia(:,ista+1:3:iend+1) + TV.Ia(:,ista+2:3:iend+2))/9;
  if MG.zshift(1)
   zref1(:,ksta:kend) = zref1(:,ista:3:iend);
  else
   zref1(:,ksta:kend) = (zref1(:,ista-2:3:iend-2) + 2*zref1(:,ista-1:3:iend-1) + 3*zref1(:,ista:3:iend) + 2*zref1(:,ista+1:3:iend+1) + zref1(:,ista+2:3:iend+2))/9;
  end
  zref2(:,ksta:kend) = (zref2(:,ista-2:3:iend-2) + 2*zref2(:,ista-1:3:iend-1) + 3*zref2(:,ista:3:iend) + 2*zref2(:,ista+1:3:iend+1) + zref2(:,ista+2:3:iend+2))/9;
  zref3(:,ksta:kend) = (zref3(:,ista-2:3:iend-2) + 2*zref3(:,ista-1:3:iend-1) + 3*zref3(:,ista:3:iend) + 2*zref3(:,ista+1:3:iend+1) + zref3(:,ista+2:3:iend+2))/9;
  zbref(:,ksta:kend) = (zbref(:,ista-2:3:iend-2) + 2*zbref(:,ista-1:3:iend-1) + 3*zbref(:,ista:3:iend) + 2*zbref(:,ista+1:3:iend+1) + zbref(:,ista+2:3:iend+2))/9;
  Rref (:,ksta:kend) = (Rref (:,ista-2:3:iend-2) + 2*Rref (:,ista-1:3:iend-1) + 3*Rref (:,ista:3:iend) + 2*Rref (:,ista+1:3:iend+1) + Rref (:,ista+2:3:iend+2))/9;
  Rbref(:,ksta:kend) = (Rbref(:,ista-2:3:iend-2) + 2*Rbref(:,ista-1:3:iend-1) + 3*Rbref(:,ista:3:iend) + 2*Rbref(:,ista+1:3:iend+1) + Rbref(:,ista+2:3:iend+2))/9;
  kref (:,ksta:kend) = (kref (:,ista-2:3:iend-2) + 2*kref (:,ista-1:3:iend-1) + 3*kref (:,ista:3:iend) + 2*kref (:,ista+1:3:iend+1) + kref (:,ista+2:3:iend+2))/9;
  Ipref(:,ksta:kend) = (Ipref(:,ista-2:3:iend-2) + 2*Ipref(:,ista-1:3:iend-1) + 3*Ipref(:,ista:3:iend) + 2*Ipref(:,ista+1:3:iend+1) + Ipref(:,ista+2:3:iend+2))/9;
  % Stop traces
  RF.t (  kend+1:end) = [];
  RF.Ua(:,kend+1:end) = [];
  RF.Ia(:,kend+1:end) = [];
  zref1(:,kend+1:end) = [];
  zref2(:,kend+1:end) = [];
  zref3(:,kend+1:end) = [];
  zbref(:,kend+1:end) = [];
  Rref (:,kend+1:end) = [];
  Rbref(:,kend+1:end) = [];
  kref (:,kend+1:end) = [];
  Ipref(:,kend+1:end) = [];
 else
  RF.t  = TV.t;
  RF.Ua = TV.Ua; % RF.Ua = f.*TV.Ua;
  RF.Ia = TV.Ia; % RF.Ia = f.*TV.Ia;
 end
 
 Ref = zeros(numel(RF.t),24);           Ref(:,EF) = (Nca * RF.Ia(EF,:))';
 % doublet EFWAVE reference
 if midplan==8,                         Ref(:,EF) = (Tca*RF.Ia(EF,:))'; end
                                        Ref(:,18) = Rref;
                                        Ref(:,19) = Ipref;
                                        Ref(:,20) = zref1;
 if      midplan == 8,                  Ref(:,22) = zbref;
                                        Ref(:,23) = Rbref;
 elseif MG.ikriz == 2 || MG.nfast == 0, Ref(:,23) = zref2;
  if    MG.inova == 4,                  Ref(:, 3) = zref3; end
 end
 
 %% M matrix (MMATRIX)
 nM = 20; mM = 40; nsw = 10;
 Mmx = zeros(nM,mM);
 
 % Gains below give maximum load_mat elements of 0.495
 % i[e,f,oh]_gain(2) = 0.495*[650.7,1250.1,1398.6]/(G.Maa([1,9,17],[1,9,17])*[2020,2020,6868]/50
 gM = [repmat(MG.egain,1,8) repmat(MG.fgain,1,8) MG.ohgain MG.ohgain]*MG.emgain;
 if MG.iscramb>0
  switch MG.mvloop, case {1   11    21},  gM( 1     ) = 0;
                    case {  2    12   },  gM([1 2  ]) = 0;
                    case 6,               gM([  2 3]) = 0;     % mvloop = 6, elongation control by ECRH
                                          gM([4 8  ]) = 0;
  end
  if MG.ierat == 4,                       gM( 1     ) = 0; end
  if mod(MG.mvloop, 10) ==  5,            gM( 1     ) = 0; end % mvloop = 5, plasma current control by ECRH
  if mod(MG.mvloop,100) == 25,            gM(   2   ) = 0; end
  
  if MG.inova == 4,                       gM(     3 ) = 0; end
 end
  
 Mmx(1:G.na,1:G.na)          = G.Maa.*gM;
 Mmx(1:G.na,G.na+5:2*G.na+4) = (~inores)*diag(G.Ra);
 if MG.nfast == -5, MG.fastm = 0.85*MG.fastm; end % account for difference in observer, Nuno Cruz calculates 85% reduction need
 if MG.nfast ~=  0, Mmx(20,22) = MG.fastm; end %FPS
 Mmx(19,21) = 6; %GAS
 
 Mmx = repmat(Mmx,[1,1,nsw]);
  
 %% A matrix (AMATRIX)
 mA = 24; nA = 120; Amx = zeros(mA,nA);
 kAf = 1:G.nf; kAm = G.nf+(1:G.nm); kAa = G.nf+G.nm+(1:G.na); kAvloop = 119; kADOH = 115; kAne = 116;
 
 if midplan ~= 8
  % Radial observer
  Amx(18,kAf) = -Radf;
  Amx(18,kAm) = -Radm;
  
  % Vertical observer
  if MG.nfast == 0,    Amx(20,kAf)          =  Ipzf;
                       Amx(20,kAm)          =  Ipzm;
  else,                Amx(20,kAf)          =  Ipzf*(  2*MG.rampt);
                       Amx(20,kAm)          =  Ipzm*(2-2*MG.rampt);
                       Amx(20,kAm([14 26])) =  MG.boost*[1 -1];
  end
  if MG.ikriz == 2 || MG.nfast == 0
   if MG.inova == 4
    if MG.nfast == 0,  Amx(23,kAf)          =  Ipzf;
                       Amx(23,kAm)          =  Ipzm;
    else,              Amx(23,kAf)          =  Ipzf*(  2*MG.rampt);
                       Amx(23,kAm)          =  Ipzm*(2-2*MG.rampt);
    end
    if MG.iscramb ~=0
     switch MG.nfast
      case  -2    ,    Amx( 3,kAm)          =  Ipzmm;
      case {-3 -5},    Amx( 3,kAm)          =  Ipzm;
      end
    end
   else % MG.inova ~= 4
    switch MG.nfast
     case  -2    ,     Amx(23,kAm)          =  Ipzmm;
     case {-3 -5},     Amx(23,kAm)          =  Ipzm;
     case   0    ,     Amx(23,kAm)          =  Ipzm; % fort_equiv (seems new addition in .f)
    end
   end
   switch MG.nfast
    case {-3 -5},      Amx(24,95:114)       =  Ipzdm([2:6 8 11 17:19 21:23 29 32 34:38]);
   end
  end
  if MG.nfast == -3,   Amx(23,kAm([14 26])) =  MG.uc1*[-1 1];
   if MG.inova == 4 && MG.iscramb~=0
                       Amx( 3,kAm([14 26])) =  MG.uc1*[-1 1];
   end
  end
  switch MG.inova
   case {3 4},         Amx(20,kAa)          = -Ipza;
  end
  
  % Kappa
  Amx(22,kAf) =  Ipkf;
  Amx(22,kAm) =  Ipkm;
  if ismember(MG.ifour,[10:19 30:39]), Amx(22,kAa) = -Ipka; end
  %TODO: Add Kappa observer for mvloop=5,6 - iscramb>0 - MGA.F#L5971
  
  % E and F
  if MG.iscramb~=0
   if ismember(MG.mvloop,[1 2 11 12 21]) || MG.ierat == 4, Nca(E1,E1      ) = 0; Tca(E1,E1) = 0; end
   if ismember(MG.mvloop,[  2    12   ])                 , Nca(E2,E2      ) = 0; Tca(E2,E2) = 0; end
   if MG.inova == 4                                      , Nca(E3,E3      ) = 0; Tca(E3,E3) = 0; end
   if ismember(MG.mvloop,[1 2    12   ]) || MG.ierat == 4, Amx(E1,kAa(OH1)) = 1; end
   if ismember(MG.mvloop,[    11    21])                 , Amx(E1,kAvloop ) = 1; end
   %TODO: Add mvloop=5,6 - MGA.F#L5244
  end
  Amx(1:16,kAa(EF)) = Nca;
  
  % iscramb<0 case for external control tests
  if MG.iscramb < 0 || MG.ierat == -1, Amx(22,       :) = 0;
                                       Amx(22,kAa(OH1)) = 1;end
  % iscramb=0 case
  if MG.iscramb == 0, Amx(18,       :) = 0;
                      Amx(18,kAa(OH1)) = 1;end
  
 else % midplane == 8
  if MG.ikriz == 1 %choose b/w new vertical observer or F.Hofmann's observer
   % New vertical observer for doublet
   Amx(20,   1:38) =    2*MG.rampt *Ipzf_t(1:38);
   Amx(20,39:2*38) = 2*(1-MG.rampt)*Ipzm_t(1:38);
   Amx(20,  77:94) = -Ipza_t;
   Amx(22,   1:38) =    2*MG.rampt *Ipzf_b(1:38);
   Amx(22,39:2*38) = 2*(1-MG.rampt)*Ipzm_b(1:38);
   Amx(22,  77:94) = -Ipza_b;
  else
   % Vertical observer
   Amx(20,   1:38) =    MG.rampt *IpDz_xf_t(1:38);
   Amx(20,39:2*38) = (1-MG.rampt)*IpDz_xm_t(1:38);
   Amx(22,   1:38) =    MG.rampt *IpDz_xf_b(1:38);
   Amx(22,39:2*38) = (1-MG.rampt)*IpDz_xm_b(1:38);
   Amx(24,   1:38) =             -Ipzf(1:38);
   Amx(24,39:2*38) =              Ipzm(1:38);
  end
  % E and F coils
  Tca = eye(16);
  % z=+-0.4m actuators
  % Tca([9 12 13 16],[9 12 13 16]) = ...
  %     [ 0.540 -0.580  0.264 -0.064;
  %       0.226  0.908 -0.765  0.212;
  %       0.212 -0.765  0.908  0.226;
  %       0.064 -0.264  0.580 -0.540];
  % Amx(1:16,76+(1:16)) = Tca;
  % z=+-0.25m actuators
  Tca([10 12 13 15],[10 12 13 15]) = ...
      [ 0.540 -0.580  0.264 -0.064;
        0.226  0.908 -0.765  0.212;
        0.212 -0.765  0.908  0.226;
        0.064 -0.264  0.580 -0.540];
      
  % to include ierat=4 for doublet
  if MG.ierat == 4, Tca(E1,E1) = 0; end
  if MG.ierat == 4, Amx(E1,kAa(OH1)) = 1; end
  
  Amx(1:16,76+(1:16)) = Tca;
  
  % Radial observer
  Amx([18 23],kAf) = -Radf;
  Amx([18 23],kAm) = -Radm;
  
 end
 
 % Plasma current
 Amx(19,kAm)   = Ipm;
 
 % OH1-OH2
 Amx(17,kADOH) = 1+(MG.iscramb>= 0); % Legacy factor 2 for DOH observer?
 
 % Gas
 Amx(21,kAne)  = 1;
 
 % Conversion from relative to absolute flux input, to be removed when correcting fbte18
 k = 1 + (MG.mvloop == 21);
 Amx(k:end,kAf(1)) = -sum(Amx(k:end,kAf(2:end)),2);
 
 Amx = repmat(Amx,[1,1,nsw]);
 if midplan == 8
  Amx(18,:,4) = Amx(24,:,4);
  Amx(24,:,:) = 0;
 end
 
 %% G matrix (GMATRIX)
 mG  = 22; nG  = 24; Gmx = zeros(3,mG,nG);
 GI = 1; GP = 2; GD = 3;
 kMFPS = 22; kMGAS=21;               % indices in M matrix inputs (G outputs)
 kGIp = 19; kGDOH = 17; kGGAS  = 21; % indices in G matrix inputs (A outputs)
 tauD  = 0.006/(2*5.25*pi);
 kGOH = NaN;
 if MG.ierat   ==  4, kGOH =  1;end
 if MG.ierat   == -1, kGOH = 22;end
 if MG.iscramb  <  0, kGOH = 22;end
 if MG.iscramb ==  0, kGOH = 18;end
 
 if midplan ~= 8
  % Gas
  Gmx([GP GI GD],21,21) = [MG.gain MG.ggain MG.hgain]/60;
  
  % OH coils (Ip and DOH)
  Gmx([GP GI],OH,kGIp ) = -repmat(MG.aipgain(1:2)',1,2);
  Gmx([GP GI],OH,kGDOH) = [MG.ohsame ; MG.ohsami]*[1 -1]*0.25;
  
  % EF coils (coil combinations, radial, vertical, kappa)
  Gmx(GP,EF,EF) = Tca'; Gmx(GP,E1:E8,E1:E8) = 0.25*Gmx(GP,E1:E8,E1:E8);
  if nr>0
   Gmx(GP,EF,18) = MG.gainr            * Tra(1,:);
   if ismember(MG.inova,[2 3 4]), k = 23;
   else,                          k = 20; end
   Gmx(GP,EF,20) = MG.gainz*0.05       * Tra(2,:);
   Gmx(GD,EF,k ) = MG.gainvz*0.05*tauD * Tra(2,:);
   if nr >= 3
    Gmx(GP,EF,22) = MG.gainext(1)      * Tra(3,:);
   end
  end
  if MG.iscramb == 0
   Gmx(GP,E1:E8,E1:E8) = 0.35*eye(8);
   Gmx(GP,F1:F8,F1:F8) = 1.40*eye(8);
   Gmx(GI,E1:E8,E1:E8) = 0.10*eye(8);
   Gmx(GI,F1:F8,F1:F8) = 0.40*eye(8);
  end
  
  if fort_equiv
   if ismember(MG.mvloop,[ 1  2 11 12 21]) || MG.ierat == 4, Gmx(GP,1,1) = 0;end
   if ismember(MG.mvloop,[    2    12   ]),                  Gmx(GP,2,2) = 0;end
  end
  
  % FPS
  if MG.nfast ~= 0
   if     MG.ikriz ~=  2, kGFPS = 20;
   elseif MG.nfast == -5, kGFPS = 24;
   elseif MG.inova ==  4, kGFPS =  3;
   else                 , kGFPS = 23;
   end
   Gmx(GD,kMFPS,kGFPS) = -MG.gainvz*0.1*tauD;
  end
  
  % Switching
  Gmx = repmat(Gmx,[1,1,1,nsw]);
  
  % OH
  if MG.iscramb ~= 0
   if ismember(MG.ierat,[-1,4]) || MG.iscramb < 0
    if ~MG.aipgain(3),MG.aipgain(3)=0.39;end
    Gmx( GP    ,OH,kGIp, 1   ) = 0;
    Gmx([GP GI],OH,kGOH, 1   ) = repmat([1;5]*MG.aipgain(3),1,2);
   end
   if MG.iscramb<0
    Gmx(GI,OH,kGOH,1) = 0; % hard set integral gains to 0
   end
   if MG.ierat == 4
    Gmx( GP    ,OH,kGIp,[2,8]) = repmat(-min([abs(MG.aipgain(1)),0.01]),[1 2 1 2]);
   elseif MG.ierat ~= -1
    Gmx( GP    ,OH,kGIp,[1,2]) = repmat(-min([abs(MG.aipgain(1)),0.01]),[1 2 1 2]);
   end
  else %MG.iscramb == 0
   Gmx([GP GI],OH,kGIp ,[1,7:10]) = 0;
   Gmx( GP    ,OH,kGOH ,[1,7:10]) = 0.39;
  end
  Gmx( GI    ,OH,kGIp,[1:2,8]) = 0; % No integral gain on current in first two matrices
  
  % EF, .f keeps GD on 20
  if nr>0, Gmx(:      ,EF,kc(1:2), 2:6 ) = 0;end
  % NOTE: I can only assume the following is the result of a .f bug
  if fort_equiv && MG.iscramb == 5, Gmx(GP,12,kc(1),2:6) = 1; end 
  Gmx(:      ,EF,     18,[1 7]) = 0;
  Gmx([GP GI],EF,     20,[1 7]) = 0;
  Gmx( GD    ,EF,     20, 1   ) = 0;
  Gmx(:      ,EF,     23, 1   ) = 0;
  if nr >= 3, Gmx(:,EF,kc(3)  , 5:6   ) = 0;
              Gmx(:,EF,22     ,[1:4 7]) = 0; end
  % No gains for vertical, radial and kappa control for matrix 8 at t=0s
  Gmx(:      ,EF,     18,8) = 0;
  Gmx([GP GI],EF,     20,8) = 0;
  Gmx( GD    ,EF,     20,8) = 0;
  Gmx(:      ,EF,     22,8) = 0;
  Gmx(:      ,EF,     23,8) = 0;
  % Open-loop for all coils for matrix 7 in the case iscramb<0
  if MG.iscramb < 0
   Gmx(:    ,EF,      :,7) = 0;
   Gmx(:    ,OH,      :,7) = 0;
  end
  % FPS
  Gmx(:,kMFPS,[3 20 23],[1 2]) = 0;
  % FPS no gains at t=0s matrix 8
  Gmx(:,kMFPS,[3 20 23],8) = 0;
  
  % Gas
  Gmx(: ,kMGAS,kGGAS,1) = 0  ;
  Gmx(GP,kMGAS,kGGAS,2) = 0.1;
  % Gas no gains at t=0s matrix 8
  Gmx(: ,kMGAS,kGGAS,8) = 0  ;
  % No gas feedback for iscramb = -2
  if MG.iscramb < -1
   Gmx(: ,kMGAS,kGGAS,:) = 0;
  end
  
  % ECCD
  switch MG.mvloop
   case   1    , Gmx(GP,OH, 1   ,4:5) =   0.39;
   case {11 21}, Gmx(GP,OH, 1   ,4:5) = 100   ;
   case { 2 12}, Gmx(GP,OH,[1 2],4:5) =   0.39;
  end
  Gmx(:,OH,kGIp,4) = 0;
  switch MG.mvloop
   case {1 2 11 12 21}, Gmx(:,OH,kGIp,5) = 0;
  end
  %TODO: mvloop = 5,6 and derived (incompatible with iscramb=0) - MGA.F#L4974
   
  % Special matrices
  % Matrix 7 (???)
  Gmx(:,kMGAS,kGGAS,7) = Gmx(:,kMGAS,kGGAS,3);
  Gmx(GD,:,[3,20,23,24],7) = Gmx(GD,:,[3,20,23,24],3)*0.45;
  % Matrix 9 (3 + Gas FF)
  Gmx(:,:,:,9) = Gmx(:,:,:,3);
  Gmx(:,kMGAS,kGGAS,9) = 0;
  % Matrix 2 (1 + Gas FB P only)
  if MG.inova == 5
   Gmx(:,:,:,2) = Gmx(:,:,:,1);
   Gmx(GP,kMGAS,kGGAS,2) = Gmx(GP,kMGAS,kGGAS,3);
  end
  % Matrix 10 (3 + No FPS)
  Gmx(:,:,:,10) = Gmx(:,:,:,3);
  Gmx(:,kMFPS,:,10) = 0;
  
 else % midplan == 8
  % init Gmx
  Gmx = zeros(3,mG,nG,nsw);   
  % Gas
  Gmx([GP GI GD],kMGAS,kGGAS) = [MG.gain MG.ggain MG.hgain]/60;
  
  % OH
  Gmx([GP GI],OH,17) = [MG.ohsame ; MG.ohsami]*[1 -1]*0.25;
  
  % Plasma current
  Gmx(GP,OH,kGIp) = -repmat(MG.aipgain(1)',1,2);
  
  % EF
  Gmx(GP,EF,EF) = Tca';
  
  
  % Gas
  Gmx(: ,kMGAS,kGGAS,[1 8]) = 0  ;
  Gmx(GP,kMGAS,kGGAS,2) = MG.gain/60;
  
  % delta OH
  Gmx([GP GI],OH,kGDOH,2) = [MG.ohsame ; MG.ohsami]*[1 -1]*0.25;
  % OH
  Gmx([GP GI],OH,kGIp,1:2) = repmat(-[min(MG.aipgain(1),0.01);0],[1 2 1 2]);
  % 3 and 4 matrix for doublet
  Gmx([GP GI],OH,kGIp,3:4) = repmat(-[min(MG.aipgain(1),0.01);0],[1 2 1 2]);
  % OH Ip control at t=0s for matrix 8
  Gmx([GP GI],OH,kGIp,8) = repmat(-[max(MG.aipgain(1),0.01);0],[1 2 1 1]);
  
  % EF
  Gmx(GP,E1:E8,E1:E8,1:4)=repmat(0.25*squeeze(Gmx(GP,E1:E8,E1:E8)),[1,1,4]);
  
  % F coils not used for position control F1 F3 F6 F8 z=+-0.25m
  Gmx(GP,[9 11 14 16],[9 11 14 16],1:4) = repmat(eye(4),[1,1,4]);
  
  % IOH control
  if MG.ierat == 4, Gmx(:      ,OH,19,1) = 0;
                    Gmx([GP GI],OH, 1,1) = repmat([1;40]*1.3,1,2); end
     
  % % F coils used for position control F1 F4 F5 F8 coils
  % Gmx(GP, 9,[9 12 13 16],1) =  1*[ 1.167  1.177  0.350  0.472];
  % Gmx(GP,12,[9 12 13 16],1) =  1*[-1.069  1.448  0.856  1.056];
  % Gmx(GP,13,[9 12 13 16],1) =  1*[-1.056  0.856  1.448  1.069];
  % Gmx(GP,16,[9 12 13 16],1) =  1*[-0.472  0.350  1.177 -1.167];

  % z=+-0.25m F2 F4 F5 F7 to be used for the position control
  Gmx(GP,10,[10 12 13 15],1) = 1*[ 1.167  1.177  0.350  0.472];
  Gmx(GP,12,[10 12 13 15],1) = 1*[-1.069  1.448  0.856  1.056];
  Gmx(GP,13,[10 12 13 15],1) = 1*[-1.056  0.856  1.448  1.069];
  Gmx(GP,15,[10 12 13 15],1) = 1*[-0.472  0.350  1.177 -1.167];
  
  % z=+-0.4m
  %     c = [ 0.770  0.150  0.091  0.141;
  %          -0.827  0.603 -0.377 -0.508;
  %           0.377 -0.508  0.827  0.603;
  %          -0.091  0.141 -0.770  0.150];
  % Using F2 F4 F5 F7  as actuators for z=+-0.25m
  c = [0.6938 -0.1349  0.0769  0.3117;
      -0.6618  0.6845 -0.2734 -0.6451;
       0.2734 -0.6451  0.6618  0.6845;
      -0.0769  0.3317 -0.6938 -0.1349];
   
   
  switch iplatop
   case 1
    c2 = c(:,[1 2 3 4]).*[1 1        1      MG.gainr];
    c3 = c(:,[1 2 4 3]).*[1 1 MG.gainr -0.1*MG.gainz];
    cd = c(:,       3 )              *(-0.1*MG.gainvz*tauD);
    Gmx(GP,[9 12 13 16],[9 12 16 18   ],2) = c2;
    Gmx(GP,[9 12 13 16],[9 12    18 20],3) = c3;
    Gmx(GD,[9 12 13 16],            20 ,3) = cd;
   case -1
    c2 = c(:,[1 4 3 2]).*[1 1             1 MG.gainr];
    c3 = c(:,[4 3 1 2]).*[1 1 -0.1*MG.gainz MG.gainr];
    cd = c(:,       1 )     *(-0.1*MG.gainvz*tauD);
    Gmx(GP,[9 12 13 16],[9 13 16    23],2) = c2;
    Gmx(GP,[9 12 13 16],[  13 16 22 23],3) = c3;
    Gmx(GD,[9 12 13 16],         22    ,3) = cd;
   case 0
    % c2 = c(:,[1 3 4 2]).*[       1             1      MG.gainr MG.gainr];
    c3 = c(:,[4 3 1 2]).*[MG.gainr -0.1*MG.gainz  0.1*MG.gainz MG.gainr];
    cd = c(:,[3 1])    .*[-1 1]*0.1*MG.gainvz*tauD;
    % c3 = c(:,[4 3 1 2]).*[MG.gainr  0.1*MG.gainz -0.1*MG.gainz MG.gainr];
    % cd = c(:,[3 1])    .*[1 -1]*0.1*MG.gainvz*tauD;
    % Gmx(GP,[9 12 13 16],[9 16 18          23],2) = c2;
    % Gmx(GP,[9 12 13 16],[        18 20 22 23],3) = c3;
    % Gmx(GD,[9 12 13 16],[           20 22   ],3) = cd;
    % Gmx(GP,[9 12 13 16],[        18 20 22 23],2) = c3;
    % Gmx(GD,[9 12 13 16],[           20 22   ],2) = cd;
    % z=+0.25m
    % Gmx(GP,[10 12 13 15],[10 15 18          23],2) = c2;
    Gmx(GP,[10 12 13 15],[        18 20 22 23],2) = c3;
    Gmx(GD,[10 12 13 15],[           20 22   ],2) = cd;
   end
  Gmx(:,:,:,3) = Gmx(:,:,:,2);
  Gmx(:,:,:,4) = Gmx(:,:,:,2);
 end
  
 %% WRITE MDS
 % To comply with .f,, to be removed
 if nr == 3                           , k = [1 4 5 2 6 3];
 elseif ismember(MG.iscramb,[8 10 11]), k = [3 2 1 4];
 elseif length(kc) == 4               , k = [3 1 2 4];
 elseif isempty(kc)                   , k = [];
 else                                 , k = [3 4 1 5 6 2];
 end
 Amx(  kc,:,:) = Amx(  kc(k),:,:);
 Gmx(:,:,kc,:) = Gmx(:,:,kc(k),:);
 Ref(:,kc)     = Ref(:,kc(k));
 if midplan ~= 8, Ref(:,[19 20]) = 0; end
 if ~fort_equiv, kref(RF.t <= 0) = 0; end
 
 MGA.AAASW   = Amx                         ;
 MGA.AIPOFT  = SC.Ip(:)*MG.curfac+MG.aipdel;
 MGA.AMLSW   = Mmx                         ;
 MGA.BZEROFT = SC.Bt(:)                    ;
 MGA.CAYREF  = kref'                       ;
 MGA.EFCUR   = RF.Ia'                      ;
 MGA.EFVOLT  = RF.Ua'                      ;
 MGA.EFWAVE  = Ref                         ;
 MGA.G1SW    = squeeze(Gmx(GI,:,:,:))      ;
 MGA.G2SW    = squeeze(Gmx(GP,:,:,:))      ;
 MGA.G3SW    = squeeze(Gmx(GD,:,:,:))      ;
 MGA.MIDPLAN = midplan                     ;
 MGA.OHCUR   = 0*SC.t(:)                   ;
 MGA.TOFT    = SC.t(:)                     ;
 MGA.NTOFT   = SC.nt                       ; % to be removed for fbte18.m
 MGA.TOH     = RF.t(:)                     ;
 MGA.ZEREF   = zref1'                      ;
 
 for node = fieldnames(MGA)'
  data = MGA.(node{1});
  assert(all(isfinite(data(:))),'TCVCS:MGA:NoConvergence','NaNs found in MGA result field "%s"',node{1});
 end
 
 if nargout == 0
  for node = fieldnames(MGA)'
   mdsw(['\MGA.' node{1}],MGA.(node{1}))
  end
 end
end

function [Iyf,Iym,Iya,Iyv] = Iy(Wf,Wm,Wp,Mfa,Mfv,Mfy,Bma,Bmv,Bmy,Ipm,Ipy,D)
 Ayf = Mfy'*diag(Wf);
 Bym = Bmy'*diag(Wm);
 Aym = Bym     + Ipy'*Wp(1)*Ipm;
 Byy = Bym*Bmy + Ipy'*Wp(2)*Ipy;
 Ayy = Ayf*Mfy + Byy;
 Aya = Ayf*Mfa + Aym*Bma;
 Ayv = Ayf*Mfv + Aym*Bmv;
 Iyf = D*(Ayy\Ayf);
 Iym = D*(Ayy\Aym);
 Iya = D*(Ayy\Aya);
 Iyv = D*(Ayy\Ayv);
end

function mdso(shot)
 [~,s] = mdsopen('PCS',shot);
 assert(isodd(s),'Cannot open PCS tree for shot %d',shot);
end

function mdsw(node,x)
 assert(isodd(mdsput(node,mdscvt(double(x),'f'))),'TCVCS:MGA:mdsputFailed','Can not write MDS node %s',node)
end

function x = mdsr(e,varargin)
 [x,s] = mdsdata(e,varargin{:});
 assert(isodd(s),'Cannot evaluate MDS expression %s',e)
 if iscellstr(x), x = deblank(x); end
end
