diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m index 3c8d3fd6520a1eccac019cabdadba149bd51de08..9e781fde81a724b5522e47d908102d6382509e1b 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]);