Skip to content
Snippets Groups Projects
Commit 6f3f31d3 authored by Olivier Sauter's avatar Olivier Sauter
Browse files

start add core_profiles ids 

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11789 d63d8f72-b253-0410-a779-e742ad2e26cf
parent bff7f955
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
crpptbx/TCV_IMAS/tcv_get_ids_core_profiles.m 0 → 100644
+ 299
− 0
View file @ 6f3f31d3
function [ids_cores_profiles,ids_cores_profiles_description,varargout] = tcv_get_ids_core_profiles(shot,ids_equil_empty,varargin);
%
% [ids_cores_profiles,ids_cores_profiles_description,varargout] = tcv_get_ids_cores_profiles(shot,ids_equil_empty,varargin);
%
%
machine = 'tcv';
tens_time = -1;
tens_rho = -0.1;
[ids_cores_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',varargin{:});
ids_cores_profiles_description = [];
% base all from times fro nete_rho (which should be conf by default)
temp_1d.ne_rho = gdat(params_cores_profiles.shot,'ne_rho','machine',machine,'fit',1);
temp_1d.te_rho = gdat(params_cores_profiles.shot,'te_rho','machine',machine,'fit',1);
% compute grids_1d for fit rhopol grid, need a gdat_data like structure
temp_1d.fit.te_rho = temp_1d.te_rho.fit;
temp_1d.fit.te_rho.gdat_params = temp_1d.te_rho.gdat_params;
temp_1d.fit.te_rho.shot = temp_1d.te_rho.shot;
temp_1d.fit.te_rho = get_grids_1d(temp_1d.fit.te_rho,1,1);
temp_1d.fit.ne_rho = temp_1d.ne_rho.fit;
ids_cores_profiles.time = temp_1d.te_rho.fit.t;
ids_cores_profiles_description.time = ['.t subfield from: [''te_rho'',''subnode'',''fit.t''] thus from ' temp_1d.te_rho.fit.data_fullpath];
% make empty cell arrays for subnodes not filled in before copying default structure to all times
ids_cores_profiles.profiles_1d{1}.ion{1}.state = {};
ids_cores_profiles.profiles_1d{1}.neutral{1}.state = {};
ids_cores_profiles.profiles_1d(1:length(ids_cores_profiles.time)) = ids_cores_profiles.profiles_1d(1);
% As a general rule, for a new substructure under the main ids, construct a local structure like:
% "global_quantities" with subfields being the relevant data to get and a local structure:
% "global_quantities_desc" which contains the same subfields themselves containing the gdat string aftre shot used
%
% vacuum_toroidal_field and time, using homogeneous
%
vacuum_toroidal_field.b0=gdat(params_cores_profiles.shot,'b0','source','liuqe','machine',machine); % to get on liuqe time array
vacuum_toroidal_field_desc.b0 = '''b0'',''source'',''liuqe''';
vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'',''source'',''liuqe'']';
ids_cores_profiles.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
ids_cores_profiles.vacuum_toroidal_field.b0 = interpos(63,vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data,ids_cores_profiles.time,-1);
ids_cores_profiles_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
% global_quantities data into local global_quantities.* structure with correct end names and global_quantities_desc.* with description. Use temp.* and temp_desc.* structures for temporary data
global_quantities.ip = gdat(params_cores_profiles.shot,'ip','machine',machine);
global_quantities_desc.ip = 'ip';
current_non_inductive = gdat(params_cores_profiles.shot,'ohm_data','machine',machine);
global_quantities.current_non_inductive.data = current_non_inductive.ohm.ohm_data.cd_tot.data;
global_quantities.current_non_inductive.t = current_non_inductive.ohm.ohm_data.cd_tot.t;
global_quantities_desc.current_non_inductive = current_non_inductive.ohm.ohm_data.cd_tot.label;
current_bootstrap = gdat(params_cores_profiles.shot,'bs_data','machine',machine);
global_quantities.current_bootstrap.data = current_bootstrap.bs.bs_data.cd_tot.data;
global_quantities.current_bootstrap.t = current_bootstrap.bs.bs_data.cd_tot.t;
global_quantities_desc.current_bootstrap = current_bootstrap.bs.bs_data.cd_tot.label;
global_quantities.v_loop = gdat(params_cores_profiles.shot,'vloop','machine',machine);
global_quantities_desc.v_loop = 'vloop';
global_quantities.li_3 = gdat(params_cores_profiles.shot,'li','machine',machine);
global_quantities_desc.li_3 = 'li';
global_quantities.beta_tor = gdat(params_cores_profiles.shot,'beta','machine',machine);
global_quantities_desc.beta_tor = 'beta';
global_quantities.beta_tor_norm = gdat(params_cores_profiles.shot,'betan','machine',machine);
global_quantities_desc.beta_tor_norm = 'betan';
global_quantities.beta_pol = gdat(params_cores_profiles.shot,'betap','machine',machine);
global_quantities_desc.beta_pol = 'betap';
global_quantities.energy_diamagnetic = gdat(params_cores_profiles.shot,'w_mhd','machine',machine);
global_quantities_desc.energy_diamagnetic = 'w_mhd';
global_quantities_fieldnames = fieldnames(global_quantities);
special_fields = {''}; % fields needing non-automatic treatments
for i=1:length(global_quantities_fieldnames)
if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
if ~isstruct(ids_cores_profiles.global_quantities.(global_quantities_fieldnames{i}))
ids_cores_profiles.global_quantities.(global_quantities_fieldnames{i}) = ...
interpos(global_quantities.(global_quantities_fieldnames{i}).t, ...
global_quantities.(global_quantities_fieldnames{i}).data,ids_cores_profiles.time,tens_time);
else
special_fields{end+1} = global_quantities_fieldnames{i};
end
end
end
% special case
%
%% profiles_1d (cannot use eqdsk since not same radial mesh)
temp_1d.area = gdat(shot,'area_rho','machine',machine);
for ir=1:length(temp_1d.area.x)
aa(ir,:) = interpos(temp_1d.area.t,temp_1d.area.data(ir,:),ids_cores_profiles.time,tens_time);
end
it_thom = iround_os(temp_1d.te_rho.t,ids_cores_profiles.time);
for it=1:length(ids_cores_profiles.time)
ids_cores_profiles.profiles_1d{it}.grid.rho_tor_norm = temp_1d.fit.te_rho.grids_1d.rhotornorm(:,it);
ids_cores_profiles.profiles_1d{it}.grid.rho_tor = temp_1d.fit.te_rho.grids_1d.rhotornorm(:,it) ...
.* temp_1d.fit.te_rho.grids_1d.rhotor_edge(it);
ids_cores_profiles.profiles_1d{it}.grid.psi = temp_1d.fit.te_rho.grids_1d.psi(:,it);
ids_cores_profiles.profiles_1d{it}.grid.volume = temp_1d.fit.te_rho.grids_1d.rhovolnorm(:,it).^2 ...
.* temp_1d.fit.te_rho.grids_1d.volume_edge(it);
ids_cores_profiles.profiles_1d{it}.grid.area = interpos(temp_1d.area.x,aa(:,it),temp_1d.fit.te_rho.grids_1d.rhopolnorm, ...
tens_rho,[1 2],[0 aa(end,it)]);
ids_cores_profiles.profiles_1d{it}.time = ids_cores_profiles.time(it);
ids_cores_profiles.profiles_1d{it}.electrons.temperature = temp_1d.fit.te_rho.data(:,it);
ids_cores_profiles.profiles_1d{it}.electrons.temperature_fit.measured = temp_1d.te_rho.data(:,it_thom(it));
ids_cores_profiles.profiles_1d{it}.electrons.temperature_fit.measured_error_upper = temp_1d.te_rho.error_bar(:,it_thom(it));
ids_cores_profiles.profiles_1d{it}.electrons.temperature_fit.time_measurement = temp_1d.te_rho.t(it_thom(it));
ids_cores_profiles.profiles_1d{it}.electrons.temperature_fit.rho_tor_norm = temp_1d.te_rho.grids_1d.rhotornorm(:,it_thom(it));
ids_cores_profiles.profiles_1d{it}.electrons.temperature_fit.source = 'Thomson, interpos fit';
ids_cores_profiles.profiles_1d{it}.electrons.density = temp_1d.fit.ne_rho.data(:,it);
ids_cores_profiles.profiles_1d{it}.electrons.density_fit.measured = temp_1d.ne_rho.data(:,it_thom(it));
ids_cores_profiles.profiles_1d{it}.electrons.density_fit.measured_error_upper = temp_1d.ne_rho.error_bar(:,it_thom(it));
ids_cores_profiles.profiles_1d{it}.electrons.density_fit.time_measurement = temp_1d.ne_rho.t(it_thom(it));
ids_cores_profiles.profiles_1d{it}.electrons.density_fit.rho_tor_norm = temp_1d.ne_rho.grids_1d.rhotornorm(:,it_thom(it));
ids_cores_profiles.profiles_1d{it}.electrons.density_fit.source = 'Thomson, interpos fit';
ids_cores_profiles.profiles_1d{it}.electrons.pressure_thermal = 1.6022e-19.*ids_cores_profiles.profiles_1d{it}.electrons.density ...
.* ids_cores_profiles.profiles_1d{it}.electrons.temperature;
end
% ion, assume only D if no CXRS (need to ask how to check if H...)
temp_1d.cxrs_rho = gdat(params_cores_profiles.shot,'cxrs','machine',machine,'fit',1);
temp_1d.ti_conf_rho = gdat(params_cores_profiles.shot,'results.conf:ti','machine',machine,'fit',1);
temp_1d.ni_conf_rho = gdat(params_cores_profiles.shot,'results.conf:ni','machine',machine,'fit',1);
temp_1d.z_eff_conf = gdat(params_cores_profiles.shot,'results.conf:z_eff','machine',machine,'fit',1);
data_fullpath_raw = 'Ti(C sometimes B) from cxrs system 1 to 3';
temp_1d.ti.raw = temp_1d.cxrs_rho.ti.raw;
temp_1d.ti.raw.shot = temp_1d.cxrs_rho.shot;temp_1d.ti.raw.gdat_params = temp_1d.cxrs_rho.gdat_params;
temp_1d.ti.raw.x =temp_1d.cxrs_rho.ti.raw.rho; temp_1d.ti.raw.t =temp_1d.cxrs_rho.t;
if ~isempty(temp_1d.cxrs_rho.ti.raw.data)
data_fullpath_fit = 'Ti from fit from cxrs thus Ti(C)';
temp_1d.ti.raw =get_grids_1d(temp_1d.ti.raw,2,1);
else
data_fullpath_fit = 'Ti from fit in CONF node, using Te/Ti and equilibrium Wmhd';
end
temp_1d.ti.fit = temp_1d.ti_conf_rho;
temp_1d.ti.fit =get_grids_1d(temp_1d.ti.fit,1,1);
temp_1d.ni.fit = temp_1d.ni_conf_rho;
temp_1d.ni.fit =get_grids_1d(temp_1d.ni.fit,1,1);
it_ti = iround_os(temp_1d.ti.fit.t,ids_cores_profiles.time);
for it=1:length(ids_cores_profiles.time)
ids_cores_profiles.profiles_1d{it}.ion{1}.element{1}.a = 2.;
ids_cores_profiles.profiles_1d{it}.ion{1}.element{1}.z_n = 1;
ids_cores_profiles.profiles_1d{it}.ion{1}.element{1}.atoms_n = 1;
ids_cores_profiles.profiles_1d{it}.ion{1}.element{1}.z_ion = 1;
ids_cores_profiles.profiles_1d{it}.ion{1}.element{1}.label = 'D+';
ids_cores_profiles.profiles_1d{it}.ion{1}.multiple_states_flag = 0;
ids_cores_profiles.profiles_1d{it}.ion{1}.temperature = temp_1d.ti.fit.data(:,it_ti(it));
ids_cores_profiles.profiles_1d{it}.ion{1}.density = temp_1d.ni.fit.data(:,it_ti(it));
ids_cores_profiles.profiles_1d{it}.ion{1}.pressure_thermal = 1.6022e-19.*ids_cores_profiles.profiles_1d{it}.ion{1}.density ...
.* ids_cores_profiles.profiles_1d{it}.ion{1}.temperature;
%
ids_cores_profiles.profiles_1d{it}.t_i_average = ids_cores_profiles.profiles_1d{it}.ion{1}.temperature;
ids_cores_profiles.profiles_1d{it}.n_i_thermal_total = ids_cores_profiles.profiles_1d{it}.ion{1}.density;
ids_cores_profiles.profiles_1d{it}.pressure_ion_total = 1.6022e-19 .* ids_cores_profiles.profiles_1d{it}.n_i_thermal_total ...
.* ids_cores_profiles.profiles_1d{it}.t_i_average;
ids_cores_profiles.profiles_1d{it}.pressure_thermal = ids_cores_profiles.profiles_1d{it}.pressure_ion_total ...
+ ids_cores_profiles.profiles_1d{it}.electrons.pressure_thermal;
end
if ~isempty(temp_1d.cxrs_rho.ti.fit.data)
it_raw = iround_os(temp_1d.ti.raw.t,ids_cores_profiles.time);
for it=1:length(ids_cores_profiles.time)
% ids_cores_profiles.profiles_1d{it}.ion{1}.temperature_fit = temp_1d.ti.fit(:,it_ti(it));
ids_cores_profiles.profiles_1d{it}.ion{1}.density_fit.source = 'from Zeff and ne profile';
ids_cores_profiles.profiles_1d{it}.t_i_average_fit.measured = temp_1d.ti.raw.data(:,it_raw(it));
ids_cores_profiles.profiles_1d{it}.t_i_average_fit.source = 'from CXRS on C usually';
end
end
keyboard
%
profiles_1d.dpressure_dpsi = gdat(params_cores_profiles.shot,'pprime','machine',machine);
% dpsi_drho_tor = gdat(params_cores_profiles.shot,'dpsi_drho_tor','machine',machine);
% dvolume_dpsi = gdat(params_cores_profiles.shot,'dvolume_dpsi','machine',machine);
% dvolume_drho_tor = gdat(params_cores_profiles.shot,'dvolume_drho_tor','machine',machine);
% elongation = gdat(params_cores_profiles.shot,'elongation','machine',machine);
profiles_1d.f_df_dpsi = gdat(params_cores_profiles.shot,'ttprime','machine',machine);
profiles_1d.f = gdat(params_cores_profiles.shot,'rbphi_rho','machine',machine);
% geometric_axis = gdat(params_cores_profiles.shot,'geometric_axis','machine',machine);
% gm1 = gdat(params_cores_profiles.shot,'gm1','machine',machine);
% gm2 = gdat(params_cores_profiles.shot,'gm2','machine',machine);
% gm3 = gdat(params_cores_profiles.shot,'gm3','machine',machine);
% gm4 = gdat(params_cores_profiles.shot,'gm4','machine',machine);
% gm5 = gdat(params_cores_profiles.shot,'gm5','machine',machine);
% gm6 = gdat(params_cores_profiles.shot,'gm6','machine',machine);
% gm7 = gdat(params_cores_profiles.shot,'gm7','machine',machine);
% gm8 = gdat(params_cores_profiles.shot,'gm8','machine',machine);
% gm9 = gdat(params_cores_profiles.shot,'gm9','machine',machine);
% j_parallel = gdat(params_cores_profiles.shot,'j_parallel','machine',machine);
% j_tor = gdat(params_cores_profiles.shot,'j_tor','machine',machine);
% magnetic_shear = gdat(params_cores_profiles.shot,'magnetic_shear','machine',machine);
% mass_density = gdat(params_cores_profiles.shot,'mass_density','machine',machine);
profiles_1d.phi = gdat(params_cores_profiles.shot,'phi_tor','machine',machine);
profiles_1d.pressure = gdat(params_cores_profiles.shot,'pressure','machine',machine);
% psi = gdat(params_cores_profiles.shot,'psi_rho','machine',machine); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities)
profiles_1d.q = gdat(params_cores_profiles.shot,'q_rho','machine',machine);
profiles_1d.rho_tor = gdat(params_cores_profiles.shot,'rhotor','machine',machine);
%rho_tor_norm = gdat(params_cores_profiles.shot,'rhotor_norm','machine',machine); % from rho_tor
profiles_1d.rho_volume_norm = gdat(params_cores_profiles.shot,'rhovol','machine',machine);
% r_inboard = gdat(params_cores_profiles.shot,'r_inboard','machine',machine);
% r_outboard = gdat(params_cores_profiles.shot,'r_outboard','machine',machine);
% squareness_lower_inner = gdat(params_cores_profiles.shot,'squareness_lower_inner','machine',machine);
% squareness_lower_outer = gdat(params_cores_profiles.shot,'squareness_lower_outer','machine',machine);
% squareness_upper_inner = gdat(params_cores_profiles.shot,'squareness_upper_inner','machine',machine);
% squareness_upper_outer = gdat(params_cores_profiles.shot,'squareness_upper_outer','machine',machine);
% surface = gdat(params_cores_profiles.shot,'surface','machine',machine);
% trapped_fraction = gdat(params_cores_profiles.shot,'trapped_fraction','machine',machine);
% triangularity_lower = gdat(params_cores_profiles.shot,'triangularity_lower','machine',machine);
% triangularity_upper = gdat(params_cores_profiles.shot,'triangularity_upper','machine',machine);
profiles_1d.volume = gdat(params_cores_profiles.shot,'volume_rho','machine',machine);
profiles_1d_fieldnames = fieldnames(profiles_1d);
special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
for it=1:length(ids_cores_profiles.time)
for i=1:length(profiles_1d_fieldnames)
if ~any(strcmp(profiles_1d_fieldnames{i},special_fields))
if ~isstruct(ids_cores_profiles.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}))
if ~ischar(profiles_1d.(profiles_1d_fieldnames{i}).data) && ~isempty(profiles_1d.(profiles_1d_fieldnames{i}).data) ...
&& size(profiles_1d.(profiles_1d_fieldnames{i}).data,2)>=it
ids_cores_profiles.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}) = ...
profiles_1d.(profiles_1d_fieldnames{i}).data(:,it);
end
else
special_fields{end+1} = profiles_1d_fieldnames{i};
end
end
end
end
% special cases
nrho = length(profiles_1d.rho_tor.x);
ntime = length(temp.psi_axis.data);
for it=1:length(ids_cores_profiles.time)
ids_cores_profiles.time_slice{it}.profiles_1d.rho_tor_norm = ids_cores_profiles.time_slice{it}.profiles_1d.rho_tor./ ...
ids_cores_profiles.time_slice{it}.profiles_1d.rho_tor(end);
ids_cores_profiles.time_slice{it}.profiles_1d.psi = temp.psi_axis.data(it) + ...
profiles_1d.rho_tor.x.^2.*(global_quantities.psi_boundary.data(it)-temp.psi_axis.data(it));
end
%
%% profiles_2d{1} ala eqdsk, only this one thus grid_type=1
%
% b_field_r = gdat(params_cores_profiles.shot,'b_field_r','machine',machine);
% b_field_tor = gdat(params_cores_profiles.shot,'b_field_tor','machine',machine);
% b_field_z = gdat(params_cores_profiles.shot,'b_field_z','machine',machine);
% b_r = gdat(params_cores_profiles.shot,'b_r','machine',machine);
% b_tor = gdat(params_cores_profiles.shot,'b_tor','machine',machine);
% b_z = gdat(params_cores_profiles.shot,'b_z','machine',machine);
% grid = gdat(params_cores_profiles.shot,'grid','machine',machine); % special
profiles_2d.grid_type.name = 'rectangular';
profiles_2d.grid_type.index = 1;
profiles_2d.grid_type.description = 'Cylindrical R,Z ala eqdsk';
% j_parallel = gdat(params_cores_profiles.shot,'j_parallel','machine',machine);
% j_tor = gdat(params_cores_profiles.shot,'j_tor','machine',machine);
% phi = gdat(params_cores_profiles.shot,'phi','machine',machine);
profiles_2d.psi = gdat(params_cores_profiles.shot,'psi','machine',machine); % add psi_bound in a second step in special cases
% r = gdat(params_cores_profiles.shot,'r','machine',machine); % not to be filled since in grid.dim1
% theta = gdat(params_cores_profiles.shot,'theta','machine',machine);
% z = gdat(params_cores_profiles.shot,'z','machine',machine); % not to be filled since in grid.dim2
profiles_2d_fieldnames = fieldnames(profiles_2d);
special_fields = {'grid', 'grid_type'}; % fields needing non-automatic treatments
for it=1:length(ids_cores_profiles.time)
for i=1:length(profiles_2d_fieldnames)
if ~any(strcmp(profiles_2d_fieldnames{i},special_fields))
if ~isstruct(ids_cores_profiles.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}))
if ~ischar(profiles_2d.(profiles_2d_fieldnames{i}).data) && ~isempty(profiles_2d.(profiles_2d_fieldnames{i}).data) ...
&& size(profiles_2d.(profiles_2d_fieldnames{i}).data,3)>=it
ids_cores_profiles.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}) = ...
profiles_2d.(profiles_2d_fieldnames{i}).data(:,:,it);
end
else
special_fields{end+1} = profiles_2d_fieldnames{i};
end
end
end
end
% special cases
for it=1:length(ids_cores_profiles.time)
ids_cores_profiles.time_slice{it}.profiles_2d{1}.grid_type.name = profiles_2d.grid_type.name;
ids_cores_profiles.time_slice{it}.profiles_2d{1}.grid_type.index = profiles_2d.grid_type.index;
ids_cores_profiles.time_slice{it}.profiles_2d{1}.grid_type.description = profiles_2d.grid_type.description;
ids_cores_profiles.time_slice{it}.profiles_2d{1}.grid.dim1 = profiles_2d.psi.dim{1};
ids_cores_profiles.time_slice{it}.profiles_2d{1}.grid.dim2 = profiles_2d.psi.dim{2};
ids_cores_profiles.time_slice{it}.profiles_2d{1}.psi(:,:) = ids_cores_profiles.time_slice{it}.profiles_2d{1}.psi(:,:) + ...
global_quantities.psi_boundary.data(it);
end
% make arrays not filled in empty:
ids_cores_profiles.grids_ggd = {};
for it=1:length(ids_cores_profiles.time_slice)
ids_cores_profiles.time_slice{it}.ggd = {};
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment