From 4503a18bfe3bafbd348dfdbf48adba6891c8bb99 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Thu, 2 Jun 2016 16:48:28 +0000
Subject: [PATCH] add cxrs_rho, thus grids_1d, call a function for grids_1d,
add li, wmhd, psi_edge with vsurf for time variation. pgyro not added since
in powers.ec
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5870 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/TCV/gdat_tcv.m | 122 ++++++++++++++++++-----------
crpptbx/TCV/tcv_requests_mapping.m | 15 +++-
2 files changed, 89 insertions(+), 48 deletions(-)
diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
index 106b268f..6360949d 100644
--- a/crpptbx/TCV/gdat_tcv.m
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -576,7 +576,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.dimunits{1} = 's';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- case {'cxrs'}
+ case {'cxrs','cxrs_rho'}
%not yet finished, just started
% load typical data from cxrs, Ti, ni, vtori and vpoli (if available), as well as zeff from cxrs
% if 'fit' option is added: 'fit',1, then the fitted profiles are returned
@@ -709,7 +709,13 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
end
gdat_data.dimunits{1} = '';
gdat_data.dimunits{2} = 's';
-
+ % add grids_1d to have rhotor, etc if cxrs_rho was asked for
+ if strcmp(data_request_eff,'cxrs_rho')
+ gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);
+ else
+ gdat_data = get_grids_1d(gdat_data,2,0,gdat_params.nverbose);
+ end
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'eqdsk'}
%
@@ -936,30 +942,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.dim{1}=rho;
gdat_data.dimunits=[{'sqrt(psi_norm)'} ; {'time [s]'}];
gdat_data.x=rho;
- % add various rhos in grids_1d
- params_eff = gdat_data.gdat_params;params_eff.doplot=0;
- params_eff.data_request='rhotor';
- rhotor = gdat(gdat_data.shot,params_eff);
- params_eff.data_request='rhovol';
- rhovol = gdat(gdat_data.shot,params_eff);
- gdat_data.grids_1d.rhopolnorm = gdat_data.x;
- psi0 = interpos(rhotor.t',rhotor.psi_axis,gdat_data.t,-0.01);
- gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2.*repmat(reshape(psi0,1,length(psi0)),size(gdat_data.grids_1d.rhopolnorm,1),1);
- gdat_data.grids_1d.rhotornorm = NaN*ones(size(gdat_data.data));
- gdat_data.grids_1d.rhovolnorm = NaN*ones(size(gdat_data.data));
- it_rt = iround_os(rhotor.t,gdat_data.t);
- it_vol = iround_os(rhovol.t,gdat_data.t);
- for it=1:length(gdat_data.t)
- % do an interpolation on closest point to avoid 2D interp
- it_rt_eff = it_rt(it);
- it_vol_eff = it_vol(it);
- ii=find(~isnan(gdat_data.grids_1d.rhopolnorm(:,it)));
- if length(ii)==length(gdat_data.grids_1d.rhopolnorm)
- gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor.x,rhotor.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
- gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it));
- end
- end
- gdat_data.grids_1d.rhotor_edge=interpos(rhotor.t',rhotor.rhotor_edge,gdat_data.t',-0.1);
+ % add grids_1d to have rhotor, etc if cxrs_rho was asked for
+ gdat_data = get_grids_1d(gdat_data,2,1,gdat_params.nverbose);
%%%%%%%%%%%
% if nete_rho, copy data as .ne, get .te and put pe=e ne Te in data:
@@ -1307,31 +1291,13 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.dimunits{2} = 's';
gdat_data.units = '';
gdat_data.request_description = 'q(rhopol\_norm)';
- % add various rhos in grids_1d
- params_eff = gdat_data.gdat_params;params_eff.doplot=0;
- params_eff.data_request='rhotor';
- rhotor = gdat(gdat_data.shot,params_eff);
- params_eff.data_request='rhovol';
- rhovol = gdat(gdat_data.shot,params_eff);
- gdat_data.grids_1d.rhopolnorm = gdat_data.x;
- psi0 = interpos(rhotor.t',rhotor.psi_axis,gdat_data.t,-0.01);
- gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2*reshape(psi0,1,length(psi0));
- gdat_data.grids_1d.rhotornorm = NaN*ones(size(gdat_data.data));
- gdat_data.grids_1d.rhovolnorm = NaN*ones(size(gdat_data.data));
- it_rt = iround_os(rhotor.t,gdat_data.t);
- it_vol = iround_os(rhovol.t,gdat_data.t);
- for it=1:length(gdat_data.t)
- % do an interpolation on closest point to avoid 2D interp
- it_rt_eff = it_rt(it);
- it_vol_eff = it_vol(it);
- gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor.x,rhotor.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm);
- gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm);
- end
- gdat_data.grids_1d.rhotor_edge=interpos(rhotor.t',rhotor.rhotor_edge,gdat_data.t',-0.1);
+ % add grids_1d to have rhotor, etc if cxrs_rho was asked for
+ gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'psi_edge'}
% psi at edge, 0 by construction in Liuqe, thus not given
+ % add surface_psi from surface_flux and d(surface_flux)/dt = vloop
nodenameeff=['\results::psi_axis' substr_liuqe];
if liuqe_version==-1
nodenameeff=[begstr 'q_psi' substr_liuqe];
@@ -1349,6 +1315,13 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
gdat_data.dimunits = tracetdi.dimunits;
gdat_data.units = tracetdi.units;
gdat_data.request_description = '0 since LIUQE construct psi to be zero at LCFS';
+ %
+ params_eff.data_request='\results::surface_flux';
+ surface_psi=gdat_tcv([],params_eff);
+ [aa,vsurf] = interpos(surface_psi.t,surface_psi.data,-3);
+ gdat_data.surface_psi = surface_psi.data;
+ gdat_data.vsurf = vsurf;
+ gdat_data.vsurf_description = ['time derivative from surface_psi, obtained from \results::surface_flux'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'rhotor_edge','rhotor'}
@@ -1838,3 +1811,58 @@ gdat_data.t=time;
if any(strcmp(fieldnames(tracetdi),'units'))
gdat_data.units=tracetdi.units;
end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [gdat_data] = get_grids_1d(gdat_data,nbdim_x,nopt,nverbose);
+%
+% add various rhos in grids_1d
+%
+% nbdim_x = 1: same x (rhopol) for all times
+% = 2: rhopol depends on time like Thomson projection
+%
+% nopt = 0: do not fill in, just make the empty structure
+% = 1: do compute various fields
+%
+
+gdat_data.grids_1d.rhopolnorm = gdat_data.x;
+if (nopt == 0)
+ gdat_data.grids_1d.rhotornorm = [];
+ gdat_data.grids_1d.rhovolnorm = [];
+ gdat_data.grids_1d.psi = [];
+ gdat_data.grids_1d.rhotor_edge = [];
+ gdat_data.grids_1d.volume_edge = [];
+ return
+end
+
+params_eff = gdat_data.gdat_params;params_eff.doplot=0;
+params_eff.data_request='rhotor';
+rhotor = gdat(gdat_data.shot,params_eff);
+params_eff.data_request='rhovol';
+rhovol = gdat(gdat_data.shot,params_eff);
+psi0 = interpos(rhotor.t',rhotor.psi_axis,gdat_data.t,-0.01);
+if (nbdim_x == 1)
+ gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2*reshape(psi0,1,length(psi0));
+elseif (nbdim_x == 2)
+ gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2.*repmat(reshape(psi0,1,length(psi0)),size(gdat_data.grids_1d.rhopolnorm,1),1);
+else
+ if nverbose>=0; disp(['option: nbdim_x = ' numstr(nbdim_x) ' not implemented, check with O. Sauter']); end
+ return
+end
+gdat_data.grids_1d.rhotornorm = NaN*ones(size(gdat_data.data));
+gdat_data.grids_1d.rhovolnorm = NaN*ones(size(gdat_data.data));
+it_rt = iround_os(rhotor.t,gdat_data.t);
+it_vol = iround_os(rhovol.t,gdat_data.t);
+for it=1:length(gdat_data.t)
+ % do an interpolation on closest point to avoid 2D interp
+ it_rt_eff = it_rt(it);
+ it_vol_eff = it_vol(it);
+ ii=find(~isnan(gdat_data.grids_1d.rhopolnorm(:,it)));
+ if length(ii)==length(gdat_data.grids_1d.rhopolnorm)
+ gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor.x,rhotor.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+ gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+ end
+end
+gdat_data.grids_1d.rhotor_edge=interpos(rhotor.t',rhotor.rhotor_edge,gdat_data.t',-0.01);
+gdat_data.grids_1d.volume_edge=interpos(rhovol.t',rhovol.volume_edge,gdat_data.t',-0.01);
+
+return
diff --git a/crpptbx/TCV/tcv_requests_mapping.m b/crpptbx/TCV/tcv_requests_mapping.m
index 955aede6..b88a7986 100644
--- a/crpptbx/TCV/tcv_requests_mapping.m
+++ b/crpptbx/TCV/tcv_requests_mapping.m
@@ -61,6 +61,11 @@ switch lower(data_request)
mapping.label = 'cxrs';
mapping.method = 'switchcase';
mapping.expression = '';
+ case 'cxrs_rho'
+ mapping.timedim = 2;
+ mapping.label = 'cxrs';
+ mapping.method = 'switchcase';
+ mapping.expression = '';
case 'delta'
mapping.timedim = 1;
mapping.method = 'tdiliuqe';
@@ -89,7 +94,7 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.label = 'Halpha';
mapping.method = 'tdi';
- mapping.expression = '\base::pd:pd_011';
+ mapping.expression = '\base::pd:pd_001';
case 'ioh'
mapping.timedim = 1;
mapping.label = 'I ohmic transformer';
@@ -104,6 +109,10 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.method = 'tdiliuqe';
mapping.expression = '\results::kappa_edge';
+ case 'li'
+ mapping.timedim = 1;
+ mapping.method = 'tdiliuqe';
+ mapping.expression = '\results::l_i';
case 'mhd'
mapping.timedim = 1;
mapping.label = 'n=1,2, etc';
@@ -231,6 +240,10 @@ switch lower(data_request)
mapping.label = 'Volume(\rho)';
mapping.method = 'switchcase';
% mapping.expression = '\results::psitbx:vol'; (if exists for liuqe2 and 3 as well)
+ case 'wmhd'
+ mapping.timedim = 1;
+ mapping.method = 'tdiliuqe';
+ mapping.expression = '\results::total_energy';
case 'z_contour'
mapping.timedim = 1;
mapping.method = 'tdiliuqe';
--
GitLab