add core sources to tcv imas

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

+function [ids_core_sources,ids_core_sources_description,varargout] = ...
+  tcv_get_ids_core_sources(shot,ids_equil_empty,gdat_params,varargin)
+%
+% [ids_core_sources,ids_core_sources_description,varargout] = ...
+%     tcv_get_ids_core_sources(shot,ids_equil_empty,gdat_params,varargin);
+%
+%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call
+%
+
+machine = 'tcv';
+tens_time = -1;
+tens_rho = -0.1;
+
+if exist('gdat_params','var')
+  [ids_core_sources, params_core_sources] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_sources','gdat_params',gdat_params,varargin{:});
+else
+  [ids_core_sources, params_core_sources] = ...
+    tcv_ids_headpart(shot,ids_equil_empty,'cores_sources',varargin{:});
+  aa=gdat_tcv;
+  gdat_params = aa.gdat_params; % to get default params
+end
+params_eff_ref = gdat_params; params_eff_ref.doplot=0;
+try params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def)
+
+% initialize description
+ids_core_sources_description = [];
+
+%%
+% save source template to fill for differnet source options
+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 form gdat & initialize source from template
+ohm_gdat = gdat(shot,'ohm_data');
+% 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));
+
+% 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};
+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};
+% load data
+ec_gdat = gdat(shot,'ec_data'); ec_data = ec_gdat.ec.ec_data;
+ec_time = ec_gdat.ec.ec_data.cd_dens.t;
+% interpolate totals on t_grid
+ec_total_pow = interp1(ec_data.p_abs_plasma.t,ec_data.p_abs_plasma.data(12,:),t_grid);
+ec_total_pow(isnan(ec_total_pow)) = 0;
+ec_total_cur = interp1(ec_data.p_abs_plasma.t,ec_data.cd_tot.data(12,:),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;
+
+% nbi
+ids_core_sources.source{2} = source_template;
+% to be filled
+
+% total
+ids_core_sources.source{1} = source_template;
+% to be filled
+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
+  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
+
+%% add descriptions for profiles_1d
+
+desc.source = 'available by now, 1:total, 3:ec, 7:ohm, 13:bootstrap';
+desc.profiles_1d.j_parallel = 'parallel current density';
+desc.profiles_1d.current_parallel_inside = 'integrated parallel current density';
+desc.globals.power = 'total power coupled to the plasma';
+desc.globals.current_parallel = 'total parallel current driven';
+ids_core_sources_description = desc;
+
+%%
+if nargin <= 2
+  ids_core_sources.code.name = ['tcv_get_ids_core_sources, within gdat, with shot= ' num2str(params_core_sources.shot) ];
+else
+  ids_core_sources.code.name = ['tcv_get_ids_core_sources, within gdat, with shot= ' num2str(params_core_sources.shot) '; varargin: ' varargin{:}];
+end
+ids_core_sources_description.code.name = ids_core_sources.code.name;
+
+ids_core_sources.code.output_flag = zeros(size(ids_core_sources.time));
+
+% cocos automatic transform
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
+  [ids_core_sources,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_core_sources,'core_sources',gdat_params.cocos_in, ...
+          gdat_params.cocos_out,gdat_params.ipsign_out,gdat_params.b0sign_out,gdat_params.ipsign_in,gdat_params.b0sign_in, ...
+          gdat_params.error_bar,gdat_params.nverbose);
+
+end
\ No newline at end of file