diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index 00eb0980395422472d9424f9c633fd957413a432..3bbed99345c41b3d96245a5977ed06550b367167 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -208,7 +208,7 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
% 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>)
+ % 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)';
@@ -216,14 +216,17 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
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;
+ Vprime = -2*pi*q./T./Rm2_fs;
+ B2_fs = -mu0*Ip./Vprime + T.^2.*Rm2_fs;
+ % Iphi? needs to be a profile?
+ % R is also a profile in icdbseval? should not be! It's definitely R0 see
+ % eq.(30) in dok
% 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;
+ jtoray_to_jpar0 = 2*pi./B0./T./Rm2_fs.*B2_fs;
% interpolate on ec_data rho_pol grid
jtoray_to_jpar0_mapped = nan(size(ec_data.cd_dens.x));
@@ -233,6 +236,10 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
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));
+ % for testing
+ T_mapped(:,ii) = interp1(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) = interp1(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,17 +249,58 @@ 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;
+
+ % initialize variables
+ jpar0 = zeros(size(cd_dens));
+ jpar0_integrated = zeros(size(cd_dens));
+
+ cd_integrated = ec_data.cd_integrated.data;
- cd_integrated = zeros(size(ec_data.cd_integrated.data));
for i_lau = active_launchers
% convert to j// with conversion factor
- cd_dens(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
+ 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,:)));
+ jpar02(:,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);
+ [~,~,~,jpar0_integrated(:,i_lau,ii)] = interpos(vol_mapped(:,ii),jpar0(:,i_lau,ii)/(2*pi*R0));
+ [~,~,~,jpar02_integrated(:,i_lau,ii)] = interpos(vol_mapped(:,ii),jpar02(:,i_lau,ii)/(2*pi*R0));
end
end
-
+
+ figure; hold on; plot(squeeze(jpar0(:,1,50)));
+ plot(squeeze(jpar02(:,1,50)),'x');
+ legend('ids','icdbseval'); ylabel('j||0')
+
+ figure; plot(squeeze(cd_integrated(end,1,:))); hold on;
+ plot(squeeze(jpar0_integrated(end,1,:)));
+ plot(squeeze(jpar02_integrated(end,1,:)));
+ legend('I_\phi','I||,ids','I||,icdbseval');
+
+ % test to integrate but still is not the same
+ ilau = 4;
+ jtilde_test = B0./R0./T_mapped./Rm2_fs_mapped.*squeeze(jpar0(:,ilau,:));
+ jtilde_test2 = squeeze(jpar_tilde(:,ilau,:));
+ for ii = 1:nt_ec_toray
+ [~,~,~,icd_test(:,ii)] = interpos(vol_mapped(:,ii),Rm2_fs_mapped(:,ii)./T_mapped(:,ii).*jtilde_test(:,ii));
+ out(:,ii) = R0/2/pi.*T_mapped(:,ii).*icd_test(:,ii);
+ % icdbseval
+ [~,~,~,icd_test2(:,ii)] = interpos(vol_mapped(:,ii),Rm2_fs_mapped(:,ii)./T_mapped(:,ii).*jtilde_test2(:,ii));
+ out2(:,ii) = R0/2/pi.*T_mapped(:,ii).*icd_test2(:,ii);
+ end
+
+ figure;
+ plot(ec_tgrid_toray,out(end,:)); hold on;
+ plot(ec_tgrid_toray,squeeze(cd_integrated(end,ilau,:)));
+ plot(ec_tgrid_toray,out2(end,:),'x');
+ ylabel('I_\phi'); legend('from ids','from icdbseval','TORAY')
+
+ figure;
+ plot(squeeze(out(end,:))-squeeze(cd_integrated(end,1,:))'); hold on; plot(squeeze(out2(end,:))-squeeze(cd_integrated(end,1,:))');
+ ylabel('diff (I_\phi,I_TORARY)'); legend('from ids','from icdbseval')
+
%% interpolating 1d_profiles from 'ec_tgrid_toray' grid (toray tgrid) to 'interp_tgrid' (powers tgrid when powers.data>0)
ij = iround_os(ec_tgrid_out,ec_tgrid_toray);
sparse_p_ec_injected = p_ec_injected(ij); % injected ec power vals corresponding to ec_tgrid_toray
@@ -262,35 +310,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));
+ 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 +355,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 +381,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 +392,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