diff --git a/crpptbx_new/TCV/gdat_tcv.m b/crpptbx_new/TCV/gdat_tcv.m index d77241c779c88e1a592c8d1fb3a3d172804fb927..8e07dd683609b458c47c6c213631893d97b034b9 100644 --- a/crpptbx_new/TCV/gdat_tcv.m +++ b/crpptbx_new/TCV/gdat_tcv.m @@ -287,6 +287,7 @@ if strcmp(mapping_for_tcv.method(1:3),'tdi') end else if liuqe_version==-1 + mapping_for_tcv_expression_eff = mapping_for_tcv.expression; if strcmp(lower(mapping_for_tcv.expression(1:8)),'\results') mapping_for_tcv_expression_eff = mapping_for_tcv.expression(11:end); end @@ -401,7 +402,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % First the request names valid for "all" machines: % - case {'a_minor'} + case {'a_minor','rgeom'} + % compute average minor or major radius (on z=zaxis normally) nodenameeff=['\results::r_max_psi' substr_liuqe]; rmaxpsi=tdi(nodenameeff); nodenameeff2=['\results::r_min_psi' substr_liuqe]; @@ -410,8 +412,16 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') if ~isempty(ij); rmaxpsi.data(ij)=NaN; end ij=find(rminpsi.data<0.5 | rminpsi.data>1.2); if ~isempty(ij); rminpsi.data(ij)=NaN; end - gdat_data.data=0.5.*(rmaxpsi.data(end,:) - rminpsi.data(end,:)); - gdat_data.data_fullpath=[nodenameeff ' - ' nodenameeff2 ' /2']; + if strcmp(data_request_eff,'a_minor') + gdat_data.data=0.5.*(rmaxpsi.data(end,:) - rminpsi.data(end,:)); + gdat_data.data_fullpath=[nodenameeff ' - ' nodenameeff2 ' /2']; + elseif strcmp(data_request_eff,'rgeom') + gdat_data.data=0.5.*(rmaxpsi.data(end,:) + rminpsi.data(end,:)); + gdat_data.data_fullpath=[nodenameeff ' + ' nodenameeff2 ' /2']; + else + disp(['should not be in this case with data_request_eff = ' data_request_eff]) + return + end gdat_data.dim = rmaxpsi.dim(2); gdat_data.t = gdat_data.dim{1}; if any(strcmp(fieldnames(rmaxpsi),'units')) @@ -419,6 +429,24 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') end gdat_data.dimunits = rmaxpsi.dimunits(2); + case {'zgeom'} + % compute average minor or major radius (on z=zaxis normally) + nodenameeff=['\results::z_contour' substr_liuqe]; + zcontour=tdi(nodenameeff); + if strcmp(data_request_eff,'zgeom') + gdat_data.data=0.5.*(max(zcontour.data,[],1) + min(zcontour.data,[],1)); + gdat_data.data_fullpath=['(max+min)/2 of ' nodenameeff]; + gdat_data.dim{1} = zcontour.dim{2}; + gdat_data.dimunits{1} = zcontour.dimunits{2}; + else + disp(['should not be in this case with data_request_eff = ' data_request_eff]) + return + end + gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; + if any(strcmp(fieldnames(zcontour),'units')) + gdat_data.units = zcontour.units; + end + case {'b0'} % B0 at R0=0.88 R0EXP=0.88; @@ -474,7 +502,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.data = 100.*beta.data ./ abs(ip_t).*1.e6 .* abs(b0_t) .* a_minor_t; gdat_data.data_fullpath='100*beta/ip*1e6*b0*a_minor, each from gdat_tcv'; gdat_data.units = ''; - gdat_data.dimunits = 's'; + gdat_data.dimunits{1} = 's'; case {'cxrs'} %not yet finished, just started @@ -751,14 +779,18 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.ohm.data_fullpath=['-vloop(tens=' num2str(tension,'%.0e') ')*ip, from gdat']; gdat_data.ohm.label='P_{OHM}'; - % total power from all above - gdat_data.data = gdat_data.ec.data(:,10) + gdat_data.ohm.data; + % total power from each and total + gdat_data.data(:,1) = gdat_data.ohm.data; + gdat_data.data(:,2) = gdat_data.ec.data(:,10); + gdat_data.data(:,3) = gdat_data.ec.data(:,10) + gdat_data.ohm.data; + gdat_data.dim{2} = [1:3]; + gdat_data.dimunits{2} = 'Pohm;Pec;Ptot'; gdat_data.data_fullpath=['tot power from EC and ohm']; - gdat_data.label = 'P_{tot}=P_{EC}+P_{ohm}'; + gdat_data.label = 'P_{ohm};P_{EC};P_{tot}'; case {'q_rho'} % q profile on psi from liuqe - nodenameeff='\results::q_psi'; + nodenameeff=['\results::q_psi' substr_liuqe]; if liuqe_version==-1 nodenameeff=[begstr 'q_psi' substr_liuqe]; end @@ -770,11 +802,117 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') rhopol_eff = ones(size(tracetdi.dim{1})); rhopol_eff(:) = sqrt(linspace(0,1,length(tracetdi.dim{1}))); gdat_data.dim{1} = rhopol_eff; + gdat_data.x = gdat_data.dim{1}; gdat_data.dimunits{1} = ''; gdat_data.dimunits{2} = 's'; gdat_data.units = ''; gdat_data.request_description = 'q(rhopol\_norm)'; + case {'psi_edge'} + % psi at edge, 0 by construction in Liuqe, thus not given + nodenameeff=['\results::psi_axis' substr_liuqe]; + if liuqe_version==-1 + nodenameeff=[begstr 'q_psi' substr_liuqe]; + end + tracetdi=tdi(nodenameeff); + gdat_data.data = tracetdi.data.*0; + gdat_data.dim = tracetdi.dim; + gdat_data.t = gdat_data.dim{1}; + gdat_data.data_fullpath=[' zero ']; + gdat_data.dimunits = tracetdi.dimunits; + gdat_data.units = tracetdi.units; + gdat_data.request_description = '0 since LIUQE construct psi to be zero at LCFS'; + + case {'rhotor_edge','rhotor'} + % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper: + % O. Sauter, S.Yu. Medvedev, Comput. Phys. Commun. 184 (2013) 293–302 + % since cocos=17 for LIUQE we get: + % q = -dPhi/dpsi => Phi = - int(q*dpsi) which should always have the sign of B0 + params_eff = gdat_data.gdat_params; + params_eff.data_request='q_rho'; + q_rho=gdat_tcv([],params_eff); + params_eff.data_request='psi_axis'; % psi_edge=0 with LIUQE + psi_axis=gdat_tcv([],params_eff); + params_eff.data_request='b0'; % psi_edge=0 with LIUQE + b0=gdat_tcv([],params_eff); + b0tpsi = interp1(b0.t,b0.data,psi_axis.t); %q_rho on same time base as psi_axis + if isempty(psi_axis.data) || isempty(psi_axis.dim) || isempty(q_rho.data) || isempty(q_rho.dim) + warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]) + return + end + rhoequal = linspace(0,1,length(q_rho.dim{1})); + if strcmp(data_request,'rhotor_edge') + gdat_data.data = psi_axis.data; % to have the dimensions correct + gdat_data.dim = psi_axis.dim; + gdat_data.t = gdat_data.dim{1}; + 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,'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)'; + gdat_data.units = ''; + gdat_data.dimunits{1} = 'rhopol\_norm'; + gdat_data.dimunits{2} = 's'; + end + for it=1:length(psi_axis.data) + ij=find(~isnan(q_rho.data(:,it))); + if ~isempty(ij) + [qfit,~,~,phi]=interpos(q_rho.x(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 + if strcmp(data_request,'rhotor_edge') + gdat_data.data(it) = dataeff(end); + elseif strcmp(data_request,'rhotor') + gdat_data.data(:,it) = dataeff./dataeff(end); + gdat_data.rhotor_edge(it) = dataeff(end); + end + gdat_data.b0 = b0tpsi(it); + end + + case {'rhovol','volume_rho','volume'} + % volume_rho = vol(rho); volume = vol(LCFS) = vol(rho=1); + % rhovol = sqrt(vol(rho)/vol(rho=1)); + nodenameeff='\results::psitbx:vol'; + if liuqe_version==-1 + nodenameeff=[begstr 'vol' substr_liuqe]; + end + tracetdi=tdi(nodenameeff); + if isempty(tracetdi.data) || isempty(tracetdi.dim) + return + end + gdat_data.units = tracetdi.units; + if strcmp(data_request,'volume') + gdat_data.data = tracetdi.data(end,:); + gdat_data.dim{1} = tracetdi.dim{2}; + gdat_data.data_fullpath=['\results::psitbx:vol(end,:)']; + gdat_data.dimunits{1} = tracetdi.dimunits{2}; + gdat_data.request_description = 'volume(LCFS)=volume(rhopol=1)'; + else + gdat_data.data = tracetdi.data; + gdat_data.dim = tracetdi.dim; + gdat_data.dimunits = tracetdi.dimunits; + if strcmp(data_request,'volume_rho') + gdat_data.data_fullpath=['\results::psitbx:vol']; + gdat_data.request_description = 'volume(rho)'; + elseif strcmp(data_request,'rhovol') + gdat_data.volume_edge = gdat_data.data(end,:); + gdat_data.data = sqrt(gdat_data.data./repmat(reshape(gdat_data.volume_edge,1,size(gdat_data.data,2)),size(gdat_data.data,1),1)); + gdat_data.data_fullpath='sqrt(\results::psitbx:vol/vol_edge)'; + gdat_data.request_description = 'sqrt(volume(rho)/volume(edge))'; + else + disp(['should not be here in vol cases with data_request = ' data_request]); + return + end + end + gdat_data.t = gdat_data.dim{mapping_for_tcv.gdat_timedim}; + otherwise warning(['switchcase= ' data_request_eff ' not known in gdat_tcv']) error_status=901; diff --git a/crpptbx_new/TCV/tcv_requests_mapping.m b/crpptbx_new/TCV/tcv_requests_mapping.m index 5c40373ad8e585113d6a5da1d27b93e37d7e9039..ddbae25b6861b81698edfd01f56a1b97e0af3996 100644 --- a/crpptbx_new/TCV/tcv_requests_mapping.m +++ b/crpptbx_new/TCV/tcv_requests_mapping.m @@ -134,6 +134,14 @@ switch lower(data_request) mapping.timedim = 1; mapping.label = 'various powers'; mapping.method = 'switchcase'; + case 'psi_axis' + mapping.timedim = 1; + mapping.method = 'tdiliuqe'; + mapping.expression = '\results::psi_axis'; + case 'psi_edge' + mapping.timedim = 1; + mapping.method = 'switchcase'; % is set to zero, so not in tree nodes + mapping.label = 'psi\_edge'; case 'q0' mapping.timedim = 1; mapping.method = 'tdiliuqe'; @@ -191,7 +199,7 @@ switch lower(data_request) mapping.method = 'switchcase'; case 'vloop' mapping.timedim = 1; - mapping.label = ''; + mapping.label = 'Vloop'; mapping.method = 'tdi'; mapping.expression = [{'\magnetics::vloop[*,$1]'} {'001'}]; case 'volume' diff --git a/crpptbx_new/gdat_data_request_names_rho.xlsx b/crpptbx_new/gdat_data_request_names_rho.xlsx index 7816b8dab8164387a65273c5600312467169ac69..74f47e6263dadbe5c0879ab7e22d0c00441094c4 100644 Binary files a/crpptbx_new/gdat_data_request_names_rho.xlsx and b/crpptbx_new/gdat_data_request_names_rho.xlsx differ diff --git a/crpptbx_new/test_all_requestnames.m b/crpptbx_new/test_all_requestnames.m new file mode 100644 index 0000000000000000000000000000000000000000..ff449dbb6ebe8256394d963343b27e5fc153fd51 --- /dev/null +++ b/crpptbx_new/test_all_requestnames.m @@ -0,0 +1,12 @@ + +machine='TCV'; +shot=48836; +aa=gdat('machine',machine); +all_request_names = aa.gdat_request +%break +istart=1; +for irequest=istart:length(all_request_names) + request=all_request_names{irequest} + ab{irequest} = gdat(shot,request,'machine',machine,'doplot',1); + pause +end