From 3163f3460506135243e52495b4140b0a99eddd21 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Sat, 15 Jan 2022 23:03:43 +0100
Subject: [PATCH 01/15] fix gdat_plot, mex for jet, and other small

---
 matlab/JET/gdat_jet.m                     |  8 ++++
 matlab/JET/rda_jet.m                      |  6 ++-
 matlab/TCV/gdat_tcv.m                     | 12 ++---
 matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 55 +++++++++++++++++------
 matlab/gdat_plot.m                        |  2 +-
 5 files changed, 62 insertions(+), 21 deletions(-)

diff --git a/matlab/JET/gdat_jet.m b/matlab/JET/gdat_jet.m
index 2197e8f..4bb12dd 100644
--- a/matlab/JET/gdat_jet.m
+++ b/matlab/JET/gdat_jet.m
@@ -574,7 +574,15 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
         end
       else
         % assume same time for all
+        % mds_dir = '/home/amerle/public/mds-ssh';
+        mds_dir = '/home/amerle/public/mds-ssh3';
+        if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
+          addpath(mds_dir);
+        end
 	mdsconnect(['ssh://' gdat_params.jet_user '@mdsplus.jetdata.eu']);
+        if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
+          rmpath(mds_dir);
+        end
         rda_data_request = ['ppf/kk3/' num2str(i,'te%.2d')];
         bb = mdsvalue(['_bb=jet("' rda_data_request '",' num2str(shot) ')']);
         if ~isempty(bb) && numel(bb)==size(aa.data,2)
diff --git a/matlab/JET/rda_jet.m b/matlab/JET/rda_jet.m
index 1d93ae6..f80caa5 100644
--- a/matlab/JET/rda_jet.m
+++ b/matlab/JET/rda_jet.m
@@ -34,7 +34,8 @@ end
 if usemdsplus
   % use mdsplus
 
-  mds_dir = '/home/amerle/codes/mds-ssh';
+  % mds_dir = '/home/amerle/public/mds-ssh';
+  mds_dir = '/home/amerle/public/mds-ssh3';
   if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
     addpath(mds_dir);
   end
@@ -170,6 +171,9 @@ if usemdsplus
   if ~unix('test -d /home/duval/mdsplus')
     rmpath('/home/duval/mdsplus')
   end
+  if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
+    rmpath(mds_dir);
+  end
 
 else
   % use RDA
diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index c82eb8d..0bef612 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -348,11 +348,13 @@ if do_mdsopen_mdsclose
     end
     mapping_for_tcv.expression = mapping_for_tcv.expression(7:end);
   elseif shot==-1 || liuqe_version_eff==-1
-    ishot = mdsopen('pcs', shot);
+    [ishot,shot_stat] = mdsopen('pcs', shot);
   else
-    ishot = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
+    [ishot,shot_stat] = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
   end
-  if isempty(ishot) || ishot~=shot
+  if mod(shot_stat,2)==0
+    ishot
+    warning('error opening shot');
     if gdat_params.nverbose>=1; warning(['cannot open shot= ' num2str(shot)]); end
     return
   end
@@ -2515,7 +2517,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
     %
     if any(strmatch('nbi2',gdat_data.gdat_params.source))
-    % NB2  
+    % NB2
     nodenameeff = '\results::NB2:POWR_TCV';
       nb2_data_tdi = tdi(nodenameeff);
       if ~isempty(nb2_data_tdi.data) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
@@ -2543,7 +2545,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       end
     end
     if any(strmatch('dnbi',gdat_data.gdat_params.source))
-    % NB2  
+    % NB2
     nodenameeff = '\RESULTS::DNBI:POWR_TCV';
       nb2_data_tdi = tdi(nodenameeff);
       if ~isempty(nb2_data_tdi.data) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index df66620..54f872b 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -107,6 +107,7 @@ global_quantities_desc.volume = params_eff.data_request;
 params_eff.data_request = 'w_mhd';
 global_quantities.w_mhd = gdat(params_equilibrium.shot,params_eff);
 global_quantities_desc.w_mhd = params_eff.data_request;
+ids_equilibrium_description.global_quantities = global_quantities_desc;
 
 global_quantities_fieldnames = fieldnames(global_quantities);
 special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
@@ -214,9 +215,10 @@ for it=1:ntime
     ids_equilibrium.time_slice{it}.boundary.x_point = {};
   end
 end
+ids_equilibrium_description.boundary = boundary_desc;
+ids_equilibrium_description.boundary_temp = temp_desc;
 
 %% constraints
-% TODO: Add description
 
 % Measured values
 liuqe_time = ids_equilibrium.time; % Not necessarily magnetics time so far
@@ -224,9 +226,15 @@ mag_time   = mdsvalue('\magnetics::bpol_003:dim0');
 itime = iround_os(mag_time, liuqe_time);
 mag_time_req = mdscvt(mag_time(itime),'f');
 bpol = mdsvalue('\magnetics::bpol_003[$1,*]',mag_time_req);
+nbpol = size(bpol,2);
+constraints_desc.bpol_probe.measured = sprintf('from %s i=1:%d','\magnetics::bpol_003[i,*]',nbpol);
 flux = mdsvalue('tcv_idealloop("flux")[$1,*]',mag_time_req);
-diam = mdsvalue('\results::dmlcor[$1]',mag_time_req);
+nflux = size(flux,2);
+constraints_desc.flux_loop.measured = sprintf('from %s i=1:%d','tcv_idealloop("flux")[i,*]',nflux);
 ip   = mdsvalue('\magnetics::iplasma:trapeze[$1]',mag_time_req);
