Merge branch 'IMAS_codecamp_develop' into 'master'

Imas codecamp develop

See merge request spc/tcv/tbx/gdat!16

matlab/TCV/gdat_tcv.m

@@ -1991,11 +1991,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
         end
       end
       if NBH_in_TCV
-        nodenameeff = '\ATLAS::NBH.DATA.MAIN_ADC:DATA';
+        nodenameeff = '\results::NBH:POWR_TCV';
         nbh_data_tdi = tdi(nodenameeff);
         if ~isempty(nbh_data_tdi.data) && ~ischar(nbh_data_tdi.data) && ~isempty(nbh_data_tdi.dim)
-          index_for_tot_pow = 37;
-          nbi_neutral_power_tot = nbh_data_tdi.data(:,index_for_tot_pow).*1e6; % in W
+          nbi_neutral_power_tot = nbh_data_tdi.data.*1e6; % in W
           nbi_neutral_power_tot = max(nbi_neutral_power_tot,0.);
           gdat_data.nbi.data = nbi_neutral_power_tot; % at this stage p_gyro is in kW'
           gdat_data.nbi.units = 'W';
@@ -2003,11 +2002,13 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
           gdat_data.nbi.dimunits{1}=nbh_data_tdi.dimunits{1};
           gdat_data.nbi.t=gdat_data.nbi.dim{1};
           gdat_data.nbi.x=[];
-          gdat_data.nbi.data_fullpath=[nodenameeff ', index ' num2str(index_for_tot_pow) ', CLC NEUTRAL POWER'];
+          gdat_data.nbi.data_fullpath=[nodenameeff];
           gdat_data.nbi.label='P_{nbi}';
           gdat_data.nbi.help = nbh_data_tdi.help;
-          index_for_energy = 33;
-          gdat_data.nbi.energy = nbh_data_tdi.data(:,index_for_energy);
+          nodenameeff2 = '\results::nbh:energy';
+          gdat_data.nbi.data_fullpath=[nodenameeff ' *1e6 for W; and ' nodenameeff2 ' *1e3 for eV'];
+          nbh_energy_data_tdi = tdi(nodenameeff2);
+          gdat_data.nbi.energy = nbh_energy_data_tdi.data*1e3; % for eV
           % add to main with linear interpolation and 0 for extrapolated values
           ij=[isfinite(gdat_data.nbi.data)];
           gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t(ij),gdat_data.nbi.data(ij),gdat_data.t);

matlab/TCV_IMAS/ids2database.m

@@ -76,7 +76,7 @@ ids_names=fieldnames(ids2put);
 try
   ids_full_list = IDS_list;
 catch
-  ids_full_list = {'equilibrium', 'magnetics', 'tf', 'pf_active','wall','core_profiles'};
+  ids_full_list = {'equilibrium', 'magnetics', 'tf', 'pf_active','wall','core_profiles','ec_antennas','nbi'};
   warning(['IDS_list not available, quick fix introducing list of ids ready for TCV: ' fprintf('%s ',ids_full_list{:}) char(10)]);
 end
 ids_names_notok = setdiff(ids_names,ids_full_list);
@@ -100,9 +100,9 @@ try
   end
   if shot_is_new
     idx  = imas_create_env('ids',shot,run,0,0,params_ids2database.tree_user,params_ids2database.tree_tokamak, ...
-          params_ids2database.tree_majorversion); % 
+          params_ids2database.tree_majorversion); %
   else
-    idx  = imas_open_env('ids',shot,run,params_ids2database.tree_user,params_ids2database.tree_tokamak,params_ids2database.tree_majorversion); % 
+    idx  = imas_open_env('ids',shot,run,params_ids2database.tree_user,params_ids2database.tree_tokamak,params_ids2database.tree_majorversion); %
   end
 
   %% Put the field
@@ -114,7 +114,7 @@ try
       ids_put(idx,[ids_to_get '/' num2str(occurence)],ids2put.(ids_to_get));
     end
   end
-  
+
   %% Close the file
   imas_close(idx)
 catch ME

matlab/TCV_IMAS/tcv2ids.m

@@ -11,7 +11,7 @@ function [ids_from_tcv,varargout] = tcv2ids(shot,varargin);
 % mdsconnect('localhost:5555')
 %
 % varargin: 'ids_names': cell, idss to load, by default all defined so far (if empty or empty string or not given)
