Update conversion to j_parallel and integral

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

@@ -123,7 +123,19 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
     ids_core_sources.source{last_index+i_lau}.global_quantities = globals_template;
   end
 
-  % get geometrical parameters and conversion from j_V,TORAY to j// = <jdotB>/B0 
+  % 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
+    ids_core_sources.source{last_index+i_lau}.global_quantities.time = ec_powers_tgrid;
+    ids_core_sources.source{last_index+i_lau}.global_quantities.power = ec_total_pow;
+    ids_core_sources.source{last_index+i_lau}.global_quantities.current_parallel = ec_total_cur;
+  end
+
+    % get geometrical parameters and conversion from j_V,TORAY to j// = <jdotB>/B0
   % from LIUQE
   [L,~,LY] = liuqe(shot,ec_data_tgrid,'iterq',50,'ilim',3,'icsint',true);
   rho_pol_norm_liu = L.pQ; % normalized sqrt poloidal flux including axis
@@ -132,27 +144,33 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
   Vprime = -1./LY.Q1Q; % Q1Q: -dpsi/dV
   Ip_psi = LY.ItQ; % Ip(psi) = I_tor(psi)
   T      = LY.TQ;
+  vol    = LY.VQ;
   B2_fs  = mu0*Ip_psi.*2.*pi./Vprime + T.^2.*Rm2_fs;
   toray_to_parallel = 1./(T.*Rm2_fs).*B2_fs./R0./B0;
-    
-  % 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
-    ids_core_sources.source{last_index+i_lau}.global_quantities.time = ec_powers_tgrid;
-    ids_core_sources.source{last_index+i_lau}.global_quantities.power = ec_total_pow;
-    ids_core_sources.source{last_index+i_lau}.global_quantities.current_parallel = ec_total_cur;
+  % interpolate on ec_data rho_pol grid
+  toray_to_parallel_ec = nan(size(toray_to_parallel));
+  vol_ec = nan(size(toray_to_parallel));
+  for ii = 1:nt_ec_data
+    toray_to_parallel_ec(:,ii) = ...
+      interp1(rho_pol_norm_liu,toray_to_parallel(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+    vol_ec(:,ii) = ...
+      interp1(rho_pol_norm_liu,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
   end
 
   % 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;
-
+  cd_dens       = nan(size(ec_data.cd_dens.data));
+  cd_integrated = nan(size(ec_data.cd_integrated.data));
+  %%
+  for i_lau = active_launchers
+    cd_dens(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*toray_to_parallel_ec;
+    for ii = 1:nt_ec_data
+      cd_integrated(:,i_lau,ii) = cumtrapz(vol_ec(:,ii),cd_dens(:,i_lau,ii))/(2*pi*0.88);
+    end
+  end
+%%
   p_ec_injected = powers_gdat.ec.data;
   ij = iround_os(ec_powers_tgrid,ec_data_tgrid);
   sparse_p_ec_injected = p_ec_injected(ij,:); % injected ec power vals corresponding to ec_data_tgrid