From fb0df26afc70eb8d6c057f87ef1e75ea5af9f7a1 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Mon, 13 Apr 2020 15:50:53 +0200
Subject: [PATCH 1/4] start adding tcv_get_ids_summary, add global
 tokamak_agnostic ids_summary_disruption.m at top

to be modified by Alessandro
---
 matlab/TCV_IMAS/tcv_get_ids_summary.m | 664 ++++++++++++++++++++++++++
 matlab/TCV_IMAS/tcv_ids_headpart.m    |  28 +-
 matlab/ids_summary_disruption.m       |  30 ++
 3 files changed, 702 insertions(+), 20 deletions(-)
 create mode 100644 matlab/TCV_IMAS/tcv_get_ids_summary.m
 create mode 100644 matlab/ids_summary_disruption.m

diff --git a/matlab/TCV_IMAS/tcv_get_ids_summary.m b/matlab/TCV_IMAS/tcv_get_ids_summary.m
new file mode 100644
index 00000000..b804946f
--- /dev/null
+++ b/matlab/TCV_IMAS/tcv_get_ids_summary.m
@@ -0,0 +1,664 @@
+function [ids_summary,ids_summary_description,varargout] = tcv_get_ids_summary(shot,ids_equil_empty, gdat_params,varargin)
+%
+%  [ids_summary,ids_summary_description,varargout] = tcv_get_ids_summary(shot,ids_equil_empty,varargin);
+%
+%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar and cocos_out options
+%
+
+if exist('gdat_params','var')
+  [ids_summary, params_summary] = tcv_ids_headpart(shot,ids_equil_empty,'summary','gdat_params',gdat_params,varargin{:});
+else
+  [ids_summary, params_summary] = tcv_ids_headpart(shot,ids_equil_empty,'summary',varargin{:});
+  aa=gdat_tcv;
+  gdat_params = aa.gdat_params; % to get default params
+end
+
+% As a general rule, for a new substructure under the main ids, construct a local structure like:
+% "global_quantities" with subfields being the relevant data to get and a local structure:
+% "global_quantities_desc" which contains the same subfields themselves containing the gdat string aftre shot used
+%
+% vacuum_toroidal_field and time, using homogeneous
+%
+%% liuqe.m at this stage is missing correction 0.996, so use std source by time of default liuqe to make sure
+params_eff_ref = gdat_params; params_eff_ref.doplot=0;
+try params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def)
+
+% code description
+[abcpath,abcfile,abcext]=fileparts(mfilename('fullpath'));
+[aa1,aa2]=unix(['cd ' abcpath '; git describe --tags --first-parent --abbrev=11 --long --dirty --always']);
+git_release_hash='gdat git hash not found';
+if aa1==0; git_release_hash = aa2(1:end-1); end % avoid newline
+[aa1,aa2]=unix(['cd ' abcpath '; git config --get remote.origin.url']);
+git_url = 'gdat git url not found';
+if aa1==0; git_url = aa2(1:end-1); end
+ids_summary.code.name = 'gdat';
+ids_summary.code.version = git_release_hash;
+ids_summary.code.repository = git_url;
+
+% use equilibrium time as reference, so load first elongation(time)
+
+params_eff = params_eff_ref;
+params_eff.data_request='kappa'; params_eff.source='liuqe'; % to get liuqe time array
+kappa = gdat(params_summary.shot,params_eff);
+ids_summary.time = kappa.t;
+ids_summary_description.time = 'time from gdat(shot,''kappa'')';
+
+params_eff = params_eff_ref;
+% boundary:
+boundary.elongation = kappa;
+boundary_desc.elongation = kappa.gdat_params.data_request;
+params_eff.data_request = 'r_geom';
+boundary.geometric_axis_r = gdat(params_summary.shot,params_eff);
+boundary_desc.geometric_axis_r = params_eff.data_request;
+params_eff.data_request = 'z_geom';
+boundary.geometric_axis_z = gdat(params_summary.shot,params_eff);
+boundary_desc.geometric_axis_z = params_eff.data_request;
+params_eff.data_request = 'r_mag';
+boundary.magnetic_axis_r = gdat(params_summary.shot,params_eff);
+boundary_desc.magnetic_axis_r = params_eff.data_request;
+params_eff.data_request = 'z_mag';
+boundary.magnetic_axis_z = gdat(params_summary.shot,params_eff);
+boundary_desc.magnetic_axis_z = params_eff.data_request;
+params_eff.data_request = 'a_minor';
+boundary.minor_radius = gdat(params_summary.shot,params_eff);
+boundary_desc.minor_radius = params_eff.data_request;
+params_eff.data_request = 'delta_bottom';
+boundary.triangularity_lower = gdat(params_summary.shot,params_eff);
+boundary_desc.triangularity_lower = params_eff.data_request;
+params_eff.data_request = 'delta_top';
+boundary.triangularity_upper = gdat(params_summary.shot,params_eff);
+boundary_desc.triangularity_upper = params_eff.data_request;
+% $$$ params_eff.data_request = 'tcv_eq(''''n_xpts'''',''''liuqe.m'''')';
+% $$$ temp.n_x_point = gdat(params_summary.shot,params_eff);
+% $$$ temp_desc.n_x_point = params_eff.data_request;
+% $$$ params_eff.data_request = 'tcv_eq(''''r_xpts'''',''''liuqe.m'''')';
+% $$$ temp.r_x_point = gdat(params_summary.shot,params_eff);
+% $$$ temp_desc.r_x_point = params_eff.data_request;
+% $$$ params_eff.data_request = 'tcv_eq(''''z_xpts'''',''''liuqe.m'''')';
+% $$$ temp.z_x_point = gdat(params_summary.shot,params_eff);
+% $$$ temp_desc.z_x_point = params_eff.data_request;
+
+boundary_fieldnames = fieldnames(boundary);
+special_fields = {}; % fields needing non-automatic treatments
+% assume all boundary have same time base since called with same gdat_params for liuqe
+for i=1:numel(boundary_fieldnames)
+  if ~any(strcmp(boundary_fieldnames{i},special_fields))
+    if ~isstruct(ids_summary.boundary.(boundary_fieldnames{i}).value)
+      ids_summary.boundary.(boundary_fieldnames{i}).value = boundary.(boundary_fieldnames{i}).data;
+      ids_summary.boundary.(boundary_fieldnames{i}).source = ['gdat request: ' boundary_desc.(boundary_fieldnames{i})];
+    else
+      special_fields{end+1} = boundary_fieldnames{i};
+    end
+  end
+end
+
+% special case
+if ~isempty(special_fields)
+% $$$   for it=1:ntime
+% $$$     ids_summary.time_slice{it}.boundary.magnetic_axis.r = temp.r_magnetic_axis.data(it);
+% $$$     ids_summary.time_slice{it}.boundary.magnetic_axis.z = temp.z_magnetic_axis.data(it);
+% $$$     ids_summary.time_slice{it}.boundary.psi_axis = temp.psi_axis.data(it) + ...
+% $$$         ids_summary.time_slice{it}.boundary.psi_boundary;
+% $$$     [zz,izz]=min(temp.q_rho.data(:,it));
+% $$$     ids_summary.time_slice{it}.boundary.q_min.value = zz;
+% $$$     ids_summary.time_slice{it}.boundary.q_min.rho_tor_norm = temp.q_rho.grids_1d.rhotornorm(izz);
+% $$$   end
+end
+
+% disruption, use specific function
+% to be decided if changes only ids_summary_out.disruption or also elsewhere..., if ids_summary_out.time is changed,
+% then all the other data should be related, thus call this first?
+
+[ids_summary_out,ids_summary_out_description] = ids_summary_disruption(shot, ids_summary, gdat_params);
+
+return
+% to continue later...
+
+params_eff = params_eff_ref;
+params_eff.data_request='b0'; params_eff.source='liuqe'; % to get liuqe time array
+bb = gdat(params_summary.shot,params_eff);
+params_eff = rmfield(params_eff,'source'); % to get original magnetics data
+vacuum_toroidal_field.b0=gdat(params_summary.shot,params_eff);
+ij_ok = isfinite(vacuum_toroidal_field.b0.data);
+vacuum_toroidal_field.b0.data = interpos(21,vacuum_toroidal_field.b0.t(ij_ok),vacuum_toroidal_field.b0.data(ij_ok),bb.t);
+vacuum_toroidal_field.b0.t = bb.t;
+vacuum_toroidal_field.b0.dim = {vacuum_toroidal_field.b0.t};
+vacuum_toroidal_field_desc.b0 = ['''b0'',''timing source'',''liuqe=' num2str(gdat_params.liuqe) ''''];
+vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'',''source'',''liuqe'']';
+ids_summary.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
+ids_summary.vacuum_toroidal_field.b0 = vacuum_toroidal_field.b0.data;
+ids_summary_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
+
+ids_summary.time = vacuum_toroidal_field.b0.t;
+ids_summary_description.time = '.t subfield from: [''b0'',''source'',''liuqe'']';
+ntime = numel(ids_summary.time);
+
+ids_summary.time_slice(1:ntime) = ids_summary.time_slice(1);
+
+% load time array data to copy to time_slices
+
+% global_quantities data into local global_quantities.* structure with correct end names and global_quantities_desc.* with description. Use temp.* and temp_desc.* structures for temporary data
+
+% brute force solution load all eqdsks
+% $$$ for it=1:ntime
+% $$$   ids_summary.time(it)
+% $$$   temp.eqdsks{it}=gdat(params_summary.shot,'eqdsk','time',ids_summary.time(it),'write',0,'machine',gdat_params.machine);
+% $$$ end
+% $$$ temp_desc.eqdsks{1} = '''eqdsk'',''time'',ids_summary.time(it)';
+
+params_eff = params_eff_ref;
+params_eff.data_request = 'area_edge';
+global_quantities.area = gdat(params_summary.shot,params_eff);
+global_quantities_desc.area = params_eff.data_request;
+params_eff.data_request = 'betan';
+global_quantities.beta_normal = gdat(params_summary.shot,params_eff);
+global_quantities_desc.beta_normal = params_eff.data_request;
+params_eff.data_request = 'betap';
+global_quantities.beta_pol = gdat(params_summary.shot,params_eff);
+global_quantities_desc.beta_pol = params_eff.data_request;
+params_eff.data_request = 'beta';
+global_quantities.beta_tor = gdat(params_summary.shot,params_eff);
+global_quantities_desc.beta_tor = params_eff.data_request;
+params_eff.data_request = 'w_mhd';
+global_quantities.energy_mhd = gdat(params_summary.shot,params_eff);
+global_quantities_desc.energy_mhd = params_eff.data_request;
+params_eff.data_request = 'ip';
+global_quantities.ip = gdat(params_summary.shot,params_eff);
+global_quantities_desc.ip = params_eff.data_request;
+% length_pol = gdat(params_summary.shot,'length_pol','machine',gdat_params.machine); % to be added
+params_eff.data_request = 'li';
+global_quantities.li_3 = gdat(params_summary.shot,params_eff);
+global_quantities_desc.li_3 = params_eff.data_request;
+params_eff.data_request = 'r_axis';
+temp.r_magnetic_axis = gdat(params_summary.shot,params_eff);
+temp_desc.r_magnetic_axis = params_eff.data_request;
+params_eff.data_request = 'z_axis';
+temp.z_magnetic_axis = gdat(params_summary.shot,params_eff);
+temp_desc.z_magnetic_axis = params_eff.data_request;
+params_eff.data_request = 'psi_axis';
+temp.psi_axis = gdat(params_summary.shot,params_eff); % needs to add psi_edge sincepsi_axis liuqe assuming 0 dege value
+temp_desc.psi_axis = params_eff.data_request;
+params_eff.data_request = 'psi_edge';
+global_quantities.psi_boundary = gdat(params_summary.shot,params_eff);
+global_quantities_desc.psi_boundary = params_eff.data_request;
+params_eff.data_request = 'q95';
+global_quantities.q_95 = gdat(params_summary.shot,params_eff);
+global_quantities_desc.q_95 = params_eff.data_request;
+params_eff.data_request = 'q0';
+global_quantities.q_axis = gdat(params_summary.shot,params_eff); % will be checked with q_rho?
+global_quantities_desc.q_axis = params_eff.data_request;
+params_eff.data_request = 'q_rho';
+temp.q_rho = gdat(params_summary.shot,params_eff);
+temp_desc.q_rho = params_eff.data_request;
+% surface = gdat(params_summary.shot,'surface','machine',gdat_params.machine); % to be added
+params_eff.data_request = 'volume';
+global_quantities.volume = gdat(params_summary.shot,params_eff);
+global_quantities_desc.volume = params_eff.data_request;
+params_eff.data_request = 'w_mhd';
+global_quantities.w_mhd = gdat(params_summary.shot,params_eff);
+global_quantities_desc.w_mhd = params_eff.data_request;
+
+global_quantities_fieldnames = fieldnames(global_quantities);
+special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
+for it=1:ntime
+  for i=1:numel(global_quantities_fieldnames)
+    if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
+      if ~isstruct(ids_summary.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}))
+        ids_summary.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = ...
+            global_quantities.(global_quantities_fieldnames{i}).data(it);
+      else
+        special_fields{end+1} = global_quantities_fieldnames{i};
+      end
+    end
+  end
+end
+
+% special case
+for it=1:ntime
+  ids_summary.time_slice{it}.global_quantities.magnetic_axis.r = temp.r_magnetic_axis.data(it);
+  ids_summary.time_slice{it}.global_quantities.magnetic_axis.z = temp.z_magnetic_axis.data(it);
+  ids_summary.time_slice{it}.global_quantities.psi_axis = temp.psi_axis.data(it) + ...
+      ids_summary.time_slice{it}.global_quantities.psi_boundary;
+  [zz,izz]=min(temp.q_rho.data(:,it));
+  ids_summary.time_slice{it}.global_quantities.q_min.value = zz;
+  ids_summary.time_slice{it}.global_quantities.q_min.rho_tor_norm = temp.q_rho.grids_1d.rhotornorm(izz);
+end
+
+% for boundary in addition to lcfs
+% active_limiter_point = gdat(params_summary.shot,'active_limiter_point','machine',gdat_params.machine);
+params_eff.data_request = 'kappa';
+boundary.elongation = gdat(params_summary.shot,params_eff);
+boundary_desc.elongation = params_eff.data_request;
+% elongation_lower = gdat(params_summary.shot,'elongation_lower','machine',gdat_params.machine);
+% elongation_upper = gdat(params_summary.shot,'elongation_upper','machine',gdat_params.machine);
+params_eff.data_request = 'a_minor';
+boundary.minor_radius = gdat(params_summary.shot,params_eff);
+boundary_desc.minor_radius = params_eff.data_request;
+% squareness_lower_inner = gdat(params_summary.shot,'squareness_lower_inner','machine',gdat_params.machine);
+% squareness_lower_outer = gdat(params_summary.shot,'squareness_lower_outer','machine',gdat_params.machine);
+% squareness_upper_inner = gdat(params_summary.shot,'squareness_upper_inner','machine',gdat_params.machine);
+% squareness_upper_outer = gdat(params_summary.shot,'squareness_upper_outer','machine',gdat_params.machine);
+% strike_point = gdat(params_summary.shot,'strike_point','machine',gdat_params.machine);
+params_eff.data_request = 'delta';
+boundary.triangularity = gdat(params_summary.shot,params_eff);
+boundary_desc.triangularity = params_eff.data_request;
+params_eff.data_request = 'delta_bottom';
+boundary.triangularity_lower = gdat(params_summary.shot,params_eff);
+boundary_desc.triangularity_lower = params_eff.data_request;
+params_eff.data_request = 'delta_top';
+boundary.triangularity_upper = gdat(params_summary.shot,params_eff);
+boundary_desc.triangularity_upper = params_eff.data_request;
+params_eff.data_request = 'tcv_eq(''''n_xpts'''',''''liuqe.m'''')';
+temp.n_x_point = gdat(params_summary.shot,params_eff);
+temp_desc.n_x_point = params_eff.data_request;
+params_eff.data_request = 'tcv_eq(''''r_xpts'''',''''liuqe.m'''')';
+temp.r_x_point = gdat(params_summary.shot,params_eff);
+temp_desc.r_x_point = params_eff.data_request;
+params_eff.data_request = 'tcv_eq(''''z_xpts'''',''''liuqe.m'''')';
+temp.z_x_point = gdat(params_summary.shot,params_eff);
+temp_desc.z_x_point = params_eff.data_request;
+params_eff.data_request = 'rgeom';
+temp.rgeom = gdat(params_summary.shot,params_eff);
+temp_desc.rgeom = params_eff.data_request;
+params_eff.data_request = 'zgeom';
+temp.zgeom = gdat(params_summary.shot,params_eff);
+temp_desc.zgeom = params_eff.data_request;
+params_eff.data_request = 'r_contour_edge';
+temp.r_lcfs = gdat(params_summary.shot,params_eff);
+temp_desc.r_lcfs = params_eff.data_request;
+params_eff.data_request = 'z_contour_edge';
+temp.z_lcfs = gdat(params_summary.shot,params_eff);
+temp_desc.z_lcfs = params_eff.data_request;
+
+boundary_fieldnames = fieldnames(boundary);
+special_fields = {'lcfs', 'geometric_axis', 'x_point'}; % fields needing non-automatic treatments
+for it=1:ntime
+  for i=1:numel(boundary_fieldnames)
+    if ~any(strcmp(boundary_fieldnames{i},special_fields))
+      if ~isstruct(ids_summary.time_slice{it}.boundary.(boundary_fieldnames{i}))
+        ids_summary.time_slice{it}.boundary.(boundary_fieldnames{i}) = ...
+            boundary.(boundary_fieldnames{i}).data(it);
+      else
+        special_fields{end+1} = boundary_fieldnames{i};
+      end
+    end
+  end
+end
+
+% special cases
+for it=1:ntime
+  ids_summary.time_slice{it}.boundary.outline.r = temp.r_lcfs.data(:,it);
+  ids_summary.time_slice{it}.boundary.outline.z = temp.z_lcfs.data(:,it);
+  ids_summary.time_slice{it}.boundary.lcfs.r = ids_summary.time_slice{it}.boundary.outline.r;
+  ids_summary.time_slice{it}.boundary.lcfs.z = ids_summary.time_slice{it}.boundary.outline.z;
+  ids_summary.time_slice{it}.boundary.geometric_axis.r = temp.rgeom.data(it);
+  ids_summary.time_slice{it}.boundary.geometric_axis.z = temp.zgeom.data(it);
+  if temp.n_x_point.data(it) > 0
+    % not that asking for time_out may lead to interpolated nb of X-points non integer, should included piece-wise constant interpolation for equil quantities...
+    ids_summary.time_slice{it}.boundary.x_point(1:fix(temp.n_x_point.data(it))) = ids_summary.time_slice{it}.boundary.x_point(1);
+    for i=1:fix(temp.n_x_point.data(it))
+      ids_summary.time_slice{it}.boundary.x_point{i}.r = temp.r_x_point.data(i,it);
+      ids_summary.time_slice{it}.boundary.x_point{i}.z = temp.z_x_point.data(i,it);
+    end
+  else
+    ids_summary.time_slice{it}.boundary.x_point = {};
+  end
+end
+
+%% constraints
+% TODO: Add description
+
+% Measured values
+liuqe_time = ids_summary.time; % Not necessarily magnetics time so far
+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);
+flux = mdsvalue('tcv_idealloop("flux")[$1,*]',mag_time_req);
+diam = mdsvalue('\results::dmlcor[$1]',mag_time_req);
+ip   = mdsvalue('\magnetics::iplasma:trapeze[$1]',mag_time_req);
+% 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',...
+           'F_001','F_002','F_003','F_004','F_005','F_006','F_007','F_008',...
+           'G_001','G_001','G_001','G_001','G_001','G_001'};
+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
+
+% 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")');
+flux_liuqe = mdsvalue('tcv_eq("psi_loop","liuqe.m")');
+diam_liuqe = mdsvalue('tcv_eq("tor_flux_dml","liuqe.m")');
+ip_liuqe   = mdsvalue('tcv_eq("i_pl","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
+
+% Weights (using old parameters tree for now)
+bpol_err = mdsvalue('\results::parameters:berr')./mdsvalue('\results::parameters:vvv[0:37]');
+flux_err = mdsvalue('\results::parameters:ferr')./mdsvalue('\results::parameters:www[0:37]')*2*pi;
+diam_err = 0.13e-3./mdsvalue('\results::parameters:idml');
+ip_err   = mdsvalue('\results::parameters:plcerr')*1e3;
+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;
+
+if ntime > 0
+  constraints_orig = ids_summary.time_slice{1}.constraints;
+  % Remove unused arrays
+  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
+end
+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}.time_measurement = mag_time(itime(it));
+    bpol_probe{ib}.exact = 0;
+    bpol_probe{ib}.weight = 1/(bpol_err(ib)).^2;
+    bpol_probe{ib}.reconstructed = bpol_liuqe(ib,it);
+  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}.time_measurement = mag_time(itime(it));
+    flux_loop{il}.exact = 0;
+    flux_loop{il}.weight = 1/(flux_err(il)).^2;
+    flux_loop{il}.reconstructed = flux_liuqe(il,it);
+  end
+  constraints.flux_loop = flux_loop;
+  % ip
+  constraints.ip.measured = ip(it);
+  constraints.ip.source = 'IDS:magnetics/method[1]/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.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}.time_measurement = mag_time(itime(it));
+    if strcmp(dim_pol{ic},'TOR_001')
+      pf_current{ic}.source = [pf_current{ic}.source,' replaced with 0'];
+      pf_current{ic}.exact = 1;
+    else
+      pf_current{ic}.exact = 0;
+      pf_current{ic}.weight = 1/(ipol_err(ic)).^2;
+      pf_current{ic}.reconstructed = ipol_liuqe(ic,it);
+    end
+  end
+  constraints.pf_current = pf_current;
+
+  ids_summary.time_slice{it}.constraints = constraints;
+end
+
+
+%
+%% profiles_1d (cannot use eqdsk since not same radial mesh)
+%
+% area = gdat(params_summary.shot,'area','machine',gdat_params.machine);
+% b_average = gdat(params_summary.shot,'b_average','machine',gdat_params.machine);
+% beta_pol = gdat(params_summary.shot,'beta_pol','machine',gdat_params.machine);
+% b_field_average = gdat(params_summary.shot,'b_field_average','machine',gdat_params.machine);
+% b_field_max = gdat(params_summary.shot,'b_field_max','machine',gdat_params.machine);
+% b_field_min = gdat(params_summary.shot,'b_field_min','machine',gdat_params.machine);
+% b_max = gdat(params_summary.shot,'b_max','machine',gdat_params.machine);
+% b_min = gdat(params_summary.shot,'b_min','machine',gdat_params.machine);
+% darea_dpsi = gdat(params_summary.shot,'darea_dpsi','machine',gdat_params.machine);
+% darea_drho_tor = gdat(params_summary.shot,'darea_drho_tor','machine',gdat_params.machine);
+params_eff.data_request = 'pprime';
+profiles_1d.dpressure_dpsi = gdat(params_summary.shot,params_eff);
+profiles_1d_desc.dpressure_dpsi = params_eff.data_request;
+% dpsi_drho_tor = gdat(params_summary.shot,'dpsi_drho_tor','machine',gdat_params.machine);
+% dvolume_dpsi = gdat(params_summary.shot,'dvolume_dpsi','machine',gdat_params.machine);
+% dvolume_drho_tor = gdat(params_summary.shot,'dvolume_drho_tor','machine',gdat_params.machine);
+% elongation = gdat(params_summary.shot,'elongation','machine',gdat_params.machine);
+params_eff.data_request = 'ttprime';
+profiles_1d.f_df_dpsi = gdat(params_summary.shot,params_eff);
+profiles_1d_desc.f_df_dpsi = [params_eff.data_request '* 0.996^2'];
+params_eff.data_request = 'rbphi_rho';
+profiles_1d.f = gdat(params_summary.shot,params_eff);
+profiles_1d_desc.f = [params_eff.data_request '* 0.996'];
+profiles_1d.f.data = 0.996 * profiles_1d.f.data;
+profiles_1d.f_df_dpsi.data = 0.996.^2 .* profiles_1d.f_df_dpsi.data;
+% geometric_axis = gdat(params_summary.shot,'geometric_axis','machine',gdat_params.machine);
+% gm1 = gdat(params_summary.shot,'gm1','machine',gdat_params.machine);
+% gm2 = gdat(params_summary.shot,'gm2','machine',gdat_params.machine);
+% gm3 = gdat(params_summary.shot,'gm3','machine',gdat_params.machine);
+% gm4 = gdat(params_summary.shot,'gm4','machine',gdat_params.machine);
+% gm5 = gdat(params_summary.shot,'gm5','machine',gdat_params.machine);
+% gm6 = gdat(params_summary.shot,'gm6','machine',gdat_params.machine);
+% gm7 = gdat(params_summary.shot,'gm7','machine',gdat_params.machine);
+% gm8 = gdat(params_summary.shot,'gm8','machine',gdat_params.machine);
+% gm9 = gdat(params_summary.shot,'gm9','machine',gdat_params.machine);
+% j_parallel = gdat(params_summary.shot,'j_parallel','machine',gdat_params.machine);
+% j_tor = gdat(params_summary.shot,'j_tor','machine',gdat_params.machine);
+% magnetic_shear = gdat(params_summary.shot,'magnetic_shear','machine',gdat_params.machine);
+% mass_density = gdat(params_summary.shot,'mass_density','machine',gdat_params.machine);
+params_eff.data_request = 'phi_tor';
+profiles_1d.phi = gdat(params_summary.shot,params_eff);
+profiles_1d.phi.data = 0.996 * profiles_1d.phi.data;
+profiles_1d_desc.phi = [params_eff.data_request '* 0.996'];
+params_eff.data_request = 'pressure';
+profiles_1d.pressure = gdat(params_summary.shot,params_eff);
+profiles_1d_desc.pressure = params_eff.data_request;
+% psi = gdat(params_summary.shot,'psi_rho','machine',gdat_params.machine); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities)
+params_eff.data_request = 'q_rho';
+profiles_1d.q = gdat(params_summary.shot,params_eff);
+profiles_1d_desc.q = params_eff.data_request;
+params_eff.data_request = 'rhotor';
+profiles_1d.rho_tor = gdat(params_summary.shot,params_eff);
+profiles_1d_desc.rho_tor = params_eff.data_request;
+%rho_tor_norm = gdat(params_summary.shot,'rhotor_norm','machine',gdat_params.machine); % from rho_tor
+params_eff.data_request = 'rhovol';
+profiles_1d.rho_volume_norm = gdat(params_summary.shot,params_eff);
+profiles_1d_desc.rho_volume_norm = params_eff.data_request;
+% r_inboard = gdat(params_summary.shot,'r_inboard','machine',gdat_params.machine);
+% r_outboard = gdat(params_summary.shot,'r_outboard','machine',gdat_params.machine);
+% squareness_lower_inner = gdat(params_summary.shot,'squareness_lower_inner','machine',gdat_params.machine);
+% squareness_lower_outer = gdat(params_summary.shot,'squareness_lower_outer','machine',gdat_params.machine);
+% squareness_upper_inner = gdat(params_summary.shot,'squareness_upper_inner','machine',gdat_params.machine);
+% squareness_upper_outer = gdat(params_summary.shot,'squareness_upper_outer','machine',gdat_params.machine);
+% surface = gdat(params_summary.shot,'surface','machine',gdat_params.machine);
+% trapped_fraction = gdat(params_summary.shot,'trapped_fraction','machine',gdat_params.machine);
+% triangularity_lower = gdat(params_summary.shot,'triangularity_lower','machine',gdat_params.machine);
+% triangularity_upper = gdat(params_summary.shot,'triangularity_upper','machine',gdat_params.machine);
+params_eff.data_request = 'volume_rho';
+profiles_1d.volume = gdat(params_summary.shot,params_eff);
+profiles_1d_desc.volume = params_eff.data_request;
+
+liuqe_opt = gdat_params.liuqe; % default at this stage but could take from gdat params like error bar
+switch liuqe_opt
+ case {-1},   psitbx_str='FBTE';
+ case {1,21}, psitbx_str='LIUQE.M';
+ case {11},   psitbx_str='LIUQE';
+ case {2, 3, 22, 23}, psitbx_str=['LIUQE.M' num2str(mod(liuqe_opt,10))];
+ case {12,13}, psitbx_str=['LIUQE' num2str(mod(liuqe_opt,10))];
+ otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_opt)]);
+end
+fsd = psitbxtcv2(shot,profiles_1d.volume.t,'FS',psitbx_str,false); % will get automatically the correct time interval
+grho_metric_3D = metric(fsd,-1);
+% Introduced new anonymous function to compute FS average ...
+metric_FS = metric(grho_metric_3D.grid,[2,3]);
+denom=sum(metric_FS./grho_metric_3D,[2,3]);
+FS_av = @(x) sum(x.*metric_FS./grho_metric_3D,[2,3])./denom;
+R=metric(fsd,3);
+Rm2av=FS_av(1./R.^2);
+profiles_1d.gm1.data = Rm2av.x;
+profiles_1d_desc.gm1 = ['psitbxtcv2 with ' psitbx_str ' then Rm2av=FS_av(1./R.^2)'];
+%tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % this gives (grad rhopol/R)^2 not gm2 which is grad rhotor^2
+%profiles_1d.gm2.data = tmp_gm.x;
+tmp_gm = FS_av(1./R.^1);
+profiles_1d.gm9.data = tmp_gm.x;
+profiles_1d_desc.gm9 = 'FS_av(1./R.^1)';
+
+tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % grad rhopol^2 to get <grad psi^2>
+nrho = numel(profiles_1d.rho_tor.x);
+gradpsi_over_R_sq = NaN(nrho,ntime);
+for it=1:ntime
+  gradpsi_over_R_sq(:,it) = tmp_gm.x(:,it) .* 4 .* profiles_1d.volume.x.^2 .* ...
+      (ids_summary.time_slice{it}.global_quantities.psi_boundary-ids_summary.time_slice{it}.global_quantities.psi_axis).^2;
+end
+mu0 = 4.e-7 * pi;
+% Eq. (30) cocos paper cocos=17
+% j_tor=<jphi/R>/<1/R>=-sigma_Bp (2pi)^e_Bp dp/dpsi / <1/R> - sigma_Bp (2pi)^e_Bp F dF/dpsi / mu0 <1/R^2> / <1/R>
+% simaBp=-1 and eBp=1 for cocos=17 from TCV LIUQE
+profiles_1d.j_tor.data = - (-1.) .* 2.*pi .* profiles_1d.dpressure_dpsi.data ./ profiles_1d.gm9.data ...
+    - (-1.) .* 2.*pi .* profiles_1d.gm1.data ./ profiles_1d.gm9.data .* profiles_1d.f_df_dpsi.data ./ mu0;
+%
+% <j.B> = - sigma_Bp (2pi)^e_Bp dp/dpsi F - sigma_Bp F dF/dpsi / mu0 [ (2pi)^e_Bp F <1/R^2> + 1/(2pi)^e_Bp * <|grad psi|^2/R^2> / F ]
+% simaBp=-1 and eBp=1 for cocos=17 from TCV LIUQE
+%
+j_par = - (-1.) .* 2*pi .* profiles_1d.dpressure_dpsi.data .* profiles_1d.f.data ...
+        - (-1.) .* profiles_1d.f_df_dpsi.data ./ mu0 .* ...
+        ( (2.*pi) .* profiles_1d.f.data .* profiles_1d.gm1.data + 1./(2.*pi) .* gradpsi_over_R_sq ./ profiles_1d.f.data);
+profiles_1d.j_parallel.data = j_par./repmat(ids_summary.vacuum_toroidal_field.b0',size(profiles_1d.f.data,1),1);
+
+profiles_1d_fieldnames = fieldnames(profiles_1d);
+special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
+for it=1:ntime
+  for i=1:numel(profiles_1d_fieldnames)
+    if ~any(strcmp(profiles_1d_fieldnames{i},special_fields))
+      if ~isstruct(ids_summary.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}))
+        if ~ischar(profiles_1d.(profiles_1d_fieldnames{i}).data) && ~isempty(profiles_1d.(profiles_1d_fieldnames{i}).data) ...
+              && size(profiles_1d.(profiles_1d_fieldnames{i}).data,2)>=it
+          ids_summary.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}) = ...
+              profiles_1d.(profiles_1d_fieldnames{i}).data(:,it);
+        end
+      else
+        special_fields{end+1} = profiles_1d_fieldnames{i};
+      end
+    end
+  end
+end
+
+% special cases
+for it=1:ntime
+  ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_field_tor = ids_summary.time_slice{it}.profiles_1d.f(1) ...
+      ./ids_summary.time_slice{it}.global_quantities.magnetic_axis.r;
+  ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_tor = ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_field_tor;
+  ids_summary.time_slice{it}.profiles_1d.rho_tor_norm = ids_summary.time_slice{it}.profiles_1d.rho_tor./ ...
+      ids_summary.time_slice{it}.profiles_1d.rho_tor(end);
+  ids_summary.time_slice{it}.profiles_1d.psi = ids_summary.time_slice{it}.global_quantities.psi_axis + ...
+      profiles_1d.rho_tor.x.^2 .* ...
+      (ids_summary.time_slice{it}.global_quantities.psi_boundary- ids_summary.time_slice{it}.global_quantities.psi_axis);
+end
+
+%
+%% profiles_2d{1} ala eqdsk, only this one thus grid_type=1
+%
+% b_field_r = gdat(params_summary.shot,'b_field_r','machine',gdat_params.machine);
+% b_field_tor = gdat(params_summary.shot,'b_field_tor','machine',gdat_params.machine);
+% b_field_z = gdat(params_summary.shot,'b_field_z','machine',gdat_params.machine);
+% b_r = gdat(params_summary.shot,'b_r','machine',gdat_params.machine);
+% b_tor = gdat(params_summary.shot,'b_tor','machine',gdat_params.machine);
+% b_z = gdat(params_summary.shot,'b_z','machine',gdat_params.machine);
+% grid = gdat(params_summary.shot,'grid','machine',gdat_params.machine); % special
+profiles_2d.grid_type.name = 'rectangular';
+profiles_2d.grid_type.index = 1;
+profiles_2d.grid_type.description = 'Cylindrical R,Z ala eqdsk';
+% j_parallel = gdat(params_summary.shot,'j_parallel','machine',gdat_params.machine);
+% j_tor = gdat(params_summary.shot,'j_tor','machine',gdat_params.machine);
+% phi = gdat(params_summary.shot,'phi','machine',gdat_params.machine);
+params_eff.data_request = 'psi';
+profiles_2d.psi = gdat(params_summary.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_summary.shot,'r','machine',gdat_params.machine); % not to be filled since in grid.dim1
+% theta = gdat(params_summary.shot,'theta','machine',gdat_params.machine);
+% z = gdat(params_summary.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
+for it=1:ntime
+  for i=1:numel(profiles_2d_fieldnames)
+    if ~any(strcmp(profiles_2d_fieldnames{i},special_fields))
+      if ~isstruct(ids_summary.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}))
+        if ~ischar(profiles_2d.(profiles_2d_fieldnames{i}).data) && ~isempty(profiles_2d.(profiles_2d_fieldnames{i}).data) ...
+              && size(profiles_2d.(profiles_2d_fieldnames{i}).data,3)>=it
+          ids_summary.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}) = ...
+              profiles_2d.(profiles_2d_fieldnames{i}).data(:,:,it);
+        end
+      else
+        special_fields{end+1} = profiles_2d_fieldnames{i};
+      end
+    end
+  end
+end
+
+% special cases
+for it=1:ntime
+  ids_summary.time_slice{it}.profiles_2d{1}.grid_type.name = profiles_2d.grid_type.name;
+  ids_summary.time_slice{it}.profiles_2d{1}.grid_type.index = profiles_2d.grid_type.index;
+  ids_summary.time_slice{it}.profiles_2d{1}.grid_type.description = profiles_2d.grid_type.description;
+  ids_summary.time_slice{it}.profiles_2d{1}.grid.dim1 = profiles_2d.psi.dim{1};
+  ids_summary.time_slice{it}.profiles_2d{1}.grid.dim2 = profiles_2d.psi.dim{2};
+  ids_summary.time_slice{it}.profiles_2d{1}.psi(:,:) = ids_summary.time_slice{it}.profiles_2d{1}.psi(:,:) + ...
+      global_quantities.psi_boundary.data(it);
+end
+
+% make arrays not filled in empty:
+ids_summary.grids_ggd = {};
+for it=1:numel(ids_summary.time_slice)
+  ids_summary.time_slice{it}.ggd = {};
+  ids_summary.time_slice{it}.boundary.strike_point = {};
+  ids_summary.time_slice{it}.boundary_separatrix.x_point = {};
+  ids_summary.time_slice{it}.boundary_separatrix.strike_point = {};
+end
+
+% special test matrix cocos transform
+% $$$ ldim1=129;
+% $$$ ldim2=257;
+% $$$ it=1;
+% $$$ ids_summary.time_slice{it}.coordinate_system.grid_type.index = 13;
+% $$$ ids_summary.time_slice{it}.coordinate_system.grid.dim1 = linspace(0,1,ldim1)';
+% $$$ ids_summary.time_slice{it}.coordinate_system.grid.dim2 = linspace(0,2*pi,ldim2);
+% $$$ ids_summary.time_slice{it}.coordinate_system.tensor_contravariant = 2.*ones(ldim1,ldim2,3,3);
+% $$$ ids_summary.time_slice{it}.coordinate_system.tensor_covariant = 0.5*ones(ldim1,ldim2,3,3);
+% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant = 13.*ones(ldim1,ldim2,3,3);
+% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant_error_upper = 14.*ones(ldim1,ldim2,3,3);
+% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant_error_lower = 12.*ones(ldim1,ldim2,3,3);
+% $$$ for it=1:2100
+% $$$   ids_summary.time_slice{it}.coordinate_system.g11_contravariant = 11.*ones(ldim1,ldim2,3,3);
+% $$$   ids_summary.time_slice{it}.coordinate_system.tensor_covariant = 0.5*ones(ldim1,ldim2,3,3);
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid_type.name = profiles_2d.grid_type.name;
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid_type.index = 11;
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid_type.description = profiles_2d.grid_type.description;
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1 = linspace(0,1,ldim1)';
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1_error_upper = 1.2.*linspace(0,1,ldim1)';
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1_error_lower = 0.8.*linspace(0,1,ldim1)';
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid.dim2 = linspace(0,2*pi,ldim2);
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.psi(:,:) = 11.*ones(ldim1,ldim2);
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.psi_error_upper(:,:) = 12.*ones(ldim1,ldim2);
+% $$$   ids_summary.time_slice{it}.profiles_2d{2}.psi_error_lower(:,:) = 10.*ones(ldim1,ldim2);
+% $$$ end
+
+% cocos automatic transform
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
+  [ids_summary,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_summary,'summary',gdat_params.cocos_in, ...
+          gdat_params.cocos_out,gdat_params.ipsign_out,gdat_params.b0sign_out,gdat_params.ipsign_in,gdat_params.b0sign_in, ...
+          gdat_params.error_bar,gdat_params.nverbose);
+end
diff --git a/matlab/TCV_IMAS/tcv_ids_headpart.m b/matlab/TCV_IMAS/tcv_ids_headpart.m
index d1d45e6e..fa90b8cf 100644
--- a/matlab/TCV_IMAS/tcv_ids_headpart.m
+++ b/matlab/TCV_IMAS/tcv_ids_headpart.m
@@ -53,25 +53,7 @@ ids_generic = ids_in;
 % ids_properties
 %
 subcall='';
-% $$$ if nargin>=5 && ~isempty(varargin)
-% $$$   for i=1:length(varargin)
-% $$$     if ~isempty(varargin{i})
-% $$$       if ischar(varargin{i})
-% $$$         subcall_add = [',' varargin{i}];
-% $$$       elseif isnumeric(varargin{i})
-% $$$         subcall_add = [',' num2str(varargin{i})];
-% $$$       elseif isstruct(varargin{i})
-% $$$         subcall_add = [',' 'struct_in'];
-% $$$       elseif iscell(varargin{i})
-% $$$         subcall_add = [',' 'cell_in'];
-% $$$       end
-% $$$     else
-% $$$       subcall_add = [',[]'];
-% $$$     end
-% $$$     subcall(end+1:end+length(subcall_add)) = subcall_add;
-% $$$   end
-% $$$   subcall(end+1:end+2) = '; ';
-% $$$ end
+
 ids_generic.ids_properties.comment = params_ids_generic.comment;
 if ~isempty(subcall)
   add_to_comment = ['varargin: ' subcall];
@@ -82,7 +64,13 @@ if ~isempty(subcall)
   end
 end
 
-ids_generic.ids_properties.provider = ['tcv_get_ids_' ids_name];
+%get gdat version
+[abcpath,abcfile,abcext]=fileparts(mfilename('fullpath'));
+[aa1,aa2]=unix(['cd ' abcpath '; git describe --tags --first-parent --abbrev=11 --long --dirty --always']);
+git_release_hash='gdat git hash not found';
+if aa1==0; git_release_hash = aa2(1:end-1); end % avoid newline
+
+ids_generic.ids_properties.provider = ['gdat_git: ' git_release_hash '; tcv_get_ids_' ids_name];
 if ~isempty(params_ids_generic.gdat_params)
   if isfield(params_ids_generic.gdat_params,'cocos_in') && isfield(params_ids_generic.gdat_params,'cocos_out')
     comment_add = ['from TCV cocos_in=' num2str(params_ids_generic.gdat_params.cocos_in) ' transformed to cocos_out=' num2str(params_ids_generic.gdat_params.cocos_out)];
diff --git a/matlab/ids_summary_disruption.m b/matlab/ids_summary_disruption.m
new file mode 100644
index 00000000..bed33a8e
--- /dev/null
+++ b/matlab/ids_summary_disruption.m
@@ -0,0 +1,30 @@
+function [ids_summary_out,ids_summary_out_description] = ids_summary_disruption(shot, ids_summary_in, gdat_params, varargin)
+%
+%    [ids_summary_out,ids_summary_out_description] = tcv_ids_ip(shot, ids_summary_in, gdat_params, varargin);
+%
+% fill in subnode "disruption" expected by ids_summary, for machine gdat_params.machine
+%
+% error_bar: from gdat_params.error_bar
+%            'delta' (default): error_bar to be added inserted in "upper" only as mentioned in description
+%            'delta_with_lower' : error_bar (abs) inserted in both lower and upper
+%            'added': value already added to data: upper/lower = data +/- error_bar
+%
+
+error_bar = 'delta';
+if exist('gdat_params','var') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar)
+  error_bar = gdat_params.error_bar;
+end
+params_eff_ref = gdat_params; params_eff_ref.doplot=0;
+try params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def)
+
+% Get data
+switch lower(params_eff_ref.machine)
+ case 'tcv'
+ otherwise
+end
+
+% Preallocate dimension
+ids_summary_out = ids_summary_in;
+ids_summary_out_description.disruption.time = 'time of disruption obtained from...to fill in description per subnode';
+
+% Put data into ids structure
-- 
GitLab


From 42c5c46676fce5208d53b5efb84943b4febc85aa Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Wed, 15 Apr 2020 15:13:10 +0200
Subject: [PATCH 2/4] make B0 consistent with 0.996 between iphi, liuqe.f and
 liuqe.m

---
 matlab/TCV/gdat_tcv.m | 56 +++++++++++++++++++++++++++++++++++++++----
 1 file changed, 51 insertions(+), 5 deletions(-)

diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index fda7a801..6f1d4ca3 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -829,23 +829,29 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       if strcmp(lower(gdat_data.gdat_params.source),'iphi')
         nodenameeff=['\magnetics::iphi'];
         tracetdi=tdi(nodenameeff);
-        gdat_data.data=192.E-07 * 0.996 *tracetdi.data/r0exp;
+        added_correction = 0.996; % correction already in liuqe.f
+        added_correction_str = ['96*mu0/2piR0 * ' num2str(added_correction) ' * '];
+        gdat_data.data=192.E-07 * added_correction *tracetdi.data/r0exp;
       else
         if liuqe_matlab==0
+          added_correction = 1.0; % correction already in liuqe.f
+          added_correction_str = [''];
           nodenameeff = ['tcv_eq(''BZERO'',''LIUQE' substr_liuqe_tcv_eq ''')'];
         else
           if isempty(substr_liuqe); substr_liuqe = '_1'; end
           nodenameeff=['tcv_eq(''BZERO'',''LIUQE.M' substr_liuqe_tcv_eq ''')'];
+          added_correction = 0.996; % correction removed in liuqe.m at this stage
+          added_correction_str = [num2str(added_correction) ' * '];
         end
         tracetdi=tdi(nodenameeff);
-        gdat_data.data = tracetdi.data;
+        gdat_data.data = tracetdi.data .* added_correction;
       end
     end
     if isempty(tracetdi.data) || isempty(tracetdi.dim) % || ischar(tracetdi.data) (to add?)
       if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end
       return
     end
-    gdat_data.data_fullpath=[nodenameeff];
+    gdat_data.data_fullpath=[added_correction_str nodenameeff];
     gdat_data.dim = tracetdi.dim;
     gdat_data.t = gdat_data.dim{1};
     if any(strcmp(fieldnames(tracetdi),'units'))
@@ -862,7 +868,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-   case {'betan'}
+   case {'betan', 'beta_tor_norm'}
     % 100*beta / |Ip[MA] * B0[T]| * a[m]
     % get B0 from gdat_tcv, without re-opening the shot and using the same parameters except data_request
     % easily done thanks to structure call for options
@@ -1260,9 +1266,49 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     gdat_data.request_description = ['data=psi, x=(R,Z), eqdsk contains eqdsk structure with which ' ...
                     'plot_eqdsk, write_eqdsk, read_eqdsk can be used'];
 
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+   case {'gas', 'gas_flux', 'gas_request', 'gas_feedforward'}
+    params_eff = gdat_data.gdat_params;
+    params_eff.data_request = '\diagz::flux_gaz:piezo_1:flux';
+    gasflux = gdat_tcv(gdat_data.shot,params_eff);
+    gdat_data.gas_flux = gasflux;
+    params_eff.data_request = '\hybrid::mat_m_signals:output_019';
+    gasrequest = gdat_tcv(gdat_data.shot,params_eff);
+    gdat_data.gas_request_volt = gasrequest;
+    %volt_to_mlpers = max(0.,(in_volt-0.6879)*(72.41/(4.2673-0.6879))
+    gdat_data.gas_request_flux = gasrequest;
+    gdat_data.gas_request_flux.data = max(0.,72.41.*(gasrequest.data-0.6879)./(4.2673-0.6879));
+    gdat_data.gas_request_flux.units = gasflux.units;
+    params_eff.data_request = '\draw_feedfor_gas:alim_001';
+    gasrequest_ref = gdat_tcv(gdat_data.shot,params_eff); gasrequest_ref.units = 'V';
+    gdat_data.gas_feedforward_volt = gasrequest_ref;
+    gdat_data.gas_feedforward_flux = gasrequest_ref;
+    gdat_data.gas_feedforward_flux.data = max(0.,72.41.*(gasrequest_ref.data-0.6879)./(4.2673-0.6879));
+    gdat_data.label = data_request_eff;
+    switch data_request_eff
+      case {'gas', 'gas_flux'}
+       gdat_data.data = gdat_data.gas_flux.data;
+       gdat_data.units = gdat_data.gas_flux.units;
+       gdat_data.t = gdat_data.gas_flux.t;
+     case {'gas_request'}
+       gdat_data.data = gdat_data.gas_request_volt.data;
+       gdat_data.units = gdat_data.gas_request_volt.units;
+       gdat_data.t = gdat_data.gas_request_volt.t;
+     case {'gas_feedforward'}
+       gdat_data.data = gdat_data.gas_feedforward_volt.data;
+       gdat_data.units = gdat_data.gas_feedforward_volt.units;
+       gdat_data.t = gdat_data.gas_feedforward_volt.t;
+     otherwise
+      error('gas option not defined')
+    end
+    gdat_data.dim{1} = gdat_data.t;
+    gdat_data.dimunits{1} = 's';
+    if isfield(gdat_data.gdat_params,'time_out') && ~isempty(gdat_data.gdat_params.time_out)
+      [gdat_data] = gdat2time_out(gdat_data,1);
+    end
+
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'halphas'}
-  mapping.expression = '\base::pd:pd_001';
     channels = [1:18];
     if isfield(gdat_data.gdat_params,'channels') && ~isempty(gdat_data.gdat_params.channels)
       channels = gdat_data.gdat_params.channels;
-- 
GitLab


From 6e995a220cb034ea5f08a9da8066923ffdda61f4 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Tue, 21 Apr 2020 23:46:06 +0200
Subject: [PATCH 3/4] add gas, gas_flux, _request, etc

---
 matlab/TCV/tcv_requests_mapping.m | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/matlab/TCV/tcv_requests_mapping.m b/matlab/TCV/tcv_requests_mapping.m
index 90f5c9e8..d010ce02 100644
--- a/matlab/TCV/tcv_requests_mapping.m
+++ b/matlab/TCV/tcv_requests_mapping.m
@@ -69,18 +69,18 @@ switch lower(data_request)
   mapping.label = 'B_0';
   mapping.method = 'switchcase';
   mapping.expression = '';
- case 'beta'
+ case {'beta', 'beta_tor'}
   mapping.timedim = 1;
   mapping.label = '\beta';
   mapping.method = 'tdiliuqe';
   mapping.expression = '\results::beta_tor';
   mapping.expression = 'tcv_eq(''''beta_tor'''',''''LIUQE.M'''')';
- case 'betan'
+ case {'betan', 'beta_tor_norm'}
   mapping.timedim = 1;
   mapping.label = '\beta_N';
   mapping.method = 'switchcase';
   mapping.expression = '';
- case 'betap'
+ case {'beta_pol', 'betap'}
   mapping.timedim = 1;
   mapping.label = '\beta_p';
   mapping.method = 'tdiliuqe';
@@ -132,6 +132,11 @@ switch lower(data_request)
   mapping.timedim = 0;
   mapping.method = 'switchcase'; % could use function make_eqdsk directly?
   mapping.expression = '';
+ case {'gas', 'gas_flux', 'gas_request', 'gas_feedforward'}
+  mapping.timedim = 1;
+  mapping.label = 'gas flux';
+  mapping.method = 'switchcase';
+  mapping.expression = '';
  case 'halpha'
   mapping.timedim = 1;
   mapping.label = 'Halpha';
-- 
GitLab


From d9fe50d16a5b57c779c194252688d4fa5736517f Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Tue, 21 Apr 2020 23:53:24 +0200
Subject: [PATCH 4/4] add 1st version of ids_summary

---
 matlab/TCV_IMAS/tcv_get_ids_summary.m | 666 ++++++--------------------
 1 file changed, 144 insertions(+), 522 deletions(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_summary.m b/matlab/TCV_IMAS/tcv_get_ids_summary.m
index b804946f..0f543994 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_summary.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_summary.m
@@ -39,10 +39,20 @@ ids_summary.code.repository = git_url;
 % use equilibrium time as reference, so load first elongation(time)
 
 params_eff = params_eff_ref;
-params_eff.data_request='kappa'; params_eff.source='liuqe'; % to get liuqe time array
+params_eff.data_request='ip_trapeze';
+ip_trapeze = gdat(params_summary.shot,params_eff);
+ipsign = mdsvalue('\pcs::data:iohfb');
+i_t_end = find(ipsign.*ip_trapeze.data(end:-1:1) > 5e3); % an offset might cause a problem...
+i_t_end = length(ip_trapeze.t)-i_t_end(1)+2;
+t_end = ip_trapeze.t(i_t_end);
+b0sign = mdsvalue('\pcs::data:if36fb');
+params_eff = params_eff_ref;
+params_eff.data_request='kappa';
 kappa = gdat(params_summary.shot,params_eff);
-ids_summary.time = kappa.t;
-ids_summary_description.time = 'time from gdat(shot,''kappa'')';
+% to get liuqe time array for main time, but add points before and after according to ip_trapeze, to have other data when available
+dt_caract = 1e-3;
+ids_summary.time = unique([[0:dt_caract:kappa.t(1)]'; kappa.t; [kappa.t(end):dt_caract:t_end]']);
+ids_summary_description.time = 'time from gdat(shot,''kappa'') and adding from 0 up to t(ip_trapeze(end))';
 
 params_eff = params_eff_ref;
 % boundary:
@@ -81,11 +91,12 @@ boundary_desc.triangularity_upper = params_eff.data_request;
 
 boundary_fieldnames = fieldnames(boundary);
 special_fields = {}; % fields needing non-automatic treatments
-% assume all boundary have same time base since called with same gdat_params for liuqe
 for i=1:numel(boundary_fieldnames)
   if ~any(strcmp(boundary_fieldnames{i},special_fields))
     if ~isstruct(ids_summary.boundary.(boundary_fieldnames{i}).value)
       ids_summary.boundary.(boundary_fieldnames{i}).value = boundary.(boundary_fieldnames{i}).data;
+      ids_summary.boundary.(boundary_fieldnames{i}).value = interp1(boundary.(boundary_fieldnames{i}).t, ...
+          boundary.(boundary_fieldnames{i}).data,ids_summary.time,'linear',NaN);
       ids_summary.boundary.(boundary_fieldnames{i}).source = ['gdat request: ' boundary_desc.(boundary_fieldnames{i})];
     else
       special_fields{end+1} = boundary_fieldnames{i};
@@ -112,549 +123,160 @@ end
 
 [ids_summary_out,ids_summary_out_description] = ids_summary_disruption(shot, ids_summary, gdat_params);
 
-return
-% to continue later...
-
-params_eff = params_eff_ref;
-params_eff.data_request='b0'; params_eff.source='liuqe'; % to get liuqe time array
-bb = gdat(params_summary.shot,params_eff);
-params_eff = rmfield(params_eff,'source'); % to get original magnetics data
-vacuum_toroidal_field.b0=gdat(params_summary.shot,params_eff);
-ij_ok = isfinite(vacuum_toroidal_field.b0.data);
-vacuum_toroidal_field.b0.data = interpos(21,vacuum_toroidal_field.b0.t(ij_ok),vacuum_toroidal_field.b0.data(ij_ok),bb.t);
-vacuum_toroidal_field.b0.t = bb.t;
-vacuum_toroidal_field.b0.dim = {vacuum_toroidal_field.b0.t};
-vacuum_toroidal_field_desc.b0 = ['''b0'',''timing source'',''liuqe=' num2str(gdat_params.liuqe) ''''];
-vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'',''source'',''liuqe'']';
-ids_summary.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
-ids_summary.vacuum_toroidal_field.b0 = vacuum_toroidal_field.b0.data;
-ids_summary_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
-
-ids_summary.time = vacuum_toroidal_field.b0.t;
-ids_summary_description.time = '.t subfield from: [''b0'',''source'',''liuqe'']';
-ntime = numel(ids_summary.time);
-
-ids_summary.time_slice(1:ntime) = ids_summary.time_slice(1);
+% gas: bottom/top?, deuterium, helium_xx, ?
+params_eff.data_request = 'gas_flux';
+gas_injection_rates.deuterium = gdat(params_summary.shot,params_eff);
+gas_injection_rates_desc.deuterium = params_eff.data_request;
 
-% load time array data to copy to time_slices
-
-% global_quantities data into local global_quantities.* structure with correct end names and global_quantities_desc.* with description. Use temp.* and temp_desc.* structures for temporary data
-
-% brute force solution load all eqdsks
-% $$$ for it=1:ntime
-% $$$   ids_summary.time(it)
-% $$$   temp.eqdsks{it}=gdat(params_summary.shot,'eqdsk','time',ids_summary.time(it),'write',0,'machine',gdat_params.machine);
-% $$$ end
-% $$$ temp_desc.eqdsks{1} = '''eqdsk'',''time'',ids_summary.time(it)';
-
-params_eff = params_eff_ref;
-params_eff.data_request = 'area_edge';
-global_quantities.area = gdat(params_summary.shot,params_eff);
-global_quantities_desc.area = params_eff.data_request;
-params_eff.data_request = 'betan';
-global_quantities.beta_normal = gdat(params_summary.shot,params_eff);
-global_quantities_desc.beta_normal = params_eff.data_request;
-params_eff.data_request = 'betap';
+gas_injection_rates_fieldnames = fieldnames(gas_injection_rates);
+special_fields = {}; % fields needing non-automatic treatments
+% assume all gas_injection_rates have same time base since called with same gdat_params for liuqe
+for i=1:numel(gas_injection_rates_fieldnames)
+  if ~any(strcmp(gas_injection_rates_fieldnames{i},special_fields))
+    if ~isstruct(ids_summary.gas_injection_rates.(gas_injection_rates_fieldnames{i}).value)
+      ids_summary.gas_injection_rates.(gas_injection_rates_fieldnames{i}).value = interp1( ...
+          gas_injection_rates.(gas_injection_rates_fieldnames{i}).t,gas_injection_rates.(gas_injection_rates_fieldnames{i}).data, ...
+          ids_summary.time,'linear',NaN);
+      ids_summary.gas_injection_rates.(gas_injection_rates_fieldnames{i}).source = ['gdat request: ' gas_injection_rates_desc.(gas_injection_rates_fieldnames{i})];
+    else
+      special_fields{end+1} = gas_injection_rates_fieldnames{i};
+    end
+  end
+end
+%params_eff.doplot=1;
+% global_quantities:
+params_eff.data_request = 'b0';
+global_quantities.b0 = gdat(params_summary.shot,params_eff); % do not add 'source','liuqe' to get all times
+global_quantities_desc.b0 = [params_eff.data_request ' ; ' global_quantities.b0.data_fullpath];
+params_eff.data_request = 'beta_pol';
 global_quantities.beta_pol = gdat(params_summary.shot,params_eff);
-global_quantities_desc.beta_pol = params_eff.data_request;
-params_eff.data_request = 'beta';
+global_quantities_desc.beta_pol = [params_eff.data_request ' ; ' global_quantities.beta_pol.data_fullpath];
+params_eff.data_request = 'beta_pol';
+global_quantities.beta_pol_mhd = gdat(params_summary.shot,params_eff);
+global_quantities_desc.beta_pol_mhd = [params_eff.data_request ' ; ' global_quantities.beta_pol_mhd.data_fullpath];
+params_eff.data_request = 'beta_tor';
 global_quantities.beta_tor = gdat(params_summary.shot,params_eff);
-global_quantities_desc.beta_tor = params_eff.data_request;
+global_quantities_desc.beta_tor = [params_eff.data_request ' ; ' global_quantities.beta_tor.data_fullpath];
+params_eff.data_request = 'beta_tor';
+global_quantities.beta_tor_mhd = gdat(params_summary.shot,params_eff);
+global_quantities_desc.beta_tor_mhd = [params_eff.data_request ' ; ' global_quantities.beta_tor_mhd.data_fullpath];
+params_eff.data_request = 'betan';
+global_quantities.beta_tor_norm = gdat(params_summary.shot,params_eff);
+global_quantities_desc.beta_tor_norm = [params_eff.data_request ' ; ' global_quantities.beta_tor_norm.data_fullpath];
+params_eff.data_request = 'betan';
+global_quantities.beta_tor_norm_mhd = gdat(params_summary.shot,params_eff);
+global_quantities_desc.beta_tor_norm_mhd = [params_eff.data_request ' ; ' global_quantities.beta_tor_norm_mhd.data_fullpath];
+params_eff.data_request = 'betan';
+global_quantities.beta_tor_thermal_norm = gdat(params_summary.shot,params_eff);
+global_quantities_desc.beta_tor_thermal_norm = [params_eff.data_request ' ; ' global_quantities.beta_tor_thermal_norm.data_fullpath];
+params_eff.data_request = '\tcv_shot::top.results.ibs:ibs';
+global_quantities.current_bootstrap = gdat(params_summary.shot,params_eff);
+global_quantities_desc.current_bootstrap = [params_eff.data_request ' ; ' global_quantities.current_bootstrap.data_fullpath];
+params_eff.data_request = '\tcv_shot::top.results.ibs:icd';
+global_quantities.current_non_inductive = gdat(params_summary.shot,params_eff);
+global_quantities_desc.current_non_inductive = [params_eff.data_request ' ; ' global_quantities.current_non_inductive.data_fullpath];
+params_eff.data_request = 'w_mhd';
+global_quantities.energy_diamagnetic = gdat(params_summary.shot,params_eff);
+global_quantities_desc.energy_diamagnetic = [params_eff.data_request ' ; ' global_quantities.energy_diamagnetic.data_fullpath];
+params_eff.data_request = '\tcv_shot::top.results.conf:we';
+global_quantities.energy_electrons_thermal = gdat(params_summary.shot,params_eff);
+global_quantities_desc.energy_electrons_thermal = [params_eff.data_request ' ; ' global_quantities.energy_electrons_thermal.data_fullpath];
+params_eff.data_request = '\tcv_shot::top.results.conf:wi';
+global_quantities.energy_ion_total_thermal = gdat(params_summary.shot,params_eff);
+global_quantities_desc.energy_ion_total_thermal = [params_eff.data_request ' ; ' global_quantities.energy_ion_total_thermal.data_fullpath];
 params_eff.data_request = 'w_mhd';
 global_quantities.energy_mhd = gdat(params_summary.shot,params_eff);
-global_quantities_desc.energy_mhd = params_eff.data_request;
+global_quantities_desc.energy_mhd = [params_eff.data_request ' ; ' global_quantities.energy_mhd.data_fullpath];
+params_eff.data_request = '.conf:we + .conf:wi';
+global_quantities.energy_thermal = global_quantities.energy_electrons_thermal;
+global_quantities.energy_thermal.data = global_quantities.energy_thermal.data + global_quantities.energy_ion_total_thermal.data;
+global_quantities.energy_thermal.data_fullpath = [global_quantities.energy_thermal.data_fullpath ' + wi'];
+global_quantities_desc.energy_thermal = [params_eff.data_request ' ; ' global_quantities.energy_thermal.data_fullpath];
+params_eff.data_request = 'ngf';
+global_quantities.greenwald_fraction = gdat(params_summary.shot,params_eff);
+global_quantities_desc.greenwald_fraction = [params_eff.data_request ' ; ' global_quantities.greenwald_fraction.data_fullpath];
+params_eff.data_request = 'h98y2';
+global_quantities.h_98 = gdat(params_summary.shot,params_eff);
+global_quantities_desc.h_98 = [params_eff.data_request ' ; ' global_quantities.h_98.data_fullpath];
 params_eff.data_request = 'ip';
 global_quantities.ip = gdat(params_summary.shot,params_eff);
-global_quantities_desc.ip = params_eff.data_request;
-% length_pol = gdat(params_summary.shot,'length_pol','machine',gdat_params.machine); % to be added
+global_quantities_desc.ip = [params_eff.data_request ' ; ' global_quantities.ip.data_fullpath];
 params_eff.data_request = 'li';
-global_quantities.li_3 = gdat(params_summary.shot,params_eff);
-global_quantities_desc.li_3 = params_eff.data_request;
-params_eff.data_request = 'r_axis';
-temp.r_magnetic_axis = gdat(params_summary.shot,params_eff);
-temp_desc.r_magnetic_axis = params_eff.data_request;
-params_eff.data_request = 'z_axis';
-temp.z_magnetic_axis = gdat(params_summary.shot,params_eff);
-temp_desc.z_magnetic_axis = params_eff.data_request;
-params_eff.data_request = 'psi_axis';
-temp.psi_axis = gdat(params_summary.shot,params_eff); % needs to add psi_edge sincepsi_axis liuqe assuming 0 dege value
-temp_desc.psi_axis = params_eff.data_request;
-params_eff.data_request = 'psi_edge';
-global_quantities.psi_boundary = gdat(params_summary.shot,params_eff);
-global_quantities_desc.psi_boundary = params_eff.data_request;
+global_quantities.li = gdat(params_summary.shot,params_eff);
+global_quantities_desc.li = [params_eff.data_request ' ; ' global_quantities.li.data_fullpath];
+params_eff.data_request = '\tcv_shot::top.results.conf:ptot_ohm';
+global_quantities.power_ohm = gdat(params_summary.shot,params_eff);
+global_quantities_desc.power_ohm = [params_eff.data_request ' ; ' global_quantities.power_ohm.data_fullpath];
+params_eff.data_request = '\results::bolo:prad:total';
+global_quantities.power_radiated = gdat(params_summary.shot,params_eff);
+if strcmp(lower(global_quantities.power_radiated.units),'kw')
+  global_quantities.power_radiated.data = global_quantities.power_radiated.data*1e3; %in kW
+end
+global_quantities_desc.power_radiated = [params_eff.data_request ' ; ' global_quantities.power_radiated.data_fullpath];
 params_eff.data_request = 'q95';
 global_quantities.q_95 = gdat(params_summary.shot,params_eff);
-global_quantities_desc.q_95 = params_eff.data_request;
-params_eff.data_request = 'q0';
-global_quantities.q_axis = gdat(params_summary.shot,params_eff); % will be checked with q_rho?
-global_quantities_desc.q_axis = params_eff.data_request;
-params_eff.data_request = 'q_rho';
-temp.q_rho = gdat(params_summary.shot,params_eff);
-temp_desc.q_rho = params_eff.data_request;
-% surface = gdat(params_summary.shot,'surface','machine',gdat_params.machine); % to be added
+global_quantities_desc.q_95 = [params_eff.data_request ' ; ' global_quantities.q_95.data_fullpath];
+params_eff.data_request = 'b0';
+global_quantities.b0 = gdat(params_summary.shot,params_eff);
+global_quantities_desc.b0 = [params_eff.data_request ' ; ' global_quantities.b0.data_fullpath];
+params_eff.data_request = 'r0';
+global_quantities.r0 = global_quantities.b0;
+global_quantities.r0.data = global_quantities.r0.r0;
+global_quantities.r0.dim = [];global_quantities.r0.t=[]; global_quantities.r0.dimunits = [];
+global_quantities.r0.units = 'm';
+global_quantities.r0.data_fullpath = ['gdat(b0).r0'];
+global_quantities_desc.r0 = [params_eff.data_request ' ; ' global_quantities.r0.data_fullpath];
+% $$$ params_eff.data_request = 'tau_tot';
+% $$$ global_quantities.tau_energy = gdat(params_summary.shot,params_eff);
+% $$$ global_quantities_desc.tau_energy = [params_eff.data_request ' ; ' global_quantities.tau_energy.data_fullpath];
+params_eff.data_request = 'vloop';
+global_quantities.v_loop = gdat(params_summary.shot,params_eff);
+tensvloop=-1e3;
+global_quantities.v_loop.data = interpos(global_quantities.v_loop.t,-global_quantities.v_loop.data,tensvloop);
+global_quantities_desc.v_loop = [params_eff.data_request ' ; -1* interpos of ' global_quantities.v_loop.data_fullpath ...
+                    ' with tens=' num2str(tensvloop)];
 params_eff.data_request = 'volume';
 global_quantities.volume = gdat(params_summary.shot,params_eff);
-global_quantities_desc.volume = params_eff.data_request;
-params_eff.data_request = 'w_mhd';
-global_quantities.w_mhd = gdat(params_summary.shot,params_eff);
-global_quantities_desc.w_mhd = params_eff.data_request;
+global_quantities_desc.volume = [params_eff.data_request ' ; ' global_quantities.volume.data_fullpath];
 
+special_fields = {'r0'}; % fields needing non-automatic treatments
 global_quantities_fieldnames = fieldnames(global_quantities);
-special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
-for it=1:ntime
-  for i=1:numel(global_quantities_fieldnames)
-    if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
-      if ~isstruct(ids_summary.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}))
-        ids_summary.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = ...
-            global_quantities.(global_quantities_fieldnames{i}).data(it);
-      else
-        special_fields{end+1} = global_quantities_fieldnames{i};
-      end
+for i=1:numel(global_quantities_fieldnames)
+  if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
+    if ~isstruct(ids_summary.global_quantities.(global_quantities_fieldnames{i}).value)
+      ids_summary.global_quantities.(global_quantities_fieldnames{i}).value = interp1( ...
+          global_quantities.(global_quantities_fieldnames{i}).t,global_quantities.(global_quantities_fieldnames{i}).data, ...
+          ids_summary.time,'linear',NaN);
+      ids_summary.global_quantities.(global_quantities_fieldnames{i}).source = ['gdat request: ' global_quantities_desc.(global_quantities_fieldnames{i})];
+    else
+      special_fields{end+1} = global_quantities_fieldnames{i};
     end
   end
 end
 
 % special case
