diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
index 0bd11540af6521d0ae6e6694716969711712ee13..ea72c068798591e6de9ab89293965e2160a2fdfa 100644
--- a/crpptbx/TCV/gdat_tcv.m
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -1533,19 +1533,56 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
phi_tor.units = 'Wb';
end
if ~any(any(isfinite(phi_tor.data)))
+ % no phi_tor, compute it from q profile
q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
if liuqe_matlab==0
q_rho.dim{1} = sqrt(linspace(0,1,numel(q_rho.dim{1})));
end
phi_tor.dim = q_rho.dim;
phi_tor.dimunits = q_rho.dimunits;
+ params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
+ psi_axis=gdat_tcv(shot,params_eff);
for it=1:size(q_rho.data,2)
ij=find(isfinite(q_rho.data(:,it)));
if ~isempty(ij) && numel(ij)>5
- [qfit,~,~,phi]=interpos(q_rho.dim{1}(ij),q_rho.data(ij,it),phi_tor.dim{1});
+ [qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),phi_tor.dim{1}.^2);
+ % Phi=sigma_bp*sigma_rhothetaphi*(2pi)^(1-eBp)*int(q dpsi), dpsi=dpsi_norm * (psiedge-psaxis)
+ % cocos=17: sigma_Bp=-1,sigma_rhothetaphi=1, eBp=1, (psiedge-psaxis)=-psiaxis
+ phi = -1.* phi .* (-psi_axis.data(it));
phi_tor.data(:,it) = phi';
end
end
+ elseif any(any(~isfinite(phi_tor.data)))
+ % there are some non-finite values, usually edge values, so replace them with integrated values
+ only_edge = 0;
+ if any(any(~isfinite(phi_tor.data(1:end-1,:))))
+ only_edge = 1;
+ end
+ q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
+ if liuqe_matlab==0
+ q_rho.dim{1} = sqrt(linspace(0,1,numel(q_rho.dim{1})));
+ end
+ params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE
+ psi_axis=gdat_tcv(shot,params_eff);
+ % assume same dimensions, check
+ if numel(phi_tor.data) ~= numel(q_rho.data) || numel(psi_axis.data) ~= numel(phi_tor.dim{2})
+ disp(['problems in gdat_tcv with ' data_request_eff ' with unexpected dimensions'])
+ return;
+ end
+ for it=1:size(q_rho.data,2)
+ %if any(~isfinite(phi_tor.data(:,it)))
+ ij=find(isfinite(q_rho.data(:,it)));
+ if ~isempty(ij) && numel(ij)>5
+ [qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),phi_tor.dim{1}.^2);
+ phi = -1.* phi .* (-psi_axis.data(it));
+ if only_edge
+ phi_tor.data(end,it) = phi(end);
+ else
+ phi_tor.data(:,it) = phi';
+ end
+ end
+ %end
+ end
end
gdat_data.data = phi_tor.data;
gdat_data.dim = phi_tor.dim;
@@ -1622,82 +1659,80 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.data_fullpath='phi from q_rho, psi_axis and integral(-q dpsi)';
gdat_data.units = 'T m^2';
gdat_data.dimunits{1} = 's';
- elseif strcmp(data_request_eff,'rhotor')
+ elseif strcmp(data_request_eff,'rhotor_norm') || strcmp(data_request_eff,'rhotor')
gdat_data.data = q_rho.data; % to have the dimensions correct
gdat_data.dim{1} = ones(size(q_rho.dim{1}));
gdat_data.dim{1}(:) = rhoequal;
gdat_data.dim{2} = q_rho.dim{2};
gdat_data.t = gdat_data.dim{2};
- gdat_data.data_fullpath='sqrt(phitor/phitor_edge), rhotor_edge=sqrt(phitor/B0/pi)';
+ if strcmp(data_request_eff,'rhotor_norm')
+ gdat_data.data_fullpath='sqrt(phitor/phitor_edge), rhotor_edge=sqrt(phitor(1)/B0/pi)';
+ else
+ gdat_data.data_fullpath='sqrt(phitor/pi/B0), rhotor_edge=sqrt(phitor(1)/B0/pi)';
+ end
gdat_data.units = '';
gdat_data.dimunits{1} = 'rhopol\_norm';
gdat_data.dimunits{2} = 's';
end
gdat_data.x=gdat_data.dim{1};
- gdat_data.psi_axis = reshape(psi_axis.data,1,length(psi_axis.data));
+ gdat_data.psi_axis = reshape(psi_axis.data,length(psi_axis.data),1);
if length(gdat_data.psi_axis) ~= length(gdat_data.t)
disp('problems of time between qrho and psi_axis?')
return
end
+ gdat_data.b0 = b0tpsi;
for it=1:size(q_rho.data,2)
ij=find(isfinite(q_rho.data(:,it)));
if ~isempty(ij)
- [qfit,~,~,phi]=interpos(q_rho.dim{1}(ij).^2,q_rho.data(ij,it),rhoequal.^2);
+ [qfit,~,~,phi]=interpos(-43,q_rho.dim{1}(ij).^2,q_rho.data(ij,it),rhoequal.^2);
dataeff = sqrt(phi .* psi_axis.data(it) ./ b0tpsi(it) ./ pi) ; % Delta_psi = -psi_axis
else
dataeff = NaN;
end
+ gdat_data.rhotor_edge(it) = dataeff(end); % redundant with "rhotor_edge" but to always have all subfields
if strcmp(data_request_eff,'rhotor_edge')
gdat_data.data(it) = dataeff(end);
- elseif strcmp(data_request_eff,'rhotor')
+ elseif strcmp(data_request_eff,'rhotor_norm')
gdat_data.data(:,it) = dataeff./dataeff(end);
- gdat_data.rhotor_edge(it) = dataeff(end);
+ elseif strcmp(data_request_eff,'rhotor')
+ gdat_data.data(:,it) = dataeff;
+ else
+ disp(['problem in gdat_tcv rhotor with unknown data_request_eff = ' data_request_eff]);
+ return
end
- gdat_data.b0 = b0tpsi(it);
end
+ gdat_data.rhotor_edge = gdat_data.rhotor_edge';
else
params_eff = gdat_data.gdat_params;
params_eff.data_request='phi_tor';
phi_tor=gdat_tcv([],params_eff);
-% $$$ if ~any(any(isfinite(phi_tor.data)))
-% $$$ % not yet implemented, compute phi_tor
-% $$$ q_rho = tdi(['tcv_eq(''''''' liuqefortran2liuqematlab('q',1,liuqe_matlab) ''''''',''''''' psitbx_str ''''''')']);
-% $$$ for it=1:size(q_rho.data,2)
-% $$$ ij=find(isfinite(q_rho.data(:,it)));
-% $$$ if ~isempty(ij) && numel(ij)>5
-% $$$ [qfit,~,~,phi]=interpos(q_rho.dim{1}(ij),q_rho.data(ij,it),phi_tor.dim{1});
-% $$$ phi_tor.data(:,it) = phi;
-% $$$ end
-% $$$ end
-% $$$ end
params_eff = gdat_data.gdat_params;
params_eff.data_request='b0';
b0=gdat_tcv([],params_eff);
gdat_data.t = phi_tor.t;
ij=find(isfinite(b0.data));
gdat_data.b0 = interpos(b0.t(ij),b0.data(ij),gdat_data.t);
+ gdat_data.rhotor_edge = sqrt(phi_tor.data(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)))';
if strcmp(data_request_eff,'rhotor_edge')
- gdat_data.data = sqrt(phi_tor.data(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)));
+ gdat_data.data = gdat_data.rhotor_edge;
gdat_data.dim = phi_tor.dim(2); % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.dimunits = phi_tor.dimunits{2};
gdat_data.units = 'm';
gdat_data.data_fullpath=['rhotor_edge=sqrt(Phi(edge)/pi/B0) from phi_tor'];
- elseif strcmp(data_request_eff,'rhotor')
- gdat_data.data = sqrt(phi_tor.data./pi./(ones(size(phi_tor.data,1),1)*reshape(gdat_data.b0,1,numel(gdat_data.b0))));
- gdat_data.rhotor_edge = gdat_data.data(end,:);
- gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
- gdat_data.x=gdat_data.dim{1};
- gdat_data.units = 'm';
- gdat_data.dimunits = phi_tor.dimunits;
- gdat_data.data_fullpath=['rhotor=sqrt(Phi/pi/B0) from phi_tor'];
elseif strcmp(data_request_eff,'rhotor_norm')
gdat_data.data = sqrt(phi_tor.data./(ones(size(phi_tor.data,1),1)*reshape(phi_tor.data(end,:),1,size(phi_tor.data,2))));
- gdat_data.rhotor_edge = sqrt(phi_tor.data(end,:)./pi./reshape(gdat_data.b0,1,numel(gdat_data.b0)));
gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
gdat_data.x=gdat_data.dim{1};
gdat_data.units = '';
gdat_data.dimunits = phi_tor.dimunits;
gdat_data.data_fullpath=['rhotor_norm=sqrt(Phi/Phi(edge)) from phi_tor'];
+ elseif strcmp(data_request_eff,'rhotor')
+ gdat_data.data = sqrt(phi_tor.data./pi./(ones(size(phi_tor.data,1),1)*reshape(gdat_data.b0,1,numel(gdat_data.b0))));
+ gdat_data.dim = phi_tor.dim; % while there is a problem with \tcv_shot::top.results.equil... dim nodes
+ gdat_data.x=gdat_data.dim{1};
+ gdat_data.units = 'm';
+ gdat_data.dimunits = phi_tor.dimunits;
+ gdat_data.data_fullpath=['rhotor=sqrt(Phi/pi/B0) from phi_tor'];
else
disp(['data_request_eff = ' data_request_eff ' not defined within rhotor block in gdat_tcv.m']);
return