From f55087333a7c7849d9cc4440b7374e10d7b694df Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Mon, 13 May 2019 05:35:02 +0000
Subject: [PATCH] updates for running torbeam throughout
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11871 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/gdat_aug.m | 78 ++++++++++++++++++++++++++++++++-------
crpptbx/AUG/geteqdskAUG.m | 5 ++-
crpptbx/get_grids_1d.m | 33 ++++++++++++-----
3 files changed, 92 insertions(+), 24 deletions(-)
diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
index a08ab492..172c6e05 100644
--- a/crpptbx/AUG/gdat_aug.m
+++ b/crpptbx/AUG/gdat_aug.m
@@ -728,9 +728,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.fit.t = gdat_data.t;
for it=1:length(gdat_data.t)
% make rhotor->rhopol transformation for each time since equilibrium might have changed
- irhook=find(gdat_data.rhotornorm(:,it)>0 & gdat_data.rhotornorm(:,it)<1); % no need for ~isnan
- [rhoeff isort]=sort(gdat_data.rhotornorm(irhook,it));
- gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.rhopolnorm(irhook(isort),it); 1],rhotornormfit,-0.1,[2 2],[0 1]);
+ gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
+ gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
idata = find(gdat_data.ti.data(:,it)>0 & gdat_data.rhotornorm(:,it)<1.01);
if length(idata)>0
gdat_data.fit.ti.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it);
@@ -956,6 +955,13 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.gdat_params.zshift = zshift;
end
gdat_data.gdat_params.zshift = zshift;
+ zwrite = 0.;
+ if isfield(gdat_data.gdat_params,'write') && ~isempty(gdat_data.gdat_params.write)
+ zwrite = gdat_data.gdat_params.write;
+ else
+ gdat_data.gdat_params.write = zwrite;
+ end
+ gdat_data.gdat_params.write = zwrite;
params_equil = gdat_data.gdat_params;
params_equil.data_request = 'equil';
[equil,params_equil,error_status] = gdat_aug(shot,params_equil);
@@ -965,10 +971,12 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
end
gdat_data.gdat_params = params_equil;
gdat_data.equil = equil;
+ figure(10);
+ fignb_handle = gcf;
for itime=1:length(time_eqdsks)
time_eff = time_eqdsks(itime);
% use read_results updated to effectively obtain an eqdsk with sign correct with COCOS=2
- [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(equil,time_eff,gdat_data.gdat_params.zshift,'source',gdat_data.gdat_params.equil,'extra_arg_sf2sig',gdat_data.gdat_params.extra_arg_sf2sig);
+ [eqdskAUG, equil_all_t, equil_t_index]=geteqdskAUG(equil,time_eff,gdat_data.gdat_params.zshift,'source',gdat_data.gdat_params.equil,'extra_arg_sf2sig',gdat_data.gdat_params.extra_arg_sf2sig,'fignb',fignb_handle);
eqdskAUG.fnamefull = fullfile(['/tmp/' getenv('USER')],['EQDSK_' num2str(shot) 't' num2str(time_eff,'%.4f')]);
cocos_out = equil.cocos;
if isfield(gdat_data.gdat_params,'cocos') && ~isempty(gdat_data.gdat_params.cocos)
@@ -985,7 +993,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
% cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times,
% so project psi(R,Z) on Rmesh, Zmesh of 1st time
if length(time_eqdsks) > 1
- gdat_data.eqdsk{itime} = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout);
+ gdat_data.eqdsk{itime} = eqdsk_cocosout;
+ if gdat_data.gdat_params.write; gdat_data.eqdsk{itime} = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout); end
if itime==1
gdat_data.data(:,:,itime) = gdat_data.eqdsk{itime}.psi;
gdat_data.dim{1} = gdat_data.eqdsk{itime}.rmesh;
@@ -998,7 +1007,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.data(:,:,itime) = aa;
end
else
- gdat_data.eqdsk = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout);
+ gdat_data.eqdsk = eqdsk_cocosout;
+ if gdat_data.gdat_params.write; gdat_data.eqdsk = write_eqdsk(eqdskAUG.fnamefull,eqdsk_cocosout); end
gdat_data.data = gdat_data.eqdsk.psi;
gdat_data.dim{1} = gdat_data.eqdsk.rmesh;
gdat_data.dim{2} = gdat_data.eqdsk.zmesh;
@@ -1454,7 +1464,6 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.data_fullpath=['pe=1.6e-19*ne*Te in data, .ne, .te ' gdat_data.data_fullpath(3:end)];
gdat_data.label = 'pe';
end
-
% defaults for fits, so user always gets std structure
gdat_data.fit.rhotornorm = []; % same for both ne and te
gdat_data.fit.rhopolnorm = [];
@@ -1511,9 +1520,10 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
end
for it=1:length(gdat_data.t)
% make rhotor->rhopol transformation for each time since equilibrium might have changed
- irhook=find(gdat_data.grids_1d.rhotornorm(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<1); % no need for ~isnan
- [rhoeff isort]=sort(gdat_data.grids_1d.rhotornorm(irhook,it));
- gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.grids_1d.rhopolnorm(irhook(isort),it); 1],rhotornormfit,-0.1,[2 2],[0 1]);
+ gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
+ gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
+ gdat_data.fit.rhopolnorm(:,it)=interpos(gdat_data.grids_1d.rhotornorm_equil(:,it), ...
+ gdat_data.grids_1d.rhopolnorm_equil(:,it),rhotornormfit);
if any(strfind(data_request_eff,'te'))
idatate = find(gdat_data.te.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
if length(idatate)>0
@@ -1526,11 +1536,32 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
te_err_eff = gdat_data.te.error_bar(idatate(irhoeff),it);
% they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
ij=find(teeff./te_err_eff>10.);
- if ~isempty(ij); te_err_eff(ij) = teeff(ij)./0.1; end
+ if ~isempty(ij); te_err_eff(ij) = max(max(te_err_eff),teeff(ij)./0.1); end
%
teeff = [teeff(1); teeff];
te_err_eff = [1e4; te_err_eff];
+ % add some points to reach 0 at 1.05
+ if max(rhoeff)<1.03
+ rhoeff = [rhoeff; 1.05];
+ teeff = [teeff; 0.];
+ te_err_eff(end) = 10.*max(te_err_eff(2:end));
+ te_err_eff = [te_err_eff; min(te_err_eff)]; % to give more weight for new last point
+ else
+ teeff(end) = 0.;
+ te_err_eff(end-1) = 10.*max(te_err_eff(2:end));
+ te_err_eff(end) = min(te_err_eff); % to give more weight for new value of last point
+ end
[gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhotornorm(:,it)] = interpos(rhoeff,teeff,rhotornormfit,fit_tension.te,[1 0],[0 0],te_err_eff);
+ % avoid neg points or positive edge gradient
+ if ~isempty(find(gdat_data.fit.te.data(:,it)<0)) || gdat_data.fit.te.drhotornorm(end,it) > 0
+ ij= find(rhoeff>=0.85 & rhoeff<1.04);
+ if ~isempty(ij)
+ te_err_eff(ij) = 100.*min(te_err_eff);
+ else
+ te_err_eff(end-min(5,length(te_err_eff)-3):end-1) = 100.*min(te_err_eff);
+ end
+ [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhotornorm(:,it)] = interpos(rhoeff,teeff,rhotornormfit,fit_tension.te,[1 0],[0 0],te_err_eff);
+ end
end
end
if any(strfind(data_request_eff,'ne'))
@@ -1545,11 +1576,32 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
neeff = gdat_data.ne.data(idatane(irhoeff),it);
ne_err_eff = gdat_data.ne.error_bar(idatane(irhoeff),it);
ij=find(neeff./ne_err_eff>10.);
- if ~isempty(ij); ne_err_eff(ij) = neeff(ij)./0.1; end
+ if ~isempty(ij); ne_err_eff(ij) = max(max(ne_err_eff),neeff(ij)./0.1); end
%
neeff = [neeff(1); neeff];
ne_err_eff = [1e21; ne_err_eff];
+ % add some points to reach 0 at 1.05
+ if max(rhoeff)<1.03
+ rhoeff = [rhoeff; 1.05];
+ neeff = [neeff; 0.];
+ ne_err_eff(end) = 10.*max(ne_err_eff(2:end)); % to give more weight for new last point
+ ne_err_eff = [ne_err_eff; min(ne_err_eff)];
+ else
+ neeff(end) = 0.;
+ ne_err_eff(end-1) = 10.*max(ne_err_eff(2:end));
+ ne_err_eff(end) = min(ne_err_eff); % to give more weight for new value of last point
+ end
[gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhotornorm(:,it)] = interpos(rhoeff,neeff,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff);
+ % avoid neg points or positive edge gradient
+ if ~isempty(find(gdat_data.fit.ne.data(:,it)<0)) || gdat_data.fit.ne.drhotornorm(end,it) > 0
+ ij= find(rhoeff>=0.85 & rhoeff<1.04);
+ if ~isempty(ij)
+ ne_err_eff(ij) = 100.*min(ne_err_eff);
+ else
+ ne_err_eff(end-min(5,length(ne_err_eff)-3):end-1) = 100.*min(ne_err_eff);
+ end
+ [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhotornorm(:,it)] = interpos(rhoeff,neeff,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff);
+ end
end
end
end
@@ -2193,7 +2245,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
if strcmp(data_request_eff,'volume_rho')
gdat_data.data = gdat_data.vol;
gdat_data.units = Vol.units;
- else strcmp(data_request_eff,'rhovol')
+ else strcmp(data_request_eff,'rhovol');
gdat_data.data = gdat_data.rhovolnorm;
gdat_data.units = ' ';
end
diff --git a/crpptbx/AUG/geteqdskAUG.m b/crpptbx/AUG/geteqdskAUG.m
index d991b104..7abfccb8 100644
--- a/crpptbx/AUG/geteqdskAUG.m
+++ b/crpptbx/AUG/geteqdskAUG.m
@@ -135,10 +135,13 @@ zmag=gdat_aug(shot,'zmag','source',equilpar_source,'extra_arg_sf2sig',extra_arg_
eqdsk.zaxis = zmag.data(itrmag) - eqdsk.zshift;
% get plasma boundary
-if isempty(fignb) || ~isnumeric(fignb) || fignb < 1
+if isempty(fignb) || (isnumeric(fignb) && fignb < 1) || (~isnumeric(fignb) && ~ishandle(fignb))
figure
+ fignb_handle = gcf;
+ set(gcf,'Name','from geteqdskAUG');
else
figure(fignb)
+ fignb_handle = gcf;
clf
set(gcf,'Name','from geteqdskAUG');
end
diff --git a/crpptbx/get_grids_1d.m b/crpptbx/get_grids_1d.m
index 96c59a1f..e6a14999 100644
--- a/crpptbx/get_grids_1d.m
+++ b/crpptbx/get_grids_1d.m
@@ -27,7 +27,6 @@ if (nopt == 0) || isempty(gdat_data.x) || isempty(gdat_data.t) || isempty(gdat_d
gdat_data.grids_1d.b0 = [];
return
end
-
params_eff = gdat_data.gdat_params;params_eff.doplot=0;
params_eff.data_request='rhotor_norm';
rhotor_norm = gdat(gdat_data.shot,params_eff);
@@ -35,6 +34,12 @@ ndim_x_rhotor = length(find(size(rhotor_norm.x)>1));
params_eff.data_request='rhovol';
rhovol = gdat(gdat_data.shot,params_eff);
ndim_x_rhovol = length(find(size(rhotor_norm.x)>1));
+% check that rhotor and rhovol come from same equil psi
+if sum(abs(size(rhotor_norm.x)-size(rhovol.x)))~=0 || sum(sum(abs(rhotor_norm.x-rhovol.x)))>1e-10 ...
+ || sum(abs(rhotor_norm.t-rhovol.t))>1e-10;
+ error(['get_grids_1d: x and t arrays for rhotor_norm and rhovol are different, although should refer to same equilibrium' char(10) ...
+ 'size(rhotor_norm.x) = ' num2str(size(rhotor_norm.x)) '; size(rhovol.x) = ' num2str(size(rhovol.x))]);
+end
params_eff.data_request='psi_axis';
psi_axis = gdat(gdat_data.shot,params_eff);
params_eff.data_request='psi_edge';
@@ -62,27 +67,35 @@ 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);
+ if ndim_x_rhotor==1
+ gdat_data.grids_1d.rhopolnorm_equil(:,it) = rhotor_norm.x; % to have same sizes for all 3 reference rhos from equil
+ else
+ gdat_data.grids_1d.rhopolnorm_equil(:,it) = rhotor_norm.x(:,it_rt_eff);
+ end
+ gdat_data.grids_1d.rhotornorm_equil(:,it) = rhotor_norm.data(:,it_rt_eff);
+ % can assume same .x for rhotor and rhovol since called with same params thus equil (checked with sum before)
+ gdat_data.grids_1d.rhovolnorm_equil(:,it) = rhovol.data(:,it_vol_eff);
+ %
if (nbdim_x == 1)
ii=find(isfinite(gdat_data.grids_1d.rhopolnorm));
else
ii=find(isfinite(gdat_data.grids_1d.rhopolnorm(:,it)));
end
if (nbdim_x == 1)
-% $$$ if length(ii)==length(gdat_data.grids_1d.rhopolnorm)
if length(ii) >0
nb_ii = length(ii);
if ~isempty(rhotor_norm.x) && ~isempty(rhotor_norm.data)
if ndim_x_rhotor==1
- gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii));
+ gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii));
else
- gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii));
+ gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii));
end
end
if ~isempty(rhovol.x) && ~isempty(rhovol.data)
if ndim_x_rhovol==1
- gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii));
+ gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii));
else
- gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(-13,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii));
+ gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii));
end
end
end
@@ -92,16 +105,16 @@ for it=1:length(gdat_data.t)
if ~isempty(rhotor_norm.x) && ~isempty(rhotor_norm.data)
nb_ii = length(ii);
if ndim_x_rhotor==1
- gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
+ gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
else
- gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
+ gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
end
end
if ~isempty(rhovol.x) && ~isempty(rhovol.data)
if ndim_x_rhovol==1
- gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
+ gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
else
- gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(-13,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
+ gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
end
end
end
--
GitLab