-for it=1:ntime
-  ids_summary.time_slice{it}.global_quantities.magnetic_axis.r = temp.r_magnetic_axis.data(it);
-  ids_summary.time_slice{it}.global_quantities.magnetic_axis.z = temp.z_magnetic_axis.data(it);
-  ids_summary.time_slice{it}.global_quantities.psi_axis = temp.psi_axis.data(it) + ...
-      ids_summary.time_slice{it}.global_quantities.psi_boundary;
-  [zz,izz]=min(temp.q_rho.data(:,it));
-  ids_summary.time_slice{it}.global_quantities.q_min.value = zz;
-  ids_summary.time_slice{it}.global_quantities.q_min.rho_tor_norm = temp.q_rho.grids_1d.rhotornorm(izz);
-end
-
-% for boundary in addition to lcfs
-% active_limiter_point = gdat(params_summary.shot,'active_limiter_point','machine',gdat_params.machine);
-params_eff.data_request = 'kappa';
-boundary.elongation = gdat(params_summary.shot,params_eff);
-boundary_desc.elongation = params_eff.data_request;
-% elongation_lower = gdat(params_summary.shot,'elongation_lower','machine',gdat_params.machine);
-% elongation_upper = gdat(params_summary.shot,'elongation_upper','machine',gdat_params.machine);
-params_eff.data_request = 'a_minor';
-boundary.minor_radius = gdat(params_summary.shot,params_eff);
-boundary_desc.minor_radius = params_eff.data_request;
-% squareness_lower_inner = gdat(params_summary.shot,'squareness_lower_inner','machine',gdat_params.machine);
-% squareness_lower_outer = gdat(params_summary.shot,'squareness_lower_outer','machine',gdat_params.machine);
-% squareness_upper_inner = gdat(params_summary.shot,'squareness_upper_inner','machine',gdat_params.machine);
-% squareness_upper_outer = gdat(params_summary.shot,'squareness_upper_outer','machine',gdat_params.machine);
-% strike_point = gdat(params_summary.shot,'strike_point','machine',gdat_params.machine);
-params_eff.data_request = 'delta';
-boundary.triangularity = gdat(params_summary.shot,params_eff);
-boundary_desc.triangularity = params_eff.data_request;
-params_eff.data_request = 'delta_bottom';
-boundary.triangularity_lower = gdat(params_summary.shot,params_eff);
-boundary_desc.triangularity_lower = params_eff.data_request;
-params_eff.data_request = 'delta_top';
-boundary.triangularity_upper = gdat(params_summary.shot,params_eff);
-boundary_desc.triangularity_upper = params_eff.data_request;
-params_eff.data_request = 'tcv_eq(''''n_xpts'''',''''liuqe.m'''')';
-temp.n_x_point = gdat(params_summary.shot,params_eff);
-temp_desc.n_x_point = params_eff.data_request;
-params_eff.data_request = 'tcv_eq(''''r_xpts'''',''''liuqe.m'''')';
-temp.r_x_point = gdat(params_summary.shot,params_eff);
-temp_desc.r_x_point = params_eff.data_request;
-params_eff.data_request = 'tcv_eq(''''z_xpts'''',''''liuqe.m'''')';
-temp.z_x_point = gdat(params_summary.shot,params_eff);
-temp_desc.z_x_point = params_eff.data_request;
-params_eff.data_request = 'rgeom';
-temp.rgeom = gdat(params_summary.shot,params_eff);
-temp_desc.rgeom = params_eff.data_request;
-params_eff.data_request = 'zgeom';
-temp.zgeom = gdat(params_summary.shot,params_eff);
-temp_desc.zgeom = params_eff.data_request;
-params_eff.data_request = 'r_contour_edge';
-temp.r_lcfs = gdat(params_summary.shot,params_eff);
-temp_desc.r_lcfs = params_eff.data_request;
-params_eff.data_request = 'z_contour_edge';
-temp.z_lcfs = gdat(params_summary.shot,params_eff);
-temp_desc.z_lcfs = params_eff.data_request;
-
-boundary_fieldnames = fieldnames(boundary);
-special_fields = {'lcfs', 'geometric_axis', 'x_point'}; % fields needing non-automatic treatments
-for it=1:ntime
-  for i=1:numel(boundary_fieldnames)
-    if ~any(strcmp(boundary_fieldnames{i},special_fields))
-      if ~isstruct(ids_summary.time_slice{it}.boundary.(boundary_fieldnames{i}))
-        ids_summary.time_slice{it}.boundary.(boundary_fieldnames{i}) = ...
-            boundary.(boundary_fieldnames{i}).data(it);
-      else
-        special_fields{end+1} = boundary_fieldnames{i};
-      end
-    end
-  end
-end
-
-% special cases
-for it=1:ntime
-  ids_summary.time_slice{it}.boundary.outline.r = temp.r_lcfs.data(:,it);
-  ids_summary.time_slice{it}.boundary.outline.z = temp.z_lcfs.data(:,it);
-  ids_summary.time_slice{it}.boundary.lcfs.r = ids_summary.time_slice{it}.boundary.outline.r;
-  ids_summary.time_slice{it}.boundary.lcfs.z = ids_summary.time_slice{it}.boundary.outline.z;
-  ids_summary.time_slice{it}.boundary.geometric_axis.r = temp.rgeom.data(it);
-  ids_summary.time_slice{it}.boundary.geometric_axis.z = temp.zgeom.data(it);
-  if temp.n_x_point.data(it) > 0
-    % not that asking for time_out may lead to interpolated nb of X-points non integer, should included piece-wise constant interpolation for equil quantities...
-    ids_summary.time_slice{it}.boundary.x_point(1:fix(temp.n_x_point.data(it))) = ids_summary.time_slice{it}.boundary.x_point(1);
-    for i=1:fix(temp.n_x_point.data(it))
-      ids_summary.time_slice{it}.boundary.x_point{i}.r = temp.r_x_point.data(i,it);
-      ids_summary.time_slice{it}.boundary.x_point{i}.z = temp.z_x_point.data(i,it);
-    end
-  else
-    ids_summary.time_slice{it}.boundary.x_point = {};
-  end
-end
-
-%% constraints
-% TODO: Add description
-
-% Measured values
-liuqe_time = ids_summary.time; % Not necessarily magnetics time so far
-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);
-flux = mdsvalue('tcv_idealloop("flux")[$1,*]',mag_time_req);
-diam = mdsvalue('\results::dmlcor[$1]',mag_time_req);
-ip   = mdsvalue('\magnetics::iplasma:trapeze[$1]',mag_time_req);
-% 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',...
-           'F_001','F_002','F_003','F_004','F_005','F_006','F_007','F_008',...
-           'G_001','G_001','G_001','G_001','G_001','G_001'};
-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
+ids_summary.global_quantities.global_quantities.r0.value = global_quantities.r0.data;
+ids_summary.global_quantities.global_quantities.r0.source = ['gdat request: ' global_quantities_desc.r0];
 
