From 4e0aceb19987286bc5ea6715966d561be4ca2ff1 Mon Sep 17 00:00:00 2001
From: Antonia Frank <antonia.frank@epfl.ch>
Date: Tue, 13 Aug 2024 13:47:17 +0200
Subject: [PATCH] Add DNBI similarily to NBI

---
 matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 61 ++++++++++++----------
 1 file changed, 34 insertions(+), 27 deletions(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index afeb88d1..d34d1a14 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -321,8 +321,7 @@ if numel(active_nbi)>0
   nbi_powers_tgrid = powers_gdat.(nbi_names{active_nbi(1)}).t;
   % find times where NBI is on to define time grid with extra time slice just
   % before/after  NBI power and at start/end of shot
-  itime_nbi = find((powers_gdat.(nbi_names{active_nbi(1)}).data>0));% & ...
-  %  powers_gdat.(nbi_names{active_nbi(1)}).t<(ohm_t_grid(end)));
+  itime_nbi = find((powers_gdat.(nbi_names{active_nbi(1)}).data>0));
   if nbi_powers_tgrid(itime_nbi(end))>=ohm_t_grid(end)
     i_time_end = iround(nbi_powers_tgrid,ohm_t_grid(end));
     nbi_tgrid_out = [ohm_t_grid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
@@ -342,8 +341,8 @@ if numel(active_nbi)>0
     ids_core_sources.source{last_index+i_nbi}.global_quantities(1:nt_nbi_out) = {globals_template};
   end
 
-  sprintf('Loading of current & power densities from ASTRA not implemented yet. \n')
-  sprintf('Checking if ASTRA run available on partition /Lac8_D:')
+  disp('Loading of current & power densities from ASTRA not implemented yet.')
+  disp('Checking if ASTRA run available on partition /Lac8_D:')
   [~,hostname] = unix('hostname');
   if strcmp(hostname,'lac8.epfl.ch')
     unix(sprintf('ls /Lac8_D/ASTRA/ | grep ''%i'' && echo File for shotnumber exists! || echo File for shotnumber does not exist!',shot));
@@ -369,29 +368,36 @@ if numel(active_nbi)>0
 end
 
 %% DNBI has it's own time grid
-% if ~isempty(powers_gdat.nbi1)
-%
-% id_nbi.description = 'Source from Neutral Beam Injection';
-% id_nbi.index = 2; id_nbi.name = 'nbi';
-% ids_core_sources.source{last_index+1} = source_template;
-% ids_core_sources.source{last_index+1}.identifier = id_nbi;
-% 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};
-%
-%   % find times of injected EC power to interpolate power & current densities
-%   % on p_ec_injected tgrid
-%   ec_powers_tgrid = powers_gdat.ec.t; %nt_ec_powers = numel(ec_powers_tgrid);
-%   % find times where EC is on to define time grid with extra time slice just
-%   % before/after  EC power and at start/end of shot
-%   itime_ec = find(powers_gdat.ec.data(:,end)>0);
-%   ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_t_grid(end)];
-%   nt_ec_out = numel(ec_tgrid_out);
-%   p_ec_injected = interpos(ec_powers_tgrid,powers_gdat.ec.data(:,end),ec_tgrid_out);
-%
-%   last_index = last_index+nb_launchers;
-% end
-
+if ~isempty(powers_gdat.dnbi)
+  % get tgrid (same for NBI1 & 2 if both active) from active_nbi(1)
+  dnbi_powers_tgrid = powers_gdat.dnbi.t;
+  % find times where DNBI is on to define time grid with extra time slice just
+  % before/after  DNBI power and at start/end of shot
+  itime_dnbi = find((powers_gdat.dnbi.data>0));
+  if dnbi_powers_tgrid(itime_dnbi(end))>=ohm_t_grid(end)
+    i_time_end = iround(dnbi_powers_tgrid,ohm_t_grid(end));
+    dnbi_tgrid_out = [ohm_t_grid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:i_time_end)'];
+  else
+    dnbi_tgrid_out = [ohm_t_grid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:itime_dnbi(end)+1)',ohm_t_grid(end)];
+  end
+  nt_dnbi_out = numel(dnbi_tgrid_out);
+
+  id_dnbi.description = 'DNBI Source from Neutral Beam Injection';
+  id_dnbi.index = 2; id_dnbi.name = 'nbi';
+  ids_core_sources.source{last_index+1} = source_template;
+  ids_core_sources.source{last_index+1}.identifier = id_dnbi;
+  ids_core_sources.source{last_index+1}.profiles_1d(1:nt_dnbi_out) = {profiles_template};
+  ids_core_sources.source{last_index+1}.global_quantities(1:nt_dnbi_out) = {globals_template};
+
+  p_dnbi_injected = interpos(dnbi_powers_tgrid,powers_gdat.dnbi.data,dnbi_tgrid_out);
+  for ii = 1:nt_dnbi_out
+    % globals
+    ids_core_sources.source{last_index+i_nbi}.global_quantities{ii}.time = nbi_tgrid_out(ii);
+    ids_core_sources.source{last_index+i_nbi}.global_quantities{ii}.power = p_dnbi_injected(ii);
+  end
 
+  last_index = last_index+1;
+end
 
 %% total
 % id_total.description = 'Total source; combines all sources';
@@ -403,7 +409,8 @@ end
 
 %% add descriptions for profiles_1d
 
-desc.source = 'available by now {source_index}: ohmic {1}, bootstrap {2}, ec {3+} for individual launchers (if present in shot)';
+desc.source = ...
+  'available by now {source_index}: ohmic {1}, bootstrap {2}, ec {3+} for individual launchers (if present in shot), nbi {3+|14+} if available, dnbi {3|14|16} if available';
 desc.profiles_1d.electrons.energy = 'absorbed power density';
 desc.profiles_1d.electrons.power_inside = 'integrated absorbed power density';
 desc.profiles_1d.j_parallel = 'parallel current density';
-- 
GitLab