From 8d39fd16b143dad243390851beca4f566307d1d1 Mon Sep 17 00:00:00 2001
From: Antonia Frank <antonia.frank@epfl.ch>
Date: Tue, 17 Oct 2023 14:47:19 +0200
Subject: [PATCH] Clean up core_profiles ids and add jbs and johm
---
matlab/TCV_IMAS/tcv_get_ids_core_profiles.m | 75 ++++++++++++++-------
1 file changed, 49 insertions(+), 26 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
index 2891e6dd..d618d754 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
@@ -1,10 +1,12 @@
-function [ids_core_profiles,ids_core_profiles_description,varargout] = tcv_get_ids_core_profiles(shot,ids_equil_empty, gdat_params,varargin)
+function [ids_core_profiles,ids_core_profiles_description,varargout] = ...
+ tcv_get_ids_core_profiles(shot,ids_equil_empty,gdat_params,varargin)
%
-% [ids_core_profiles,ids_core_profiles_description,varargout] = tcv_get_ids_core_profiles(shot,ids_equil_empty,varargin);
+% [ids_core_profiles,ids_core_profiles_description,varargout] = ...
+% tcv_get_ids_core_profiles(shot,ids_equil_empty,gdat_params,varargin);
%
%
-% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar options
-%
+% gdat_params: gdat_data.gdat_params to get all params passed from original call,
+% in particular error_bar options
%
% error_bar: from gdat_params.error_bar
% 'delta' (default): error_bar to be added inserted in "upper" only as mentioned in description
@@ -22,18 +24,22 @@ tens_time = -1;
tens_rho = -0.1;
if exist('gdat_params','var')
- [ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles','gdat_params',gdat_params,varargin{:});
+ [ids_core_profiles, params_cores_profiles] = ...
+ tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles','gdat_params',gdat_params,varargin{:});
else
- [ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',varargin{:});
+ [ids_core_profiles, params_cores_profiles] = ...
+ tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',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_profiles_description = [];
-% base all from times fro nete_rho (which should be conf by default)
+%% setup profiles_1d
+% base all from times for nete_rho (which should be conf by default)
params_eff = params_eff_ref;
params_eff.data_request='ne_rho'; params_eff.fit = 1;
temp_1d.ne_rho = gdat(params_cores_profiles.shot,params_eff);
@@ -70,13 +76,11 @@ ids_core_profiles.profiles_1d{1}.ion{1}.state = {};
ids_core_profiles.profiles_1d{1}.neutral{1}.state = {};
ids_core_profiles.profiles_1d(1:length(ids_core_profiles.time)) = ids_core_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 and time, using homogeneous
params_eff = params_eff_ref;
params_eff.data_request='b0';
vacuum_toroidal_field.b0=gdat(params_cores_profiles.shot,params_eff);
@@ -85,7 +89,10 @@ ids_core_profiles.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
ids_core_profiles.vacuum_toroidal_field.b0 = interpos(63,vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data,ids_core_profiles.time,-1);
ids_core_profiles_description.vacuum_toroidal_field = [vacuum_toroidal_field_desc.b0 ' on ids_core_profiles.time, with interpos(63)'];
-% 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
+% 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
params_eff.data_request='ip';
global_quantities.ip = gdat(params_cores_profiles.shot,'ip','machine',machine);
@@ -141,25 +148,20 @@ for i=1:length(global_quantities_fieldnames_eff)
end
end
end
-% special case
-
-%
%% profiles_1d (cannot use eqdsk since not same radial mesh)
+
+% prepare area for ids_core_profiles.profiles_1d{*}.grid.area
params_eff.data_request = 'area_rho';
temp_1d.area = gdat(params_cores_profiles.shot,params_eff);
temp_1d_desc.area = params_eff.data_request;
-for ir=1:length(temp_1d.area.x)
+for ir=1:length(temp_1d.area.x) % map tmp_1d.area to core_profiles.time
area_cpt(ir,:) = interpos(temp_1d.area.t,temp_1d.area.data(ir,:),ids_core_profiles.time,tens_time);
end
-params_eff.data_request = 'q_rho';
-temp_1d.q = gdat(params_cores_profiles.shot,params_eff);
-temp_1d_desc.q = params_eff.data_request;
-for ir=1:length(temp_1d.q.x)
- q_cpt(ir,:) = interpos(temp_1d.q.t,temp_1d.q.data(ir,:),ids_core_profiles.time,tens_time);
-end
+
it_thom = iround_os(temp_1d.te_rho.t,ids_core_profiles.time);
for it=1:length(ids_core_profiles.time)
+ % fill grid
ids_core_profiles.profiles_1d{it}.grid.rho_tor_norm = temp_1d.fit.te_rho.grids_1d.rhotornorm(:,it);
ids_core_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);
@@ -171,7 +173,9 @@ for it=1:length(ids_core_profiles.time)
.* temp_1d.fit.te_rho.grids_1d.volume_edge(it);
ids_core_profiles.profiles_1d{it}.grid.area = interpos(temp_1d.area.x,area_cpt(:,it),temp_1d.fit.te_rho.grids_1d.rhopolnorm, ...
tens_rho,[1 2],[0 area_cpt(end,it)]);
+ % fill time
ids_core_profiles.profiles_1d{it}.time = ids_core_profiles.time(it);
+ % fill electrons struct
ids_core_profiles.profiles_1d{it}.electrons.temperature = temp_1d.fit.te_rho.data(:,it);
ids_core_profiles.profiles_1d{it}.electrons.temperature_fit.measured = temp_1d.te_rho.data(:,it_thom(it));
ids_core_profiles.profiles_1d{it}.electrons.temperature_fit.time_measurement = temp_1d.te_rho.t(it_thom(it));
@@ -185,6 +189,7 @@ for it=1:length(ids_core_profiles.time)
ids_core_profiles.profiles_1d{it}.electrons.density_fit.source = {'Thomson, interpos fit'};
ids_core_profiles.profiles_1d{it}.electrons.pressure_thermal = 1.6022e-19.*ids_core_profiles.profiles_1d{it}.electrons.density_thermal ...
.* ids_core_profiles.profiles_1d{it}.electrons.temperature;
+ % fill zeff
ids_core_profiles.profiles_1d{it}.zeff = global_quantities.z_eff_resistive.data(it) .* ...
ones(size(ids_core_profiles.profiles_1d{it}.electrons.density));
end
@@ -231,7 +236,8 @@ switch error_bar
error(['tcv_ids_bpol_loop: error_bar option not known: ' error_bar])
end
-% ion, assume only D if no CXRS (need to ask how to check if H...)
+%% ion struct
+% assume only D if no CXRS (need to ask how to check if H...)
params_eff_fit1.data_request = 'cxrs';
temp_1d.cxrs_rho = gdat(params_cores_profiles.shot,params_eff_fit1);
temp_1d_desc.cxrs_rho = params_eff_fit1.data_request;
@@ -316,9 +322,14 @@ if ~isempty(temp_1d.cxrs_rho.ti.fit.data)
end
end
-for ir=1:length(temp_1d.q.x)
+%% q profile and magnetic shear
+params_eff.data_request = 'q_rho';
+temp_1d.q = gdat(params_cores_profiles.shot,params_eff);
+temp_1d_desc.q = params_eff.data_request;
+for ir=1:length(temp_1d.q.x) % map tmp_1d.q to core_profiles.time
q_cpt(ir,:) = interpos(temp_1d.q.t,temp_1d.q.data(ir,:),ids_core_profiles.time,tens_time);
end
+
for it=1:length(ids_core_profiles.time)
ij=isfinite(q_cpt(:,it));
if sum(ij) >= 5
@@ -335,9 +346,21 @@ for it=1:length(ids_core_profiles.time)
ids_core_profiles.profiles_1d{it}.magnetic_shear = -9.e40;
end
end
+temp_1d_desc.magnetic_shear = 'from interpos with rhotor';
+
+%% Current densities
+
+for it=1:length(ids_core_profiles.time)
+ ids_core_profiles.profiles_1d{it}.j_bootstrap = ...
+ current_bootstrap.bs.bs_data.cd_dens.data(:,it);
+ temp_1d_desc.j_bootstrap = current_bootstrap.bs.bs_data.cd_dens.label;
+ ids_core_profiles.profiles_1d{it}.j_ohmic = ...
+ current_non_inductive.ohm.ohm_data.cd_dens.data(:,it);
+ temp_1d_desc.j_ohmic = current_non_inductive.ohm.ohm_data.cd_dens.label;
+end
-ids_core_profiles_description.temp_1d_desc = temp_1d_desc;
-ids_core_profiles_description.profiles_1d.magnetic_shear = 'from interpos with rhotor';
+%% add descriptions for profiles_1d
+ids_core_profiles_description.profiles_1d = temp_1d_desc;
if nargin <= 2
ids_core_profiles.code.name = ['tcv_get_ids_core_profiles, within gdat, with shot= ' num2str(params_cores_profiles.shot) ];
--
GitLab