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
%

if exist('gdat_params')
  [ids_nbi, params_nbi] = tcv_ids_headpart(shot,ids_nbi_empty,'nbi','homogeneous_time',0,'gdat_params',gdat_params,varargin{:});
else
  [ids_nbi, params_nbi] = tcv_ids_headpart(shot,ids_nbi_empty,'nbi','homogeneous_time',0,varargin{:});
  aa=gdat_tcv;
  gdat_params = aa.gdat_params; % to get default params
end
params_eff_ref = gdat_params; params_eff_ref.doplot=0;
try;params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def)

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;
beamlets_group.position(1).z = 0.;
beamlets_group.position(2).phi = 295.2416*pi/180.;
beamlets_group.position(2).r = 4.9274;
beamlets_group.position(2).z = 0.;

params_eff = params_eff_ref;
for iunit=1:nb_units
  ids_nbi.unit{iunit}.identifier = unit_identifier{iunit};
  ids_nbi.unit{iunit}.name = unit_name{iunit};
  %% power
  params_eff.data_request = ['\results::' results_subname{iunit} ':powr_tcv'];
  pow=gdat_tcv(shot,params_eff);
  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 = params_eff.data_request;
  %% energy
  params_eff.data_request = ['\results::' results_subname{iunit} ':energy'];
  en=gdat_tcv(shot,params_eff);
  ids_nbi.unit{iunit}.energy.data = en.data*1e3;
  ids_nbi.unit{iunit}.energy.time = en.t;
  ids_nbi_description.unit{iunit}.energy = params_eff.data_request;
  %% power & current fractions
  params_eff.data_request = ['\results::' results_subname{iunit} ':fraction'];
  p_frac=gdat(shot,params_eff);
  ids_nbi.unit{iunit}.beam_power_fraction.time = p_frac.t;
  ids_nbi_description.unit{iunit}.beam_power_fraction = params_eff.data_request;
  if ~isempty(p_frac.data) && size(p_frac.data,2)>=3
    ids_nbi.unit{iunit}.beam_power_fraction.data = p_frac.data(:,1:3)'*0.01;
    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);
  else
    ids_nbi.unit{iunit}.beam_power_fraction.data = p_frac.data;
    i_frac = p_frac.data;
  end
  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


% cocos automatic transform
if exist('ids_generic_cocos_nodes_transformation_symbolic') == 2
  [ids_nbi,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_nbi,'nbi',gdat_params.cocos_in, ...
          gdat_params.cocos_out,gdat_params.ipsign_out,gdat_params.b0sign_out,gdat_params.ipsign_in,gdat_params.b0sign_in, ...
          gdat_params.error_bar,gdat_params.nverbose);
end