From 5d8ffe686755dd16c946a74d22de3022da1488dd Mon Sep 17 00:00:00 2001
From: Antonia Frank <antonia.frank@epfl.ch>
Date: Thu, 23 May 2024 11:18:14 +0200
Subject: [PATCH] Cleanup code and spacing, add comments to make more clear
 what is happening

---
 matlab/TCV_IMAS/tcv_get_ids_core_sources.m | 159 +++++++++++----------
 1 file changed, 83 insertions(+), 76 deletions(-)

diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index f98b10b2..3504b3ac 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -88,122 +88,129 @@ last_index = last_index+1;  % add if statement to only increment if bs source ha
 
 %% ec
 % load data
-ec_gdat = gdat(shot,'ec_data'); ec_data = ec_gdat.ec.ec_data;
-ec_time = ec_data.p_dens.t;
-n_ec_time = numel(ec_time);
-ec_t_grid = powers_gdat.ec.t; ec_n_t = numel(ec_t_grid);
-ec_inputs = ec_gdat.ec.ec_inputs;
+ec_gdat = gdat(shot,'ec_data'); 
+ec_data = ec_gdat.ec.ec_data; ec_inputs = ec_gdat.ec.ec_inputs;
+% get tgrid from gdat ec_data
+ec_data_tgrid = ec_data.p_dens.t; nt_ec_data = numel(ec_data_tgrid);
 
+% get tgrid from gdat powers
+ec_powers_tgrid = powers_gdat.ec.t; nt_ec_powers = numel(ec_powers_tgrid);
+
+% retrieve active launcher information from ec_inputs
 nb_launchers = numel(ec_inputs.launchers_active.data);
 active_launchers = find(ec_inputs.launchers_active.data==1)';
 n_active_launchers = sum(ec_inputs.launchers_active.data);
 
+% Setup structures for active launchers from template
 main_desc = 'Source from electron cyclotron heating and current drive';
 id_ec.index = 3; id_ec.name = 'ec';
 for i_lau = active_launchers
   id_ec.description = sprintf('L%i/G%i, %s',i_lau,ec_inputs.gyro2launcher.data(i_lau),main_desc);
   ids_core_sources.source{last_index+i_lau} = source_template;
   ids_core_sources.source{last_index+i_lau}.identifier = id_ec;
-  ids_core_sources.source{last_index+i_lau}.profiles_1d(1:ec_n_t) = {profiles_template};
-  ids_core_sources.source{last_index+i_lau}.global_quantities(1:ec_n_t) = {globals_template};
+  ids_core_sources.source{last_index+i_lau}.profiles_1d(1:nt_ec_powers) = {profiles_template};
+  ids_core_sources.source{last_index+i_lau}.global_quantities(1:nt_ec_powers) = {globals_template};
 end
  
+% get data for globals from gdat powers and fill in ids_structure
 ec_total_pow = transpose(powers_gdat.ec.data(:,nb_launchers+1)); %use power from powers_gdat(injected ec power) instead of ec_data power
 ec_total_pow(isnan(ec_total_pow)) = 0;
 ec_total_cur = ec_data.cd_tot.data(nb_launchers+1,:);
 ec_total_cur(isnan(ec_total_cur)) = 0;
 
 for i_lau = active_launchers
-  for ii = 1:ec_n_t
-    ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = ec_t_grid(ii); % profiles time 
-    ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.time = ec_t_grid(ii); % globals time
-
-    % globals
-    ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.time = ec_t_grid(ii);
+  for ii = 1:nt_ec_powers
+    ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.time = ec_powers_tgrid(ii);
     ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.power = ec_total_pow(ii);
     ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.current_parallel = ec_total_cur(ii);
   end
 end
 
