function fbte18(shot,index,pcsconfig,ioh1values)
%% Mlint exceptions
%#ok<*STRCLFH> % avoid contains()

%% Constants
MU0   = 4e-7*pi;
R0    = 0.88;
NWTOR = 16*6; % total windings in TF coils

%% Housekeeping
narginchk(0,4)
if nargin < 1 || isempty(shot)     , shot      = -1      ; end
if nargin < 2 || isempty(index)    , index     = 108     ; end
if nargin < 3 || isempty(pcsconfig), pcsconfig = 'hybrid'; end

pcsconfig = lower(pcsconfig);
mdsopen('PCS',shot)
mdsw('\TOP:ID',mfilename);

%% read MGA
MGA = mgaget(shot,false); % Convert to doubles
MG = mgp(shot);

Ax      = MGA.AAASW;
aipoft  = MGA.AIPOFT;
Mx      = MGA.AMLSW;
bzeroft = MGA.BZEROFT;
cayref  = MGA.CAYREF;
efcur   = MGA.EFCUR;   efcur(:,18) = efcur(:,17);
efvolt  = MGA.EFVOLT;
efwave  = MGA.EFWAVE;
G1x     = MGA.G1SW;
G2x     = MGA.G2SW;
G3x     = MGA.G3SW;
midplan = MGA.MIDPLAN;
ohcur   = MGA.OHCUR;   dioh = ohcur; % MGA convention
toft    = MGA.TOFT;
toh     = MGA.TOH;
zeref   = MGA.ZEREF;

ierat   = MG.ierat;
ikriz   = MG.ikriz;
inova   = MG.inova;
iohfb   = MG.iohfb;
iscramb = MG.iscramb; iscramb=rem(iscramb,100);
mvloop  = MG.mvloop;
nfast   = MG.nfast;
relovo  = MG.relovo;
istop   = MG.istop;

%% Convert to struct with additional parameters from debug_p3
MG = struct_with_debug(MG,3);

%% Local variables

load_user_refs     = iscramb  <  0;
rm_den_ip_ipz_refs = iscramb == -2;

%% A matrix
% inputs and outputs
Ain = lstmk('CTLINT:VLOOP_001',iii('CTLINT:RVLOOP_'  ,2:38),iii('CTLINT:BPOL_AVG_',1:38),...
 iii('CTLAMP:IPOL_E_'  ,1:8),iii('CTLAMP:IPOL_F_'  ,1:8),iii('CTLAMP:IPOL_OH_' ,[1 2]),...
 iii('CTLAMP:BPOL_AVG_',[2 3 4 5 6 8 11 17 18 19 21 22 23 29 32 34 35 36 37 38]),...
 'CTLAMP:IPOL_DOH_001','FIR_FRINGE','CTLAMP:IPOL_G_001','1MM_FRINGE','CTLAMP:VLOOP_001','CTLAMP:IPHI_001');
Aoutk = [1:11 20 19 14:18 13 12 21 22:24];
Aout = lstmk(iii('CP',Aoutk));

% gates, no more in use
if midplan ~= 8
 if ismember(iscramb,[12 15 16 19 25 26 36]) || mvloop == 6
  switch nfast
   case  0,      tmp = lstmk(repmat({'_P'},1,16)                 ,'_DELTA_OH','_DELTA_R','IPA','IP*ZA','DENSA','A','!','!');
   case -5,      tmp = lstmk(repmat({'_P'},1,16)                 ,'_DELTA_OH','_DELTA_R','IPA','IP*ZA','DENSA','A','A','A');
   otherwise,    tmp = lstmk(repmat({'_P'},1,16)                 ,'_DELTA_OH','_DELTA_R','IPA','IP*ZA','DENSA','A','A','!');
  end
 else
  switch nfast
   case  0,      tmp = lstmk(repmat({'_P'},1,16)                 ,'_DELTA_OH','_DELTA_R','IPA','IP*ZA','DENSA','!','!','!');
   case -5,      tmp = lstmk(repmat({'_P'},1,16)                 ,'_DELTA_OH','_DELTA_R','IPA','IP*ZA','DENSA','!','A','A');
   otherwise,    tmp = lstmk(repmat({'_P'},1,16)                 ,'_DELTA_OH','_DELTA_R','IPA','IP*ZA','DENSA','!','A','!');
  end
 end
else,            tmp = lstmk(repmat({'_P'},1,16)                 ,'_DELTA_OH','_DELTA_R','IPA','IP*ZA','DENSA','A','A','!');
end
if mvloop == 18, tmp = lstmk(repmat({'!'},1,8),repmat({'_P'},1,8),'_DELTA_OH','_DELTA_R','IPA','IP*ZA','DENSA','A','A','!');
end
Aout = strcat(Aout,tmp);

% output scales in MKSA/V
%                                                     A     A   Ip          Amfr     A     Am^2
if ~ismember(mvloop,[1 11 21 2 12]) && rem(mvloop,10) ~= 5 && ierat ~= 4
                     Aoutf     = [     repmat(2200,1,16) 7480 2200 1.1e5 2.7e4 4 1.4e4 1.8e4 2.4e6];
 if mvloop == 6,     Aoutf( 1)  = 1.4e4;
 end
 if ierat == -1 || iscramb < 0
                     Aoutf(22)  = 7480;
 end
elseif mvloop == 1 || rem(mvloop,10) == 5 || ierat == 4
                     Aoutf     = [7480 repmat(2200,1,15) 7480 2200 1.1e5 2.7e4 4 1.0e4 1.8e4 1.2e6];
 if rem(mvloop,10) == 5
  if mvloop < 100,   Aoutf(19) =  1.1e5/1.5;
  else,              Aoutf( 4) =  1.1e5/3;
  end
 end
elseif mvloop == 11, Aoutf     = [  10 repmat(2200,1,15) 7480 2200 1.1e5 2.7e4 4 1.0e4 1.7e4 1.2e6];
elseif mvloop == 21, Aoutf     = [   1 repmat(2200,1,15) 7480 2200 1.1e5 2.7e4 4 1.0e4 1.7e4 1.2e6];
elseif mvloop == 2 || mvloop == 12
                     Aoutf     = [7480 repmat(2200,1,15) 7480 2200 1.1e5 2.7e4 4 1.0e4 1.7e4 1.2e6];
                     Aoutf( 2) =  max(abs(ioh1values))/.85; % based on actual max current required, to maximize resolution (0.85 V -> 8.5 V at G2 input)
