diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index 1a3c531d034e868b0a595929cbe2b9c3591108ac..6f630945d23ad81bdf831e41bc201d4ec6e3602f 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -103,7 +103,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 = ...
@@ -206,33 +206,44 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
ids_core_sources.source{last_index+i_lau}.global_quantities(1:nt_ec_out) = {globals_template};
end
- % load LIUQE data to convert
- % \tilde{V,TORAY}_// = dI_\phi/dV = (1/2\pi)<j_\phi/R> to j_//0 = <jdotB>/B0;
- % conversion using j_//,cd = <j_cd/B> B, thus j_// = 2\pi <B^2> / (B_0 <B_\phi/R>)
- mu0 = 4*pi*10^-7;
- % interpolate liuqe outputs on ohm_tgrid
- T = interp1(liuqe_time,T_liu.data.', ec_tgrid_toray)';
- Rm2_fs = interp1(liuqe_time,Rm2_fs_liu.data.',ec_tgrid_toray)';
- vol = interp1(liuqe_time,vol_liu.data.', ec_tgrid_toray)';
- Ip = interp1(liuqe_time,Ip_liu.data, ec_tgrid_toray);
- q = interp1(liuqe_time,q_liu.data.', ec_tgrid_toray)';
- mVprime = -2*pi*q./T./Rm2_fs;
- B2_fs = mu0*Ip./mVprime + T.^2.*Rm2_fs;
-
% get vacuum field data from ids
R0 = ids_core_sources.vacuum_toroidal_field.r0;
B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ec_tgrid_toray);
-
- jtoray_to_jpar0 = 2*pi./B0./(T.*Rm2_fs).*B2_fs;
+
+ % do conversion j_{V,TORAY} = dI_\phi/dV = (1/2\pi)<j_\phi/R> to j_//0 = <jdotB>/B0;
+ % via j_tilde// = <jdotB>/(R0 <Bphi/R>) = 2*pi/R0/Rm2_fs[j_{V,TORAY}-dT/dV/T*I_{phi,cd}]
+
+ % interpolate liuqe outputs on ec_tgrid_toray
+ T = interp1(liuqe_time,T_liu.data.', ec_tgrid_toray)';
+ Rm2_fs = interp1(liuqe_time,Rm2_fs_liu.data.',ec_tgrid_toray)';
+ vol = interp1(liuqe_time,vol_liu.data.', ec_tgrid_toray)';
+
+% % alternative conversion, slightly mismatching with the conversion
+% % above, to be clarified
+% % using j_//,cd = <j_cd/B> B, thus j_// = 2\pi <B^2> / (B_0 <B_\phi/R>)
+% mu0 = 4*pi*10^-7;
+% q = interp1(liuqe_time,q_liu.data.', ec_tgrid_toray)';
+% Ip = interp1(liuqe_time,Ip_liu.data, ec_tgrid_toray);
+% Vprime = -2*pi*q./T./Rm2_fs;
+% B2_fs = -mu0*Ip./Vprime + T.^2.*Rm2_fs;
+% jtoray_to_jpar0 = 2*pi./B0./T./Rm2_fs.*B2_fs;
+% jtoray_to_jpar0_mapped = nan(size(ec_data.cd_dens.x));
% interpolate on ec_data rho_pol grid
- jtoray_to_jpar0_mapped = nan(size(ec_data.cd_dens.x));
- vol_mapped = nan(size(ec_data.cd_dens.x));
+ vol_mapped = nan(size(ec_data.cd_dens.x));
+ T_mapped = nan(size(ec_data.cd_dens.x));
+ dTdV_mapped = nan(size(ec_data.cd_dens.x));
+ Rm2_fs_mapped = nan(size(ec_data.cd_dens.x));
for ii = 1:nt_ec_toray
- jtoray_to_jpar0_mapped(:,ii) = ...
- interp1(rho_pol_norm_liu.data,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
- vol_mapped(:,ii) = ...
- interp1(rho_pol_norm_liu.data,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+% jtoray_to_jpar0_mapped(:,ii) = ...
+% interp1(rho_pol_norm_liu.data,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+ vol_mapped(:,ii) = interpos(...
+ rho_pol_norm_liu.data,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+ T_mapped(:,ii) = interpos(...
+ rho_pol_norm_liu.data,T(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+ [~,dTdV_mapped(:,ii)] = interpos(vol_mapped(:,ii),T_mapped(:,ii));
+ Rm2_fs_mapped(:,ii) = interpos(...
+ rho_pol_norm_liu.data,Rm2_fs(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
end
% load 1d_profiles from ec_data
@@ -242,14 +253,26 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
p_integrated(isnan(p_integrated)) = 0;
cd_dens = ec_data.cd_dens_doublewidth.data; % use double width to have realistic deposition broadening
cd_dens(isnan(cd_dens)) = 0;
-
- cd_integrated = zeros(size(ec_data.cd_integrated.data));
+ cd_integrated = ec_data.cd_integrated.data;
+ cd_integrated(isnan(cd_integrated)) = 0;
+
+ % initialize variables
+ jpar_tilde = zeros(size(cd_dens));
+ jpar0 = zeros(size(cd_dens));
+ jpar0_integrated = zeros(size(cd_dens));
+
+ % convert and integrate j_TORAY to j//0
for i_lau = active_launchers
- % convert to j// with conversion factor
- cd_dens(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
+% % alternative conversion
+% jpar0(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
+
+ jpar_tilde(:,i_lau,:) = 2*pi/R0./Rm2_fs_mapped.*(squeeze(cd_dens(:,i_lau,:))...
+ -dTdV_mapped./T_mapped.*squeeze(cd_integrated(:,i_lau,:)));
+ jpar0(:,i_lau,:) = squeeze(jpar_tilde(:,i_lau,:)).*R0.*T_mapped.*Rm2_fs_mapped./B0;
+
for ii = 1:nt_ec_toray
- % integrate j//
- cd_integrated(:,i_lau,ii) = cumtrapz(vol_mapped(:,ii),cd_dens(:,i_lau,ii))/(2*pi*R0);
+ % integrate j//0
+ [~,~,~,jpar0_integrated(:,i_lau,ii)] = interpos(vol_mapped(:,ii),jpar0(:,i_lau,ii)/(2*pi*R0));
end
end
@@ -262,35 +285,35 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
rho_tor_norm = ec_data.cd_dens.grids.rho_tor_norm;
% calculate normalised profiles on ec_tgrid_toray grid
- norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_toray);
- norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_toray);
- norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_toray);
- norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+ norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+ norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+ norm_jpar0 = zeros(n_rho,nb_launchers+1,nt_ec_toray);
+ norm_jpar0_integrated = zeros(n_rho,nb_launchers+1,nt_ec_toray);
for it = 1:nt_ec_toray
- norm_p_dens(:,:,it) = p_dens(:,:,it) ./sparse_p_ec_injected(it);
- norm_p_integrated(:,:,it) = p_integrated(:,:,it) ./sparse_p_ec_injected(it);
- norm_cd_dens(:,:,it) = cd_dens(:,:,it) ./sparse_p_ec_injected(it);
- norm_cd_integrated(:,:,it) = cd_integrated(:,:,it)./sparse_p_ec_injected(it);
+ norm_p_dens(:,:,it) = p_dens(:,:,it) ./sparse_p_ec_injected(it);
+ norm_p_integrated(:,:,it) = p_integrated(:,:,it) ./sparse_p_ec_injected(it);
+ norm_jpar0(:,:,it) = jpar0(:,:,it) ./sparse_p_ec_injected(it);
+ norm_jpar0_integrated(:,:,it) = jpar0_integrated(:,:,it)./sparse_p_ec_injected(it);
end
% interpolate normalised p_dens profiles on ec_powers_tgrid
- interp_norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
- interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
- interp_norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
- interp_norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_norm_jpar0 = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_norm_jpar0_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
for i_lau = active_launchers
for irho = 1:n_rho
% get power and current density at each rho
trace_p_dens = squeeze(norm_p_dens(irho,i_lau,:));
trace_p_integrated = squeeze(norm_p_integrated(irho,i_lau,:));
- trace_cd_dens = squeeze(norm_cd_dens(irho,i_lau,:));
- trace_cd_integrated = squeeze(norm_cd_integrated(irho,i_lau,:));
+ trace_jpar0 = squeeze(norm_jpar0(irho,i_lau,:));
+ trace_jpar0_integrated = squeeze(norm_jpar0_integrated(irho,i_lau,:));
% interpolate on gdat powers tgrid, interpos 21 to extrapolate,
- % taking constant value y_edge; needed for time window between start of EC and first TS time
- interp_norm_p_dens( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_dens, ec_tgrid_out(3:end-2));
- interp_norm_p_integrated( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_integrated, ec_tgrid_out(3:end-2));
- interp_norm_cd_dens( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_cd_dens, ec_tgrid_out(3:end-2));
- interp_norm_cd_integrated(irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_cd_integrated,ec_tgrid_out(3:end-2));
+ % taking constant value y_edge; needed for time window between start of EC and first TS time
+ interp_norm_p_dens( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_dens, ec_tgrid_out(3:end-2));
+ interp_norm_p_integrated( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_p_integrated, ec_tgrid_out(3:end-2));
+ interp_norm_jpar0( irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_jpar0, ec_tgrid_out(3:end-2));
+ interp_norm_jpar0_integrated(irho,i_lau,3:end-2) = interpos(21,ec_tgrid_toray,trace_jpar0_integrated,ec_tgrid_out(3:end-2));
end
end
@@ -307,13 +330,13 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
% normalised & interpolated profiles * p_ec_injected on interp_tgrid
interp_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
interp_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
- interp_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_out);
- interp_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_jpar0 = zeros(n_rho,nb_launchers+1,nt_ec_out);
+ interp_jpar0_integrated = zeros(n_rho,nb_launchers+1,nt_ec_out);
for it = 1:nt_ec_out % fill profiles only for time slices when EC is nonzero, thus it+2 assignment
- interp_p_dens(:,:,it) = interp_norm_p_dens(:,:,it) .*p_ec_injected(it);
- interp_p_integrated(:,:,it) = interp_norm_p_integrated(:,:,it) .*p_ec_injected(it);
- interp_cd_dens(:,:,it) = interp_norm_cd_dens(:,:,it) .*p_ec_injected(it);
- interp_cd_integrated(:,:,it) = interp_norm_cd_integrated(:,:,it).*p_ec_injected(it);
+ interp_p_dens(:,:,it) = interp_norm_p_dens(:,:,it) .*p_ec_injected(it);
+ interp_p_integrated(:,:,it) = interp_norm_p_integrated(:,:,it) .*p_ec_injected(it);
+ interp_jpar0(:,:,it) = interp_norm_jpar0(:,:,it) .*p_ec_injected(it);
+ interp_jpar0_integrated(:,:,it) = interp_norm_jpar0_integrated(:,:,it).*p_ec_injected(it);
end
% fill the ids entries
@@ -333,10 +356,10 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
interp_p_integrated(:,i_lau,ii);
% current density
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ...
- interp_cd_dens(:,i_lau,ii);
+ interp_jpar0(:,i_lau,ii);
% integrated current density
ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ...
- interp_cd_integrated(:,i_lau,ii);
+ interp_jpar0_integrated(:,i_lau,ii);
end
% globals
@@ -344,7 +367,7 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.power = ...
interp_p_integrated(end,i_lau,ii);
ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.current_parallel = ...
- interp_cd_integrated(end,i_lau,ii);
+ interp_jpar0_integrated(end,i_lau,ii);
end
end