Clean up core_profiles ids and add jbs and johm

matlab/TCV_IMAS/tcv_get_ids_core_profiles.m

-function [ids_core_profiles,ids_core_profiles_description,varargout] = tcv_get_ids_core_profiles(shot,ids_equil_empty, gdat_params,varargin)
+function [ids_core_profiles,ids_core_profiles_description,varargout] = ...
+  tcv_get_ids_core_profiles(shot,ids_equil_empty,gdat_params,varargin)
 %
-%  [ids_core_profiles,ids_core_profiles_description,varargout] = tcv_get_ids_core_profiles(shot,ids_equil_empty,varargin);
+% [ids_core_profiles,ids_core_profiles_description,varargout] = ...
+%     tcv_get_ids_core_profiles(shot,ids_equil_empty,gdat_params,varargin);
 %
 %
-% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar options
-%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call, 
+%              in particular error_bar options
 %
 % error_bar: from gdat_params.error_bar
 %            'delta' (default): error_bar to be added inserted in "upper" only as mentioned in description
@@ -22,18 +24,22 @@ tens_time = -1;
 tens_rho = -0.1;
 
 if exist('gdat_params','var')
-  [ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles','gdat_params',gdat_params,varargin{:});
+  [ids_core_profiles, params_cores_profiles] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles','gdat_params',gdat_params,varargin{:});
 else
-  [ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',varargin{:});
+  [ids_core_profiles, params_cores_profiles] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',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)
 
+% initialize description
 ids_core_profiles_description = [];
 
-% base all from times fro nete_rho (which should be conf by default)
+%% setup profiles_1d
+% base all from times for nete_rho (which should be conf by default)
 params_eff = params_eff_ref;
 params_eff.data_request='ne_rho'; params_eff.fit = 1;
 temp_1d.ne_rho = gdat(params_cores_profiles.shot,params_eff);
@@ -70,13 +76,11 @@ ids_core_profiles.profiles_1d{1}.ion{1}.state = {};
 ids_core_profiles.profiles_1d{1}.neutral{1}.state = {};
 ids_core_profiles.profiles_1d(1:length(ids_core_profiles.time)) = ids_core_profiles.profiles_1d(1);
 
-
 % 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
-%
-% vacuum_toroidal_field and time, using homogeneous
-%
+
+%% vacuum_toroidal_field and time, using homogeneous
 params_eff = params_eff_ref;
 params_eff.data_request='b0';
 vacuum_toroidal_field.b0=gdat(params_cores_profiles.shot,params_eff);
@@ -85,7 +89,10 @@ ids_core_profiles.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
 ids_core_profiles.vacuum_toroidal_field.b0 = interpos(63,vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data,ids_core_profiles.time,-1);
 ids_core_profiles_description.vacuum_toroidal_field = [vacuum_toroidal_field_desc.b0 ' on ids_core_profiles.time, with interpos(63)'];
 
-% global_quantities data into local global_quantities.* structure with correct end names and global_quantities_desc.* with description. Use temp.* and temp_desc.* structures for temporary data
+%% global_quantities 
+% data into local global_quantities.* structure with correct end names and 
+% global_quantities_desc.* with description. 
+% Use temp.* and temp_desc.* structures for temporary data
 
 params_eff.data_request='ip';
 global_quantities.ip = gdat(params_cores_profiles.shot,'ip','machine',machine);
@@ -141,25 +148,20 @@ for i=1:length(global_quantities_fieldnames_eff)
     end
   end
 end
-% special case
-
 
-%
 %% profiles_1d (cannot use eqdsk since not same radial mesh)
+
+% prepare area for ids_core_profiles.profiles_1d{*}.grid.area
 params_eff.data_request = 'area_rho';
 temp_1d.area = gdat(params_cores_profiles.shot,params_eff);
 temp_1d_desc.area = params_eff.data_request;
-for ir=1:length(temp_1d.area.x)
+for ir=1:length(temp_1d.area.x) % map tmp_1d.area to core_profiles.time
   area_cpt(ir,:) = interpos(temp_1d.area.t,temp_1d.area.data(ir,:),ids_core_profiles.time,tens_time);
 end
-params_eff.data_request = 'q_rho';
-temp_1d.q = gdat(params_cores_profiles.shot,params_eff);
-temp_1d_desc.q = params_eff.data_request;
-for ir=1:length(temp_1d.q.x)
-  q_cpt(ir,:) = interpos(temp_1d.q.t,temp_1d.q.data(ir,:),ids_core_profiles.time,tens_time);
-end
+
 it_thom = iround_os(temp_1d.te_rho.t,ids_core_profiles.time);
 for it=1:length(ids_core_profiles.time)
+  % fill grid
   ids_core_profiles.profiles_1d{it}.grid.rho_tor_norm = temp_1d.fit.te_rho.grids_1d.rhotornorm(:,it);
   ids_core_profiles.profiles_1d{it}.grid.rho_tor = temp_1d.fit.te_rho.grids_1d.rhotornorm(:,it) ...
       .* temp_1d.fit.te_rho.grids_1d.rhotor_edge(it);
@@ -171,7 +173,9 @@ for it=1:length(ids_core_profiles.time)
       .* temp_1d.fit.te_rho.grids_1d.volume_edge(it);
   ids_core_profiles.profiles_1d{it}.grid.area = interpos(temp_1d.area.x,area_cpt(:,it),temp_1d.fit.te_rho.grids_1d.rhopolnorm, ...
           tens_rho,[1 2],[0 area_cpt(end,it)]);
+  % fill time
   ids_core_profiles.profiles_1d{it}.time = ids_core_profiles.time(it);
+  % fill electrons struct
   ids_core_profiles.profiles_1d{it}.electrons.temperature = temp_1d.fit.te_rho.data(:,it);
   ids_core_profiles.profiles_1d{it}.electrons.temperature_fit.measured = temp_1d.te_rho.data(:,it_thom(it));
   ids_core_profiles.profiles_1d{it}.electrons.temperature_fit.time_measurement = temp_1d.te_rho.t(it_thom(it));
@@ -185,6 +189,7 @@ for it=1:length(ids_core_profiles.time)
   ids_core_profiles.profiles_1d{it}.electrons.density_fit.source = {'Thomson, interpos fit'};
   ids_core_profiles.profiles_1d{it}.electrons.pressure_thermal = 1.6022e-19.*ids_core_profiles.profiles_1d{it}.electrons.density_thermal ...
       .* ids_core_profiles.profiles_1d{it}.electrons.temperature;
+  % fill zeff
   ids_core_profiles.profiles_1d{it}.zeff = global_quantities.z_eff_resistive.data(it) .* ...
       ones(size(ids_core_profiles.profiles_1d{it}.electrons.density));
 end
@@ -231,7 +236,8 @@ switch error_bar
   error(['tcv_ids_bpol_loop: error_bar option not known: ' error_bar])
 end
 
-% ion, assume only D if no CXRS (need to ask how to check if H...)
+%% ion struct
+% assume only D if no CXRS (need to ask how to check if H...)
 params_eff_fit1.data_request = 'cxrs';
 temp_1d.cxrs_rho = gdat(params_cores_profiles.shot,params_eff_fit1);
 temp_1d_desc.cxrs_rho = params_eff_fit1.data_request;
@@ -316,9 +322,14 @@ if ~isempty(temp_1d.cxrs_rho.ti.fit.data)
   end
 end
 
-for ir=1:length(temp_1d.q.x)
+%% q profile and magnetic shear
+params_eff.data_request = 'q_rho';
+temp_1d.q = gdat(params_cores_profiles.shot,params_eff);
+temp_1d_desc.q = params_eff.data_request;
+for ir=1:length(temp_1d.q.x) % map tmp_1d.q to core_profiles.time
   q_cpt(ir,:) = interpos(temp_1d.q.t,temp_1d.q.data(ir,:),ids_core_profiles.time,tens_time);
 end
+
 for it=1:length(ids_core_profiles.time)
   ij=isfinite(q_cpt(:,it));
   if sum(ij) >= 5
@@ -335,9 +346,21 @@ for it=1:length(ids_core_profiles.time)
     ids_core_profiles.profiles_1d{it}.magnetic_shear = -9.e40;
   end
 end
+temp_1d_desc.magnetic_shear = 'from interpos with rhotor';
+
+%% Current densities
+
+for it=1:length(ids_core_profiles.time)
+  ids_core_profiles.profiles_1d{it}.j_bootstrap = ...
+    current_bootstrap.bs.bs_data.cd_dens.data(:,it);
+  temp_1d_desc.j_bootstrap = current_bootstrap.bs.bs_data.cd_dens.label;
+  ids_core_profiles.profiles_1d{it}.j_ohmic = ...
+    current_non_inductive.ohm.ohm_data.cd_dens.data(:,it);
+  temp_1d_desc.j_ohmic = current_non_inductive.ohm.ohm_data.cd_dens.label;
+end
 
-ids_core_profiles_description.temp_1d_desc = temp_1d_desc;
-ids_core_profiles_description.profiles_1d.magnetic_shear = 'from interpos with rhotor';
+%% add descriptions for profiles_1d
+ids_core_profiles_description.profiles_1d = temp_1d_desc;
 
 if nargin <= 2
   ids_core_profiles.code.name = ['tcv_get_ids_core_profiles, within gdat, with shot= ' num2str(params_cores_profiles.shot) ];