From d577e1d130228a8c1895c7157825c9ed5eb273d7 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Tue, 4 Jun 2019 08:07:21 +0000
Subject: [PATCH] update further torbeam run start to have a shotfile
substructure
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11965 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/aug_help_parameters.m | 1 +
crpptbx/AUG/gdat_aug.m | 82 ++++++-
crpptbx/AUG/geteqdskAUG.m | 11 +-
crpptbx/AUG/rdaAUG_eff.m | 48 +++-
crpptbx/AUG/read_ech_py.m | 232 +++++++++++++------
crpptbx/AUG/torbeam_plot_output_singlerun.m | 5 +-
crpptbx/AUG/torbeam_prepare_inputs_and_run.m | 95 ++++++--
7 files changed, 369 insertions(+), 105 deletions(-)
diff --git a/crpptbx/AUG/aug_help_parameters.m b/crpptbx/AUG/aug_help_parameters.m
index dde25b8d..f936f614 100644
--- a/crpptbx/AUG/aug_help_parameters.m
+++ b/crpptbx/AUG/aug_help_parameters.m
@@ -43,6 +43,7 @@ help_struct_all.fit_nb_rho_points = 'nb of points for the radial mesh over which
help_struct_all.source = ['sxr: ''G'' (default, with ssx), camera name ''J'', ''G'', ...[F-M], case insensitive;' char(10) ...
'cxrs: ''CEZ'' (default), ''CMZ'',''CUZ'',''COZ'',''all'';' char(10) ...
'raptor: ''observer'', ''predictive'' (or ''obs'', ''pre'') to restrict the ouput to these signals'];
+help_struct_all.source_exp_name = ['ne_rho, te_rho, nete_rho for fit from TRA source typically or IDA when expname not the same as formain signal'];
help_struct_all.camera = ['[] (default, all), [i1 i2 ...] chord nbs ([1 3 5] if only chords 1, 3 and 5 are desired), ''central'' uses J_049'];
help_struct_all.freq = '''slow'', default (means ssx, 500kHz), lower sampling; ''fast'' full samples 2MHz; integer value nn for downsampling every nn''th points';
%help_struct_all. = '';
diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
index c281799c..02e766b8 100644
--- a/crpptbx/AUG/gdat_aug.m
+++ b/crpptbx/AUG/gdat_aug.m
@@ -227,6 +227,9 @@ if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'ide')
|| any(strcmp(lower(gdat_params.source),'ida')))
gdat_params.equil = 'IDE';
end
+% $$$ if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'tre')) || any(strcmp(lower(gdat_params.source),'tra')))
+% $$$ gdat_params.equil = 'TRE';
+% $$$ end
extra_arg_sf2sig = '[]';
if isfield(gdat_params,'extra_arg_sf2sig') && ~isempty(gdat_params.extra_arg_sf2sig)
extra_arg_sf2sig = gdat_params.extra_arg_sf2sig;
@@ -314,9 +317,13 @@ if strcmp(mapping_for_aug.method,'signal')
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~ischar(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = mapping_for_aug.expression{1};
else
- % special case for IDE instead of IDG when FPG is present value
+ % special case for EQI,EQH,IDE instead of IDG when FPG is present value
if strcmp(lower(gdat_data.gdat_params.source),'ide') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
gdat_data.gdat_params.source = 'IDG';
+ elseif strcmp(lower(gdat_data.gdat_params.source),'eqi') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
+ gdat_data.gdat_params.source = 'GQI';
+ elseif strcmp(lower(gdat_data.gdat_params.source),'eqh') && strcmp(lower(mapping_for_aug.expression{1}),'fpg')
+ gdat_data.gdat_params.source = 'GQH';
end
mapping_for_aug.expression{1} = gdat_data.gdat_params.source;
end
@@ -1156,7 +1163,11 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.jpol = [];
gdat_data.djpolpsi = [];
end
- gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
+ if ~strcmp(DIAG,'TRE')
+ gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)';
+ else
+ gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)'/2/pi;
+ end
gdat_data.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1);
gdat_data.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1);
gdat_data.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1);
@@ -1189,10 +1200,15 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
end
gdat_data.x = gdat_data.rhopolnorm;
%
- [equil_Rcoil,e]=rdaAUG_eff(shot,DIAG,'Rcl',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
- gdat_data.Rcoils=equil_Rcoil.value;
- [equil_Zcoil,e]=rdaAUG_eff(shot,DIAG,'Zcl',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
- gdat_data.Zcoils=equil_Zcoil.value;
+ try
+ [equil_Rcoil,e]=rdaAUG_eff(shot,DIAG,'Rcl',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
+ gdat_data.Rcoils = equil_Rcoil.value;
+ [equil_Zcoil,e]=rdaAUG_eff(shot,DIAG,'Zcl',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
+ gdat_data.Zcoils = equil_Zcoil.value;
+ catch
+ gdat_data.Rcoils = [];
+ gdat_data.Zcoils = [];
+ end
%
if strcmp(DIAG,'IDE')
[IpiPSI,e]=rdaAUG_eff(shot,'IDG','Itor',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig);
@@ -1207,7 +1223,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
%
gdat_data.data = gdat_data.qvalue; % put q in data
gdat_data.units=[]; % not applicable
- gdat_data.data_fullpath = [DIAG ' from expname: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
+ gdat_data.data_fullpath = [DIAG ' from exp_name: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
gdat_data.cocos = 17; % should check FPP
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
@@ -1335,6 +1351,11 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
gdat_data.gdat_params.source = 'VTA'; % default source
end
+ sources_exp_name_available = 'AUGD';
+ if isfield(gdat_data.gdat_params,'source_exp_name') && ~isempty(gdat_data.gdat_params.source_exp_name)
+ sources_exp_name_available = gdat_data.gdat_params.source_exp_name;
+ end
+ gdat_data.gdat_params.source_exp_name = sources_exp_name_available;
gdat_data.gdat_params.source = upper(gdat_data.gdat_params.source);
if strcmp(gdat_data.gdat_params.source(1:2),'EQ');
warning('set equilibrium choice in ''equil'' parameter and not source, moved to equil')
@@ -1663,7 +1684,41 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.fit.ne.error_bar = gdat_data_ne_ida_err.data; % time as 2nd dim
end
case 'TRA'
+ gdat_data.fit.full_path = 'from TRA: TE, NE, rhop from PLFLX and XB (.x is rhotornorm mid-points, XB end-point)';
+ % (tra.XB.data(1:end,it)'+[0 ,tra.XB.data(1:end-1,it)'])/2 = tra.TE.x(1:end,it)'
+ if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'te')
+ % Te
+ params_eff.data_request={'TRA','TE',gdat_data.gdat_params.source_exp_name};
+ [gdat_data_te_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
+ gdat_data.fit.te.data = gdat_data_te_tra.data;
+ gdat_data.fit.t = gdat_data_te_tra.t;
+ gdat_data.fit.rhotornorm = gdat_data_te_tra.x;
+ params_eff.data_request={'TRA','PLFLX',gdat_data.gdat_params.source_exp_name};
+ [gdat_data_plflx_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
+ rhopol_endpoints = sqrt((gdat_data_plflx_tra.data-0)./(repmat(gdat_data_plflx_tra.data(end,:),size(gdat_data_plflx_tra.data,1),1)-0.));
+ for it=1:size(gdat_data_plflx_tra.data,2)
+ gdat_data.fit.rhopolnorm(:,it) = interpos([0. ; gdat_data_plflx_tra.x(:,it)],[0. ; rhopol_endpoints(:,it)],gdat_data.fit.rhotornorm(:,it));
+ end
+ end
+ if strcmp(data_request_eff(1:4),'nete') || strcmp(data_request_eff(1:2),'ne')
+ % ne
+ params_eff.data_request={'TRA','NE',gdat_data.gdat_params.source_exp_name};
+ [gdat_data_ne_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
+ gdat_data.fit.ne.data = gdat_data_ne_tra.data .* 1e6; % ne in cm^3 !!
+ if ~strcmp(data_request_eff(1:4),'nete')
+ gdat_data.fit.t = gdat_data.t;
+ gdat_data.fit.rhotornorm = gdat_data_ne_tra.x;
+ params_eff.data_request={'TRA','PLFLX',gdat_data.gdat_params.source_exp_name};
+ [gdat_data_plflx_tra,params_kin,error_status]=gdat_aug(shot,params_eff);
+ rhopol_endpoints = sqrt((gdat_data_plflx_tra.data-0)./(repmat(gdat_data_plflx_tra.data(end,:),size(gdat_data_plflx_tra.data,1),1)-0.));
+ for it=1:size(gdat_data_plflx_tra.data,2)
+ gdat_data.fit.rhopolnorm(:,it) = interpos([0. ; gdat_data_plflx_tra.x(:,it)],[0. ; rhopol_endpoints(:,it)],gdat_data.fit.rhotornorm(:,it));
+ end
+ end
+ end
otherwise
+ disp(['profiles ''fit'' from ' gdat_data.gdat_params.source ' not defined. At this stage only: ' sources_available{:} ...
+ ' are available. Ask O. Suater if need be']);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -2097,7 +2152,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
% get time values
gdat_data.data = gdat_data.qvalue; % put q in data
gdat_data.units=[]; % not applicable
- gdat_data.data_fullpath = [DIAG ' from expname: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
+ gdat_data.data_fullpath = [DIAG ' from exp_name: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
gdat_data.cocos = 17; % should check FPP
gdat_data.dim{1} = gdat_data.x;
gdat_data.dim{2} = gdat_data.t;
@@ -2476,6 +2531,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
% most of the times the exp for the shotfile should be provided
shotfile_exp_eff = gdat_params.exp_name;
diagname='TRA';
+ gdat_params.equil = 'TRE';
+ gdat_data.gdat_params = gdat_params;
TRANSP_signals;
for i=1:size(transp_sig,1)
if strcmp(lower(transp_sig{i,2}),'signal') || strcmp(lower(transp_sig{i,2}),'signal-group')
@@ -2490,6 +2547,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
clear adata_area
try
eval(['[adata_area]=sf2ab(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']);
+ adata_area.data = adata_area.value{1};
catch
adata_area.value = cell(0);
end
@@ -2498,8 +2556,9 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
clear adata_time
try
eval(['[adata_time]=sf2tb(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff,shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']);
+ adata_time.data = adata_time.value;
catch
- adata_time.value = cell(0);
+ adata_time.value = [];
end
eval(['gdat_data.' transp_sig{i,1} '=adata_time;']);
end
@@ -2511,6 +2570,11 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.dimunits{1} = gdat_data.TIME.unit;
end
+ % add equilibrium
+ params_eff = gdat_params;
+ params_eff.data_request = 'equil';
+ gdat_data.equil = gdat_aug(shot,params_eff);
+
otherwise
if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_aug']); end
error_status=901;
diff --git a/crpptbx/AUG/geteqdskAUG.m b/crpptbx/AUG/geteqdskAUG.m
index 7abfccb8..f0624981 100644
--- a/crpptbx/AUG/geteqdskAUG.m
+++ b/crpptbx/AUG/geteqdskAUG.m
@@ -45,7 +45,10 @@ else
equil_source = 'EQI';
end
equilpar_source = 'FPG';
-if strcmp(lower(equil_source),'ide'); equilpar_source = 'IDG'; end
+if strcmp(upper(equil_source),'IDE'); equilpar_source = 'IDG'; end
+if strcmp(upper(equil_source),'EQI'); equilpar_source = 'GQI'; end
+if strcmp(upper(equil_source),'EQH'); equilpar_source = 'GQH'; end
+if strcmp(upper(equil_source),'TRA'); equilpar_source = 'TRE'; end
if nargin >= 7 && ~isempty(varargin{3}) && strcmp(lower(varargin{3}),'extra_arg_sf2sig')
if ~isempty(varargin{4})
@@ -112,7 +115,11 @@ if abs(equil.qvalue(end,it)) > 25
else
eqdsk.q = interpos(equil.rhopolnorm(:,it),equil.qvalue(:,it),eqdsk.rhopsi,-0.01,[1 2],[0 equil.qvalue(end,it)]);
end
-
+if strcmp(upper(equil_source),'IDE');
+ % q of IDE ids abs(q), add sign from EQI
+ q95_eqi=gdat_aug(shot,{'GQI','q95'},'extra_arg_sf2sig',extra_arg_sf2sig);
+ eqdsk.q = abs(eqdsk.q) * sign(q95_eqi.data(round(length(q95_eqi.data)/2)));
+end
eqdsk.ip = equil.Ip(it);
eqdsk.stitle=['#' num2str(shot) ' t=' num2str(time_eff) ' from ' equil.gdat_params.exp_name '/' equil.gdat_params.equil];
eqdsk.ind1=1;
diff --git a/crpptbx/AUG/rdaAUG_eff.m b/crpptbx/AUG/rdaAUG_eff.m
index 558f9907..954f4d66 100644
--- a/crpptbx/AUG/rdaAUG_eff.m
+++ b/crpptbx/AUG/rdaAUG_eff.m
@@ -117,6 +117,7 @@ if nargin_eff>=7 && ~isempty(varargin_eff{3}) && ischar(varargin_eff{3}) ...
end
end
+did_transpose = 0;
if usemdsplus
% a_remote=mdsremotelist; does not seem sufficient if connection broken, test with 1+2
is3 = mdsvalue('1+2');
@@ -247,16 +248,24 @@ if usemdsplus
eval(['tunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',' num2str(idim0) '))'');']);
adata.dimunits = {tunits};
case 2
- adata.data = adata.data';
- did_transpose = 1;
+ if strcmp(upper(diagname),'IDA')
+ % rho, time in original dimension, so do not transpose
+ idim_x = 0;
+ idim_t = 1;
+ else
+ adata.data = adata.data';
+ did_transpose = 1;
+ idim_x = 1;
+ idim_t = 0;
+ end
% transposed because of C relation and backward compatibility with sf2sig part
- eval(['x=mdsvalue(''dim_of(_rdaeff' user diagname ',1)'');']);
+ eval(['x=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(idim_x) ')'');']);
if prod(size(x))==length(x); x = reshape(x,1,length(x)); end
- eval(['time=mdsvalue(''dim_of(_rdaeff' user diagname ',0)'');']);
+ eval(['time=mdsvalue(''dim_of(_rdaeff' user diagname ',' num2str(idim_t) ')'');']);
time = reshape(time,1,length(time));
adata.dim = {x, time};
- eval(['xunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',1))'');']);
- eval(['tunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',0))'');']);
+ eval(['xunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',' num2str(idim_x) '))'');']);
+ eval(['tunits=mdsvalue(''units_of(dim_of(_rdaeff' user diagname ',' num2str(idim_t) '))'');']);
adata.dimunits = {xunits, tunits};
% $$$ case 3
@@ -387,7 +396,14 @@ else
adata.t=[1:size(adata.data,2)];
end
else
- if length(size(adata.data))<=2; adata.data = adata.data'; end % cannot transpose Nd>2 matrix
+ if strcmp(upper(diagname),'IDA')
+ % rho, time in original dimension, so do not transpose
+ else
+ if length(size(adata.data))<=2;
+ adata.data = adata.data';
+ did_transpose = 1;
+ end % cannot transpose Nd>2 matrix
+ end
if ~isempty(adata.time_aug)
if length(size(adata.data))<=2;
adata.x=[1:prod(size(adata.data))/length(adata_time.value)];
@@ -399,6 +415,9 @@ else
adata.x=[1:size(adata.data,1)];
adata.t=[1:size(adata.data,2)];
end
+ if ~isempty(adata.area)
+ adata.x = adata.area.value{1};
+ end
end
adata.units = adata.unit;
elseif isempty(param_set_name) && ~area_base && ~time_base
@@ -429,9 +448,9 @@ else
% area-base
try
if ~isempty(extra_arg_sf2sig) && ~strcmp(extra_arg_sf2sig,'[]')
- eval(['[adata]=sf2ab(diagname,shot,sigtype,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
+ eval(['[adata]=sf2ab(diagname,shot,area_base_name,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
else
- [adata]=sf2ab(diagname,shot,sigtype,'-exp',shotfile_exp);
+ [adata]=sf2ab(diagname,shot,area_base_name,'-exp',shotfile_exp);
end
adata.data = adata.value{1};
adata.value = adata.data;
@@ -442,9 +461,9 @@ else
% time-base
try
if ~isempty(extra_arg_sf2sig) && ~strcmp(extra_arg_sf2sig,'[]')
- eval(['[adata]=sf2tb(diagname,shot,sigtype,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
+ eval(['[adata]=sf2tb(diagname,shot,time_base_name,''-exp'',shotfile_exp,' extra_arg_sf2sig ');']);
else
- [adata]=sf2tb(diagname,shot,sigtype,'-exp',shotfile_exp);
+ [adata]=sf2tb(diagname,shot,time_base_name,'-exp',shotfile_exp);
end
adata.data = adata.value;
catch ME
@@ -472,4 +491,11 @@ if strcmp(upper(sigtype),'PNIQ')
end
end
+if strcmp(diagname,'IDA') && strcmp(sigtype,'rhot')
+ % In some shots, rhot was multiplied by 1e6, so fix here:
+ if max(max(adata.data)) > 1e5
+ adata.data = adata.data / 1e6;
+ end
+end
+
error=0;
diff --git a/crpptbx/AUG/read_ech_py.m b/crpptbx/AUG/read_ech_py.m
index 56018c79..da6f7305 100644
--- a/crpptbx/AUG/read_ech_py.m
+++ b/crpptbx/AUG/read_ech_py.m
@@ -7,30 +7,29 @@ function [torbeam] = read_ech_py(shot,varargin);
%
% varargin{1}: time array for Torbeam relevant data output
% varargin{2}: min_power_on to decide if the beam was on at a given time (default P>=2e3)
+% varargin{3}: source for TB angles, etc: 0 (default, aug_setval2tp_mex function useing setmirrors library); 1 (TBP/xxx/git shotfile)
%
-torbeam.rfmod = - ones(1,8); % 1 for O-mode, -1 for X-mode
-torbeam.tb_par.('xxb') = [238., 238., 231.1, 231.1, 236.4, 236.4, 236.4 ,236.4];
-torbeam.tb_par.('xzb') = [0., 0., 0., 0., 32., 32., -32., -32.];
-torbeam.tb_par.('xryyb') = [129.4,129.4,129.4, 129.4, 85.4, 85.4, 85.4, 85.4]; % 85.0 in GT's py, 85.4 in tb_gui saved files
-torbeam.tb_par.('xrzzb') = [129.4,129.4,129.4, 129.4, 85.4, 85.4, 85.4, 85.4];
-torbeam.tb_par.('xwyyb') = [2.97, 2.97, 2.97, 2.97, 2.55, 2.55, 2.55, 2.55];
-torbeam.tb_par.('xwzzb') = [2.97, 2.97, 2.97, 2.97, 2.55, 2.55, 2.55, 2.55];
-torbeam.tr2as = { ...
- 'XECECH', 'RR_n'; 'ZECECH', 'ZZ_n'; ...
- 'ECB_WIDTH_HOR', 'xwyyb'; 'ECB_WIDTH_VERT', 'xwzzb'; ...
- 'ECB_CURV_HOR' , 'xryyb'; 'ECB_CURV_VERT' , 'xrzzb' };
+if nargin==0
+ error('read_ech_py requires at least shot number as input');
+else
+ torbeam.shot = shot;
+end
+
+ecsystems=[4,4];
+nbgy = sum(ecsystems);
+igy_seq=[1:ecsystems(1);ecsystems(1)+1:ecsystems(1)+ecsystems(2)]';
-%aa=gdat(shot,'pgyro'); % loads ECS and ECN stuff
+load_option = 0,
+if nargin >= 4 && ~isempty(varargin{3})
+ load_option = varargin{3};
+end
+torbeam.read_ech_options.load_option = load_option;
+% Load parameters independent of angles, like power etc:
+%
tec_min = NaN;
tec_max = NaN;
-
-% following python script:
-
-% WARNING seems alpha, beta inverted in read_ech.py script, so here use alpha for poloidal related angle and beta for tor
-% Thus use alpha_pol and beta_tor to avoid confusion (inverted with respect to .py alpha, beta script)
-
ecsystems=[4,4];
nbgy = sum(ecsystems);
igy_seq=[1:ecsystems(1);ecsystems(1)+1:ecsystems(1)+ecsystems(2)]';
@@ -39,16 +38,6 @@ for isys=1:length(ecsystems)
if isys==1 && shot > 33725; isyseff = 3; end
for igy=1:ecsystems(isys)
ps_source = ['P_sy' num2str(isyseff) '_g' num2str(igy)];
- aa=gdat(shot,{'ECS',ps_source,'','param:GPolPos'},'machine','aug');
- torbeam.alpha_pol(igy_seq(igy,isys),1) = aa.data;
- aa=gdat(shot,{'ECS',ps_source,'','param:GTorPos'},'machine','aug');
- torbeam.beta_tor(igy_seq(igy,isys),1) = aa.data;
- aa=gdat(shot,{'ECS',ps_source,'','param:gyr_freq'},'machine','aug');
- if aa.data > 1000
- torbeam.freq(igy_seq(igy,isys),1) = aa.data;
- else
- torbeam.freq(igy_seq(igy,isys),1) = 1.398765E+11;
- end
if isyseff == 1
aaon=gdat(shot,{'ECS',ps_source,'','param:gy_on'},'machine','aug');
else
@@ -76,29 +65,6 @@ for isys=1:length(ecsystems)
end
end
-% Compute theta,phi for TORBEAM
-torbeam.theta = NaN(nbgy,1);
-torbeam.phi = NaN(nbgy,1);
-torbeam.err_status = -ones(nbgy,1);
-torbeam.isystunit = NaN(nbgy,1);
-torbeam.datum = NaN(nbgy,1);
-for isys=1:length(ecsystems)
- for igy=1:ecsystems(isys)
- igy_1_8 = igy_seq(igy,isys);
- if torbeam.is_gyro_on(igy_seq(igy,isys))
- [torbeam.theta(igy_seq(igy,isys),1),torbeam.phi(igy_seq(igy,isys),1),torbeam.err_status(igy_seq(igy,isys),1),torbeam.isystunit(igy_seq(igy,isys),1),torbeam.datum(igy_seq(igy,isys),1)] = ...
- aug_setval2tp_mex(torbeam.alpha_pol(igy_seq(igy,isys),1),torbeam.beta_tor(igy_seq(igy,isys),1),igy_1_8,shot);
- end
- end
-end
-
-% get time dependent data
-
-if shot < 23187
- % Before ECRH2 installed
- return
-end
-
if ~isfinite(tec_min) || ~isfinite(tec_max)
warning('no good time interval')
tec_min
@@ -114,9 +80,6 @@ if tec_max <= tec_min
return
end
-isys=length(ecsystems);
-isyseff = isys;
-diagname = 'ECN';
if nargin>=2 && ~isempty(varargin{1})
time_for_torbeam = varargin{1};
else
@@ -126,31 +89,19 @@ torbeam.t = time_for_torbeam;
torbeam.theta_t = NaN(nbgy,length(time_for_torbeam));
torbeam.phi_t = NaN(nbgy,length(time_for_torbeam));
torbeam.err_status_t = NaN(nbgy,length(time_for_torbeam));
-% 1st system angles fixed so use above values
+
isys = 1;
for igy=1:ecsystems(isys)
igy_1_8 = igy_seq(igy,isys);
if torbeam.is_gyro_on(igy_seq(igy,isys))
- torbeam.theta_t(igy_1_8,:) = torbeam.theta(igy_seq(igy,isys),1);
- torbeam.phi_t(igy_1_8,:) = torbeam.phi(igy_seq(igy,isys),1);
- torbeam.err_status_t(igy_1_8,:) = torbeam.err_status(igy_seq(igy,isys),1);
- % linear interpolation for power to avoid over/under shoot
+ % linear interpolation for power to avoid over/under shoot
torbeam.power_t(igy_1_8,:) = interpos(-21,torbeam.pow{igy_1_8}.t,torbeam.pow{igy_1_8}.data,torbeam.t);
end
end
isys = 2;
for igy=1:ecsystems(isys)
- torbeam.alpha_polN{igy} = gdat(shot,{diagname,['G' num2str(igy) 'POL']},'machine','aug');
- torbeam.alpha_polN{igy}.data = 0.01 .* torbeam.alpha_polN{igy}.data;
- torbeam.alpha_pol_t = interpos(torbeam.alpha_polN{igy}.t,torbeam.alpha_polN{igy}.data,torbeam.t,-1e2);
-% $$$ torbeam.beta_torN{igy} = gdat(shot,{diagname,['G' num2str(igy) 'TOR']},'machine','aug');
igy_1_8 = igy_seq(igy,isys);
if torbeam.is_gyro_on(igy_seq(igy,isys))
- [theta_t,phi_t,err_status_t,torbeam.isystunit(igy_seq(igy,isys),1),torbeam.datum(igy_seq(igy,isys),1)] = ...
- aug_setval2tp_mex(torbeam.alpha_pol_t,torbeam.beta_tor(igy_seq(igy,isys),1).*ones(size(torbeam.alpha_pol_t)),igy_1_8,shot);
- torbeam.theta_t(igy_1_8,:) = theta_t';
- torbeam.phi_t(igy_1_8,:) = phi_t';
- torbeam.err_status_t(igy_1_8,:) = err_status_t';
torbeam.power_t(igy_1_8,:) = interpos(-21,torbeam.pow{igy_1_8}.t,torbeam.pow{igy_1_8}.data,torbeam.t);
end
end
@@ -160,6 +111,151 @@ if nargin>=3 && ~isempty(varargin{2})
else
min_power_on = 2e3;
end
+torbeam.read_ech_options.min_power_on = min_power_on;
+
torbeam.beam_on = false.*ones(size(torbeam.power_t));
torbeam.beam_on(torbeam.power_t >= min_power_on) = true;
torbeam.t_some_beam_on = torbeam.t(any(torbeam.beam_on,1));
+
+switch load_option
+ case 0
+ torbeam.rfmod = - ones(1,8); % 1 for O-mode, -1 for X-mode
+ torbeam.tb_par.('xxb') = [238., 238., 231.1, 231.1, 236.4, 236.4, 236.4 ,236.4];
+ torbeam.tb_par.('xzb') = [0., 0., 0., 0., 32., 32., -32., -32.];
+ torbeam.tb_par.('xryyb') = [129.4,129.4,129.4, 129.4, 85.4, 85.4, 85.4, 85.4]; % 85.0 in GT's py, 85.4 in tb_gui saved files
+ torbeam.tb_par.('xrzzb') = [129.4,129.4,129.4, 129.4, 85.4, 85.4, 85.4, 85.4];
+ torbeam.tb_par.('xwyyb') = [2.97, 2.97, 2.97, 2.97, 2.55, 2.55, 2.55, 2.55];
+ torbeam.tb_par.('xwzzb') = [2.97, 2.97, 2.97, 2.97, 2.55, 2.55, 2.55, 2.55];
+ % from Ringberg 2018 Martin Schubert
+ torbeam.tb_par.('xryyb') = [129.4,129.4,129.4, 129.4, 122.7, 122.7, 139.7, 139.7]; % 85.0 in GT's py, 85.4 in tb_gui saved files
+ torbeam.tb_par.('xrzzb') = [129.4,129.4,129.4, 129.4, 122.7, 122.7, 139.7, 139.7];
+ torbeam.tb_par.('xwyyb') = [2.97, 2.97, 2.97, 2.97, 3.589, 3.589, 3.196, 3.196];
+ torbeam.tb_par.('xwzzb') = [2.97, 2.97, 2.97, 2.97, 3.589, 3.589, 3.196, 3.196];
+
+% $$$ torbeam.tr2as = { ...
+% $$$ 'XECECH', 'RR_n'; 'ZECECH', 'ZZ_n'; ...
+% $$$ 'ECB_WIDTH_HOR', 'xwyyb'; 'ECB_WIDTH_VERT', 'xwzzb'; ...
+% $$$ 'ECB_CURV_HOR' , 'xryyb'; 'ECB_CURV_VERT' , 'xrzzb' };
+
+ % following python script:
+ % WARNING seems alpha, beta inverted in read_ech.py script, so here use alpha for poloidal related angle and beta for tor
+ % Thus use alpha_pol and beta_tor to avoid confusion (inverted with respect to .py alpha, beta script)
+ for isys=1:length(ecsystems)
+ isyseff = isys;
+ if isys==1 && shot > 33725; isyseff = 3; end
+ for igy=1:ecsystems(isys)
+ ps_source = ['P_sy' num2str(isyseff) '_g' num2str(igy)];
+ aa=gdat(shot,{'ECS',ps_source,'','param:GPolPos'},'machine','aug');
+ torbeam.alpha_pol(igy_seq(igy,isys),1) = aa.data;
+ aa=gdat(shot,{'ECS',ps_source,'','param:GTorPos'},'machine','aug');
+ torbeam.beta_tor(igy_seq(igy,isys),1) = aa.data;
+ if load_option == 0
+ aa=gdat(shot,{'ECS',ps_source,'','param:gyr_freq'},'machine','aug');
+ if aa.data > 1000
+ torbeam.freq(igy_seq(igy,isys),1) = aa.data;
+ else
+ torbeam.freq(igy_seq(igy,isys),1) = 1.398765E+11;
+ end
+ else
+ % assume already loaded in above switch case
+ end
+ end
+ end
+
+ % Compute theta,phi for TORBEAM
+ torbeam.theta = NaN(nbgy,1);
+ torbeam.phi = NaN(nbgy,1);
+ torbeam.err_status = -ones(nbgy,1);
+ torbeam.isystunit = NaN(nbgy,1);
+ torbeam.datum = NaN(nbgy,1);
+ for isys=1:length(ecsystems)
+ for igy=1:ecsystems(isys)
+ igy_1_8 = igy_seq(igy,isys);
+ if torbeam.is_gyro_on(igy_seq(igy,isys))
+ [torbeam.theta(igy_seq(igy,isys),1),torbeam.phi(igy_seq(igy,isys),1),torbeam.err_status(igy_seq(igy,isys),1),torbeam.isystunit(igy_seq(igy,isys),1),torbeam.datum(igy_seq(igy,isys),1)] = ...
+ aug_setval2tp_mex(torbeam.alpha_pol(igy_seq(igy,isys),1),torbeam.beta_tor(igy_seq(igy,isys),1),igy_1_8,shot);
+ end
+ end
+ end
+
+ % get time dependent data
+ if shot < 23187
+ % Before ECRH2 installed
+ return
+ end
+
+ isys=length(ecsystems);
+ isyseff = isys;
+ diagname = 'ECN';
+ % 1st system angles fixed so use above values
+ isys = 1;
+ for igy=1:ecsystems(isys)
+ igy_1_8 = igy_seq(igy,isys);
+ if torbeam.is_gyro_on(igy_seq(igy,isys))
+ torbeam.theta_t(igy_1_8,:) = torbeam.theta(igy_seq(igy,isys));
+ torbeam.phi_t(igy_1_8,:) = torbeam.phi(igy_seq(igy,isys));
+ torbeam.err_status_t(igy_1_8,:) = torbeam.err_status(igy_seq(igy,isys),1);
+ end
+ end
+ isys = 2;
+ for igy=1:ecsystems(isys)
+ torbeam.alpha_polN{igy} = gdat(shot,{diagname,['G' num2str(igy) 'POL']},'machine','aug');
+ torbeam.alpha_polN{igy}.data = 0.01 .* torbeam.alpha_polN{igy}.data;
+ torbeam.alpha_pol_t = interpos(-21,torbeam.alpha_polN{igy}.t,torbeam.alpha_polN{igy}.data,torbeam.t,-1e2);
+% $$$ torbeam.beta_torN{igy} = gdat(shot,{diagname,['G' num2str(igy) 'TOR']},'machine','aug');
+ igy_1_8 = igy_seq(igy,isys);
+ if torbeam.is_gyro_on(igy_seq(igy,isys))
+ [theta_t,phi_t,err_status_t,torbeam.isystunit(igy_seq(igy,isys),1),torbeam.datum(igy_seq(igy,isys),1)] = ...
+ aug_setval2tp_mex(torbeam.alpha_pol_t,torbeam.beta_tor(igy_seq(igy,isys),1).*ones(size(torbeam.alpha_pol_t)),igy_1_8,shot);
+ torbeam.theta_t(igy_1_8,:) = theta_t';
+ torbeam.phi_t(igy_1_8,:) = phi_t';
+ torbeam.err_status_t(igy_1_8,:) = err_status_t';
+ end
+ end
+
+ case 1
+ torbeam.rfmod = - ones(1,8); % 1 for O-mode, -1 for X-mode
+ tmp=gdat(shot,{'TBP','TB_par','git','param:Rlaunch'},'machine','aug');
+ torbeam.tb_par.('xxb') = reshape(tmp.data,1,length(tmp.data)) * 100.;
+ tmp=gdat(shot,{'TBP','TB_par','git','param:Zlaunch'},'machine','aug');
+ torbeam.tb_par.('xzb') = reshape(tmp.data,1,length(tmp.data)) * 100.;
+ tmp=gdat(shot,{'TBP','TB_par','git','param:wid_hor'},'machine','aug');
+ torbeam.tb_par.('xwyyb') = reshape(tmp.data,1,length(tmp.data)) * 100.;
+ tmp=gdat(shot,{'TBP','TB_par','git','param:wid_ver'},'machine','aug');
+ torbeam.tb_par.('xwzzb') = reshape(tmp.data,1,length(tmp.data)) * 100.;
+ tmp=gdat(shot,{'TBP','TB_par','git','param:curv_hor'},'machine','aug');
+ torbeam.tb_par.('xryyb') = reshape(tmp.data,1,length(tmp.data)) * 100.;
+ tmp=gdat(shot,{'TBP','TB_par','git','param:curv_ver'},'machine','aug');
+ torbeam.tb_par.('xrzzb') = reshape(tmp.data,1,length(tmp.data)) * 100.;
+ tmp=gdat(shot,{'TBP','TB_par','git','param:the_ps'},'machine','aug');
+ torbeam.tb_par_sf.('the_ps') = reshape(tmp.data,1,length(tmp.data));
+ tmp=gdat(shot,{'TBP','TB_par','git','param:phi_ps'},'machine','aug');
+ torbeam.tb_par_sf.('phi_ps') = reshape(tmp.data,1,length(tmp.data));
+ tmp=gdat(shot,{'TBP','TB_par','git','param:freq'},'machine','aug');
+ torbeam.tb_par_sf.('freq') = reshape(tmp.data,1,length(tmp.data));
+ tmp=gdat(shot,{'TBP','the_ec2','git'},'machine','aug');
+ torbeam.tbp.the_ec2 = tmp;
+ tmp=gdat(shot,{'TBP','phi_ec2','git'},'machine','aug');
+ torbeam.tbp.phi_ec2 = tmp;
+
+ % fill in relevant output tables
+ torbeam.freq = torbeam.tb_par_sf.freq';
+ torbeam.theta = torbeam.tb_par_sf.the_ps';
+ torbeam.phi = torbeam.tb_par_sf.phi_ps';
+ isys = 1;
+ igy_1_8 = igy_seq(1:ecsystems(isys),isys); % sys1 has fixed angles vs time
+ torbeam.theta_t(igy_1_8,:) = torbeam.theta(igy_1_8)*ones(1,size(torbeam.theta_t,2));
+ torbeam.phi_t(igy_1_8,:) = torbeam.phi(igy_1_8)*ones(1,size(torbeam.phi_t,2));
+ isys = 2;
+ igy_1_8 = igy_seq(1:ecsystems(isys),isys); % sys1 has fixed angles vs time
+ for igy=1:size(igy_seq,1)
+ igy_1_8 = igy_seq(igy,isys);
+ torbeam.theta_t(igy_1_8,:) = interpos(-21,torbeam.tbp.the_ec2.t,torbeam.tbp.the_ec2.data(igy,:),torbeam.t,-1e2);
+ torbeam.phi_t(igy_1_8,:) = interpos(-21,torbeam.tbp.phi_ec2.t,torbeam.tbp.phi_ec2.data(igy,:),torbeam.t,-1e2);
+ end
+
+ otherwise
+ error(['option for Torbeam parameters and angles = ' num2str(load_option) ' not known nor implemented yet. Ask O. Sauter']);
+end
+
+torbeam = orderfields(torbeam);
diff --git a/crpptbx/AUG/torbeam_plot_output_singlerun.m b/crpptbx/AUG/torbeam_plot_output_singlerun.m
index 654e9d0e..a7a34ec8 100644
--- a/crpptbx/AUG/torbeam_plot_output_singlerun.m
+++ b/crpptbx/AUG/torbeam_plot_output_singlerun.m
@@ -122,12 +122,15 @@ xlabel('rhopol')
figure;
subplot(2,1,1)
-plot(zzz.t2_new_LIB(:,1),zzz.t2_new_LIB(:,2))
+try
+ plot(zzz.t2_new_LIB(:,1),zzz.t2_new_LIB(:,2))
ylabel('Pdens')
subplot(2,1,2)
plot(zzz.t2_new_LIB(:,1),zzz.t2_new_LIB(:,3))
ylabel('jcd')
xlabel('rhopol')
+catch
+end
figure
theta=linspace(0,2.*pi,300);
diff --git a/crpptbx/AUG/torbeam_prepare_inputs_and_run.m b/crpptbx/AUG/torbeam_prepare_inputs_and_run.m
index 20be22b0..4f081cf6 100644
--- a/crpptbx/AUG/torbeam_prepare_inputs_and_run.m
+++ b/crpptbx/AUG/torbeam_prepare_inputs_and_run.m
@@ -1,20 +1,25 @@
function [torbeam_out,torbeam_in,varargout] = torbeam_prepare_inputs_and_run(shot,times_in,varargin);
%
-% [torbeam_out,torbeam_in,varargout] = aaa(shot,times_in,varargin);
+% [torbeam_out,torbeam_in,varargout] = torbeam_prepare_inputs_and_run(shot,times_in,varargin);
%
% varargin{1}: torbeam_in (from a previous run, so avoids to reload all the data)
% leave then times_in empty to rerun at all times in torbeam_in.inbeam.t_some_beam_on or select times points from that array
+% varargin{2}: source for gdat profiles/equilibrium related: [] (default sources), 'IDA' (IDA/IDE sources), {'TRA','expname'} for TRANSP
+% varargin{3}: source for torbeam angles: 0 (default) from conversion from ECS/ECN with aug_setval2tp_mex; 1 from TBP/git shotfile
+% varargin{4}: nb of radial points for outputs in shotfile structure (default 401)
%
% Examples:
% [torbeam_out,torbeam_in] = torbeam_prepare_inputs_and_run(30594,linpace(0,10,101));
% [torbeam_out2] = torbeam_prepare_inputs_and_run(30594,torbeam_in.inbeam.t_some_beam_on(1:30:end),torbeam_in); % to rerun a few cases
%
+torbeam_in = [];
+torbeam_out = [];
+
if nargin==0 || isempty(shot)
error('requires shot number')
return
end
-
times_for_torbeam_eff = [];
if nargin>=2 && ~isempty(times_in)
times_for_torbeam_eff = times_in;
@@ -25,10 +30,44 @@ if nargin>=3 && ~isempty(varargin{1})
torbeam_in = varargin{1};
doload = 0;
end
+torbeam_in.prepare_run_options.doload = doload;
+
+source_gdat_profiles = [];
+source_gdat_profiles_exp_name = [];
+if nargin>=4 && ~isempty(varargin{2})
+ if ~iscell( varargin{2})
+ source_gdat_profiles = varargin{2};
+ else
+ if length(varargin{2}) >= 2
+ source_gdat_profiles = varargin{2}{1};
+ source_gdat_profiles_exp_name = varargin{2}{2};
+ else
+ error('expects cell(2) if for source and exp_name if a cell is provided for profiles');
+ end
+ end
+end
+torbeam_in.prepare_run_options.source_gdat_profiles = source_gdat_profiles;
+torbeam_in.prepare_run_options.source_gdat_profiles_exp_name = source_gdat_profiles_exp_name;
+
+tbangles_source = 0;
+if nargin>=5 && ~isempty(varargin{3})
+ tbangles_source = varargin{3};
+end
+torbeam_in.prepare_run_options.tbangles_source = tbangles_source;
+
+nb_points_shotfile = 401;
+if nargin>=6 && ~isempty(varargin{4})
+ nb_points_shotfile = varargin{4};
+end
+torbeam_in.prepare_run_options.nb_points_shotfile = nb_points_shotfile;
if doload
tic
- torbeam_in.nete=gdat(shot,'nete_rho','machine','aug');
+ if isempty(source_gdat_profiles_exp_name)
+ torbeam_in.nete=gdat(shot,'nete_rho','machine','aug','source',source_gdat_profiles);
+ else
+ torbeam_in.nete=gdat(shot,'nete_rho','machine','aug','source',source_gdat_profiles,'source_exp_name',source_gdat_profiles_exp_name);
+ end
% qrho_and_grid=gdat(shot,'q_rho','machine','aug');
else
if ~isfield(torbeam_in,'nete') || isempty(torbeam_in.nete)
@@ -42,7 +81,7 @@ if doload && (~exist('times_for_torbeam_eff') || isempty(times_for_torbeam_eff))
end
% load beam releated data including theta,phi in torbeam reference
if doload,
- [torbeam_in.inbeam] = read_ech_py(shot,times_for_torbeam_eff);
+ [torbeam_in.inbeam] = read_ech_py(shot,times_for_torbeam_eff,[],tbangles_source);
times_for_torbeam_eff = torbeam_in.inbeam.t_some_beam_on;
else
if ~isfield(torbeam_in,'inbeam') || isempty(torbeam_in.inbeam)
@@ -60,7 +99,7 @@ if doload==0
end
end
if doload,
- torbeam_in.eqd_all=gdat(shot,'eqdsk','time',torbeam_in.inbeam.t_some_beam_on,'write',0);
+ torbeam_in.eqd_all=gdat(shot,'eqdsk','time',torbeam_in.inbeam.t_some_beam_on,'write',0,'source',source_gdat_profiles);
torbeam_in.dir_torbeam = ['/tmp/' getenv('USER') '/' num2str(shot) '_torbeam'];
unix(['mkdir -p ' torbeam_in.dir_torbeam]);
else
@@ -81,6 +120,10 @@ if doload,
itime_ne = iround_os(torbeam_in.nete.fit.t,times_for_torbeam_eff);
it_inb = iround_os(torbeam_in.inbeam.t,times_for_torbeam_eff);
cocos_out = 17;
+ cocos_in = 17;
+ if strcmp(upper(source_gdat_profiles),'TRA')
+ cocos_in = 17;
+ end
tic
for it=1:length(torbeam_in.inbeam.t_some_beam_on)
% Te
@@ -103,8 +146,8 @@ if doload,
fprintf(fid,'%.6f %.6f\n',aa');
fclose(fid);
% topfile
- % eqd_with_b = read_eqdsk(torbeam_in.eqd_all.eqdsk{it},17,1);
- eqd_with_b2=eqdsk_cocos_transform(torbeam_in.eqd_all.eqdsk{it},[17 cocos_out]);
+ % eqd_with_b = read_eqdsk(torbeam_in.eqd_all.eqdsk{it},cocos_in,1);
+ eqd_with_b2=eqdsk_cocos_transform(torbeam_in.eqd_all.eqdsk{it},[cocos_in cocos_out]);
eqd_with_b1=eqdsk_transform(eqd_with_b2,65,129,-0.1);
eqd_with_b = read_eqdsk(eqd_with_b1,cocos_out,1); % to compute BR, etc
torbeam_in.fname{it}.topfile = fullfile(torbeam_in.dir_torbeam,['topfile' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it))]);
@@ -171,6 +214,10 @@ torbeam_out.run_torbeam = false .* ones(size(torbeam_in.inbeam.power_t,1),length
it_inb_beam_on = iround_os(torbeam_in.inbeam.t_some_beam_on,times_for_torbeam_eff);
it_inb = iround_os(torbeam_in.inbeam.t,times_for_torbeam_eff);
torbeam_out.t = torbeam_in.inbeam.t_some_beam_on(it_inb_beam_on);
+torbeam_out.shotfile.rhopol = linspace(0,1,nb_points_shotfile);
+torbeam_out.shot = shot;
+torbeam_out.shotfile.shot = shot;
+torbeam_out.shotfile.t = torbeam_out.t;
tic
icount = 0;
for it_eff=1:length(times_for_torbeam_eff)
@@ -196,11 +243,31 @@ for it_eff=1:length(times_for_torbeam_eff)
torbeam_out.(tb_out_files_to_cp_varname{j}){ibeam,it_eff} = load(fullfile(torbeam_in.dir_torbeam,[tb_out_files_to_cp{j}]));
end
nnn = load('cntpr.dat');
- torbeam_out.cntpr.tbegabs(ibeam,it_eff) = nnn(1);
- torbeam_out.cntpr.dtabs(ibeam,it_eff) = nnn(2);
+ torbeam_out.cntpr.tbegabs(ibeam,it_eff) = nnn(2);
+ torbeam_out.cntpr.dtabs(ibeam,it_eff) = nnn(1);
torbeam_out.pdens_jcd_rho(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,1);
- torbeam_out.pdens(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,2);
- torbeam_out.jcd(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,3);
+ torbeam_out.pdens(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,2) * 1e6;
+ torbeam_out.jcd(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,3) * 1e6;
+ % structure for shot file, add rho=0 point, use a generic nb_points_shotfile
+ x2 = [0.; torbeam_out.pdens_jcd_rho(:,ibeam,it_eff)];
+ pdens2 = [0.; torbeam_out.pdens(:,ibeam,it_eff)];
+ jcd2 = [0.; torbeam_out.jcd(:,ibeam,it_eff)];
+ torbeam_out.shotfile.volume(:,it_eff) = ...
+ interpos([0.;torbeam_out.volumes{ibeam,it_eff}(3:end,1)],[0.;torbeam_out.volumes{ibeam,it_eff}(3:end,2)],torbeam_out.shotfile.rhopol);
+ torbeam_out.shotfile.volume(1,it_eff) = 0.;
+ torbeam_out.shotfile.pdens(:,ibeam,it_eff) = interpos(x2,pdens2,torbeam_out.shotfile.rhopol);
+ torbeam_out.shotfile.jcd(:,ibeam,it_eff) = interpos(x2,jcd2,torbeam_out.shotfile.rhopol);
+ torbeam_out.shotfile.rz_r_ray(1:size(torbeam_out.t1_LIB{ibeam,it_eff},1),1:size(torbeam_out.t1_LIB{ibeam,it_eff},2)/2,ibeam,it_eff) = torbeam_out.t1_LIB{ibeam,it_eff}(:,1:2:end-1);
+ torbeam_out.shotfile.rz_z_ray(1:size(torbeam_out.t1_LIB{ibeam,it_eff},1),1:size(torbeam_out.t1_LIB{ibeam,it_eff},2)/2,ibeam,it_eff) = torbeam_out.t1_LIB{ibeam,it_eff}(:,2:2:end);
+ torbeam_out.shotfile.ry_r_ray(1:size(torbeam_out.t1_LIB{ibeam,it_eff},1),1:size(torbeam_out.t1_LIB{ibeam,it_eff},2)/2,ibeam,it_eff) = torbeam_out.t1tor_LIB{ibeam,it_eff}(:,1:2:end-1);
+ torbeam_out.shotfile.ry_y_ray(1:size(torbeam_out.t1_LIB{ibeam,it_eff},1),1:size(torbeam_out.t1_LIB{ibeam,it_eff},2)/2,ibeam,it_eff) = torbeam_out.t1tor_LIB{ibeam,it_eff}(:,2:2:end);
+ torbeam_out.shotfile.index_abs_ray(1,ibeam,it_eff) = nnn(2) + 1;
+ torbeam_out.shotfile.index_abs_ray(2,ibeam,it_eff) = nnn(2) + 1 + nnn(1);
+ torbeam_out.shotfile.torbeam_was_run = torbeam_out.run_torbeam;
+ torbeam_out.shotfile.pabs = torbeam_out.pabs;
+ torbeam_out.shotfile.icd = torbeam_out.icd;
+ torbeam_out.shotfile.pabs_rhopolmax = torbeam_out.pabs_rhopolmax;
+ torbeam_out.shotfile.icd_rhopolmax = torbeam_out.icd_rhopolmax;
else
disp(['TORBEAM did not run properly, no MW result'])
end
@@ -223,17 +290,17 @@ end
figure
subplot(3,1,1)
-plot(torbeam_out.t,torbeam_out.pabs/1e6,'*-')
+plot(torbeam_out.t,torbeam_out.pabs/1e6,'-')
ylabel('Pabs [MW]')
title(['AUG shot #' num2str(shot)])
subplot(3,1,2)
-plot(torbeam_out.t,torbeam_out.icd/1e6,'*-')
+plot(torbeam_out.t,torbeam_out.icd/1e6,'-')
ylabel('Icd [MA]')
subplot(3,1,3)
-plot(torbeam_out.t,torbeam_out.pabs_rhopolmax,'*-')
+plot(torbeam_out.t,torbeam_out.pabs_rhopolmax,'-')
ylabel('rhopol(max Pdens)')
xlabel('t [s]')
--
GitLab