From 372705b5ac04dda6377c054681cee8f39de56c85 Mon Sep 17 00:00:00 2001
From: Aliki Chatzilouka <chatzilouka.aliki@epfl.ch>
Date: Sun, 17 Mar 2024 21:52:54 +0100
Subject: [PATCH] adding electron energy profile to structure
---
matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 110 +++++++++++----------
1 file changed, 57 insertions(+), 53 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index 652a2339..d30dd1a6 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/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';
--
GitLab