+constraints_desc.ip.measured = sprintf('from %s','\magnetics::iplasma:trapeze');
+diam = mdsvalue('\results::dmlcor[$1]',mag_time_req);
+constraints_desc.diamagnetic_flux.measured = sprintf('from %s','\results::dmlcor');
 % Coil currents since dim of constraints pf_current is IDS:pf_active/coil
 dim_pol = {'OH_001','OH_002','OH_002','OH_002','OH_002','OH_002','OH_002',...
            'E_001','E_002','E_003','E_004','E_005','E_006','E_007','E_008',...
@@ -236,26 +244,38 @@ ipol = mdsvalue('\magnetics::ipol[$1,$2]',mag_time_req,dim_pol);
 ipol(:,27:29) = -ipol(:,27:29); % Bottom G-coils
 dim_pol(30:32) = {'TOR_001'};
 ipol(:,30:32) = 0; % TOR_001 is not used in LIUQE
+nipol = size(ipol,2);
+constraints_desc.pf_current.measured = sprintf('from %s, ii from %s','\magnetics::ipol[t,ii]',cell2str(dim_pol,15));
 
 % Reconstructed values
 ipol_liuqe_order = [18,19*ones(1,6),1:16,17*ones(1,6)]; % LIUQE order is E F G OH
 bpol_liuqe = mdsvalue('tcv_eq("b_probe","liuqe.m")');
+constraints_desc.bpol_probe.reconstructed = sprintf('%s','tcv_eq("b_probe","liuqe.m")');
 flux_liuqe = mdsvalue('tcv_eq("psi_loop","liuqe.m")');
-diam_liuqe = mdsvalue('tcv_eq("tor_flux_dml","liuqe.m")');
+constraints_desc.flux_liuqe.reconstructed = sprintf('%s','tcv_eq("psi_loop","liuqe.m")');
 ip_liuqe   = mdsvalue('tcv_eq("i_pl","liuqe.m")');
+constraints_desc.ip.reconstructed = sprintf('%s','tcv_eq("i_pl","liuqe.m")');
+diam_liuqe = mdsvalue('tcv_eq("tor_flux_dml","liuqe.m")');
+constraints_desc.diamagnetic_flux.reconstructed = sprintf('%s','tcv_eq("tor_flux_dml","liuqe.m")');
 ipol_liuqe = mdsvalue('tcv_eq("i_pol","liuqe.m")');
 ipol_liuqe = ipol_liuqe(ipol_liuqe_order,:);
 ipol_liuqe(27:29,:) = -ipol_liuqe(27:29,:); % Bottom G-coils
 ipol_liuqe(30:32,:) = 0; % ... TOR
+constraints_desc.ipol_liuqe.reconstructed = sprintf('%s','tcv_eq("i_pol","liuqe.m")');
 
 % Weights (using old parameters tree for now)
 bpol_err = mdsvalue('\results::parameters:berr')./mdsvalue('\results::parameters:vvv[0:37]');
+constraints_desc.bpol_probe.weight = sprintf('(%s/%s)^2','\results::parameters:vvv[0:37]','\results::parameters:berr');
 flux_err = mdsvalue('\results::parameters:ferr')./mdsvalue('\results::parameters:www[0:37]')*2*pi;
-diam_err = 0.13e-3./mdsvalue('\results::parameters:idml');
+constraints_desc.flux_loop.weight = sprintf('(%s/%s)^2','\results::parameters:www[0:37]','\results::parameters:ferr');
 ip_err   = mdsvalue('\results::parameters:plcerr')*1e3;
+constraints_desc.ip.weight = sprintf('(1e-3/%s)^2','\results::parameters:plcerr');
+diam_err = 0.13e-3./mdsvalue('\results::parameters:idml');
+constraints_desc.diamagnetic_flux.weight = sprintf('(%s/0.13e-3)^2','\results::parameters:idml');
 ipol_err = mdsvalue('\results::parameters:cerr')./mdsvalue('\results::parameters:uuu[0:18]')*1e3;
 ipol_err = ipol_err(ipol_liuqe_order);
 ipol_err(30:32) = NaN;
+constraints_desc.pf_current.weight = sprintf('(1e-3*%s/%s)^2','\results::parameters:uuu[0:18]','\results::parameters:cerr');
 
 if ntime > 0
   constraints_orig = ids_equilibrium.time_slice{1}.constraints;
@@ -263,15 +283,16 @@ if ntime > 0
   ununsed_constraints = {'faraday_angle','mse_polarisation_angle','iron_core_segment',...
                          'n_e','n_e_line','pressure','q','x_point'};
   for name = ununsed_constraints, constraints_orig.(name{1})={}; end
+  constraints_desc.ununsed_constraints = ununsed_constraints;
 end
+keyboard
 for it = 1:ntime
   constraints = constraints_orig;
   % bpol_probe
-  nbpol = size(bpol,2);
   bpol_probe(1:nbpol) = constraints.bpol_probe(1);
   for ib = 1:nbpol
     bpol_probe{ib}.measured = bpol(it,ib);
-    bpol_probe{ib}.source = sprintf('IDS:magnetics/bpol_probe[%02d]/field',ib);
+    bpol_probe{ib}.source = sprintf('IDS:equilibrium/../constraints/bpol_probe[%02d]',ib);
     bpol_probe{ib}.time_measurement = mag_time(itime(it));
     bpol_probe{ib}.exact = 0;
     bpol_probe{ib}.weight = 1/(bpol_err(ib)).^2;
@@ -279,11 +300,10 @@ for it = 1:ntime
   end
   constraints.bpol_probe = bpol_probe;
   % flux_loop
-  nflux = size(flux,2);
   flux_loop(1:nflux) = constraints.flux_loop(1);
   for il = 1:nflux
     flux_loop{il}.measured = flux(it,il);
-    flux_loop{il}.source = sprintf('IDS:magnetics/flux_loop[%02d]/flux',il);
+    flux_loop{il}.source = sprintf('IDS:equilibrium/../constraints/flux_loop[%02d]',il);
     flux_loop{il}.time_measurement = mag_time(itime(it));
     flux_loop{il}.exact = 0;
     flux_loop{il}.weight = 1/(flux_err(il)).^2;
@@ -292,24 +312,23 @@ for it = 1:ntime
   constraints.flux_loop = flux_loop;
   % ip
   constraints.ip.measured = ip(it);
-  constraints.ip.source = 'IDS:magnetics/method[1]/ip';
+  constraints.ip.source = 'IDS:equilibrium/../constraints/ip';
   constraints.ip.time_measurement = mag_time(itime(it));
   constraints.ip.exact = 0;
   constraints.ip.weight = 1/(ip_err).^2;
   constraints.ip.reconstructed = ip_liuqe(it);
   % diamagnetic_flux
   constraints.diamagnetic_flux.measured = diam(it);
-  constraints.diamagnetic_flux.source = 'IDS:magnetics/method[1]/diamagnetic_flux';
+  constraints.diamagnetic_flux.source = 'IDS:equilibrium/../constraints/diamagnetic_flux';
   constraints.diamagnetic_flux.time_measurement = mag_time(itime(it));
   constraints.diamagnetic_flux.exact = 0;
   constraints.diamagnetic_flux.weight = 1/(diam_err).^2;
   constraints.diamagnetic_flux.reconstructed = diam_liuqe(it);
   % pf_current
-  nipol = size(ipol,2);
   pf_current(1:nipol) = constraints.pf_current(1);
   for ic = 1:nipol
     pf_current{ic}.measured = ipol(it,ic);
-    pf_current{ic}.source = sprintf('IDS:pf_active/coil[%02d]/current',ic);
+    pf_current{ic}.source = sprintf('IDS:equilibrium/../constraints/pf_current[%02d]',ic);
     pf_current{ic}.time_measurement = mag_time(itime(it));
     if strcmp(dim_pol{ic},'TOR_001')
       pf_current{ic}.source = [pf_current{ic}.source,' replaced with 0'];
@@ -465,6 +484,8 @@ for it=1:ntime
   end
 end
 
+ids_equilibrium_description.profiles_1d = profiles_1d_desc;
+
 % special cases
 for it=1:ntime
   ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.b_field_tor = ids_equilibrium.time_slice{it}.profiles_1d.f(1) ...
@@ -496,9 +517,13 @@ profiles_2d.grid_type.description = 'Cylindrical R,Z ala eqdsk';
 params_eff.data_request = 'psi';
 profiles_2d.psi = gdat(params_equilibrium.shot,params_eff); % add psi_bound in a second step in special cases
 profiles_2d_desc.psi = [params_eff.data_request ' adding psi_bound in a 2nd step'];
-% r = gdat(params_equilibrium.shot,'r','machine',gdat_params.machine); % not to be filled since in grid.dim1
+%profiles_2d.r = profiles_2d.psi;
+profiles_2d.r.data = repmat(repmat(profiles_2d.psi.dim{1},1,65),1,1,1299);
+profiles_2d_desc.r = 'from dim{1} of ''psi'' repeated';
+profiles_2d.z.data = repmat(repmat(profiles_2d.psi.dim{2}',28,1),1,1,1299);
+profiles_2d_desc.z = 'from dim{2} of ''psi'' repeated';
+
 % theta = gdat(params_equilibrium.shot,'theta','machine',gdat_params.machine);
-% z = gdat(params_equilibrium.shot,'z','machine',gdat_params.machine); % not to be filled since in grid.dim2
 
 profiles_2d_fieldnames = fieldnames(profiles_2d);
 special_fields = {'grid', 'grid_type'}; % fields needing non-automatic treatments
@@ -529,6 +554,8 @@ for it=1:ntime
       global_quantities.psi_boundary.data(it);
 end
 
+ids_equilibrium_description.profiles_2d = profiles_2d_desc;
+
 % make arrays not filled in empty:
 ids_equilibrium.grids_ggd = {};
 for it=1:numel(ids_equilibrium.time_slice)
diff --git a/matlab/gdat_plot.m b/matlab/gdat_plot.m
index 3334870..aba5eac 100644
--- a/matlab/gdat_plot.m
+++ b/matlab/gdat_plot.m
@@ -124,7 +124,7 @@ if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(
     zoom on;
   end
   maxnblines = 1;
-  if strcmp(gdat_data.gdat_params.data_request,'powers') && length(ab)==length(gdat_data.label)
+  if any(strcmp(gdat_data.gdat_params.data_request,'powers')) && length(ab)==length(gdat_data.label)
     % keep legend from label
   else
     for i=1:numel(linehandles)
-- 
GitLab


From 6b4736512550ff690c26fc3ccc07d9789d2fd2ed Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Sun, 16 Jan 2022 09:00:35 +0100
Subject: [PATCH 02/15] forgot an mdsopen call

---
 matlab/TCV/gdat_tcv.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 0bef612..86f4461 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -342,7 +342,7 @@ if do_mdsopen_mdsclose
   %%%  if liuqe_version_eff==-1
   if ~iscell(mapping_for_tcv.expression) && any(strfind(mapping_for_tcv.expression,'\rtc::'))
     mdsdisconnect;
-    ishot = mdsopen('rtc',shot); % sub-tree
+    [ishot,shot_stat]  = mdsopen('rtc',shot); % sub-tree
     if any(strfind(data_request_eff,'\rtc::'))
       data_request_eff = data_request_eff(7:end);
     end
-- 
GitLab


From df1666ab464748ab15e3b205d8e295b8bbe0a026 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 15:22:09 +0100
Subject: [PATCH 03/15] add gas to AUG, add Tags to gdat_plot

---
 matlab/AUG/aug_requests_mapping.m |  3 ++
 matlab/AUG/gdat_aug.m             | 80 +++++++++++++++++++++++++++++++
 matlab/gdat_plot.m                | 40 ++++++++++++++--
 3 files changed, 118 insertions(+), 5 deletions(-)

diff --git a/matlab/AUG/aug_requests_mapping.m b/matlab/AUG/aug_requests_mapping.m
index 9b4e16a..4112db7 100644
--- a/matlab/AUG/aug_requests_mapping.m
+++ b/matlab/AUG/aug_requests_mapping.m
@@ -130,6 +130,9 @@ switch lower(data_request)
   mapping.gdat_timedim = 2;
   mapping.method = 'switchcase'; % could use function make_eqdsk directly?
   mapping.expression = '';
+ case {'gas', 'gas_valve'}
+  mapping.gdat_timedim = 2;
+  mapping.method = 'switchcase';
  case 'halpha'
   mapping.timedim = 1;
   mapping.label = 'Halpha';
diff --git a/matlab/AUG/gdat_aug.m b/matlab/AUG/gdat_aug.m
index 3cc7d98..58d292a 100644
--- a/matlab/AUG/gdat_aug.m
+++ b/matlab/AUG/gdat_aug.m
@@ -1297,6 +1297,86 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
     gdat_data.dimunits{1} = 'rhopolnorm';
     gdat_data.dimunits{2} = 'time [s]';
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+   case {'gas', 'gas_valve', 'gas_valves'}
+    %
+    uvs_signals_source = [{'D_tot','Sum of calibrated D2 gas-fluxes (= D/s)'}; ...
+                    {'N_tot','Sum of calibrated N2 gas-fluxes (= N*7/s)'}; ...
+                    {'H_tot','Sum of calibrated H gas-fluxes (= H/s)'}; ...
+                    {'He_tot','Sum of calibrated He gas-fluxes (= He*2/s)'}; ...
+                    {'Ne_tot','Sum of calibrated Ne gas-fluxes (= Ne*10/s)'}; ...
+                    {'Ar_tot','Sum of calibrated Ar gas-fluxes (= Ar*18/s)'}];
+    for i=1:size(uvs_signals_source,1)
+      uvs_signals.(uvs_signals_source{i,1}) = gdat_aug(shot,{'UVS',uvs_signals_source{i,1}});
+      if isempty(gdat_data.t) && prod(isfinite(uvs_signals.(uvs_signals_source{i,1}).t)) && isvector(uvs_signals.(uvs_signals_source{i,1}).t) ...
+          && numel(uvs_signals.(uvs_signals_source{i,1}).t)>10
+        gdat_data.t = uvs_signals.(uvs_signals_source{i,1}).t; % all have same time base, could use time signal
+      end
+      if prod(isfinite(uvs_signals.(uvs_signals_source{i,1}).data)) && isvector(uvs_signals.(uvs_signals_source{i,1}).data) ...
+          && numel(uvs_signals.(uvs_signals_source{i,1}).data) == numel(gdat_data.t);
+        gdat_data.data(i,:) = uvs_signals.(uvs_signals_source{i,1}).data; % all have same time base
+      else
+        gdat_data.data(i,:) = NaN;
+      end
+      gdat_data.units{i} = uvs_signals_source{i,2};
+    end
+    gdat_data.data_fullpath = 'from UVS and uvs_signals_source';
+    gdat_data.uvs_signals_source = uvs_signals_source;
+    gdat_data.label = uvs_signals_source(:,1)';
+    gdat_data.dim{2} = gdat_data.t;
+    gdat_data.dimunits{2} = 's';
+    gdat_data.dim{1} = [1:size(gdat_data.data,1)];
+    gdat_data.dimunits{1} = 'uvs_signals_source index';
+    gdat_data.x = gdat_data.dim{1};
+    gdat_data.mapping_for.aug.timedim = 1;
+
+    % note also available in UVS, not commented the same way on ISIS
+    uvd_signals_source=[{'CFCu01T','PV01','empty comment'}; ...
+                        {'CFCo02A','PV02','valve flux, normally for prefill'}; ...
+                        {'CFFo02A','PV03','valve flux,'}; ...
+                        {'CFA13B','PV04','valve flux, in-vessel pipe for gas imaging'}; ...
+                        {'CFA03B','PV05','valve flux,'}; ...
+                        {'CFFo07A','PV06','valve flux,'}; ...
+                        {'CFDu05X','PV07','valve flux, normally with N2'}; ...
+                        {'CFDu05B','PV08','valve flux, normally with D2'}; ...
+                        {'CFCo10A','PV09','valve flux, in-vessel pipe into ICRH limiter'}; ...
+                        {'CFDu01B','PV10','valve flux, normally with D2'}; ...
+                        {'CFDu09B','PV11','valve flux, normally with D2'}; ...
+                        {'CFFo10A','PV12','valve flux,'}; ...
+                        {'CFFo14A','PV13','valve flux,'}; ...
+                        {'CFA03C','PV14','valve flux, short feeding pipe for rare gases'}; ...
+                        {'CFA13A','PV15','valve flux,'}; ...
+                        {'CFA03A','PV16','valve flux,'}; ...
+                        {'CFDu09X','PV17','valve flux, normally with N2'}; ...
+                        {'CFDu13X','PV18','valve flux, normally with N2'}; ...
+                        {'CFDu13B','PV19','valve flux, normally with D2'}; ...
+                        {'CFDu01X','PV20','valve flux, normally with N2'}];
+    gdat_data.pvxx.t = [];
+    gdat_data.pvxx.t = [];
+
+    for i=1:size(uvd_signals_source,1)
+      uvd_signals.(uvd_signals_source{i,2}) = gdat_aug(shot,{'UVD',uvd_signals_source{i,1}});
+      if isempty(gdat_data.pvxx.t) && prod(isfinite(uvd_signals.(uvd_signals_source{i,2}).t)) && isvector(uvd_signals.(uvd_signals_source{i,2}).t) ...
+          && numel(uvd_signals.(uvd_signals_source{i,2}).t)>10
+        gdat_data.pvxx.t = uvd_signals.(uvd_signals_source{i,2}).t; % all have same time base
+      end
+      if prod(isfinite(uvd_signals.(uvd_signals_source{i,2}).data)) && isvector(uvd_signals.(uvd_signals_source{i,2}).data) ...
+          && numel(uvd_signals.(uvd_signals_source{i,2}).data) == numel(gdat_data.pvxx.t);
+        gdat_data.pvxx.data(i,:) = uvd_signals.(uvd_signals_source{i,2}).data; % all have same time base
+      else
+        gdat_data.pvxx.data(i,:) = NaN;
+      end
+      gdat_data.pvxx.units{i} = [uvd_signals_source{i,2} ': ' uvd_signals_source{i,3}];
+    end
+    gdat_data.pvxx.data_fullpath = 'from UVD and uvd_signals_source';
+    gdat_data.pvxx.uvd_signals_source = uvd_signals_source;
+    gdat_data.pvxx.label = uvd_signals_source(:,2)';
+    gdat_data.pvxx.dim{2} = gdat_data.t;
+    gdat_data.pvxx.dimunits{2} = 's';
+    gdat_data.pvxx.dim{1} = [1:size(gdat_data.pvxx.data,1)];
+    gdat_data.pvxx.dimunits{1} = 'PVxx index';
+    gdat_data.pvxx.x = gdat_data.dim{1};
+
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ids'}
     params_eff = gdat_data.gdat_params;
diff --git a/matlab/gdat_plot.m b/matlab/gdat_plot.m
index aba5eac..acd7768 100644
--- a/matlab/gdat_plot.m
+++ b/matlab/gdat_plot.m
@@ -37,6 +37,7 @@ elseif isfield(gdat_data,'gdat_params') && isfield(gdat_data.gdat_params,'doplot
   end
 end
 if doplot==0; return; end
+redo_legend_from_Tags = 0;
 
 if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(gdat_data.t)
   fighandle = get(0,'CurrentFigure');
@@ -106,25 +107,54 @@ if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(
     else
       hylabel=ylabel(ylabel_eff,'interpreter','none');
     end
-    if iscell(gdat_data.label) && length(gdat_data.label)>=2;
-      ab=get(gca,'children');
+    ab=get(gca,'children');
+    if iscell(gdat_data.label) && numel(gdat_data.label) > 1;
+      if numel(ab)==numel(gdat_data.label) || mod(numel(ab),numel(gdat_data.label))==0
+        % Assumes a single shot with several lines and labels, if mod==0 then this is a new shot, only Tag present lines
+        for iab=1:numel(gdat_data.label)
+          set(ab(iab),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label{end-iab+1}]);
+        end
+        if numel(ab)>numel(gdat_data.label), redo_legend_from_Tags = 1; end
+      end
       if iscell(gdat_data.label) && length(ab) < length(gdat_data.label)
         if length(ab) == 1
           % assume combine label is best
           tempaaa = sprintf('%s/',gdat_data.label{:});
           hhhleg=legend(tempaaa(1:end-1));
+          set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' tempaaa(1:min(10,numel(tempaaa)-1))]);
         else
+          % not sure why would get there
           hhhleg=legend(gdat_data.label{1:length(ab)});
         end
-      else
+      elseif numel(ab)==numel(gdat_data.label)
         hhhleg=legend(gdat_data.label);
       end
-      set(hhhleg,'Interpreter','none');
+      if exist('hhhleg'), set(hhhleg,'Interpreter','none'); end
+    else
+      % assume one line per call
+      if isempty(gdat_data.label)
+        set(ab(1),'Tag',[num2str(gdat_data.shot)]);
+      elseif ischar(gdat_data.label)
+        % assume one signal, take max 10 chars
+        set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label(1:min(10,numel(gdat_data.label)))]);
+      elseif iscell(gdat_data.label) && numel(gdat_data.label) == 1
+        if isempty(gdat_data.label{1})
+          set(ab(1),'Tag',[num2str(gdat_data.shot)]);
+        else
+          set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label{1}(1:min(10,numel(gdat_data.label{1})))]);
+        end
+      end
+    end
+    if redo_legend_from_Tags
+      for iab=1:numel(ab)
+        % Use Tag for DisplayName, when overlay plots of multiple data at this stage
+        set(ab(iab),'DisplayName',strrep(get(ab(iab),'Tag'),'_','\_'));
+      end
     end
     zoom on;
   end
   maxnblines = 1;
-  if any(strcmp(gdat_data.gdat_params.data_request,'powers')) && length(ab)==length(gdat_data.label)
+  if redo_legend_from_Tags || any(strcmp(gdat_data.gdat_params.data_request,'powers')) || (numel(ab)==numel(gdat_data.label) && numel(ab)>1)
     % keep legend from label
   else
     for i=1:numel(linehandles)
-- 
GitLab


From e030263300f5e3c965a6b76bacaaafc7b079ee22 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 16:41:14 +0100
Subject: [PATCH 04/15] include changes for NBI from Matteo to ease rebase

---
 matlab/TCV/gdat_tcv.m             | 39 ++++++++++++++++++++---
 matlab/TCV_IMAS/tcv_get_ids_nbi.m | 53 ++++++++++++++++++++-----------
 2 files changed, 69 insertions(+), 23 deletions(-)

diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 86f4461..0b760d6 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1,4 +1,4 @@
-function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_request,varargin)
+function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv_OS(shot,data_request,varargin)
 %
 % function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
 %
