Clean, modularize and add interpolation when the time array is completely different

matlab/TCV_IMAS/tcv_get_ids_core_profiles.m

@@ -98,6 +98,9 @@ params_eff.data_request='ip';
 global_quantities.ip = gdat(params_cores_profiles.shot,'ip','machine',machine);
 global_quantities_desc.ip = params_eff.data_request;
 
+params_eff.data_request = '\results::ibs:signeo';
+signeo = gdat(params_cores_profiles.shot,params_eff);
+
 params_eff.data_request = 'ohm_data';
 current_non_inductive = gdat(params_cores_profiles.shot,params_eff);
 global_quantities.current_non_inductive.data = current_non_inductive.ohm.ohm_data.cd_tot.data;
@@ -359,95 +362,46 @@ for it=1:length(ids_core_profiles.time)
 end
 temp_1d_desc.magnetic_shear = 'from interpos with rhotor';
 
-
 %% Current densities
-data_unavailable_indexes = [];
-for it=1:numel(ids_core_profiles.time)
-  % ids_core_profiles.time might be different from current_bootstrap.ohm.ohm_data.cd_tot.t.
-  % If dimensions are different, a warning will be issued and interpolation will be performed
-  if it > size(current_non_inductive.ohm.ohm_data.cd_dens.data,2)
-    data_unavailable_indexes = [data_unavailable_indexes, it];
-  % Do the same when cd_dens.data is filled with nans. Cannot be an or on previous statement
-  elseif all(isnan(current_non_inductive.ohm.ohm_data.cd_dens.data(:,it)))
-    data_unavailable_indexes = [data_unavailable_indexes, it];
-    %nans as place holder, will be substituted later
-    ids_core_profiles.profiles_1d{it}.j_ohmic = current_non_inductive.ohm.ohm_data.cd_dens.data(:,it);
-  else
-    ids_core_profiles.profiles_1d{it}.j_ohmic = current_non_inductive.ohm.ohm_data.cd_dens.data(:,it);
-  end
-end
-temp_1d_desc.j_ohmic = current_non_inductive.ohm.ohm_data.cd_dens.label;
 
-if ~isempty(data_unavailable_indexes)
-  % Create new data interpolating on time.
-  new_x_slices = create_new_x_slices(data_unavailable_indexes, current_non_inductive.ohm.ohm_data.cd_dens.t, ...
-                                     current_non_inductive.ohm.ohm_data.cd_dens.data, ids_core_profiles.time, tens_time);
-  for ii=1:numel(data_unavailable_indexes)
-    ids_core_profiles.profiles_1d{data_unavailable_indexes(ii)}.j_ohmic = new_x_slices(ii,:);
-  end
-  warning('Ohmic current density data not available for time(s) %s interpolating instead', sprintf(' %3f,', ids_core_profiles.time(data_unavailable_indexes)));
-  temp_1d_desc.j_ohmic = string(temp_1d_desc.j_ohmic) + ' and interpolation on time(s)' + sprintf(' %3f,', ids_core_profiles.time(data_unavailable_indexes));
-end
+% Take care of Nans and differences of one element if present
+% If the time arrays are completely different, rebase. Still assume the same xgrid
+is_almost_equal = compare_arrays_tolength(current_non_inductive.ohm.ohm_data.cd_dens.t, ids_core_profiles.time);
+if ~is_almost_equal
+  [ids_core_profiles, temp_1d_desc] = insert_current_density( ...
+      ids_core_profiles, current_non_inductive.ohm.ohm_data, ...
+      'cd_dens', 'j_ohmic', 'Ohmic current', tens_time, temp_1d_desc);
+else
+  current_non_inductive.ohm.ohm_data.cd_dens.data = rebase_by_time(current_non_inductive.ohm.ohm_data.cd_dens.t, ...
+        current_non_inductive.ohm.ohm_data.cd_dens.data, ids_core_profiles.time, tens_time);
 
-data_unavailable_indexes = [];
-for it=1:numel(ids_core_profiles.time)
-  % ids_core_profiles.time might be different from current_bootstrap.bs.bs_data.cd_tot.t.
-  % If dimensions are different, a warning will be issued and interpolation will be performed
-  if it > size(current_bootstrap.bs.bs_data.cd_dens.data,2)
-    data_unavailable_indexes = [data_unavailable_indexes, it];
-  % Do the same when cd_dens.data is filled with nans. Cannot be an or on previous statement
-  elseif all(isnan(current_bootstrap.bs.bs_data.cd_dens.data(:,it)))
-    data_unavailable_indexes = [data_unavailable_indexes, it];
-    %nans as place holder, will be substituted later
-    ids_core_profiles.profiles_1d{it}.j_bootstrap = current_bootstrap.bs.bs_data.cd_dens.data(:,it);
-  else
-    ids_core_profiles.profiles_1d{it}.j_bootstrap = current_bootstrap.bs.bs_data.cd_dens.data(:,it);
-  end
+  [ids_core_profiles, temp_1d_desc] = fill_rebased_data(ids_core_profiles, temp_1d_desc, 'j_ohmic', current_non_inductive.ohm.ohm_data.cd_dens);
 end
