diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 56d436ce8f5e95f0e72dac0763d3cb3aca775138..9e781fde81a724b5522e47d908102d6382509e1b 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1652,64 +1652,31 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
end
if liuqe_matlab==1 && strcmp(substr_liuqe,'_1'); substr_liuqe = ''; end
% psiscatvol obtained from linear interpolation in fact so not quite ok near axis where psi is quadratic
+ % NOTE: Since 2019 psi_scatvol uses psitbx with cubic-spline interpolation, unsure what shots have new data or not ...
recompute_psiscatvol_always = 1;
if liuqe_version==-1; recompute_psiscatvol_always = 1; end
+ % NOTE: Since oct.2019 psi_scatvol stored in nodes now uses LIUQE.M (shot range to check if any ...)
if all(abs(zshift)<1e-5) && liuqe_matlab==0 && recompute_psiscatvol_always==0
psi_max=gdat_tcv([],['\results::thomson' edge_str_dot ':psi_max' substr_liuqe],'nverbose',gdat_params.nverbose);
psiscatvol=gdat_tcv([],['\results::thomson' edge_str_dot ':psiscatvol' substr_liuqe],'nverbose',gdat_params.nverbose);
else
% calculate new psiscatvol
- if liuqe_matlab==0
- psitdi = tdi(['tcv_eq(''psi'',''LIUQE' substr_liuqe_tcv_eq ''')']);
- psiaxis = tdi(['tcv_eq(''psi_axis'',''LIUQE' substr_liuqe_tcv_eq ''')']);
- else
- psitdi = tdi(['tcv_eq(''psi'',''LIUQE.M' substr_liuqe_tcv_eq ''')']);
- psiaxis = tdi(['tcv_eq(''psi_axis'',''LIUQE.M' substr_liuqe_tcv_eq ''')']);
- end
- if numel(psitdi.dim)<1
- warning('problem with psitdi in gdat_tcv ')
- psitdi
- psiaxis
- return
- end
- rmesh=psitdi.dim{1};
- zmesh=psitdi.dim{2};
- zthom=mdsdata(['dim_of(\thomson' edge_str_dot ':te,1)']);
- zeffshift=zshift;
- % set zeffshift time array same as psitdi
- switch length(zeffshift)
- case 1
- zeffshift=zeffshift * ones(size(psitdi.dim{3}));
- case length(psitdi.dim{3})
- % ok
- case length(gdat_data.t)
- zeffshift=interp1(gdat_data.t,zeffshift,psitdi.dim{3});
- otherwise
- if (gdat_params.nverbose>=1);
- disp(' bad time dimension for zshift')
- disp(['it should be 1 or ' num2str(length(gdat_data.t)) ' or ' num2str(length(psitdi.dim{3}))])
- end
- end
- for it=1:length(gdat_data.t)
- itpsitdi=iround_os(psitdi.dim{3},gdat_data.t(it));
- psirz=psitdi.data(:,:,itpsitdi);
- %psiscatvol0=interp2(rmesh,zmesh,psirz',0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi),'spline'); % faster with interpos
- psiscatvol0=interpos2Dcartesian(rmesh,zmesh,psirz,0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi),-0.1,-0.1);
- %psiscatvol0=interp2(rmesh,zmesh,psirz',0.9*ones(size(zthom)),zthom-zeffshift(itpsitdi),'linear');
- psiscatvol.data(it,:)=psiscatvol0;
- % since take closest psi(R,Z) from psitdi, should also take closest for psi_max and not interpolating
- itpsiaxis = iround_os(psiaxis.dim{1},gdat_data.t(it));
- psi_max.data(it,1) = psiaxis.data(itpsiaxis);
- end
- psiscatvol.dim{1} = gdat_data.t;
- psiscatvol.dim{2} = gdat_data.x;
+ [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, edge_str_dot, psitbx_str, gdat_params);
+
+ % Add the results to the output of gdat
+ gdat_data.psiscatvol = psiscatvol;
+ gdat_data.psi_max = psi_max;
end
if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
- for ir=1:length(psiscatvol.dim{2})
- rho2 = max(1.-psiscatvol.data(:,ir)./psi_max.data,0);
- rho(ir,:)= sqrt(rho2');
- % rho(ir,:)= sqrt(1.-psiscatvol.data(:,ir)./psi_max.data(:))';
- 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]);
@@ -3123,3 +3090,83 @@ gdat_data.t=time;
if any(strcmp(fieldnames(tracetdi),'units'))
gdat_data.units=tracetdi.units;
end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, system_str_dot, psitbx_str, gdat_params )
+% Input arguments:
+% gdat_data: current structure containing TS data
+% zshift: vertical shift of the equilibrium, can be a vector for time-dependent data
+% system_str_dot: '' or '.edge' indicating thomson syste,
+% psitbx_str: source for the equilibrium map, psitbx-style
+% gdat_params: parameters for the gdat call % TODO: Why is it different from gdat_data.gdat_params?
+
+% Use psitbx
+t_th = gdat_data.t;
+th_z = mdsdata(['dim_of(\thomson' system_str_dot ':te,1)']); % TODO: Isn't this gdat_data.dim{1}?
+th_r = 0.9*ones(size(th_z));
+ntt = numel(t_th);
+nch = numel(th_z); % number of chords
+
+[t_psi,status] = mdsvalue('dim_of(tcv_eq("i_p",$))',psitbx_str); % TODO: replace with time_psi when tcv_eq supports it
+if ~rem(status,2)
+ warning('problem with getting equilibrium time base in gdat_tcv')
+ display(t_psi);
+ return
+end
+
+zshifteff=zshift;
+% set zshifteff time array same as t_psi
+switch numel(zshifteff)
+ case 1
+ zshifteff=zshifteff * ones(size(t_th));
+ case numel(t_psi)
+ zshifteff=interp1(t_psi,zshifteff,t_th);
+ case numel(t_th)
+ % ok
+ otherwise
+ if (gdat_params.nverbose>=1)
+ disp(' bad time dimension for zshift')
+ disp(['it should be 1 or ' num2str(numel(t_th)) ' or ' num2str(numel(t_psi))])
+ return
+ end
+end
+
+% Get LIUQE times corresponding to TS times
+t_eq = t_psi(iround(t_psi,t_th));
+
+% Get PSI map
+psi = psitbxtcv(gdat_data.shot,t_eq,'*0',psitbx_str);
+% PSITBXTCV will have removed duplicate times
+i_psi = iround(psi.t,t_eq); % Mapping from psi times to TS times
+
+% Define grid points where to evaluate PSI
+if ~any(diff(zshifteff))
+ % Constant vector
+ pgrid = psitbxgrid('cylindrical','points',{th_r,th_z - zshifteff(1),0});
+else
+ % We need as many time points in psi as TS times
+ psi = subs_time(psi,i_psi); % Now has ntt times
+ i_psi = 1:ntt; % Trivial mapping from psi times to TS times
+ % Time-varying vector
+ th_z_eff = repmat(th_z(:),1,ntt) - repmat(reshape(zshifteff,1,ntt),nch,1);
+ th_r_eff = repmat(th_r(:),1,ntt);
+ pgrid = psitbxgrid('cylindrical','time-points',{th_r_eff,th_z_eff,0});
+end
+
+% Compute psi at TS positions
+psi_ts = psitbxf2f(psi,pgrid);
+psiscatvol.data = squeeze(psi_ts.x(:,i_psi));
+psiscatvol.dim{1} = gdat_data.x;
+psiscatvol.dim{2} = t_th;
+% NOTE: we should probably not include time points where equilibrium time is far from TS time.
+
+% Compute psi_axis at TS times
+% Do not use psi_axis node because:
+% - For legacy LIUQE, psi_axis is not computed by psitbx, so we could get a
+% different answer and cause complex numbers computing rho_psi
+% - For MATLAB LIUQE, psi_axis is only the result of asxymex whereas the
+% psitbx adds some Newton iterations so again complex numbers are possible
+psi_norm = psitbxp2p(psi,'01');
+psi_max.data = psi_norm.psimag(i_psi);
+psi_max.dim = {t_th};
+