diff --git a/matlab/TCV_IMAS/tcv2ids.m b/matlab/TCV_IMAS/tcv2ids.m
index 312cb4bc04e85f2f8b17efb95ded310f19b8234d..83662207267327b831b59dafa8fcb1c4cbcb589f 100644
--- a/matlab/TCV_IMAS/tcv2ids.m
+++ b/matlab/TCV_IMAS/tcv2ids.m
@@ -16,7 +16,7 @@ function [ids_from_tcv,varargout] = tcv2ids(shot,varargin);
 %                empty or 'delta' (default): meaning difference in upper is set (standard error_bar
 %                'added': errorbar is added: upper=data+delta and lower=data-delta
 %                'delta_with_lower': as 'delta' but lower also set
-%           'cocos_out': (default 11) cocos to transform ids from TCV cocos_in=17 to cocos_out
+%           'cocos_out': (DD3 and before default 11, DD4 default 17) cocos to transform ids from TCV cocos_in=17 to cocos_out
 %           'ipsign_out': if a specific sign fo Ip is desired in output within the cocos_out system (default 0=no specific sign)
 %           'b0sign_out': if a specific sign fo B0 is desired in output within the cocos_out system (default 0=no specific sign)
 %           'nverbose': (default 1), set it to 3 to have more messages, for example about not fully valid nodes when doing transformation (empty or Nans)
@@ -29,6 +29,16 @@ function [ids_from_tcv,varargout] = tcv2ids(shot,varargin);
 % varargout{1}: return also the ids in array of structure with the names, to allow easy use of plotallids
 %
 
+imas_version_number=getenv('IMAS_VERSION');
+cocos_out_default=17;
+try
+    if str2num(imas_version_number(1)) < 4
+        cocos_out_default=11;
+    end
+catch
+    disp('IMAS_VERSION undefined, assuming cocos_out_default=17');
+end
+
 % initialize input parser
 p = inputParser;
 % no required inputs here so one can call with empty args and get defaults parameters
@@ -37,7 +47,7 @@ ids_names = tcv_available_ids;
 p.addOptional('shot', [], @(x) (isnumeric(x) && isscalar(x) && (x == round(x)))); % integer
 p.addOptional('ids_names', ids_names, @(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('cocos_out', 11, @(x) isempty(x) || isnumeric(x) ); % char
+p.addOptional('cocos_out', cocos_out_default, @(x) isempty(x) || isnumeric(x) ); % char
 p.addOptional('ipsign_out', 0, @(x) isempty(x) || (isscalar(x) && (x==0 || x==-1 || x==+1)) );
 p.addOptional('b0sign_out', 0, @(x) isempty(x) || (isscalar(x) && (x==0 || x==-1 || x==+1)) );
 p.addOptional('nverbose', 1, @(x) isempty(x) || isnumeric(x) );
diff --git a/matlab/TCV_IMAS/tcv_get_ids_nbi.m b/matlab/TCV_IMAS/tcv_get_ids_nbi.m
index 9d5cf6a56080a327c02766f4cba32082639886f6..a9baab2b5a4bed15c9bef56f0de900878d7fb1af 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_nbi.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_nbi.m
@@ -5,6 +5,10 @@ function [ids_nbi,ids_nbi_description,varargout] = tcv_get_ids_nbi(shot,ids_nbi_
 %
 % gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar options
 %
+% Most information taken from NBI and NB_Model wiki pages:
+% https://spcwiki.epfl.ch/wiki/NBI
+% https://spcwiki.epfl.ch/wiki/NB_Model
+i_nbi1=1; i_nbi2=2; i_dnbi=3;
 
 if exist('gdat_params')
   [ids_nbi, params_nbi] = tcv_ids_headpart(shot,ids_nbi_empty,'nbi','homogeneous_time',0,'gdat_params',gdat_params,varargin{:});
@@ -22,17 +26,21 @@ ids_nbi_description='';
 % "global_quantities_desc" which contains the same subfields themselves containing the gdat string aftre shot used
 %
 
-
-nb_units = 3; % 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
-
+nb_units = 3; % assume 3 units: 1st NB1. 2nd NB2 and DNBI
 % create lists of what is different for each units so that can scan through units
 unit_identifier = {'NB1', 'NB2', 'DNBI'};
 unit_name = {'25keV 1st NBH source', '50keV 2nd NBH source', 'diagnostic NBI'};
 results_subname = {'nb1', 'nb2', 'dnbi'};
+
 if shot<70811
-  results_subname = {'nbh', 'nb2', 'dnbi'};
+  unit_identifier = {'NB1', 'DNBI'};
+  unit_name = {'25keV 1st NBH source', 'diagnostic NBI'};
+  results_subname = {'nbh','dnbi'};
+  nb_units = 2; % assume 2 units: 1st NBI and DNBI
 end
+nbi_indices=1:nb_units;
+ids_nbi.unit(nbi_indices) = ids_nbi.unit(1); % copy empty structure for all units, then fill in
+
 species.a = [2., 2., 1.];
 species.z_n = [1., 1., 1.];
 species.label = {'D', 'D','H'};
@@ -43,45 +51,45 @@ beamlets_group.width_horizontal = [250, 250, 87.2]*1e-3;
 beamlets_group.width_vertical = [250, 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(i_nbi1).focal_length_horizontal = 3.76;
+beamlets_group.focus(i_nbi1).focal_length_vertical = 3.98;
+beamlets_group.focus(i_nbi1).width_min_horizontal = 21.6*1e-2;
+beamlets_group.focus(i_nbi1).width_min_vertical = 9.4*1e-2;
 % So far NB2 parameters are merely a copy of NB1 parameters
-beamlets_group.focus(2).focal_length_horizontal = 3.76;
-beamlets_group.focus(2).focal_length_vertical = 3.98;
-beamlets_group.focus(2).width_min_horizontal = 21.6*1e-2;
-beamlets_group.focus(2).width_min_vertical = 9.4*1e-2;
-beamlets_group.focus(3).focal_length_horizontal = 1.8;
-beamlets_group.focus(3).focal_length_vertical = 1.8;
-beamlets_group.focus(3).width_min_horizontal = 12.1*1e-2;
-beamlets_group.focus(3).width_min_vertical = 12.1*1e-2;
+beamlets_group.focus(i_nbi2).focal_length_horizontal = 3.76;
+beamlets_group.focus(i_nbi2).focal_length_vertical = 3.98;
+beamlets_group.focus(i_nbi2).width_min_horizontal = 21.6*1e-2;
+beamlets_group.focus(i_nbi2).width_min_vertical = 9.4*1e-2;
+beamlets_group.focus(i_dnbi).focal_length_horizontal = 1.8;
+beamlets_group.focus(i_dnbi).focal_length_vertical = 1.8;
+beamlets_group.focus(i_dnbi).width_min_horizontal = 12.1*1e-2;
+beamlets_group.focus(i_dnbi).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.59 *pi/180.;
-beamlets_group.divergence(2).horizontal = 1.4 *pi/180.;
-beamlets_group.divergence(3).particle_fraction = 1.;
-beamlets_group.divergence(3).vertical = 0.53 *pi/180.;
-beamlets_group.divergence(3).horizontal = 0.53 *pi/180.;
+beamlets_group.divergence(i_nbi1).particle_fraction = 1.;
+beamlets_group.divergence(i_nbi1).vertical = 0.59 *pi/180.;
+beamlets_group.divergence(i_nbi1).horizontal = 1.4 *pi/180.;
+beamlets_group.divergence(i_nbi2).particle_fraction = 1.;
+beamlets_group.divergence(i_nbi2).vertical = 0.59 *pi/180.;
+beamlets_group.divergence(i_nbi2).horizontal = 1.4 *pi/180.;
+beamlets_group.divergence(i_dnbi).particle_fraction = 1.;
+beamlets_group.divergence(i_dnbi).vertical = 0.53 *pi/180.;
+beamlets_group.divergence(i_dnbi).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 = 58.8255*pi/180.;
-beamlets_group.position(2).r = 4.5889;
-beamlets_group.position(2).z = 0.;
-beamlets_group.position(2).phi = 295.2416*pi/180.;
-beamlets_group.position(2).r = 4.9274;
-beamlets_group.position(2).z = 0.;
+beamlets_group.position(i_nbi1).phi = 211.9535*pi/180.;
+beamlets_group.position(i_nbi1).r = 4.5889;
+beamlets_group.position(i_nbi1).z = 0.;
+beamlets_group.position(i_nbi2).phi = 58.8255*pi/180.;
+beamlets_group.position(i_nbi2).r = 4.5889;
+beamlets_group.position(i_nbi2).z = 0.;
+beamlets_group.position(i_dnbi).phi = 295.2416*pi/180.;
+beamlets_group.position(i_dnbi).r = 4.9274;
+beamlets_group.position(i_dnbi).z = 0.;
 
 params_eff = params_eff_ref;
-for iunit=1:nb_units
+for iunit=nbi_indices
   ids_nbi.unit{iunit}.identifier = unit_identifier{iunit};
   ids_nbi.unit{iunit}.name = unit_name{iunit};
   %% power
@@ -128,21 +136,9 @@ for iunit=1:nb_units
   % 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
 
   % Should always copy "leaves" only to avoid deleting non-filled in values like error_bars
   fields_to_copy = fieldnames(beamlets_group.focus(iunit));