fixing time grid for ec sources

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

@@ -36,51 +36,58 @@ globals_template = source_template.global_quantities{1};
 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
+ohm_t_grid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_t_grid); t_spacing = mean(diff(ohm_t_grid));
+
 id_ohm.description = 'Source from ohmic heating';
 id_ohm.index = 7; id_ohm.name = 'ohmic';
-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};
+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;
 
 %fill profiles for times n t_grid
-for ii = 1:n_t
+
+for ii = 1:ohm_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 = ...
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = ohm_t_grid(ii);
+  ids_core_sources.source{last_index+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);
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.time = ohm_t_grid(ii);
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii);
 end
-last_index = last_index+1;
+
+last_index = last_index+1;  % add if statement to only increment if ohmic source has been added
+
 %% bs
 id_bs.description = 'Bootstrap current';
-id_bs.index = 13; id_bs.name = 'bootstrap_cuurent';
+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:n_t) = {profiles_template};
-ids_core_sources.source{last_index+1}.global_quantities(1:n_t) = {globals_template};
+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};
 % 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
+for ii = 1:ohm_n_t
    ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = ...
-    bs_data.cd_dens.data(:,ii)'; 
+     bs_data.cd_dens.data(:,ii)'; 
 end
-last_index = last_index+1;
+ 
+last_index = last_index+1;  % add if statement to only increment if bs source has been added
+
 %% ec
 % load data
 ec_gdat = gdat(shot,'ec_data'); ec_data = ec_gdat.ec.ec_data;
 ec_time = ec_data.cd_dens.t;
+ec_t_grid = powers_gdat.ec.t; ec_n_t = numel(ec_t_grid); t_spacing = mean(diff(ec_t_grid));
 
 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';
@@ -88,22 +95,23 @@ 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};
+    ids_core_sources.source{last_index+i}.profiles_1d(1:ec_n_t) = {profiles_template};
+    ids_core_sources.source{last_index+i}.global_quantities(1:ec_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(nb_launchers+1,:),t_grid); %the index 10 works for shot 56044
+
+% interpolate totals on t_grid ??
+ec_total_pow = ec_data.p_abs_plasma.data(nb_launchers+1,:); 
 ec_total_pow(isnan(ec_total_pow)) = 0;
-ec_total_cur = interp1(ec_data.cd_tot.t,ec_data.cd_tot.data(nb_launchers+1,:),t_grid);
+ec_total_cur = ec_data.cd_tot.data(nb_launchers+1,:);
 ec_total_cur(isnan(ec_total_cur)) = 0;
 
 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 
+    for ii = 1:ec_n_t
+    % ec, need to find matching time index
+    ids_core_sources.source{last_index+i}.profiles_1d{ii}.time = ec_t_grid(ii); % profiles time (should profiles time be ec_time?)
+    ids_core_sources.source{last_index+i}.global_quantities{ii}.time = ec_t_grid(ii); % globals time
+    [~, ec_idx] = min(abs(ec_time-ec_t_grid(ii)));
+    if ec_time(ec_idx) == ec_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)';
@@ -118,12 +126,12 @@ for i = 1:nb_launchers
           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}.time = ec_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
-
+last_index = last_index+nb_launchers;
 
 % nbi  
 id_nbi.description = 'Source from Neutral Beam Injection';