Update inputs for ohmic & bs source, do conversion from differen j_parallel...

Update inputs for ohmic & bs source, do conversion from differen j_parallel and add integrated current densities

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;