end
if inova == 4,       Aoutf( 3) =  1.7e4; end
if     midplan == 8, Aoutf     = [     repmat(2200,1,16) 7480 1200 1.0e5 8.0e3 4 8.0e3 1.2e3 2.2e5];
elseif midplan == 9, Aoutf     = [     repmat(2200,1,16) 7480 2200 1.0e5 2.3e4 4 2.2e3 2.2e3 2.2e3];
end
if iscramb == 9,     Aoutf(18) =  7480; end
% Adjust disruption detector sensitivity to Ip
if istop(1)<0,       Aoutf(19) = (abs(istop(1))*1000-2e4)/.67;end
% Adjust Ipz scale
if rm_den_ip_ipz_refs
                     Aoutf(20) = Aoutf(20)*4;
end

% 3D to 2D, real to virtual flux loop transformation
Ax = Ax(:,:,1);
Ax(:,1:38) = Ax(:,1:38)*mdsu('MAGNETICS','\T_F_R'); 
Ax(:,1) = 0; % keep relative fluxes only, discard absolute flux term

% Remove Ipz observer
if rm_den_ip_ipz_refs
  Ax(20,:) = 0;
end

%% read voltages and references from file if iscramb=-1 or -2
if load_user_refs
  refdir = '~/mgams_ref/'; fname = fullfile(refdir,'tcv_coil_traces.mat');
  if isfield(MG,'fname'), fname = MG.fname; end
  LL=load(fname); fprintf('loaded %s\n',fname)
  efcur = LL.efcur;
  efvolt = LL.efvolt;
  toh = LL.toh; % times for efcur, efvolt
  tC = LL.tC; % times for switching between FF and FB
  % Load custom Ip trace if it exists
  if isfield(LL,'aipoft')
    aipcustom = LL.aipoft;
  else
    aipcustom = 0*toh;
  end
  assert(numel(cell2mat(tC))>1,'must have more than 1 matrix switching time in tC')
  dioh = efcur(:,17) - efcur(:,18);
  
  assert(              iscolumn(toh),'loaded toh must be a column vector');
  assert(size(efcur,1) == numel(toh),'loaded efcur must have as many rows as numel(toh)')
  assert(size(efvolt,1)== numel(toh),'loaded efvolt must have as many rows as numel(toh)')
  assert(size(aipcustom,1)== numel(toh),'loaded aipcustom must have as many rows as numel(toh)')
end

% matrix switching
teq = mdsdata('fbte_time("FBTE")');
teq = max(ceil(teq(end)/.01)*.01,0.42);
if midplan ~= 8
 if ismember(mvloop,[1 2 11 12 21]) || rem(mvloop,10) == 5,   tI  = 0.40;
  if ismember(iscramb,[12 15 16 19 25 26 36]),                tsw = [-0.30 0.01 0.07 tI teq];     ksw = [0 1 2 4 0];
  else,                                                       tsw = [-0.30 0.01 0.07 tI teq];     ksw = [0 1 2 3 0];
  end
 else
  if ismember(iscramb,[12 15 16 19 25 26 36]) || mvloop == 6, tI = 0.30;
                                                              tsw = [-0.30 0.01 0.07 tI teq]; ksw = [0 1 2 4 0];
  else,                                                       tsw = [-0.30 0.01 0.07 teq];    ksw = [0 1 2 0];
  end
 end
else,                                                         tsw = [-0.30 0.01 teq];         ksw = [0 2 0];
end
if iscramb == 0,                                              tsw = [-0.30 0.01 0.07 3.0 3.10];    ksw = [0 0 0 0 0];
end
% Temporary for IOH control, for runaways and more generally
if mvloop(1)==1 && (istop(1)<0 || istop(1)==666)
  assert(istop(3)<teq,'IOH control window extends past planned shot termination, set istop(3)<%.3f',teq);
  offpre = find(ksw==3|(tsw>=istop(2)&tsw<=istop(3))); tsw(offpre)=[]; ksw(offpre)=[];
  [tsw,iremap] = sort([tsw,istop(2),istop(3)]);
  ksw = [ksw,3,2];ksw = ksw(iremap);
end
if load_user_refs
  % switch matrices based on timings specified in tC
  [tsw,isw] = sort(cell2mat(tC));
  ksw = [zeros(1,length(tC{1})), 6*ones(1,length(tC{2}))];
  ksw = ksw(isw);
end
wave_put('\PHYS_MAT_ADDRESSES:A',zeros(size(ksw)),tsw)
wave_put('\PHYS_MAT_ADDRESSES:G',           ksw  ,tsw)
wave_put('\PHYS_MAT_ADDRESSES:M',zeros(size(ksw)),tsw)

% amplifier gains
Again = [1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 1 2 2 2 1 4 8 8];
if         ierat      ==   4, Again(1)  = 4; end
if	       ierat      ==  -1 || iscramb < 0
                              Again(22) = 4;
end
if         mvloop     ==  11, Again( 1) = 8;
elseif rem(mvloop,10) ==   5, Again(19) = 8;
 if        mvloop     >  100, Again( 4) = 8; end
elseif     mvloop     ==   6, Again( 1) = 8;
end
if         inova      ==   4, Again( 3) = 2; end
% Adjust Ipz amplifier gain
if rm_den_ip_ipz_refs,        Again(20) = Again(20)*2; end

% ### Put the A Matrix
pcsmatput('A',Ax,Ain,Aout,Aoutf,Again)

%% G matrix
Gin = Aout;
Gout = lstmk(iii('GO',1:22));
if nfast == -5, Goutf    = [repmat(44,1,16) 149.6 149.6 2000 240 0.15 88];
else,           Goutf    = [repmat(44,1,16) 149.6 149.6 2000 240 0.15 30];
end
if rem(mvloop, 10) ==  5, Goutf(1) = 1.e4; end % RHVPS
if rem(mvloop,100) == 25, Goutf(2) = 1.e4; end % RHVPS
if mvloop == 6, error('UNSTRANLATED'), end

% ### PID gain, not used
mdsw(strcat('\LOAD_TAU:',{'D' 'I' 'P'}),0)
% ### Put the G Matrix
pcsmatput('G1',G1x,Gin,Gout,Goutf)
pcsmatput('G2',G2x,Gin,Gout,Goutf)
pcsmatput('G3',G3x,Gin,Gout,Goutf)

