From 5e7a913fd3d1953de4ba8ef364f8c16ddef62895 Mon Sep 17 00:00:00 2001
From: Antonia Frank <antonia.frank@epfl.ch>
Date: Thu, 22 Aug 2024 11:05:04 +0200
Subject: [PATCH] make sure ohm tgrids of gdat powers and ohm data are the same

---
 matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index bd6a9e3b..4866c59c 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -43,9 +43,7 @@ globals_template = source_template.global_quantities{1};
 
 last_index = 0; % index to be incremented when sources are added
 
-%% get time grid and data from gdat
-params_eff = params_eff_ref; params_eff.data_request='ohm_data';
-ohm_gdat = gdat(params_core_sources.shot,params_eff);
+%% get powers information from gdat
 params_eff = params_eff_ref; params_eff.data_request='powers';
 powers_gdat = gdat(params_core_sources.shot,params_eff);
 
@@ -82,7 +80,12 @@ else
 end
 %% initialize source from template
 % ohm
-ohm_tgrid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_tgrid);
+params_eff = params_eff_ref; params_eff.data_request='ohm_data';
+ohm_gdat = gdat(params_core_sources.shot,params_eff); ohm_data = ohm_gdat.ohm.ohm_data;
+ohm_tgrid = ohm_gdat.ohm.t; ohm_n_t = numel(ohm_tgrid);
+
+% make sure that ohm data is on correct tgrid
+ohm_power = interpos(powers_gdat.ohm.t,powers_gdat.ohm.data,ohm_tgrid);
 
 main_desc = 'Source from ohmic heating'; production = 'IBS nodes';
 id_ohm.description = sprintf('%s from %s',main_desc,production);
@@ -91,7 +94,6 @@ ids_core_sources.source{last_index+1} = source_template;
 ids_core_sources.source{last_index+1}.identifier = id_ohm;
 ids_core_sources.source{last_index+1}.profiles_1d(1:ohm_n_t) = {profiles_template};
 ids_core_sources.source{last_index+1}.global_quantities(1:ohm_n_t) = {globals_template};
-ohm_data = ohm_gdat.ohm.ohm_data;
 
 % load LIUQE data to convert
 % \tilde{j}_// = <jdotB>/(R0<Bphi/R>) to j_//0 = <jdotB>/B0; Bphi=T(psi)/R
@@ -103,7 +105,7 @@ vol    = interp1(liuqe_time,vol_liu.data.',   ohm_tgrid)';
 R0 = ids_core_sources.vacuum_toroidal_field.r0;
 B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ohm_tgrid);
 
-jpar_tilde_to_jpar_0 = R0*T.*Rm2_fs./B0';
+jpar_tilde_to_jpar_0 = R0*T.*Rm2_fs./B0;
 
 % map volume and conversion factor on rho_pol grid of ohm_data cd_dens
 jpar_tilde_to_jpar_0_mapped = ...
@@ -124,7 +126,7 @@ for ii = 1:ohm_n_t
 
   % globals
   ids_core_sources.source{last_index+1}.global_quantities{ii}.time = ohm_tgrid(ii);
-  ids_core_sources.source{last_index+1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii);
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.power = ohm_power(ii);
   ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
 end
 
-- 
GitLab