From a3b6f7ab51e955199dc7501e8da3f8448bbaa389 Mon Sep 17 00:00:00 2001
From: Aliki Chatzilouka <chatzilouka.aliki@epfl.ch>
Date: Thu, 14 Mar 2024 18:17:26 +0100
Subject: [PATCH] fixing time grid for ec sources
---
matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 68 ++++++++++++----------
1 file changed, 38 insertions(+), 30 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index 0a4a3963..652a2339 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/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';
--
GitLab