function [varargout] = tcv_recipe(recipetype,varargin)
% TODO: message and time of start (whole program, individual subroutines)
% TODO: message and time of completion (whole program, individual subroutines)
% TODO: error handling [review old recipes; errors for wrong number of
%   arguments, etc. - or only more serious ones?]
%       "Gates, bad vcs_bput"

% Master switchyard
% TODO: revise next line depending on scheme chosen for daemon
global tcv_recipe_error
if nargin == 0, recipetype = 'fbte'; end
if strcmp(recipetype,'fbte') && nargin>1 && varargin{1}~=-1,return;end
tcv_recipe_error=0;

recipes = {'fbte','stray','ecrhonly','nbhonly','backoff', ...
        'fpsonly','gasonly','diagonly','pcs_test','timer_only'};
recipetype = validatestring(recipetype,recipes);

% TODO: understand why "mdsdisconnect" is necessary after running the recipe
%   or the following one hangs
%   indefinitely at the matrix r_download: this is true for tcv_oper but *not*
%   for tcv_root
mdsdisconnect;
mdsopen('PCS',-1)
vdbw('tcv_protectdb::shot_design:pending','on')
vdbw('tcv_protectdb::shot_create:command_1',recipetype)
switch recipetype
  case 'timer_only'
    pcstimeronly;
  case 'pcs_test'
    pcspcsdesign;
  case 'stray'
    pcsstraydesign(varargin{:});
  case 'backoff'
    pcsbackoff(varargin{:});
  case 'fbte'
    disp('Blank recipe started');
end
pcslimits
% must do this better
ecrhwasntready=false;
ecrhready = vdbr('TCV_ECRH_X2DB::ECRH:READINESS');
ecrhreq = ...
     strcmpi(deblank(vdbr('TCV_ECRH_X2DB::ECRH_ACQ:SELECT')),'ON') | ...
     strcmpi(deblank(vdbr('TCV_ECRH_X3DB::ECRH_ACQ:SELECT')),'ON') ;
while ecrhreq && (ecrhready<1 || ecrhready>3)
  ecrhwasntready=true;
  vdbw('TCV_PROTECTDB::CYCLE:STATE','CREATE WAITING ECRH')
  disp('Waiting for ECRH')
  pause(2.5)
  ecrhready = vdbr('TCV_ECRH_X2DB::ECRH:READINESS');
  ecrhreq = ...
     strcmpi(deblank(vdbr('TCV_ECRH_X2DB::ECRH_ACQ:SELECT')),'ON') | ...
     strcmpi(deblank(vdbr('TCV_ECRH_X3DB::ECRH_ACQ:SELECT')),'ON') ;
end
if ecrhwasntready
  disp('ECRH is now ready')
  vdbw('TCV_PROTECTDB::CYCLE:STATE','')
end
nbh1wasntready=false;
nbh2wasntready=false;
nbh1ready = vdbr('TCV_NBHDB::NBI:READINESS');
nbh2ready = vdbr('TCV_NBHDB::NB2:READINESS');
nbh1req = strncmpi(vdbr('tcvpcdb::nbh:reqfordoshot'),'ON',2) && ...
  vdbr('TCV_NBHDB::NBI:LCS_MODE') == 9;
nbh2req = strncmpi(vdbr('tcvpcdb::nb2:reqfordoshot'),'ON',2) && ...
  vdbr('TCV_NBHDB::NB2:LCS_MODE') == 9;
while nbh1req && (nbh1ready<1 || nbh1ready>3)
  nbh1wasntready=true;
  vdbw('TCV_PROTECTDB::CYCLE:STATE','CREATE WAITING NBH 1')
  disp('Waiting for NBH 1')
  pause(2.5)
  nbh1ready = vdbr('TCV_NBHDB::NBI:READINESS');
  nbh1req = strncmpi(vdbr('tcvpcdb::nbh:reqfordoshot'),'ON',2) && ...
    vdbr('TCV_NBHDB::NBI:LCS_MODE') == 9;
end
if nbh1wasntready
  disp('NBH 1 is now ready')
  vdbw('TCV_PROTECTDB::CYCLE:STATE','')
end
while nbh2req && (nbh2ready<1 || nbh2ready>3)
  nbh2wasntready=true;
  vdbw('TCV_PROTECTDB::CYCLE:STATE','CREATE WAITING NBH 2')
  disp('Waiting for NBH 2')
  pause(2.5)
  nbh2ready = vdbr('TCV_NBHDB::NB2:READINESS');
  nbh2req = strncmpi(vdbr('tcvpcdb::nb2:reqfordoshot'),'ON',2) && ...
    vdbr('TCV_NBHDB::NB2:LCS_MODE') == 9;
end
if nbh2wasntready
  disp('NBH 2 is now ready')
  vdbw('TCV_PROTECTDB::CYCLE:STATE','')
end
pcssetup
pcsecrhsetup
pcsnbhsetup
%pcseclaunrecsetup
pcscondsettings
pcsalimgates
mdsvalue('timing_setup()');
pcsmatswitching
pcsspectimers
pcsmdsprepare
% temporary to avoid problems with PREF/ECE
wave_put('\DRAW_FEEDFOR_ECRH:POL_MIR_003',[0 0],[-1 3])
timenow=[' ',char(datetime('now','Format','HH:mm')),' '];
%if strcmp(recipetype,'stray'),
% close primary valves for NEG test
% button_with_timeout('tcv_protectdb','pompe_no:command_8', ...
%                     'tcv_protectdb::pompe_no:state_b7',0,10,1.5)
% button_with_timeout('tcv_protectdb','pompe_ne:command_8', ...
%                     'tcv_protectdb::pompe_ne:state_b7',0,10,1.5)
% button_with_timeout('tcv_protectdb','pompe_se:command_8', ...
%                     'tcv_protectdb::pompe_se:state_b7',0,10,1.5)
% button_with_timeout('tcv_protectdb','pompe_so:command_8', ...
%                     'tcv_protectdb::pompe_so:state_b7',0,10,1.5)
%else
% button_with_timeout('tcv_protectdb','pompe_no:command_7', ...
%                     'tcv_protectdb::pompe_no:state_b7',1,10,1.5)
% button_with_timeout('tcv_protectdb','pompe_ne:command_7', ...
%                     'tcv_protectdb::pompe_ne:state_b7',1,10,1.5)
% button_with_timeout('tcv_protectdb','pompe_se:command_7', ...
%                     'tcv_protectdb::pompe_se:state_b7',1,10,1.5)
% button_with_timeout('tcv_protectdb','pompe_so:command_7', ...
%                     'tcv_protectdb::pompe_so:state_b7',1,10,1.5)
%end
recipestring=recipetype;
vdbw('tcv_protectdb::comment:text',recipestring)
vdbw('tcv_protectdb::shot_design:recipe',[timenow,recipestring])
if ~tcv_recipe_error
  vdbw('tcv_protectdb::shot_design:pending','off')% Switch OFF Pending
  vdbw('tcv_protectdb::shot_create:command_1',' ')
else
  disp(['Exited with final error code ',int2str(tcv_recipe_error)])
end
varargout{1}=tcv_recipe_error;
disp('Recipe finished')

%%%%%***** function separator #####

% Timer only
function pcstimeronly

