diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m index 2d9d588f2fdaeda748ce5ad0e1ed0b10cf861cc5..0a4a3963561d525534f70d06fad32f4b9fce50a3 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_core_sources.m +++ b/matlab/TCV_IMAS/tcv_get_ids_core_sources.m @@ -32,106 +32,158 @@ ids_core_sources_description = []; source_template = ids_core_sources.source{1}; profiles_template = source_template.profiles_1d{1}; globals_template = source_template.global_quantities{1}; -%% get time grid and data from gdat & initialize source from template +%% get time grid and data from gdat ohm_gdat = gdat(shot,'ohm_data'); +powers_gdat = gdat(shot,'powers'); % fix time grid to ohm time grid, same as bs t_grid = ohm_gdat.t; n_t = numel(t_grid); t_spacing = mean(diff(t_grid)); +last_index = 0; +%% initialize source from template % ohm id_ohm.description = 'Source from ohmic heating'; id_ohm.index = 7; id_ohm.name = 'ohmic'; -ids_core_sources.source{7} = source_template; -ids_core_sources.source{7}.identifier = id_ohm; -ids_core_sources.source{7}.profiles_1d(1:n_t) = {profiles_template}; -ids_core_sources.source{7}.global_quantities(1:n_t) = {globals_template}; +ids_core_sources.source{1} = source_template; +ids_core_sources.source{1}.identifier = id_ohm; +ids_core_sources.source{1}.profiles_1d(1:n_t) = {profiles_template}; +ids_core_sources.source{1}.global_quantities(1:n_t) = {globals_template}; ohm_data = ohm_gdat.ohm.ohm_data; -% ec -id_ec.description = 'Sources from electron cyclotron heating and current drive'; -id_ec.index = 3; id_ec.name = 'ec'; -ids_core_sources.source{3} = source_template; -ids_core_sources.source{3}.identifier = id_ec; -ids_core_sources.source{3}.profiles_1d(1:n_t) = {profiles_template}; -ids_core_sources.source{3}.global_quantities(1:n_t) = {globals_template}; +%fill profiles for times n t_grid +for ii = 1:n_t + % ohm + % 1d_profiles + ids_core_sources.source{1}.profiles_1d{ii}.time = t_grid(ii); + ids_core_sources.source{1}.profiles_1d{ii}.j_parallel = ... + ohm_data.cd_dens.data(:,ii)'; + % globals + ids_core_sources.source{1}.global_quantities{ii}.time = t_grid(ii); + ids_core_sources.source{1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii); +end +last_index = last_index+1; +%% bs +id_bs.description = 'Bootstrap current'; +id_bs.index = 13; id_bs.name = 'bootstrap_cuurent'; +ids_core_sources.source{last_index+1} = source_template; +ids_core_sources.source{last_index+1}.identifier = id_bs; +ids_core_sources.source{last_index+1}.profiles_1d(1:n_t) = {profiles_template}; +ids_core_sources.source{last_index+1}.global_quantities(1:n_t) = {globals_template}; +% load data +bs_gdat = gdat(shot,'bs_data'); bs_data = bs_gdat.bs.bs_data; + +% fill profiles for times n t_grid +for ii = 1:n_t + ids_core_sources.source{last_index+1}.profiles_1d{ii}.j_parallel = ... + bs_data.cd_dens.data(:,ii)'; +end +last_index = last_index+1; +%% ec % load data -%nb_launchers = size(aa.ec.ec_data.rho_max_pow_dens.data,1); ec_gdat = gdat(shot,'ec_data'); ec_data = ec_gdat.ec.ec_data; ec_time = ec_data.cd_dens.t; + +nb_launchers = size(ec_gdat.ec.ec_data.rho_max_pow_dens.data,1); +id_ec.description = 'Sources from electron cyclotron heating and current drive'; +id_ec.index = 3; id_ec.name = 'ec'; +for i = 1:nb_launchers + ids_core_sources.source{last_index+i} = source_template; + ids_core_sources.source{last_index+i}.identifier = id_ec; + ids_core_sources.source{last_index+i}.profiles_1d(1:n_t) = {profiles_template}; + ids_core_sources.source{last_index+i}.global_quantities(1:n_t) = {globals_template}; +end % interpolate totals on t_grid -ec_total_pow = interp1(ec_data.p_abs_plasma.t,ec_data.p_abs_plasma.data(10,:),t_grid); %the index 10 works for shot 56044 +ec_total_pow = interp1(ec_data.p_abs_plasma.t,ec_data.p_abs_plasma.data(nb_launchers+1,:),t_grid); %the index 10 works for shot 56044 ec_total_pow(isnan(ec_total_pow)) = 0; -ec_total_cur = interp1(ec_data.cd_tot.t,ec_data.cd_tot.data(10,:),t_grid); +ec_total_cur = interp1(ec_data.cd_tot.t,ec_data.cd_tot.data(nb_launchers+1,:),t_grid); ec_total_cur(isnan(ec_total_cur)) = 0; -% bs -id_bs.description = 'Bootstrap current'; -id_bs.index = 13; id_bs.name = 'bootstrap_cuurent'; -ids_core_sources.source{13} = source_template; -ids_core_sources.source{13}.identifier = id_bs; -ids_core_sources.source{13}.profiles_1d(1:n_t) = {profiles_template}; -ids_core_sources.source{13}.global_quantities(1:n_t) = {globals_template}; -% load data -bs_gdat = gdat(shot,'bs_data'); bs_data = bs_gdat.bs.bs_data; +for i = 1:nb_launchers + for ii = 1:n_t + % ec, need to find matching time index\ + ids_core_sources.source{last_index+i}.profiles_1d{ii}.time = t_grid(ii); % profiles time (should profiles time be ec_time?) + ids_core_sources.source{last_index+i}.global_quantities{ii}.time = t_grid(ii); % globals time + [~, ec_idx] = min(abs(ec_time-t_grid(ii))); + if ec_time(ec_idx) == t_grid(ii) % same grid + % current density + ids_core_sources.source{last_index+i}.profiles_1d{ii}.j_parallel = ... + ec_data.cd_dens.data(:,end,ec_idx)'; + % integrated current density + ids_core_sources.source{last_index+i}.profiles_1d{ii}.current_parallel_inside = ... + ec_data.cd_integrated.data(:,ii)'; + + else + ids_core_sources.source{last_index+i}.profiles_1d{ii}.j_parallel = ... + zeros(1,numel(ec_data.cd_dens.x(1,:)))'; + ids_core_sources.source{last_index+i}.profiles_1d{ii}.current_parallel_inside = ... + zeros(1,numel(ec_data.cd_integrated.x(1,:)))'; + end + % globals + ids_core_sources.source{last_index+i}.global_quantities{ii}.time = t_grid(ii); + ids_core_sources.source{last_index+i}.global_quantities{ii}.power = ec_total_pow(ii); + ids_core_sources.source{last_index+i}.global_quantities{ii}.current_parallel = ec_total_cur(ii); + end +end + % nbi id_nbi.description = 'Source from Neutral Beam Injection'; id_nbi.index = 2; id_nbi.name = 'nbi'; -ids_core_sources.source{2} = source_template; -ids_core_sources.source{2}.identifier = id_nbi; -ids_core_sources.source{2}.profiles_1d(1:n_t) = {profiles_template}; -ids_core_sources.source{2}.global_quantities(1:n_t) = {globals_template}; +%ids_core_sources.source{} = source_template; +%ids_core_sources.source{}.identifier = id_nbi; +%ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template}; +%ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template}; % total id_total.description = 'Total source; combines all sources'; id_total.index = 1; id_total.name = 'total'; -ids_core_sources.source{1} = source_template; -ids_core_sources.source{1}.identifier = id_total; -ids_core_sources.source{1}.profiles_1d(1:n_t) = {profiles_template}; -ids_core_sources.source{1}.global_quantities(1:n_t) = {globals_template}; +%ids_core_sources.source{} = source_template; +%ids_core_sources.source{}.identifier = id_total; +%ids_core_sources.source{}.profiles_1d(1:n_t) = {profiles_template}; +%ids_core_sources.source{}.global_quantities(1:n_t) = {globals_template}; + -powers_gdat = gdat(shot,'powers'); %% fill profiles for times n t_grid -for ii = 1:n_t - % ohm - % 1d_profiles - ids_core_sources.source{7}.profiles_1d{ii}.time = t_grid(ii); - ids_core_sources.source{7}.profiles_1d{ii}.j_parallel = ... - ohm_data.cd_dens.data(:,ii)'; - % globals - ids_core_sources.source{7}.global_quantities{ii}.time = t_grid(ii); - ids_core_sources.source{7}.global_quantities{ii}.power = powers_gdat.ohm.data(ii); - - % ec, need to find matching time index\ - ids_core_sources.source{3}.profiles_1d{ii}.time = t_grid(ii); % profiles time (should profiles time be ec_time?) - ids_core_sources.source{3}.global_quantities{ii}.time = t_grid(ii); % globals time - [~, ec_idx] = min(abs(ec_time-t_grid(ii))); - if ec_time(ec_idx) == t_grid(ii) % same grid - % current density - ids_core_sources.source{3}.profiles_1d{ii}.j_parallel = ... - ec_data.cd_dens.data(:,end,ec_idx)'; - % integrated current density - ids_core_sources.source{3}.profiles_1d{ii}.current_parallel_inside = ... - ec_data.cd_integrated.data(:,ii)'; - - else - ids_core_sources.source{3}.profiles_1d{ii}.j_parallel = ... - zeros(1,numel(ec_data.cd_dens.x(1,:)))'; - ids_core_sources.source{3}.profiles_1d{ii}.current_parallel_inside = ... - zeros(1,numel(ec_data.cd_integrated.x(1,:)))'; - end - % globals - ids_core_sources.source{3}.global_quantities{ii}.time = t_grid(ii); - ids_core_sources.source{3}.global_quantities{ii}.power = ec_total_pow(ii); - ids_core_sources.source{3}.global_quantities{ii}.current_parallel = ec_total_cur(ii); - - % bs - ids_core_sources.source{13}.profiles_1d{ii}.j_parallel = ... - bs_data.cd_dens.data(:,ii)'; - -end +% for ii = 1:n_t +% % ohm +% % 1d_profiles +% ids_core_sources.source{1}.profiles_1d{ii}.time = t_grid(ii); +% ids_core_sources.source{1}.profiles_1d{ii}.j_parallel = ... +% ohm_data.cd_dens.data(:,ii)'; +% % globals +% ids_core_sources.source{1}.global_quantities{ii}.time = t_grid(ii); +% ids_core_sources.source{1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii); +% +% % ec, need to find matching time index\ +% ids_core_sources.source{4}.profiles_1d{ii}.time = t_grid(ii); % profiles time (should profiles time be ec_time?) +% ids_core_sources.source{4}.global_quantities{ii}.time = t_grid(ii); % globals time +% [~, ec_idx] = min(abs(ec_time-t_grid(ii))); +% if ec_time(ec_idx) == t_grid(ii) % same grid +% % current density +% ids_core_sources.source{4}.profiles_1d{ii}.j_parallel = ... +% ec_data.cd_dens.data(:,end,ec_idx)'; +% % integrated current density +% ids_core_sources.source{4}.profiles_1d{ii}.current_parallel_inside = ... +% ec_data.cd_integrated.data(:,ii)'; +% +% else +% ids_core_sources.source{4}.profiles_1d{ii}.j_parallel = ... +% zeros(1,numel(ec_data.cd_dens.x(1,:)))'; +% ids_core_sources.source{4}.profiles_1d{ii}.current_parallel_inside = ... +% zeros(1,numel(ec_data.cd_integrated.x(1,:)))'; +% end +% % globals +% ids_core_sources.source{4}.global_quantities{ii}.time = t_grid(ii); +% ids_core_sources.source{4}.global_quantities{ii}.power = ec_total_pow(ii); +% ids_core_sources.source{4}.global_quantities{ii}.current_parallel = ec_total_cur(ii); +% +% % bs +% ids_core_sources.source{2}.profiles_1d{ii}.j_parallel = ... +% bs_data.cd_dens.data(:,ii)'; +% +% end %% add descriptions for profiles_1d