diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index e6c59afd1c112564f0b6857ec5d693afd6ea10e8..9f81305df8cda80c047fc7d4bd2d1a50d2397fce 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -40,14 +40,14 @@ ids_core_sources.time = b0_gdat.t;
source_template = ids_core_sources.source{1};
profiles_template = source_template.profiles_1d{1};
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);
params_eff = params_eff_ref; params_eff.data_request='powers';
powers_gdat = gdat(params_core_sources.shot,params_eff);
-% fix time grid to ohm time grid, same as bs
-
-last_index = 0;
%% initialize source from template
% ohm
ohm_t_grid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_t_grid);
@@ -60,40 +60,79 @@ ids_core_sources.source{last_index+1}.profiles_1d(1:ohm_n_t) = {profiles_templat
ids_core_sources.source{last_index+1}.global_quantities = globals_template;
ohm_data = ohm_gdat.ohm.ohm_data;
-%fill profiles for times n t_grid
-
+% load LIUQE data to convert
+% \tilde{j}_// = <jdotB>/(R0<Bphi/R>) to j_//0 = <jdotB>/B0; Bphi=T(psi)/R
+[L,~,LY] = liuqe(shot,ohm_t_grid,'iterq',50,'ilim',3,'icsint',true);
+rho_pol_norm_liu = L.pQ; % normalized sqrt poloidal flux including axis
+T = LY.TQ;
+Rm2_fs = LY.Q2Q;
+vol = LY.VQ;
+R0 = ids_core_sources.vacuum_toroidal_field.r0;
+B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ohm_t_grid);
+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 = interp1(rho_pol_norm_liu,jpar_tilde_to_jpar_0,ohm_data.cd_dens.rhopol_norm');
+vol_mapped = interp1(rho_pol_norm_liu,vol,ohm_data.cd_dens.rhopol_norm');
+
+% profiles_1d
+current_parallel_tmp = nan(1,ohm_n_t);
for ii = 1:ohm_n_t
- % ohm
- % 1d_profiles
ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = ohm_t_grid(ii);
- ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = ...
- ohm_data.cd_dens.data(:,ii)';
+ ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_pol_norm = ...
+ ohm_data.cd_dens.rhopol_norm';
+ ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_tor_norm = ...
+ ohm_data.cd_dens.rhotor_norm(:,ii)';
+ jpar_0_tmp = ohm_data.cd_dens.data(:,ii)'.*jpar_tilde_to_jpar_0_mapped(:,ii)';
+ ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = jpar_0_tmp;
+ integrated_jpar_0_tmp = cumtrapz(vol_mapped(:,ii),jpar_0_tmp)/2/pi/R0;
+ ids_core_sources.source{last_index+1}.profiles_1d{ii}.current_parallel_inside = integrated_jpar_0_tmp;
+ current_parallel_tmp(ii) = integrated_jpar_0_tmp(end);
end
- % globals
+
+% globals
ids_core_sources.source{last_index+1}.global_quantities.time = ohm_t_grid;
ids_core_sources.source{last_index+1}.global_quantities.power = powers_gdat.ohm.data;
+ids_core_sources.source{last_index+1}.global_quantities.current_parallel = current_parallel_tmp;
last_index = last_index+1; % add if statement to only increment if ohmic source has been added
%% bs
-bs_gdat = gdat(shot,'bs_data'); bs_data = bs_gdat.bs.bs_data;
+params_eff = params_eff_ref; params_eff.data_request='bs_data';
+bs_gdat = gdat(params_core_sources.shot,params_eff); bs_data = bs_gdat.bs.bs_data;
bs_t_grid = bs_gdat.bs.t; bs_n_t = numel(bs_t_grid);
+if any(bs_t_grid-ohm_t_grid')>0
+ warning('Bootstrap and ohmic time grids are not the same! Interpolation needed.')
+end
+
id_bs.description = 'Bootstrap current';
id_bs.index = 13; id_bs.name = 'bootstrap_current';
ids_core_sources.source{last_index+1} = source_template;
ids_core_sources.source{last_index+1}.identifier = id_bs;
-ids_core_sources.source{last_index+1}.profiles_1d(1:ohm_n_t) = {profiles_template};
+ids_core_sources.source{last_index+1}.profiles_1d(1:bs_n_t) = {profiles_template};
ids_core_sources.source{last_index+1}.global_quantities = globals_template;
-
-% fill profiles for times n t_grid
+% profiles_1d
+current_parallel_tmp = nan(1,bs_n_t);
for ii = 1:bs_n_t
ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = bs_t_grid(ii);
- ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = ...
- bs_data.cd_dens.data(:,ii)';
+ ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_pol_norm = ...
+ bs_data.cd_dens.rhopol_norm';
+ ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_tor_norm = ...
+ bs_data.cd_dens.rhotor_norm(:,ii)';
+ % use same conversion as for ohmic data from LIUQE
+ jpar_0_tmp = bs_data.cd_dens.data(:,ii)'.*jpar_tilde_to_jpar_0_mapped(:,ii)';
+ ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = jpar_0_tmp;
+ integrated_jpar_0_tmp = cumtrapz(vol_mapped(:,ii),jpar_0_tmp)/2/pi/R0;
+ ids_core_sources.source{last_index+1}.profiles_1d{ii}.current_parallel_inside = integrated_jpar_0_tmp;
+ current_parallel_tmp(ii) = integrated_jpar_0_tmp(end);
end
+% globals
+ids_core_sources.source{last_index+1}.global_quantities.time = bs_t_grid;
+ids_core_sources.source{last_index+1}.global_quantities.current_parallel = current_parallel_tmp;
+
last_index = last_index+1; % add if statement to only increment if bs source has been added
%% ec
@@ -135,8 +174,8 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
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
+ % load LIUQE data to convert
+ % \tilde{j}_// = <jdotB>/(R0<Bphi/R>) to j_//0 = <jdotB>/B0; Bphi=T(psi)/R
[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
mu0 = 4*pi*10^-7;