-% 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")');
-flux_liuqe = mdsvalue('tcv_eq("psi_loop","liuqe.m")');
-diam_liuqe = mdsvalue('tcv_eq("tor_flux_dml","liuqe.m")');
-ip_liuqe   = mdsvalue('tcv_eq("i_pl","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
+% cst zeff
+params_eff.data_request = '\tcv_shot::top.results.ibs:z_eff';
+volume_average.zeff = gdat(params_summary.shot,params_eff);
+volume_average_desc.zeff = [params_eff.data_request ' ; ' volume_average.zeff.data_fullpath ' (not vol avrg but to match Iohm)'];
 
-% Weights (using old parameters tree for now)
-bpol_err = mdsvalue('\results::parameters:berr')./mdsvalue('\results::parameters:vvv[0:37]');
-flux_err = mdsvalue('\results::parameters:ferr')./mdsvalue('\results::parameters:www[0:37]')*2*pi;
-diam_err = 0.13e-3./mdsvalue('\results::parameters:idml');
-ip_err   = mdsvalue('\results::parameters:plcerr')*1e3;
-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;
-
-if ntime > 0
-  constraints_orig = ids_summary.time_slice{1}.constraints;
-  % Remove unused arrays
-  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
-end
-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}.time_measurement = mag_time(itime(it));
-    bpol_probe{ib}.exact = 0;
-    bpol_probe{ib}.weight = 1/(bpol_err(ib)).^2;
-    bpol_probe{ib}.reconstructed = bpol_liuqe(ib,it);
-  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}.time_measurement = mag_time(itime(it));
-    flux_loop{il}.exact = 0;
-    flux_loop{il}.weight = 1/(flux_err(il)).^2;
-    flux_loop{il}.reconstructed = flux_liuqe(il,it);
-  end
-  constraints.flux_loop = flux_loop;
-  % ip
-  constraints.ip.measured = ip(it);
-  constraints.ip.source = 'IDS:magnetics/method[1]/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.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}.time_measurement = mag_time(itime(it));
-    if strcmp(dim_pol{ic},'TOR_001')
-      pf_current{ic}.source = [pf_current{ic}.source,' replaced with 0'];
-      pf_current{ic}.exact = 1;
+special_fields = {}; % fields needing non-automatic treatments
+volume_average_fieldnames = fieldnames(volume_average);
+for i=1:numel(volume_average_fieldnames)
+  if ~any(strcmp(volume_average_fieldnames{i},special_fields))
+    if ~isstruct(ids_summary.volume_average.(volume_average_fieldnames{i}).value)
+      ids_summary.volume_average.(volume_average_fieldnames{i}).value = interp1( ...
+          volume_average.(volume_average_fieldnames{i}).t,volume_average.(volume_average_fieldnames{i}).data, ...
+          ids_summary.time,'linear',NaN);
+      ids_summary.volume_average.(volume_average_fieldnames{i}).source = ['gdat request: ' volume_average_desc.(volume_average_fieldnames{i})];
     else
