fixing indexing of sources

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

@@ -32,106 +32,158 @@ ids_core_sources_description = [];
 source_template = ids_core_sources.source{1};
 profiles_template = source_template.profiles_1d{1};
 globals_template = source_template.global_quantities{1};
-%% get time grid and data from gdat & initialize source from template
+%% get time grid and data from gdat 
 ohm_gdat = gdat(shot,'ohm_data');
+powers_gdat = gdat(shot,'powers');
 % fix time grid to ohm time grid, same as bs
 t_grid = ohm_gdat.t; n_t = numel(t_grid); t_spacing = mean(diff(t_grid));
 
+last_index = 0;
+%% initialize source from template
 % ohm
 id_ohm.description = 'Source from ohmic heating';
 id_ohm.index = 7; id_ohm.name = 'ohmic';
-ids_core_sources.source{7} = source_template;
-ids_core_sources.source{7}.identifier = id_ohm;
-ids_core_sources.source{7}.profiles_1d(1:n_t) = {profiles_template};
-ids_core_sources.source{7}.global_quantities(1:n_t) = {globals_template};
+ids_core_sources.source{1} = source_template;
+ids_core_sources.source{1}.identifier = id_ohm;
+ids_core_sources.source{1}.profiles_1d(1:n_t) = {profiles_template};
+ids_core_sources.source{1}.global_quantities(1:n_t) = {globals_template};
 ohm_data = ohm_gdat.ohm.ohm_data;
 
-% ec
-id_ec.description = 'Sources from electron cyclotron heating and current drive';
-id_ec.index = 3; id_ec.name = 'ec';
-ids_core_sources.source{3} = source_template;
-ids_core_sources.source{3}.identifier = id_ec;
-ids_core_sources.source{3}.profiles_1d(1:n_t) = {profiles_template};
-ids_core_sources.source{3}.global_quantities(1:n_t) = {globals_template};
+%fill profiles for times n t_grid
+for ii = 1:n_t
+  % ohm
+  % 1d_profiles
+  ids_core_sources.source{1}.profiles_1d{ii}.time = t_grid(ii);
+  ids_core_sources.source{1}.profiles_1d{ii}.j_parallel = ...
+    ohm_data.cd_dens.data(:,ii)';
+  % globals
+  ids_core_sources.source{1}.global_quantities{ii}.time = t_grid(ii);
+  ids_core_sources.source{1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii);
+end
+last_index = last_index+1;
+%% bs
+id_bs.description = 'Bootstrap current';
+id_bs.index = 13; id_bs.name = 'bootstrap_cuurent';
+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:n_t) = {profiles_template};
+ids_core_sources.source{last_index+1}.global_quantities(1:n_t) = {globals_template};
+% load data
+bs_gdat = gdat(shot,'bs_data'); bs_data = bs_gdat.bs.bs_data;
+
+% fill profiles for times n t_grid
+for ii = 1:n_t
+   ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = ...
+    bs_data.cd_dens.data(:,ii)'; 
+end
+last_index = last_index+1;
+%% ec
 % load data
-%nb_launchers = size(aa.ec.ec_data.rho_max_pow_dens.data,1);
 ec_gdat = gdat(shot,'ec_data'); ec_data = ec_gdat.ec.ec_data;
 ec_time = ec_data.cd_dens.t;
+
+nb_launchers = size(ec_gdat.ec.ec_data.rho_max_pow_dens.data,1);
+id_ec.description = 'Sources from electron cyclotron heating and current drive';
+id_ec.index = 3; id_ec.name = 'ec';
+for i = 1:nb_launchers
+    ids_core_sources.source{last_index+i} = source_template;
+    ids_core_sources.source{last_index+i}.identifier = id_ec;
+    ids_core_sources.source{last_index+i}.profiles_1d(1:n_t) = {profiles_template};
+    ids_core_sources.source{last_index+i}.global_quantities(1:n_t) = {globals_template};
+end
 % interpolate totals on t_grid
