diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
index a025571dda5ea5668f80618d6d3c39f9c29804d2..4104ea8a1e223eed7447e09aaff562f87f0149f6 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m
@@ -49,31 +49,40 @@ 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);
-%% load liuqe data from nodes for conversions
-params_eff = params_eff_ref;
-% normalized sqrt poloidal flux including axis
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:rho';
-rho_pol_norm_liu = gdat(shot,params_eff);
-% <1/R^2>
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
-Rm2_fs = gdat(shot,params_eff);
-% Volume
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:vol';
-vol= gdat(shot,params_eff);
-% T(\psi) = R*B_tor
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:rbtor_rho';
-T = gdat(shot,params_eff);
-% Q1Q: -dpsi/dV
-%mVprime = 1./LY.Q1Q; % Q1Q: -dpsi/dV
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:DPSI_VOL';
-Vprime = gdat(shot,params_eff);
-% I_tor = I_p, total plasma current
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:i_pl';
-Ip = gdat(shot,params_eff);
-
+%% load liuqe data from nodes for conversions
+if params_eff_ref.liuqe == 1
+ nodes_liuqe = 'equil';
+else
+ nodes_liuqe = ['equil',num2str(params_eff_ref.liuqe)];
+end
+if check_nodes_filled(params_core_sources.shot,nodes_liuqe)
+ params_eff = params_eff_ref;
+ % normalized sqrt poloidal flux including axis
+ params_eff.data_request='\tcv_shot::top.results.equil_1.results:rho';
+ rho_pol_norm_liu = gdat(params_core_sources.shot,params_eff);
+ % <1/R^2>
+ params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
+ Rm2_fs_liu = gdat(params_core_sources.shot,params_eff);
+ % Volume
+ params_eff.data_request='\tcv_shot::top.results.equil_1.results:vol';
+ vol_liu= gdat(params_core_sources.shot,params_eff);
+ % T(\psi) = R*B_tor
+ params_eff.data_request='\tcv_shot::top.results.equil_1.results:rbtor_rho';
+ T_liu = gdat(params_core_sources.shot,params_eff);
+ % safety factor q = -dPhi/dPsi
+ params_eff.data_request='\tcv_shot::top.results.equil_1.results:q';
+ q_liu = gdat(params_core_sources.shot,params_eff);
+ % I_tor = I_p, total plasma current
+ params_eff.data_request='\tcv_shot::top.results.equil_1.results:i_pl';
+ Ip_liu = gdat(params_core_sources.shot,params_eff);
+ % define liuqe_time from Ip_liu
+ liuqe_time = Ip_liu.t;
+else
+ error('Liuqe nodes %s not filled. Contact O. Sauter.')
+end
%% initialize source from template
% ohm
-ohm_t_grid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_t_grid);
+ohm_tgrid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_tgrid);
main_desc = 'Source from ohmic heating'; production = 'IBS nodes';
id_ohm.description = sprintf('%s from %s',main_desc,production);
@@ -86,35 +95,24 @@ ohm_data = ohm_gdat.ohm.ohm_data;
% 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);
-params_eff = params_eff_ref;
-% normalized sqrt poloidal flux including axis
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:rho';
-rho_pol_norm_liu = gdat(shot,params_eff);
-% <1/R^2>
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
-Rm2_fs = gdat(shot,params_eff);
-% Volume
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:vol';
-vol= gdat(shot,params_eff);
-% T(\psi) = R*B_tor
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:rbtor_rho';
-T = gdat(shot,params_eff);
-%rho_pol_norm_liu = L.pQ; % normalized sqrt poloidal flux including axis
-%T = LY.TQ;
-%Rm2_fs = LY.Q2Q;
-%vol = LY.VQ;
+% interpolate liuqe outputs on ohm_tgrid
+T = interp1(liuqe_time,T_liu.data.', ohm_tgrid)';
+Rm2_fs = interp1(liuqe_time,Rm2_fs_liu.data.',ohm_tgrid)';
+vol = interp1(liuqe_time,vol_liu.data.', ohm_tgrid)';
+% 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,ohm_t_grid);
-jpar_tilde_to_jpar_0 = R0*T.data.*Rm2_fs.data./B0';
+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';
% 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');
+jpar_tilde_to_jpar_0_mapped = ...
+ interp1(rho_pol_norm_liu.data,jpar_tilde_to_jpar_0,ohm_data.cd_dens.rhopol_norm);
+vol_mapped = interp1(rho_pol_norm_liu.data,vol,ohm_data.cd_dens.rhopol_norm);
for ii = 1:ohm_n_t
% profiles_1d
- 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}.time = ohm_tgrid(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 = ...
@@ -125,7 +123,7 @@ for ii = 1:ohm_n_t
ids_core_sources.source{last_index+1}.profiles_1d{ii}.current_parallel_inside = integrated_jpar_0_tmp;
% globals
- ids_core_sources.source{last_index+1}.global_quantities{ii}.time = ohm_t_grid(ii);
+ ids_core_sources.source{last_index+1}.global_quantities{ii}.time = ohm_tgrid(ii);
ids_core_sources.source{last_index+1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii);
ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
end
@@ -135,7 +133,7 @@ last_index = last_index+1; % add if statement to only increment if ohmic source
%% bs
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);
+bs_tgrid = bs_gdat.bs.t; bs_n_t = numel(bs_tgrid);
main_desc = 'Bootstrap current'; production = 'IBS nodes';
id_bs.description = sprintf('%s from %s',main_desc,production);
@@ -147,7 +145,7 @@ ids_core_sources.source{last_index+1}.global_quantities(1:bs_n_t) = {globals_tem
for ii = 1:bs_n_t
% profiles_1d
- 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}.time = bs_tgrid(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 = ...
@@ -159,7 +157,7 @@ for ii = 1:bs_n_t
ids_core_sources.source{last_index+1}.profiles_1d{ii}.current_parallel_inside = integrated_jpar_0_tmp;
% globals
- ids_core_sources.source{last_index+1}.global_quantities{ii}.time = bs_t_grid(ii);
+ ids_core_sources.source{last_index+1}.global_quantities{ii}.time = bs_tgrid(ii);
ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
end
@@ -187,11 +185,11 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
% find times where EC is on to define time grid with extra time slice just
% before/after EC power and at start/end of shot
itime_ec = find(powers_gdat.ec.data(:,end)>0);
- if ec_powers_tgrid(itime_ec(end))>=ohm_t_grid(end)
- i_time_end = iround(ec_powers_tgrid,ohm_t_grid(end));
- ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:i_time_end)'];
+ if ec_powers_tgrid(itime_ec(end))>=ohm_tgrid(end)
+ i_time_end = iround(ec_powers_tgrid,ohm_tgrid(end));
+ ec_tgrid_out = [ohm_tgrid(1),ec_powers_tgrid(itime_ec(1)-1:i_time_end)'];
else
- ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_t_grid(end)];
+ ec_tgrid_out = [ohm_tgrid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_tgrid(end)];
end
nt_ec_out = numel(ec_tgrid_out);
p_ec_injected = interpos(ec_powers_tgrid,powers_gdat.ec.data(:,end),ec_tgrid_out);
@@ -210,23 +208,21 @@ 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> / (R_0 <B_\phi/R>^2)
- [L,~,LY] = liuqe(shot,ec_tgrid_toray,'iterq',50,'ilim',3,'icsint',true);
- rho_pol_norm_liu = L.pQ; % normalized sqrt poloidal flux including axis
+ % conversion using j_//,cd = <j_cd/B> B, thus j_// = 2\pi <B^2> / (B_0 <B_\phi/R>)
mu0 = 4*pi*10^-7;
- Rm2_fs = LY.Q2Q;
- params_eff = params_eff_ref;
- params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
- Rm2_fs = gdat(shot,params_eff);
+ % 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;
- mVprime = 1./LY.Q1Q; % Q1Q: -dpsi/dV
- Ip_psi = LY.ItQ; % Ip(psi) = I_tor(psi)
- T = LY.TQ;
- vol = LY.VQ;
- B2_fs = mu0*Ip_psi./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;
% interpolate on ec_data rho_pol grid
@@ -234,9 +230,9 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
vol_mapped = nan(size(ec_data.cd_dens.x));
for ii = 1:nt_ec_toray
jtoray_to_jpar0_mapped(:,ii) = ...
- interp1(rho_pol_norm_liu,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,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,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+ interp1(rho_pol_norm_liu.data,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
end
% load 1d_profiles from ec_data
@@ -373,11 +369,11 @@ if numel(active_nbi)>0
% find times where NBI is on to define time grid with extra time slice just
% before/after NBI power and at start/end of shot
itime_nbi = find((powers_gdat.(nbi_names{active_nbi(1)}).data>0));
- if nbi_powers_tgrid(itime_nbi(end))>=ohm_t_grid(end)
- i_time_end = iround(nbi_powers_tgrid,ohm_t_grid(end));
- nbi_tgrid_out = [ohm_t_grid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
+ if nbi_powers_tgrid(itime_nbi(end))>=ohm_tgrid(end)
+ i_time_end = iround(nbi_powers_tgrid,ohm_tgrid(end));
+ nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
else
- nbi_tgrid_out = [ohm_t_grid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_t_grid(end)];
+ nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_tgrid(end)];
end
nt_nbi_out = numel(nbi_tgrid_out);
@@ -425,11 +421,11 @@ if ~isempty(powers_gdat.dnbi)
% find times where DNBI is on to define time grid with extra time slice just
% before/after DNBI power and at start/end of shot
itime_dnbi = find((powers_gdat.dnbi.data>0));
- if dnbi_powers_tgrid(itime_dnbi(end))>=ohm_t_grid(end)
- i_time_end = iround(dnbi_powers_tgrid,ohm_t_grid(end));
- dnbi_tgrid_out = [ohm_t_grid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:i_time_end)'];
+ if dnbi_powers_tgrid(itime_dnbi(end))>=ohm_tgrid(end)
+ i_time_end = iround(dnbi_powers_tgrid,ohm_tgrid(end));
+ dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:i_time_end)'];
else
- dnbi_tgrid_out = [ohm_t_grid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:itime_dnbi(end)+1)',ohm_t_grid(end)];
+ dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:itime_dnbi(end)+1)',ohm_tgrid(end)];
end
nt_dnbi_out = numel(dnbi_tgrid_out);