-% interpoating p_dens profiles from 'ec_time' grid (toray time nodes) to 'ec_t_grid' (injected power time)
-
-p_dens = ec_data.p_dens.data;
-p_integrated = ec_data.p_integrated.data;
-cd_dens = ec_data.cd_dens.data;
+% interpoating p_dens profiles from 'ec_data_tgrid' grid (toray tgrid) to 'ec_powers_tgrid' (powers tgrid)
+p_dens        = ec_data.p_dens.data;
+p_integrated  = ec_data.p_integrated.data;
+cd_dens       = ec_data.cd_dens.data;
 cd_integrated = ec_data.cd_integrated.data;
 
 p_ec_injected = powers_gdat.ec.data;
-it = iround_os(ec_t_grid,ec_time);
-sparse_p_ec_injected = p_ec_injected(it,:); % injected ec power vals corresponding to ec_time grid
-rho_grid = size(p_dens, 1);
-% calculate normalised profiles on ec_time grid
-norm_p_dens = zeros(rho_grid, n_active_launchers, n_ec_time);
-norm_p_integrated = zeros(rho_grid, n_active_launchers, n_ec_time);
-norm_cd_dens = zeros(rho_grid, n_active_launchers, n_ec_time);
-norm_cd_integrated = zeros(rho_grid, n_active_launchers, n_ec_time);
-for t = 1:n_ec_time
-  norm_p_dens_temp = zeros(rho_grid, n_active_launchers);
-  norm_p_integrated_temp = zeros(rho_grid, n_active_launchers);
-  norm_cd_dens_temp = zeros(rho_grid, n_active_launchers);
-  norm_cd_integrated_temp = zeros(rho_grid, n_active_launchers);
+it = iround_os(ec_powers_tgrid,ec_data_tgrid);
+sparse_p_ec_injected = p_ec_injected(it,:); % injected ec power vals corresponding to ec_data_tgrid
+n_rho = size(p_dens, 1);
+
+% calculate normalised profiles on ec_data_tgrid grid
+norm_p_dens        = zeros(n_rho,n_active_launchers,nt_ec_data);
+norm_p_integrated  = zeros(n_rho,n_active_launchers,nt_ec_data);
+norm_cd_dens       = zeros(n_rho,n_active_launchers,nt_ec_data);
+norm_cd_integrated = zeros(n_rho,n_active_launchers,nt_ec_data);
+for it = 1:nt_ec_data
+  norm_p_dens_temp        = zeros(n_rho,n_active_launchers);
+  norm_p_integrated_temp  = zeros(n_rho,n_active_launchers);
+  norm_cd_dens_temp       = zeros(n_rho,n_active_launchers);
+  norm_cd_integrated_temp = zeros(n_rho,n_active_launchers);
   for i_lau = n_active_launchers
-    norm_p_dens_temp(:, i_lau) = p_dens(:, active_launchers(i_lau), t) ./ sparse_p_ec_injected(t, active_launchers(i_lau));
-    norm_p_integrated_temp(:, i_lau) = p_integrated(:, active_launchers(i_lau)) ./ sparse_p_ec_injected(t, active_launchers(i_lau));
-    norm_cd_dens_temp(:, i_lau) = cd_dens(:, active_launchers(i_lau), t) ./ sparse_p_ec_injected(t, i_lau);
-    norm_cd_integrated_temp(:, i_lau) = cd_integrated(:, active_launchers(i_lau), t) ./ sparse_p_ec_injected(t, active_launchers(i_lau));
+    norm_p_dens_temp(:,i_lau)        = p_dens(:,active_launchers(i_lau),it) ./ sparse_p_ec_injected(it, active_launchers(i_lau));
+    norm_p_integrated_temp(:,i_lau)  = p_integrated(:,active_launchers(i_lau)) ./ sparse_p_ec_injected(it, active_launchers(i_lau));
+    norm_cd_dens_temp(:,i_lau)       = cd_dens(:,active_launchers(i_lau),it) ./ sparse_p_ec_injected(it, i_lau);
+    norm_cd_integrated_temp(:,i_lau) = cd_integrated(:,active_launchers(i_lau),it) ./ sparse_p_ec_injected(it, active_launchers(i_lau));
   end