-      pf_current{ic}.exact = 0;
-      pf_current{ic}.weight = 1/(ipol_err(ic)).^2;
-      pf_current{ic}.reconstructed = ipol_liuqe(ic,it);
+      special_fields{end+1} = volume_average_fieldnames{i};
     end
   end
-  constraints.pf_current = pf_current;
-
-  ids_summary.time_slice{it}.constraints = constraints;
-end
-
-
-%
-%% profiles_1d (cannot use eqdsk since not same radial mesh)
-%
-% area = gdat(params_summary.shot,'area','machine',gdat_params.machine);
-% b_average = gdat(params_summary.shot,'b_average','machine',gdat_params.machine);
-% beta_pol = gdat(params_summary.shot,'beta_pol','machine',gdat_params.machine);
-% b_field_average = gdat(params_summary.shot,'b_field_average','machine',gdat_params.machine);
-% b_field_max = gdat(params_summary.shot,'b_field_max','machine',gdat_params.machine);
-% b_field_min = gdat(params_summary.shot,'b_field_min','machine',gdat_params.machine);
-% b_max = gdat(params_summary.shot,'b_max','machine',gdat_params.machine);
-% b_min = gdat(params_summary.shot,'b_min','machine',gdat_params.machine);
-% darea_dpsi = gdat(params_summary.shot,'darea_dpsi','machine',gdat_params.machine);
-% darea_drho_tor = gdat(params_summary.shot,'darea_drho_tor','machine',gdat_params.machine);
-params_eff.data_request = 'pprime';
-profiles_1d.dpressure_dpsi = gdat(params_summary.shot,params_eff);
-profiles_1d_desc.dpressure_dpsi = params_eff.data_request;
-% dpsi_drho_tor = gdat(params_summary.shot,'dpsi_drho_tor','machine',gdat_params.machine);
-% dvolume_dpsi = gdat(params_summary.shot,'dvolume_dpsi','machine',gdat_params.machine);
-% dvolume_drho_tor = gdat(params_summary.shot,'dvolume_drho_tor','machine',gdat_params.machine);
-% elongation = gdat(params_summary.shot,'elongation','machine',gdat_params.machine);
-params_eff.data_request = 'ttprime';
-profiles_1d.f_df_dpsi = gdat(params_summary.shot,params_eff);
-profiles_1d_desc.f_df_dpsi = [params_eff.data_request '* 0.996^2'];
-params_eff.data_request = 'rbphi_rho';
-profiles_1d.f = gdat(params_summary.shot,params_eff);
-profiles_1d_desc.f = [params_eff.data_request '* 0.996'];
-profiles_1d.f.data = 0.996 * profiles_1d.f.data;
-profiles_1d.f_df_dpsi.data = 0.996.^2 .* profiles_1d.f_df_dpsi.data;
-% geometric_axis = gdat(params_summary.shot,'geometric_axis','machine',gdat_params.machine);
-% gm1 = gdat(params_summary.shot,'gm1','machine',gdat_params.machine);
-% gm2 = gdat(params_summary.shot,'gm2','machine',gdat_params.machine);
-% gm3 = gdat(params_summary.shot,'gm3','machine',gdat_params.machine);
-% gm4 = gdat(params_summary.shot,'gm4','machine',gdat_params.machine);
-% gm5 = gdat(params_summary.shot,'gm5','machine',gdat_params.machine);
-% gm6 = gdat(params_summary.shot,'gm6','machine',gdat_params.machine);
-% gm7 = gdat(params_summary.shot,'gm7','machine',gdat_params.machine);
-% gm8 = gdat(params_summary.shot,'gm8','machine',gdat_params.machine);
-% gm9 = gdat(params_summary.shot,'gm9','machine',gdat_params.machine);
-% j_parallel = gdat(params_summary.shot,'j_parallel','machine',gdat_params.machine);
-% j_tor = gdat(params_summary.shot,'j_tor','machine',gdat_params.machine);
-% magnetic_shear = gdat(params_summary.shot,'magnetic_shear','machine',gdat_params.machine);
-% mass_density = gdat(params_summary.shot,'mass_density','machine',gdat_params.machine);
-params_eff.data_request = 'phi_tor';
-profiles_1d.phi = gdat(params_summary.shot,params_eff);
-profiles_1d.phi.data = 0.996 * profiles_1d.phi.data;
-profiles_1d_desc.phi = [params_eff.data_request '* 0.996'];
-params_eff.data_request = 'pressure';
-profiles_1d.pressure = gdat(params_summary.shot,params_eff);
-profiles_1d_desc.pressure = params_eff.data_request;
-% psi = gdat(params_summary.shot,'psi_rho','machine',gdat_params.machine); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities)
-params_eff.data_request = 'q_rho';
-profiles_1d.q = gdat(params_summary.shot,params_eff);
-profiles_1d_desc.q = params_eff.data_request;
-params_eff.data_request = 'rhotor';
-profiles_1d.rho_tor = gdat(params_summary.shot,params_eff);
-profiles_1d_desc.rho_tor = params_eff.data_request;
-%rho_tor_norm = gdat(params_summary.shot,'rhotor_norm','machine',gdat_params.machine); % from rho_tor
-params_eff.data_request = 'rhovol';
-profiles_1d.rho_volume_norm = gdat(params_summary.shot,params_eff);
-profiles_1d_desc.rho_volume_norm = params_eff.data_request;
-% r_inboard = gdat(params_summary.shot,'r_inboard','machine',gdat_params.machine);
-% r_outboard = gdat(params_summary.shot,'r_outboard','machine',gdat_params.machine);
-% squareness_lower_inner = gdat(params_summary.shot,'squareness_lower_inner','machine',gdat_params.machine);
-% squareness_lower_outer = gdat(params_summary.shot,'squareness_lower_outer','machine',gdat_params.machine);
-% squareness_upper_inner = gdat(params_summary.shot,'squareness_upper_inner','machine',gdat_params.machine);
-% squareness_upper_outer = gdat(params_summary.shot,'squareness_upper_outer','machine',gdat_params.machine);
-% surface = gdat(params_summary.shot,'surface','machine',gdat_params.machine);
-% trapped_fraction = gdat(params_summary.shot,'trapped_fraction','machine',gdat_params.machine);
-% triangularity_lower = gdat(params_summary.shot,'triangularity_lower','machine',gdat_params.machine);
-% triangularity_upper = gdat(params_summary.shot,'triangularity_upper','machine',gdat_params.machine);
-params_eff.data_request = 'volume_rho';
-profiles_1d.volume = gdat(params_summary.shot,params_eff);
-profiles_1d_desc.volume = params_eff.data_request;
-
-liuqe_opt = gdat_params.liuqe; % default at this stage but could take from gdat params like error bar
-switch liuqe_opt
- case {-1},   psitbx_str='FBTE';
- case {1,21}, psitbx_str='LIUQE.M';
- case {11},   psitbx_str='LIUQE';
- case {2, 3, 22, 23}, psitbx_str=['LIUQE.M' num2str(mod(liuqe_opt,10))];
- case {12,13}, psitbx_str=['LIUQE' num2str(mod(liuqe_opt,10))];
- otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_opt)]);
 end