@@ -2466,6 +2466,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     %
     if any(strmatch('nbi1',gdat_data.gdat_params.source))
       nodenameeff = '\results::NBH:POWR_TCV';
+      if shot<70811
+        nodenameeff = '\results::NBH:POWR_TCV';
+      end
       nbh_data_tdi = tdi(nodenameeff);
       if ~isempty(nbh_data_tdi.data) && ~ischar(nbh_data_tdi.data) && ~isempty(nbh_data_tdi.dim)
         nbi_neutral_power_tot = nbh_data_tdi.data.*1e6; % in W
@@ -2489,6 +2492,20 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi1.t(ij),gdat_data.nbi1.data(ij),gdat_data.t);
         gdat_data.x(end+1) = size(gdat_data.data,2);
         gdat_data.label{end+1}=gdat_data.nbi1.label;
+      else
+        if gdat_params.nverbose>=3
+          fprintf('No NB1 data in MDS+ nodes for shot %i\n', shot);
+          model_signame = '\ATLAS::NB1.DATA.MODEL:POWR_NEUTRAL';
+          if shot<70811
+            nodenameeff = '\ATLAS::NB1.DATA.MODEL:POWR_NEUTRAL';
+          end
+          model_data=tdi(model_signame);
+          if isempty(model_data.dim)
+            fprintf('And it was NOT requested \n');
+          else
+            fprintf('And it was requested. Check if plasma was a blip or if NB1 did not work \n');
+          end
+        end
       end
     end
     %
@@ -2517,10 +2534,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
     %
     if any(strmatch('nbi2',gdat_data.gdat_params.source))
