Fix bug in loading of NBI data, NBI zero for all times

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

@@ -380,61 +380,64 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
 end
 
 %% NBI. same tgrid for NBI1 and NBI2
-
 check_nbi = [~isempty(powers_gdat.nbi1.data),~isempty(powers_gdat.nbi2.data)];
 active_nbi = find(check_nbi==1);
 nb_nbi = numel(check_nbi);
 nbi_names = {'nbi1','nbi2'};
 
-if numel(active_nbi)>0
-  % get tgrid (same for NBI1 & 2 if both active) from active_nbi(1)
-  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));
-  if nbi_powers_tgrid(itime_nbi(end))>=ohm_tgrid(end)
-    i_time_end = iround(nbi_powers_tgrid,ohm_tgrid(end));
-    nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
-  else
-    nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_tgrid(end)];
-  end
-  nt_nbi_out = numel(nbi_tgrid_out);
-
-  % Setup source structs for active nbi from template
-  main_desc = 'Source from Neutral Beam Injection';
-  id_nbi.index = 2; id_nbi.name = 'nbi';
-  for i_nbi = active_nbi
-    id_nbi.description = sprintf('NBI%i %s',i_nbi,main_desc);
-    ids_core_sources.source{last_index+i_nbi} = source_template;
-    ids_core_sources.source{last_index+i_nbi}.identifier = id_nbi;
-    ids_core_sources.source{last_index+i_nbi}.profiles_1d(1:nt_nbi_out) = {profiles_template};
-    ids_core_sources.source{last_index+i_nbi}.global_quantities(1:nt_nbi_out) = {globals_template};
-  end
+% Setup source structs for active nbi from template
+main_desc = 'Source from Neutral Beam Injection';
+id_nbi.index = 2; id_nbi.name = 'nbi';
 
-  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));
-  else
-    unix(sprintf('ssh $(whoami)@lac8 "ls /Lac8_D/ASTRA/ | grep ''%i'' && echo File for shotnumber exists! || echo File for shotnumber does not exist! && exit"',shot));
-  end
+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));
+else
+  unix(sprintf('ssh $(whoami)@lac8 "ls /Lac8_D/ASTRA/ | grep ''%i'' && echo File for shotnumber exists! || echo File for shotnumber does not exist! && exit"',shot));
+end
 
+if numel(active_nbi)>0
   for i_nbi = active_nbi
-    p_nbi_injected_tmp = interpos(nbi_powers_tgrid,powers_gdat.(nbi_names{i_nbi}).data,nbi_tgrid_out);
-    for ii = 1:nt_nbi_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_nbi_injected_tmp(ii);
+    % get tgrid (same for NBI1 & 2 if both active) from active_nbi(1)
+    nbi_powers_tgrid = powers_gdat.(nbi_names{i_nbi}).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{i_nbi}).data>0));
+    if ~isempty(itime_nbi)
+      if nbi_powers_tgrid(itime_nbi(end))>=ohm_tgrid(end)
+        i_time_end = iround(nbi_powers_tgrid,ohm_tgrid(end));
+        nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
+      else
+        nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_tgrid(end)];
+      end
+      nt_nbi_out = numel(nbi_tgrid_out);
+
+      % Setup source structs for active nbi from template
+      id_nbi.description = sprintf('NBI%i %s',i_nbi,main_desc);
+      ids_core_sources.source{last_index+i_nbi} = source_template;
+      ids_core_sources.source{last_index+i_nbi}.identifier = id_nbi;
+      ids_core_sources.source{last_index+i_nbi}.profiles_1d(1:nt_nbi_out) = {profiles_template};
+      ids_core_sources.source{last_index+i_nbi}.global_quantities(1:nt_nbi_out) = {globals_template};
+
+      p_nbi_injected_tmp = interpos(nbi_powers_tgrid,powers_gdat.(nbi_names{i_nbi}).data,nbi_tgrid_out);
+      for ii = 1:nt_nbi_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_nbi_injected_tmp(ii);
+      end
+    else
+      disp('NBI power is zero for all times, skip source.')
     end
   end
-
+  
   % add empty for unused NBI
-  if numel(ids_core_sources.source)-last_index ~= nb_nbi
+  if numel(ids_core_sources.source)-last_index == 1 % if 1 then add second empty source for NBI2
     ids_core_sources.source{last_index+nb_nbi} = [];
   end
 
-  last_index = last_index+nb_nbi;
+  last_index = last_index+(numel(ids_core_sources.source)-last_index);
 end
 
 %% DNBI has it's own time grid
@@ -444,29 +447,32 @@ if ~isempty(powers_gdat.dnbi)
   % 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_tgrid(end)
-    i_time_end = iround(dnbi_powers_tgrid,ohm_tgrid(end));
-    dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:i_time_end)'];
+  if ~isempty(itime_dnbi)
+    if dnbi_powers_tgrid(itime_dnbi(end))>=ohm_tgrid(end)
+      i_time_end = iround(dnbi_powers_tgrid,ohm_tgrid(end));
+      dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:i_time_end)'];
+    else
+      dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:itime_dnbi(end)+1)',ohm_tgrid(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+1}.global_quantities{ii}.time = dnbi_tgrid_out(ii);
+      ids_core_sources.source{last_index+1}.global_quantities{ii}.power = p_dnbi_injected(ii);
+    end   
+    last_index = last_index+1;  
   else
-    dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:itime_dnbi(end)+1)',ohm_tgrid(end)];
+    disp('DNBI power is zero for all times, skip source.')
   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+1}.global_quantities{ii}.time = nbi_tgrid_out(ii);
-    ids_core_sources.source{last_index+1}.global_quantities{ii}.power = p_dnbi_injected(ii);
-  end
-
-  last_index = last_index+1;
 end
 
 %% total