adding electron energy profile to structure

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

@@ -77,6 +77,7 @@ bs_gdat = gdat(shot,'bs_data'); bs_data = bs_gdat.bs.bs_data;
 
 % fill profiles for times n t_grid
 for ii = 1:ohm_n_t
+   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 = ...
      bs_data.cd_dens.data(:,ii)'; 
 end
@@ -86,7 +87,7 @@ last_index = last_index+1;  % add if statement to only increment if bs source ha
 %% ec
 % load data
 ec_gdat = gdat(shot,'ec_data'); ec_data = ec_gdat.ec.ec_data;
-ec_time = ec_data.cd_dens.t;
+ec_time = ec_data.p_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);
@@ -98,9 +99,8 @@ for i = 1:nb_launchers
     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 = ec_data.p_abs_plasma.data(nb_launchers+1,:); 
+ 
+ec_total_pow = transpose(powers_gdat.ec.data(:,nb_launchers+1)); %use power from powers_gdat(injected ec power) instead of ec_data power
 ec_total_pow(isnan(ec_total_pow)) = 0;
 ec_total_cur = ec_data.cd_tot.data(nb_launchers+1,:);
 ec_total_cur(isnan(ec_total_cur)) = 0;
@@ -108,22 +108,22 @@ ec_total_cur(isnan(ec_total_cur)) = 0;
 for i = 1:nb_launchers
     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}.profiles_1d{ii}.time = ec_t_grid(ii); % profiles 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
+        % current density (to adapt to <J.B>/B0)
         ids_core_sources.source{last_index+i}.profiles_1d{ii}.j_parallel = ...
           ec_data.cd_dens.data(:,end,ec_idx)';
-        % integrated current density
+        % integrated current density (to adapt to INTEGRAL(<J.B>/B0)*ds_phi)
         ids_core_sources.source{last_index+i}.profiles_1d{ii}.current_parallel_inside = ...
-          ec_data.cd_integrated.data(:,ii)';
+          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,:)))';
+        %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 = ec_t_grid(ii);
@@ -131,9 +131,55 @@ for i = 1:nb_launchers
     ids_core_sources.source{last_index+i}.global_quantities{ii}.current_parallel = ec_total_cur(ii);
     end
 end
+
+% interpoating p_dens profiles from 'ec_time' grid (toray time nodes) to 'ec_t_grid' (injected power time)
+
+p_dens = ec_data.p_dens.data;
+p_ec_injected = powers_gdat.ec.data;
+it = iround_os(ec_t_grid,ec_time);
+sparse_p_ec_injected = p_ec_injected(it,:); % injected ec power vals corresponding to ec_time grid
+
+% calculate normalised profiles on ec_time grid
+norm_p_dens = zeros(size(p_dens, 1), nb_launchers+1, length(ec_time));
+for t = 1:length(ec_time)
+    normalised = zeros(size(p_dens, 1), nb_launchers+1);
+    for launcher = 1:nb_launchers+1
+        normalised(:, launcher) = p_dens(:, launcher, t) ./ sparse_p_ec_injected(t, launcher);
+    end
+    norm_p_dens(:,:,t) = normalised;
+end
+
+% interpolate normalised p_dens profiles on ec_t_grid 
+interp_norm_p_dens = zeros(size(norm_p_dens, 1), size(norm_p_dens, 2), ec_n_t);
+for rho = 1:size(norm_p_dens, 1)
+    for launcher = 1:size(norm_p_dens, 2)
+        profile = squeeze(norm_p_dens(rho, launcher, :));
+        interp_profile = interp1(ec_time, profile, ec_t_grid);
+        interp_norm_p_dens(rho, launcher, :) = interp_profile;
+    end
+end
+
+% normalised & interpolated p_dens profiles * p_ec_injected on ec_t_grid
+interp_p_dens = zeros(size(interp_norm_p_dens, 1), size(interp_norm_p_dens, 2), ec_n_t);
+for t = 1:length(ec_t_grid)
+    unnormalised = zeros(size(p_dens, 1), nb_launchers+1);
+    for launcher = 1:nb_launchers+1
+        unnormalised(:,launcher) = interp_norm_p_dens(:,launcher,t).* p_ec_injected(t, launcher);
+    end
+    interp_p_dens(:,:,t) = unnormalised;
+end
+
+for ii = 1:ec_n_t
+    for i = 1:nb_launchers
+        ids_core_sources.source{last_index+i}.profiles_1d{ii}.grid.rho_tor_norm = ...
+          ec_data.p_dens.rhotor_norm(1,:);
+        ids_core_sources.source{last_index+i}.profiles_1d{ii}.electrons.energy = ...
+          interp_p_dens(:,i,ii);
+    end
+end
 last_index = last_index+nb_launchers;
 
-% nbi  
+%% nbi  
 id_nbi.description = 'Source from Neutral Beam Injection';
 id_nbi.index = 2; id_nbi.name = 'nbi';
 %ids_core_sources.source{} = source_template;
@@ -151,48 +197,6 @@ id_total.index = 1; id_total.name = 'total';
 %ids_core_sources.source{}.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);
-%     
-%   % 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
 
 desc.source = 'available by now, 1:total, 3:ec, 7:ohm, 13:bootstrap';