Add implementation and check for ASTAR file in Lac8D for NBI

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

@@ -124,7 +124,7 @@ for ii = 1:bs_n_t
 
   % globals
   ids_core_sources.source{last_index+1}.global_quantities{ii}.time = bs_t_grid(ii);
-  ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);  
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
 end
 
 last_index = last_index+1;  % add if statement to only increment if bs source has been added
@@ -147,10 +147,10 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
 
   % 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 
+  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);  
+  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);
@@ -246,7 +246,7 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
       interp_norm_cd_integrated(irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_cd_integrated,ec_tgrid_out(3:end-2));
     end
   end
-  
+
   % interpolate the rho_pol & rho_tor on ec_powers_tgrid
   interp_rho_pol = zeros(n_rho,nt_ec_out);
   interp_rho_tor = zeros(n_rho,nt_ec_out);
@@ -254,7 +254,7 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
     interp_rho_pol(irho,3:end-2) = interpos(ec_tgrid_toray,rho_pol_norm(irho,:),ec_tgrid_out(3:end-2));
     interp_rho_tor(irho,3:end-2) = interpos(ec_tgrid_toray,rho_tor_norm(irho,:),ec_tgrid_out(3:end-2));
   end
-  % fill rho_pol just outside of TORAY times(zero power) with the known profiles one time slice after/before 
+  % fill rho_pol just outside of TORAY times(zero power) with the known profiles one time slice after/before
   interp_rho_pol(:,2) = interp_rho_pol(:,3); interp_rho_pol(:,end-1) = interp_rho_pol(:,end-2);
 
   % normalised & interpolated profiles * p_ec_injected on interp_tgrid
@@ -278,30 +278,30 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
       if ii~=1 || ii~=nt_ec_out
         ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = interp_rho_tor(:,ii);
         ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_pol_norm = interp_rho_pol(:,ii);
-      % power density
-      ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.energy = ...
-        interp_p_dens(:,i_lau,ii);
-      % integrated power density
-      ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.power_inside = ...
-        interp_p_integrated(:,i_lau,ii);
-      % current density
-      ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ...
-        interp_cd_dens(:,i_lau,ii);
-      % integrated current density
-      ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
-        interp_cd_integrated(:,i_lau,ii);
+        % power density
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.energy = ...
+          interp_p_dens(:,i_lau,ii);
+        % integrated power density
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.power_inside = ...
+          interp_p_integrated(:,i_lau,ii);
+        % current density
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ...
+          interp_cd_dens(:,i_lau,ii);
+        % integrated current density
+        ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
+          interp_cd_integrated(:,i_lau,ii);
       end
-     
+
       % globals
       ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.time = ec_tgrid_out(ii);
       ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.power = ...
         interp_p_integrated(end,i_lau,ii);
       ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.current_parallel = ...
-        interp_cd_integrated(end,i_lau,ii);  
+        interp_cd_integrated(end,i_lau,ii);
     end
   end
 
-  %add empty sources for rest of unsused launchers
+  %add empty sources for rest of unused launchers
   if numel(ids_core_sources.source)-last_index ~= nb_launchers
     ids_core_sources.source{last_index+nb_launchers} = [];
   end
@@ -309,13 +309,89 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
   last_index = last_index+nb_launchers;
 end
 
-%% nbi
+%% 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));% & ...
+  %  powers_gdat.(nbi_names{active_nbi(1)}).t<(ohm_t_grid(end)));
+  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)'];
+  else
+    nbi_tgrid_out = [ohm_t_grid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_t_grid(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
+
+  sprintf('Loading of current & power densities from ASTRA not implemented yet. \n')
+  sprintf('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
+
+  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);
+    end
+  end
+
+  % add empty for unused NBI
+  if numel(ids_core_sources.source)-last_index ~= nb_nbi
+    ids_core_sources.source{last_index+nb_nbi} = [];
+  end
+
+  last_index = last_index+nb_nbi;
+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{} = source_template;
-% ids_core_sources.source{}.identifier = id_nbi;
-% ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template};
-% ids_core_sources.source{}.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_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
+
+
 
 %% total
 % id_total.description = 'Total source; combines all sources';