-fsd = psitbxtcv2(shot,profiles_1d.volume.t,'FS',psitbx_str,false); % will get automatically the correct time interval
-grho_metric_3D = metric(fsd,-1);
-% Introduced new anonymous function to compute FS average ...
-metric_FS = metric(grho_metric_3D.grid,[2,3]);
-denom=sum(metric_FS./grho_metric_3D,[2,3]);
-FS_av = @(x) sum(x.*metric_FS./grho_metric_3D,[2,3])./denom;
-R=metric(fsd,3);
-Rm2av=FS_av(1./R.^2);
-profiles_1d.gm1.data = Rm2av.x;
-profiles_1d_desc.gm1 = ['psitbxtcv2 with ' psitbx_str ' then Rm2av=FS_av(1./R.^2)'];
-%tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % this gives (grad rhopol/R)^2 not gm2 which is grad rhotor^2
-%profiles_1d.gm2.data = tmp_gm.x;
-tmp_gm = FS_av(1./R.^1);
-profiles_1d.gm9.data = tmp_gm.x;
-profiles_1d_desc.gm9 = 'FS_av(1./R.^1)';
-
-tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % grad rhopol^2 to get <grad psi^2>
-nrho = numel(profiles_1d.rho_tor.x);
-gradpsi_over_R_sq = NaN(nrho,ntime);
-for it=1:ntime
-  gradpsi_over_R_sq(:,it) = tmp_gm.x(:,it) .* 4 .* profiles_1d.volume.x.^2 .* ...
-      (ids_summary.time_slice{it}.global_quantities.psi_boundary-ids_summary.time_slice{it}.global_quantities.psi_axis).^2;
-end
-mu0 = 4.e-7 * pi;
-% Eq. (30) cocos paper cocos=17
-% j_tor=<jphi/R>/<1/R>=-sigma_Bp (2pi)^e_Bp dp/dpsi / <1/R> - sigma_Bp (2pi)^e_Bp F dF/dpsi / mu0 <1/R^2> / <1/R>
-% simaBp=-1 and eBp=1 for cocos=17 from TCV LIUQE
-profiles_1d.j_tor.data = - (-1.) .* 2.*pi .* profiles_1d.dpressure_dpsi.data ./ profiles_1d.gm9.data ...
-    - (-1.) .* 2.*pi .* profiles_1d.gm1.data ./ profiles_1d.gm9.data .* profiles_1d.f_df_dpsi.data ./ mu0;
-%
-% <j.B> = - sigma_Bp (2pi)^e_Bp dp/dpsi F - sigma_Bp F dF/dpsi / mu0 [ (2pi)^e_Bp F <1/R^2> + 1/(2pi)^e_Bp * <|grad psi|^2/R^2> / F ]
-% simaBp=-1 and eBp=1 for cocos=17 from TCV LIUQE
-%
-j_par = - (-1.) .* 2*pi .* profiles_1d.dpressure_dpsi.data .* profiles_1d.f.data ...
-        - (-1.) .* profiles_1d.f_df_dpsi.data ./ mu0 .* ...
-        ( (2.*pi) .* profiles_1d.f.data .* profiles_1d.gm1.data + 1./(2.*pi) .* gradpsi_over_R_sq ./ profiles_1d.f.data);
-profiles_1d.j_parallel.data = j_par./repmat(ids_summary.vacuum_toroidal_field.b0',size(profiles_1d.f.data,1),1);
-
-profiles_1d_fieldnames = fieldnames(profiles_1d);
-special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
-for it=1:ntime
-  for i=1:numel(profiles_1d_fieldnames)
-    if ~any(strcmp(profiles_1d_fieldnames{i},special_fields))
-      if ~isstruct(ids_summary.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}))
-        if ~ischar(profiles_1d.(profiles_1d_fieldnames{i}).data) && ~isempty(profiles_1d.(profiles_1d_fieldnames{i}).data) ...
-              && size(profiles_1d.(profiles_1d_fieldnames{i}).data,2)>=it
-          ids_summary.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}) = ...
-              profiles_1d.(profiles_1d_fieldnames{i}).data(:,it);
-        end
-      else
-        special_fields{end+1} = profiles_1d_fieldnames{i};
-      end
-    end
-  end
-end
-
-% special cases
-for it=1:ntime
-  ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_field_tor = ids_summary.time_slice{it}.profiles_1d.f(1) ...
-      ./ids_summary.time_slice{it}.global_quantities.magnetic_axis.r;
-  ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_tor = ids_summary.time_slice{it}.global_quantities.magnetic_axis.b_field_tor;
-  ids_summary.time_slice{it}.profiles_1d.rho_tor_norm = ids_summary.time_slice{it}.profiles_1d.rho_tor./ ...
-      ids_summary.time_slice{it}.profiles_1d.rho_tor(end);
-  ids_summary.time_slice{it}.profiles_1d.psi = ids_summary.time_slice{it}.global_quantities.psi_axis + ...
-      profiles_1d.rho_tor.x.^2 .* ...
-      (ids_summary.time_slice{it}.global_quantities.psi_boundary- ids_summary.time_slice{it}.global_quantities.psi_axis);
-end
-
-%
-%% profiles_2d{1} ala eqdsk, only this one thus grid_type=1
-%
-% b_field_r = gdat(params_summary.shot,'b_field_r','machine',gdat_params.machine);
-% b_field_tor = gdat(params_summary.shot,'b_field_tor','machine',gdat_params.machine);
-% b_field_z = gdat(params_summary.shot,'b_field_z','machine',gdat_params.machine);
-% b_r = gdat(params_summary.shot,'b_r','machine',gdat_params.machine);
-% b_tor = gdat(params_summary.shot,'b_tor','machine',gdat_params.machine);
-% b_z = gdat(params_summary.shot,'b_z','machine',gdat_params.machine);
-% grid = gdat(params_summary.shot,'grid','machine',gdat_params.machine); % special
-profiles_2d.grid_type.name = 'rectangular';
-profiles_2d.grid_type.index = 1;
-profiles_2d.grid_type.description = 'Cylindrical R,Z ala eqdsk';
-% j_parallel = gdat(params_summary.shot,'j_parallel','machine',gdat_params.machine);
-% j_tor = gdat(params_summary.shot,'j_tor','machine',gdat_params.machine);
-% phi = gdat(params_summary.shot,'phi','machine',gdat_params.machine);
-params_eff.data_request = 'psi';
-profiles_2d.psi = gdat(params_summary.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_summary.shot,'r','machine',gdat_params.machine); % not to be filled since in grid.dim1
-% theta = gdat(params_summary.shot,'theta','machine',gdat_params.machine);
-% z = gdat(params_summary.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
-for it=1:ntime
-  for i=1:numel(profiles_2d_fieldnames)
-    if ~any(strcmp(profiles_2d_fieldnames{i},special_fields))
-      if ~isstruct(ids_summary.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}))
-        if ~ischar(profiles_2d.(profiles_2d_fieldnames{i}).data) && ~isempty(profiles_2d.(profiles_2d_fieldnames{i}).data) ...
-              && size(profiles_2d.(profiles_2d_fieldnames{i}).data,3)>=it
-          ids_summary.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}) = ...
-              profiles_2d.(profiles_2d_fieldnames{i}).data(:,:,it);
-        end
-      else
-        special_fields{end+1} = profiles_2d_fieldnames{i};
-      end
-    end
-  end
-end
-
-% special cases
-for it=1:ntime
-  ids_summary.time_slice{it}.profiles_2d{1}.grid_type.name = profiles_2d.grid_type.name;
-  ids_summary.time_slice{it}.profiles_2d{1}.grid_type.index = profiles_2d.grid_type.index;
-  ids_summary.time_slice{it}.profiles_2d{1}.grid_type.description = profiles_2d.grid_type.description;
-  ids_summary.time_slice{it}.profiles_2d{1}.grid.dim1 = profiles_2d.psi.dim{1};
-  ids_summary.time_slice{it}.profiles_2d{1}.grid.dim2 = profiles_2d.psi.dim{2};
-  ids_summary.time_slice{it}.profiles_2d{1}.psi(:,:) = ids_summary.time_slice{it}.profiles_2d{1}.psi(:,:) + ...
-      global_quantities.psi_boundary.data(it);
-end
-
-% make arrays not filled in empty:
-ids_summary.grids_ggd = {};
-for it=1:numel(ids_summary.time_slice)
-  ids_summary.time_slice{it}.ggd = {};
-  ids_summary.time_slice{it}.boundary.strike_point = {};
-  ids_summary.time_slice{it}.boundary_separatrix.x_point = {};
-  ids_summary.time_slice{it}.boundary_separatrix.strike_point = {};
-end
-
-% special test matrix cocos transform
-% $$$ ldim1=129;
-% $$$ ldim2=257;
-% $$$ it=1;
-% $$$ ids_summary.time_slice{it}.coordinate_system.grid_type.index = 13;
-% $$$ ids_summary.time_slice{it}.coordinate_system.grid.dim1 = linspace(0,1,ldim1)';
-% $$$ ids_summary.time_slice{it}.coordinate_system.grid.dim2 = linspace(0,2*pi,ldim2);
-% $$$ ids_summary.time_slice{it}.coordinate_system.tensor_contravariant = 2.*ones(ldim1,ldim2,3,3);
-% $$$ ids_summary.time_slice{it}.coordinate_system.tensor_covariant = 0.5*ones(ldim1,ldim2,3,3);
-% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant = 13.*ones(ldim1,ldim2,3,3);
-% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant_error_upper = 14.*ones(ldim1,ldim2,3,3);
-% $$$ ids_summary.time_slice{it}.coordinate_system.g13_contravariant_error_lower = 12.*ones(ldim1,ldim2,3,3);
-% $$$ for it=1:2100
-% $$$   ids_summary.time_slice{it}.coordinate_system.g11_contravariant = 11.*ones(ldim1,ldim2,3,3);
-% $$$   ids_summary.time_slice{it}.coordinate_system.tensor_covariant = 0.5*ones(ldim1,ldim2,3,3);
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid_type.name = profiles_2d.grid_type.name;
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid_type.index = 11;
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid_type.description = profiles_2d.grid_type.description;
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1 = linspace(0,1,ldim1)';
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1_error_upper = 1.2.*linspace(0,1,ldim1)';
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid.dim1_error_lower = 0.8.*linspace(0,1,ldim1)';
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.grid.dim2 = linspace(0,2*pi,ldim2);
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.psi(:,:) = 11.*ones(ldim1,ldim2);
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.psi_error_upper(:,:) = 12.*ones(ldim1,ldim2);
-% $$$   ids_summary.time_slice{it}.profiles_2d{2}.psi_error_lower(:,:) = 10.*ones(ldim1,ldim2);
-% $$$ end
 
 % cocos automatic transform
 if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
-- 
GitLab