-ec_total_pow = interp1(ec_data.p_abs_plasma.t,ec_data.p_abs_plasma.data(10,:),t_grid); %the index 10 works for shot 56044
+ec_total_pow = interp1(ec_data.p_abs_plasma.t,ec_data.p_abs_plasma.data(nb_launchers+1,:),t_grid); %the index 10 works for shot 56044
 ec_total_pow(isnan(ec_total_pow)) = 0;
-ec_total_cur = interp1(ec_data.cd_tot.t,ec_data.cd_tot.data(10,:),t_grid);
+ec_total_cur = interp1(ec_data.cd_tot.t,ec_data.cd_tot.data(nb_launchers+1,:),t_grid);
 ec_total_cur(isnan(ec_total_cur)) = 0;
 
-% bs
-id_bs.description = 'Bootstrap current';
-id_bs.index = 13; id_bs.name = 'bootstrap_cuurent';
-ids_core_sources.source{13} = source_template;
-ids_core_sources.source{13}.identifier = id_bs;
-ids_core_sources.source{13}.profiles_1d(1:n_t) = {profiles_template};
-ids_core_sources.source{13}.global_quantities(1:n_t) = {globals_template};
-% load data
-bs_gdat = gdat(shot,'bs_data'); bs_data = bs_gdat.bs.bs_data;
+for i = 1:nb_launchers
+    for ii = 1:n_t
+    % ec, need to find matching time index\
+    ids_core_sources.source{last_index+i}.profiles_1d{ii}.time = t_grid(ii); % profiles time (should profiles time be ec_time?)
+    ids_core_sources.source{last_index+i}.global_quantities{ii}.time = t_grid(ii); % globals time
+    [~, ec_idx] = min(abs(ec_time-t_grid(ii)));
+    if ec_time(ec_idx) == t_grid(ii) % same grid 
+        % current density
+        ids_core_sources.source{last_index+i}.profiles_1d{ii}.j_parallel = ...
+          ec_data.cd_dens.data(:,end,ec_idx)';
+        % integrated current density
+        ids_core_sources.source{last_index+i}.profiles_1d{ii}.current_parallel_inside = ...
+          ec_data.cd_integrated.data(:,ii)';
+    
+    else
+        ids_core_sources.source{last_index+i}.profiles_1d{ii}.j_parallel = ...
+          zeros(1,numel(ec_data.cd_dens.x(1,:)))';
+        ids_core_sources.source{last_index+i}.profiles_1d{ii}.current_parallel_inside = ...
+          zeros(1,numel(ec_data.cd_integrated.x(1,:)))';
+    end
+    % globals
+    ids_core_sources.source{last_index+i}.global_quantities{ii}.time = t_grid(ii);
+    ids_core_sources.source{last_index+i}.global_quantities{ii}.power = ec_total_pow(ii);
+    ids_core_sources.source{last_index+i}.global_quantities{ii}.current_parallel = ec_total_cur(ii);
+    end
+end
+
 
 % nbi  
 id_nbi.description = 'Source from Neutral Beam Injection';
 id_nbi.index = 2; id_nbi.name = 'nbi';
-ids_core_sources.source{2} = source_template;
-ids_core_sources.source{2}.identifier = id_nbi;
-ids_core_sources.source{2}.profiles_1d(1:n_t) = {profiles_template};
-ids_core_sources.source{2}.global_quantities(1:n_t) = {globals_template};
+%ids_core_sources.source{} = source_template;
+%ids_core_sources.source{}.identifier = id_nbi;
+%ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
+%ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template};
 
 
 % total  
 id_total.description = 'Total source; combines all sources';
 id_total.index = 1; id_total.name = 'total';
-ids_core_sources.source{1} = source_template;
-ids_core_sources.source{1}.identifier = id_total;
-ids_core_sources.source{1}.profiles_1d(1:n_t) = {profiles_template};
-ids_core_sources.source{1}.global_quantities(1:n_t) = {globals_template};
+%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};
+
 
 
-powers_gdat = gdat(shot,'powers');
 
 %% fill profiles for times n t_grid