-%             {'equilibrium', 'magnetics', 'pf_active','wall','core_profiles'} or a subset
+%             {'equilibrium', 'magnetics', 'pf_active','wall','core_profiles','ec_antennas','nbi'} or a subset
 % varargin{2}: 'error_bar': type (string)
 %                empty or 'delta' (default): meaning difference in upper is set (standard error_bar
 %                'added': errorbar is added: upper=data+delta and lower=data-delta
@@ -26,7 +26,7 @@ p = inputParser;
 % no required inputs here so one can call with empty args and get defaults parameters
 % effectively required inputs should have defaults as empty
 p.addOptional('shot', [], @(x) (isnumeric(x) && isscalar(x) && (x == round(x)))); % integer
-p.addOptional('ids_names', {'equilibrium', 'magnetics', 'pf_active','wall', 'tf','core_profiles'}, @(x) isempty(x) | iscell(x) ); % char or cell array
+p.addOptional('ids_names', {'equilibrium', 'magnetics', 'pf_active','wall', 'tf','core_profiles','ec_antennas','nbi'}, @(x) isempty(x) | iscell(x) ); % char or cell array
 p.addOptional('error_bar', 'delta', @(x) isempty(x) | ischar(x) ); % char or cell array
 p.parse;
 defaults_tcv2ids = p.Results; % to keep track of defaults

matlab/TCV_IMAS/tcv2ids2database.m

 function [ids_from_tcv,varargout] = tcv2ids2database(shot,run_out,varargin);
-% 
+%
 % [ids_from_tcv,varargout] = tcv2ids2database(shot,run_out,varargin)
-% 
+%
 % script to get TCV data as ids and then write them on the database
 %
 % Uses function tcv2ids and ids2database
@@ -9,7 +9,7 @@ function [ids_from_tcv,varargout] = tcv2ids2database(shot,run_out,varargin);
 %
 % varargin: 'occurence': occurence value
 % varargin: 'ids_names': cell, ids to load, by default all defined so far (if empty or empty string or not given)
-%             {'equilibrium', 'magnetics', 'tf', 'pf_active','wall','core_profiles'} or a subset
+%             {'equilibrium', 'magnetics', 'tf', 'pf_active','wall','core_profiles','ec_antennas','nbi'} or a subset
 % varargin: 'error_bar': type (string)
 %                'delta' or empty (default): meaning difference in upper is set (standard error_bar
 %                'added': errorbar is added: upper=data+delta and lower=data-delta
@@ -42,7 +42,7 @@ p = inputParser;
 p.addOptional('shot', [], @(x) (isnumeric(x) && isscalar(x) && (x == round(x)))); % integer
 p.addOptional('run_out', [], @(x) (isnumeric(x) && isscalar(x) && (x == round(x)))); % integer
 p.addOptional('occurence', 0, @(x) (isnumeric(x) && isscalar(x) && (x == round(x)))); % integer
-p.addOptional('ids_names', {'equilibrium', 'magnetics', 'pf_active','wall', 'tf','core_profiles'}, @(x) isempty(x) | iscell(x)); % char or cell array
+p.addOptional('ids_names', {'equilibrium', 'magnetics', 'pf_active','wall', 'tf','core_profiles','ec_antennas','nbi'}, @(x) isempty(x) | iscell(x)); % char or cell array
 p.addOptional('error_bar', 'delta', @(x) isempty(x) | ischar(x) ); % char or cell array
 p.addOptional('tree_user', getenv('USER'), @(x) isempty(x) | ischar(x) ); % char
 p.addOptional('tree_tokamak', 'tcv', @(x) isempty(x) | ischar(x) ); % char
@@ -115,3 +115,5 @@ ids_from_tcv.ids_put_status = ids_put_status;
 % tic;[TTTp,idsp]=evalc('ids_get(expIdx, ''pf_active'')');toc
 % tic;[TTTt,idst]=evalc('ids_get(expIdx, ''tf'')');toc
 % tic;[TTTw,idsw]=evalc('ids_get(expIdx, ''wall'')');toc
+% tic;[TTTw,idsw]=evalc('ids_get(expIdx, ''ec_antennas'')');toc
+% tic;[TTTw,idsw]=evalc('ids_get(expIdx, ''nbi'')');toc

matlab/TCV_IMAS/tcv_get_ids_ec_antennas.m

+function [ids_ec_antennas,ids_ec_antennas_description,varargout] = tcv_get_ids_ec_antennas(shot,ids_ec_antennas_empty, gdat_params,varargin);
+%
+%  [ids_ec_antennas,ids_ec_antennas_description,varargout] = tcv_get_ids_ec_antennas(shot,ids_ec_antennas_empty,varargin);
+%
+%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar options
+%
+
+[ids_ec_antennas, params_ec_antennas] = tcv_ids_headpart(shot,ids_ec_antennas_empty,'ec_antennas','homogeneous_time',0,varargin{:});
+
+% As a general rule, for a new substructure under the main ids, construct a local structure like:
+% "global_quantities" with subfields being the relevant data to get and a local structure:
+% "global_quantities_desc" which contains the same subfields themselves containing the gdat string aftre shot used
+%
+
+pow=gdat_tcv(shot,'powers');
+
+if shot>=63000
+  warning('not sure if set-up correct for new launchers, beams, gyrotrons. Check with O. Sauter')
+  return
+end
+
+nb_antennas = 9; % at this stage keep same numbering as standard TCV launchers before mid-2019
+
+ids_ec_antennas.antenna(1:nb_antennas) = ids_ec_antennas.antenna(1); % copy empty structure for all units, then fill in
+% $$$ ids_ec_antennas.time = pow.ec.t; if homogeneous time
+% $$$ ids_ec_antennas_description.time = 'time from power';
+[r_phi0,z_phi0,r1,z1,rc,zc,thetphi_L,thetphi_tor,raysbylauncher,powers,fname,z_axis,launch_params] = ...
+    toray_raygeom_TCV(shot,pow.ec.t(1:10:end));
+
+for iant=1:nb_antennas
+  ids_ec_antennas.antenna{iant}.identifier = ['L' num2str(iant)];
+  switch iant
+   case {1,4}
+    ids_ec_antennas.antenna{iant}.name = 'LFS equatorial launcher';
+   case {2, 3, 5, 6}
+    ids_ec_antennas.antenna{iant}.name = 'LFS top launcher';
+   otherwise
+    ids_ec_antennas.antenna{iant}.name = 'TOP launcher';
+  end
+  if any(isfinite(pow.ec.data(:,iant)))
+    ids_ec_antennas.antenna{iant}.power_launched.data = pow.ec.data(:,iant);
+    ids_ec_antennas.antenna{iant}.power_launched.time = pow.ec.t;
+    ids_ec_antennas_description.antenna{iant}.power_launched = 'from gdat powers .ec';
+    ids_ec_antennas_description.antenna{iant}.launching_position = 'from launch_params obtained from toray_raygeom_TCV';
+  else
+    ids_ec_antennas_description.antenna{iant}.power_launched = 'no power for this launcher from gdat powers .ec';
+  end
+  if length(launch_params)>=iant && ~isempty(launch_params{iant})
+    % find 1st non-empty data
+    it_ok{iant} = [];
+    for it=1:length(launch_params{iant})
+      if ~isempty(launch_params{iant}{it})
+        it_ok{iant}(end+1) = it;
+      end
+    end
+    % non time-dependent quantities, take 1st ok values
+    ids_ec_antennas.antenna{iant}.frequency =launch_params{iant}{it_ok{iant}(1)}.freq;
+    ids_ec_antennas.antenna{iant}.mode.time = [pow.ec.t(1) pow.ec.t(end)];
+    ids_ec_antennas.antenna{iant}.mode.data = [-1 -1]; % at this stage assume X mode always, to change when available
+    for i=1:length(it_ok{iant})
+      r0 = sqrt(launch_params{iant}{it_ok{iant}(i)}.x0.^2 + launch_params{iant}{it_ok{iant}(i)}.y0.^2) / 100.; % in [m]
+      time_launch = launch_params{iant}{it_ok{iant}(i)}.time;
+      ids_ec_antennas.antenna{iant}.launching_position.r.data(i) = r0;
+      ids_ec_antennas.antenna{iant}.launching_position.z.data(i) = launch_params{iant}{it_ok{iant}(i)}.z0/100.;
+      ids_ec_antennas.antenna{iant}.launching_position.phi.data(i) = atan2(launch_params{iant}{it_ok{iant}(i)}.y0/100,r0);
+      ids_ec_antennas.antenna{iant}.launching_position.r.time(i) = time_launch;
+      ids_ec_antennas.antenna{iant}.launching_position.z.time(i) = time_launch;
+      ids_ec_antennas.antenna{iant}.launching_position.phi.time(i) = time_launch;
+      ids_ec_antennas.antenna{iant}.launching_angle_pol.data(i) = launch_params{iant}{it_ok{iant}(i)}.theta_toray * pi/180.;
+      ids_ec_antennas.antenna{iant}.launching_angle_pol.time(i) = time_launch;
+      ids_ec_antennas.antenna{iant}.launching_angle_tor.data(i) = launch_params{iant}{it_ok{iant}(i)}.phi_toray * pi/180.;
+      ids_ec_antennas.antenna{iant}.launching_angle_tor.time(i) = time_launch;
+    end
+  end
+end

matlab/TCV_IMAS/tcv_get_ids_nbi.m

+function [ids_nbi,ids_nbi_description,varargout] = tcv_get_ids_nbi(shot,ids_nbi_empty, gdat_params,varargin);
+%
+%  [ids_nbi,ids_nbi_description,varargout] = tcv_get_ids_nbi(shot,ids_nbi_empty,varargin);
+%
+%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar options
+%
+
+[ids_nbi, params_nbi] = tcv_ids_headpart(shot,ids_nbi_empty,'nbi','homogeneous_time',0,varargin{:});
+ids_nbi_description='';
+% As a general rule, for a new substructure under the main ids, construct a local structure like:
+% "global_quantities" with subfields being the relevant data to get and a local structure:
+% "global_quantities_desc" which contains the same subfields themselves containing the gdat string aftre shot used
+%
+
+
+nb_units = 2; % assume 2 units: 1st NBH and DNBI
+ids_nbi.unit(1:nb_units) = ids_nbi.unit(1); % copy empty structure for all units, then fill in
+
+% create lists of what is different for each units so that can scan through units
+unit_identifier = {'NBH1', 'DNBI'};
+unit_name = {'25keV 1st NBH source', 'diagnostic NBI'};
+results_subname = {'nbh', 'dnbi'};
+species.a = [2., 1.];
+species.z_n = [1., 1.];
+species.label = {'D', 'H'};
+beamlets_group.direction = [-1, 1];
+beamlets_group.tangency_radius = [736, 235.3]*1e-3;
+beamlets_group.angle = [0., 0.];
+beamlets_group.width_horizontal = [250, 87.2]*1e-3;
+beamlets_group.width_vertical = [250, 87.2]*1e-3;
+
+beamlets_group.focus(1:nb_units)=struct('focal_length_horizontal',[],'focal_length_vertical',[],'width_min_horizontal',[],'width_min_vertical',[]);
+beamlets_group.focus(1).focal_length_horizontal = 3.76;
+beamlets_group.focus(1).focal_length_vertical = 3.98;
+beamlets_group.focus(1).width_min_horizontal = 21.6*1e-2;
+beamlets_group.focus(1).width_min_vertical = 9.4*1e-2;
+beamlets_group.focus(2).focal_length_horizontal = 1.8;
+beamlets_group.focus(2).focal_length_vertical = 1.8;
+beamlets_group.focus(2).width_min_horizontal = 12.1*1e-2;
+beamlets_group.focus(2).width_min_vertical = 12.1*1e-2;
+
+beamlets_group.divergence(1:nb_units) = struct('particle_fraction',[],'vertical',[],'horizontal',[]);
+beamlets_group.divergence(1).particle_fraction = 1.;
+beamlets_group.divergence(1).vertical = 0.59 *pi/180.;
+beamlets_group.divergence(1).horizontal = 1.4 *pi/180.;
+beamlets_group.divergence(2).particle_fraction = 1.;
+beamlets_group.divergence(2).vertical = 0.53 *pi/180.;
+beamlets_group.divergence(2).horizontal = 0.53 *pi/180.;
+
+%dcd_NBH = psitbxdcd(4.5889, 0.0, 211.9535*pi/180, 0.0, -9.2308*pi/180);
+beamlets_group.position(1:nb_units) = struct('phi',[],'r',[],'z',[]);
+beamlets_group.position(1).phi = 211.9535*pi/180.;
+beamlets_group.position(1).r = 4.5889*1e-2; % in [m]?? OS
+beamlets_group.position(1).z = 0.; % in [m]?? OS
+beamlets_group.position(2).phi = 295.2416*pi/180.;
+beamlets_group.position(2).r = 4.9274*1e-2; % in [m]?? OS
+beamlets_group.position(2).z = 0.; % in [m]?? OS
+
+for iunit=1:nb_units
+  ids_nbi.unit{iunit}.identifier = unit_identifier{iunit};
+  ids_nbi.unit{iunit}.name = unit_name{iunit};
+  %% power
+  pow=gdat_tcv(shot,['\results::' results_subname{iunit} ':powr_tcv']);
+  ids_nbi.unit{iunit}.power_launched.data = pow.data*1e6;
+  ids_nbi.unit{iunit}.power_launched.time = pow.t;
+  ids_nbi_description.unit{iunit}.power_launched = ['from \results::' results_subname{iunit} ':powr_tcv'];
+  %% energy
+  en=gdat_tcv(shot,['\results::' results_subname{iunit} ':energy']);
+  ids_nbi.unit{iunit}.energy.data = en.data*1e3;
+  ids_nbi.unit{iunit}.energy.time = en.t;
+  ids_nbi_description.unit{iunit}.energy = ['from \results::' results_subname{iunit} ':energy'];
+  %% power & current fractions
+  p_frac=gdat(shot,['\results::' results_subname{iunit} ':fraction']);
+  ids_nbi.unit{iunit}.beam_power_fraction.data = p_frac.data(:,1:3)*0.01;
+  ids_nbi.unit{iunit}.beam_power_fraction.time = p_frac.t;
+  ids_nbi_description.unit{iunit}.beam_power_fraction = ['from \results::' results_subname{iunit} ':fraction'];
+  i_frac = p_frac.data(:,1:3).*repmat([1 2 3],size(p_frac.data,1),1); % to be compatible with older matlab version .*[1 2 3] not ok
+  i_frac = i_frac.*1./repmat(sum(i_frac),size(p_frac.data,1),1);
+  ids_nbi.unit{iunit}.beam_current_fraction.data = i_frac;
+  ids_nbi.unit{iunit}.beam_current_fraction.time = p_frac.t;
+  %% species
+  ids_nbi.unit{iunit}.species.a   = species.a(iunit);
+  ids_nbi.unit{iunit}.species.z_n = species.z_n(iunit);
+  ids_nbi.unit{iunit}.species.label = species.label{iunit};
+  %% beamlets group, now only one single beamlet
+  % https://spcwiki.epfl.ch/wiki/NB_Model
+  ids_nbi.unit{iunit}.beamlets_group{1}.direction = beamlets_group.direction(iunit); %clockwise
+  ids_nbi.unit{iunit}.beamlets_group{1}.tangency_radius = beamlets_group.tangency_radius(iunit);
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.tangency_radius_error_index: -999999999
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.tangency_radius_error_lower: -9.0000e+40
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.tangency_radius_error_upper: -9.0000e+40
+  ids_nbi.unit{iunit}.beamlets_group{1}.angle = beamlets_group.angle(iunit); %injection parallel to midplane
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.angle_error_index: -999999999
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.angle_error_lower: -9.0000e+40
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.angle_error_upper: -9.0000e+40
+  ids_nbi.unit{iunit}.beamlets_group{1}.width_horizontal = beamlets_group.width_horizontal(iunit);
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.width_horizontal_error_index: -999999999
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.width_horizontal_error_lower: -9.0000e+40
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.width_horizontal_error_upper: -9.0000e+40
+  ids_nbi.unit{iunit}.beamlets_group{1}.width_vertical = beamlets_group.width_vertical(iunit);
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.width_vertical_error_index: -999999999
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.width_vertical_error_lower: -9.0000e+40
+  %     ids_nbi.unit{iunit}.beamlets_group{1}.width_vertical_error_upper: -9.0000e+40
+  ids_nbi.unit{iunit}.beamlets_group{1}.focus = beamlets_group.focus(iunit);
+  %% divergence component struct
+  ids_nbi.unit{iunit}.beamlets_group{1}.divergence_component{1} = beamlets_group.divergence(iunit);
+  %% tilting
+  % it is fixed in time. what should we do about it?
+  %% position
+  ids_nbi.unit{iunit}.beamlets_group{1}.position = beamlets_group.position(iunit);
+
+  %% beamlets
+  % M. Vallar thinks it is useless now
+
+end