-  norm_p_dens(:,:,t) = norm_p_dens_temp;
-  norm_p_integrated(:,:,t) = norm_p_integrated_temp;
-  norm_cd_dens(:,:,t) = norm_cd_dens_temp;
-  norm_cd_integrated(:,:,t) = norm_cd_integrated_temp;
+  norm_p_dens(:,:,it)        = norm_p_dens_temp;
+  norm_p_integrated(:,:,it)  = norm_p_integrated_temp;
+  norm_cd_dens(:,:,it)       = norm_cd_dens_temp;
+  norm_cd_integrated(:,:,it) = norm_cd_integrated_temp;
 end
 
-% interpolate normalised p_dens profiles on ec_t_grid 
-interp_norm_p_dens = zeros(rho_grid, n_active_launchers, ec_n_t);
-interp_norm_p_integrated = zeros(rho_grid, n_active_launchers, ec_n_t);
-interp_norm_cd_dens = zeros(rho_grid, n_active_launchers, ec_n_t);
-interp_norm_cd_integrated = zeros(rho_grid, n_active_launchers, ec_n_t);
-for rho = 1:rho_grid
+% interpolate normalised p_dens profiles on ec_powers_tgrid 
+interp_norm_p_dens        = zeros(n_rho,n_active_launchers,nt_ec_powers);
+interp_norm_p_integrated  = zeros(n_rho,n_active_launchers,nt_ec_powers);
+interp_norm_cd_dens       = zeros(n_rho,n_active_launchers,nt_ec_powers);
+interp_norm_cd_integrated = zeros(n_rho,n_active_launchers,nt_ec_powers);
+for irho = 1:n_rho
   for i_lau = 1:n_active_launchers
-    profile_p_dens = squeeze(norm_p_dens(rho, i_lau, :));
-    profile_p_integrated = squeeze(norm_p_integrated(rho, i_lau, :));
-    profile_cd_dens = squeeze(norm_cd_dens(rho, i_lau, :));
-    profile_cd_integrated = squeeze(norm_cd_integrated(rho, i_lau, :));
-    interp_profile_p_dens = interp1(ec_time, profile_p_dens, ec_t_grid);
-    interp_profile_p_integrated = interp1(ec_time, profile_p_integrated, ec_t_grid);
-    interp_profile_cd_dens = interp1(ec_time, profile_cd_dens, ec_t_grid);
-    interp_profile_cd_integrated = interp1(ec_time, profile_cd_integrated, ec_t_grid);
-    interp_norm_p_dens(rho, i_lau, :) = interp_profile_p_dens;
-    interp_norm_p_integrated(rho, i_lau, :) = interp_profile_p_integrated;
-    interp_norm_cd_dens(rho, i_lau, :) = interp_profile_cd_dens;
-    interp_norm_cd_integrated(rho, i_lau, :) = interp_profile_cd_integrated;
+    profile_p_dens        = squeeze(norm_p_dens(irho,i_lau,:));
+    profile_p_integrated  = squeeze(norm_p_integrated(irho,i_lau,:));
+    profile_cd_dens       = squeeze(norm_cd_dens(irho,i_lau,:));
+    profile_cd_integrated = squeeze(norm_cd_integrated(irho,i_lau,:));
+    
+    interp_profile_p_dens        = interp1(ec_data_tgrid, profile_p_dens, ec_powers_tgrid);
+    interp_profile_p_integrated  = interp1(ec_data_tgrid, profile_p_integrated, ec_powers_tgrid);
+    interp_profile_cd_dens       = interp1(ec_data_tgrid, profile_cd_dens, ec_powers_tgrid);
+    interp_profile_cd_integrated = interp1(ec_data_tgrid, profile_cd_integrated, ec_powers_tgrid);
+    
+    interp_norm_p_dens(irho,i_lau,:)        = interp_profile_p_dens;
+    interp_norm_p_integrated(irho,i_lau,:)  = interp_profile_p_integrated;
+    interp_norm_cd_dens(irho,i_lau,:)       = interp_profile_cd_dens;
+    interp_norm_cd_integrated(irho,i_lau,:) = interp_profile_cd_integrated;
   end
 end
 