-    % NB2
-    nodenameeff = '\results::NB2:POWR_TCV';
+      % NB2
+      nodenameeff = '\results::NB2:POWR_TCV';
       nb2_data_tdi = tdi(nodenameeff);
-      if ~isempty(nb2_data_tdi.data) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
+      if ~isempty(nb2_data_tdi.dim) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
         nbi_neutral_power_tot = nb2_data_tdi.data.*1e6; % in W
         ij = nbi_neutral_power_tot<100;
         nbi_neutral_power_tot(ij) = 0.;
@@ -2542,8 +2559,20 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi2.t(ij),gdat_data.nbi2.data(ij),gdat_data.t);
         gdat_data.x(end+1) = size(gdat_data.data,2);
         gdat_data.label{end+1}=gdat_data.nbi2.label;
+      else
+        if gdat_params.nverbose>=3
+          fprintf('No NB2 data in MDS+ nodes for shot %i\n', shot);
+          model_signame = '\ATLAS::NB2.DATA.MODEL:POWR_NEUTRAL';
+          model_data=tdi(model_signame);
+          if isempty(model_data.dim)
+            fprintf('And it was NOT requested \n');
+          else
+            fprintf('And it was requested. Check if plasma was a blip or if NB1 did not work \n');
+          end
+        end
       end
     end
+
     if any(strmatch('dnbi',gdat_data.gdat_params.source))
     % NB2
     nodenameeff = '\RESULTS::DNBI:POWR_TCV';
@@ -3286,7 +3315,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         gdat_data.x = [1:size(sig,1)];
         gdat_data.dimunits = {'20 chords per camera'; 's'};
       else
-        keyboard
+        % keyboard
         % extract only given channels
         gdat_data.data = sig(channel_xtomo,:);
         gdat_data.x = channel_xtomo;
diff --git a/matlab/TCV_IMAS/tcv_get_ids_nbi.m b/matlab/TCV_IMAS/tcv_get_ids_nbi.m
index 3542cfb..eb6777a 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_nbi.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_nbi.m
@@ -23,45 +23,59 @@ ids_nbi_description='';
 %
 
 
-nb_units = 2; % assume 2 units: 1st NBH and DNBI
+nb_units = 3; % assume 2 units: 1st NBH and DNBI
 ids_nbi.unit(1:nb_units) = ids_nbi.unit(1); % copy empty structure for all units, then fill in
 
 % create lists of what is different for each units so that can scan through units
-unit_identifier = {'NBH1', 'DNBI'};
-unit_name = {'25keV 1st NBH source', 'diagnostic NBI'};
-results_subname = {'nbh', 'dnbi'};
-species.a = [2., 1.];
-species.z_n = [1., 1.];
-species.label = {'D', 'H'};
-beamlets_group.direction = [-1, 1];
-beamlets_group.tangency_radius = [736, 235.3]*1e-3;
-beamlets_group.angle = [0., 0.];
-beamlets_group.width_horizontal = [250, 87.2]*1e-3;
-beamlets_group.width_vertical = [250, 87.2]*1e-3;
+unit_identifier = {'NB1', 'NB2', 'DNBI'};
+unit_name = {'25keV 1st NBH source', '50keV 2nd NBH source', 'diagnostic NBI'};
+results_subname = {'nb1', 'nb2', 'dnbi'};
+if shot<70811
+  results_subname = {'nbh', 'nb2', 'dnbi'};
+end
+species.a = [2., 2., 1.];
+species.z_n = [1., 1., 1.];
+species.label = {'D', 'D','H'};
+beamlets_group.direction = [-1, 1, 1];
+beamlets_group.tangency_radius = [736, 736, 235.3]*1e-3;
+beamlets_group.angle = [0., 0., 0.];
+beamlets_group.width_horizontal = [250, 250, 87.2]*1e-3;
+beamlets_group.width_vertical = [250, 250, 87.2]*1e-3;
 
 beamlets_group.focus(1:nb_units)=struct('focal_length_horizontal',[],'focal_length_vertical',[],'width_min_horizontal',[],'width_min_vertical',[]);
 beamlets_group.focus(1).focal_length_horizontal = 3.76;
 beamlets_group.focus(1).focal_length_vertical = 3.98;
 beamlets_group.focus(1).width_min_horizontal = 21.6*1e-2;
 beamlets_group.focus(1).width_min_vertical = 9.4*1e-2;
-beamlets_group.focus(2).focal_length_horizontal = 1.8;
-beamlets_group.focus(2).focal_length_vertical = 1.8;
-beamlets_group.focus(2).width_min_horizontal = 12.1*1e-2;
-beamlets_group.focus(2).width_min_vertical = 12.1*1e-2;
+% So far NB2 parameters are merely a copy of NB1 parameters
+beamlets_group.focus(2).focal_length_horizontal = 3.76;
+beamlets_group.focus(2).focal_length_vertical = 3.98;
+beamlets_group.focus(2).width_min_horizontal = 21.6*1e-2;
+beamlets_group.focus(2).width_min_vertical = 9.4*1e-2;
+beamlets_group.focus(3).focal_length_horizontal = 1.8;
+beamlets_group.focus(3).focal_length_vertical = 1.8;
+beamlets_group.focus(3).width_min_horizontal = 12.1*1e-2;
+beamlets_group.focus(3).width_min_vertical = 12.1*1e-2;
 
 beamlets_group.divergence(1:nb_units) = struct('particle_fraction',[],'vertical',[],'horizontal',[]);
 beamlets_group.divergence(1).particle_fraction = 1.;
 beamlets_group.divergence(1).vertical = 0.59 *pi/180.;
 beamlets_group.divergence(1).horizontal = 1.4 *pi/180.;
 beamlets_group.divergence(2).particle_fraction = 1.;
-beamlets_group.divergence(2).vertical = 0.53 *pi/180.;
-beamlets_group.divergence(2).horizontal = 0.53 *pi/180.;
+beamlets_group.divergence(2).vertical = 0.59 *pi/180.;
+beamlets_group.divergence(2).horizontal = 1.4 *pi/180.;
+beamlets_group.divergence(3).particle_fraction = 1.;
+beamlets_group.divergence(3).vertical = 0.53 *pi/180.;
+beamlets_group.divergence(3).horizontal = 0.53 *pi/180.;
 
 %dcd_NBH = psitbxdcd(4.5889, 0.0, 211.9535*pi/180, 0.0, -9.2308*pi/180);
 beamlets_group.position(1:nb_units) = struct('phi',[],'r',[],'z',[]);
 beamlets_group.position(1).phi = 211.9535*pi/180.;
 beamlets_group.position(1).r = 4.5889;
 beamlets_group.position(1).z = 0.;
+beamlets_group.position(2).phi = 58.8255*pi/180.;
+beamlets_group.position(2).r = 4.5889;
+beamlets_group.position(2).z = 0.;
 beamlets_group.position(2).phi = 295.2416*pi/180.;
 beamlets_group.position(2).r = 4.9274;
 beamlets_group.position(2).z = 0.;
@@ -73,6 +87,9 @@ for iunit=1:nb_units
   %% power
   params_eff.data_request = ['\results::' results_subname{iunit} ':powr_tcv'];
   pow=gdat_tcv(shot,params_eff);
+  if ischar(pow.data)
+    pow.data=0;
+  end
   ids_nbi.unit{iunit}.power_launched.data = pow.data*1e6;
   ids_nbi.unit{iunit}.power_launched.time = pow.t;
   ids_nbi_description.unit{iunit}.power_launched = params_eff.data_request;
-- 
GitLab


From 36d4bfe56646366714bea7299ee7a2ef105b7a8a Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 17:19:15 +0100
Subject: [PATCH 05/15] try to finalize tcv_get_ids_equilibrium changes

to be tested
---
 matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index 54f872b..10dcc0c 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -340,7 +340,7 @@ for it = 1:ntime
     end
   end
   constraints.pf_current = pf_current;
-  
+
   ids_equilibrium.time_slice{it}.constraints = constraints;
 end
 
@@ -518,9 +518,9 @@ params_eff.data_request = 'psi';
 profiles_2d.psi = gdat(params_equilibrium.shot,params_eff); % add psi_bound in a second step in special cases
 profiles_2d_desc.psi = [params_eff.data_request ' adding psi_bound in a 2nd step'];
 %profiles_2d.r = profiles_2d.psi;
-profiles_2d.r.data = repmat(repmat(profiles_2d.psi.dim{1},1,65),1,1,1299);
+profiles_2d.r.data = repmat(repmat(profiles_2d.psi.dim{1},1,numel(profiles_2d.psi.dim{2})),1,1,numel(profiles_2d.psi.dim{3}));
 profiles_2d_desc.r = 'from dim{1} of ''psi'' repeated';
