From 3f9f872ed0e40b779eaa2a2d690fb2750259ce62 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Wed, 19 Aug 2020 19:11:34 +0200
Subject: [PATCH] move private flux points psi_norm to max to be outside

---
 matlab/TCV/gdat_tcv.m | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 619fb239..55aaa180 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1997,6 +1997,20 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         warning('Found negative values of psi_norm at TS positions, replacing them with 0');
         psi_norm(mask) = 0;
       end
+      params_eff = gdat_data.gdat_params;
+      params_eff.doplot=0;
+      params_eff.data_request='z_axis';
+      z_axis = gdat(gdat_data.shot,params_eff);
+      ii=iround_os(psiscatvol.dim{1},z_axis.data);
+      is_psi_plasma = psi_norm < 1;
+      is_in_lcfs = zeros(size(psi_norm));
+      for it=1:size(is_psi_plasma,2)
+        is_in_lcfs(ii(it):-1:1,it)=cumprod(is_psi_plasma(ii(it):-1:1,it));
+        is_in_lcfs(ii(it):end,it)=cumprod(is_psi_plasma(ii(it):end,it));
+      end
+      is_private_flux = is_psi_plasma & ~is_in_lcfs;
+      % at this stage set private flux "rho" to max, could use 1.3^2
+      psi_norm(is_private_flux) = max(max(psi_norm));
       rho = sqrt(psi_norm);
     else
       if gdat_params.nverbose>=1 && gdat_data.gdat_params.edge==0
-- 
GitLab