From b103596e81e98ec3c1b7fe0028f4d6555e3e2b37 Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Fri, 11 Oct 2019 14:52:52 +0200
Subject: [PATCH] Create separate function to compute psi_scatvol with zshift.
---
matlab/TCV/gdat_tcv.m | 150 +++++++++++++++++++++++-------------------
1 file changed, 81 insertions(+), 69 deletions(-)
diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m
index b332400d..3c8d3fd6 100644
--- a/matlab/TCV/gdat_tcv.m
+++ b/matlab/TCV/gdat_tcv.m
@@ -1661,75 +1661,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
psiscatvol=gdat_tcv([],['\results::thomson' edge_str_dot ':psiscatvol' substr_liuqe],'nverbose',gdat_params.nverbose);
else
% calculate new psiscatvol
- % 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
- 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))])
- end
- end
-
-
- % 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} = 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};
+ [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;
@@ -3154,3 +3086,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};
+
--
GitLab