-for ii = 1:n_t
-  % ohm
-  % 1d_profiles
-  ids_core_sources.source{7}.profiles_1d{ii}.time = t_grid(ii);
-  ids_core_sources.source{7}.profiles_1d{ii}.j_parallel = ...
-    ohm_data.cd_dens.data(:,ii)';
-  % globals
-  ids_core_sources.source{7}.global_quantities{ii}.time = t_grid(ii);
-  ids_core_sources.source{7}.global_quantities{ii}.power = powers_gdat.ohm.data(ii);
-    
-  % ec, need to find matching time index\
-  ids_core_sources.source{3}.profiles_1d{ii}.time = t_grid(ii); % profiles time (should profiles time be ec_time?)
-  ids_core_sources.source{3}.global_quantities{ii}.time = t_grid(ii); % globals time
-  [~, ec_idx] = min(abs(ec_time-t_grid(ii)));
-  if ec_time(ec_idx) == t_grid(ii) % same grid 
-    % current density
-    ids_core_sources.source{3}.profiles_1d{ii}.j_parallel = ...
-      ec_data.cd_dens.data(:,end,ec_idx)';
-    % integrated current density
-    ids_core_sources.source{3}.profiles_1d{ii}.current_parallel_inside = ...
-      ec_data.cd_integrated.data(:,ii)';
-    
-  else
-    ids_core_sources.source{3}.profiles_1d{ii}.j_parallel = ...
-      zeros(1,numel(ec_data.cd_dens.x(1,:)))';
-    ids_core_sources.source{3}.profiles_1d{ii}.current_parallel_inside = ...
-      zeros(1,numel(ec_data.cd_integrated.x(1,:)))';
-  end
-  % globals
-  ids_core_sources.source{3}.global_quantities{ii}.time = t_grid(ii);
-  ids_core_sources.source{3}.global_quantities{ii}.power = ec_total_pow(ii);
-  ids_core_sources.source{3}.global_quantities{ii}.current_parallel = ec_total_cur(ii);
-  
-  % bs
-  ids_core_sources.source{13}.profiles_1d{ii}.j_parallel = ...
-    bs_data.cd_dens.data(:,ii)';
-    
-end
+% for ii = 1:n_t
+%   % ohm
+%   % 1d_profiles
+%   ids_core_sources.source{1}.profiles_1d{ii}.time = t_grid(ii);
+%   ids_core_sources.source{1}.profiles_1d{ii}.j_parallel = ...
+%     ohm_data.cd_dens.data(:,ii)';
+%   % globals
+%   ids_core_sources.source{1}.global_quantities{ii}.time = t_grid(ii);
+%   ids_core_sources.source{1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii);
+%     
+%   % ec, need to find matching time index\
+%   ids_core_sources.source{4}.profiles_1d{ii}.time = t_grid(ii); % profiles time (should profiles time be ec_time?)
+%   ids_core_sources.source{4}.global_quantities{ii}.time = t_grid(ii); % globals time
+%   [~, ec_idx] = min(abs(ec_time-t_grid(ii)));
+%   if ec_time(ec_idx) == t_grid(ii) % same grid 
+%     % current density
+%     ids_core_sources.source{4}.profiles_1d{ii}.j_parallel = ...
+%       ec_data.cd_dens.data(:,end,ec_idx)';
+%     % integrated current density
+%     ids_core_sources.source{4}.profiles_1d{ii}.current_parallel_inside = ...
+%       ec_data.cd_integrated.data(:,ii)';
+%     
+%   else
+%     ids_core_sources.source{4}.profiles_1d{ii}.j_parallel = ...
+%       zeros(1,numel(ec_data.cd_dens.x(1,:)))';
+%     ids_core_sources.source{4}.profiles_1d{ii}.current_parallel_inside = ...
+%       zeros(1,numel(ec_data.cd_integrated.x(1,:)))';
+%   end
+%   % globals
+%   ids_core_sources.source{4}.global_quantities{ii}.time = t_grid(ii);
+%   ids_core_sources.source{4}.global_quantities{ii}.power = ec_total_pow(ii);
+%   ids_core_sources.source{4}.global_quantities{ii}.current_parallel = ec_total_cur(ii);
+%   
+%   % bs
+%   ids_core_sources.source{2}.profiles_1d{ii}.j_parallel = ...
+%     bs_data.cd_dens.data(:,ii)';
+%     
+% end
 
 %% add descriptions for profiles_1d