diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m index 83bdd9a85a73f77d6d29834ddabf4582c179d41a..a5d65d69ea9098d35f72817593d516c9a2c29d4f 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m +++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m @@ -87,129 +87,132 @@ end last_index = last_index+1; % add if statement to only increment if bs source has been added %% ec -% load data -ec_gdat = gdat(shot,'ec_data'); -ec_data = ec_gdat.ec.ec_data; ec_inputs = ec_gdat.ec.ec_inputs; -% get tgrid from gdat ec_data -ec_data_tgrid = ec_data.p_dens.t; nt_ec_data = numel(ec_data_tgrid); - -% get tgrid from gdat powers -ec_powers_tgrid = powers_gdat.ec.t; nt_ec_powers = numel(ec_powers_tgrid); - -% retrieve active launcher information from ec_inputs -nb_launchers = numel(ec_inputs.launchers_active.data); -active_launchers = find(ec_inputs.launchers_active.data==1)'; -n_active_launchers = sum(ec_inputs.launchers_active.data); - -% Setup structures for active launchers from template -main_desc = 'Source from electron cyclotron heating and current drive'; -id_ec.index = 3; id_ec.name = 'ec'; -for i_lau = active_launchers - id_ec.description = sprintf('L%i/G%i, %s',i_lau,ec_inputs.gyro2launcher.data(i_lau),main_desc); - ids_core_sources.source{last_index+i_lau} = source_template; - ids_core_sources.source{last_index+i_lau}.identifier = id_ec; - ids_core_sources.source{last_index+i_lau}.profiles_1d(1:nt_ec_powers) = {profiles_template}; - ids_core_sources.source{last_index+i_lau}.global_quantities(1:nt_ec_powers) = {globals_template}; -end - -% get data for globals from gdat powers and fill in ids_structure -ec_total_pow = transpose(powers_gdat.ec.data(:,nb_launchers+1)); %use power from powers_gdat(injected ec power) instead of ec_data power -ec_total_pow(isnan(ec_total_pow)) = 0; -ec_total_cur = ec_data.cd_tot.data(nb_launchers+1,:); -ec_total_cur(isnan(ec_total_cur)) = 0; - -for i_lau = active_launchers - for ii = 1:nt_ec_powers - ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.time = ec_powers_tgrid(ii); - ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.power = ec_total_pow(ii); - ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.current_parallel = ec_total_cur(ii); +ec_gdat = gdat(shot,'ec_data'); + +if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources + + ec_data = ec_gdat.ec.ec_data; ec_inputs = ec_gdat.ec.ec_inputs; + % get tgrid from gdat ec_data + ec_data_tgrid = ec_data.p_dens.t; nt_ec_data = numel(ec_data_tgrid); + + % get tgrid from gdat powers + ec_powers_tgrid = powers_gdat.ec.t; nt_ec_powers = numel(ec_powers_tgrid); + + % retrieve active launcher information from ec_inputs + nb_launchers = numel(ec_inputs.launchers_active.data); + active_launchers = find(ec_inputs.launchers_active.data==1)'; + + % Setup structures for active launchers from template + main_desc = 'Source from electron cyclotron heating and current drive'; + id_ec.index = 3; id_ec.name = 'ec'; + for i_lau = active_launchers + id_ec.description = sprintf('L%i/G%i, %s',i_lau,ec_inputs.gyro2launcher.data(i_lau),main_desc); + ids_core_sources.source{last_index+i_lau} = source_template; + ids_core_sources.source{last_index+i_lau}.identifier = id_ec; + ids_core_sources.source{last_index+i_lau}.profiles_1d(1:nt_ec_powers) = {profiles_template}; + ids_core_sources.source{last_index+i_lau}.global_quantities(1:nt_ec_powers) = {globals_template}; end -end -% interpoating p_dens profiles from 'ec_data_tgrid' grid (toray tgrid) to 'ec_powers_tgrid' (powers tgrid) -p_dens = ec_data.p_dens.data; -p_integrated = ec_data.p_integrated.data; -cd_dens = ec_data.cd_dens.data; -cd_integrated = ec_data.cd_integrated.data; - -p_ec_injected = powers_gdat.ec.data; -ij = iround_os(ec_powers_tgrid,ec_data_tgrid); -sparse_p_ec_injected = p_ec_injected(ij,:); % injected ec power vals corresponding to ec_data_tgrid -n_rho = size(p_dens, 1); - -% calculate normalised profiles on ec_data_tgrid grid -norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_data); -norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_data); -norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_data); -norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_data); -for it = 1:nt_ec_data - 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,:); -end + % get data for globals from gdat powers and fill in ids_structure + ec_total_pow = transpose(powers_gdat.ec.data(:,nb_launchers+1)); %use power from powers_gdat(injected ec power) instead of ec_data power + ec_total_pow(isnan(ec_total_pow)) = 0; + ec_total_cur = ec_data.cd_tot.data(nb_launchers+1,:); + ec_total_cur(isnan(ec_total_cur)) = 0; + + for i_lau = active_launchers + for ii = 1:nt_ec_powers + ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.time = ec_powers_tgrid(ii); + ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.power = ec_total_pow(ii); + ids_core_sources.source{last_index+i_lau}.global_quantities{ii}.current_parallel = ec_total_cur(ii); + end + end -% interpolate normalised p_dens profiles on ec_powers_tgrid -interp_norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); -interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); -interp_norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); -interp_norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); -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,:)); - % interpolate on gdat powers tgrid - interp_norm_p_dens(irho,i_lau,:) = interp1(ec_data_tgrid,trace_p_dens,ec_powers_tgrid); - interp_norm_p_integrated(irho,i_lau,:) = interp1(ec_data_tgrid,trace_p_integrated,ec_powers_tgrid); - interp_norm_cd_dens(irho,i_lau,:) = interp1(ec_data_tgrid,trace_cd_dens,ec_powers_tgrid); - interp_norm_cd_integrated(irho,i_lau,:) = interp1(ec_data_tgrid,trace_cd_integrated,ec_powers_tgrid); + % interpoating p_dens profiles from 'ec_data_tgrid' grid (toray tgrid) to 'ec_powers_tgrid' (powers tgrid) + p_dens = ec_data.p_dens.data; + p_integrated = ec_data.p_integrated.data; + cd_dens = ec_data.cd_dens.data; + cd_integrated = ec_data.cd_integrated.data; + + p_ec_injected = powers_gdat.ec.data; + ij = iround_os(ec_powers_tgrid,ec_data_tgrid); + sparse_p_ec_injected = p_ec_injected(ij,:); % injected ec power vals corresponding to ec_data_tgrid + n_rho = size(p_dens, 1); + + % calculate normalised profiles on ec_data_tgrid grid + norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_data); + norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_data); + norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_data); + norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_data); + for it = 1:nt_ec_data + 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,:); end -end -% normalised & interpolated profiles * p_ec_injected on ec_powers_tgrid -interp_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); -interp_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); -interp_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); -interp_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); -for it = 1:nt_ec_powers - 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,:); -end + % interpolate normalised p_dens profiles on ec_powers_tgrid + interp_norm_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); + interp_norm_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); + interp_norm_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); + interp_norm_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); + 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,:)); + % interpolate on gdat powers tgrid + interp_norm_p_dens(irho,i_lau,:) = interp1(ec_data_tgrid,trace_p_dens,ec_powers_tgrid); + interp_norm_p_integrated(irho,i_lau,:) = interp1(ec_data_tgrid,trace_p_integrated,ec_powers_tgrid); + interp_norm_cd_dens(irho,i_lau,:) = interp1(ec_data_tgrid,trace_cd_dens,ec_powers_tgrid); + interp_norm_cd_integrated(irho,i_lau,:) = interp1(ec_data_tgrid,trace_cd_integrated,ec_powers_tgrid); + end + end -% fill profiles_1d with interpolated profiles -for i_lau = active_launchers - for ii = 1:nt_ec_powers - % time - ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = ec_powers_tgrid(ii); - % rhotor grid - ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = ... - ec_data.p_dens.rhotor_norm(1,:); % constant in time - % power density - ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.energy = ... - interp_p_dens(:,i_lau,ii); - % integrated power density - ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.power_inside = ... - interp_p_integrated(:,i_lau,ii); - % current density (to adapt to <J.B>/B0) - ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ... - interp_cd_dens(:,i_lau,ii); - % integrated current density (to adapt to INTEGRAL(<J.B>/B0)*ds_phi) - ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ... - interp_cd_integrated(:,i_lau,ii); + % normalised & interpolated profiles * p_ec_injected on ec_powers_tgrid + interp_p_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); + interp_p_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); + interp_cd_dens = zeros(n_rho,nb_launchers+1,nt_ec_powers); + interp_cd_integrated = zeros(n_rho,nb_launchers+1,nt_ec_powers); + for it = 1:nt_ec_powers + 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,:); + end + + % fill profiles_1d with interpolated profiles + for i_lau = active_launchers + for ii = 1:nt_ec_powers + % time + ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.time = ec_powers_tgrid(ii); + % rhotor grid + ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.grid.rho_tor_norm = ... + ec_data.p_dens.rhotor_norm(1,:); % constant in time + % power density + ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.energy = ... + interp_p_dens(:,i_lau,ii); + % integrated power density + ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.electrons.power_inside = ... + interp_p_integrated(:,i_lau,ii); + % current density (to adapt to <J.B>/B0) + ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.j_parallel = ... + interp_cd_dens(:,i_lau,ii); + % integrated current density (to adapt to INTEGRAL(<J.B>/B0)*ds_phi) + ids_core_sources.source{last_index+i_lau}.profiles_1d{ii}.current_parallel_inside = ... + interp_cd_integrated(:,i_lau,ii); + end end -end -%add empty sources for rest of unsused launchers -if numel(ids_core_sources.source)-last_index ~= nb_launchers - ids_core_sources.source{last_index+nb_launchers} = []; -end - -last_index = last_index+nb_launchers; + %add empty sources for rest of unsused launchers + if numel(ids_core_sources.source)-last_index ~= nb_launchers + ids_core_sources.source{last_index+nb_launchers} = []; + end + + last_index = last_index+nb_launchers; + +end %% nbi id_nbi.description = 'Source from Neutral Beam Injection';