diff --git a/crpptbx/TCV_IMAS/tcv_get_ids_core_profiles.m b/crpptbx/TCV_IMAS/tcv_get_ids_core_profiles.m new file mode 100644 index 0000000000000000000000000000000000000000..1f6d82cfbb3fc038993e36dd0aeafcdc242bf1bd --- /dev/null +++ b/crpptbx/TCV_IMAS/tcv_get_ids_core_profiles.m @@ -0,0 +1,299 @@ +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