-temp_1d_desc.j_bootstrap = current_bootstrap.bs.bs_data.cd_dens.label;
 
-if ~isempty(data_unavailable_indexes)
-  % Create new data interpolating on time.
-  new_x_slices = create_new_x_slices(data_unavailable_indexes, current_bootstrap.bs.bs_data.cd_dens.t, ...
-                                     current_bootstrap.bs.bs_data.cd_dens.data, ids_core_profiles.time, tens_time);
-  for ii=1:numel(data_unavailable_indexes)
-    ids_core_profiles.profiles_1d{data_unavailable_indexes(ii)}.j_bootstrap = new_x_slices(ii,:);
-  end
-  warning('Bootstrap current density data not available for time(s) %s interpolating instead', sprintf(' %3f,', ids_core_profiles.time(data_unavailable_indexes)));
-  temp_1d_desc.j_bootstrap = string(temp_1d_desc.j_bootstrap) + ' and interpolation on time(s)' + sprintf(' %3f,', ids_core_profiles.time(data_unavailable_indexes));
+is_almost_equal = compare_arrays_tolength(current_bootstrap.bs.bs_data.cd_dens.t, ids_core_profiles.time);
+if ~is_almost_equal
+  [ids_core_profiles, temp_1d_desc] = insert_current_density( ...
+      ids_core_profiles, current_non_inductive.ohm.ohm_data, ...
+      'cd_dens', 'j_bootstrap', 'Bootstrap current', tens_time, temp_1d_desc);
+else
+  current_bootstrap.bs.bs_data.cd_dens.data = rebase_by_time(current_bootstrap.bs.bs_data.cd_dens.t, ...
+        current_bootstrap.bs.bs_data.cd_dens.data, ids_core_profiles.time, tens_time);
+
+  [ids_core_profiles, temp_1d_desc] = fill_rebased_data(ids_core_profiles, temp_1d_desc, 'j_bootstrap', current_bootstrap.bs.bs_data.cd_dens);
 end
 
 %% parallel conductivity (neoclassical)
-signeo = gdat(params_cores_profiles.shot,'\results::ibs:signeo');
-data_unavailable_indexes = [];
-for it=1:numel(ids_core_profiles.time)
-  % ids_core_profiles.time might be different from signeo.t.
-  % If dimensions are different, a warning will be issued and interpolation will be performed
-  if it > size(signeo.data,2)
-    data_unavailable_indexes = [data_unavailable_indexes, it];
-  % Do the same when cd_dens.data is filled with nans. Cannot be an or on previous statement
-  elseif all(isnan(signeo.data(:,it)))
-    data_unavailable_indexes = [data_unavailable_indexes, it];
-    %nans as place holder, will be substituted later
-    ids_core_profiles.profiles_1d{it}.conductivity_parallel = signeo.data(:,it);
-  else
-    ids_core_profiles.profiles_1d{it}.conductivity_parallel = signeo.data(:,it);
-  end
-end
-temp_1d_desc.conductivity_parallel = signeo.label;
+is_almost_equal = compare_arrays_tolength(signeo.t, ids_core_profiles.time);
+if ~is_almost_equal
 
-if ~isempty(data_unavailable_indexes)
-  % Create new data interpolating on time.
-  new_x_slices = create_new_x_slices(data_unavailable_indexes, signeo.t', ...
-                                     signeo.data, ids_core_profiles.time, tens_time);
-  for ii=1:numel(data_unavailable_indexes)
-    ids_core_profiles.profiles_1d{data_unavailable_indexes(ii)}.conductivity_parallel = new_x_slices(ii,:);
-  end
-  warning('Parallel conductivity data not available for time(s) %s interpolating instead', sprintf(' %3f,', ids_core_profiles.time(data_unavailable_indexes)));
-  temp_1d_desc.conductivity_parallel = string(temp_1d_desc.conductivity_parallel) + ...
-                                       ' and interpolation on time(s)' + sprintf(' %3f,', ids_core_profiles.time(data_unavailable_indexes));
-end
+  [ids_core_profiles, temp_1d_desc] = insert_current_density( ...
+      ids_core_profiles, signeo, '', 'conductivity_parallel', 'Parallel conductivity', tens_time, temp_1d_desc);
+else
+  signeo.data = rebase_by_time(signeo.t', ...
+        signeo.data, ids_core_profiles.time, tens_time);
 
