matlab/TCV/gdat_tcv.m

@@ -1639,7 +1639,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
 
           % power deposition related:
           ec_data.p_abs_plasma.data = pabs_gyro.data * 1e6;
-          ec_data.p_abs_plasma.data(end+1,:) = sum(ec_data.p_abs_plasma.data,1,'omitnan');
+          ec_data.p_abs_plasma.data(end+1,:) = sum(ec_data.p_abs_plasma.data,1,'omitnan');  % add total
           ec_data.p_abs_plasma.label = [strrep(pabs_gyro.comment,'MW','W') ' ; last index is total'];
           ec_data.p_abs_plasma.units = 'W';
           ec_data.p_abs_plasma.x = launchers_grid;
@@ -1648,21 +1648,21 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           ec_data.p_abs_plasma.dimunits = {launchers_label, 's'};
           %
           ec_data.p_dens.data = pow_dens.data * 1e6;
-          ec_data.p_dens.data(:,end+1,:) = sum(ec_data.p_dens.data,2,'omitnan');
+          ec_data.p_dens.data(:,end+1,:) = sum(ec_data.p_dens.data,2,'omitnan'); % add total
           ec_data.p_dens.label = [strrep(pow_dens.comment,'MW','W') ' ; last index is total'];
           ec_data.p_dens.units = 'W/m^3';
-          ec_data.p_dens.x = pow_dens.rgrid';
-          ec_data.p_dens.rhotor_norm = ec_data.p_dens.x;
+          ec_data.p_dens.x = pow_dens.rgrid;
+          ec_data.p_dens.grids = pow_dens.grids;
           ec_data.p_dens.t = pow_dens.tgrid;
           ec_data.p_dens.dim = {ec_data.p_dens.x, launchers_grid, ec_data.p_dens.t};
           ec_data.p_dens.dimunits = {'rhotor_norm', launchers_label, 's'};
           %
           ec_data.p_integrated.data = power_integrated.data * 1e6;
-          ec_data.p_integrated.data(:,end+1,:) = sum(ec_data.p_integrated.data,2,'omitnan');
+          ec_data.p_integrated.data(:,end+1,:) = sum(ec_data.p_integrated.data,2,'omitnan');  % add total
           ec_data.p_integrated.label = [strrep(power_integrated.comment,'MW','W') ' ; last index is total'];
           ec_data.p_integrated.units = 'W';
-          ec_data.p_integrated.x = power_integrated.rgrid';
-          ec_data.p_integrated.rhotor_norm = ec_data.p_integrated.x;
+          ec_data.p_integrated.x = power_integrated.rgrid;
+          ec_data.p_integrated.grids = power_integrated.grids;
           ec_data.p_integrated.t = power_integrated.tgrid;
           ec_data.p_integrated.dim = {ec_data.p_integrated.x, launchers_grid, ec_data.p_integrated.t};
           ec_data.p_integrated.dimunits = {'rhotor_norm', launchers_label, 's'};
@@ -1693,7 +1693,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
 
           % current drive deposition related:
           ec_data.cd_tot.data = icdtot.data * 1e6;
-          ec_data.cd_tot.data(end+1,:) = sum(ec_data.cd_tot.data,1,'omitnan');
+          ec_data.cd_tot.data(end+1,:) = sum(ec_data.cd_tot.data,1,'omitnan'); % add total
           ec_data.cd_tot.label = [strrep(icdtot.comment,'MA','A') ' ; last index is total'];
           ec_data.cd_tot.units = 'A';
           ec_data.cd_tot.x = launchers_grid;
@@ -1702,28 +1702,28 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           ec_data.cd_tot.dimunits = {launchers_label, 's'};
           %
           ec_data.cd_dens.data = currentdrive_dens.data * 1e6;
-          ec_data.cd_dens.data(:,end+1,:) = sum(ec_data.cd_dens.data,2,'omitnan');
+          ec_data.cd_dens.data(:,end+1,:) = sum(ec_data.cd_dens.data,2,'omitnan');  % add total
           ec_data.cd_dens.label = [strrep(currentdrive_dens.comment,'MA','A') ' ; last index is total'];
-          ec_data.cd_dens.units = 'A/m^2';
-          ec_data.cd_dens.x = currentdrive_dens.rgrid';
-          ec_data.cd_dens.rhotor_norm = ec_data.cd_dens.x;
+          ec_data.cd_dens.units = 'A/m^3';
+          ec_data.cd_dens.x = currentdrive_dens.rgrid;
+          ec_data.cd_dens.grids = currentdrive_dens.grids;
           ec_data.cd_dens.t = currentdrive_dens.tgrid;
           ec_data.cd_dens.dim = {ec_data.cd_dens.x, launchers_grid, ec_data.cd_dens.t};
           ec_data.cd_dens.dimunits = {'rhotor_norm', launchers_label, 's'};
           %
           ec_data.cd_integrated.data = currentdrive_integrated.data * 1e6;
-          ec_data.cd_integrated.data(:,end+1,:) = sum(ec_data.cd_integrated.data,2,'omitnan');
+          ec_data.cd_integrated.data(:,end+1,:) = sum(ec_data.cd_integrated.data,2,'omitnan'); % add total
           ec_data.cd_integrated.label = [strrep(currentdrive_integrated.comment,'MA','A') ' ; last index is total'];
           ec_data.cd_integrated.units = 'A';
-          ec_data.cd_integrated.x = currentdrive_integrated.rgrid';
-          ec_data.cd_integrated.rhotor_norm = ec_data.cd_integrated.x;
+          ec_data.cd_integrated.x = currentdrive_integrated.rgrid;
+          ec_data.cd_integrated.grids = currentdrive_integrated.grids;
           ec_data.cd_integrated.t = currentdrive_integrated.tgrid;
           ec_data.cd_integrated.dim = {ec_data.cd_integrated.x, launchers_grid, ec_data.cd_integrated.t};
           ec_data.cd_integrated.dimunits = {'rhotor_norm', launchers_label, 's'};
           %
           ec_data.max_cd_dens.data = icdmax.data * 1e6;
           ec_data.max_cd_dens.label = strrep(icdmax.comment,'MA','A');
-          ec_data.max_cd_dens.units = 'A/m^2';
+          ec_data.max_cd_dens.units = 'A/m^3';
           ec_data.max_cd_dens.x = [];
           ec_data.max_cd_dens.t = icdmax.tgrid;
           ec_data.max_cd_dens.dim = {ec_data.max_cd_dens.t};
@@ -1746,8 +1746,9 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           ec_data.width_cd_dens.dimunits = {'s'};
           %
           ec_data.cd_dens_doublewidth.data = currentdrive_dens_w2.data * 1e6;
+          ec_data.cd_dens_doublewidth.data(:,end+1,:) = sum(ec_data.cd_dens_doublewidth.data,2,'omitnan');  % add total
           ec_data.cd_dens_doublewidth.label = [strrep(currentdrive_dens_w2.comment,'MA','A') ' ; last index is total'];
-          for subfields={'x','rhotor_norm','t','dim','dimunits','units'}
+          for subfields={'x','grids','t','dim','dimunits','units'}
             ec_data.cd_dens_doublewidth.(subfields{1}) = ec_data.cd_dens.(subfields{1});
           end
 