%% M matrix
Min  = lstmk(Gout,strcat('CTLAMP:IPOL_',lstmk(iii('E_',1:8),iii('F_',1:8),'OH_001','OH_002')));
Mout = lstmk(                                 iii('E_',1:8),iii('F_',1:8),'OH_001','OH_002','GAS','FAST');
switch rem(mvloop,100)
 case { 5 25}, Mout{[1 2]} = {'RHVPSA' 'RHVPSB'};
 case  15    , Mout{ 1   } =           'RHVPSB';
end
if mvloop == 6, error('UNSTRANLATED'), end
Moutf = zeros(length(Mout));
Mx = Mx(:,:,1);
% ### Put the M matrix
pcsmatput('M',Mx,Min,Mout,Moutf)

%% Loadable matrices
% ### Prepare loadable matrices
pcsmatdo(index,pcsconfig)
pcsmatoverflow(shot)

%% 3-way summator
if shot == -1
 assert(logical(mdsr('JP_SET_SUM_SYS($)',pcsconfig)),'THREE-WAY SUMMATOR SET-UP FAILED')
end

%% Power supply settings
% read currents from file if iscramb=0 and inova=5
if iscramb == 0 && inova == 5, error('UNTRANSLATED'), end

if     iscramb <= 0
	alimmode = 'BBBBBBBBBBBBBBBBBBBA';
elseif ismember(mvloop,[2 12])    || rem(mvloop,100) == 25
	alimmode = 'CCBBBBBBBBBBBBBBBBBA';
elseif ismember(mvloop,[1 11 21]) || rem(mvloop, 10) ==  5 || ierat == 4
	alimmode = 'CBBBBBBBBBBBBBBBBBBA';
elseif mvloop == 18
	alimmode = 'CCCCCCCCBBBBBBBBBBBA';
else
	alimmode = 'BBBBBBBBBBBBBBBBBBBA';