-% normalised & interpolated profiles * p_ec_injected on ec_t_grid
-interp_p_dens = zeros(rho_grid, n_active_launchers, ec_n_t);
-interp_p_integrated = zeros(rho_grid, n_active_launchers, ec_n_t);
-interp_cd_dens = zeros(rho_grid, n_active_launchers, ec_n_t);
-interp_cd_integrated = zeros(rho_grid, n_active_launchers, ec_n_t);
-for t = 1:numel(ec_t_grid)
-  unnormalised_p_dens = zeros(rho_grid, n_active_launchers);
-  unnormalised_p_integrated = zeros(rho_grid, n_active_launchers);
-  unnormalised_cd_dens = zeros(rho_grid, n_active_launchers);
-  unnormalised_cd_integrated = zeros(rho_grid, n_active_launchers);
+% normalised & interpolated profiles * p_ec_injected on ec_powers_tgrid
+interp_p_dens        = zeros(n_rho,n_active_launchers,nt_ec_powers);
+interp_p_integrated  = zeros(n_rho,n_active_launchers,nt_ec_powers);
+interp_cd_dens       = zeros(n_rho,n_active_launchers,nt_ec_powers);
+interp_cd_integrated = zeros(n_rho,n_active_launchers,nt_ec_powers);
+for it = 1:numel(ec_powers_tgrid)
+  unnormalised_p_dens        = zeros(n_rho,n_active_launchers);
+  unnormalised_p_integrated  = zeros(n_rho,n_active_launchers);
+  unnormalised_cd_dens       = zeros(n_rho,n_active_launchers);
+  unnormalised_cd_integrated = zeros(n_rho,n_active_launchers);
   for i_lau = 1:n_active_launchers
-    unnormalised_p_dens(:,i_lau) = interp_norm_p_dens(:,i_lau,t).* p_ec_injected(t, active_launchers(i_lau));
-    unnormalised_p_integrated(:,i_lau) = interp_norm_p_integrated(:,i_lau,t).* p_ec_injected(t, active_launchers(i_lau));
-    unnormalised_cd_dens(:,i_lau) = interp_norm_cd_dens(:,i_lau,t).* p_ec_injected(t, active_launchers(i_lau));
-    unnormalised_cd_integrated(:,i_lau) = interp_norm_cd_integrated(:,i_lau,t).* p_ec_injected(t, active_launchers(i_lau));
+    unnormalised_p_dens(:,i_lau)        = interp_norm_p_dens(:,i_lau,it).* p_ec_injected(it, active_launchers(i_lau));
+    unnormalised_p_integrated(:,i_lau)  = interp_norm_p_integrated(:,i_lau,it).* p_ec_injected(it, active_launchers(i_lau));
+    unnormalised_cd_dens(:,i_lau)       = interp_norm_cd_dens(:,i_lau,it).* p_ec_injected(it, active_launchers(i_lau));
+    unnormalised_cd_integrated(:,i_lau) = interp_norm_cd_integrated(:,i_lau,it).* p_ec_injected(it, active_launchers(i_lau));
   end
-  interp_p_dens(:,:,t) = unnormalised_p_dens;
-  interp_p_integrated(:,:,t) = unnormalised_p_integrated;
-  interp_cd_dens(:,:,t) = unnormalised_cd_dens;
-  interp_cd_integrated(:,:,t) = unnormalised_cd_integrated;
+  interp_p_dens(:,:,it)        = unnormalised_p_dens;
+  interp_p_integrated(:,:,it)  = unnormalised_p_integrated;
+  interp_cd_dens(:,:,it)       = unnormalised_cd_dens;
+  interp_cd_integrated(:,:,it) = unnormalised_cd_integrated;
 end
 
-for ii = 1:ec_n_t
+% fill profiles_1d with interpolated profiles
+for ii = 1:nt_ec_powers
+  
   for i_lau = active_launchers
+    ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = ec_powers_tgrid(ii); % profiles time 
+    
     launcher_index = find(active_launchers == i_lau);
     ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = ...
       ec_data.p_dens.rhotor_norm(1,:);
-- 
GitLab