@@ -1755,7 +1756,6 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           filled_successfully = true; %set flag to true
         catch ME
           warning('Problem retrieving TORAY data. \nError message: %s',ME.message);
-          getReport(ME) % without ; to get output
           % try to identify cause of the error
           if ~check_nodes_filled(shot,'toray')
             if ~isempty(ec_inputs.launchers_active)
@@ -1764,6 +1764,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
               else
                 msg = 'write_pyro(shot) found active launchers, but TORAY nodes are not filled, check hldsi(shot), and (ask to) relaunch TORAY.';
               end
+            else
+              msg = 'Run gdat(shot,''ec_data'',''ec_inputs'',1) to check whether there were active launchers in this shot or not.';
             end
           elseif ~isempty(gdat_data.gdat_params.trialindx)
             msg = 'Is the trial index filled? Check TORAY nodes with hdlsi(shot).';
@@ -1866,7 +1868,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       data_fullpath = '';
       ohm_help = '';
       if strcmp(lower(source_icd.ohm),'ibs')
-        ohm_data.cd_tot = gdat_tcv(shot,'\results::ibs:iohm');
+        ohm_data.cd_tot = gdat_tcv(shot,'\results::ibs:iohm','trialindx',gdat_params.trialindx);
         if gdat_data.mapping_for.tcv.gdat_timedim ==2
           tgrid_to_change = {'ohm_data.cd_tot.t','ohm_data.cd_tot.dim{1}'};
           for i=1:length(tgrid_to_change)
@@ -1875,7 +1877,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         end
         ohm_data.cd_tot.data = ohm_data.cd_tot.data';
         ohm_data.cd_tot.units = strrep(ohm_data.cd_tot.units,'kA','A');