% Deselect all the power supplies (mode and level are not important)
nps=20;
pcssetalim(repmat({'A'},1,nps),repmat({'20'},1,nps),repmat({'OFF'},1,nps))
vdbw('tcv_protectdb::alims:select','OFF')
% Turn off all RHVPS and mirror motors
pcsecrhmode('OFF')
% Select and deselect
vdbw(strcat('tcv_protectdb::',{'timer','acq','pcs','gas','vide'},':select'), ...
     {'ON','OFF','OFF','OFF','OFF'}')
vdbw('tcv_publicdb::diags:select','OFF')
% Deselect auto-glow
vdbw('tcvpcdb::doglow:select_auto','OFF')

%%%%%***** function separator #####

% PCS design
function pcspcsdesign

% Common actions to no-plasma shots (PCS test and stray)
pcsnoplasma

% Deselect all the power supplies (mode and level are not important)
nps=20;
pcssetalim(repmat({'A'},1,nps),repmat({'20'},1,nps),repmat({'OFF'},1,nps))
vdbw('tcv_protectdb::alims:select','OFF')

% sinusoids in references
refs_osc

%%%%%***** function separator #####

% Stray design
function pcsstraydesign(tor_sign)
if nargin<1 || isempty(tor_sign),tor_sign=-1;end
% set consistent toroidal field sign for observers
mdsw('\data:if36fb',tor_sign)
pcsnoplasma
% toroidal field sign
mdsw('\TOR_SIGN',tor_sign)
nps=20;
pcssetalim([repmat({'C'},1,nps-1),'B'],[repmat({'50'},1,nps-1),'60'], ...
           repmat({'ON'},1,nps),ones(1,nps-2))
% zero references
for k = 1:24, wave_put(['\DRAW_REFS:REF_',iii(k)],[0,0],[0,3]); end


%%%%%***** function separator #####

% Prepare shot with no plasma
function pcsnoplasma

mu0=4e-7*pi;
bphirad=0.865;
bphi2iphi=2*pi*bphirad/mu0/16/6;
%ABA zeroing all perts
% Set up fbte18 observers for testing
% - Load standard shot to avoid errors
try
  fbte18(-1,108);
catch
  MG=mgp(100000);
  MG.save(-1);
  fbte18(-1,108);
  disp('tcv_recipe: Error in loading observers for current MGAMS shot: loaded standard shot observers instead')
end
mdsopen('pcs',-1)
% toroidal field
bphi = [0,0.4,0.4,0];
tbphi = [0,0.5,0.7,1.05];
wave_put('\DRAW_FEEDFOR_TOR_I:ALIM_001',bphi*bphi2iphi,tbphi)

n_oh=2;n_e=8;n_f=8;
% poloidal field
ped_oh = 800;amp_oh = 2500;
ped_ef = 400;amp_ef = 1000;
% if the sign of the first nonzero value (sign bit) changes, modify signbit in
%   straydesign accordingly
ohcur = [0,700,0,0,0,ped_oh,ped_oh+[0,1,1,-1,0]*amp_oh,ped_oh,0];
efcur = [0,300,0,0,0,ped_ef,ped_ef+[0,1,1,-1,0,0]*amp_ef,ped_ef,0];
fpsv = [0,0,-1,-1,0,0,1,1,0,0]*27;
fpsdcv = [100,[1,1]*250,100];
fpscur = [0,0,0,0];
tpol_start = 0;
tpol_init = tpol_start+[-20,-16,-12,-10,0,10]*1e-3;
toh_rel = [0,30,60,120,150]*1e-3;
toh_start = tbphi(end)+0.1;
te_rel = [0,10,44,64,70,76]*1e-3;
tf_rel = [0,14,60,84,100,106]*1e-3;
tef_start = toh_start+2*toh_rel(end)+0.04;
tfpsv = [0,1,11,12,17,18,28,29]*1e-3;
tfpsdcv = [-0.1,tfpsv(end)];
tfpscur = [0,tfpsv(end)];
tpol_end = tef_start+8*(te_rel(end)+tf_rel(end))+0.1;
tfps_start = tpol_end-0.07;
% Complete time base
tpol_tot = cell(1,n_e+n_f+2);
for k = 1:n_e/2
  tpol_tot{k}     = [tpol_init,tef_start+te_rel+(k+3)*te_rel(end)+ ...
                        8*tf_rel(end),tpol_end-0.03,tpol_end-0.01];
end
for k = n_e/2+1:n_e
  tpol_tot{k}     = [tpol_init,tef_start+te_rel+(k-5)*te_rel(end), ...
                        tpol_end-0.03,tpol_end-0.01];
end
for k = 1:n_f
  tpol_tot{n_e+k} = [tpol_init,tef_start+tf_rel+4*te_rel(end)+ ...
                        (8-k)*tf_rel(end),tpol_end-0.03,tpol_end-0.01];
end
tpol_tot{n_e+n_f+1} = [tpol_init,toh_start+toh_rel,tpol_end-0.03,tpol_end-0.01];
tpol_tot{n_e+n_f+2} = [tpol_init,toh_start+toh_rel+0.17,tpol_end-0.03,tpol_end-0.01];
%
ip_ref = mdsr('_ipref=\DRAW_REFS:REF_013');
[ya(:,1:n_e+n_f+n_oh),tc]= ...
      pcsalignipol([repmat({efcur},1,n_e+n_f),repmat({ohcur},1,n_oh)],tpol_tot);
ya(:,end+1) = interp1(mdsr('dim_of(_ipref)')',ip_ref',tc,'linear','extrap');
wave_put('\DRAW_REFS:REF_013',ya(:,end),tc)
drawstr='\DRAW_FEEDFOR_';
for k=1:n_e+n_f+n_oh
  if k <= n_e,          wpstr = [drawstr,'E_I:ALIM_',iii(k)];
  elseif k <=n_e+n_f,   wpstr = [drawstr,'F_I:ALIM_',iii(k-n_e)];
  else,                 wpstr = [drawstr,'OH_I:ALIM_',iii(k-n_e-n_f)];
  end
  wave_put(wpstr,ya(:,k),tc)
end
wave_put('\DRAW_FPS:UREF_001',fpsv, ...
         [tpol_init(1),tfps_start+tfpsv,tpol_end+0.02])
wave_put('\DRAW_FPS:UDCREF_001',fpsdcv, ...
         [tpol_init(1),tfps_start+tfpsdcv,tpol_end+0.02])
wave_put('\DRAW_FPS:IREF_001',fpscur, ...
         [tpol_init(1),tfps_start+tfpscur,tpol_end+0.02])
% remove sign bit part
ya(tc<tc(1)+0.021,:)=0;
%
pcshybcurrtovolt({'E','F','OH';n_e,n_f,n_oh},ya(:,1:n_e+n_f+n_oh),tc,1.0)
% gas
ngas=5;
gas = repmat([0,3.2,3.2,0],1,ngas);
tgas = repmat([0,1,30,31]*1e-3,1,ngas)+ ...
       reshape(repmat([0.1,0.3,0.5,0.7,0.9],length(gas)/ngas,1),1,length(gas));

% NEG test
% recipetype = lower(deblank(vdbr('tcv_protectdb::shot_create:command_1')));
%if contains(recipetype,'stray'),
%  tgas=tgas-3.5;
%end
% temporary test for FF
% recipetype = lower(deblank(vdbr('tcv_protectdb::shot_create:command_1')));
%if contains(recipetype,'stray'),
%  gas = [0.000 0.000 2.000 2.000 0.000 0.000 0.000 2.000 2.000 0.000 0.000 1.000 1.000 1.500 1.500 0.500 0.500 0.500 0.500 0.500 1.500 0.500 0.500 1.500 0.500 0.500 0.500 0.500 0.500 1.500 1.500 0.500 0.500 0.500 0.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 0.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 1.500 0.000];
%  tgas = [-3.970 -3.951 -3.950 -3.930 -3.929 -3.000 -2.981 -2.980 -2.960 -2.959 -1.500 -1.499 -1.400 -1.390 -1.380 -1.370 -1.360 -1.350 -1.340 -1.330 -1.320 -1.310 -1.300 -1.290 -1.280 -1.270 -1.260 -1.250 -1.240 -1.230 -1.220 -1.210 -1.200 -1.190 -1.180 -1.170 -1.160 -1.150 -1.140 -1.130 -1.120 -1.110 -1.100 -1.090 -1.080 -1.070 -1.060 -1.050 -1.040 -1.030 -1.020 -1.010 -1.000 -0.990 -0.980 -0.970 -0.960 -0.950 -0.940 -0.930 -0.920 -0.910 -0.900 -0.890 -0.880 -0.870 -0.860 -0.850 -0.840 -0.830 -0.820 -0.810 -0.800 -0.790 -0.780 -0.770 -0.760 -0.750 -0.740 -0.730 -0.720 -0.710 -0.700 -0.690 -0.680 -0.670 -0.660 -0.650 -0.640 -0.630 -0.620 -0.610 -0.600 -0.590 -0.580 -0.570 -0.560 -0.550 -0.540 -0.530 -0.520 -0.510 -0.500 -0.490 -0.480 -0.470 -0.460 -0.450 -0.440 -0.430 -0.420 -0.410 -0.400 -0.300 -0.299];
%end
wave_put('\DRAW_FEEDFOR_GAS:ALIM_001',gas,tgas)
wave_put('\DRAW_FEEDFOR_GAS:ALIM_002',[0 0],tgas([1 end]))
wave_put('\DRAW_FEEDFOR_GAS:ALIM_003',[0 0],tgas([1 end]))
% timing references
t_start = min(tbphi(1),tpol_start);
t_stop = max(tbphi(end),tpol_end);
vdbw('tcv_protectdb::shot_design:ip_end_time',1.0)
vdbw(strcat('tcv_protectdb::shot_design:', ...
      {'ip','pol','bphi','ecrh','ecrh_a','ecrh_b','ecrh_c','flat','ramp'}, ...
        '_start_ref'),[t_start,t_start,t_start,0,0,0,0,0,0])
vdbw(strcat('tcv_protectdb::shot_design:', ...
      {'ip','pol','bphi','ecrh','ecrh_a','ecrh_b','ecrh_c','flat','ramp'}, ...
        '_stop_ref'),[t_stop,t_stop,t_stop,0,0,0,0,0,0]')
% reset acquisition time for PCS and magnetics
dt_start = -0.04;
dt_stop = 0.04;
acq_ch = strcat('tcv_publicdb::', ...
                [strcat('trcf_mag_',{'bpol','flux','loop','curr'},'_001'), ...
                 strcat('cadf_',{'wav_a_001','ain_a1_1'})]);
vdbw(strcat(acq_ch,':delta_t_start'),repmat(dt_start,1,length(acq_ch)))
vdbw(strcat(acq_ch(4:5),':delta_t_stop'),repmat(dt_stop,1,2))
% reset acquisition frequency for magnetics
vdbw(strcat('tcv_publicdb::trcf_mag_',{'bpol','flux','loop','curr','ffps'}, ...
            '_001:frequence'),[0,0,0,0,0])
% integrator clamp/unclamp
vdbw(strcat('timers_operdb::timer_01',{'9','A','B','C'},':absolute'), ...
     [0.05,1.0,3.0,4.0])
% set up matrix addresses, no switching
tsw=[-0.05,0.3,1.9];
wave_put('\PHYS_MAT_ADDRESSES:A',zeros(size(tsw)),tsw)
wave_put('\PHYS_MAT_ADDRESSES:G',zeros(size(tsw)),tsw)
wave_put('\PHYS_MAT_ADDRESSES:M',zeros(size(tsw)),tsw)
% zero G matrices
for kG = '123'
  mdsw(['\LOAD_MAT_G',kG],zeros(24,32,10))
  mdsw(['\PHYS_MAT_G',kG],zeros(22,24,10))
end
% Deselect auto-glow
%vdbw('tcvpcdb::doglow:select_auto','OFF')
% Turn off all RHVPS and mirror motors
pcsecrhmode('OFF')

%%%%%***** function separator #####
function pcsbackoff(mode,coilparam,tor_sign)
if nargin<2
  disp('Insufficient number of arguments')
  return
end
if nargin<3,tor_sign=-1;end
if ~iscell(coilparam)
  disp('The second argument must be a cell array')
end
if ~iscell(coilparam{1})
  coilparam={coilparam};% deal with case in which there is only one coilparam set
end
if mode==2
 if length(coilparam)>1
  disp(['With mode=2 only one coilparam component is accepted; the first ', ...
        'component will be retained']);
 end
 coilparam = coilparam{1};
end
switch mode
 case 1
  for icp=length(coilparam):-1:1
   if length(coilparam{icp})<5
    disp(['Insufficient number of components in coilparam element #', ...
          int2str(icp),': ignored'])
    coilparam(icp)=[];
   else
    if length(coilparam{icp})<6,coilparam{icp}{6} = 0;end
   end
  end
 case 2
  if length(coilparam)<5
   disp('Insufficient number of components in coilparam: aborted')
  else
   if length(coilparam)<6,coilparam{6} = 0;end
  end
end
if isempty(coilparam)
  disp('No components left in coilparam: aborted')
  return
end

% Set up fbte18 observers for testing
%   First set consistent toroidal field sign for observers
MG = mgp(-1);
MG.if36fb = tor_sign;
MG.save(-1);
%  Now run fbte18
fbte18(-1,108)

mdsopen('pcs',-1)
n_e=8;n_f=8;n_oh=2;
coils = {'E1','E2','E3','E4','E5','E6','E7','E8', ...
         'F1','F2','F3','F4','F5','F6','F7','F8', ...
         'OH1','OH2','TOR','FPS'};
nps = length(coils)-1;
tdef = [0,2]'; sdef = [0,0]';% default traces
drawstr='\DRAW_FEEDFOR_';
[wpstr,sig,tsig,alimsel] = deal(cell(1,nps));
for k=1:nps
  if k <= n_e,              wpstr{k} = [drawstr,'E_I:ALIM_',iii(k)];
  elseif k <=n_e+n_f,       wpstr{k} = [drawstr,'F_I:ALIM_',iii(k-n_e)];
  elseif k <=n_e+n_f+n_oh,  wpstr{k} = [drawstr,'OH_I:ALIM_',iii(k-n_e-n_f)];
  else,                     wpstr{k} = [drawstr,'TOR_I:ALIM_001'];
  end
  sig{k} = sdef';tsig{k} = tdef';
  alimsel{k} = 'OFF';
end

t0 = [-20,-15,-12,-10]*1e-3;
sig0 = [0,40,0,0];
amplim = [repmat(7,1,16),repmat(26,1,2),72,2]*1e3;
sloplim = [repmat(100,1,16),repmat(150,1,2),70,1000]*1e3;
tlim = 4.5;
t_start_pol = Inf;
t_start_bphi = Inf;
t_stop_pol = -Inf;
t_stop_bphi = -Inf;
switch mode
  case 1 % one coil; requested amplitude reached in a specified number of steps with a specified slope within a specified total amount of time
   for icp=1:length(coilparam)
    coil = deblank(upper(coilparam{icp}{1}));
    amp = coilparam{icp}{2};
    slop = coilparam{icp}{3};
    nper = coilparam{icp}{4};
    t_tot = coilparam{icp}{5};
    tstart = coilparam{icp}{6};
    coilind = find(strcmpi(coil,coils));
    amp = max(min(amp,amplim(coilind)),-amplim(coilind));
    slop = max(min(slop,sloplim(coilind)),-sloplim(coilind));
    t_tot = max(min(t_tot,tlim),0);
    t_slop = 2*abs(amp)/slop;
    t_rise = t_slop/2/nper;
    t_flat = (t_tot-t_slop)/nper;
    if (t_slop>t_tot)
      disp('Slope is insufficient for sequence requested. Aborted')
      return
    end
    t = t_rise+[0,t_flat];
    t = repmat(t',1,nper)+repmat((0:nper-1)*(t_flat+t_rise),2,1);
    t = tstart+[t0,[0,t(:)',t_tot]];
    if strcmp(coil,'TOR')
      t_start_bphi=min(t_start_bphi,tstart-0.00005);
      t_stop_bphi = max(t_stop_bphi,tstart+t_tot);
    else
      t_start_pol=min(t_start_pol,tstart-0.00005);
      t_stop_pol = max(t_stop_pol,tstart+t_tot);
    end
    sig{coilind} = repmat(amp/nper*(1:nper),2,1);
    sig{coilind} = [sig0,0,sig{coilind}(:)',0];
    tsig{coilind} = t;
    alimsel{coilind} = 'ON';
   end
  case 2 % multiple coils in same group (TOR-FPS excluded) energized sequentially; requested amplitude reached with maximum allowed slope within a specified total amount of time,
%            with three possible waveforms defined by varargin{4}: 1=up only (square), 2=square up and spike down, 3=square up and down
    coil = deblank(upper(coilparam{1}));
    ncoil = coilparam{2};
    waveform = coilparam{3};
    amp = coilparam{4};
    t_tot = coilparam{5};
    tstart = coilparam{6};
    coilind = find(strcmpi(coil,coils));
    if coilind>n_e+n_f+n_oh
      disp('Coil name is outside allowed set')
      return
    end
    if coilind<=n_e
      ncoil=min(ncoil,n_e-coilind+1);
    elseif coilind<=n_e+n_f
      ncoil=min(ncoil,n_e+n_f-coilind+1);
    else
      ncoil=min(ncoil,n_e+n_f+n_oh-coilind+1);
    end
    amp = max(min(amp,amplim(coilind)),-amplim(coilind));
    slop = sloplim(coilind);
    t_tot = max(min(t_tot,tlim),0);
    t_slop = abs(amp)/slop;
    t_slopes = t_slop*2*ncoil;
    if waveform>1,t_slopes = 2*t_slopes;end
    t_flat = (t_tot-t_slopes)/ncoil;
    switch waveform
      case 1
        t = [-0.1,0,t_slop,t_slop+t_flat,2*t_slop+t_flat,2*t_slop+t_flat+0.1];
        sig1 = [0,0,amp,amp,0,0];
        dt = 2*t_slop+t_flat;
      case 2
        t = [-0.1,0,t_slop,t_slop+t_flat,3*t_slop+t_flat,4*t_slop+t_flat, ...
             4*t_slop+t_flat+0.1];
        sig1 = [0,0,amp,amp,-amp,0,0];
        dt = 4*t_slop+t_flat;
      case 3
        t = [-0.1,0,t_slop,t_slop+t_flat/2,3*t_slop+t_flat/2, ...
             3*t_slop+t_flat,4*t_slop+t_flat,4*t_slop+t_flat+0.1];
        sig1 = [0,0,amp,amp,-amp,-amp,0,0];
        dt = 4*t_slop+t_flat;
    end
    for k=1:ncoil
      sig{coilind+k-1} = [sig0,sig1];
      tsig{coilind+k-1} = tstart+[t0,t+dt*(k-1)];
      alimsel{coilind+k-1} = 'ON';
    end
    t_start_pol=tstart-0.00005;
    t_stop_pol=tstart+ttot;
end
[ya(:,1:nps),tc]=pcsalignipol(sig,tsig);
for k=1:nps
  wave_put(wpstr{k},ya(:,k),tc)
end

tstart = min(t_start_pol,t_start_bphi);
tstop = max(t_stop_pol,t_stop_bphi);
if isinf(t_start_bphi),t_start_bphi = t_start_pol;end
if isinf(t_start_pol),t_start_pol = t_start_bphi;end
if isinf(t_stop_bphi),t_stop_bphi = t_stop_pol;end
if isinf(t_stop_pol),t_stop_pol = t_stop_bphi;end
% timing references
vdbw(strcat('tcv_protectdb::shot_design:', ...
      {'ip','pol','bphi','ecrh','ecrh_a','ecrh_b','ecrh_c','flat','ramp'}, ...
        '_start_ref'),[t_start_pol,t_start_pol,t_start_bphi,0,0,0,0,0,0])
vdbw(strcat('tcv_protectdb::shot_design:', ...
      {'ip','pol','bphi','ecrh','ecrh_a','ecrh_b','ecrh_c','flat','ramp'}, ...
        '_stop_ref'),[t_stop_pol,t_stop_pol,t_stop_bphi,0,0,0,0,0,0]')
wave_put('\DRAW_REFS:REF_013',sdef,tdef)
wave_put('\DRAW_FEEDFOR_GAS:ALIM_001',sdef,tdef)
% the next line appears unnecessary since we operate in current mode (not in
%   hybrid mode)
pcshybcurrtovolt({'E','F','OH',;n_e,n_f,n_oh},ya(:,1:n_e+n_f+n_oh),tc,1.0)
% toroidal field sign
mdsw('\TOR_SIGN',tor_sign)
% Select necessary power supplies
nps = nps+1;
alimsel{nps} = 'ON';%force FPS
pcssetalim(repmat({'A'},1,nps),repmat({'100'},1,nps),alimsel)
% FPS default
if strcmpi(coil,'FPS')
  disp('Using default trace for FPS')
end
tfps = [tstart-0.005,tstart,tstop,tstop+0.005];
wave_put('\DRAW_FPS:UDCREF_001',250*ones(size(tfps)),tfps)
wave_put('\DRAW_FPS:UREF_001',zeros(size(tfps)),tfps)
wave_put('\DRAW_FPS:IREF_001',zeros(size(tfps)),tfps)
% obsolete? non-existent?
vdbw('tcv_protectdb::fps:modea','ON')
% integrator clamp/unclamp
vdbw(strcat('timers_operdb::timer_01',{'9','A','B','C'},':absolute'), ...
     [0.05,1.0,3.0,4.0])
% set up matrix addresses, no switching
tsw=[-0.05,0.3,1.9];
wave_put('\PHYS_MAT_ADDRESSES:A',zeros(size(tsw)),tsw)
wave_put('\PHYS_MAT_ADDRESSES:G',zeros(size(tsw)),tsw)
wave_put('\PHYS_MAT_ADDRESSES:M',zeros(size(tsw)),tsw)
% zero G matrices
for kG = '123'
  mdsw(['\LOAD_MAT_G',kG],zeros(24,32,10))
  mdsw(['\PHYS_MAT_G',kG],zeros(22,24,10))
end
% Turn off all RHVPS and mirror motors
pcsecrhmode('OFF')

%%%%%***** function separator #####

% Check limits and take corrective actions
function pcslimits
global tcv_recipe_error
recipetype = lower(deblank(vdbr('tcv_protectdb::shot_create:command_1')));
mdsopen('HYBRID',-1)
if any(strcmp(mdsr('\ALIM_SELECTS[16:17]'),'ON')) && ...
   contains(recipetype,'fbte')
% disruption detector Ip limit = 0.67 V
  curforbhi = 10e3+0.67*mdsr('\MAT_A_OUT_SCALES[12]');
  curforblo = 10e3;
  curtimmax = 0.25;
% special for runaway experiments
  if mdsr('\PCS::MGAMS.DATA:ISTOP[0]')<0
    curtimmax = 2.25;
  end
  curderlim = 20e6;%A/s
  % dt=1e-3;
  dtramp=10e-3;
% Ip reference
  ip_ref = mdsr('_ipref=\DRAW_REFS:REF_013');
  tip_ref = mdsr('dim_of(_ipref)'); 
  tip_ref(1)=min(tip_ref(1),-0.31);
  ip_sign = 2*(abs(max(ip_ref))>abs(min(ip_ref)))-1;
  if (abs(ip_ref(1))>curforblo || abs(ip_ref(end))>curforblo)
    disp(['FATAL ERROR: Plasma current at start or end exceeds ', ...
          num2str(curforblo,'%5.2e'),' A'])
    tcv_recipe_error=2;
    return
  end
% M inputs for OH coils
  [~,oh_ind] = ...
        max(abs(mdsr('\LOAD_MAT_M[$1,*]',find(contains( ...
         upper(mdsr('\HYBRID::MAT_M:OUT_CABLING')),'OH'))-1)),[],2);
% G2 elements from Ip observer to OH coil voltages
  g2mat = mdsr('\LOAD_MAT_G2[$1,12]',str2double(cellfun(@(x) x(4:5), ...
            mdsr('\HYBRID::MAT_M:IN_CABLING[$1]',oh_ind-1), ...
                'UniformOutput',false))-1,12);
  g2addr = mdsr('_add=\PHYS_MAT_ADDRESSES:G')+1;
  g2times = mdsr('dim_of(_add)')';
  ip_ref = [ip_ref;interp1(tip_ref,ip_ref,g2times','linear',0)];% add matrix switching times to Ip vector
  tip_ref = [tip_ref;g2times'];
  [ip_ref,tip_ref] = pcsalignipol({ip_ref},{tip_ref});
  tip_ref = [tip_ref;Inf];% this time vector has an extra point
  ip_ref = [ip_ref;0];
  g2times = [g2times,Inf];% extra time at +Inf
  nfi = find(~sum(abs(g2mat(:,:,g2addr))))';% indices of no-feedback matrices
  if isempty(nfi)
   t_int = [-Inf,g2times(end-1);g2times(1),Inf];
  else
   inc = find(diff(nfi)>1);% indices of non-consecutive no-feedback matrices
   t_int = [g2times(nfi(1)),g2times(nfi(inc+1)); ...
            g2times(nfi(inc)+1),g2times(nfi(end)+1)];% no-feedback time ranges
  end
  ipi = 0;t_extra = [];ip_extra = [];remi = [];
  for k = 1:length(t_int(:))-1 % alternate no-feedback and feedback ranges
   ipi = (ipi(end)+1):find(tip_ref<=t_int(k+1),1,'last');% relevant Ip indices
   if ~rem(k,2) %derivative check only for Ip feedback ranges
    if any(abs(diff(ip_ref(ipi))./diff(tip_ref(ipi)))>curderlim)
     disp(['FATAL ERROR: Plasma current derivative exceeds ', ...
         num2str(curderlim,'%5.2e'),' A/s'])
     tcv_recipe_error=2;
     return
    end
   end
   forbi = find(ip_ref(ipi)*ip_sign>=curforblo&ip_ref(ipi)*ip_sign<curforbhi)';% points in forbidden zone within Ip subvector
   if ~isempty(forbi)
    iinc = find(diff(forbi)>1);% indices of non-consecutive forbidden-zone points
    inti = [ipi(forbi(1)),ipi(forbi(iinc+1));ipi(forbi(iinc)),ipi(forbi(end))];% alternate forbidden and allowed index ranges
    inti1 = [inti(1,:)-1;inti(2,:)+1];% indices adjacent to inti on the outside
    t_iint = tip_ref(inti);
    t_iint1 = tip_ref(inti1);
    ip_iint = ip_ref(inti)*ip_sign;
    ip_iint1 = ip_ref(inti1)*ip_sign;
    curlim = curforblo*(ip_iint>ip_iint1)+curforbhi*(ip_iint1>ip_iint);% relevant current limit for each crossing
    tcr = (t_iint.*(ip_iint1-curlim)-t_iint1.*(ip_iint-curlim))./ ...
          (ip_iint1-ip_iint);% boundary crossing times
    tcr(1,:) = max(tcr(1,:),t_int(k));
    tcr(2,:) = min(tcr(2,:),t_int(k+1));
    if any(diff(tcr)>curtimmax)
     if rem(k,2) % no-feedback case: clamp
      for kj=find(diff(tcr)>curtimmax)
       t_extra = [t_extra,tcr(1,kj),tcr(1,kj)+dtramp, ...
                  tcr(2,kj)-dtramp,tcr(2,kj)];% add an extra trapeze to clamp current to maximum
       ip_extra = [ip_extra,curlim(1,kj),curforbhi,curforbhi,curlim(2,kj)];
       remi = [remi,inti(1,kj):inti(2,kj)];% all points inside interval are scheduled for removal
      end
     else% Ip-feedback case: fatal error
      disp(['FATAL ERROR: Plasma current in forbidden zone (', ...
            num2str(curforblo,'%5.2e'),' to ',num2str(curforbhi,'%5.2e'), ...
            ' A) for over ',num2str(curtimmax,'%6.3e'), ...
            ' s with Ip feedback on'])
      tcv_recipe_error=2;
      return
     end
    end
   end
  end
  if ~isempty(t_extra)
   tip_ref([remi,end]) = [];% remove scheduled points + additional end point
   ip_ref([remi,end]) = [];
   ip_ref = [ip_ref;ip_extra(:)];
   [tip_ref,k] = sort([tip_ref;t_extra(:)]);
   [ya,tc]=pcsalignipol({ip_ref(k)},{tip_ref});
   wave_put('\DRAW_REFS:REF_013',ya,tc)
  end
end

%%%%%***** function separator #####

function pcssetup
global tcv_recipe_error
recipetype = lower(deblank(vdbr('tcv_protectdb::shot_create:command_1')));

% Book-keeping of A matrix outputs (feedback channels)
Aout = upper(mdsdata('\PHYS_MAT_A_OUTPUTS'));
loadi = str2double(cellfun(@(x){x(3:5)},Aout));% mapping of physics to loadable signals
chty = cellfun(@(x){x(7:end)},Aout(loadi));% channel type
vdbw(strcat('tcv_protectdb::fback:chan_',cellstr(iii(1:size(Aout,1)))),chty)

% Book-keeping of PCS switches (on actuator side)
mout = upper(mdsdata('\PHYS_MAT_M_OUTPUTS'));
nmout = numel(mout); % Do not assume size is 20
nsw = nmout+12+16;
swn = cell(nsw,1);
swn(1:nmout) = mout;% 1-20: M matrix outputs
for k=5:16
  swn{nmout+k-4} = upper(mdsdata([ ...
      'decompile(`getnci("\\load_wave_gen_a:channel_',iii(k),'","record"))']));% NMOUT+(1:12): Wavegen A 5-16
end
for k=1:16
  swn{nmout+12+k} = upper(mdsdata('\summator.summer_001.output_names[$1]',k-1));% NMOUT+12+(1:16): summator 1
end
swn = deblank(swn);
% Coil Power Supply Voltages
ki = startsWith(swn,{'E_','F_','OH_','FAST'});
swn(ki) = strcat('U_',swn(ki));% precede by U_ to signify voltage
% GAS
ki = strncmp(swn,'GAS',3)&cellfun(@length,swn)==3;
swn(ki) = {'GAS1'};% uniform notation
% Remove conditional part
swn = regexprep(swn,'[^?]+\? ','','once');
swn = regexprep(swn,' : \\PCS::DRAW_ZERO','');
% Unused channels
ki = contains(swn,'DRAW_ZERO');
swn(ki) = {'unused'};
% keep everything after DRAW_FEEDFOR_ or DRAW_FPS: or DRAW_THRESH:
swn = regexprep(swn,'.*DRAW_(FEEDFOR_|FPS:|THRESH:)','');
% keep everything after DRAW_ if no subnode
swn = regexprep(swn,'.*DRAW_([^:]*).*','$1');
% chuck everything after * (multiplication sign) or / (division sign)
swn = regexprep(swn,'[\*\/].*','');

% put names into Vsystem
vdbw(strcat('tcv_protectdb::name:chan_',cellstr(iii(1:nsw))),swn)

% close or open switches
gates = repmat({'ON'},nsw,1);
if ~contains(recipetype,{'timer_only','diagonly'})
  gates(~strcmp(swn,'unused')) = {'OFF'};
end
vdbw(strcat('tcv_protectdb::gates:chan_',cellstr(iii(1:nsw))),gates)

% reset summators except for *fbte* recipe
if ~contains(recipetype,'fbte')
  mdsdata('JP_SET_SUM_SYS($)','hybrid');
  mdsw('\LOAD_WAVE_MASK',int8(ones(1,128)))
end

% programmable amplifier gains
pag = mdsdata('\LOAD_FEED_GAINS');
if (length(pag)==1),pag=repmat(pag,size(Aout,1),1);end
if any(~ismember(fix(pag),[1 2 4 8]))
  disp('FATAL ERROR: illegal programmable amplifier gain values')
  tcv_recipe_error=2;
  return
end
vdbw(strcat('tcv_protectdb::gain:chan_',cellstr(iii(1:size(Aout,1)))), ...
        log2(pag))

% types of references
ngr = 4;
nty = 6;
refty = mdsdata('\LOAD_WAVE_GEN_B:REFS');
ipi = find(strcmp(refty,'Ip'));
reftyg = repmat({'ON'},nty,1);
if ~contains(recipetype,'pcs_test')% force to 10 V for PCS test
  reftyg(unique(ceil(ipi/ngr)))={'OFF'};% Ip by request, in groups of four (if at least one in the group is requested)
end
vdbw(strcat('tcv_protectdb::wavegen_b:source_', ...
       cellstr(iii(((1:nty)-1)*ngr+1)),'_',cellstr(iii((1:nty)*ngr))),reftyg)

pcsalimsetup

%%%%%***** function separator #####

% Set up the power supplies
function pcsalimsetup

global tcv_recipe_error
recipetype = lower(deblank(vdbr('tcv_protectdb::shot_create:command_1')));

n_e = 8; n_f = 8; n_oh = 2;
alim = [strcat({'e'},num2str(1:n_e,'%i')');strcat({'f'},num2str(1:n_f,'%i')'); ...
        strcat({'oh'},num2str(1:n_oh,'%i')');'tor';'fps'];

mode = mdsdata('\ALIM_MODES');
lev  = str2double(mdsdata('\ALIM_LEVELS'));
sel  = mdsdata('\ALIM_SELECTS');
sbit = mdsdata('\LOAD_SIGN_BITS');

vdbw(strcat('tcv_protectdb::',alim,':select'),sel)

% FPS
lev_fps = fix(lev(20))*2.3;
vdbw('prot_232db::fps_1:limit_set',lev_fps)
vdbw(['prot_232db::fps_1:mode',mode{20}],'ON')
%TODO: replce this with an on_off (PAL?)
pause(1)
vdbw(['prot_232db::fps_1:mode',mode{20}],'OFF')

% E, F, OH
efoh_ctrl = ones(1,18)*2;
efoh_ctrl(16) = 1; %1=OPCUA (Imperix), 2=traditional
efoh_db = {'opcua_rect','wago_jsi'};
efoh_mod = {'cmd_mode','io_mod_crt'};
efoh_lim = {'idc_lim','io_lim_crt'};
i_ctrl = find(efoh_ctrl==1)';
efoh_mods(i_ctrl) = strcat(efoh_db(1),'_',alim(i_ctrl), ...
			   'db::',efoh_mod(1));
efoh_modr(i_ctrl) = strcat(efoh_db(1),'_',alim(i_ctrl), ...
			   'db::',efoh_mod(1),'_R');
efoh_lims(i_ctrl) = strcat(efoh_db(1),'_',alim(i_ctrl), ...
			   'db::',efoh_lim(1));
efoh_limr(i_ctrl) = strcat(efoh_db(1),'_',alim(i_ctrl), ...
			   'db::',efoh_lim(1),'_R');
lev(i_ctrl) = lev(i_ctrl)/100*7700;
i_ctrl = find(efoh_ctrl==2)';
efoh_mods(i_ctrl) = strcat(efoh_db(2)','db::',alim(i_ctrl), ...
			   ':',efoh_mod(2)');
efoh_modr(i_ctrl) = strcat(efoh_db(2)','db::',alim(i_ctrl), ...
			   ':o_r_mod_crt');
efoh_lims(i_ctrl) = strcat(efoh_db(2)','db::',alim(i_ctrl), ...
			   ':',efoh_lim(2)');
efoh_limr(i_ctrl) = strcat(efoh_db(2)','db::',alim(i_ctrl), ...
			   ':o_r_lim_crt');
lev(i_ctrl) = lev(i_ctrl)*10;

rsel = find(strcmp(deblank(sel(1:n_e+n_f+n_oh)),'ON'));
act = rsel;
actsel{1} = intersect(act,find(efoh_ctrl==1));
actsel{2} = intersect(act,find(efoh_ctrl==2));
for j=1:2% two attempts
  for ity=1:2
    if ~isempty(actsel{ity})
      moden{ity} = double(char(mode(actsel{ity}))-'A');
      vdbw(efoh_mods(actsel{ity}),moden{ity})
      vdbw(efoh_lims(actsel{ity}),lev(actsel{ity}))
    end
  end
  pause(3)
  for ity=1:2
    if ~isempty(actsel{ity})
      moder = vdbr(efoh_modr(actsel{ity}));
      levr = vdbr(efoh_limr(actsel{ity}));
      actsel{ity} = actsel{ity}(moder~=moden{ity}| ...
        abs(levr-fix(lev(actsel{ity})))>0.2*fix(lev(actsel{ity})));
    end
  end
  if isempty([actsel{1};actsel{2}]),break;end
  if j==2
    act = union(actsel{1},actsel{2});
    fprintf('FATAL ERROR: wrong levels or modes on some power supplies:\n');
    fprintf('\talim=%s, requested mode=%s, requested level=%d\n',alim{act},mode{act},lev(act));
    tcv_recipe_error=2;
    return
  end
  pause(10)
end

% need motor generator and power supplies if any coil is requested or ECRH (not
%   on grid) or NBH; only motor generator for FPS
coilson = any(strcmp(deblank(sel),'ON'));
if contains(recipetype,{'ecrhonly','nbhonly'}) || coilson
  vdbw(strcat({'tcv_protectdb::alims','prot_232db::mt'},':select'),{'ON';'ON'})
else
  vdbw(strcat({'tcv_protectdb::alims','prot_232db::mt'},':select'),{'OFF';'OFF'})
end
if contains(recipetype,'fpsonly')
  vdbw('prot_232db::mt:select','ON')
end

% ensure last finite-voltage FPS point is after JSI end gate
tstop = vdbr('tcv_protectdb::shot_design:pol_stop_ref');
fpsv = mdsdata('_f=\DRAW_FPS:UDCREF_001');
tfpsv = mdsdata('dim_of(_f)');
if max(tfpsv)<tstop
  wave_put('\DRAW_FPS:UDCREF_001',[fpsv;fpsv(end)],[tfpsv;tstop+0.05])
end

% toroidal field sign
if strcmpi(deblank(vdbr('tcv_protectdb::tor:select')),'ON')
  n_to = 10;
  condch = 'wago_jsidb::inv:o_etat';
  tor_sign = mdsr('\TOR_SIGN');
  if vdbr(condch)~=1 && vdbr(condch)~=3
    disp('FATAL ERROR: toroidal switch position undefined')
    tcv_recipe_error=2;
    return
  end
  i_to = 0;
  while (tor_sign>0 && vdbr(condch)~=3) || (tor_sign<0 && vdbr(condch)~=1)
    i_to = i_to+1;
    if i_to>=n_to,return;end
    if tor_sign>0
      mdsdata('SWITCH_1_0($1,$2)','wago_jsidb','inv:io_cmd_inv');
    else
      mdsdata('SWITCH_1_0($1,$2)','wago_jsidb','inv:io_cmd_not_inv');
    end
    pause(1.5)
  end
end

% sign bits
sbitv = repmat({'OFF'},n_e+n_f+n_oh,1);
sbitv(sbit>0) = {'ON'};
vdbw(strcat('tcv_protectdb::',alim(1:n_e+n_f+n_oh),':signe'),sbitv)

%%%%%***** function separator #####

% Set up power supply gates
function pcsalimgates

recipetype = lower(deblank(vdbr('tcv_protectdb::shot_create:command_1')));

magncard = 32;
invaltime = -10;
curstart = vdbr(strcat('TCV_PROTECTDB::SHOT_DESIGN:',{'POL','BPHI'}, ...
                       '_START_REF'))-0.01;
curstop = vdbr(strcat('TCV_PROTECTDB::SHOT_DESIGN:',{'POL','BPHI'}, ...
                       '_STOP_REF'))+[0.01;0.11];
alims = {'E1','E2','E3','E4','E5','E6','E7','E8', ...
         'F1','F2','F3','F4','F5','F6','F7','F8','OH1','OH2', ...
         'FPS','FPS','TOR','TOR'};
% E coils, F coils, OH coils, FPSDC, FPSFOS, TOR, magnetics integrator, DML
timers_on = [31+2*(0:3),41+2*(0:3),51+2*(0:3),61+2*(0:3),71+2*(0:1), ...
             431,434,075,325,323];
timers_off = timers_on+1;
idrift = find(floor(timers_on/10)==magncard);%magnetics card has special requirements
time_on = invaltime*ones(1,length(timers_on));
time_off = invaltime*ones(1,length(timers_off));
ntime = length(time_on);
coilstart = [repmat(curstart(1),1,18), ...
             curstart(1)-[0.5,0],min(curstart),curstart(2)-0.25];
coilstop = [repmat(curstop(1),1,18), ...
             curstop(1)+[0.01,0],max(curstop),curstop(2)+0.25];
recipes = {'fbte','stray','ecrhonly','nbhonly','backoff', ...
        'fpsonly','gasonly','diagonly','pcs_test','timer_only'};
recind = find(contains(recipes,recipetype));
if ismember(recind,[1,2,5,6])
  alimsel = find(strcmpi(deblank(...
                vdbr(strcat('TCV_PROTECTDB::',alims,':SELECT'))),'ON'));
  time_on(alimsel) = coilstart(alimsel);
  time_off(alimsel) = coilstop(alimsel);
  if (any(time_on(idrift)==invaltime))
    time_on(idrift)=min(curstart)-0.25;
    time_off(idrift)=max(curstop)+0.85;
  end
end
alimnodes = strcat([repmat('TIMERS_OPERDB::TIMER_', ...
                length(timers_on)-2,1);repmat('TIMERS_DIAGDB::TIMER_',2,1); ...
                repmat('TIMERS_OPERDB::TIMER_', ...
               length(timers_off)-2,1);repmat('TIMERS_DIAGDB::TIMER_',2,1)], ...
                reshape(num2str([timers_on,timers_off],'%03d')',3, ...
                        length(timers_on)+length(timers_off))',':ABSOLUTE');
vdbw(alimnodes([1:ntime-2,ntime+1:2*ntime-2],:), ...
     invaltime*ones(2*ntime-2,1))
vdbw(alimnodes([1:ntime-2,ntime+1:2*ntime-2],:), ...
     fix([time_on(1:end-2),time_off(1:end-2)]'*500)*0.002)
vdbw(alimnodes([ntime-1:ntime,2*ntime-1:2*ntime],:),fix([time_on(end-1:end),time_off(end-1:end)]'*500)*0.002)
% Motor generator
if ismember(recind,[1,2,3,4,5,6])
  if recind==4
    curstart = [-1.6,curstart];
    curstop = [curstop,2];
  end
  timers_mt = 330+(1:6);
  time_mt = [min(curstart)-[5.7,4.7,5.2,4.2],max(curstop)+[0.5,1.5]];
  mtnodes = strcat('TIMERS_OPERDB::TIMER_', ...
                reshape(num2str(timers_mt,'%03d')',3, ...
                        length(timers_mt))',':ABSOLUTE');
  vdbw(mtnodes,fix(time_mt'*500)*.002)
% TODO: these writes don't actually work (protected?) - anything to do here?
  vdbw('PROT_232DB::MT:SET_SPEED',117.0)
  vdbw('PROT_232DB::MT:SET_IDW',537.8)
  vdbw('PROT_232DB::MT:SET_VOLTS',10.18)
end

%%%%%***** function separator #####

% Protect ECRH when not needed
function pcsecrhmode(mode)

if strcmpi(mode,'OFF')
% turn off launchers
  nlaun=6;
  vdbw(strcat('tcv_ecrh_x2db::x2_',cellstr(num2str(1:nlaun,'%1i')')', ...
              ':wr_etat'),cellstr(repmat('off',nlaun,1)))
  button_with_timeout('tcv_ecrh_x3db','gr_c:bt_0_cm4', ...
                      'tcv_ecrh_x3db::gr_c:bt_0_a1',0,10,1.5)
% turn off RHVPS
% Q: are these channels defined??
  vdbw(strcat('tcv_ecrh_x',{'2','2','3'},'db::rhv_',{'a','b','c'},':pre_on'), ...
       {'ON','ON','ON'}')
% Q: third one exists in both X2 and X3??
  vdbw(strcat('tcv_ecrh_x',{'2','2','3'},'db::gr_',{'a','b','c'},':status'), ...
       {'set to preset mode','set to preset mode','set to preset mode'}')
end

%%%%%***** function separator #####

% Set up ECRH
function pcsecrhsetup

global tcv_recipe_error
% ECHCS_AUX_PARAM_05 parameter handles exception allowing ECRH injection with
%   no plasma current
mdsopen('HYBRID',-1)
if mdsr('\ECHCS_AUX_PARAM_05')
  ip_thres = 0;
else
  ip_thres = 20000;% current threshold
end
tforb = 0.035;% maximum time in Ip forbidden zone
gymax = 11;laumax = 11;
cl = {'A','B','C'};% clusters
gy = {[1,2],[10,11],[7,8,9]};% gyrotrons
lau = {[1,2],[4,5],[],[10,11]};% launchers
ltog = mdsr('\LINETOGYRO');
gtol = zeros(1,gymax);
for ic=1:length(gy)
  [~,gtol(gy{ic})] = ismember(gy{ic},ltog);
end
gtol(~gtol)=NaN;
timstac = {'451','C11','631'};% cathode start timers
timstoc = {'452','C12','632'};% cathode stop timers
timbods = {'503','505','501','507'};% body power supply start timers
timbode = {'504','506','502','508'};% body power supply stop timers
timstai(gy{2}) = {'C04','C54'};% XU anode inhibit on timers
timstoi(gy{2}) = {'C05','C55'};% XU anode inhibit off timers
timstag(gy{2}) = {'C01','C51'};% XU anode start timers
timstog(gy{2}) = {'C02','C52'};% XU anode stop timers
timstag(gy{3}) = {'651','653','655'};% X3 anode start timers
timstog(gy{3}) = {'652','654','656'};% X3 anode stop timers
timstagc = '657';% prep gate for X3 anodes
timstogc = '658';
tactr = [0,0.005]; %empirical activation delay for RHVPS
tdeact = [0.008,0,0];% empirical deactivation delay for RHVPS
tgates = [-0.009,-0.05,-0.005];% starting gate time shift for cathodes
tgatee = [0.007,0.05,0.01];% stopping gate time shift for cathodes
tbodgates = [-0.008,-0.008,-0.038,-0.038];% starting gate time shift for body supplies
tbodgatee = [-0.001,0.007,0.032,0.032];% stopping gate time shift for body supplies
tact(gy{2}) = [.022,.022];% empirical activation delay for anodes
tgate(gy{2}) = [-.1,-.1];% gate time shift for anodes
tact(gy{3}) = [.0,.0,.0];% empirical activation delay for anodes
tgate(gy{3}) = [-.002,-.002,-.002];% gate time shift for anodes
tdef = [0,2]'; sdef = [0,0]';% default traces
timcdef = -5;% default timer
% nt_to = 20;% time-out limit
nsumm = 3;% summator branches
nscat = 1;nsan = 14;nslx2 = 4;%first index of summator channels for cathodes, anodes, X2 mirrors
terr = 1;%time interval used to calculate minimum significant slope from signal accuracy
dbx2 = 'tcv_ecrh_x2db';
dbx3 = 'tcv_ecrh_x3db';
dbxu = 'tcv_ecrh_xudb';
dbrx2t = 'timers_x2db::timer_';
dbrx3t = 'timers_x3db::timer_';
dbrxut = 'timers_xudb::timer_';
suft = ':absolute';
densloroot = '\DRAW_THRESH:EC_DENS_LOW';
denshiroot = '\DRAW_THRESH:EC_DENS_HIGH';
if strcmpi(deblank(vdbr([dbxu,'::launcher_10:select'])),'ON') || ...
   strcmpi(deblank(vdbr([dbxu,'::launcher_11:select'])),'ON')
     vdbw([dbxu,'::ecrh_acq:select'],'ON');
end
x2req = strcmpi(deblank(vdbr([dbx2,'::ecrh_acq:select'])),'ON');
xureq = strcmpi(deblank(vdbr([dbxu,'::ecrh_acq:select'])),'ON');
x3req = strcmpi(deblank(vdbr([dbx3,'::ecrh_acq:select'])),'ON');
greq (1:gymax) = false;
ecintcv (1:gymax) = false;
dbrhm = 'tcv_protectdb::h_mere';
hm_t0 = vdbr([dbrhm,':t0']);
hm_period_f = vdbr([dbrhm,':period_f']);
hm_n_slow = vdbr([dbrhm,':n_slow']);
hm_period_s = hm_period_f*hm_n_slow*1.e-6 ;
if x2req || x3req || xureq
% master clock parameters
% get planned current
  ip = mdsdata('_x=tcv_eq("i_p","fbte")');
  tip = mdsdata('dim_of(_x)');
% add two extra points
  ip = [0;ip;0];
  tip = [-Inf;tip;Inf];
% get periods with current below threshold
  forbi = find(abs(ip)<ip_thres)';% points in forbidden zone
  inc = find(diff(forbi)>1);% indices of non-consecutive forbidden-zone points
  inti = [forbi(1),forbi(inc+1);forbi(inc),forbi(end)];% alternate forbidden and allowed index ranges
  inti1 = max(min([inti(1,:)-1;inti(2,:)+1],length(ip)),1);% indices adjacent to inti on the outside
  t_int = tip(inti);
  t_int1 = tip(inti1);
  ip_int = ip(inti);
  ip_int1 = ip(inti1);
  tcr = (t_int.*(ip_int1-ip_thres)-t_int1.*(ip_int-ip_thres))./ ...
        (ip_int1-ip_int);% boundary crossing times
  if isnan(tcr(1)),tcr(1)=-Inf;end
  if isnan(tcr(end)),tcr(end)=Inf;end
% summator settings
  mdsdata('_summn=\SUMMATOR:INPUT_NAMES');
  summ = mdsdata('\SUMMATOR.SUMMER_001.INPUTS');
end

% X3 power
avsl(7) = mdsdata(['_avsl=0.5/USING(DATA(\X3_LAU_SETUP:ANGLE_SLOPE),', ...
                '\TOP,-1,"ECRH")']);
avof(7) = mdsdata(['_avof=2.*_avsl*(USING(DATA(\X3_LAU_SETUP:ANGLE_ZERO),', ...
                 '\TOP,-1,"ECRH")-90.)']);
ic = 3;%cluster number
[va,tva] = deal(cell(1,gymax));
[tsa,tea] = deal(zeros(1,gymax));
if x3req
% X3 requested into TCV?
% at this stage, greq and ecintcv do not incorporate the information that the
%   power traces are nonzero
  for k=gy{3}
   kl = gtol(k);
   greq(k) = strcmpi(mdsdata(['\online_gyro_',num2str(k,'%02d')]),'ON');
   ecintcv(k) = ~isnan(kl) && (greq(k) && ...
                strcmpi(mdsdata(['\switch_line_',num2str(kl,'%02d')]),'TCV'));
  end
% power traces
  vc{3} = mdsdata('_vc=min(\RHVPS_VOLT_C,maxval(\LIM_RHVPS_C:VMAX))');
  vc{3}(vc{3}~=0) = max(vc{3}(vc{3}~=0),mdsdata('maxval(\LIM_RHVPS_C:VMIN)'));
  tvc{3} = mdsdata('dim_of(_vc)');
  if length(vc{3})<2,tvc{3} = tdef;vc{3} = sdef;end
  [vc{3},tvc{3},~,~] = massage_traces(vc{3},tvc{3},[],[],0.003, ...
        tdeact(3),[],mdsdata('\LIM_RHVPS_C:DV_DT_ERR'), ...
        mdsdata('\LIM_RHVPS_C:DV_DT_MIN'),mdsdata('\LIM_RHVPS_C:DV_DT_MAX'), ...
        mdsdata('minval(\LIM_RHVPS_C:VMIN)'), ...
                mdsdata('maxval(\LIM_RHVPS_C:VMAX)'));
  for k=gy{3}
   ks = num2str(k,'%02d');
   va{k} = mdsdata(['_va=min(\ANODE_VOLT_GYRO_',ks, ...
                          ',\LIM_ANODE_',ks,':VMAX)']);
   va{k}(va{k}~=0) = max(va{k}(va{k}~=0),mdsdata(['\LIM_ANODE_',ks,':VMIN']));
   tva{k} = mdsdata('dim_of(_va)');
   if length(va{k})<2,tva{k} = tdef;va{k} = sdef;end
   [va{k},tva{k},tsa(k),tea(k)] = massage_traces(va{k},tva{k},[],[],tact(k), ...
                             [],[],mdsdata(['\LIM_ANODE_',ks,':DV_DT_ERR']), ...
                                   mdsdata(['\LIM_ANODE_',ks,':DV_DT_MIN']), ...
                                   mdsdata(['\LIM_ANODE_',ks,':DV_DT_MAX']), ...
                                   mdsdata(['\LIM_ANODE_',ks,':VMIN']), ...
                                   mdsdata(['\LIM_ANODE_',ks,':VMAX']));
   greq(k) = greq(k)&any(va{k});
  end
% timers
  if any(greq(gy{3}))
% greq and ecintcv now incorporate the information that the anode trace is
%   nonzero
   for i=gy{3}
     ecintcv(i) = ecintcv(i)&any(va{k});
   end
   greqx3 = greq(gy{3});
   ts(ic) = min(tsa(gy{3}(greqx3)));
   te(ic) = max(tea(gy{3}(greqx3)));
   vdbw(strcat(dbrx3t,timstag(gy{3}(greqx3)),suft), ...
        tsa(gy{3}(greqx3))+tgate(gy{3}(greqx3)))% anode gates
   vdbw(strcat(dbrx3t,timstog(gy{3}(greqx3)),suft),tea(gy{3}(greqx3))+0.003)
   vdbw(strcat(dbrx3t,{timstagc,timstogc},suft),[ts(ic)-1,te(ic)+0.2])% anodes' prep gates
   vcmod = mdsdata('\RHVPS_MOD_C:MODULATOR[0]');
   if ~strcmpi(vcmod(1),'INTERNAL')
    if ~any(vc{3}) || tvc{3}(1)+.002>ts(ic)
     disp('FATAL ERROR: RHVPS C trace starts too late for anodes')
     tcv_recipe_error=2;
     return
    end
    if tvc{3}(end)-.002<te(ic)
     disp('FATAL ERROR: RHVPS C trace stops too early for anodes')
     tcv_recipe_error=2;
     return
    end
    ts(ic) = min(ts(ic),tvc{3}(2)-.002);
    te(ic) = max(te(ic),tvc{3}(end-1)+.002);
    vdbw(strcat(dbrx3t,{timstac{ic},timstoc{ic}},suft), ...
         [ts(ic)+tgates(ic),te(ic)+tgatee(ic)])% cathode gates
% pulse length check vs limit
    if te(ic)-ts(ic)>mdsdata('\LIM_RHVPS_C:TMAX')-0.09
      disp('FATAL ERROR: cluster C power traces exceed the maximum allowed duration')
      tcv_recipe_error=2;
      return
    end
   end
  else
    ts(ic) = 0;
    te(ic) = 0;
    vdbw(strcat(dbrx3t,{timstagc,timstogc,timstag{gy{3}},timstog{gy{3}}, ...
                        timstac{ic},timstoc{ic}}),timcdef)
  end
% pulse length check vs Ip
  if any(ecintcv(gy{3}))
   is = find(tcr(1,:)<=ts(ic),1,'last');
   ie = find(tcr(2,:)>=te(ic),1,'first');
   if min(tcr(2,is),te(ic))-ts(ic)+te(ic)-max(tcr(1,ie),ts(ic))+ ...
      sum(diff(tcr(:,is+1:ie-1)))>tforb
     disp('FATAL ERROR: finite cluster C power requested in TCV for too long at too low a plasma current')
     tcv_recipe_error=2;
     return
   end
  end
% summator settings
  mdsdata('_mod=\RHVPS_MOD_C:MODULATOR');
  for kk=1:nsumm
   summ(nscat-1+ic,kk) = mdsdata(['(_mod[0]=="PCS"&&', ...
                                 'any(_mod[1:3]==_summn[$1]))'],kk-1);
  end
  for k=gy{3}
   mdsdata(['_mod=\ANODE_MOD_GYRO_',num2str(k,'%02d'),':MODULATOR']);
   for kk=1:nsumm
    summ(nsan-k(1)+k,kk) = mdsdata(['(_mod[0]=="PCS"&&', ...
                                   'any(_mod[1:3]==_summn[$1]))'],kk-1);
   end
  end
else
  tvc{3} = tdef;vc{3} = sdef;
  for k=gy{3}
    tva{k} = tdef;va{k} = sdef;
  end
  tal{7} = tdef;al{7} = (sdef-avof(7))/avsl(7);
  ts(ic) = 0;te(ic) = 0;
  ecintcv(gy{3}) = false;
end
% write traces
wave_put('\DRAW_FEEDFOR_ECRH:ALIM_C',vc{3},tvc{3})
for k=gy{3}
  wave_put(['\DRAW_FEEDFOR_ECRH:X3_ANODE_',iii(k)],-va{k},tva{k});%with sign inversion
end

% reference times
vdbw(strcat('tcv_protectdb::shot_design:ecrh_c_',{'start','stop'},'_ref'), ...
     [ts(ic),te(ic)])

% X2 power
ic = 1;%cluster number
if x2req
% X2 requested into TCV?
% at this stage, greq and ecintcv do not incorporate the  information that the
%   RHVPS trace is nonzero
   for il = gy{ic}
    kl = gtol(il);
    greq(il) = strcmpi(mdsdata(['\online_gyro_',num2str(il,'%02d')]),'ON');
    ecintcv(il) = ~isnan(kl)&&(greq(il) && ...
                  strcmpi(mdsdata(['\switch_line_',num2str(kl,'%02d')]),'TCV'));
   end
% power traces
   power = mdsdata(['\POWER_CLUSTER_',cl{ic}]);
   vc{ic} = mdsdata(['_vc=max(min(\RHVPS_VOLT_',cl{ic}, ...
             ',maxval(\LIM_RHVPS_', ...
             cl{ic},':VMAX)),minval(\LIM_RHVPS_',cl{ic},':VMIN))']);
   vc{ic}(power<=0)=0;
   tvc{ic} = mdsdata('dim_of(_vc)');
   if length(vc{ic})<2,tvc{ic} = tdef;vc{ic} = sdef;end
   [vc{ic},tvc{ic},ts(ic),te(ic)] = massage_traces(vc{ic},tvc{ic},[],[], ...
        tact(ic),tdeact(ic),[],mdsdata(['\LIM_RHVPS_',cl{ic},':DV_DT_ERR']), ...
                mdsdata(['\LIM_RHVPS_',cl{ic},':DV_DT_MIN']), ...
                mdsdata(['\LIM_RHVPS_',cl{ic},':DV_DT_MAX']), ...
                mdsdata(['minval(\LIM_RHVPS_',cl{ic},':VMIN)']), ...
                mdsdata(['maxval(\LIM_RHVPS_',cl{ic},':VMAX)']));
   if ~strcmpi(mdsdata(['\RHVPS_MOD_',cl{ic},':MODULATOR[0]']),'INTERNAL')
     ts(ic) = tvc{ic}(1);
   end
% now greq and ecintcv incorporate the information that the power trace is
%   nonzero
   greq(gy{ic}) = greq(gy{ic})&any(vc{ic});
% timers
   if any(greq(gy{ic}))
    ecintcv(gy{ic}) = ecintcv(gy{ic})&any(vc{ic});
    % greqx2 = greq(gy{ic});
    vdbw(strcat(dbrx2t,{timstac{ic},timstoc{ic}},suft), ...
         [ts(ic)+tgates(ic),te(ic)+tgatee(ic)])% cathode gates
    vdbw(strcat(dbrx2t,[timbods,timbode],suft), ...
         [ts(ic)+tbodgates,te(ic)+tbodgatee])% gates for body power supplies
 % pulse length check vs limit
    if te(ic)-ts(ic)>mdsdata(['\LIM_RHVPS_',cl{ic},':TMAX'])-0.09
      disp(['FATAL ERROR: cluster ',cl{ic}, ...
            ' power traces exceed the maximum allowed duration'])
      tcv_recipe_error=2;
      return
    end
   else
    ts(ic) = 0;
    te(ic) = 0;
    vdbw(strcat(dbrx2t,[timstac(ic),timstoc(ic), ...
                        timbods,timbode]),timcdef)
   end
% pulse length check vs Ip
   if any(ecintcv(gy{ic}))
    is = find(tcr(1,:)<=ts(ic),1,'last');
    ie = find(tcr(2,:)>=te(ic),1,'first');
    if min(tcr(2,is),te(ic))-ts(ic)+te(ic)-max(tcr(1,ie),ts(ic))+ ...
       sum(diff(tcr(:,is+1:ie-1)))>tforb
     disp(['FATAL ERROR: finite cluster ',cl{ic}, ...
           ' power requested in TCV for too long at too low a plasma current'])
     tcv_recipe_error=2;
     return
    end
   end
% summator settings
   mdsdata(['_mod=\RHVPS_MOD_',cl{ic},':MODULATOR']);
   for kk=1:nsumm
    summ(nscat-1+ic,kk) = mdsdata(['(_mod[0]=="PCS"&&', ...
                                  'any(_mod[1:3]==_summn[$1]))'],kk-1);
   end
   for il=lau{ic}
    ils = num2str(il,'%02d');
    mdsdata(['_mod=\MIR_MOD_',ils,':MODULATOR']);
    for kk=1:nsumm
     summ(nslx2-1+il,kk) = mdsdata(['(_mod[0]=="PCS"&&', ...
                                   'any(_mod[1:3]==_summn[$1]))'],kk-1);
    end
   end
else
  tvc{ic} = tdef;vc{ic} = sdef;
  for il=gy{ic}
   tal{il} = tdef;al{il} = (sdef-avof(il))/avsl(il);
  end
  ts(ic) = 0;te(ic) = 0;
  ecintcv(gy{ic}) = false;
end
% write traces
wave_put(['\DRAW_FEEDFOR_ECRH:ALIM_',cl{ic}],vc{ic},tvc{ic})
% reference times
vdbw(strcat('tcv_protectdb::shot_design:ecrh_',cl{ic},'_',{'start','stop'}, ...
                '_ref'),[ts(ic),te(ic)])
%
% XU power
ic = 2;%cluster number
if xureq
% XU requested into TCV?
% at this stage, greq and ecintcv do not incorporate the information that the
%   RHVPS trace is nonzero
   for il = gy{ic}
    kl = gtol(il);
    greq(il) = strcmpi(mdsdata(['\online_gyro_',num2str(il,'%02d')]),'ON');
    ecintcv(il) = ~isnan(kl)&&(greq(il) && ...
                  strcmpi(mdsdata(['\switch_line_',num2str(kl,'%02d')]),'TCV'));
   end
% power traces
   power = mdsdata(['\POWER_CLUSTER_',cl{ic}]);
   vc{ic} = mdsdata(['_vc=max(min(\RHVPS_VOLT_',cl{ic}, ...
                ',maxval(\LIM_RHVPS_', ...
             cl{ic},':VMAX)),minval(\LIM_RHVPS_',cl{ic},':VMIN))']);
   vc{ic}(power<=0)=0;
   tvc{ic} = mdsdata('dim_of(_vc)');
   if length(vc{ic})<2,tvc{ic} = tdef;vc{ic} = sdef;end
   [vc{ic},tvc{ic},ts(ic),te(ic)] = massage_traces(vc{ic},tvc{ic},[],[], ...
        [],tdeact(ic),[],mdsdata(['\LIM_RHVPS_',cl{ic},':DV_DT_ERR']), ...
                mdsdata(['\LIM_RHVPS_',cl{ic},':DV_DT_MIN']), ...
                mdsdata(['\LIM_RHVPS_',cl{ic},':DV_DT_MAX']), ...
                mdsdata(['minval(\LIM_RHVPS_',cl{ic},':VMIN)']), ...
                mdsdata(['maxval(\LIM_RHVPS_',cl{ic},':VMAX)']));
   if ~strcmpi(mdsdata(['\RHVPS_MOD_',cl{ic},':MODULATOR[0]']),'INTERNAL')
     ts(ic) = tvc{ic}(1);
   end
   tsa(10) = max(ts(ic),mdsdata('\power_gate_10[0]'));
   tsa(11) = max(ts(ic),mdsdata('\power_gate_11[0]'));
   tea(10) = min(te(ic),mdsdata('\power_gate_10[1]'));
   tea(11) = min(te(ic),mdsdata('\power_gate_11[1]'));
   ts(ic) = min(tsa(10:11));
   te(ic) = max(tea(10:11));
% now greq and ecintcv incorporate the information that the power trace is
%   nonzero
   greq(gy{ic}) = greq(gy{ic})&any(vc{ic});
% timers
   if any(greq(gy{ic}))
    ecintcv(gy{ic}) = ecintcv(gy{ic})&any(vc{ic});
    greqxu = greq(gy{ic});
    vdbw(strcat(dbrxut,timstai(gy{2}(greqxu)),suft), ...
        tsa(gy{2}(greqxu))-tact(gy{2}(greqxu)))% anode inhibit on gates
    vdbw(strcat(dbrxut,timstoi(gy{2}(greqxu)),suft), ...
        tea(gy{2}(greqxu)))% anode inhibit off gates
    vdbw(strcat(dbrxut,timstag(gy{2}(greqxu)),suft), ...
        tsa(gy{2}(greqxu))-tact(gy{2}(greqxu))+tgate(gy{2}(greqxu)))% anode HV on gates
    vdbw(strcat(dbrxut,timstog(gy{2}(greqxu)),suft), ...
        tea(gy{2}(greqxu))-tgate(gy{2}(greqxu)))% anode HV off gates
    vdbw(strcat(dbrxut,{timstac{ic},timstoc{ic}},suft), ...
         [ts(ic)-max(tact(gy{2}(greqxu)))+tgates(ic),te(ic)+tgatee(ic)])% cathode gates
% pulse length check vs limit
    if te(ic)-ts(ic)>mdsdata(['\LIM_RHVPS_',cl{ic},':TMAX'])-0.09
      disp(['FATAL ERROR: cluster ',cl{ic}, ...
            ' power traces exceed the maximum allowed duration'])
      tcv_recipe_error=2;
      return
    end
% extend RHVPS trace
    tvc{ic}(1:2)=tvc{ic}(1:2)-max(tact(gy{2}(greqxu)))-tactr(2);
    tvc{ic}(end-1:end)=tvc{ic}(end-1:end)+max(tact(gy{2}(greqxu)))+tactr(2);
    ts(ic) = tvc{ic}(1);
    te(ic) = tvc{ic}(end);
   else
    ts(ic) = 0;
    te(ic) = 0;
    vdbw(strcat(dbrx2t,[timstac(ic),timstoc(ic), ...
                        timbods,timbode]),timcdef)
   end
% pulse length check vs Ip
   if any(ecintcv(gy{ic}))
    is = find(tcr(1,:)<=ts(ic),1,'last');
    ie = find(tcr(2,:)>=te(ic),1,'first');
    if min(tcr(2,is),te(ic))-ts(ic)+te(ic)-max(tcr(1,ie),ts(ic))+ ...
       sum(diff(tcr(:,is+1:ie-1)))>tforb
     disp(['FATAL ERROR: finite cluster ',cl{ic}, ...
           ' power requested in TCV for too long at too low a plasma current'])
     tcv_recipe_error=2;
     return
    end
   end
% summator settings
   mdsdata(['_mod=\RHVPS_MOD_',cl{ic},':MODULATOR']);
   for kk=1:nsumm
    summ(nscat-1+ic,kk) = mdsdata(['(_mod[0]=="PCS"&&', ...
                                  'any(_mod[1:3]==_summn[$1]))'],kk-1);
   end
   for il=lau{ic}
    ils = num2str(il,'%02d');
    mdsdata(['_mod=\MIR_MOD_',ils,':MODULATOR']);
    for kk=1:nsumm
     summ(nslx2-1+il,kk) = mdsdata(['(_mod[0]=="PCS"&&', ...
                                   'any(_mod[1:3]==_summn[$1]))'],kk-1);
    end
   end
else
  tvc{ic} = tdef;vc{ic} = sdef;
  for il=lau{ic}
   tal{il} = tdef;al{il} = (sdef-avof(il))/avsl(il);
  end
  ts(ic) = 0;te(ic) = 0;
  ecintcv(gy{ic}) = false;
end
% write traces
wave_put(['\DRAW_FEEDFOR_ECRH:ALIM_',cl{ic}],vc{ic},tvc{ic})
% reference times
vdbw(strcat('tcv_protectdb::shot_design:ecrh_',cl{ic},'_',{'start','stop'}, ...
                '_ref'),[ts(ic),te(ic)])
%
% X2 launchers
[aof,phbl,phi,phios] = deal(zeros(1,laumax));
for ic = [1,2]
 for il=lau{ic}
  ils1 = num2str(il,'%01d');
  avsl(il) = mdsdata(['_avsl=-2/USING(DATA(\CLUSTER_',cl{ic},'.LINE_', ...
                              ils1,':LAU_POL_SLOP),\TOP,-1,"ECRH")']);
  aof(il) = mdsdata(['_aof=USING(DATA(\CLUSTER_',cl{ic},'.LINE_',ils1, ...
                      ':LAU_POL_ZERO),\TOP,-1,"ECRH")']);
  avof(il) = mdsdata('_avof=-1*_avsl*_aof');
 end
end
% X2 toroidal angles (preliminary overshoot = 5 times the backlash)
for ic = [1,2]
 for il = lau{ic}
  kg = ltog(il);
  ils = num2str(il,'%02d');
  if ecintcv(kg)
   phbl(il) = abs(mdsdata(['\CALIB_LAUNCHER_',ils,':BACKLASH']));
   phi(il) = mdsdata(['max(min(\LAU_ANGLE_',ils, ...
               ',\LIM_LAU_ROT_',ils,':PHI_MAX),', ...
                '\LIM_LAU_ROT_',ils,':PHI_MIN)']);
   phios(il) = phi(il)-5*phbl(il)*mdsr(['\CALIB_LAUNCHER_',ils,':REF_DIR']);
% move if farther than twice the error (using backlash instead as error is
%   currently undefined)
   phir = vdbr([dbx2,'::x2_',num2str(il,'%01d'),':mes_tor']);
   if min(abs(phir-phi(il)),abs(phir-phios(il))) > 2*abs(phbl(il))
     vdbw([dbx2,'::x2_',num2str(il,'%01d'),':tor_val'],phios(il))
   end
  end
% mirror traces
  aerr = mdsdata(['\LIM_MIR_',ils,':TH_ERR']);
  almin = mdsdata(['_almin=\LIM_MIR_',ils,':TH_MIN']);
  almax = mdsdata(['_almax=\LIM_MIR_',ils,':TH_MAX']);
  dalmin = mdsdata(['\LIM_MIR_',ils,':DTH_DT_MIN']);
  dalmax = mdsdata(['\LIM_MIR_',ils,':DTH_DT_MAX']);
  al0 = mdsdata(['_al=max(min(\DESIGN_MIRROR_',ils,',_almax),_almin)']);
  tal0 = mdsdata('dim_of(_al)');
  if length(al0)<2,tal0 = tdef;al0 = aof(il)*ones(size(tal0));end
% ensure slow motion from rest position to starting position to avoid sudden
%   jerks (starts at -4 s)
  al0 = [-avof(il)/avsl(il);al0(1);al0;-avof(il)/avsl(il)];
  tal01 = hm_t0+4+5*hm_period_s;
  tal02 = max([tal0(1)-1,tal01+(al0(2)-al0(1))/dalmin, ...
              tal01+(al0(2)-al0(1))/dalmax]);
  if tal02>tal0(1)-0.001
    tal01 = max(tal01-tal02+tal0(1)-0.001,hm_t0+4+hm_period_s);
    tal02 = tal0(1)-0.001;
  end
  tal0 = [tal01;tal02;tal0;tal0(end)+1];
  [al{il},tal{il}] = massage_traces(al0,tal0,-avof(il)/avsl(il),[],[], ...
                      [],[],aerr/terr,dalmin,dalmax,almin,almax);
 end
end
% X2 toroidal angles (final setting)
for ic = [1,2]
 for il = lau{ic}
  kg = ltog(il);
  if ecintcv(kg)
    ils = num2str(il,'%01d');
    if abs(vdbr([dbx2,'::x2_',ils,':mes_tor'])-phi(il)) > 2*abs(phbl(il))
     n_to = 10;
     i_to = 0;
     while abs(vdbr([dbx2,'::x2_',ils,':mes_tor'])-phios(il)) > 2*abs(phbl(il))
      i_to = i_to+1;
      if i_to>=n_to
       disp(['FATAL ERROR: timed out on toroidal mirror ',int2str(il), ...
             ' movement to overshoot position'])
       tcv_recipe_error=2;
       return
      end
      pause(2)
     end
     vdbw([dbx2,'::x2_',ils,':tor_val'],phi(il))
     n_to = 20;
     i_to = 0;
     while abs(vdbr([dbx2,'::x2_',ils,':mes_tor'])-phi(il)) > 2*abs(phbl(il))
      i_to = i_to+1;
      if i_to>=n_to
       disp(['FATAL ERROR: timed out on toroidal mirror ',int2str(il), ...
             ' movement to final position'])
       tcv_recipe_error=2;
       return
      end
      pause(2)
     end
    end
  end
 end
end
%
% X3 top launchers
ic = 4;
for il = lau{ic}
 kg = ltog(il);
 if ~isnan(kg)&&ecintcv(kg)
  ils = num2str(il,'%02d');
% mirror traces
  aerr = mdsdata(['\LIM_MIR_',ils,':TH_ERR']);
  almin = mdsdata(['_almin=\LIM_MIR_',ils,':TH_MIN']);
  almax = mdsdata(['_almax=\LIM_MIR_',ils,':TH_MAX']);
  dalmin = mdsdata(['\LIM_MIR_',ils,':DTH_DT_MIN']);
  dalmax = mdsdata(['\LIM_MIR_',ils,':DTH_DT_MAX']);
  al0 = mdsdata(['_al=max(min(\MIR_ANGLE_',ils,',_almax),_almin)']);
  tal0 = mdsdata('dim_of(_al)');
  if length(al0)<2,tal0 = tdef;al0 = (almin+almax)/2*ones(size(tal0));end
  [al{il},tal{il}] = massage_traces(al0,tal0,(almin+almax)/2,0.5,[], ...
                      [],[],aerr/terr,dalmin,dalmax,almin,almax);
 else
  al{il}=[45.;45.];tal{il}=tdef;
 end
end
% write traces
for ic = [1,2]
 for il=lau{ic}
  wave_put(['\DRAW_FEEDFOR_ECRH:POL_MIR_0',num2str(il,'%02d')], ...
           avof(il)+avsl(il)*al{il},tal{il})
 end
end
ic=4;
for il=lau{ic}
 kg = ltog(il);
 if ~isnan(kg)&&ecintcv(kg)
  wave_put(['\DRAW_FEEDFOR_ECRH:POL_MIR_0',num2str(il,'%02d')],al{il},tal{il})
 end
end
%
ts0=Inf;te0=-Inf;
for ic=1:3
 if any(greq(gy{ic}))
  ts0=min(ts0,ts(ic));
  te0=max(te0,te(ic));
 end
end
if ts0>te0
  ts0=0;te0=0;
end
vdbw(strcat('tcv_protectdb::shot_design:ecrh_',{'start','stop'},'_ref'), ...
     [ts0,te0])

% write summator
if x2req || x3req
  mdsw('\SUMMATOR.SUMMER_001.INPUTS',int32(summ))
end

ts0=Inf;te0=-Inf;
for ic=1:3
 if any(ecintcv(gy{ic}))
  ts0=min(ts0,ts(ic));
  te0=max(te0,te(ic));
 end
end
if ts0>te0
  ts0=0;te0=0;
end

% density thresholds
if any(ecintcv)
  ipstart = vdbr('tcv_protectdb::shot_design:ip_start_ref');
  interlocktest = mdsr('\ECHCS_AUX_PARAM_04');
  tdt = [min([ipstart,ts0])-0.1,ipstart+0.04,ipstart+0.07,ipstart+0.1, ...
        ipstart+0.2,max([te0+0.1,ipstart+0.3])]';
  dt = [0.0,0.15,0.3,0.5,0.5,0.5];% lower
  wave_put('\DRAW_PHYS:PHYS_007',dt,tdt)
  dt = [40.0,20.0,10.0,2.0,1.7,1.7];% upper
  wave_put('\DRAW_PHYS:PHYS_008',dt,tdt)
%
  pcsdothresh(densloroot,'\DRAW_PHYS:PHYS_007',ts0,-1.,-1.,...
              interlocktest,[-1.,-1.,2.,2.],[-3.99,0.03,0.05,4.92]);
%
  if any(ecintcv([gy{1},gy{2}]))
    vign = 20.; vmin = 2.5;
  else% X3 only
    vign = 20.; vmin = 7.0;
  end
  pcsdothresh(denshiroot,'\DRAW_PHYS:PHYS_008',ts0,vign,vmin,...
              interlocktest,[20.,20.,1.5,1.5],[-3.99,0.03,0.05,4.92]);
% Output gain for density threshold channels (this must be consistent with gain
%   in pcsinscale
  threshgain = 0.1*2^(8-vdbr('diag2db::fir:fringe_gain'));
  mdsw('\HYBRID::OUTPUT_SCALES_THR_DENS','$',threshgain,'x')
% Timer (end gate)
  if interlocktest
    vdbw('timers_operdb::timer_D14:absolute',1.0)
  else
    vdbw('timers_operdb::timer_D14:absolute',te0+0.02)
  end
else
  wave_put('\DRAW_PHYS:PHYS_007',sdef,tdef)
  wave_put('\DRAW_PHYS:PHYS_008',sdef+1,tdef)
  wave_put(densloroot,[-2,-2]',[-3,4]')
  wave_put(denshiroot,[20,20]',[-3,4]')
  vdbw('timers_operdb::timer_D14:absolute',-0.5)
end
% stray detector gains
if any(ecintcv) && (te0>=-0.05 || te0-ts0<=0.003)
  vdbw([dbx2,'::x2_2:ch_003_pregain'],2)
  vdbw([dbx2,'::x2_4:ch_003_pregain'],1)
  vdbw([dbx2,'::x2_5:ch_003_pregain'],1)
  vdbw([dbx2,'::x2_6:ch_003_pregain'],1)
else
  vdbw([dbx2,'::x2_2:ch_003_pregain'],3)
  vdbw([dbx2,'::x2_4:ch_003_pregain'],3)
  vdbw([dbx2,'::x2_5:ch_003_pregain'],3)
  vdbw([dbx2,'::x2_6:ch_003_pregain'],3)
end

%%%%%***** function separator #####

% Set up NBH
function pcsnbhsetup
tdef = [0,2]'; sdef = [0,0]';% default traces
densloroot = '\DRAW_THRESH:NB_DENS_LOW';
denshiroot = '\DRAW_THRESH:NB_DENS_HIGH';

if strncmpi(vdbr('tcvpcdb::nbh:reqfordoshot'),'ON',2) && ...
   vdbr('TCV_NBHDB::NBI:LCS_MODE') == 9
  ts1 = vdbr('TIMERS_OPERDB::TIMER_439:ABSOLUTE')+1;
% TODO: could replace te with a true end of trace
  te1 = ts1+2;
else
  ts1 = Inf;
  te1 = -Inf;
end
if strncmpi(vdbr('tcvpcdb::nb2:reqfordoshot'),'ON',2) && ...
   vdbr('TCV_NBHDB::NB2:LCS_MODE') == 9
  ts2 = vdbr('TIMERS_OPERDB::TIMER_43A:ABSOLUTE')+1;
  te2 = ts2+2;
else
  ts2 = Inf;
  te2 = -Inf;
end
ts = min(ts1,ts2);
te = max(te1,te2);
if ~isinf(ts)
% allow 0.5 s at low density as wall can tolerate this
  tsd = ts+0.5;
% TODO: could do Ip check as in ECRH case once have full time trace of NBH
  ipstart = vdbr('tcv_protectdb::shot_design:ip_start_ref');
  tdt = [min([ipstart,tsd])-0.1,ipstart+0.04,ipstart+0.07,ipstart+0.1, ...
        ipstart+0.2,max([te+0.1,ipstart+0.3])]';
  dt = [0.0,0.15,0.3,0.5,0.5,0.5];% lower
  wave_put('\DRAW_PHYS:PHYS_009',dt,tdt)
  dt = [40.0,20.0,10.0,8.0,8.0,8.0];% upper
  wave_put('\DRAW_PHYS:PHYS_010',dt,tdt)
%
  pcsdothresh(densloroot,'\DRAW_PHYS:PHYS_009',tsd,-1.,0.5);
%
  pcsdothresh(denshiroot,'\DRAW_PHYS:PHYS_010',tsd,20.,2.5);
% Timers (gates)
  vdbw('timers_diagdb::timer_971:absolute',-3.5)
  vdbw('timers_diagdb::timer_972:absolute',3.5)
else
  wave_put('\DRAW_PHYS:PHYS_009',sdef,tdef)
  wave_put('\DRAW_PHYS:PHYS_010',sdef+1,tdef)
  wave_put(densloroot,[-2,-2]',[-3,4]')
  wave_put(denshiroot,[20,20]',[-3,4]')
  vdbw('timers_operdb::timer_971:absolute',-3.5)
end

%%%%%***** function separator #####

% Set up EC launcher/receiver (diagnostic actuator)
function pcseclaunrecsetup
global tcv_recipe_error
tdef = [0,2]'; sdef = [0,0]';% default traces
nsumm = 3;% summator branches
nsecl = 6;%index of summator channel
dbx2 = 'tcv_ecrh_x2db';
avsl = mdsdata(['_avsl=-2/USING(DATA(\CLUSTER_B.LINE_7:LAU_POL_SLOP),', ...
                     '\TOP,-1,"ECRH")']);
mdsdata(['_aof=USING(DATA(\CLUSTER_B.LINE_7:LAU_POL_ZERO),', ...
                     '\TOP,-1,"ECRH")']);
avof = mdsdata('_avof=-1*_avsl*_aof');
if vdbr([dbx2,'::x2_7:ask_select'])~=3
% master clock parameters
  dbrhm = 'tcv_protectdb::h_mere';
  hm_t0 = vdbr([dbrhm,':t0']);
  hm_period_f = vdbr([dbrhm,':period_f']);
  hm_n_slow = vdbr([dbrhm,':n_slow']);
  hm_period_s = hm_period_f*hm_n_slow*1.e-6 ;
% summator settings
  mdsdata('_summn=\SUMMATOR:INPUT_NAMES');
  summ = mdsdata('\SUMMATOR.SUMMER_001.INPUTS');
% toroidal angle (preliminary overshoot = 5 times the backlash)
  phbl = abs(mdsdata(['USING(DATA(\CALIB_LAUNCHER_D01:BACKLASH),', ...
                             '\TOP,-1,"DIAGZ")']));
  phi = mdsdata(['max(min(USING(DATA(\LAU_ANGLE_D01),\TOP,-1,"DIAGZ"),', ...
                  'USING(DATA(\LIM_LAU_ROT_D01:PHI_MAX),\TOP,-1,"DIAGZ")),', ...
                   'USING(DATA(\LIM_LAU_ROT_D01:PHI_MIN),\TOP,-1,"DIAGZ"))']);
  phios = phi-5*phbl* ...
           mdsr('USING(DATA(\CALIB_LAUNCHER_D01:REF_DIR),\TOP,-1,"DIAGZ")');
% move if farther than twice the error (using backlash instead as error is
%   currently undefined)
  phir = vdbr([dbx2,'::x2_7:mes_tor']);
  if min(abs(phir-phi),abs(phir-phios)) > 2*abs(phbl)
   m_to = 10;
   j_to = 0;
   while ~strcmpi(deblank(vdbr([dbx2,'::X2_7:MOT_TOR'])),'ON')
    n_to = 10;
    i_to = 0;
    while strcmpi(deblank(vdbr([dbx2,'::X2_7:ETAT_MOT'])),'ON')
% turn off poloidal motor and turn on toroidal motor
      i_to = i_to+1;
      if i_to>=n_to
        disp('FATAL ERROR: timed out on diagnostic poloidal mirror turn-off')
        tcv_recipe_error=2;
        return
      end
      vdbw([dbx2,'X2_7:WR_ETAT'],'OFF')
      pause(3)
    end
    j_to = j_to+1;
    if j_to>=m_to
      disp('FATAL ERROR: timed out on diagnostic toroidal mirror turn-on')
      tcv_recipe_error=2;
      return
    end
    vdbw([dbx2,'::X2_7:MOT_TOR_ON'],'ON')
    pause(3)
   end
   vdbw([dbx2,'::X2_7:TOR_VAL'],phios)
  end
% mirror trace
  aerr = mdsdata('USING(DATA(\LIM_MIR_D01:TH_ERR),\TOP,-1,"DIAGZ")');
  almin = mdsdata('_almin=USING(DATA(\LIM_MIR_D01:TH_MIN),\TOP,-1,"DIAGZ")');
  almax = mdsdata('_almax=USING(DATA(\LIM_MIR_D01:TH_MAX),\TOP,-1,"DIAGZ")');
  dalmin = mdsdata('USING(DATA(\LIM_MIR_D01:DTH_DT_MIN),\TOP,-1,"DIAGZ")');
  dalmax = mdsdata('USING(DATA(\LIM_MIR_D01:DTH_DT_MAX),\TOP,-1,"DIAGZ")');
  al0 = mdsdata(['max(min(USING(DATA(\DESIGN_MIRROR_D01),\TOP,-1,"DIAGZ"),', ...
                       '_almax),_almin)']);
  tal0 = mdsdata('USING(dim_of(\DESIGN_MIRROR_D01),\TOP,-1,"DIAGZ")');
% ensure slow motion from rest position to starting position to avoid sudden
%   jerks (starts at -4 s)
  al0 = [-avof/avsl;al0(1);al0;-avof/avsl];
  tal01 = hm_t0+4+5*hm_period_s;
  tal02 = max([tal0(1)-1,tal01+(al0(2)-al0(1))/dalmin, ...
              tal01+(al0(2)-al0(1))/dalmax]);
  if tal02>tal0(1)-0.001
    tal01 = max(tal01-tal02+tal0(1)-0.001,hm_t0+4+hm_period_s);
    tal02 = tal0(1)-0.001;
  end
  tal0 = [tal01;tal02;tal0;tal0(end)+1];
  [al,tal] = massage_traces(al0,tal0,-avof/avsl,[],[], ...
                      [],[],aerr/terr,dalmin,dalmax,almin,almax);
% summator settings
  mdsdata('_mod=USING(DATA(\MIR_MOD_D01:MODULATOR),\TOP,-1,"DIAGZ")');
  for kk=1:nsumm
   summ(nsecl,kk) = mdsdata(['(_mod[0]=="PCS"&&', ...
                                  'any(_mod[1:3]==_summn[$1]))'],kk-1);
  end
% toroidal angle (final setting)
  if abs(vdbr([dbx2,'::x2_7:mes_tor'])-phi) > 2*abs(phbl)
   m_to = 10;
   j_to = 0;
   while ~strcmpi(deblank(vdbr([dbx2,'::X2_7:MOT_TOR'])),'ON')
    n_to = 10;
    i_to = 0;
    while strcmpi(deblank(vdbr([dbx2,'::X2_7:ETAT_MOT'])),'ON')
% turn off poloidal motor and turn on toroidal motor
      i_to = i_to+1;
      if i_to>=n_to
        disp('FATAL ERROR: timed out on diagnostic poloidal mirror turn-off')
        tcv_recipe_error=2;
        return
      end
      vdbw([dbx2,'X2_7:WR_ETAT'],'OFF')
      pause(3)
    end
    j_to = j_to+1;
    if j_to>=m_to
      disp('FATAL ERROR: timed out on diagnostic toroidal mirror turn-on')
      tcv_recipe_error=2;
      return
    end
    vdbw([dbx2,'::X2_7:MOT_TOR_ON'],'ON')
    pause(3)
   end
   n_to = 10;
   i_to = 0;
   while abs(vdbr([dbx2,'::x2_7:mes_tor'])-phios) > 2*abs(phbl)
    i_to = i_to+1;
    if i_to>=n_to
     disp(['FATAL ERROR: timed out on diagnostic toroidal mirror ', ...
           'movement to overshoot position'])
     tcv_recipe_error=2;
     return
    end
    pause(2)
   end
   vdbw([dbx2,'::X2_7:TOR_VAL'],phi)
  end
  n_to = 20;
  i_to = 0;
  while abs(vdbr([dbx2,'::x2_7:mes_tor'])-phi) > 2*abs(phbl)
   i_to = i_to+1;
   if i_to>=n_to
    disp(['FATAL ERROR: timed out on diagnostic toroidal mirror ', ...
          'movement to final position'])
    tcv_recipe_error=2;
    return
   end
   pause(2)
  end
  n_to = 10;
  i_to = 0;
  while strcmpi(deblank(vdbr([dbx2,'::X2_7:MOT_TOR'])),'ON')
% turn off toroidal motor
   i_to = i_to+1;
   if i_to>=n_to
     disp('FATAL ERROR: timed out on diagnostic toroidal mirror turn-off')
     tcv_recipe_error=2;
     return
   end
   vdbw([dbx2,'X2_7:MOT_TOR_ON'],'OFF')
   pause(3)
  end
  n_to = 10;
  i_to = 0;
  while ~strcmpi(deblank(vdbr([dbx2,'::X2_7:ETAT_MOT'])),'ON')
% turn on poloidal motor
   i_to = i_to+1;
   if i_to>=n_to
     disp('FATAL ERROR: timed out on diagnostic poloidal mirror turn-on')
     tcv_recipe_error=2;
     return
   end
   vdbw([dbx2,'X2_7:WR_ETAT'],'ON')
   pause(3)
  end
% write summator
  mdsw('\SUMMATOR.SUMMER_001.INPUTS',summ)
else
  tal = tdef;al = (sdef-avof)/avsl;
end
% write trace
wave_put('\DRAW_FEEDFOR_DIAG:ECLR_MIR_001',avof+avsl*al,tal)

%%%%%***** function separator #####

% Load matrix switching sequences to real-time node
function pcsmatswitching

phroot = '\phys_mat_addresses:';
mats = {'a','g','m'};
mats_hard = {'a1','g1','m'};
bb_addr = [1,3,5];
pb_addr = {'2','3','4'};% private bus addresses
dt = 1e-2;%granularity

nmats = numel(mats);
[switch_times,mat_hard_addr_list,mat_hard_addr] = deal(cell(1,nmats));
[last_time,mhalen] = deal(zeros(1,nmats));

for i=1:nmats
% add one sample at start
  switch_times{i} = [-2;dt*fix((mdsdata('dim_of(getnci($1//$2,"RECORD"))', ...
                phroot,mats{i},mats{2},mats{3})+dt/2)/dt)];
  last_time(i) = max(switch_times{i});
  mdsdata('_mao=floor(0.5+getnci($1//$2,"RECORD"))',phroot,mats{i});
% physical addresses (only valid for low matrix indices)
  mat_hard_addr_list{i} = mdsdata('_mhal=8+getnci($1//$3//$2,"RECORD")', ...
                '\LOAD_MAT_','_ADDRESS',mats_hard{i});
  mat_hard_addr{i} = mdsdata('_mhal[_mao]');
  mat_hard_addr{i} = [mat_hard_addr{i}(1);mat_hard_addr{i}];
  mhalen(i) = length(mat_hard_addr{i});
end
timebase = dt*(-1:ceil(min(max(last_time/dt),396)/2)*2);% timebase starts at -dt; clamp to 398; should be even
tblen = length(timebase);
seq = 1:2*tblen;
seq_byte = 4*(0:fix((tblen+1)/2)-1)+1;
for i=1:nmats
  mat_seq = mat_hard_addr{i}(sum(repmat(timebase,mhalen(i),1) > ...
                                      repmat(switch_times{i},1,tblen)-dt/6))';
  mat_seq = [fix(mat_seq/256);mat_seq];
  mat_seq_byte = int32(uint8(mat_seq));% byte-switching
  mdsr('VDBPUT($1,$2)', ...
        ['TIMERS_OPERDB::SWITCH:LIST_',mats{i}], ...
        int32(mat_seq_byte(seq_byte)+ ...
                mat_seq_byte(seq_byte+1)*256+ ...
    (mat_seq_byte(seq_byte+2)+mat_seq_byte(seq_byte+3)*256)*65536));
  mdsr('r_download(3,$1,$2)',mat_seq(seq),4096*bb_addr(i));
  vdbw(['TCV_PROTECTDB::RT_NODE:PRIVATE_BUS_',pb_addr{i}], ...
        mat_hard_addr_list{i}(1))
end
vdbw('TIMERS_OPERDB::TIMER_013:ABSOLUTE',-dt)
mdsr('VDBPUT($1,$2)', ...
        strcat('TIMERS_OPERDB::SWITCH:',{'INTERVAL','SWITCHES','POINTER'}), ...
        [int32(dt*1e6),tblen-2,4098]');

%%%%%***** function separator #####

% Conditional settings
function pcscondsettings
recipetype = lower(deblank(vdbr('tcv_protectdb::shot_create:command_1')));
% these kinds of look-up tables could be globals, run at beginning of Matlab session?
recipes = {'fbte','stray','ecrhonly','nbhonly','backoff', ...
        'fpsonly','gasonly','diagonly','pcs_test', ...
        'timer_only'};
recind = find(contains(recipes,recipetype));
acq = {'ON','ON','OFF','ON','ON','ON','ON','ON','ON','OFF'};
pcs = {'ON','ON','OFF','OFF','ON','ON','OFF','ON','ON','OFF'};
magnspeed = {'FAST','FAST','FAST','FAST','FAST','FAST','FAST','FAST', ...
             'FAST','FAST'};
magnrate = [10,10,10,10,10,10,10,10,10,10];
gasvacuum = {'ON','OFF','OFF','OFF','OFF','OFF','ON','OFF','OFF','OFF'};
diags = {'ON','ON','OFF','OFF','OFF','OFF','OFF','ON','OFF','OFF'};
gas = [1,1,0,0,0,0,1,0,0,0];
gastimers = {[-8,-8,-8,-8,-8,-8],[-5.9,4.0,-5.9,4.0,-5.9,4.0]};
vdbw('TCV_PROTECTDB::TIMER:SELECT','ON')
vdbw('TCV_PROTECTDB::ACQ:SELECT',acq(recind))
vdbw('TCV_PROTECTDB::PCS:SELECT',pcs(recind))
vdbw(strcat('TIMERS_OPERDB::TIMER_58',{'1','2','3','4','5','6'}, ...
                ':ABSOLUTE'),gastimers{gas(recind)+1})
vdbw('TCV_PROTECTDB::ACQ:BASE_RATE',magnspeed(recind))
vdbw('TCV_PROTECTDB::SEQ:COMMAND', ...
     ['Acquisition ',int2str(magnrate(recind)),'kHz'])
vdbw(strcat('TCV_PROTECTDB::',{'GAS','VIDE'},':SELECT'), ...
     repmat(gasvacuum(recind),2,1))
vdbw('TCV_PUBLICDB::DIAGS:SELECT',diags(recind))
% special cases
% bolometer gains in stray
if recind==2
 vdbw(strcat('TCV_PUBLICDB::BOLO_GAIN:CH_',iii(1:64)),repmat(5,64,1))
end
% magnetic acquisition end
vdbw(strcat('TCV_PUBLICDB::',{'TRCF_MAG_BPOL_001','TRCF_MAG_FLUX_001', ...
            'TRCF_MAG_LOOP_001','TRCF_MAG_CURR_001','TRCH_MAG_FFPS_001'}), ...
     repmat(0.4,1,5))
% MHD acq for non-plasma shot (reset elsewhere for plasma shot)
if recind~=1
  magnsuff = {'DELTA_T_START','DELTA_T_STOP','REF_START','REF_STOP'};
  magnval = {-0.04,0.04,'POL_START_REF','POL_STOP_REF'};
  vdbw(strcat('TCV_PUBLICDB::TRCH_MHD_VFAST_001:',magnsuff(1:2)), ...
        cell2mat(magnval(1:2))')
  vdbw(strcat('TCV_PUBLICDB::TRCH_MHD_VFAST_001:',magnsuff(3:4)),magnval(3:4)')
  vdbw(strcat('TCV_PUBLICDB::TRCH_MHD_FAST_001:',magnsuff(1:2)), ...
        cell2mat(magnval(1:2))')
  vdbw(strcat('TCV_PUBLICDB::TRCH_MHD_FAST_001:',magnsuff(3:4)),magnval(3:4)')
  vdbw(strcat('TCV_PUBLICDB::TRCH_MHD_SLOW_001:',magnsuff(1:2)), ...
        cell2mat(magnval(1:2))')
  vdbw(strcat('TCV_PUBLICDB::TRCH_MHD_SLOW_001:',magnsuff(3:4)),magnval(3:4)')
  vdbw('TCV_PUBLICDB::TRCH_MHD_SLOW_001:FREQUENCE',5000)
end
% PD
vdbw(strcat('TCV_PUBLICDB::TRCF_PD_002:',{'REF_START','REF_STOP'}), ...
     {'POL_START_REF','POL_STOP_REF'}')
vdbw('TIMERS_DIAGDB::TIMER_114:ABSOLUTE',-0.1)

%%%%%***** function separator #####

% Special timers
function pcsspectimers

mdsopen('BASE',-1)
% pressure slaves
tim=mdsr('getnci("\\BASE::PRESS_SLAVES.TCV_*", "MINPATH")')';
for i=1:length(tim)
  t0 = mdsr('getnci($//":T0","RECORD")',tim{i});
  timn = mdsr('getnci($//":TIMER","RECORD")',tim{i});
  vdbw(['TIMERS_OPERDB::TIMER_',timn,':ABSOLUTE'],t0)
end
mdsclose;

%%%%%***** function separator #####

% Prepare data acquisition
function pcsmdsprepare

global tcv_recipe_error
recipetype = lower(deblank(vdbr('tcv_protectdb::shot_create:command_1')));
% these kinds of look-up tables could be globals, run at beginning of Matlab session?
recipes = {'fbte','stray','ecrhonly','nbhonly','backoff', ...
        'fpsonly','gasonly','diagonly','pcs_test', ...
        'timer_only'};
recind = find(contains(recipes,recipetype),1);
mdsopen('TCV_SHOT',-1)
vdbw('TCV_PROTECTDB::CYCLE:PREPARE_MDS_RECEIVED','ON')
vdbw('TCV_PROTECTDB::CYCLE:PREPARE_MDS_COMPLETE','OFF')
recacq = logical([1,1,0,0,1,1,0,1,0,0]);
exps = {'AXUV',     'LANG',     'MPX',      '',         '', ...
        '',         '',         '',         '',         '',         '', ...
        'BOLO',     'XTOMO',    '',         '',         'TWOMM', ...
        'DIAG_NEU', 'BREMSIR',  '',         '',         '', ...
        '',         '',         ''};% empty string = sel not used
ndts = {'axuv',     'lang',     'mpx',      'ece',      'se', ...
        'xtep',     'rthom',    'tpci',     'hxrs',     'mag',      'mhd', ...
        'bolo',     'xtomo',    'xte',      'xmode',    '1mm', ...
        'neutron',  'pm',       'msi',      'mix',      'ltcc', ...
        'frp',      'cher',     'hyb'};
trees = {'atlas','base','diagz'};
branches = {'tcpip:dt_','camac:tf_','camac:th_','tcpip.'};
actions = {{'INIT','VERIF'},{'ARM'},{'STORE'},{'FIRE'}};
dbs = {'TCV_PROTECT','TCVPC','TCV_PUBLIC'};
sels = {'SELECT','ACQ_SELECT'};
tree = [1,      1,      1,      1,      1, ...
        1,      1,      1,      1,      1,      1, ...
        2,      1,      2,      2,      2, ...
        2,      2,      3,      1,      1, ...
        1,      1,      1];
bra =  [1,      1,      1,      1,      1, ...
        1,      1,      1,      1,      1,      1, ...
        2,      1,      2,      3,      3, ...
        3,      2,      4,      1,      1, ...
        1,      1,      1];
act =  {1,      1,      1,      1,      1, ...
        1,      1,      1,      1,      1,      1, ...
        [2,3],  1,      [2,3],  [2,3],  [2,3], ...
        [2,3],  [2,3],  4,      1,      1, ...
        1,      1,      1};% index into actions array
db =   [2,      2,      3,      3,      3, ...
        3,      0,      0,      0,      0,      0, ...
        2,      2,      0,      0,      2, ...
        3,      2,      0,      0,      0, ...
        0,      0,      0];
sel =  [1,      1,      2,      0,      0, ...
        0,      0,      0,      0,      0,      0, ...
        1,      1,      0,      0,      1, ...
        1,      1,      0,      0,      0, ...
        0,      0,      0];% index into sels array
dgsel= [0,      0,      0,      0,      0, ...
        0,      0,      0,      0,      0,      0, ...
        1,      1,      0,      0,      1, ...
        1,      1,      1,      0,      0, ...
        0,      0,      0];% dependence on general diagnostic selection
acq1 = [1,      1,      1,      1,      1, ...
        1,      1,      1,      1,      1,      1, ...
        1,      1,      1,      1,      1, ...
        1,      1,      1,      1,      1, ...
        1,      1,      1];% selection defaults for recacq recipes when sel not used (-1 or 2 to override Vsystem selection towards 0 or 1, respectively)
acq2 = [0,      0,      0,      0,      0, ...
        0,      0,      0,      0,      1,      1, ...
        0,      0,      0,      0,      0, ...
        0,      0,      0,      1,      1, ...
        0,      0,      1];% selection defaults for ~recacq recipes when sel not used (-1 or 2 to override Vsystem selection towards 0 or 1, respectively)
isel = find(~cellfun('isempty',exps));
selstr = strcmp(sels(sel(isel)),'SELECT');
% turn off non-power-supply ATLAS nodes for certain recipes (with exceptions
%   defined by acq2), turn them on for all others (before dealing with
%   individual named diagnostics)
dts = deblank(mdsdata( ...
        '_dts=getnci("\\atlas::top.hardware.tcpip:dt_*","FULLPATH")'));
if ~recacq(recind)
 acq=acq2;
 for j=1:length(dts)
   if ~contains(upper(dts{j}),{'DT_MT','DT_FPS','DT_RED','DT_X2','DT_X3',...
                               'DT_AFPS','DT_GYRO','DT_NBHP'}) && ...
      mdsdata('getnci(_dts[$1],"ON")[0]',j-1)
     mdstcl(['set node/off ',dts{j},':*_ACTION']);
   end
 end
 acq(strcmp(trees(tree),'atlas')&~acq) = -1;%override Vsystem selection
else
 acq=acq1;
 for j=1:length(dts)
% We don't touch the ECRH acquisition at all (this is handled in ecrh_select)
% We also don't touch the NBH protection acquisition (handled separately below)
  if ~contains(upper(dts{j}),{'DT_X2','DT_X3','DT_AFPS','DT_GYRO','DT_NBHP'})
   if mdsdata('getnci(_dts[$1],"ON")[0]',j-1)
     mdstcl(['set node/on ',dts{j},':*_ACTION']);
   end
  end
 end
end
% have to deal separately with different channels so 'ON'/'OFF' is turned to 1/0
acqtmp(selstr) = strncmpi(vdbr(strcat(dbs(db(isel(selstr))),'DB::', ...
                       exps(isel(selstr)),':',sels(sel(isel(selstr))))),'ON',2);%acqtmp is a logical array
acqtmp(~selstr) = vdbr(strcat(dbs(db(isel(~selstr))),'DB::', ...
                       exps(isel(~selstr)),':',sels(sel(isel(~selstr)))));%values different from zero are turned into true
acqtmp(acq(isel)==-1) = 0;
acqtmp(acq(isel)==2) = 1;
acq(isel) = acqtmp;%logical turns into double
for i=1:length(ndts)
  dts = cellstr(mdsr(['_dts=getnci("\\',trees{tree(i)}, ...
                 '::top.hardware.',branches{bra(i)},ndts{i},'*","MINPATH")']));
  if acq(i)==1 && ...
     (~dgsel(i) || strncmpi(vdbr('TCV_PUBLICDB::DIAGS:SELECT'),'ON',2))% acquisition
    for j=1:length(dts)
      if mdsdata('getnci(_dts[$1],"ON")[0]',j-1)
       for k=act{i}
        for k1=1:length(actions{k})
         nodtmp=[dts{j},':',actions{k}{k1},'_ACTION'];
         if mdsdata(['node_exists("\',nodtmp,'")'])
          mdstcl(['set node/on ',nodtmp]);
         end
        end
       end
      end
    end
  else% no acquisition
    for j=1:length(dts)
      if mdsdata('getnci(_dts[$1],"ON")[0]',j-1)
        mdstcl(['set node/off ',dts{j},':*_ACTION']);
      end
    end
  end
end
% special cases
mdstcl('set node/off \atlas::dt4g_mt_001:verif_action');
mdstcl('set node/off \atlas::dt4g_thcpl_001:verif_action')
mdstcl('set node/off \atlas::dt4g_test_001:verif_action');
if vdbr('TCV_NBHDB::NBI:LCS_MODE') == 9 || vdbr('TCV_NBHDB::NB2:LCS_MODE') == 9
  mdstcl('set node/on \ATLAS::DT4G_NBHP_001:INIT_ACTION');
  mdstcl('set node/on \ATLAS::DT4G_NBHP_001:VERIF_ACTION');
else
  mdstcl('set node/off \ATLAS::DT4G_NBHP_001:INIT_ACTION');
  mdstcl('set node/off \ATLAS::DT4G_NBHP_001:VERIF_ACTION');
end
ecrhacq=mdsr('ecrh_select()');% ECRH
if ecrhacq~=1
  disp('FATAL ERROR: ECRH acquisition setup failed')
  tcv_recipe_error=2;
  return
end
mdsevent('prepare_done');
vdbw('TCV_PROTECTDB::CYCLE:PREPARE_MDS_COMPLETE','ON')


%%%%%***** function separator #####

% Align timebases
function [ya,tc] = pcsalignipol(y,t)
dt = 2e-3;
tc = unique(round(sort(cell2mat(t))/dt)*dt);
ntraces = numel(y);
ya = zeros(numel(tc),ntraces);
for i=1:ntraces
  [ti,js] = unique(t{i});
  ya(:,i) = interp1(ti,y{i}(js),tc,'linear',0);
end

%%%%%***** function separator #####
% Calculate feedforward voltages based on the feedforward currents in hybrid
%   mode
function pcshybcurrtovolt(coil_list,curr,t_curr,scalfac)
% Q: make dt into a global?
dt = 1e-3;
% mutuals
n_oh=2;n_e=8;n_f=8;
mutuals = mdsr('static("\\mut_a_a")');
mutuals = mutuals(1:n_e+n_f+n_oh,1:n_e+n_f+n_oh);
coils = [repmat({'OH'},1,n_oh),repmat({'E'},1,n_e),repmat({'F'},1,n_f)];
coilsort=cell(1,size(coil_list,2));
for k=1:size(coil_list,2)
  sorttemp=find(strcmp(coils,coil_list{1,k}));
  coilsort{k}=sorttemp(1:coil_list{2,k});
end
coilsort = horzcat(coilsort{:}); % Avoids growing inside a loop
n_coil = length(coilsort);
selfcorrfac = [0.75,0.87, ...
               0.76,0.76,0.76,0.75,0.75,0.82,0.80,0.72, ...
               0.83,0.78,0.80,0.85,0.76,0.80,0.79,0.89];
oh1corrfac = [1,0.50,repmat(0.60,1,n_e),repmat(0.30,1,n_f)];
oh2corrfac = [0.72,1,[0.65,0.70,0.70,0.70,0.70,0.70,0.70,0.65], ...
              repmat(0.50,1,n_f)];
efcorrfac = repmat(0.7,16,18);
corrfac = [oh1corrfac;oh2corrfac;efcorrfac];
mutuals = mutuals.*(1-eye(size(mutuals))+diag(selfcorrfac)).*corrfac;
mutuals = mutuals(coilsort,coilsort);
curr_dot = diff(curr)./repmat(diff(t_curr'),1,n_coil);
volt = curr_dot*mutuals'*scalfac;
volt = [zeros(1,n_coil);repelem(volt,2,1);zeros(1,n_coil)];
t_volt = repelem(t_curr',2,1);
t_volt(1:2:end-1,:)=t_volt(2:2:end,:)-dt;
drawstr='\DRAW_FEEDFOR_';
k1=0;
for ck=1:size(coil_list,2)
  for k=1:coil_list{2,ck}
    wpstr = [drawstr,coil_list{1,ck},'_U:ALIM_',iii(k)];
    wave_put(wpstr,volt(:,k1+k),t_volt)
  end
  k1=k1+coil_list{2,ck};
end

%%%%%***** function separator #####

% Utility functions

function [x,tx,tsta,tend] = massage_traces(x,tx,xlo,tramp,tact,tdeact,ttail,sigder,dermin,dermax,xmin,xmax)
% general trace massaging, tweaking, plus checks on limits
if nargin<12 || isempty(xmax)   || isnan(xmax),   xmax  = max(x);end %signal maximum
if nargin<11 || isempty(xmin)   || isnan(xmin),   xmin  = min(x);end %signal minimum
if nargin<10 || isempty(dermax) || isnan(dermax), dermax= Inf;   end %maximum acceptable derivative
if nargin<9  || isempty(dermin) || isnan(dermin), dermin=-Inf;   end %minimum acceptable derivative
if nargin<8  || isempty(sigder) || isnan(sigder), sigder= 0;     end %minimum significant derivative
if nargin<7  || isempty(ttail)  || isnan(ttail),  ttail = 0;     end %extra time to "off" trigger
if nargin<6  || isempty(tdeact) || isnan(tdeact), tdeact= 0;     end %empirical deactivation delay
if nargin<5  || isempty(tact)   || isnan(tact),   tact  = 0;     end %empirical activation delay
if nargin<4  || isempty(tramp)  || isnan(tramp),  tramp = 0.001; end %ramp-up and ramp-down time
if nargin<3  || isempty(xlo)    || isnan(xlo),    xlo   = 0;     end %baseline signal value
limfac = 0.2;% compression factor of limit derivatives in the proximity of signal limits
limlev = 0.1;% definition of "proximity"
tx = [tx(1)-tramp;tx;tx(end)+tramp];
x = [xlo;x;xlo];
iv = find(x~=xlo);
if ~isempty(iv)
 iv = [min(iv)-1,min(iv):max(iv),max(iv)+1]';
else
 iv = (2:length(x)-1)';
 tact = 0;
 tdeact = 0;
end
x = x(iv);
tx = tx(iv);
tsta = tx(2);
tend = tx(end)+ttail;
tx(1:2) = tx(1:2)-tact;
tx(end-1:end) = tx(end-1:end)+tdeact;
[tx,iu] = unique(tx);
x = x(iu);
dt = diff(tx);
dxdt = diff(x)./dt;
dermax = repmat(dermax,size(dxdt));
dermin = repmat(dermin,size(dxdt));
ilim = find(x(2:end)-xmin<limlev*(xmax-xmin)|xmax-x(2:end)<limlev*(xmax-xmin));
dermax(ilim) = dermax(ilim)*limfac;
dermin(ilim) = dermin(ilim)*limfac;
if any(dxdt>dermax | dxdt<dermin)
 for k=2:length(x)
   if (x(k)-x(k-1))/dt(k-1)>dermax(k-1),x(k)=x(k-1)+dt(k-1)*dermax(k-1);end
   if (x(k)-x(k-1))/dt(k-1)<dermin(k-1),x(k)=x(k-1)+dt(k-1)*dermin(k-1);end
 end
end
if sigder
  iv = mdsdata('erp($1,$2,$3,1)',tx,x,sigder)+1;
  x = x(iv);
  tx = tx(iv);
end

%%%%%***** function separator #####

function refs_osc (ch,amp)
if nargin<1,ch=[];end
if nargin<2 || isempty(amp),amp=0.001;end
times = 0:0.05:2;
if isempty(ch)
% get all channels
  nodes = mdsr('GETNCI("\\DRAW_REFS:REF_*","PATH")')';
  ch = 1:length(nodes);
else
% SMO
  nodes{ch} = ['\DRAW_REFS:REF_',iii(ch)];
end
% SMO
for i=ch
  toplim = mdsr([nodes{i},':LIMITS']);
  wave_put(nodes{i}, ...
       (toplim(1)+toplim(2))/2+cos(2*i*times)*(toplim(2)-toplim(1))*amp, ...
       times);
end

%%%%%***** function separator #####

function pcsdothresh(root,sig,tmin,vign,vmin,interlocktest,ilockdata,ilocktime)
% root is the node containing the absolute threshold values
% sig is the node containing the relative threshold values
if nargin > 5 && interlocktest
  if nargin < 8, error('pcsdothresh:missingarg','If interlocktest is true, values and time base must be provided'); end
  expr = 'MAKE_SIGNAL($1,*,$2)';
  args = {mdscvt(ilockdata,'f'),mdscvt(ilocktime,'f')};
else
% The expression depends on 3 values:
%  $1: minimum starting time
%  $2: values outside of EC/NB activation values (vign)
%  $3: minimum accepted value (vmin)
  expr = ['_th = ',sig,';', ...
          '_tr = DATA(DIM_OF(_th));', ...
          '_de = \DRAW_REFS:REF_021;', ...
          '_td = DATA(DIM_OF(_de));', ...
          '_te = min(_tr[SIZE(_tr)-1],_td[SIZE(_td)-1]);', ...
          '_ts = max(_tr[0],_td[0],0.035,$1);', ...
          '_tn = union(min(max([_tr,_td],_ts),_te));', ...
          '_th = resample_cn(_th, _tn);_de = resample_cn(_de, _tn);', ...
          'Make_Signal([$2,$2,max(DATA(_th) * DATA(_de),$3),$2,$2], *, ', ...
                      '[-3.99,max(0.03,$1-0.01),_tn,_tn[SIZE(_tn)-1]+0.01,max(_tn[SIZE(_tn)-1]+0.02,4.92)])'];
  args = {mdscvt(tmin,'f'),mdscvt(vign,'f'),mdscvt(vmin,'f')};
end
for it=0:7
  mdsw([root,':TRACK_',num2str(it,'%03d')],expr,args{:},'x')
end
mdsw([root,':TIME'],['dim_of(',root,')'],'x')
mdsw([root,':CONTROL'],'Build_Signal(DRAW_KNOTS, *, TIME)','x')
mdsw([root,':DRAW_KNOTS'],'$',uint8(ones(1001,1)),'x')

%%%%%***** function separator #####

% Data server utility functions
function pal_on_off(db,chan)

% function pal_on_off(db,chan)
%
% "Presses" PAL bouton (ON followed by OFF, allowing for appropriate prompter
%    and poller cycle); returns only when readback can be guaranteed to be
%    correct or times out after 1 s

%       S. Coda, 19/04/04
%       updated  08/06/16

chixp = mdsr('vdbchix($)',[db,'::poller:loops']);
vdbw([db,'::',chan],'ON');
% clear and then book wake (must do at each call to guarantee change has come
%   *after* any actuation command given)
% wait one poller cycle or 200 ms prompter cycle before giving PAL off command
mdsvalue(['vdbsync3(-',int2str(chixp),'),vdbsync3(',int2str(chixp), ...
               '),vdbsync3(,250)']);
vdbw([db,'::',chan],'OFF');
% wait one poller cycle before attempting readback, time out after 1 sec
mdsvalue(['vdbsync3(-',int2str(chixp),'),vdbsync3(',int2str(chixp),')']);
mdsvalue('vdbsync3(,1000)');

% clear wake
mdsvalue(['vdbsync3(-',int2str(chixp),')']);

%%%%%***** function separator #####

function button_with_timeout(db,ch,condch,cond,n_to,t_pa)

condst = {'OFF','ON'};% for cond=0 or 1
i_to = 0;
while ~strcmpi(deblank(vdbr(condch)),condst{cond+1})
  i_to = i_to+1;
  if i_to>=n_to,return;end
  pal_on_off(db,ch)
  pause(t_pa)
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

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

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

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)

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

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

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

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
