From 969009fc3a5f57c1758007d74a3c886a0c547f67 Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Fri, 11 Oct 2019 11:29:13 +0200
Subject: [PATCH] Use psitbx to compute psi at the TS volume positions.
---
matlab/TCV/gdat_tcv.m | 94 +++++++++++++++++++++++++++----------------
1 file changed, 60 insertions(+), 34 deletions(-)
diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index 516e2f19..165876bd 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1652,57 +1652,83 @@ 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
+ % Use psitbx
+ t_th = gdat_data.t;
+ th_z = mdsdata(['dim_of(\thomson' edge_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
- rmesh=psitdi.dim{1};
- zmesh=psitdi.dim{2};
- zthom=mdsdata(['dim_of(\thomson' edge_str_dot ':te,1)']);
zshifteff=zshift;
- % set zshifteff time array same as psitdi
- switch length(zshifteff)
+ % set zshifteff time array same as t_psi
+ switch numel(zshifteff)
case 1
- zshifteff=zshifteff * ones(size(psitdi.dim{3}));
- case length(psitdi.dim{3})
+ zshifteff=zshifteff * ones(size(t_th));
+ case numel(t_psi)
+ zshifteff=interp1(t_psi,zshifteff,t_th);
+ case numel(t_th)
% ok
- case length(gdat_data.t)
- zshifteff=interp1(gdat_data.t,zshifteff,psitdi.dim{3});
otherwise
- if (gdat_params.nverbose>=1);
+ 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}))])
+ disp(['it should be 1 or ' num2str(numel(t_th)) ' or ' num2str(numel(t_psi))])
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-zshifteff(itpsitdi),'spline'); % faster with interpos
- psiscatvol0=interpos2Dcartesian(rmesh,zmesh,psirz,0.9*ones(size(zthom)),zthom-zshifteff(itpsitdi),-0.1,-0.1);
- %psiscatvol0=interp2(rmesh,zmesh,psirz',0.9*ones(size(zthom)),zthom-zshifteff(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;
+
+
+ % Get LIUQE times corresponding to TS times
+ t_eq = t_psi(iround(t_psi,t_th));
+
+ % Get PSI map
+ psi = psitbxtcv(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} = t_th;
psiscatvol.dim{2} = gdat_data.x;
+ % 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);
end
if ~isempty(psiscatvol.data) && ~ischar(psiscatvol.data) && ~isempty(psi_max.data) && ~ischar(psi_max.data)
for ir=1:length(psiscatvol.dim{2})
--
GitLab