-        ohm_data.cd_dens = gdat_tcv(shot,'\results::ibs:johmav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
+        ohm_data.cd_dens = gdat_tcv(shot,'\results::ibs:johmav','trialindx',gdat_params.trialindx); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
         if gdat_data.mapping_for.tcv.gdat_timedim ==2
           tgrid_to_change = {'ohm_data.cd_dens.t','ohm_data.cd_dens.dim{2}'};
           for i=1:length(tgrid_to_change)
@@ -1913,7 +1915,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       data_fullpath = '';
       bs_help = '';
       if strcmp(lower(source_icd.bs),'ibs')
-        bs_data.cd_tot = gdat_tcv(shot,'\results::ibs:ibs');
+        bs_data.cd_tot = gdat_tcv(shot,'\results::ibs:ibs','trialindx',gdat_params.trialindx);
         if gdat_data.mapping_for.tcv.gdat_timedim ==2
           tgrid_to_change = {'bs_data.cd_tot.t','bs_data.cd_tot.dim{1}'};
           for i=1:length(tgrid_to_change)
@@ -1922,7 +1924,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         end
         bs_data.cd_tot.data = bs_data.cd_tot.data';
         bs_data.cd_tot.units = strrep(bs_data.cd_tot.units,'kA','A');
-        bs_data.cd_dens = gdat_tcv(shot,'\results::ibs:jbsav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
+        bs_data.cd_dens = gdat_tcv(shot,'\results::ibs:jbsav','trialindx',gdat_params.trialindx); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
         if gdat_data.mapping_for.tcv.gdat_timedim ==2
           tgrid_to_change = {'bs_data.cd_dens.t','bs_data.cd_dens.dim{2}'};
           for i=1:length(tgrid_to_change)

matlab/TCV/tcv_requests_mapping.m

@@ -473,6 +473,11 @@ switch lower(data_request)
                     'bb=interpos(aa.t,aa.data(ix,:),gdat_tmp.t,-100);' ...
                     'gdat_tmp.data = 4e-7*pi.*gdat_tmp.data.^2.*bb./1.22; gdat_tmp.label=''\tau_R'';gdat_tmp.gdat_request=''tau_r'';' ...
                     'gdat_tmp.units=''s'';gdat_tmp.help=''tauR=mu0 a^2 / 1.22 etaneo(rhopol=0.4)'''];
+ case 'tau_tot'
+  mapping.timedim = 1;
+  mapping.label = 'total energy confinment time';
+  mapping.method = 'tdi';
+  mapping.expression = '\results::conf:tau';
  case 'te'
   mapping.timedim = 2;
   mapping.label = 'Te';

matlab/TCV_IMAS/tcv_get_ids_core_profiles.m

-function [ids_core_profiles,ids_core_profiles_description,varargout] = tcv_get_ids_core_profiles(shot,ids_equil_empty, gdat_params,varargin)
+function [ids_core_profiles,ids_core_profiles_description,varargout] = ...
+  tcv_get_ids_core_profiles(shot,ids_equil_empty,gdat_params,varargin)
 %
-%  [ids_core_profiles,ids_core_profiles_description,varargout] = tcv_get_ids_core_profiles(shot,ids_equil_empty,varargin);
+% [ids_core_profiles,ids_core_profiles_description,varargout] = ...
+%     tcv_get_ids_core_profiles(shot,ids_equil_empty,gdat_params,varargin);
 %
 %
-% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar options
-%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call,
+%              in particular error_bar options
 %
 % error_bar: from gdat_params.error_bar
 %            'delta' (default): error_bar to be added inserted in "upper" only as mentioned in description
@@ -22,18 +24,22 @@ tens_time = -1;
 tens_rho = -0.1;
 
 if exist('gdat_params','var')
-  [ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles','gdat_params',gdat_params,varargin{:});
+  [ids_core_profiles, params_cores_profiles] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles','gdat_params',gdat_params,varargin{:});
 else
-  [ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',varargin{:});
+  [ids_core_profiles, params_cores_profiles] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',varargin{:});
   aa=gdat_tcv;
   gdat_params = aa.gdat_params; % to get default params
 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)
 
+% initialize description
 ids_core_profiles_description = [];
 
-% base all from times fro nete_rho (which should be conf by default)
+%% setup profiles_1d
+% base all from times for nete_rho (which should be conf by default)
 params_eff = params_eff_ref;
 params_eff.data_request='ne_rho'; params_eff.fit = 1;
 temp_1d.ne_rho = gdat(params_cores_profiles.shot,params_eff);
@@ -70,13 +76,11 @@ ids_core_profiles.profiles_1d{1}.ion{1}.state = {};
 ids_core_profiles.profiles_1d{1}.neutral{1}.state = {};
 ids_core_profiles.profiles_1d(1:length(ids_core_profiles.time)) = ids_core_profiles.profiles_1d(1);
 
-
 % 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
-%
+
+%% vacuum_toroidal_field and time, using homogeneous
 params_eff = params_eff_ref;
 params_eff.data_request='b0';
 vacuum_toroidal_field.b0=gdat(params_cores_profiles.shot,params_eff);
@@ -85,7 +89,10 @@ ids_core_profiles.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
 ids_core_profiles.vacuum_toroidal_field.b0 = interpos(63,vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data,ids_core_profiles.time,-1);
 ids_core_profiles_description.vacuum_toroidal_field = [vacuum_toroidal_field_desc.b0 ' on ids_core_profiles.time, with interpos(63)'];
 
-% 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
+%% 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
 
 params_eff.data_request='ip';
 global_quantities.ip = gdat(params_cores_profiles.shot,'ip','machine',machine);
@@ -141,25 +148,20 @@ for i=1:length(global_quantities_fieldnames_eff)
     end
   end
 end
-% special case
-
 
-%
 %% profiles_1d (cannot use eqdsk since not same radial mesh)
+
+% prepare area for ids_core_profiles.profiles_1d{*}.grid.area
 params_eff.data_request = 'area_rho';
 temp_1d.area = gdat(params_cores_profiles.shot,params_eff);
 temp_1d_desc.area = params_eff.data_request;
-for ir=1:length(temp_1d.area.x)
+for ir=1:length(temp_1d.area.x) % map tmp_1d.area to core_profiles.time
   area_cpt(ir,:) = interpos(temp_1d.area.t,temp_1d.area.data(ir,:),ids_core_profiles.time,tens_time);
 end
-params_eff.data_request = 'q_rho';
-temp_1d.q = gdat(params_cores_profiles.shot,params_eff);
-temp_1d_desc.q = params_eff.data_request;
-for ir=1:length(temp_1d.q.x)
-  q_cpt(ir,:) = interpos(temp_1d.q.t,temp_1d.q.data(ir,:),ids_core_profiles.time,tens_time);
-end
+
 it_thom = iround_os(temp_1d.te_rho.t,ids_core_profiles.time);
 for it=1:length(ids_core_profiles.time)
+  % fill grid
   ids_core_profiles.profiles_1d{it}.grid.rho_tor_norm = temp_1d.fit.te_rho.grids_1d.rhotornorm(:,it);
   ids_core_profiles.profiles_1d{it}.grid.rho_tor = temp_1d.fit.te_rho.grids_1d.rhotornorm(:,it) ...
       .* temp_1d.fit.te_rho.grids_1d.rhotor_edge(it);
@@ -171,7 +173,9 @@ for it=1:length(ids_core_profiles.time)
       .* temp_1d.fit.te_rho.grids_1d.volume_edge(it);
   ids_core_profiles.profiles_1d{it}.grid.area = interpos(temp_1d.area.x,area_cpt(:,it),temp_1d.fit.te_rho.grids_1d.rhopolnorm, ...
           tens_rho,[1 2],[0 area_cpt(end,it)]);
+  % fill time
   ids_core_profiles.profiles_1d{it}.time = ids_core_profiles.time(it);
+  % fill electrons struct
   ids_core_profiles.profiles_1d{it}.electrons.temperature = temp_1d.fit.te_rho.data(:,it);
   ids_core_profiles.profiles_1d{it}.electrons.temperature_fit.measured = temp_1d.te_rho.data(:,it_thom(it));
   ids_core_profiles.profiles_1d{it}.electrons.temperature_fit.time_measurement = temp_1d.te_rho.t(it_thom(it));
@@ -185,6 +189,7 @@ for it=1:length(ids_core_profiles.time)
   ids_core_profiles.profiles_1d{it}.electrons.density_fit.source = {'Thomson, interpos fit'};
   ids_core_profiles.profiles_1d{it}.electrons.pressure_thermal = 1.6022e-19.*ids_core_profiles.profiles_1d{it}.electrons.density_thermal ...
       .* ids_core_profiles.profiles_1d{it}.electrons.temperature;
+  % fill zeff
   ids_core_profiles.profiles_1d{it}.zeff = global_quantities.z_eff_resistive.data(it) .* ...
       ones(size(ids_core_profiles.profiles_1d{it}.electrons.density));
 end
@@ -231,7 +236,8 @@ switch error_bar
   error(['tcv_ids_bpol_loop: error_bar option not known: ' error_bar])
 end
 
-% ion, assume only D if no CXRS (need to ask how to check if H...)
+%% ion struct
+% assume only D if no CXRS (need to ask how to check if H...)
 params_eff_fit1.data_request = 'cxrs';
 temp_1d.cxrs_rho = gdat(params_cores_profiles.shot,params_eff_fit1);
 temp_1d_desc.cxrs_rho = params_eff_fit1.data_request;
@@ -316,9 +322,14 @@ if ~isempty(temp_1d.cxrs_rho.ti.fit.data)
   end
 end
 
-for ir=1:length(temp_1d.q.x)
+%% q profile and magnetic shear
+params_eff.data_request = 'q_rho';
+temp_1d.q = gdat(params_cores_profiles.shot,params_eff);
+temp_1d_desc.q = params_eff.data_request;
+for ir=1:length(temp_1d.q.x) % map tmp_1d.q to core_profiles.time
   q_cpt(ir,:) = interpos(temp_1d.q.t,temp_1d.q.data(ir,:),ids_core_profiles.time,tens_time);
 end
+
 for it=1:length(ids_core_profiles.time)
   ij=isfinite(q_cpt(:,it));
   if sum(ij) >= 5
@@ -335,9 +346,28 @@ for it=1:length(ids_core_profiles.time)
     ids_core_profiles.profiles_1d{it}.magnetic_shear = -9.e40;
   end
 end
+temp_1d_desc.magnetic_shear = 'from interpos with rhotor';
+
+%% Current densities
+
+for it=1:length(ids_core_profiles.time)
+  ids_core_profiles.profiles_1d{it}.j_bootstrap = ...
+    current_bootstrap.bs.bs_data.cd_dens.data(:,it);
+  temp_1d_desc.j_bootstrap = current_bootstrap.bs.bs_data.cd_dens.label;
+  ids_core_profiles.profiles_1d{it}.j_ohmic = ...
+    current_non_inductive.ohm.ohm_data.cd_dens.data(:,it);
+  temp_1d_desc.j_ohmic = current_non_inductive.ohm.ohm_data.cd_dens.label;
+end
+
+%% parallel conductivity (neoclassical)
+signeo = gdat(params_cores_profiles.shot,'\results::ibs:signeo');
+for it=1:length(ids_core_profiles.time)
+  ids_core_profiles.profiles_1d{it}.conductivity_parallel = signeo.data(:,it);
+  temp_1d_desc.conductivity_parallel = signeo.label;
+end
 
-ids_core_profiles_description.temp_1d_desc = temp_1d_desc;
-ids_core_profiles_description.profiles_1d.magnetic_shear = 'from interpos with rhotor';
+%% add descriptions for profiles_1d
+ids_core_profiles_description.profiles_1d = temp_1d_desc;
 
 if nargin <= 2
   ids_core_profiles.code.name = ['tcv_get_ids_core_profiles, within gdat, with shot= ' num2str(params_cores_profiles.shot) ];

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

+function [ids_core_sources,ids_core_sources_description,varargout] = ...
+  tcv_get_ids_core_sources(shot,ids_equil_empty,gdat_params,varargin)
+%
+% [ids_core_sources,ids_core_sources_description,varargout] = ...
+%     tcv_get_ids_core_sources(shot,ids_equil_empty,gdat_params,varargin);
+%
+%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call
+%
+
+machine = 'tcv';
+tens_time = -1;
+tens_rho = -0.1;
+
+if exist('gdat_params','var')
+  [ids_core_sources, params_core_sources] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_sources','gdat_params',gdat_params,...
+    'homogeneous_time',0,varargin{:});
+else
+  [ids_core_sources, params_core_sources] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_sources','homogeneous_time',0,varargin{:});
+  aa=gdat_tcv;
+  gdat_params = aa.gdat_params; % to get default params
+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)
+
+% initialize description
+ids_core_sources_description = [];
+
+%% fill vacuum_toroidal field
+params_eff = params_eff_ref; params_eff.data_request='b0';
+b0_gdat=gdat(params_core_sources.shot,params_eff);
+ids_core_sources.vacuum_toroidal_field.r0 = b0_gdat.r0;
+ids_core_sources.vacuum_toroidal_field.b0 = b0_gdat.data;
+ids_core_sources_description.vacuum_toroidal_field = b0_gdat.request_description;
+ids_core_sources.time = b0_gdat.t;
+%%
+% save source template to fill for differnet source options
+source_template = ids_core_sources.source{1};
+profiles_template = source_template.profiles_1d{1};
+globals_template = source_template.global_quantities{1};
+
+last_index = 0; % index to be incremented when sources are added
+
+%% get time grid and data from gdat
+params_eff = params_eff_ref; params_eff.data_request='ohm_data';
+ohm_gdat = gdat(params_core_sources.shot,params_eff);
+params_eff = params_eff_ref; params_eff.data_request='powers';
+powers_gdat = gdat(params_core_sources.shot,params_eff);
+
+%% load liuqe data from nodes for conversions
+if params_eff_ref.liuqe == 1
+  nodes_liuqe = 'equil';
+else
+  nodes_liuqe = ['equil',num2str(params_eff_ref.liuqe)];
+end
+if check_nodes_filled(params_core_sources.shot,nodes_liuqe)
+  params_eff = params_eff_ref;
+  % normalized sqrt poloidal flux including axis
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:rho';
+  rho_pol_norm_liu = gdat(params_core_sources.shot,params_eff);
+  % <1/R^2>
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
+  Rm2_fs_liu = gdat(params_core_sources.shot,params_eff);
+  % Volume
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:vol';
+  vol_liu= gdat(params_core_sources.shot,params_eff);
+  % T(\psi) = R*B_tor
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:rbtor_rho';
+  T_liu = gdat(params_core_sources.shot,params_eff);
+  % safety factor q = -dPhi/dPsi
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:q';
+  q_liu = gdat(params_core_sources.shot,params_eff);
+  % I_tor = I_p, total plasma current
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:i_pl';
+  Ip_liu = gdat(params_core_sources.shot,params_eff);
+  % define liuqe_time from Ip_liu
+  liuqe_time = Ip_liu.t;
+else
+  error('Liuqe nodes %s not filled. Contact O. Sauter.',nodes_liuqe)
+end
+%% initialize source from template
+% ohm
+ohm_tgrid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_tgrid);
+
+main_desc = 'Source from ohmic heating'; production = 'IBS nodes';
+id_ohm.description = sprintf('%s from %s',main_desc,production);
+id_ohm.index = 7; id_ohm.name = 'ohmic';
+ids_core_sources.source{last_index+1} = source_template;
+ids_core_sources.source{last_index+1}.identifier = id_ohm;
+ids_core_sources.source{last_index+1}.profiles_1d(1:ohm_n_t) = {profiles_template};
+ids_core_sources.source{last_index+1}.global_quantities(1:ohm_n_t) = {globals_template};
+ohm_data = ohm_gdat.ohm.ohm_data;
+
+% load LIUQE data to convert
+% \tilde{j}_// = <jdotB>/(R0<Bphi/R>) to j_//0 = <jdotB>/B0; Bphi=T(psi)/R
+% interpolate liuqe outputs on ohm_tgrid
+T      = interp1(liuqe_time,T_liu.data.',     ohm_tgrid)';
+Rm2_fs = interp1(liuqe_time,Rm2_fs_liu.data.',ohm_tgrid)';
+vol    = interp1(liuqe_time,vol_liu.data.',   ohm_tgrid)';
+% get vacuum field data from ids
+R0 = ids_core_sources.vacuum_toroidal_field.r0;
+B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ohm_tgrid);
+
+jpar_tilde_to_jpar_0 = R0*T.*Rm2_fs./B0';
+
+% map volume and conversion factor on rho_pol grid of ohm_data cd_dens
+jpar_tilde_to_jpar_0_mapped = ...
+  interp1(rho_pol_norm_liu.data,jpar_tilde_to_jpar_0,ohm_data.cd_dens.rhopol_norm);
+vol_mapped = interp1(rho_pol_norm_liu.data,vol,ohm_data.cd_dens.rhopol_norm);
+
+for ii = 1:ohm_n_t
+  % profiles_1d
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = ohm_tgrid(ii);
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_pol_norm = ...
+    ohm_data.cd_dens.rhopol_norm';
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_tor_norm = ...
+    ohm_data.cd_dens.rhotor_norm(:,ii)';
+  jpar_0_tmp = ohm_data.cd_dens.data(:,ii)'.*jpar_tilde_to_jpar_0_mapped(:,ii)';
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = jpar_0_tmp;
+  integrated_jpar_0_tmp = cumtrapz(vol_mapped(:,ii),jpar_0_tmp)/2/pi/R0;
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.current_parallel_inside = integrated_jpar_0_tmp;
+
+  % globals
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.time = ohm_tgrid(ii);
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii);
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
+end
+
+last_index = last_index+1;  % add if statement to only increment if ohmic source has been added
+
+%% bs
+params_eff = params_eff_ref; params_eff.data_request='bs_data';
+bs_gdat = gdat(params_core_sources.shot,params_eff); bs_data = bs_gdat.bs.bs_data;
+bs_tgrid = bs_gdat.bs.t; bs_n_t = numel(bs_tgrid);
+
+main_desc = 'Bootstrap current'; production = 'IBS nodes';
+id_bs.description = sprintf('%s from %s',main_desc,production);
+id_bs.index = 13; id_bs.name = 'bootstrap_current';
+ids_core_sources.source{last_index+1} = source_template;
+ids_core_sources.source{last_index+1}.identifier = id_bs;
+ids_core_sources.source{last_index+1}.profiles_1d(1:bs_n_t) = {profiles_template};
+ids_core_sources.source{last_index+1}.global_quantities(1:bs_n_t) = {globals_template};
+
+for ii = 1:bs_n_t
+  % profiles_1d
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = bs_tgrid(ii);
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_pol_norm = ...
+    bs_data.cd_dens.rhopol_norm';
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_tor_norm = ...
+    bs_data.cd_dens.rhotor_norm(:,ii)';
+  % use same conversion as for ohmic data from LIUQE
+  jpar_0_tmp = bs_data.cd_dens.data(:,ii)'.*jpar_tilde_to_jpar_0_mapped(:,ii)';
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = jpar_0_tmp;
+  integrated_jpar_0_tmp = cumtrapz(vol_mapped(:,ii),jpar_0_tmp)/2/pi/R0;
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.current_parallel_inside = integrated_jpar_0_tmp;
+
+  % globals
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.time = bs_tgrid(ii);
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
+end
+
+last_index = last_index+1;  % add if statement to only increment if bs source has been added
+
+%% ec
+params_eff = params_eff_ref;
+params_eff.data_request='ec_data';
+params_eff.ec_inputs = 1; % load EC input information
+ec_gdat = gdat(params_core_sources.shot,params_eff);
+
+if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
+
+  ec_data = ec_gdat.ec.ec_data; ec_inputs = ec_gdat.ec.ec_inputs;
+  % get tgrid from gdat ec_data
+  ec_tgrid_toray = ec_data.p_dens.t; nt_ec_toray = numel(ec_tgrid_toray);
+
+  % retrieve active launcher information from ec_inputs
+  nb_launchers = numel(ec_inputs.launchers_active.data);
+  active_launchers = find(ec_inputs.launchers_active.data==1)';
+
+  % find times of injected EC power to interpolate power & current densities
+  % on p_ec_injected tgrid
+  ec_powers_tgrid = powers_gdat.ec.t; %nt_ec_powers = numel(ec_powers_tgrid);
+  % find times where EC is on to define time grid with extra time slice just
+  % before/after  EC power and at start/end of shot
+  itime_ec = find(powers_gdat.ec.data(:,end)>0);
+  if ec_powers_tgrid(itime_ec(end))>=ohm_tgrid(end)
+    i_time_end = iround(ec_powers_tgrid,ohm_tgrid(end));
+    ec_tgrid_out = [ohm_tgrid(1),ec_powers_tgrid(itime_ec(1)-1:i_time_end)'];
+  else
+    ec_tgrid_out = [ohm_tgrid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_tgrid(end)];
+  end
+  nt_ec_out = numel(ec_tgrid_out);
+  p_ec_injected = interpos(ec_powers_tgrid,powers_gdat.ec.data(:,end),ec_tgrid_out);
+
+  % Setup source structs for active launchers from template
+  main_desc = 'Source from electron cyclotron heating and current drive';
+  production = 'TORAY';
+  id_ec.index = 3; id_ec.name = 'ec';
+  for i_lau = active_launchers
+    id_ec.description = sprintf('L%i/G%i, %s from %s double width CD profiles',i_lau,ec_inputs.gyro2launcher.data(i_lau),main_desc,production);
+    ids_core_sources.source{last_index+i_lau} = source_template;
+    ids_core_sources.source{last_index+i_lau}.identifier = id_ec;
+    ids_core_sources.source{last_index+i_lau}.profiles_1d(1:nt_ec_out) = {profiles_template};
+    ids_core_sources.source{last_index+i_lau}.global_quantities(1:nt_ec_out) = {globals_template};
+  end
+
+  % get vacuum field data from ids
+  R0 = ids_core_sources.vacuum_toroidal_field.r0;
+  B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ec_tgrid_toray);
+
+  % do conversion j_{V,TORAY} = dI_\phi/dV = (1/2\pi)<j_\phi/R> to j_//0 = <jdotB>/B0;
+  % via j_tilde// = <jdotB>/(R0 <Bphi/R>) = 2*pi/R0/Rm2_fs[j_{V,TORAY}-dT/dV/T*I_{phi,cd}]
+
+  % interpolate liuqe outputs on ec_tgrid_toray
+  T      = interp1(liuqe_time,T_liu.data.',     ec_tgrid_toray)';
+  Rm2_fs = interp1(liuqe_time,Rm2_fs_liu.data.',ec_tgrid_toray)';
+  vol    = interp1(liuqe_time,vol_liu.data.',   ec_tgrid_toray)';
+
+%  % alternative conversion, slightly mismatching with the conversion
+%  % above, to be clarified
+%  % using j_//,cd = <j_cd/B> B, thus j_// = 2\pi <B^2> / (B_0 <B_\phi/R>)
+%   mu0     = 4*pi*10^-7;
+%   q      = interp1(liuqe_time,q_liu.data.',     ec_tgrid_toray)';
+%   Ip     = interp1(liuqe_time,Ip_liu.data,      ec_tgrid_toray);
+%   Vprime = -2*pi*q./T./Rm2_fs;
+%   B2_fs   = -mu0*Ip./Vprime + T.^2.*Rm2_fs;
+%   jtoray_to_jpar0 = 2*pi./B0./T./Rm2_fs.*B2_fs;
+%   jtoray_to_jpar0_mapped = nan(size(ec_data.cd_dens.x));
+
+  % interpolate on ec_data rho_pol grid
+  vol_mapped    = nan(size(ec_data.cd_dens.x));
+  T_mapped      = nan(size(ec_data.cd_dens.x));
+  dTdV_mapped   = nan(size(ec_data.cd_dens.x));
+  Rm2_fs_mapped = nan(size(ec_data.cd_dens.x));
+  for ii = 1:nt_ec_toray
+%     jtoray_to_jpar0_mapped(:,ii) = ...
+%       interp1(rho_pol_norm_liu.data,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+    vol_mapped(:,ii)      = interpos(...
+      rho_pol_norm_liu.data,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+    T_mapped(:,ii)        = interpos(...
+      rho_pol_norm_liu.data,T(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+    [~,dTdV_mapped(:,ii)] = interpos(vol_mapped(:,ii),T_mapped(:,ii));
+    Rm2_fs_mapped(:,ii)   = interpos(...
+      rho_pol_norm_liu.data,Rm2_fs(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+  end
+
+  % load 1d_profiles from ec_data
+  p_dens        = ec_data.p_dens.data;
+  p_dens(isnan(p_dens)) = 0;
+  p_integrated  = ec_data.p_integrated.data;
+  p_integrated(isnan(p_integrated)) = 0;
+  cd_dens       = ec_data.cd_dens_doublewidth.data; % use double width to have realistic deposition broadening
+  cd_dens(isnan(cd_dens)) = 0;
+  cd_integrated = ec_data.cd_integrated.data;
+  cd_integrated(isnan(cd_integrated)) = 0;
+
+  % initialize variables
+  jpar_tilde       = zeros(size(cd_dens));
+  jpar0            = zeros(size(cd_dens));
+  jpar0_integrated = zeros(size(cd_dens));
+
+  % convert and integrate j_TORAY to j//0
+  for i_lau = active_launchers
+%     % alternative conversion
+%     jpar0(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
+
+    jpar_tilde(:,i_lau,:) = 2*pi/R0./Rm2_fs_mapped.*(squeeze(cd_dens(:,i_lau,:))...
+      -dTdV_mapped./T_mapped.*squeeze(cd_integrated(:,i_lau,:)));
+    jpar0(:,i_lau,:) = squeeze(jpar_tilde(:,i_lau,:)).*R0.*T_mapped.*Rm2_fs_mapped./B0;
+
+    for ii = 1:nt_ec_toray
+      % integrate j//0
+      [~,~,~,jpar0_integrated(:,i_lau,ii)] = interpos(vol_mapped(:,ii),jpar0(:,i_lau,ii)/(2*pi*R0));
+    end
+  end
+
+%% interpolating 1d_profiles from 'ec_tgrid_toray' grid (toray tgrid) to 'interp_tgrid' (powers tgrid when powers.data>0)
+  ij = iround_os(ec_tgrid_out,ec_tgrid_toray);
+  sparse_p_ec_injected = p_ec_injected(ij); % injected ec power vals corresponding to ec_tgrid_toray
+  n_rho = size(p_dens, 1);
+  %rho grids to be interpolated on power time grid
+  rho_pol_norm = ec_data.cd_dens.grids.rho_pol_norm;
+  rho_tor_norm = ec_data.cd_dens.grids.rho_tor_norm;
+
+  % calculate normalised profiles on ec_tgrid_toray grid
+  norm_p_dens           = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+  norm_p_integrated     = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+  norm_jpar0            = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+  norm_jpar0_integrated = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+  for it = 1:nt_ec_toray
+    norm_p_dens(:,:,it)           = p_dens(:,:,it)          ./sparse_p_ec_injected(it);
+    norm_p_integrated(:,:,it)     = p_integrated(:,:,it)    ./sparse_p_ec_injected(it);
+    norm_jpar0(:,:,it)            = jpar0(:,:,it)           ./sparse_p_ec_injected(it);
+    norm_jpar0_integrated(:,:,it) = jpar0_integrated(:,:,it)./sparse_p_ec_injected(it);
+  end
+
+  % interpolate normalised p_dens profiles on ec_powers_tgrid
+  interp_norm_p_dens           = zeros(n_rho,nb_launchers+1,nt_ec_out);
+  interp_norm_p_integrated     = zeros(n_rho,nb_launchers+1,nt_ec_out);
+  interp_norm_jpar0            = zeros(n_rho,nb_launchers+1,nt_ec_out);
+  interp_norm_jpar0_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
+  for i_lau = active_launchers
+    for irho = 1:n_rho
+      % get power and current density at each rho
+      trace_p_dens        = squeeze(norm_p_dens(irho,i_lau,:));
+      trace_p_integrated  = squeeze(norm_p_integrated(irho,i_lau,:));
+      trace_jpar0       = squeeze(norm_jpar0(irho,i_lau,:));
+      trace_jpar0_integrated = squeeze(norm_jpar0_integrated(irho,i_lau,:));
+      % interpolate on gdat powers tgrid, interpos 21 to extrapolate,
+      % taking constant value y_edge; needed for time window between start of EC and first TS time
+      interp_norm_p_dens(          irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_dens,          ec_tgrid_out(3:end-2));
+      interp_norm_p_integrated(    irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_integrated,    ec_tgrid_out(3:end-2));
+      interp_norm_jpar0(           irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_jpar0,           ec_tgrid_out(3:end-2));
+      interp_norm_jpar0_integrated(irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_jpar0_integrated,ec_tgrid_out(3:end-2));
+    end
+  end
+
+  % interpolate the rho_pol & rho_tor on ec_powers_tgrid
+  interp_rho_pol = zeros(n_rho,nt_ec_out);
+  interp_rho_tor = zeros(n_rho,nt_ec_out);
+  for irho = 1:n_rho
+    interp_rho_pol(irho,3:end-2) = interpos(21,ec_tgrid_toray,rho_pol_norm(irho,:),ec_tgrid_out(3:end-2));
+    interp_rho_tor(irho,3:end-2) = interpos(21,ec_tgrid_toray,rho_tor_norm(irho,:),ec_tgrid_out(3:end-2));
+  end
+  % fill rho_pol just outside of TORAY times(zero power) with the known profiles one time slice after/before
+  interp_rho_pol(:,2) = interp_rho_pol(:,3); interp_rho_pol(:,end-1) = interp_rho_pol(:,end-2);
+
+  % normalised & interpolated profiles * p_ec_injected on interp_tgrid
+  interp_p_dens        = zeros(n_rho,nb_launchers+1,nt_ec_out);
+  interp_p_integrated  = zeros(n_rho,nb_launchers+1,nt_ec_out);
+  interp_jpar0         = zeros(n_rho,nb_launchers+1,nt_ec_out);
+  interp_jpar0_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
+  for it = 1:nt_ec_out % fill profiles only for time slices when EC is nonzero, thus it+2 assignment
+    interp_p_dens(:,:,it)           = interp_norm_p_dens(:,:,it)          .*p_ec_injected(it);
+    interp_p_integrated(:,:,it)     = interp_norm_p_integrated(:,:,it)    .*p_ec_injected(it);
+    interp_jpar0(:,:,it)            = interp_norm_jpar0(:,:,it)           .*p_ec_injected(it);
+    interp_jpar0_integrated(:,:,it) = interp_norm_jpar0_integrated(:,:,it).*p_ec_injected(it);
+  end
+
+  % fill the ids entries
+  for i_lau = active_launchers
+    for ii = 1:nt_ec_out
+      % 1d_profiles
+      ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = ec_tgrid_out(ii);
+      % grids (do not fill for first and last time slice)
+      if ii~=1 || ii~=nt_ec_out
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = interp_rho_tor(:,ii);
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_pol_norm = interp_rho_pol(:,ii);
+        % power density
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.energy = ...
+          interp_p_dens(:,i_lau,ii);
+        % integrated power density
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.power_inside = ...
+          interp_p_integrated(:,i_lau,ii);
+        % current density
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ...
+          interp_jpar0(:,i_lau,ii);
+        % integrated current density
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
+          interp_jpar0_integrated(:,i_lau,ii);
+      end
+
+      % globals
+      ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.time = ec_tgrid_out(ii);
+      ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.power = ...
+        interp_p_integrated(end,i_lau,ii);
+      ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.current_parallel = ...
+        interp_jpar0_integrated(end,i_lau,ii);
+    end
+  end
+
+  %add empty sources for rest of unused launchers
+  if numel(ids_core_sources.source)-last_index ~= nb_launchers
+    ids_core_sources.source{last_index+nb_launchers} = [];
+  end
+
+  last_index = last_index+nb_launchers;
+end
+
+%% NBI. same tgrid for NBI1 and NBI2
+
+check_nbi = [~isempty(powers_gdat.nbi1.data),~isempty(powers_gdat.nbi2.data)];
+active_nbi = find(check_nbi==1);
+nb_nbi = numel(check_nbi);
+nbi_names = {'nbi1','nbi2'};
+
+if numel(active_nbi)>0
+  % get tgrid (same for NBI1 & 2 if both active) from active_nbi(1)
+  nbi_powers_tgrid = powers_gdat.(nbi_names{active_nbi(1)}).t;
+  % find times where NBI is on to define time grid with extra time slice just
+  % before/after  NBI power and at start/end of shot
+  itime_nbi = find((powers_gdat.(nbi_names{active_nbi(1)}).data>0));
+  if nbi_powers_tgrid(itime_nbi(end))>=ohm_tgrid(end)
+    i_time_end = iround(nbi_powers_tgrid,ohm_tgrid(end));
+    nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
+  else
+    nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_tgrid(end)];
+  end
+  nt_nbi_out = numel(nbi_tgrid_out);
+
+  % Setup source structs for active nbi from template
+  main_desc = 'Source from Neutral Beam Injection';
+  id_nbi.index = 2; id_nbi.name = 'nbi';
+  for i_nbi = active_nbi
+    id_nbi.description = sprintf('NBI%i %s',i_nbi,main_desc);
+    ids_core_sources.source{last_index+i_nbi} = source_template;
+    ids_core_sources.source{last_index+i_nbi}.identifier = id_nbi;
+    ids_core_sources.source{last_index+i_nbi}.profiles_1d(1:nt_nbi_out) = {profiles_template};
+    ids_core_sources.source{last_index+i_nbi}.global_quantities(1:nt_nbi_out) = {globals_template};
+  end
+
+  disp('Loading of current & power densities from ASTRA not implemented yet.')
+  disp('Checking if ASTRA run available on partition /Lac8_D:')
+  [~,hostname] = unix('hostname');
+  if strcmp(hostname,'lac8.epfl.ch')
+    unix(sprintf('ls /Lac8_D/ASTRA/ | grep ''%i'' && echo File for shotnumber exists! || echo File for shotnumber does not exist!',shot));
+  else
+    unix(sprintf('ssh $(whoami)@lac8 "ls /Lac8_D/ASTRA/ | grep ''%i'' && echo File for shotnumber exists! || echo File for shotnumber does not exist! && exit"',shot));
+  end
+
+  for i_nbi = active_nbi
+    p_nbi_injected_tmp = interpos(nbi_powers_tgrid,powers_gdat.(nbi_names{i_nbi}).data,nbi_tgrid_out);
+    for ii = 1:nt_nbi_out
+      % globals
+      ids_core_sources.source{last_index+i_nbi}.global_quantities{ii}.time = nbi_tgrid_out(ii);
+      ids_core_sources.source{last_index+i_nbi}.global_quantities{ii}.power = p_nbi_injected_tmp(ii);
+    end
+  end
+
+  % add empty for unused NBI
+  if numel(ids_core_sources.source)-last_index ~= nb_nbi
+    ids_core_sources.source{last_index+nb_nbi} = [];
+  end
+
+  last_index = last_index+nb_nbi;
+end
+
+%% DNBI has it's own time grid
+if ~isempty(powers_gdat.dnbi)
+  % get tgrid (same for NBI1 & 2 if both active) from active_nbi(1)
+  dnbi_powers_tgrid = powers_gdat.dnbi.t;
+  % find times where DNBI is on to define time grid with extra time slice just
+  % before/after  DNBI power and at start/end of shot
+  itime_dnbi = find((powers_gdat.dnbi.data>0));
+  if dnbi_powers_tgrid(itime_dnbi(end))>=ohm_tgrid(end)
+    i_time_end = iround(dnbi_powers_tgrid,ohm_tgrid(end));
+    dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:i_time_end)'];
+  else
+    dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:itime_dnbi(end)+1)',ohm_tgrid(end)];
+  end
+  nt_dnbi_out = numel(dnbi_tgrid_out);
+
+  id_dnbi.description = 'DNBI Source from Neutral Beam Injection';
+  id_dnbi.index = 2; id_dnbi.name = 'nbi';
+  ids_core_sources.source{last_index+1} = source_template;
+  ids_core_sources.source{last_index+1}.identifier = id_dnbi;
+  ids_core_sources.source{last_index+1}.profiles_1d(1:nt_dnbi_out) = {profiles_template};
+  ids_core_sources.source{last_index+1}.global_quantities(1:nt_dnbi_out) = {globals_template};
+
+  p_dnbi_injected = interpos(dnbi_powers_tgrid,powers_gdat.dnbi.data,dnbi_tgrid_out);
+  for ii = 1:nt_dnbi_out
+    % globals
+    ids_core_sources.source{last_index+1}.global_quantities{ii}.time = nbi_tgrid_out(ii);
+    ids_core_sources.source{last_index+1}.global_quantities{ii}.power = p_dnbi_injected(ii);
+  end
+
+  last_index = last_index+1;
+end
+
+%% total
+% id_total.description = 'Total source; combines all sources';
+% id_total.index = 1; id_total.name = 'total';
+% ids_core_sources.source{} = source_template;
+% ids_core_sources.source{}.identifier = id_total;
+% ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
+% ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
+
+%% add descriptions for profiles_1d
+
+desc.source = ...
+  'available by now {source_index}: ohmic {1}, bootstrap {2}, ec {3+} for individual launchers (if present in shot), nbi {3+|14+} if available, dnbi {3|14|16} if available';
+desc.profiles_1d.electrons.energy = 'absorbed power density';
+desc.profiles_1d.electrons.power_inside = 'integrated absorbed power density';
+desc.profiles_1d.j_parallel = 'parallel current density';
+desc.profiles_1d.current_parallel_inside = 'integrated parallel current density';
+desc.globals.power = 'total power coupled to the plasma';
+desc.globals.current_parallel = 'total parallel current driven';
+ids_core_sources_description = desc;
+
+%%
+if nargin <= 2
+  ids_core_sources.code.name = ['tcv_get_ids_core_sources, within gdat, with shot= ' num2str(params_core_sources.shot) ];
+else
+  ids_core_sources.code.name = ['tcv_get_ids_core_sources, within gdat, with shot= ' num2str(params_core_sources.shot) '; varargin: ' varargin{:}];
+end
+ids_core_sources_description.code.name = ids_core_sources.code.name;
+
+ids_core_sources.code.output_flag = zeros(size(ids_core_sources.time));
+
+% cocos automatic transform
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
+  [ids_core_sources,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_core_sources,'core_sources',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

matlab/TCV_IMAS/tcv_get_ids_core_transport.m

+function [ids_core_transport,ids_core_transport_description,varargout] = ...
+  tcv_get_ids_core_transport(shot,ids_equil_empty,gdat_params,varargin)
+%
+% [ids_core_transport,ids_core_transport_description,varargout] = ...
+%     tcv_get_ids_core_transport(shot,ids_equil_empty,gdat_params,varargin);
+%
+%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call
+%
+
+machine = 'tcv';
+tens_time = -1;
+tens_rho = -0.1;
+
+if exist('gdat_params','var')
+  [ids_core_transport, params_core_transport] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles','gdat_params',gdat_params,varargin{:});
+else
+  [ids_core_transport, params_core_transport] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',varargin{:});
+  aa=gdat_tcv;
+  gdat_params = aa.gdat_params; % to get default params
+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)
+
+% initialize description
+ids_core_transport_description = [];
+
+%%
+last_index = 0;
+
+% fill model [name: transport solver, desc: output from transport solver]
+% setup model and profiles_1d
+comment = 'Output from a transport solver';
+ids_core_transport.model{last_index+1}.comment = comment;
+ids_core_transport.model{last_index+1}.name = comment;
+ids_core_transport.model{last_index+1}.identifier.index = 2;
+ids_core_transport.model{last_index+1}.name = 'transport_solver';
+
+% read data and setup time
+params_eff.data_request='\results::conf:chie';
+temp_1d.chie = gdat(params_core_transport.shot,params_eff.data_request);
+temp_1d_desc.chie = params_eff.data_request;
+
+if isempty(temp_1d.chie.t)
+  warning('no data in chie node, might need to rerun anaprofs')
+  return
+end
+
+ids_core_transport.time = temp_1d.chie.t;
+ids_core_transport_description.time = ['from node' params_eff.data_request];
+
+ids_core_transport.model{last_index+1}.profiles_1d(1:length(ids_core_transport.time)) = ...
+  ids_core_transport.model{last_index+1}.profiles_1d(1);
+
+% fill profiles_1d
+for it=1:length(ids_core_transport.time)
+  ids_core_transport.model{last_index+1}.profiles_1d{it}.time = ids_core_transport.time(it);
+  ids_core_transport.model{last_index+1}.profiles_1d{it}.electrons.energy.d = temp_1d.chie.data(:,it);
+  temp_1d_desc.electrons.energy.d = temp_1d.chie.label;
+end
+
+%% add descriptions for profiles_1d
+ids_core_transport_description.profiles_1d = temp_1d_desc;
+
+%%
+if nargin <= 2
+  ids_core_transport.code.name = ['tcv_get_ids_core_transport, within gdat, with shot= ' num2str(params_core_transport.shot) ];
+else
+  ids_core_transport.code.name = ['tcv_get_ids_core_transport, within gdat, with shot= ' num2str(params_core_transport.shot) '; varargin: ' varargin{:}];
+end
+ids_core_transport_description.code.name = ids_core_transport.code.name;
+
+ids_core_transport.code.output_flag = zeros(size(ids_core_transport.time));
+
+% cocos automatic transform
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
+  [ids_core_transport,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_core_transport,'core_transport',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

matlab/TCV_IMAS/tcv_get_ids_equilibrium.m

-function [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,ids_equil_empty, gdat_params,varargin)
+function [ids_equilibrium,ids_equilibrium_description,varargout] = ...
+  tcv_get_ids_equilibrium(shot,ids_equil_empty,gdat_params,varargin)
 %
-%  [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,ids_equil_empty,varargin);
+% [ids_equilibrium,ids_equilibrium_description,varargout] = ...
+%     tcv_get_ids_equilibrium(shot,ids_equil_empty,gdat_params,varargin);
 %
 %
-% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar and cocos_out options
+% 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_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium','gdat_params',gdat_params,varargin{:});
+  [ids_equilibrium, params_equilibrium] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'equilibrium','gdat_params',gdat_params,varargin{:});
 else
-  [ids_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium',varargin{:});
+  [ids_equilibrium, params_equilibrium] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'equilibrium',varargin{:});
   aa=gdat_tcv;
   gdat_params = aa.gdat_params; % to get default params
 end

matlab/TCV_IMAS/tcv_get_ids_summary.m

@@ -16,7 +16,7 @@ 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
+% "global_quantities_desc" which contains the same subfields themselves containing the gdat string after shot used
 %
 % vacuum_toroidal_field and time, using homogeneous
 %
@@ -36,6 +36,7 @@ ids_summary.code.name = 'gdat';
 ids_summary.code.version = git_release_hash;
 ids_summary.code.repository = git_url;
 
+%% boundary
 % use equilibrium time as reference, so load first elongation(time)
 
 params_eff = params_eff_ref;
@@ -59,7 +60,6 @@ ids_summary.time = unique([t_before; [kappa.t(1):max((kappa.t(end)-kappa.t(1))/(
 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:
 boundary.elongation = kappa;
 boundary_desc.elongation = kappa.gdat_params.data_request;
 params_eff.data_request = 'r_geom';
@@ -121,13 +121,13 @@ if ~isempty(special_fields)
 % $$$   end
 end
 
-% disruption, use specific function
+%% 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);
 
-% gas: bottom/top?, deuterium, helium_xx, ?
+%% 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;
@@ -148,7 +148,8 @@ for i=1:numel(gas_injection_rates_fieldnames)
   end
 end
 %params_eff.doplot=1;
-% global_quantities:
+
+%% 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];
@@ -182,6 +183,9 @@ global_quantities_desc.current_non_inductive = [params_eff.data_request ' ; ' gl
 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 = '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 = '\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];
@@ -263,6 +267,7 @@ end
 ids_summary.global_quantities.r0.value = global_quantities.r0.data;
 ids_summary.global_quantities.r0.source = ['gdat request: ' global_quantities_desc.r0];
 
+%% volume average
 % cst zeff
 params_eff.data_request = '\tcv_shot::top.results.ibs:z_eff';
 volume_average.zeff = gdat(params_summary.shot,params_eff);
@@ -283,6 +288,29 @@ for i=1:numel(volume_average_fieldnames)
   end
 end
 
+%% line average
+params_eff.data_request = 'nel';
+line_average.n_e = gdat(params_summary.shot,params_eff);
+line_average_desc.n_e = [params_eff.data_request ' ; ' volume_average.zeff.data_fullpath ' (not vol avrg but to match Iohm)'];
+
+special_fields = {}; % fields needing non-automatic treatments
+line_average_fieldnames = fieldnames(line_average);
+for i=1:numel(line_average_fieldnames)
+  if ~any(strcmp(line_average_fieldnames{i},special_fields))
+    if ~isstruct(ids_summary.line_average.(line_average_fieldnames{i}).value)
+      ids_summary.line_average.(line_average_fieldnames{i}).value = interp1( ...
+          line_average.(line_average_fieldnames{i}).t,line_average.(line_average_fieldnames{i}).data, ...
+          ids_summary.time,'linear',NaN);
+      ids_summary.line_average.(line_average_fieldnames{i}).source = ['gdat request: ' line_average_desc.(line_average_fieldnames{i})];
+    else
+      special_fields{end+1} = line_average_fieldnames{i};
+    end
+  end
+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, ...

matlab/TCV_IMAS/tcv_ids_headpart.m

-function [ids_generic, params_ids_generic] = tcv_ids_headpart(shot,ids_in,ids_name,varargin);
+function [ids_generic, params_ids_generic] = tcv_ids_headpart(shot,ids_in,ids_name,varargin)
 %
 % [ids_generic, params_ids_generic] = tcv_ids_headpart(shot,ids_in,ids_name,varargin);
 %
 % parses inputs and fill in ids_properties
 %
-%
 % varargin options:
-%
-% 'comment': comment to include in ids_properties, using gdat_params for example cocos_in and cocos_out
-% 'homogeneous_time': homogeneous_time in ids_properties: 1 (default) if the whole ids has same time, 0 otherwise
-% 'gdat_params': gdat params structure
+% 'comment':          comment to include in ids_properties, using gdat_params 
+%                     for example cocos_in and cocos_out
+% 'homogeneous_time': homogeneous_time in ids_properties: 
+%                     1 (default) if the whole ids has same time, 0 otherwise
+% 'gdat_params':      gdat params structure
 %
 %
 % example:
 %   [ids_equilibrium, params_ids_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium','comment','your comment');
 %
-%
 
 % initialize input parser
 p = inputParser;