-profiles_2d.z.data = repmat(repmat(profiles_2d.psi.dim{2}',28,1),1,1,1299);
+profiles_2d.z.data = repmat(repmat(profiles_2d.psi.dim{2}',numel(profiles_2d.psi.dim{1}),1),1,1,numel(profiles_2d.psi.dim{3}));
 profiles_2d_desc.z = 'from dim{2} of ''psi'' repeated';
 
 % theta = gdat(params_equilibrium.shot,'theta','machine',gdat_params.machine);
-- 
GitLab


From 2e85bf0c09ec580fa73cfccda2fb32274379bd23 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 17:43:35 +0100
Subject: [PATCH 06/15] add constraints description

tested on TCV
---
 matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index 10dcc0c..2638760 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -285,7 +285,6 @@ if ntime > 0
   for name = ununsed_constraints, constraints_orig.(name{1})={}; end
   constraints_desc.ununsed_constraints = ununsed_constraints;
 end
-keyboard
 for it = 1:ntime
   constraints = constraints_orig;
   % bpol_probe
@@ -344,6 +343,7 @@ for it = 1:ntime
   ids_equilibrium.time_slice{it}.constraints = constraints;
 end
 
+ids_equilibrium_description.constraints = constraints_desc;
 
 %
 %% profiles_1d (cannot use eqdsk since not same radial mesh)
-- 
GitLab


From 0a7fb2d6df4a9c3c4efbaa77f3e5c5f9d8d3c36e Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 17:50:33 +0100
Subject: [PATCH 07/15] temp

---
 matlab/TCV/gdat_tcv.m | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 0b760d6..b2a647e 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1,4 +1,4 @@
-function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv_OS(shot,data_request,varargin)
+function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_request,varargin)
 %
 % function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
 %
@@ -2574,8 +2574,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
 
     if any(strmatch('dnbi',gdat_data.gdat_params.source))
-    % NB2
-    nodenameeff = '\RESULTS::DNBI:POWR_TCV';
+      % NB2
+      nodenameeff = '\RESULTS::DNBI:POWR_TCV';
       nb2_data_tdi = tdi(nodenameeff);
       if ~isempty(nb2_data_tdi.data) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
         nbi_neutral_power_tot = nb2_data_tdi.data.*1e6; % in W
@@ -2706,7 +2706,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
    case {'phi_tor', 'phitor', 'toroidal_flux'}
     % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper:
-    % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293–302
+    % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293???302
     % since cocos=17 for LIUQE we get:
     % q = -dPhi/dpsi => Phi = - int(q*dpsi) which should always have the sign of B0
     % need to get q_rho but to avoid loop for rhotor in grids_1d, get q_rho explicitely here
@@ -2928,7 +2928,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'rhotor_edge', 'rhotor', 'rhotor_norm'}
     % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper:
-    % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293–302
+    % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293???302
     % since cocos=17 for LIUQE we get:
     % q = -dPhi/dpsi => Phi = - int(q*dpsi) which should always have the sign of B0
     % need to get q_rho but to avoid loop for rhotor in grids_1d, get q_rho explicitely here
@@ -3315,7 +3315,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         gdat_data.x = [1:size(sig,1)];
         gdat_data.dimunits = {'20 chords per camera'; 's'};
       else
-        % keyboard
+        keyboard
         % extract only given channels
         gdat_data.data = sig(channel_xtomo,:);
         gdat_data.x = channel_xtomo;
-- 
GitLab


From c52416df6c91d9f5008c92a2238930d574eb2dc7 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Sat, 15 Jan 2022 23:03:43 +0100
Subject: [PATCH 08/15] fix gdat_plot, mex for jet, and other small

---
 matlab/JET/gdat_jet.m                     |  8 ++++
 matlab/JET/rda_jet.m                      |  6 ++-
 matlab/TCV/gdat_tcv.m                     | 20 +++++----
 matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 55 +++++++++++++++++------
 matlab/gdat_plot.m                        |  2 +-
 5 files changed, 66 insertions(+), 25 deletions(-)

diff --git a/matlab/JET/gdat_jet.m b/matlab/JET/gdat_jet.m
index 2197e8f..4bb12dd 100644
--- a/matlab/JET/gdat_jet.m
+++ b/matlab/JET/gdat_jet.m
@@ -574,7 +574,15 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
         end
       else
         % assume same time for all
+        % mds_dir = '/home/amerle/public/mds-ssh';
+        mds_dir = '/home/amerle/public/mds-ssh3';
+        if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
+          addpath(mds_dir);
+        end
 	mdsconnect(['ssh://' gdat_params.jet_user '@mdsplus.jetdata.eu']);
+        if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
+          rmpath(mds_dir);
+        end
         rda_data_request = ['ppf/kk3/' num2str(i,'te%.2d')];
         bb = mdsvalue(['_bb=jet("' rda_data_request '",' num2str(shot) ')']);
         if ~isempty(bb) && numel(bb)==size(aa.data,2)
diff --git a/matlab/JET/rda_jet.m b/matlab/JET/rda_jet.m
index 1d93ae6..f80caa5 100644
--- a/matlab/JET/rda_jet.m
+++ b/matlab/JET/rda_jet.m
@@ -34,7 +34,8 @@ end
 if usemdsplus
   % use mdsplus
 
-  mds_dir = '/home/amerle/codes/mds-ssh';
+  % mds_dir = '/home/amerle/public/mds-ssh';
+  mds_dir = '/home/amerle/public/mds-ssh3';
   if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
     addpath(mds_dir);
   end
@@ -170,6 +171,9 @@ if usemdsplus
   if ~unix('test -d /home/duval/mdsplus')
     rmpath('/home/duval/mdsplus')
   end
+  if exist(mds_dir,'dir') && ~strcmp(which('mdsipmex'),fullfile(mds_dir,'mdsipmex.mexa64'))
+    rmpath(mds_dir);
+  end
 
 else
   % use RDA
diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 8b5bb3b..23935d1 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -348,11 +348,13 @@ if do_mdsopen_mdsclose
     end
     mapping_for_tcv.expression = mapping_for_tcv.expression(7:end);
   elseif shot==-1 || liuqe_version_eff==-1
-    ishot = mdsopen('pcs', shot);
+    [ishot,shot_stat] = mdsopen('pcs', shot);
   else
-    ishot = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
+    [ishot,shot_stat] = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
   end
-  if isempty(ishot) || ishot~=shot
+  if mod(shot_stat,2)==0
+    ishot
+    warning('error opening shot');
     if gdat_params.nverbose>=1; warning(['cannot open shot= ' num2str(shot)]); end
     return
   end
@@ -2501,7 +2503,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           if isempty(model_data.dim)
             fprintf('And it was NOT requested \n');
           else
-            fprintf('And it was requested. Check if plasma was a blip or if NB1 did not work \n');   
+            fprintf('And it was requested. Check if plasma was a blip or if NB1 did not work \n');
           end
         end
       end
@@ -2532,7 +2534,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
     %
     if any(strmatch('nbi2',gdat_data.gdat_params.source))
-      % NB2  
+      % NB2
       nodenameeff = '\results::NB2:POWR_TCV';
       nb2_data_tdi = tdi(nodenameeff);
       if ~isempty(nb2_data_tdi.dim) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
@@ -2565,14 +2567,14 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           if isempty(model_data.dim)
             fprintf('And it was NOT requested \n');
           else
-            fprintf('And it was requested. Check if plasma was a blip or if NB1 did not work \n');   
+            fprintf('And it was requested. Check if plasma was a blip or if NB1 did not work \n');
           end
         end
-      end        
+      end
     end
-    
+
     if any(strmatch('dnbi',gdat_data.gdat_params.source))
-    % NB2  
+    % NB2
     nodenameeff = '\RESULTS::DNBI:POWR_TCV';
       nb2_data_tdi = tdi(nodenameeff);
       if ~isempty(nb2_data_tdi.data) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index df66620..54f872b 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -107,6 +107,7 @@ global_quantities_desc.volume = params_eff.data_request;
 params_eff.data_request = 'w_mhd';
 global_quantities.w_mhd = gdat(params_equilibrium.shot,params_eff);
 global_quantities_desc.w_mhd = params_eff.data_request;
+ids_equilibrium_description.global_quantities = global_quantities_desc;
 
 global_quantities_fieldnames = fieldnames(global_quantities);
 special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
@@ -214,9 +215,10 @@ for it=1:ntime
     ids_equilibrium.time_slice{it}.boundary.x_point = {};
   end
 end
+ids_equilibrium_description.boundary = boundary_desc;
+ids_equilibrium_description.boundary_temp = temp_desc;
 
 %% constraints
-% TODO: Add description
 
 % Measured values
 liuqe_time = ids_equilibrium.time; % Not necessarily magnetics time so far
@@ -224,9 +226,15 @@ mag_time   = mdsvalue('\magnetics::bpol_003:dim0');
 itime = iround_os(mag_time, liuqe_time);
 mag_time_req = mdscvt(mag_time(itime),'f');
 bpol = mdsvalue('\magnetics::bpol_003[$1,*]',mag_time_req);
+nbpol = size(bpol,2);
+constraints_desc.bpol_probe.measured = sprintf('from %s i=1:%d','\magnetics::bpol_003[i,*]',nbpol);
 flux = mdsvalue('tcv_idealloop("flux")[$1,*]',mag_time_req);
-diam = mdsvalue('\results::dmlcor[$1]',mag_time_req);
+nflux = size(flux,2);
+constraints_desc.flux_loop.measured = sprintf('from %s i=1:%d','tcv_idealloop("flux")[i,*]',nflux);
 ip   = mdsvalue('\magnetics::iplasma:trapeze[$1]',mag_time_req);
+constraints_desc.ip.measured = sprintf('from %s','\magnetics::iplasma:trapeze');
+diam = mdsvalue('\results::dmlcor[$1]',mag_time_req);
+constraints_desc.diamagnetic_flux.measured = sprintf('from %s','\results::dmlcor');
 % Coil currents since dim of constraints pf_current is IDS:pf_active/coil
 dim_pol = {'OH_001','OH_002','OH_002','OH_002','OH_002','OH_002','OH_002',...
            'E_001','E_002','E_003','E_004','E_005','E_006','E_007','E_008',...
@@ -236,26 +244,38 @@ ipol = mdsvalue('\magnetics::ipol[$1,$2]',mag_time_req,dim_pol);
 ipol(:,27:29) = -ipol(:,27:29); % Bottom G-coils
 dim_pol(30:32) = {'TOR_001'};
 ipol(:,30:32) = 0; % TOR_001 is not used in LIUQE
+nipol = size(ipol,2);
+constraints_desc.pf_current.measured = sprintf('from %s, ii from %s','\magnetics::ipol[t,ii]',cell2str(dim_pol,15));
 
 % Reconstructed values
 ipol_liuqe_order = [18,19*ones(1,6),1:16,17*ones(1,6)]; % LIUQE order is E F G OH
 bpol_liuqe = mdsvalue('tcv_eq("b_probe","liuqe.m")');
+constraints_desc.bpol_probe.reconstructed = sprintf('%s','tcv_eq("b_probe","liuqe.m")');
 flux_liuqe = mdsvalue('tcv_eq("psi_loop","liuqe.m")');
-diam_liuqe = mdsvalue('tcv_eq("tor_flux_dml","liuqe.m")');
+constraints_desc.flux_liuqe.reconstructed = sprintf('%s','tcv_eq("psi_loop","liuqe.m")');
 ip_liuqe   = mdsvalue('tcv_eq("i_pl","liuqe.m")');
+constraints_desc.ip.reconstructed = sprintf('%s','tcv_eq("i_pl","liuqe.m")');
+diam_liuqe = mdsvalue('tcv_eq("tor_flux_dml","liuqe.m")');
+constraints_desc.diamagnetic_flux.reconstructed = sprintf('%s','tcv_eq("tor_flux_dml","liuqe.m")');
 ipol_liuqe = mdsvalue('tcv_eq("i_pol","liuqe.m")');
 ipol_liuqe = ipol_liuqe(ipol_liuqe_order,:);
 ipol_liuqe(27:29,:) = -ipol_liuqe(27:29,:); % Bottom G-coils
 ipol_liuqe(30:32,:) = 0; % ... TOR
+constraints_desc.ipol_liuqe.reconstructed = sprintf('%s','tcv_eq("i_pol","liuqe.m")');
 
 % Weights (using old parameters tree for now)
 bpol_err = mdsvalue('\results::parameters:berr')./mdsvalue('\results::parameters:vvv[0:37]');
+constraints_desc.bpol_probe.weight = sprintf('(%s/%s)^2','\results::parameters:vvv[0:37]','\results::parameters:berr');
 flux_err = mdsvalue('\results::parameters:ferr')./mdsvalue('\results::parameters:www[0:37]')*2*pi;
-diam_err = 0.13e-3./mdsvalue('\results::parameters:idml');
+constraints_desc.flux_loop.weight = sprintf('(%s/%s)^2','\results::parameters:www[0:37]','\results::parameters:ferr');
 ip_err   = mdsvalue('\results::parameters:plcerr')*1e3;
+constraints_desc.ip.weight = sprintf('(1e-3/%s)^2','\results::parameters:plcerr');
+diam_err = 0.13e-3./mdsvalue('\results::parameters:idml');
+constraints_desc.diamagnetic_flux.weight = sprintf('(%s/0.13e-3)^2','\results::parameters:idml');
 ipol_err = mdsvalue('\results::parameters:cerr')./mdsvalue('\results::parameters:uuu[0:18]')*1e3;
 ipol_err = ipol_err(ipol_liuqe_order);
 ipol_err(30:32) = NaN;
+constraints_desc.pf_current.weight = sprintf('(1e-3*%s/%s)^2','\results::parameters:uuu[0:18]','\results::parameters:cerr');
 
 if ntime > 0
   constraints_orig = ids_equilibrium.time_slice{1}.constraints;
@@ -263,15 +283,16 @@ if ntime > 0
   ununsed_constraints = {'faraday_angle','mse_polarisation_angle','iron_core_segment',...
                          'n_e','n_e_line','pressure','q','x_point'};
   for name = ununsed_constraints, constraints_orig.(name{1})={}; end
+  constraints_desc.ununsed_constraints = ununsed_constraints;
 end
+keyboard
 for it = 1:ntime
   constraints = constraints_orig;
   % bpol_probe
-  nbpol = size(bpol,2);
   bpol_probe(1:nbpol) = constraints.bpol_probe(1);
   for ib = 1:nbpol
     bpol_probe{ib}.measured = bpol(it,ib);
-    bpol_probe{ib}.source = sprintf('IDS:magnetics/bpol_probe[%02d]/field',ib);
+    bpol_probe{ib}.source = sprintf('IDS:equilibrium/../constraints/bpol_probe[%02d]',ib);
     bpol_probe{ib}.time_measurement = mag_time(itime(it));
     bpol_probe{ib}.exact = 0;
     bpol_probe{ib}.weight = 1/(bpol_err(ib)).^2;
@@ -279,11 +300,10 @@ for it = 1:ntime
   end
   constraints.bpol_probe = bpol_probe;
   % flux_loop
-  nflux = size(flux,2);
   flux_loop(1:nflux) = constraints.flux_loop(1);
   for il = 1:nflux
     flux_loop{il}.measured = flux(it,il);
-    flux_loop{il}.source = sprintf('IDS:magnetics/flux_loop[%02d]/flux',il);
+    flux_loop{il}.source = sprintf('IDS:equilibrium/../constraints/flux_loop[%02d]',il);
     flux_loop{il}.time_measurement = mag_time(itime(it));
     flux_loop{il}.exact = 0;
     flux_loop{il}.weight = 1/(flux_err(il)).^2;
@@ -292,24 +312,23 @@ for it = 1:ntime
   constraints.flux_loop = flux_loop;
   % ip
   constraints.ip.measured = ip(it);
-  constraints.ip.source = 'IDS:magnetics/method[1]/ip';
+  constraints.ip.source = 'IDS:equilibrium/../constraints/ip';
   constraints.ip.time_measurement = mag_time(itime(it));
   constraints.ip.exact = 0;
   constraints.ip.weight = 1/(ip_err).^2;
   constraints.ip.reconstructed = ip_liuqe(it);
   % diamagnetic_flux
   constraints.diamagnetic_flux.measured = diam(it);
-  constraints.diamagnetic_flux.source = 'IDS:magnetics/method[1]/diamagnetic_flux';
+  constraints.diamagnetic_flux.source = 'IDS:equilibrium/../constraints/diamagnetic_flux';
   constraints.diamagnetic_flux.time_measurement = mag_time(itime(it));
   constraints.diamagnetic_flux.exact = 0;
   constraints.diamagnetic_flux.weight = 1/(diam_err).^2;
   constraints.diamagnetic_flux.reconstructed = diam_liuqe(it);
   % pf_current
-  nipol = size(ipol,2);
   pf_current(1:nipol) = constraints.pf_current(1);
   for ic = 1:nipol
     pf_current{ic}.measured = ipol(it,ic);
-    pf_current{ic}.source = sprintf('IDS:pf_active/coil[%02d]/current',ic);
+    pf_current{ic}.source = sprintf('IDS:equilibrium/../constraints/pf_current[%02d]',ic);
     pf_current{ic}.time_measurement = mag_time(itime(it));
     if strcmp(dim_pol{ic},'TOR_001')
       pf_current{ic}.source = [pf_current{ic}.source,' replaced with 0'];
@@ -465,6 +484,8 @@ for it=1:ntime
   end
 end
 
+ids_equilibrium_description.profiles_1d = profiles_1d_desc;
+
 % special cases
 for it=1:ntime
   ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.b_field_tor = ids_equilibrium.time_slice{it}.profiles_1d.f(1) ...
@@ -496,9 +517,13 @@ profiles_2d.grid_type.description = 'Cylindrical R,Z ala eqdsk';
 params_eff.data_request = 'psi';
 profiles_2d.psi = gdat(params_equilibrium.shot,params_eff); % add psi_bound in a second step in special cases
 profiles_2d_desc.psi = [params_eff.data_request ' adding psi_bound in a 2nd step'];
-% r = gdat(params_equilibrium.shot,'r','machine',gdat_params.machine); % not to be filled since in grid.dim1
+%profiles_2d.r = profiles_2d.psi;
+profiles_2d.r.data = repmat(repmat(profiles_2d.psi.dim{1},1,65),1,1,1299);
+profiles_2d_desc.r = 'from dim{1} of ''psi'' repeated';
+profiles_2d.z.data = repmat(repmat(profiles_2d.psi.dim{2}',28,1),1,1,1299);
+profiles_2d_desc.z = 'from dim{2} of ''psi'' repeated';
+
 % theta = gdat(params_equilibrium.shot,'theta','machine',gdat_params.machine);
-% z = gdat(params_equilibrium.shot,'z','machine',gdat_params.machine); % not to be filled since in grid.dim2
 
 profiles_2d_fieldnames = fieldnames(profiles_2d);
 special_fields = {'grid', 'grid_type'}; % fields needing non-automatic treatments
@@ -529,6 +554,8 @@ for it=1:ntime
       global_quantities.psi_boundary.data(it);
 end
 
+ids_equilibrium_description.profiles_2d = profiles_2d_desc;
+
 % make arrays not filled in empty:
 ids_equilibrium.grids_ggd = {};
 for it=1:numel(ids_equilibrium.time_slice)
diff --git a/matlab/gdat_plot.m b/matlab/gdat_plot.m
index 3334870..aba5eac 100644
--- a/matlab/gdat_plot.m
+++ b/matlab/gdat_plot.m
@@ -124,7 +124,7 @@ if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(
     zoom on;
   end
   maxnblines = 1;
-  if strcmp(gdat_data.gdat_params.data_request,'powers') && length(ab)==length(gdat_data.label)
+  if any(strcmp(gdat_data.gdat_params.data_request,'powers')) && length(ab)==length(gdat_data.label)
     % keep legend from label
   else
     for i=1:numel(linehandles)
-- 
GitLab


From 47396e7586a4e0116b72b210877e215436a4ed7c Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Sun, 16 Jan 2022 09:00:35 +0100
Subject: [PATCH 09/15] forgot an mdsopen call

---
 matlab/TCV/gdat_tcv.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 23935d1..22dbcbf 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -342,7 +342,7 @@ if do_mdsopen_mdsclose
   %%%  if liuqe_version_eff==-1
   if ~iscell(mapping_for_tcv.expression) && any(strfind(mapping_for_tcv.expression,'\rtc::'))
     mdsdisconnect;
-    ishot = mdsopen('rtc',shot); % sub-tree
+    [ishot,shot_stat]  = mdsopen('rtc',shot); % sub-tree
     if any(strfind(data_request_eff,'\rtc::'))
       data_request_eff = data_request_eff(7:end);
     end
-- 
GitLab


From ceadcebd8d18d16a46db557a6956b344ebd2349d Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 15:22:09 +0100
Subject: [PATCH 10/15] add gas to AUG, add Tags to gdat_plot

---
 matlab/AUG/aug_requests_mapping.m |  3 ++
 matlab/AUG/gdat_aug.m             | 80 +++++++++++++++++++++++++++++++
 matlab/gdat_plot.m                | 40 ++++++++++++++--
 3 files changed, 118 insertions(+), 5 deletions(-)

diff --git a/matlab/AUG/aug_requests_mapping.m b/matlab/AUG/aug_requests_mapping.m
index 9b4e16a..4112db7 100644
--- a/matlab/AUG/aug_requests_mapping.m
+++ b/matlab/AUG/aug_requests_mapping.m
@@ -130,6 +130,9 @@ switch lower(data_request)
   mapping.gdat_timedim = 2;
   mapping.method = 'switchcase'; % could use function make_eqdsk directly?
   mapping.expression = '';
+ case {'gas', 'gas_valve'}
+  mapping.gdat_timedim = 2;
+  mapping.method = 'switchcase';
  case 'halpha'
   mapping.timedim = 1;
   mapping.label = 'Halpha';
diff --git a/matlab/AUG/gdat_aug.m b/matlab/AUG/gdat_aug.m
index 3cc7d98..58d292a 100644
--- a/matlab/AUG/gdat_aug.m
+++ b/matlab/AUG/gdat_aug.m
@@ -1297,6 +1297,86 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
     gdat_data.dimunits{1} = 'rhopolnorm';
     gdat_data.dimunits{2} = 'time [s]';
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+   case {'gas', 'gas_valve', 'gas_valves'}
+    %
+    uvs_signals_source = [{'D_tot','Sum of calibrated D2 gas-fluxes (= D/s)'}; ...
+                    {'N_tot','Sum of calibrated N2 gas-fluxes (= N*7/s)'}; ...
+                    {'H_tot','Sum of calibrated H gas-fluxes (= H/s)'}; ...
+                    {'He_tot','Sum of calibrated He gas-fluxes (= He*2/s)'}; ...
+                    {'Ne_tot','Sum of calibrated Ne gas-fluxes (= Ne*10/s)'}; ...
+                    {'Ar_tot','Sum of calibrated Ar gas-fluxes (= Ar*18/s)'}];
+    for i=1:size(uvs_signals_source,1)
+      uvs_signals.(uvs_signals_source{i,1}) = gdat_aug(shot,{'UVS',uvs_signals_source{i,1}});
+      if isempty(gdat_data.t) && prod(isfinite(uvs_signals.(uvs_signals_source{i,1}).t)) && isvector(uvs_signals.(uvs_signals_source{i,1}).t) ...
+          && numel(uvs_signals.(uvs_signals_source{i,1}).t)>10
+        gdat_data.t = uvs_signals.(uvs_signals_source{i,1}).t; % all have same time base, could use time signal
+      end
+      if prod(isfinite(uvs_signals.(uvs_signals_source{i,1}).data)) && isvector(uvs_signals.(uvs_signals_source{i,1}).data) ...
+          && numel(uvs_signals.(uvs_signals_source{i,1}).data) == numel(gdat_data.t);
+        gdat_data.data(i,:) = uvs_signals.(uvs_signals_source{i,1}).data; % all have same time base
+      else
+        gdat_data.data(i,:) = NaN;
+      end
+      gdat_data.units{i} = uvs_signals_source{i,2};
+    end
+    gdat_data.data_fullpath = 'from UVS and uvs_signals_source';
+    gdat_data.uvs_signals_source = uvs_signals_source;
+    gdat_data.label = uvs_signals_source(:,1)';
+    gdat_data.dim{2} = gdat_data.t;
+    gdat_data.dimunits{2} = 's';
+    gdat_data.dim{1} = [1:size(gdat_data.data,1)];
+    gdat_data.dimunits{1} = 'uvs_signals_source index';
+    gdat_data.x = gdat_data.dim{1};
+    gdat_data.mapping_for.aug.timedim = 1;
+
+    % note also available in UVS, not commented the same way on ISIS
+    uvd_signals_source=[{'CFCu01T','PV01','empty comment'}; ...
+                        {'CFCo02A','PV02','valve flux, normally for prefill'}; ...
+                        {'CFFo02A','PV03','valve flux,'}; ...
+                        {'CFA13B','PV04','valve flux, in-vessel pipe for gas imaging'}; ...
+                        {'CFA03B','PV05','valve flux,'}; ...
+                        {'CFFo07A','PV06','valve flux,'}; ...
+                        {'CFDu05X','PV07','valve flux, normally with N2'}; ...
+                        {'CFDu05B','PV08','valve flux, normally with D2'}; ...
+                        {'CFCo10A','PV09','valve flux, in-vessel pipe into ICRH limiter'}; ...
+                        {'CFDu01B','PV10','valve flux, normally with D2'}; ...
+                        {'CFDu09B','PV11','valve flux, normally with D2'}; ...
+                        {'CFFo10A','PV12','valve flux,'}; ...
+                        {'CFFo14A','PV13','valve flux,'}; ...
+                        {'CFA03C','PV14','valve flux, short feeding pipe for rare gases'}; ...
+                        {'CFA13A','PV15','valve flux,'}; ...
+                        {'CFA03A','PV16','valve flux,'}; ...
+                        {'CFDu09X','PV17','valve flux, normally with N2'}; ...
+                        {'CFDu13X','PV18','valve flux, normally with N2'}; ...
+                        {'CFDu13B','PV19','valve flux, normally with D2'}; ...
+                        {'CFDu01X','PV20','valve flux, normally with N2'}];
+    gdat_data.pvxx.t = [];
+    gdat_data.pvxx.t = [];
+
+    for i=1:size(uvd_signals_source,1)
+      uvd_signals.(uvd_signals_source{i,2}) = gdat_aug(shot,{'UVD',uvd_signals_source{i,1}});
+      if isempty(gdat_data.pvxx.t) && prod(isfinite(uvd_signals.(uvd_signals_source{i,2}).t)) && isvector(uvd_signals.(uvd_signals_source{i,2}).t) ...
+          && numel(uvd_signals.(uvd_signals_source{i,2}).t)>10
+        gdat_data.pvxx.t = uvd_signals.(uvd_signals_source{i,2}).t; % all have same time base
+      end
+      if prod(isfinite(uvd_signals.(uvd_signals_source{i,2}).data)) && isvector(uvd_signals.(uvd_signals_source{i,2}).data) ...
+          && numel(uvd_signals.(uvd_signals_source{i,2}).data) == numel(gdat_data.pvxx.t);
+        gdat_data.pvxx.data(i,:) = uvd_signals.(uvd_signals_source{i,2}).data; % all have same time base
+      else
+        gdat_data.pvxx.data(i,:) = NaN;
+      end
+      gdat_data.pvxx.units{i} = [uvd_signals_source{i,2} ': ' uvd_signals_source{i,3}];
+    end
+    gdat_data.pvxx.data_fullpath = 'from UVD and uvd_signals_source';
+    gdat_data.pvxx.uvd_signals_source = uvd_signals_source;
+    gdat_data.pvxx.label = uvd_signals_source(:,2)';
+    gdat_data.pvxx.dim{2} = gdat_data.t;
+    gdat_data.pvxx.dimunits{2} = 's';
+    gdat_data.pvxx.dim{1} = [1:size(gdat_data.pvxx.data,1)];
+    gdat_data.pvxx.dimunits{1} = 'PVxx index';
+    gdat_data.pvxx.x = gdat_data.dim{1};
+
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ids'}
     params_eff = gdat_data.gdat_params;
diff --git a/matlab/gdat_plot.m b/matlab/gdat_plot.m
index aba5eac..acd7768 100644
--- a/matlab/gdat_plot.m
+++ b/matlab/gdat_plot.m
@@ -37,6 +37,7 @@ elseif isfield(gdat_data,'gdat_params') && isfield(gdat_data.gdat_params,'doplot
   end
 end
 if doplot==0; return; end
+redo_legend_from_Tags = 0;
 
 if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(gdat_data.t)
   fighandle = get(0,'CurrentFigure');
@@ -106,25 +107,54 @@ if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(
     else
       hylabel=ylabel(ylabel_eff,'interpreter','none');
     end
-    if iscell(gdat_data.label) && length(gdat_data.label)>=2;
-      ab=get(gca,'children');
+    ab=get(gca,'children');
+    if iscell(gdat_data.label) && numel(gdat_data.label) > 1;
+      if numel(ab)==numel(gdat_data.label) || mod(numel(ab),numel(gdat_data.label))==0
+        % Assumes a single shot with several lines and labels, if mod==0 then this is a new shot, only Tag present lines
+        for iab=1:numel(gdat_data.label)
+          set(ab(iab),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label{end-iab+1}]);
+        end
+        if numel(ab)>numel(gdat_data.label), redo_legend_from_Tags = 1; end
+      end
       if iscell(gdat_data.label) && length(ab) < length(gdat_data.label)
         if length(ab) == 1
           % assume combine label is best
           tempaaa = sprintf('%s/',gdat_data.label{:});
           hhhleg=legend(tempaaa(1:end-1));
+          set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' tempaaa(1:min(10,numel(tempaaa)-1))]);
         else
+          % not sure why would get there
           hhhleg=legend(gdat_data.label{1:length(ab)});
         end
-      else
+      elseif numel(ab)==numel(gdat_data.label)
         hhhleg=legend(gdat_data.label);
       end
-      set(hhhleg,'Interpreter','none');
+      if exist('hhhleg'), set(hhhleg,'Interpreter','none'); end
+    else
+      % assume one line per call
+      if isempty(gdat_data.label)
+        set(ab(1),'Tag',[num2str(gdat_data.shot)]);
+      elseif ischar(gdat_data.label)
+        % assume one signal, take max 10 chars
+        set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label(1:min(10,numel(gdat_data.label)))]);
+      elseif iscell(gdat_data.label) && numel(gdat_data.label) == 1
+        if isempty(gdat_data.label{1})
+          set(ab(1),'Tag',[num2str(gdat_data.shot)]);
+        else
+          set(ab(1),'Tag',[num2str(gdat_data.shot) ' ' gdat_data.label{1}(1:min(10,numel(gdat_data.label{1})))]);
+        end
+      end
+    end
+    if redo_legend_from_Tags
+      for iab=1:numel(ab)
+        % Use Tag for DisplayName, when overlay plots of multiple data at this stage
+        set(ab(iab),'DisplayName',strrep(get(ab(iab),'Tag'),'_','\_'));
+      end
     end
     zoom on;
   end
   maxnblines = 1;
-  if any(strcmp(gdat_data.gdat_params.data_request,'powers')) && length(ab)==length(gdat_data.label)
+  if redo_legend_from_Tags || any(strcmp(gdat_data.gdat_params.data_request,'powers')) || (numel(ab)==numel(gdat_data.label) && numel(ab)>1)
     % keep legend from label
   else
     for i=1:numel(linehandles)
-- 
GitLab


From 34e583439b66725af9b19158f3506ebd3569533a Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 16:41:14 +0100
Subject: [PATCH 11/15] include changes for NBI from Matteo to ease rebase

---
 matlab/TCV/gdat_tcv.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 22dbcbf..ea437b8 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1,4 +1,4 @@
-function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_request,varargin)
+function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv_OS(shot,data_request,varargin)
 %
 % function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
 %
@@ -3315,7 +3315,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         gdat_data.x = [1:size(sig,1)];
         gdat_data.dimunits = {'20 chords per camera'; 's'};
       else
-        keyboard
+        % keyboard
         % extract only given channels
         gdat_data.data = sig(channel_xtomo,:);
         gdat_data.x = channel_xtomo;
-- 
GitLab


From d9fdb2aaf0822141e4d3a23069c1ae807c5de937 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 17:19:15 +0100
Subject: [PATCH 12/15] try to finalize tcv_get_ids_equilibrium changes

to be tested
---
 matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index 54f872b..10dcc0c 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -340,7 +340,7 @@ for it = 1:ntime
     end
   end
   constraints.pf_current = pf_current;
-  
+
   ids_equilibrium.time_slice{it}.constraints = constraints;
 end
 
@@ -518,9 +518,9 @@ params_eff.data_request = 'psi';
 profiles_2d.psi = gdat(params_equilibrium.shot,params_eff); % add psi_bound in a second step in special cases
 profiles_2d_desc.psi = [params_eff.data_request ' adding psi_bound in a 2nd step'];
 %profiles_2d.r = profiles_2d.psi;
-profiles_2d.r.data = repmat(repmat(profiles_2d.psi.dim{1},1,65),1,1,1299);
+profiles_2d.r.data = repmat(repmat(profiles_2d.psi.dim{1},1,numel(profiles_2d.psi.dim{2})),1,1,numel(profiles_2d.psi.dim{3}));
 profiles_2d_desc.r = 'from dim{1} of ''psi'' repeated';
-profiles_2d.z.data = repmat(repmat(profiles_2d.psi.dim{2}',28,1),1,1,1299);
+profiles_2d.z.data = repmat(repmat(profiles_2d.psi.dim{2}',numel(profiles_2d.psi.dim{1}),1),1,1,numel(profiles_2d.psi.dim{3}));
 profiles_2d_desc.z = 'from dim{2} of ''psi'' repeated';
 
 % theta = gdat(params_equilibrium.shot,'theta','machine',gdat_params.machine);
-- 
GitLab


From 48034ef913094b1e585c95a98bbcb00a1870c78a Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 17:43:35 +0100
Subject: [PATCH 13/15] add constraints description

tested on TCV
---
 matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index 10dcc0c..2638760 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -285,7 +285,6 @@ if ntime > 0
   for name = ununsed_constraints, constraints_orig.(name{1})={}; end
   constraints_desc.ununsed_constraints = ununsed_constraints;
 end
-keyboard
 for it = 1:ntime
   constraints = constraints_orig;
   % bpol_probe
@@ -344,6 +343,7 @@ for it = 1:ntime
   ids_equilibrium.time_slice{it}.constraints = constraints;
 end
 
+ids_equilibrium_description.constraints = constraints_desc;
 
 %
 %% profiles_1d (cannot use eqdsk since not same radial mesh)
-- 
GitLab


From a58c6ef127baf3c0c891ef3fca9d258c0a9253bb Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 17:50:33 +0100
Subject: [PATCH 14/15] temp

---
 matlab/TCV/gdat_tcv.m | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index ea437b8..b2a647e 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1,4 +1,4 @@
-function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv_OS(shot,data_request,varargin)
+function [gdat_data,gdat_params,error_status,varargout] = gdat_tcv(shot,data_request,varargin)
 %
 % function [gdat_data,gdat_params,error_status,varargout] = gdat(shot,data_request,varargin)
 %
@@ -2574,8 +2574,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
 
     if any(strmatch('dnbi',gdat_data.gdat_params.source))
-    % NB2
-    nodenameeff = '\RESULTS::DNBI:POWR_TCV';
+      % NB2
+      nodenameeff = '\RESULTS::DNBI:POWR_TCV';
       nb2_data_tdi = tdi(nodenameeff);
       if ~isempty(nb2_data_tdi.data) && ~ischar(nb2_data_tdi.data) && ~isempty(nb2_data_tdi.dim)
         nbi_neutral_power_tot = nb2_data_tdi.data.*1e6; % in W
@@ -3315,7 +3315,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         gdat_data.x = [1:size(sig,1)];
         gdat_data.dimunits = {'20 chords per camera'; 's'};
       else
-        % keyboard
+        keyboard
         % extract only given channels
         gdat_data.data = sig(channel_xtomo,:);
         gdat_data.x = channel_xtomo;
-- 
GitLab


From f5bc4eb7d98a2cba68c6de5e05283f87e4f8f34f Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 17 Jan 2022 19:21:28 +0100
Subject: [PATCH 15/15] simplifications

---
 matlab/AUG/gdat_aug.m | 12 ++++--------
 matlab/TCV/gdat_tcv.m |  4 +---
 matlab/gdat_plot.m    |  2 +-
 3 files changed, 6 insertions(+), 12 deletions(-)

diff --git a/matlab/AUG/gdat_aug.m b/matlab/AUG/gdat_aug.m
index 58d292a..93b3ba7 100644
--- a/matlab/AUG/gdat_aug.m
+++ b/matlab/AUG/gdat_aug.m
@@ -1308,12 +1308,10 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
                     {'Ar_tot','Sum of calibrated Ar gas-fluxes (= Ar*18/s)'}];
     for i=1:size(uvs_signals_source,1)
       uvs_signals.(uvs_signals_source{i,1}) = gdat_aug(shot,{'UVS',uvs_signals_source{i,1}});
-      if isempty(gdat_data.t) && prod(isfinite(uvs_signals.(uvs_signals_source{i,1}).t)) && isvector(uvs_signals.(uvs_signals_source{i,1}).t) ...
-          && numel(uvs_signals.(uvs_signals_source{i,1}).t)>10
+      if isempty(gdat_data.t) && numel(uvs_signals.(uvs_signals_source{i,1}).t)>10
         gdat_data.t = uvs_signals.(uvs_signals_source{i,1}).t; % all have same time base, could use time signal
       end
-      if prod(isfinite(uvs_signals.(uvs_signals_source{i,1}).data)) && isvector(uvs_signals.(uvs_signals_source{i,1}).data) ...
-          && numel(uvs_signals.(uvs_signals_source{i,1}).data) == numel(gdat_data.t);
+      if ~isempty(uvs_signals.(uvs_signals_source{i,1}).data) && numel(uvs_signals.(uvs_signals_source{i,1}).data) == numel(gdat_data.t);
         gdat_data.data(i,:) = uvs_signals.(uvs_signals_source{i,1}).data; % all have same time base
       else
         gdat_data.data(i,:) = NaN;
@@ -1356,12 +1354,10 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
 
     for i=1:size(uvd_signals_source,1)
       uvd_signals.(uvd_signals_source{i,2}) = gdat_aug(shot,{'UVD',uvd_signals_source{i,1}});
-      if isempty(gdat_data.pvxx.t) && prod(isfinite(uvd_signals.(uvd_signals_source{i,2}).t)) && isvector(uvd_signals.(uvd_signals_source{i,2}).t) ...
-          && numel(uvd_signals.(uvd_signals_source{i,2}).t)>10
+      if isempty(gdat_data.pvxx.t) && numel(uvd_signals.(uvd_signals_source{i,2}).t)>10
         gdat_data.pvxx.t = uvd_signals.(uvd_signals_source{i,2}).t; % all have same time base
       end
-      if prod(isfinite(uvd_signals.(uvd_signals_source{i,2}).data)) && isvector(uvd_signals.(uvd_signals_source{i,2}).data) ...
-          && numel(uvd_signals.(uvd_signals_source{i,2}).data) == numel(gdat_data.pvxx.t);
+      if ~isempty(uvd_signals.(uvd_signals_source{i,2}).data) && numel(uvd_signals.(uvd_signals_source{i,2}).data) == numel(gdat_data.pvxx.t);
         gdat_data.pvxx.data(i,:) = uvd_signals.(uvd_signals_source{i,2}).data; % all have same time base
       else
         gdat_data.pvxx.data(i,:) = NaN;
diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index b2a647e..f84e5af 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -353,9 +353,7 @@ if do_mdsopen_mdsclose
     [ishot,shot_stat] = mdsopen(shot); % if ishot equal to shot, then mdsclose at the end
   end
   if mod(shot_stat,2)==0
-    ishot
-    warning('error opening shot');
-    if gdat_params.nverbose>=1; warning(['cannot open shot= ' num2str(shot)]); end
+    if gdat_params.nverbose>=1; warning('error opening shot %d. Value returned:%d',shot,ishot); end
     return
   end
 end
diff --git a/matlab/gdat_plot.m b/matlab/gdat_plot.m
index acd7768..d2a1ede 100644
--- a/matlab/gdat_plot.m
+++ b/matlab/gdat_plot.m
@@ -129,7 +129,7 @@ if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(
       elseif numel(ab)==numel(gdat_data.label)
         hhhleg=legend(gdat_data.label);
       end
-      if exist('hhhleg'), set(hhhleg,'Interpreter','none'); end
+      if exist('hhhleg','var'), set(hhhleg,'Interpreter','none'); end
     else
       % assume one line per call
       if isempty(gdat_data.label)
-- 
GitLab