diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m index aa08b0dc62481a45b82fc16fb839e71d1341960d..e31c6fad1d3c8ebd29c07cfbd71aebcc8fbef9a0 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m +++ b/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 + + + +