end
if inova == 4 && iscramb > 0 && ierat ~=-1           , alimmode( 3) = 'C'; end
if rem(mvloop,10) == 5 && mvloop > 100 && iscramb > 0, alimmode( 4) = 'C'; end
if nfast ~= 0                                        , alimmode(20) = 'B'; end
if mvloop == 6, error('UNTRANSLATED'), end
alimmode = lstmk(alimmode');
alimsel = repmat({'ON'},1,20);
alimlim = repmat({'50'},1,18);
lim50 = [repmat(0.75*3850,1,16) repmat(27000*0.5,1,2)];
lim80 = [repmat(0.75*5544,1,16) repmat(27000*0.8,1,2)];
tmp = max(abs(efcur));
alimlim(tmp > lim50) = { '80'};
alimlim(tmp > lim80) = {'100'};
alimlim([19 20]) = {'100' '75'};
alimsign( 1:16 ) = 1;
alimsign(any(efcur < -100)) = -1;
[~,im] = max(abs(aipoft));
alimsign([17 18]) = sign(aipoft(im));

% ### Put the alim settings
pcssetalim(alimmode,alimlim,alimsel,alimsign)

%% Wavegen
switch pcsconfig
 case 'hybrid'    , tmp = ones(1,128);
 case 'rtdtacq'   , tmp = [0 0 0 ones(1,13) zeros(1,16) ones(1,11) zeros(1,2) ones(1,3) zeros(1,16) ones(1,64)];
 case 'rtdtacqfps', tmp = [0 0 0 0 ones(1,12) zeros(1,16) ones(1,11) zeros(1,2) ones(1,3) zeros(1,16) ones(1,64)];
 case 'rtdtacqden', tmp = [1 1 0 ones(1,125)];
 otherwise, error('Invalid recipe selected %s',pcsconfig);
end
mdsw('\LOAD_WAVE_MASK',int8(tmp))

%% Feed-forward
if iscramb>0
  efvolt(:,17) = efvolt(:,17) + 20; % obscure correction
end

% Voltage feed forward
switch pcsconfig
 case 'hybrid'
  for i = 1:8, wave_put(['\DRAW_FEEDFOR_E_U:ALIM_' iii(i)],efvolt(:,i  ),toh,8), end
  for i = 1:8, wave_put(['\DRAW_FEEDFOR_F_U:ALIM_' iii(i)],efvolt(:,i+8),toh,8), end
               wave_put( '\DRAW_FEEDFOR_OH_U:ALIM_001'    ,efvolt(:,17 ),toh,8)
               wave_put( '\DRAW_FEEDFOR_OH_U:ALIM_002'    ,efvolt(:,18 ),toh,8)
 otherwise
  for i = 1:8, wave_put_scd(['\DRAW_FEEDFOR_E_U:ALIM_' iii(i)],efvolt(:,i  ),toh), end
  for i = 1:8, wave_put_scd(['\DRAW_FEEDFOR_F_U:ALIM_' iii(i)],efvolt(:,i+8),toh), end
               wave_put_scd( '\DRAW_FEEDFOR_OH_U:ALIM_001'    ,efvolt(:,17 ),toh)
               wave_put_scd( '\DRAW_FEEDFOR_OH_U:ALIM_002'    ,efvolt(:,18 ),toh)
end

% Compute sign bits
[tbit,vbit] = deal(cell(1,size(efcur,2)));
for k = 1:size(efcur,2)
 [tbit{k},vbit{k}] = pcssignbit(toh,efcur(:,k));
 assert(all(diff(tbit{k}) > 0),'pcssignbit failed to provide monotonic timebase for sign bits of alim #%d',k);
end

% OH current (always in voltage - or B - mode, containts sign bits)
t = min(toft) - [0.036 0.030 0.024];
nodes = mdsr('GETNCI("\\DRAW_FEEDFOR_OH_I:ALIM*","PATH")');
for ii = 1:numel(nodes)
 if iscramb >= 0
  % Use default sign bit sequence
  wave_put(nodes{ii},0.8*sign(iohfb)*[0 1000 0 0 -1000 0 0],[t t+0.500 toh(end)])
 else % iscramb<0
  % Use computed sign bit sequence based on custom IOH reference
  k = 16+ii;
  wave_put(nodes{ii},vbit{k},tbit{k});
 end
end

% E and F current (contains sign bits if in voltage - or B - mode)
for k = 1:16
  if k <= 8, n = ['\DRAW_FEEDFOR_E_I:ALIM_00' int2str(k  )];
  else,      n = ['\DRAW_FEEDFOR_F_I:ALIM_00' int2str(k-8)]; end
  if strcmp(alimmode{k},'B'), wave_put(n,vbit{k},tbit{k},8)
  else,                       assert(tbit{k}(3) < toh(1),'First sign bit overlaps with current trace');
                              % NOTE: This should never happen if pcssignbit works correctly.
                              wave_put(n,[vbit{k}(1:3);efcur(:,k)],[tbit{k}(1:3);toh])
  end
end

% TOR current
% sign
if mdsr('\pcs::mgams.data:if36fb') > 0
  mdsw('\pcs::tor_sign','1','x')
else
  mdsw('\pcs::tor_sign','-1','x')
end

% current trace
k = find(toft(3:end-2) > -0.103) + 2;
bphi = [0;bzeroft(2)/4;bzeroft([1;2;k;end-1;end])]*1.02; % correction for power supply error
tphi = [-0.353-[bphi(4)+2*bphi(2);bphi(4)+bphi(2);bphi(4);0]/1.2;toft([k;end-1]);toft(end-1)+bphi(end-1)/1.2];
Iphi = 2*pi*R0/MU0/NWTOR*bphi;
wave_put('\DRAW_FEEDFOR_TOR_I:ALIM_001',Iphi,tphi)

% FPS
t = toh([1 end]);
wave_put('\DRAW_FPS:UDCREF_001',ifc(nfast == 0 || iscramb <= 0,350,270)*[1 1],t-[0.1;0]  )
wave_put('\DRAW_FPS:UREF_001'  ,                                        [0 0],t        ,8)
wave_put('\DRAW_FPS:IREF_001'  ,                                        [0 0],t        ,8)

% GAS
wave_put('\DRAW_FEEDFOR_GAS:ALIM_001',mdsr('\DRAW_FEEDFOR_GAS:ALIM_001'),mdsr('DIM_OF(\DRAW_FEEDFOR_GAS:ALIM_001)'),8)

%% References
% specify whether the references are 'Ip' or '10V'
pcssetref(Aout,repmat({'10V'},1,24))
t = [min(toft(1),toh(1)) max(toft(end),toh(end))];
for k = [1:20 22:24], wave_put(['\DRAW_REFS:REF_' iii(k)],[0 0],t); end

% E and F
for k = find(strcmp(alimmode(1:16),'B'))
  if iscramb>=0
	  wave_put(iii('\DRAW_REFS:REF_',Aoutk(k)),efwave(:,k),toh )
  else
	  wave_put(iii('\DRAW_REFS:REF_',Aoutk(k)),efcur(:,k),toh )
  end
end

% Special options for E and F
if     ismember(mvloop,[1 2 12]) || rem(mvloop,10) == 5 && iscramb > 0
	wave_put('\DRAW_REFS:REF_001',repmat(-1500,1,2),t)
elseif ismember(mvloop,[11 21]) && iscramb > 0
	wave_put('\DRAW_REFS:REF_001',repmat(relovo,1,2),t)
elseif ierat == 4 && iscramb > 0
	wave_put('\DRAW_REFS:REF_001',efcur(:,17),toh)
end
if mvloop == 2 && iscramb > 0
	wave_put('\DRAW_REFS:REF_002',zeros(1,2),t)
end
if midplan ~= 8 && ikriz == 2 && nfast ~= 0 && inova == 4 && ierat~=-1 && iscramb >= 0
	wave_put('\DRAW_REFS:REF_003',efwave(:,3),toh)
end
if rem(mvloop,10) == 5 && mvloop > 100
	wave_put('\DRAW_REFS:REF_004',aipoft,toft)
end
if mvloop == 18 && iscramb > 0
	for k = 1:8, wave_put(iii('\DRAW_REFS:REF_',Aoutk(k)),zeros(1,2),t), end
end
if mvloop == 6 && iscramb > 0
 error('UNTRANSLATED')
end

% Temporary for IOH control, for runaways and more generally
if mvloop(1)==1&&(istop(1)<0||istop(1)==666)
 if ierat==4
  i1 = find(toh<istop(2)-0.001); i2 = find(toh>istop(3)+0.001);
  wave_put('\DRAW_REFS:REF_001',[efcur(i1,17);istop(4);istop(5);efcur(i2,17)],[toh(i1);istop(2);istop(3);toh(i2)])
 else
  wave_put('\DRAW_REFS:REF_001',[istop(4);istop(5)],[istop(2);istop(3)])
 end
end

if iscramb == 0,   wave_put('\DRAW_REFS:REF_018',efcur (:,17)  ,toh );
elseif iscramb >0, wave_put('\DRAW_REFS:REF_018',efwave(:,18)  ,toh ); % radial
                   wave_put('\DRAW_REFS:REF_023',efwave(:,23)  ,toh ); % vertical
 switch midplan
  case 8,          wave_put('\DRAW_REFS:REF_012',efwave(:,20)  ,toh );
                   wave_put('\DRAW_REFS:REF_022',efwave(:,22)  ,toh ); 
  case 9,          wave_put('\DRAW_REFS:REF_012',zeref         ,toh );
                   wave_put('\DRAW_REFS:REF_022',efwave(:,22)  ,toh );
                   wave_put('\DRAW_REFS:REF_024',efwave(:,24)  ,toh );
  otherwise,       wave_put('\DRAW_REFS:REF_012',zeref         ,toh ); % vertical
                   wave_put('\DRAW_REFS:REF_022',cayref        ,toh ); % elongation
                   wave_put('\DRAW_REFS:REF_024',efwave(:,24)  ,toh ); % vertical (fast pre-differentiated)
 end
end

ntoft = numel(toft);
% Ip and DIOH references
if load_user_refs
  % custom loaded values
  wave_put('\DRAW_REFS:REF_013',aipcustom     ,toh ); % Ip
  wave_put('\DRAW_REFS:REF_017',dioh          ,toh ); % DOH
else
  % default
  wave_put('\DRAW_REFS:REF_013',aipoft        ,toft); % Ip
  wave_put('\DRAW_REFS:REF_017',dioh          ,toft); % DOH
end


if ierat == -1 || load_user_refs
  % With IERAT=-1 or ISCRAMB<0, OH1 current is on observer 22
  wave_put('\DRAW_REFS:REF_022',efcur (:,17)  ,toh ); % OH current
  if load_user_refs
    if rm_den_ip_ipz_refs
      % zeroing density and Ipz references for disr detector (non-plasma shot)
      wave_put('\DRAW_REFS:REF_012',zeros(ntoft,1),toft); % zIp
      wave_put('\DRAW_REFS:REF_021',zeros(ntoft,1),toft); % density
    end
  end
end

%% Timing and LIUQE
dow = (shot == -1); % actually do the write only if this is true

% Acquisition references
k = find(abs(diff(aipoft)) < 10 & abs(aipoft(1:end-1)) > max(abs(aipoft))/2);
if isempty(k), [~,k1] = max(aipoft); k2 = k1;
else,             k1  = min(k)     ; k2 = max(k)+1;
end
vdbw(dow,strcat('SHOTDB::SHOT_DESIGN:',{'IP' 'POL' 'BPHI' 'FLAT' 'RAMP'},'_START_REF'),[        0 toh(  1) tphi(  1) toft(k1) toft( 4)])
vdbw(dow,strcat('SHOTDB::SHOT_DESIGN:',{'IP' 'POL' 'BPHI' 'FLAT' 'RAMP'},'_STOP_REF' ),[toft(end) toh(end) tphi(end) toft(k2) toft(k1)])

% PID integrator unclamp
vdbw(dow,strcat('TIMER3_OPERDB::TIMER_01',{'9' 'A'},':ABSOLUTE'),[-5 -4])
if     ierat  == 4 || ierat  == -1, tI = toh([1,end]); % Full powered window [POL_START_REF, POL_STOP_REF]
elseif mvloop ~= 6, tI = [ 0.3 2.5]; end
vdbw(dow,strcat('TIMER3_OPERDB::TIMER_01',{'B' 'C'},':ABSOLUTE'),tI     )

% LIUQE setting
pcsliuqe(dow);

end

%% PCS utils
function pcsmatput(id,X,Xin,Xout,Xscale,Xg)
s = size(X);
assert(s(2) == length(Xin) && s(1) == length(Xout) && s(1) == length(Xscale),'pcsmatput %s errors',id)
mdsw(['\PHYS_MAT_' id              ],X)
mdsw(['\PHYS_MAT_' id '_INPUTS'    ],Xin)
mdsw(['\PHYS_MAT_' id '_OUTPUTS'   ],Xout)
mdsw(['\PHYS_MAT_' id '_OUT_SCALES'],Xscale)
if strcmp(id,'A')
 mdsw('\FEED_GAINS',Xg)
 Xg(kkk(Xout,3:5)) = Xg;
 mdsw('\LOAD_FEED_GAINS',Xg)
end
end

function pcssetref(Aout,mode)

kout = kkk(Aout,3:5);
off = ~cellfun(@isempty,strfind(Aout,'!'));
refs = repmat({'XXX'},1,32);
refs(kout) = mode;

mdsw('\LOAD_WAVE_GEN_B:REFS','BUILD_SIGNAL([*,$1],,1..32)','x',refs,int8(1:32));
for k = 1:length(kout)
 mdstcl(['SET NODE \DRAW_REFS:REF_' iii(k) ' '       ifc(off(k),'/OFF','/ON')])  
 mdsw  ([         '\DRAW_REFS:REF_' iii(k) ':LABEL'],ifc(off(k),'UNUSED',['REF_' Aout{k} ' ' mode{k}]))
end

end

function pcsmatdo(index,pcsconfig)
if nargin < 2, pcsconfig = 'hybrid'; end
pcsconfig = lower(pcsconfig);

% A matrix
Ax    = mdsr('\PHYS_MAT_A');
Ain   = mdsr('\PHYS_MAT_A_INPUTS');
Aout  = mdsr('\PHYS_MAT_A_OUTPUTS');
Aoutf = mdsr('\PHYS_MAT_A_OUT_SCALES');

% Channel 13 must be positive for the disruption detector
% Ip   =	mdsr('\DRAW_REFS:REF_013');
% k = find(kkk(Aout,3:5) == 13,1);
% if abs(min(Ip)) > abs(max(Ip)) && Aoutf(k) > 0 || ...
%    abs(min(Ip)) < abs(max(Ip)) && Aoutf(k) < 0
%  Aoutf(k) = -Aoutf(k);
%  mdsw('\PHYS_MAT_A_OUT_SCALES',Aoutf)
% end

% add Ip and Bphi backoffs; store actual A matrix
[Ain,Ax] = pcsaddip(Ain,Aout,Ax);
[Ain,Ax] = pcsaddbo(Ain,Aout,Ax);
mdsw('\PHYS_MAT_A',Ax)
mdsw('\PHYS_MAT_A_INPUTS',Ain)

% calculate and store the inscale parameters
Ainf = pcsinscale(Ain);
if mdsr('$SHOT') == -1 % why is this in HYBRID tree, not PCS ?
 mdsopen('HYBRID',-1)
 mdsw('\hybrid::mat_a_in_scale',Ainf)
end

% loadable A matrix inputs and outputs (5 significant characters of Aout)
Yin = mdsu('HYBRID','\MAT_A1_CABLING');
Yout = lstmk(iii('CP',1:24));
Aout = char(Aout); Aout = cellstr(Aout(:,1:5));

% compute and store loadable A matrix
[Ay,outs,~,Aoutk] = pcsmatarrange(Ain,Aout,Ainf,Aoutf,Yin,Yout,Ax,'A');
matw('A1',Ay(:, 1: 48,:),index)
matw('A2',Ay(:,49: 96,:),index)
matw('A3',Ay(:,97:128,:),index)
mdsw('\MAT_A_OUT_SCALES',outs)

% modify the A-matrix output scales for the gains
Aoutf = Aoutf ./ mdsr('\FEED_GAINS');

% G matrix
Gin = 	Aout;
Yin =  lstmk(Yout,repmat('dummy for size',8,1));
Yout = lstmk(iii('GO',1:24));

for kG = '123'
 Gx = 	mdsr(['\PHYS_MAT_G' kG]);
 Gout = 	upper(mdsr(['\PHYS_MAT_G' kG '_OUTPUTS']));
 Gout = char(Gout); Gout = cellstr(Gout(:,1:5));
 Goutf = 	mdsr(['\PHYS_MAT_G' kG '_OUT_SCALES']);
 PIDg = mdsu('HYBRID',['\PID_' ifc(kG=='1','I',ifc(kG=='2','P','D'))]);
 PIDg(PIDg == 0) = 1;
 Ginf = Aoutf ./ PIDg(Aoutk);
 Yx = pcsmatarrange(Gin,Gout,Ginf,Goutf,Yin,Yout,Gx,['G' kG]);
 matw(['G' kG],-Yx,index)
end

% M Matrix
Mx   = mdsr('\PHYS_MAT_M');
Min  = mdsr('\PHYS_MAT_M_INPUTS');
Mout = mdsr('\PHYS_MAT_M_OUTPUTS');
Yin  = upper(lstmk(mdsu('HYBRID','\MAT_M:IN_CABLING'),repmat('DUMMY FOR SIZE',8,1)));
Yout = mdsu('HYBRID','\MAT_M:OUT_CABLING');
Moutf = pcsoutscale(Yout);
[~,k1,k2] = intersect(Min,Ain);  Minf(k1) = Ainf (k2);
[~,k1,k2] = intersect(Min,Gout); Minf(k1) = Goutf(k2);

switch pcsconfig
 case 'rtdtacqden', Minf(21                  ) = 1; Moutf(   19        ) = 1; Mx(19,21) = 1;
 case 'rtdtacq',    Minf(~strcmp(Min,'GO022')) = 1; Moutf([1:19 21:end]) = 1; Mx(1:19,[1:21 23:40]) = eye(19,39); Mx(19,[19 21])= [0 1];
 case 'rtdtacqfps', Minf(:                   ) = 1; Moutf(  :          ) = 1; Mx(:,[1:18 21:40]) = eye(size(Mx,1),38);
end

Yx = pcsmatarrange(Min,Mout,Minf,Moutf,Yin,Yout,Mx,'M');
matw('M',Yx,index)

function matw(name,x,index)
 mdsw(['\LOAD_MAT_' name],x)
 if ismatrix(x)
  mdsw(['\LOAD_MAT_' name '_ADDRESS'],'[$]','x',int16(index))
 else
  mdsw(['\LOAD_MAT_' name '_ADDRESS'],int16((0:size(Yx,3)-1)+index))
 end
end

end

function [Ain,Ax] = pcsaddip(Ain,Aout,Ax)

MU0 = 4e-7*pi;

Ain  = upper(Ain);
Aout = upper(Aout);
kout = find(~cellfun(@isempty,strfind(Aout,'IPA'))); 
if isempty(kout), return, end

% BPOL contributions
cin = lstmk(strcat('CTLINT:BPOL_AVG_',mdsu('MAGNETICS','DIM_OF(\IPLASMA:TRAPEZE:COEFF)')));
c   =                                 mdsu('MAGNETICS',       '\IPLASMA:TRAPEZE:COEFF') / MU0;
Ain = lstun(Ain,cin);
Ax(kout,lstix(cin,Ain),:) = repmat(c',[length(kout) 1 size(Ax,3)]);

% IPOL contributions
cin = lstmk(strcat('CTLAMP:IPOL_',mdsu('MAGNETICS','DIM_OF(\IPLASMA:TRAPEZE:CROSSTALK)')));
c   =                             mdsu('MAGNETICS',       '\IPLASMA:TRAPEZE:CROSSTALK');
Ain = lstun(Ain,cin);
Ax(kout,lstix(cin,Ain),:) = repmat(-c',[length(kout) 1 size(Ax,3)]);

end

function [Ain,Ax] = pcsaddbo(Ain,~,Ax)

Ain  = upper(Ain);

% add IPHI
cin = lstmk('CTLAMP:IPHI_001');
Ain = lstun(Ain,cin);

% get bo coefficients
c 	= zeros(length(Ain),1);
[k,l] = lstix(lstmk(iii('CTLINT:BPOL_AVG_',1:38)),Ain);
if ~isempty(k), c(k) = mdsu('MAGNETICS','(\BPOL_003:CROSSTALK[$1]+\BPOL_011:CROSSTALK[$1])/2',iii(l)); end
[k,l] = lstix(lstmk(iii('CTLINT:RVLOOP_',1:38)),Ain);
if ~isempty(k), c(k) = mdsu('MAGNETICS','    \RFLUX:CROSSTALK[$1]                           ',iii(l)); end
[k,l] = lstix(lstmk(iii('CTLINT:VLOOP_',1:38)),Ain);
if ~isempty(k), c(k) = mdsu('MAGNETICS','     \FLUX:CROSSTALK[$1]                           ',iii(l)); end

% apply bo
k = lstix(cin,Ain);
Ax(:,k,:) = 0;
for l = 1:size(Ax,3), Ax(:,k,l) = Ax(:,:,l)*(-c); end

end

function gains = pcsoutscale(out)

gains = zeros(length(out),1);
for k = 1:length(out)
 switch out{k}
  case lstmk(iii('E_' ,1:8)), gains(k) = mdsu('HYBRID','\OUTPUT_SCALES_E_U'   );
  case lstmk(iii('F_' ,1:8)), gains(k) = mdsu('HYBRID','\OUTPUT_SCALES_F_U'   );
  case lstmk(iii('OH_',1:2)), gains(k) = mdsu('HYBRID','\OUTPUT_SCALES_OH_U'  );
  case 'FAST'        , gains(k) = mdsu('HYBRID','\OUTPUT_SCALES_FAST_U');
  case 'GAS'         , gains(k) = 1;
  case strcat('RHVPS',{'A' 'B' 'C'}), gains(k) = mdsr(['PRODUCT(\RHVPS_GAIN_' out{k}(6) '[*,"PCS"],0)']);
  case strcat({'MIRROR'},num2str(1:6,'%1d')'), gains(k) = mdsu('TCV_SHOT',['\PCS::MIR_V_SLOPE_0' out{k}(7)]);
  otherwise, error('Invalid PCS output %s',out{k})
 end
end
end

function invgains = pcsinscale(channels)

% real-time magnetic gains
% names, relative chixes and labels from the tree, select usefull channels
magn = mdsu('VSYSTEM','\MAGDB_I:NAMES');
magn = magn(~cellfun(@isempty,regexp(magn,'(^PRE|^AMP|^INT)')));
a = regexp(mdsu('VSYSTEM','\MAGDB_L[$]',magn),'\d+/\d+/(?<chan>\d+)/(?<card>\d+)','names');
magv = vdbr(strcat('MAGDB::',magn));
magc = arrayfun(@(x)str2double(x.chan),[a{:}])';
magp = arrayfun(@(x)cellstr(x.card)   ,[a{:}])';
magn = regexprep(magn,'_...GAIN$','');

% calibrated gains
g = zeros(length(magn),1);
k = ~cellfun(@isempty,regexp(magn,'^PRE'));
for p = unique(magp(k))'
 l = k & strcmp(magp,p);
 tmp = mdsu('MAGNETICS',['\PRE:CD_' p{1}]);
 g(l) = tmp(sub2ind(size(tmp),magc(l),magv(l)+1));
end
k = ~cellfun(@isempty,regexp(magn,'^AMP'));
for p = unique(magp(k))'
 l = k & strcmp(magp,p);
 tmp = mdsu('MAGNETICS',['\AMP:CD_' p{1}]);
 g(l) = -0.9969*tmp(sub2ind(size(tmp),magc(l),magv(l)+1));
end
k = ~cellfun(@isempty,regexp(magn,'^INT'));
for p = unique(magp(k))'
 l = k & strcmp(magp,p);
 tmp = mdsu('MAGNETICS',['\INT:CD_' p{1}]);
 g(l) = -0.9969./tmp(sub2ind(size(tmp),magc(l),rem(11-magv(l),8)+1));
end

% magnetics calibration
sens_bpol_003 = mdsu('MAGNETICS','\BPOL_003:AREA*\POL_PROB_CHAR:POLARITY["003",*]');
sens_bpol_011 = mdsu('MAGNETICS','\BPOL_011:AREA*\POL_PROB_CHAR:POLARITY["011",*]');
escale        = mdsu('MAGNETICS','\IPOL:CALIBRATION["E_001" ]');
fscale        = mdsu('MAGNETICS','\IPOL:CALIBRATION["F_001" ]');
gscale        = mdsu('MAGNETICS','\IPOL:CALIBRATION["G_001" ]');
ohscale       = mdsu('MAGNETICS','\IPOL:CALIBRATION["OH_001"]');
torscale      = mdsu('MAGNETICS','\IPHI:CALIBRATION')*mdsr('\MGAMS.DATA:IF36FB');

assert(isequal(sign(sens_bpol_003),sign(sens_bpol_011)),'Polarity must be identical in sector 3 and 11')

invgains = zeros(length(channels),1);
for i = 1:length(channels)

 name = channels{i}(4:end);

 switch name(1:min(end,8))
  case { 'AMP:IPOL' 'AMP:IPHI' }
   gain = g(strcmp(name,magn));
   switch name(1:11)
    case 'AMP:IPOL_E_', gain = gain / escale;
    case 'AMP:IPOL_F_', gain = gain / fscale;
    case 'AMP:IPOL_G_', gain = gain / gscale;
    case 'AMP:UPOL_G_', gain = gain / gVscale;
    case 'AMP:IPOL_OH', gain = gain / ohscale;
    case 'AMP:IPOL_DO', gain = gain / ohscale;
    case 'AMP:IPHI_00', gain = gain / torscale;
    otherwise
     error('Invalid IPOL inside PCS_INVGAINS')
   end

  case { 'AMP:RVLO' 'AMP:VLOO'}
   gain = g(strcmp(name,magn));

  case { 'INT:RVLO' 'INT:VLOO'}
   gain = g(strcmp(name,magn));
  case { 'INT:ROG_'}
   error('UNTRANSLATED')
 
  case 'INT:BPOL'
   tmp = name(14:16);
   gain_003 = g(strcmp(['PRE:BPOL_003_' tmp],magn)) * g(strcmp(['INT:BPOL_003_' tmp],magn));
   gain_011 = g(strcmp(['PRE:BPOL_011_' tmp],magn)) * g(strcmp(['INT:BPOL_011_' tmp],magn));
   tmp = str2double(tmp);
   gain = (sens_bpol_003(tmp)*gain_003 + sens_bpol_011(tmp)*gain_011)/2;

  case 'AMP:BPOL'
   tmp = name(14:16);
   gain_003 = g(strcmp(['PRE:BPOL_003_' tmp],magn)) * g(strcmp(['AMP:BPOL_003_' tmp],magn));
   gain_011 = g(strcmp(['PRE:BPOL_011_' tmp],magn)) * g(strcmp(['AMP:BPOL_011_' tmp],magn));
   tmp = str2double(tmp);
   gain = (sens_bpol_003(tmp)*gain_003 + sens_bpol_011(tmp)*gain_011)/2;

  otherwise
   switch channels{i}
    case 'FIR_FRINGE'
     gain = 2^(vdbr('DIAG2DB::FIR:FRINGE_GAIN')-8) / 0.1075;
     
    case '1MM_FRINGE'
     gain =  2^(vdbr('DIAG2DB::TWOMM:FRINGE_GAIN')-8) / (-0.1); % mdsu('BASE','\ONEMM:FRINGE_COUNT:CONVERSION');
     
    case '2MM_FRINGE'
     error('UNTRANSLATED')
     
    otherwise
     gain = 1;
     tmp = channels{i};
     assert(strcmp(tmp(1:2),'GO'),['A-Input not found : ' channels(i)])
   end
 end
 
 invgains(i) = 1/gain;

end

end

function [Yx,Youtf,Yink,Youtk] = pcsmatarrange(Xin,Xout,Xinf,Xoutf,Yin,Yout,Xx,~)

[nout,nin,nsw] = size(Xx);
assert(length(Xin) == nin && length(Xinf) == nin && length(Xout) == nout && length(Xoutf) == nout,'MAT ARRANGE - Bad Input data dimensions')

% render case-insensitive
Xin  = deblank(upper(Xin));
Xout = deblank(upper(Xout));
Yin  = deblank(upper(Yin));
Yout = deblank(upper(Yout));

% find the indices, Xin == Yin(Yink), Xout == Yout(Youtk)
Yink  = lstix(Xin ,Yin );
Youtk = lstix(Xout,Yout);

% output scales
Youtf = NaN(length(Yout),1);
Youtf(Youtk) = Xoutf;

% arranged and scaled matrix
Yx = zeros(length(Yout),length(Yin),nsw);
Yx(Youtk,Yink,:) = Xx./Xoutf(:).*Xinf(:)';
end

function pcsmatoverflow(shot)
if shot ~= -1, return, end
for m = {'A1','A2','A3','G1','G2','G3','M'}
 tmp = mdsu('\HYBRID',['\MAT_' m{1} ':COUNTS']);
 ii = find(tmp >= 4095 | tmp <= 0 | isnan(tmp));
 if ~isempty(ii)
  % Display list of offending indices
  [i1,i2,i3] = ind2sub(size(tmp),ii);
  % Skip 3rd index if unused
  if ismatrix(tmp), tmp2 = [i1,i2];
  else,             tmp2 = [i1,i2,i3]; end
  indexlist = strjoin(cellstr(num2str(tmp2)),'\n');
  assert(isempty(ii),'MAT_OVERFLOW in %s\n  OFFENDING INDICES:\n%s',m{1},indexlist)
 end
end
end

function [t,v] = pcssignbit(t,v)
dt1 = 0.004; dt2 = 0.009;
v0 = sign(v(find(v,1))); if isempty(v0), v0 = 1; end
v(1:3) = 0;
v = sign(v); v(v == 0) = v0;
v1 = diff(v);
k  = find(v1);
t1 = t([1;k+1;end]);
v1 = [2*v0;v1(k);-2*v(end)];
% Remove pairs of sign bits closer than dt2 as they would be of opposite
% sign by design and thus have no real effect. Start with the last pair to
% ensure that in case of a triplet, the first sign bit is kept.
while numel(t1)>1
 i1 = find(diff(t1) <= dt2,1,'last');
 if isempty(i1); break; end
 t1(i1:i1+1) = [];
 v1(i1:i1+1) = [];
end
t1    = t1    - dt2;
t1(1) = t1(1) + dt2 -0.030;
t = sort([t1;t1+dt1;t1-dt1;t(end)]);
v = zeros(length(t),1); v(2:3:end-1) = 500*v1;

end

function pcssetalim(mode,lev,sel,sbit)

% force TOR
lev{19} = '100';
mode{19} = 'A';

% check requests for acceptability
mode = upper(mode);
assert(all(ismember(mode     ,{'A' 'B' 'C'}         ))                            ,'Invalid modes'  )
assert(all(ismember(lev(1:19),{'20' '50' '80' '100'})) && str2double(lev{20}) >=0 ,'Invalid levels' )
assert(all(ismember(sel      ,{'ON' 'OFF'}          ))                            ,'Invalid select')

% write to tree
                mdsw('\ALIM_MODES'    ,mode)
                mdsw('\ALIM_LEVELS'   ,lev )
                mdsw('\ALIM_SELECTS'  ,sel )
if nargin == 4, mdsw('\LOAD_SIGN_BITS',sbit), end

end

function pcsliuqe(dow)

try
 r = mdsu('HYBRID','tcv_eq("R_CONTOUR","FBTE")'); z = mdsu('HYBRID','tcv_eq("Z_CONTOUR","FBTE")');
 k = find(r == 0); r(k) = NaN; z(k) = NaN;
 r1 = min(r(:)); r2 = max(r(:)); z1 = min(z(:)); z2 = max(z(:));
 n=max(3,min(2*round(2*max((max(z)-min(z))./(max(r)-min(r)))),12));
catch
 r1 =  0.63; r2 = 1.12; z1 = -0.73; z2 = 0.73;
 n = 12;
end

% Liuqe 1,2,3
if dow
 mdsopen('RESULTS',-1)
 for k = strcat('\PARAMETERS',{':' '_2:' '_3:'})
  mdsw([k{1} 'NELEM'],int8(n))
  mdsw([k{1} 'RIP'  ],r1)
  mdsw([k{1} 'ROP'  ],r2)
  mdsw([k{1} 'ZLP'  ],z1)
  mdsw([k{1} 'ZUP'  ],z2)
 end
end

end

%% MDS utils
function x = mdsr(e,varargin)
[x,s] = mdsvalue(e,varargin{:});
assert(logical(rem(s,2)),'Cannot evaluate MDS expression %s\n%s',e,mdsdata('GETMSG($1)',s))
if iscellstr(x), x = deblank(x); end
end

function x = mdsu(tree,e,varargin)
e = ['USING(DATA(' e '),\TOP,-1,"' tree '")'];
x = mdsr(e,varargin{:});
end

function mdsw(nodes,varargin)
if iscellstr(varargin{1}) && isvector(varargin{1})
 mdsw(nodes,'SET_RANGE($1,$2)','x',length(varargin{1}),varargin{1});
elseif isfloat(varargin{1})
 mdsw(nodes,mdscvt(varargin{1},'f'))
else
 for node = cellstr(nodes)
  s = mdsput(node{1},varargin{:});
  assert(logical(rem(s,2)),'Can not put MDS node %s\n%s',node{1},mdsdata('GETMSG($1)',s))
 end
end
end

function wave_put(node,x,t,o)
if nargin == 3
 s = mdsr('WAVE_PUT($1,$2,$3)',node,x,t);
%  xx = mdsr(node); tt = mdsr(['DIM_OF(' node ')']);
%  if s == 0 && (~isequal(x(:),xx) || numel(t) ~= numel(tt) || max(abs(t(:)-tt(:))) >= 1e-4)
%   warning('WAVE_PUT failed %s',node)
%   mdsw([node ':TRACK_000'],'BUILD_SIGNAL($1,,$2)','x',x,t)
%  end
elseif o == 8
	s = mdsr('WAVE_PUT($1,[REPLICATE(SET_RANGE(,1,F_FLOAT($2)),1,7),ZERO(SHAPE($2),F_FLOAT(0.0))],F_FLOAT($3))',node,x,t);
end
assert(s==0,'WAVE_PUT failed (error %d) for node %s',s,node)
end

function wave_put_scd(node,xo,to)
[t,x] = wave_compress(to,xo); % compress points eliminating linearly interpolatable points
mdsw([node,':TIME'],t);
for ii=0:7
  mdsw([node,':TRACK_',iii(ii)],['BUILD_SIGNAL($1,,BUILD_DIM(*,',node,':TIME))'],'x',mdscvt(x,'f'));
end
end

function v = vdbr(c)
v = mdsr('VDBGET($)',c);
end

function vdbw(do,c,v)
if do, mdsr('VDBPUT($1,$2)',c,mdscvt(v,'f')); end
end

%% Utils
function a = iii(f,n)
if nargin == 1, n = f; f = ''; end
a = num2str(n(:),[strrep(f,'\','\\') '%03d']);
end

function k = kkk(a,l)
k = char(a(:)); k = str2num(k(:,l)); %#ok<ST2NM>
end

function a = lstmk(varargin)
for k = 1:nargin
 if ischar(varargin{k})
  varargin{k} = cellstr(varargin{k})';
 else
  varargin{k} = varargin{k}(:)';
 end
end
a = [varargin{:}];
end

function [k,l] = lstix(x,y)
% find k such that x == y(k) or x(l) = y(k)
[~,l,k] = intersect(x,y); [l,m] = sort(l); k = k(m);
assert(nargout == 2 | length(k) == length(x),'x not found in y')
end

function z = lstun(x,y)
[~,~,k] = union(x,y,'stable');
z = lstmk(x,y(k));
end
