From b6c0d045d3f12654da1c31ca5c91c893d4f08344 Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Fri, 11 Oct 2019 14:53:17 +0200
Subject: [PATCH] Avoid computing complex sqrt ...
---
matlab/TCV/gdat_tcv.m | 14 +++++++++-----
1 file changed, 9 insertions(+), 5 deletions(-)
diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 3c8d3fd6..9e781fde 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1668,11 +1668,15 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.psi_max = psi_max;
end
if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
- rho = sqrt(1.-psiscatvol.data./repmat(psi_max.data(:).',[size(psiscatvol.data,1),1]));
- if any(imag(rho(:)))
- warning('Found complex values when computing rho from psi TS, replacing them with 0');
- rho(imag(rho(:))) = 0;
- end
+ psi_norm = 1.-psiscatvol.data./repmat(psi_max.data(:).',[size(psiscatvol.data,1),1]);
+ mask = psi_norm < 0;
+ if any(mask(:))
+ % NOTE: I am not entirely sure there could not be a case where this
+ % is actually possible outside of the confined plasma
+ warning('Found negative values of psi_norm at TS positions, replacing them with 0');
+ psi_norm(mask) = 0;
+ end
+ rho = sqrt(psi_norm);
else
if gdat_params.nverbose>=1 && gdat_data.gdat_params.edge==0
warning(['psiscatvol empty?, no rho calculated for data_request = ' data_request_eff]);
--
GitLab