+  [ids_core_profiles, temp_1d_desc] = fill_rebased_data(ids_core_profiles, temp_1d_desc, 'conductivity_parallel', signeo);
+end
 
 %% add descriptions for profiles_1d
 ids_core_profiles_description.profiles_1d = temp_1d_desc;
@@ -482,7 +436,7 @@ function new_x_slices = create_new_x_slices(data_unavailable_indexes, time_old,
   for ix=1:numel(data(:,1));
     % Suppressing warnings for each xgrid, warnings depending on the time index will be given instead
     warning('off');
-    new_time_slice = interpos_nan(13, time_old, data(ix,:), time_new, tens_time);
+    new_time_slice = interpos_nan(13, time_old', data(ix,:), time_new', tens_time);
     warning('on');
     for ii=data_unavailable_indexes
       new_x_slices = [new_x_slices, new_time_slice(ii)];
@@ -490,3 +444,103 @@ function new_x_slices = create_new_x_slices(data_unavailable_indexes, time_old,
   end
   new_x_slices = reshape(new_x_slices, size(data_unavailable_indexes, 2), size(data, 1));
 end
+
+
+function [ids_core_profiles, temp_1d_desc] = fill_rebased_data(ids_core_profiles, temp_1d_desc, profile_field_name, struct);
+
+  for ii = 1:numel(ids_core_profiles.time)
+    ids_core_profiles.profiles_1d{ii}.(profile_field_name) = struct.data(ii,:);
+  end
+
+  temp_1d_desc.(profile_field_name) = strjoin(string(struct.label), ' interpolated in time');
+end
+
+function [ids_core_profiles, temp_1d_desc] = insert_current_density( ...
+    ids_core_profiles, current_struct, current_field_name, ...
+    profile_field_name, time_label, tens_time, temp_1d_desc);
+
+  if isempty(current_field_name)
+    data = current_struct.data;
+    time_vec = current_struct.t;
+    label = current_struct.label;
+  else
+    data = current_struct.(current_field_name).data;
+    time_vec = current_struct.(current_field_name).t;
+    label = current_struct.(current_field_name).label;
+  end
+
+  data_unavailable_indexes = [];
+  times = ids_core_profiles.time;
+
+  for it = 1:numel(times)
+    if it > size(data,2)
+      data_unavailable_indexes = [data_unavailable_indexes, it];
+    elseif all(isnan(data(:,it)))
+      data_unavailable_indexes = [data_unavailable_indexes, it];
+      ids_core_profiles.profiles_1d{it}.(profile_field_name) = data(:,it);
+    else
+      ids_core_profiles.profiles_1d{it}.(profile_field_name) = data(:,it);
+    end
+  end
+
+  temp_1d_desc.(profile_field_name) = label;
+
+  if ~isempty(data_unavailable_indexes)
+    new_x_slices = create_new_x_slices(data_unavailable_indexes, ...
+        time_vec(:), data, times, tens_time);
+
+    for ii = 1:numel(data_unavailable_indexes)
+      ids_core_profiles.profiles_1d{data_unavailable_indexes(ii)}.(profile_field_name) = ...
+          new_x_slices(ii,:);
+    end
+
+    warning('%s data not available for time(s)%s interpolating instead', ...
+        time_label, sprintf(' %3f,', times(data_unavailable_indexes)));
+
+    temp_1d_desc.(profile_field_name) = string(temp_1d_desc.(profile_field_name)) + ...
+        " and interpolation on time(s)" + sprintf(' %3f,', times(data_unavailable_indexes));
+  end
+end
+
+
+function interpolated_data = rebase_by_time(time_old, data, time_new, tens_time)
+  % Initialize new data
+  interpolated_data = zeros(size(data(:,1),1),length(time_new));
+  % Rebase the data by time, assuming constant xgrid
+  for ix = 1:numel(data(:,1))
+    warning('off');
+    new_time_slice = interpos_nan(13, time_old, data(ix,:), time_new', tens_time);
+    warning('on');
+    interpolated_data(ix,:) = new_time_slice;
+  end
+
+end
+
+% Checks if two arrays are equal or they only differ by the last value
+function is_almost_equal = compare_arrays_tolength(A, B)
+    % Ensure both are row vectors
+    A = A(:)';
+    B = B(:)';
+    
+    % Get lengths
+    lenA = length(A);
+    lenB = length(B);
+    
+    tolerance = 5e-6;
+
+    % If same length, compare all elements
+    if lenA == lenB
+        is_almost_equal = all(abs(A(:) - B(:)) < tolerance);
+        
+    % If lengths differ by 1, compare the shorter to the start of the longer
+    elseif abs(lenA - lenB) == 1
+        min_len = min(lenA, lenB);
+        is_almost_equal = all(abs(A(1:min_len) - B(1:min_len)) < tolerance);
+    else
+        is_almost_equal = false;
+    end
+end
+
+
+
+