diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index 3bbed99345c41b3d96245a5977ed06550b367167..5f480dd697fa5a097a5f0201cb1564cd003f98a4 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -205,41 +205,45 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
ids_core_sources.source{last_index+i_lau}.profiles_1d(1:nt_ec_out) = {profiles_template};
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)';
- 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;
+
+ % 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));
- % for testing
- T_mapped(:,ii) = interp1(rho_pol_norm_liu.data,T(:,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) = interp1(rho_pol_norm_liu.data,Rm2_fs(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,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
@@ -249,58 +253,29 @@ 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 = 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));
-
- cd_integrated = ec_data.cd_integrated.data;
- for i_lau = active_launchers
- % convert to j// with conversion factor
- jpar0(:,i_lau,:) = squeeze(cd_dens(:,i_lau,:)).*jtoray_to_jpar0_mapped;
+ % convert and integrate j_TORAY to j//0
+ for i_lau = active_launchers
+% % 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,:)));
- jpar02(:,i_lau,:) = squeeze(jpar_tilde(:,i_lau,:)).*R0.*T_mapped.*Rm2_fs_mapped./B0;
+ 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//
+ % integrate j//0
[~,~,~,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