From aec25717d5379444a231323ec8fc0b283d89e2a2 Mon Sep 17 00:00:00 2001
From: Antonia Frank <antonia.frank@epfl.ch>
Date: Tue, 13 Aug 2024 16:05:38 +0200
Subject: [PATCH] Resolve problem of Nans in EC sources by adding 21 to iterpos
 call for interpolation

---
 matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 29 ++++++++++++++--------
 1 file changed, 19 insertions(+), 10 deletions(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index d34d1a14..7380e42a 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -151,7 +151,12 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
   % 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)];
+  if ec_powers_tgrid(itime_ec(end))>=ohm_t_grid(end)
+    i_time_end = iround(ec_powers_tgrid,ohm_t_grid(end));
+    ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:i_time_end)'];
+  else
+    ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_t_grid(end)];
+  end
   nt_ec_out = numel(ec_tgrid_out);
   p_ec_injected = interpos(ec_powers_tgrid,powers_gdat.ec.data(:,end),ec_tgrid_out);
 
@@ -194,10 +199,13 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
 
   % load 1d_profiles from ec_data
   p_dens        = ec_data.p_dens.data;
+  p_dens(isnan(p_dens)) = 0;
   p_integrated  = ec_data.p_integrated.data;
+  p_integrated(isnan(p_integrated)) = 0;
   cd_dens       = ec_data.cd_dens_doublewidth.data; % use double width to have realistic deposition broadening
-  cd_integrated = ec_data.cd_integrated.data;
-
+  cd_dens(isnan(cd_dens)) = 0;
+  
+  cd_integrated = zeros(size(ec_data.cd_integrated.data));
   for i_lau = active_launchers
     % convert to j// with conversion factor
     cd_dens(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
@@ -239,11 +247,12 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
       trace_p_integrated  = squeeze(norm_p_integrated(irho,i_lau,:));
       trace_cd_dens       = squeeze(norm_cd_dens(irho,i_lau,:));
       trace_cd_integrated = squeeze(norm_cd_integrated(irho,i_lau,:));
-      % interpolate on gdat powers tgrid
-      interp_norm_p_dens(       irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_p_dens,       ec_tgrid_out(3:end-2));
-      interp_norm_p_integrated( irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_p_integrated, ec_tgrid_out(3:end-2));
-      interp_norm_cd_dens(      irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_cd_dens,      ec_tgrid_out(3:end-2));
-      interp_norm_cd_integrated(irho,i_lau,3:end-2) = interpos(ec_tgrid_toray,trace_cd_integrated,ec_tgrid_out(3:end-2));
+      % interpolate on gdat powers tgrid, interpos 21 to extrapolate,
+      % taking constant value y_edge; needed for time window between start of EC and first TS time 
+      interp_norm_p_dens(       irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_dens,       ec_tgrid_out(3:end-2));
+      interp_norm_p_integrated( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_integrated, ec_tgrid_out(3:end-2));
+      interp_norm_cd_dens(      irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_cd_dens,      ec_tgrid_out(3:end-2));
+      interp_norm_cd_integrated(irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_cd_integrated,ec_tgrid_out(3:end-2));
     end
   end
 
@@ -251,8 +260,8 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
   interp_rho_pol = zeros(n_rho,nt_ec_out);
   interp_rho_tor = zeros(n_rho,nt_ec_out);
   for irho = 1:n_rho
-    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));
+    interp_rho_pol(irho,3:end-2) = interpos(21,ec_tgrid_toray,rho_pol_norm(irho,:),ec_tgrid_out(3:end-2));
+    interp_rho_tor(irho,3:end-2) = interpos(21,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
   interp_rho_pol(:,2) = interp_rho_pol(:,3); interp_rho_pol(:,end-1) = interp_rho_pol(:,end-2);
-- 
GitLab