diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 4b443b3bb37f3f3fb908e8d7b9fad714eae9d104..8ce621c5708f9a2b75df0b45afb5920ee67795b5 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/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)
diff --git a/matlab/TCV/tcv_requests_mapping.m b/matlab/TCV/tcv_requests_mapping.m
index 04dcc1f5c8e483714b1882bb15b8e965bb07ed20..824348abb9c6997c13fe4206e93a67bdc80ab41f 100644
--- a/matlab/TCV/tcv_requests_mapping.m
+++ b/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';
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
index 2891e6dd8eb526d7740d301ba93a1edcc597bd27..a230c9351cf0fb4a1dadfbec902eb2903aed5d5f 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
@@ -1,10 +1,12 @@
-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) ];
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
new file mode 100644
index 0000000000000000000000000000000000000000..6f630945d23ad81bdf831e41bc201d4ec6e3602f
--- /dev/null
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -0,0 +1,508 @@
+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
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_transport.m b/matlab/TCV_IMAS/tcv_get_ids_core_transport.m
new file mode 100644
index 0000000000000000000000000000000000000000..89421f6b67f9dc1df5be26547c118a17c3614fe0
--- /dev/null
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_transport.m
@@ -0,0 +1,83 @@
+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
diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index 263876023421dd6f1178f48302443aa01b5f807b..b0b22e2d593c954f9091642d07135f68c40c2588 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -1,15 +1,20 @@
-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
diff --git a/matlab/TCV_IMAS/tcv_get_ids_summary.m b/matlab/TCV_IMAS/tcv_get_ids_summary.m
index 2a720beba4e976d3ab1443119ab95159cb7f7bf3..e24b4ee20f0e52a73dd2b1751032eda58bc3efc7 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_summary.m
+++ b/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, ...
diff --git a/matlab/TCV_IMAS/tcv_ids_headpart.m b/matlab/TCV_IMAS/tcv_ids_headpart.m
index fa90b8cf836577146f499ba9cbe0ff7e2a53d311..53216596c61a9bbccaa9c08fc501e892ad01a1af 100644
--- a/matlab/TCV_IMAS/tcv_ids_headpart.m
+++ b/matlab/TCV_IMAS/tcv_ids_headpart.m
@@ -1,21 +1,20 @@
-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;