From 0320a20446df0edabf18be0040fd5c91db56c879 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 12 Jun 2019 21:57:10 +0000
Subject: [PATCH] add some extra ids elements
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@12059 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/TCV_IMAS/tcv_get_ids_equilibrium.m | 266 ++++++++++++---------
crpptbx/TCV_IMAS/tcv_ids_bpol_probe.m | 37 ++-
crpptbx/TCV_IMAS/tcv_ids_flux_loop.m | 15 +-
3 files changed, 194 insertions(+), 124 deletions(-)
diff --git a/crpptbx/TCV_IMAS/tcv_get_ids_equilibrium.m b/crpptbx/TCV_IMAS/tcv_get_ids_equilibrium.m
index 8c8d5515..bc793c10 100644
--- a/crpptbx/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/crpptbx/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -1,13 +1,11 @@
function [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,ids_equil_empty, gdat_params,varargin);
%
% [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,ids_equil_empty,varargin);
-%
+%
%
% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar options
%
-machine = 'tcv';
-
[ids_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium',varargin{:});
% As a general rule, for a new substructure under the main ids, construct a local structure like:
@@ -16,11 +14,16 @@ machine = 'tcv';
%
% vacuum_toroidal_field and time, using homogeneous
%
-vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,'b0','source','liuqe','machine',machine); % to get on liuqe time array
-vacuum_toroidal_field_desc.b0 = '''b0'',''source'',''liuqe''';
+%% liuqe.m at this stage is missing correction 0.996, so use std source by time of default liuqe to make sure
+bb = gdat(params_equilibrium.shot,'b0','source','liuqe','machine',gdat_params.machine,'liuqe',gdat_params.liuqe); % to get liuqe time array
+vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,'b0','machine',gdat_params.machine);
+vacuum_toroidal_field.b0.data = interp1(vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data,bb.t);
+vacuum_toroidal_field.b0.t = bb.t;
+vacuum_toroidal_field.b0.dim = {vacuum_toroidal_field.b0.t};
+vacuum_toroidal_field_desc.b0 = ['''b0'',''timing source'',''liuqe=' num2str(gdat_params.liuqe) ''''];
vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'',''source'',''liuqe'']';
ids_equilibrium.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
-ids_equilibrium.vacuum_toroidal_field.b0 = vacuum_toroidal_field.b0.data;
+ids_equilibrium.vacuum_toroidal_field.b0.data = vacuum_toroidal_field.b0.data;
ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
ids_equilibrium.time = vacuum_toroidal_field.b0.t;
@@ -35,43 +38,43 @@ ids_equilibrium.time_slice(1:numel(ids_equilibrium.time)) = ids_equilibrium.time
% brute force solution load all eqdsks
% $$$ for it=1:numel(ids_equilibrium.time)
% $$$ ids_equilibrium.time(it)
-% $$$ temp.eqdsks{it}=gdat(params_equilibrium.shot,'eqdsk','time',ids_equilibrium.time(it),'write',0,'machine',machine);
+% $$$ temp.eqdsks{it}=gdat(params_equilibrium.shot,'eqdsk','time',ids_equilibrium.time(it),'write',0,'machine',gdat_params.machine);
% $$$ end
% $$$ temp_desc.eqdsks{1} = '''eqdsk'',''time'',ids_equilibrium.time(it)';
-global_quantities.area = gdat(params_equilibrium.shot,'area_edge','machine',machine);
+global_quantities.area = gdat(params_equilibrium.shot,'area_edge','machine',gdat_params.machine);
global_quantities_desc.area = 'area_edge';
-global_quantities.beta_normal = gdat(params_equilibrium.shot,'betan','machine',machine);
+global_quantities.beta_normal = gdat(params_equilibrium.shot,'betan','machine',gdat_params.machine);
global_quantities_desc.beta_normal = 'betan';
-global_quantities.beta_pol = gdat(params_equilibrium.shot,'betap','machine',machine);
+global_quantities.beta_pol = gdat(params_equilibrium.shot,'betap','machine',gdat_params.machine);
global_quantities_desc.beta_pol = 'betap';
-global_quantities.beta_tor = gdat(params_equilibrium.shot,'beta','machine',machine);
+global_quantities.beta_tor = gdat(params_equilibrium.shot,'beta','machine',gdat_params.machine);
global_quantities_desc.beta_tor = 'beta';
-global_quantities.energy_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',machine);
+global_quantities.energy_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',gdat_params.machine);
global_quantities_desc.energy_mhd = 'w_mhd';
-global_quantities.ip = gdat(params_equilibrium.shot,'ip','machine',machine);
+global_quantities.ip = gdat(params_equilibrium.shot,'ip','machine',gdat_params.machine);
global_quantities_desc.ip = 'ip';
-% length_pol = gdat(params_equilibrium.shot,'length_pol','machine',machine); % to be added
-global_quantities.li_3 = gdat(params_equilibrium.shot,'li','machine',machine);
+% length_pol = gdat(params_equilibrium.shot,'length_pol','machine',gdat_params.machine); % to be added
+global_quantities.li_3 = gdat(params_equilibrium.shot,'li','machine',gdat_params.machine);
global_quantities_desc.li_3 = 'li';
-temp.r_magnetic_axis = gdat(params_equilibrium.shot,'r_axis','machine',machine);
+temp.r_magnetic_axis = gdat(params_equilibrium.shot,'r_axis','machine',gdat_params.machine);
temp_desc.r_magnetic_axis = 'r_axis';
-temp.z_magnetic_axis = gdat(params_equilibrium.shot,'z_axis','machine',machine);
+temp.z_magnetic_axis = gdat(params_equilibrium.shot,'z_axis','machine',gdat_params.machine);
temp_desc.z_magnetic_axis = 'z_axis';
-temp.psi_axis = gdat(params_equilibrium.shot,'psi_axis','machine',machine); % needs to add psi_edge sincepsi_axis liuqe assuming 0 dege value
+temp.psi_axis = gdat(params_equilibrium.shot,'psi_axis','machine',gdat_params.machine); % needs to add psi_edge sincepsi_axis liuqe assuming 0 dege value
temp_desc.psi_axis = 'psi_axis';
-global_quantities.psi_boundary = gdat(params_equilibrium.shot,'psi_edge','machine',machine);
+global_quantities.psi_boundary = gdat(params_equilibrium.shot,'psi_edge','machine',gdat_params.machine);
global_quantities_desc.psi_boundary = 'psi_edge';
-global_quantities.q_95 = gdat(params_equilibrium.shot,'q95','machine',machine);
+global_quantities.q_95 = gdat(params_equilibrium.shot,'q95','machine',gdat_params.machine);
global_quantities_desc.q_95 = 'q95';
-global_quantities.q_axis = gdat(params_equilibrium.shot,'q0','machine',machine); % will be checked with q_rho?
+global_quantities.q_axis = gdat(params_equilibrium.shot,'q0','machine',gdat_params.machine); % will be checked with q_rho?
global_quantities_desc.q_axis = 'q0';
-temp.q_rho = gdat(params_equilibrium.shot,'q_rho','machine',machine);
+temp.q_rho = gdat(params_equilibrium.shot,'q_rho','machine',gdat_params.machine);
temp_desc.q_rho = 'q_rho';
-% surface = gdat(params_equilibrium.shot,'surface','machine',machine); % to be added
-global_quantities.volume = gdat(params_equilibrium.shot,'volume','machine',machine);
+% surface = gdat(params_equilibrium.shot,'surface','machine',gdat_params.machine); % to be added
+global_quantities.volume = gdat(params_equilibrium.shot,'volume','machine',gdat_params.machine);
global_quantities_desc.volume = 'volume';
-global_quantities.w_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',machine);
+global_quantities.w_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',gdat_params.machine);
global_quantities_desc.w_mhd = 'w_mhd';
global_quantities_fieldnames = fieldnames(global_quantities);
@@ -101,37 +104,37 @@ for it=1:numel(ids_equilibrium.time)
end
% for boundary in addition to lcfs
-% active_limiter_point = gdat(params_equilibrium.shot,'active_limiter_point','machine',machine);
-boundary.elongation = gdat(params_equilibrium.shot,'kappa','machine',machine);
+% active_limiter_point = gdat(params_equilibrium.shot,'active_limiter_point','machine',gdat_params.machine);
+boundary.elongation = gdat(params_equilibrium.shot,'kappa','machine',gdat_params.machine);
boundary_desc.elongation = 'kappa';
-% elongation_lower = gdat(params_equilibrium.shot,'elongation_lower','machine',machine);
-% elongation_upper = gdat(params_equilibrium.shot,'elongation_upper','machine',machine);
-boundary.minor_radius = gdat(params_equilibrium.shot,'a_minor','machine',machine);
+% elongation_lower = gdat(params_equilibrium.shot,'elongation_lower','machine',gdat_params.machine);
+% elongation_upper = gdat(params_equilibrium.shot,'elongation_upper','machine',gdat_params.machine);
+boundary.minor_radius = gdat(params_equilibrium.shot,'a_minor','machine',gdat_params.machine);
boundary_desc.minor_radius = 'a_minor';
-% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',machine);
-% squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer','machine',machine);
-% squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner','machine',machine);
-% squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer','machine',machine);
-% strike_point = gdat(params_equilibrium.shot,'strike_point','machine',machine);
-boundary.triangularity = gdat(params_equilibrium.shot,'delta','machine',machine);
+% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',gdat_params.machine);
+% squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer','machine',gdat_params.machine);
+% squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner','machine',gdat_params.machine);
+% squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer','machine',gdat_params.machine);
+% strike_point = gdat(params_equilibrium.shot,'strike_point','machine',gdat_params.machine);
+boundary.triangularity = gdat(params_equilibrium.shot,'delta','machine',gdat_params.machine);
boundary_desc.triangularity = 'delta';
-boundary.triangularity_lower = gdat(params_equilibrium.shot,'delta_bottom','machine',machine);
+boundary.triangularity_lower = gdat(params_equilibrium.shot,'delta_bottom','machine',gdat_params.machine);
boundary_desc.triangularity_lower = 'delta_bottom';
-boundary.triangularity_upper = gdat(params_equilibrium.shot,'delta_top','machine',machine);
+boundary.triangularity_upper = gdat(params_equilibrium.shot,'delta_top','machine',gdat_params.machine);
boundary_desc.triangularity_upper = 'delta_top';
-temp.n_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''n_xpts'''',''''liuqe.m'''')','machine',machine);
+temp.n_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''n_xpts'''',''''liuqe.m'''')','machine',gdat_params.machine);
temp_desc.n_x_point = '''tcv_eq(''''n_xpts'''',''''liuqe.m'''')''';
-temp.r_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''r_xpts'''',''''liuqe.m'''')','machine',machine);
+temp.r_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''r_xpts'''',''''liuqe.m'''')','machine',gdat_params.machine);
temp_desc.r_x_point = '''tcv_eq(''''r_xpts'''',''''liuqe.m'''')''';
-temp.z_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''z_xpts'''',''''liuqe.m'''')','machine',machine);
+temp.z_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''z_xpts'''',''''liuqe.m'''')','machine',gdat_params.machine);
temp_desc.z_x_point = '''tcv_eq(''''z_xpts'''',''''liuqe.m'''')''';
-temp.rgeom = gdat(params_equilibrium.shot,'rgeom','machine',machine);
+temp.rgeom = gdat(params_equilibrium.shot,'rgeom','machine',gdat_params.machine);
temp_desc.rgeom = 'rgeom';
-temp.zgeom = gdat(params_equilibrium.shot,'zgeom','machine',machine);
+temp.zgeom = gdat(params_equilibrium.shot,'zgeom','machine',gdat_params.machine);
temp_desc.zgeom = 'zgeom';
-temp.r_lcfs = gdat(params_equilibrium.shot,'r_contour_edge','machine',machine);
+temp.r_lcfs = gdat(params_equilibrium.shot,'r_contour_edge','machine',gdat_params.machine);
temp_desc.r_lcfs = 'r_contour_edge';
-temp.z_lcfs = gdat(params_equilibrium.shot,'z_contour_edge','machine',machine);
+temp.z_lcfs = gdat(params_equilibrium.shot,'z_contour_edge','machine',gdat_params.machine);
temp_desc.z_lcfs = 'z_contour_edge';
boundary_fieldnames = fieldnames(boundary);
@@ -171,55 +174,99 @@ end
%
%% profiles_1d (cannot use eqdsk since not same radial mesh)
%
-% area = gdat(params_equilibrium.shot,'area','machine',machine);
-% b_average = gdat(params_equilibrium.shot,'b_average','machine',machine);
-% beta_pol = gdat(params_equilibrium.shot,'beta_pol','machine',machine);
-% b_field_average = gdat(params_equilibrium.shot,'b_field_average','machine',machine);
-% b_field_max = gdat(params_equilibrium.shot,'b_field_max','machine',machine);
-% b_field_min = gdat(params_equilibrium.shot,'b_field_min','machine',machine);
-% b_max = gdat(params_equilibrium.shot,'b_max','machine',machine);
-% b_min = gdat(params_equilibrium.shot,'b_min','machine',machine);
-% darea_dpsi = gdat(params_equilibrium.shot,'darea_dpsi','machine',machine);
-% darea_drho_tor = gdat(params_equilibrium.shot,'darea_drho_tor','machine',machine);
-profiles_1d.dpressure_dpsi = gdat(params_equilibrium.shot,'pprime','machine',machine);
-% dpsi_drho_tor = gdat(params_equilibrium.shot,'dpsi_drho_tor','machine',machine);
-% dvolume_dpsi = gdat(params_equilibrium.shot,'dvolume_dpsi','machine',machine);
-% dvolume_drho_tor = gdat(params_equilibrium.shot,'dvolume_drho_tor','machine',machine);
-% elongation = gdat(params_equilibrium.shot,'elongation','machine',machine);
-profiles_1d.f_df_dpsi = gdat(params_equilibrium.shot,'ttprime','machine',machine);
-profiles_1d.f = gdat(params_equilibrium.shot,'rbphi_rho','machine',machine);
-% geometric_axis = gdat(params_equilibrium.shot,'geometric_axis','machine',machine);
-% gm1 = gdat(params_equilibrium.shot,'gm1','machine',machine);
-% gm2 = gdat(params_equilibrium.shot,'gm2','machine',machine);
-% gm3 = gdat(params_equilibrium.shot,'gm3','machine',machine);
-% gm4 = gdat(params_equilibrium.shot,'gm4','machine',machine);
-% gm5 = gdat(params_equilibrium.shot,'gm5','machine',machine);
-% gm6 = gdat(params_equilibrium.shot,'gm6','machine',machine);
-% gm7 = gdat(params_equilibrium.shot,'gm7','machine',machine);
-% gm8 = gdat(params_equilibrium.shot,'gm8','machine',machine);
-% gm9 = gdat(params_equilibrium.shot,'gm9','machine',machine);
-% j_parallel = gdat(params_equilibrium.shot,'j_parallel','machine',machine);
-% j_tor = gdat(params_equilibrium.shot,'j_tor','machine',machine);
-% magnetic_shear = gdat(params_equilibrium.shot,'magnetic_shear','machine',machine);
-% mass_density = gdat(params_equilibrium.shot,'mass_density','machine',machine);
-profiles_1d.phi = gdat(params_equilibrium.shot,'phi_tor','machine',machine);
-profiles_1d.pressure = gdat(params_equilibrium.shot,'pressure','machine',machine);
-% psi = gdat(params_equilibrium.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_equilibrium.shot,'q_rho','machine',machine);
-profiles_1d.rho_tor = gdat(params_equilibrium.shot,'rhotor','machine',machine);
-%rho_tor_norm = gdat(params_equilibrium.shot,'rhotor_norm','machine',machine); % from rho_tor
-profiles_1d.rho_volume_norm = gdat(params_equilibrium.shot,'rhovol','machine',machine);
-% r_inboard = gdat(params_equilibrium.shot,'r_inboard','machine',machine);
-% r_outboard = gdat(params_equilibrium.shot,'r_outboard','machine',machine);
-% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',machine);
-% squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer','machine',machine);
-% squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner','machine',machine);
-% squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer','machine',machine);
-% surface = gdat(params_equilibrium.shot,'surface','machine',machine);
-% trapped_fraction = gdat(params_equilibrium.shot,'trapped_fraction','machine',machine);
-% triangularity_lower = gdat(params_equilibrium.shot,'triangularity_lower','machine',machine);
-% triangularity_upper = gdat(params_equilibrium.shot,'triangularity_upper','machine',machine);
-profiles_1d.volume = gdat(params_equilibrium.shot,'volume_rho','machine',machine);
+% area = gdat(params_equilibrium.shot,'area','machine',gdat_params.machine);
+% b_average = gdat(params_equilibrium.shot,'b_average','machine',gdat_params.machine);
+% beta_pol = gdat(params_equilibrium.shot,'beta_pol','machine',gdat_params.machine);
+% b_field_average = gdat(params_equilibrium.shot,'b_field_average','machine',gdat_params.machine);
+% b_field_max = gdat(params_equilibrium.shot,'b_field_max','machine',gdat_params.machine);
+% b_field_min = gdat(params_equilibrium.shot,'b_field_min','machine',gdat_params.machine);
+% b_max = gdat(params_equilibrium.shot,'b_max','machine',gdat_params.machine);
+% b_min = gdat(params_equilibrium.shot,'b_min','machine',gdat_params.machine);
+% darea_dpsi = gdat(params_equilibrium.shot,'darea_dpsi','machine',gdat_params.machine);
+% darea_drho_tor = gdat(params_equilibrium.shot,'darea_drho_tor','machine',gdat_params.machine);
+profiles_1d.dpressure_dpsi = gdat(params_equilibrium.shot,'pprime','machine',gdat_params.machine);
+% dpsi_drho_tor = gdat(params_equilibrium.shot,'dpsi_drho_tor','machine',gdat_params.machine);
+% dvolume_dpsi = gdat(params_equilibrium.shot,'dvolume_dpsi','machine',gdat_params.machine);
+% dvolume_drho_tor = gdat(params_equilibrium.shot,'dvolume_drho_tor','machine',gdat_params.machine);
+% elongation = gdat(params_equilibrium.shot,'elongation','machine',gdat_params.machine);
+profiles_1d.f_df_dpsi = gdat(params_equilibrium.shot,'ttprime','machine',gdat_params.machine);
+profiles_1d.f = gdat(params_equilibrium.shot,'rbphi_rho','machine',gdat_params.machine);
+profiles_1d.f.data = 0.996 * profiles_1d.f.data;
+profiles_1d.f_df_dpsi.data = 0.996.^2 .* profiles_1d.f_df_dpsi.data;
+% geometric_axis = gdat(params_equilibrium.shot,'geometric_axis','machine',gdat_params.machine);
+% gm1 = gdat(params_equilibrium.shot,'gm1','machine',gdat_params.machine);
+% gm2 = gdat(params_equilibrium.shot,'gm2','machine',gdat_params.machine);
+% gm3 = gdat(params_equilibrium.shot,'gm3','machine',gdat_params.machine);
+% gm4 = gdat(params_equilibrium.shot,'gm4','machine',gdat_params.machine);
+% gm5 = gdat(params_equilibrium.shot,'gm5','machine',gdat_params.machine);
+% gm6 = gdat(params_equilibrium.shot,'gm6','machine',gdat_params.machine);
+% gm7 = gdat(params_equilibrium.shot,'gm7','machine',gdat_params.machine);
+% gm8 = gdat(params_equilibrium.shot,'gm8','machine',gdat_params.machine);
+% gm9 = gdat(params_equilibrium.shot,'gm9','machine',gdat_params.machine);
+% j_parallel = gdat(params_equilibrium.shot,'j_parallel','machine',gdat_params.machine);
+% j_tor = gdat(params_equilibrium.shot,'j_tor','machine',gdat_params.machine);
+% magnetic_shear = gdat(params_equilibrium.shot,'magnetic_shear','machine',gdat_params.machine);
+% mass_density = gdat(params_equilibrium.shot,'mass_density','machine',gdat_params.machine);
+profiles_1d.phi = gdat(params_equilibrium.shot,'phi_tor','machine',gdat_params.machine);
+profiles_1d.phi.data = 0.996 * profiles_1d.phi.data;
+profiles_1d.pressure = gdat(params_equilibrium.shot,'pressure','machine',gdat_params.machine);
+% psi = gdat(params_equilibrium.shot,'psi_rho','machine',gdat_params.machine); % (could take from .x of any like rhotor and psi_axis, psi_edge from global_quantities)
+profiles_1d.q = gdat(params_equilibrium.shot,'q_rho','machine',gdat_params.machine);
+profiles_1d.rho_tor = gdat(params_equilibrium.shot,'rhotor','machine',gdat_params.machine);
+%rho_tor_norm = gdat(params_equilibrium.shot,'rhotor_norm','machine',gdat_params.machine); % from rho_tor
+profiles_1d.rho_volume_norm = gdat(params_equilibrium.shot,'rhovol','machine',gdat_params.machine);
+% r_inboard = gdat(params_equilibrium.shot,'r_inboard','machine',gdat_params.machine);
+% r_outboard = gdat(params_equilibrium.shot,'r_outboard','machine',gdat_params.machine);
+% squareness_lower_inner = gdat(params_equilibrium.shot,'squareness_lower_inner','machine',gdat_params.machine);
+% squareness_lower_outer = gdat(params_equilibrium.shot,'squareness_lower_outer','machine',gdat_params.machine);
+% squareness_upper_inner = gdat(params_equilibrium.shot,'squareness_upper_inner','machine',gdat_params.machine);
+% squareness_upper_outer = gdat(params_equilibrium.shot,'squareness_upper_outer','machine',gdat_params.machine);
+% surface = gdat(params_equilibrium.shot,'surface','machine',gdat_params.machine);
+% trapped_fraction = gdat(params_equilibrium.shot,'trapped_fraction','machine',gdat_params.machine);
+% triangularity_lower = gdat(params_equilibrium.shot,'triangularity_lower','machine',gdat_params.machine);
+% triangularity_upper = gdat(params_equilibrium.shot,'triangularity_upper','machine',gdat_params.machine);
+profiles_1d.volume = gdat(params_equilibrium.shot,'volume_rho','machine',gdat_params.machine,'liuqe',gdat_params.liuqe);
+
+liuqe_opt = gdat_params.liuqe; % default at this stage but could take from gdat params like error bar
+switch liuqe_opt
+ case {-1}, psitbx_str='FBTE';
+ case {1,21}, psitbx_str='LIUQE.M';
+ case {11}, psitbx_str='LIUQE';
+ case {2, 3, 22, 23}, psitbx_str=['LIUQE.M' num2str(mod(liuqe_opt,10))];
+ case {12,13}, psitbx_str=['LIUQE' num2str(mod(liuqe_opt,10))];
+ otherwise, error(['Unknown LIUQE version, liuqe = ' num2str(liuqe_opt)]);
+end
+fsd = psitbxtcv2(shot,profiles_1d.volume.t,'FS',psitbx_str);
+grho_metric_3D = metric(fsd,-1);
+% Introduced new anonymous function to compute FS average ...
+metric_FS = metric(grho_metric_3D.grid,[2,3]);
+denom=sum(metric_FS./grho_metric_3D,[2,3]);
+FS_av = @(x) sum(x.*metric_FS./grho_metric_3D,[2,3])./denom;
+R=metric(fsd,3);
+Rm2av=FS_av(1./R.^2);
+profiles_1d.gm1.data = Rm2av.x;
+%tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % this gives (grad rhopol/R)^2 not gm2 which is grad rhotor^2
+%profiles_1d.gm2.data = tmp_gm.x;
+tmp_gm = FS_av(1./R.^1);
+profiles_1d.gm9.data = tmp_gm.x;
+
+tmp_gm = FS_av(grho_metric_3D.^2./R.^2); % grad rhopol^2 to get <grad psi^2>
+for it=1:numel(ids_equilibrium.time)
+ gradpsi_over_R_sq(:,it) = tmp_gm.x(:,it) .* 4 .* profiles_1d.volume.x.^2 .* ...
+ (ids_equilibrium.time_slice{it}.global_quantities.psi_boundary-ids_equilibrium.time_slice{it}.global_quantities.psi_axis).^2;
+end
+mu0 = 4.e-7 * pi;
+j_tor = -profiles_1d.dpressure_dpsi.data ./ profiles_1d.gm9.data ...
+ - profiles_1d.gm1.data ./ profiles_1d.gm9.data .* profiles_1d.f_df_dpsi.data ./ mu0;
+profiles_1d.j_tor.data = - 2.*pi.* j_tor; % 2pi sigma_bp * jtor above (Eq. (30) cocos paper cocos=17)
+% $$$ j_par = - ids_equilibrium.vacuum_toroidal_field.r0.^2 .* profiles_1d.f.data ...
+% $$$ ./repmat(profiles_1d.f.data(end,:),size(profiles_1d.f.data,1),1) ...
+% $$$ .* (profiles_1d.dpressure_dpsi.data ...
+% $$$ + profiles_1d.f_df_dpsi.data/mu0 .* (profiles_1d.gm1.data + gradpsi_sq ./ profiles_1d.f.data.^2);
+j_par = - profiles_1d.f_df_dpsi.data./profiles_1d.f.data./mu0.*gradpsi_over_R_sq./2./pi ...
+ - profiles_1d.f.data .*2*pi .*profiles_1d.dpressure_dpsi.data ...
+ - profiles_1d.f.data .*2*pi /mu0.*profiles_1d.f_df_dpsi.data.*profiles_1d.gm1.data;
+profiles_1d.j_parallel.data = j_par./repmat(ids_equilibrium.vacuum_toroidal_field.b0.data',size(profiles_1d.f.data,1),1);
profiles_1d_fieldnames = fieldnames(profiles_1d);
special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
@@ -243,6 +290,9 @@ end
nrho = numel(profiles_1d.rho_tor.x);
ntime = numel(temp.psi_axis.data);
for it=1:numel(ids_equilibrium.time)
+ ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.b_field_tor = ids_equilibrium.time_slice{it}.profiles_1d.f(1) ...
+ ./ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.r;
+ ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.b_tor = ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.b_field_tor;
ids_equilibrium.time_slice{it}.profiles_1d.rho_tor_norm = ids_equilibrium.time_slice{it}.profiles_1d.rho_tor./ ...
ids_equilibrium.time_slice{it}.profiles_1d.rho_tor(end);
ids_equilibrium.time_slice{it}.profiles_1d.psi = ids_equilibrium.time_slice{it}.global_quantities.psi_axis + ...
@@ -253,23 +303,23 @@ end
%
%% profiles_2d{1} ala eqdsk, only this one thus grid_type=1
%
-% b_field_r = gdat(params_equilibrium.shot,'b_field_r','machine',machine);
-% b_field_tor = gdat(params_equilibrium.shot,'b_field_tor','machine',machine);
-% b_field_z = gdat(params_equilibrium.shot,'b_field_z','machine',machine);
-% b_r = gdat(params_equilibrium.shot,'b_r','machine',machine);
-% b_tor = gdat(params_equilibrium.shot,'b_tor','machine',machine);
-% b_z = gdat(params_equilibrium.shot,'b_z','machine',machine);
-% grid = gdat(params_equilibrium.shot,'grid','machine',machine); % special
+% b_field_r = gdat(params_equilibrium.shot,'b_field_r','machine',gdat_params.machine);
+% b_field_tor = gdat(params_equilibrium.shot,'b_field_tor','machine',gdat_params.machine);
+% b_field_z = gdat(params_equilibrium.shot,'b_field_z','machine',gdat_params.machine);
+% b_r = gdat(params_equilibrium.shot,'b_r','machine',gdat_params.machine);
+% b_tor = gdat(params_equilibrium.shot,'b_tor','machine',gdat_params.machine);
+% b_z = gdat(params_equilibrium.shot,'b_z','machine',gdat_params.machine);
+% grid = gdat(params_equilibrium.shot,'grid','machine',gdat_params.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_equilibrium.shot,'j_parallel','machine',machine);
-% j_tor = gdat(params_equilibrium.shot,'j_tor','machine',machine);
-% phi = gdat(params_equilibrium.shot,'phi','machine',machine);
-profiles_2d.psi = gdat(params_equilibrium.shot,'psi','machine',machine); % add psi_bound in a second step in special cases
-% r = gdat(params_equilibrium.shot,'r','machine',machine); % not to be filled since in grid.dim1
-% theta = gdat(params_equilibrium.shot,'theta','machine',machine);
-% z = gdat(params_equilibrium.shot,'z','machine',machine); % not to be filled since in grid.dim2
+% j_parallel = gdat(params_equilibrium.shot,'j_parallel','machine',gdat_params.machine);
+% j_tor = gdat(params_equilibrium.shot,'j_tor','machine',gdat_params.machine);
+% phi = gdat(params_equilibrium.shot,'phi','machine',gdat_params.machine);
+profiles_2d.psi = gdat(params_equilibrium.shot,'psi','machine',gdat_params.machine); % add psi_bound in a second step in special cases
+% r = gdat(params_equilibrium.shot,'r','machine',gdat_params.machine); % not to be filled since in grid.dim1
+% theta = gdat(params_equilibrium.shot,'theta','machine',gdat_params.machine);
+% z = gdat(params_equilibrium.shot,'z','machine',gdat_params.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
diff --git a/crpptbx/TCV_IMAS/tcv_ids_bpol_probe.m b/crpptbx/TCV_IMAS/tcv_ids_bpol_probe.m
index 6a00e2b6..e0aede17 100644
--- a/crpptbx/TCV_IMAS/tcv_ids_bpol_probe.m
+++ b/crpptbx/TCV_IMAS/tcv_ids_bpol_probe.m
@@ -15,40 +15,59 @@ if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_par
error_bar = gdat_params.error_bar;
end
-% Get data
+% Get data
tmp = gdat_tcv( shot, '\MAGNETICS::BPOL_003');
names = tmp.dim{2};
time = tmp.dim{1};
data = tmp.data;
ang_m = gdat_tcv( shot, 'static("ANG_M")');
+% from mapping file from Holger:
+probes_name={ '001' '002' '003' '004' '005' '006' '007' '008' '009' '010' '011' '012' '013' '014' '015' '016' '017' '018' '019' ...
+ '020' '021' '022' '023' '024' '025' '026' '027' '028' '029' '030' '031' '032' '033' '034' '035' '036' '037' '038'};
+probes_area = [9.200040D-03, 9.126160D-03, 9.163320D-03, 9.158160D-03, 9.147700D-03, 9.141400D-03, 9.155220D-03, 9.160560D-03, 9.058100D-03, ...
+ 9.228630D-03, 9.251260D-03, 9.149670D-03, 9.271000D-03, 9.112610D-03, 9.114790D-03, 9.180870D-03, 9.131410D-03, 9.120490D-03, ...
+ 9.154670D-03, 9.115670D-03, 9.190590D-03, 9.176150D-03, 9.186540D-03, 9.207520D-03, 9.176340D-03, 9.085060D-03, 9.155950D-03, ...
+ 9.135940D-03, 9.112780D-03, 9.059920D-03, 9.316970D-03, 9.114470D-03, 9.132330D-03, 9.086280D-03, 9.122720D-03, 9.045430D-03, ...
+ 9.106670D-03, 9.124440D-03];
+probes_length = [2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, ...
+ 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, ...
+ 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, ...
+ 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, 2.400000D-02, ...
+ 2.400000D-02, 2.400000D-02];
+
% Preallocate output structure
Nprobes = size(tmp.data,2);
sub_ids_struct_out(1:Nprobes) = ids_structures;
for ii=1:Nprobes
- sub_ids_struct_out{ii}.name = [names{ii}];
+ sub_ids_struct_out{ii}.name = [names{ii}];
sub_ids_struct_out_description{ii}.name = [names{ii}];
sub_ids_struct_out{ii}.position.r = mdsvalue('STATIC("R_M" )[$1]',sub_ids_struct_out{ii}.name);
sub_ids_struct_out_description{ii}.position_r = ['from ''STATIC("R_M" )[' sub_ids_struct_out{ii}.name ']'''];
- sub_ids_struct_out{ii}.position.z = mdsvalue('STATIC("Z_M" )[$1]',sub_ids_struct_out{ii}.name);
+ sub_ids_struct_out{ii}.position.z = mdsvalue('STATIC("Z_M" )[$1]',sub_ids_struct_out{ii}.name);
sub_ids_struct_out_description{ii}.position_z = ['from ''STATIC("Z_M" )[' sub_ids_struct_out{ii}.name ']'''];
sub_ids_struct_out{ii}.field.data = data(:,ii);
sub_ids_struct_out_description{ii}.field_data = 'from ''\MAGNETICS::BPOL_003''';
- sub_ids_struct_out{ii}.field.time = time;
-
+ sub_ids_struct_out{ii}.field.time = time;
+
sub_ids_struct_out{ii}.toroidal_angle = 0.; % to see if should match sector 3 (bpol003)
- aa = gdat_tcv( shot, ['static("ANG_M")[$1]'',''' sub_ids_struct_out{ii}.name '']);
+ aa = gdat_tcv( shot, ['static("ANG_M")[$1]'',''' sub_ids_struct_out{ii}.name '']);
sub_ids_struct_out{ii}.poloidal_angle = aa.data;
+ ij=strmatch(names{ii},probes_name);
+ sub_ids_struct_out{ii}.area = probes_area(ij);
+ sub_ids_struct_out{ii}.length = probes_length(ij);
sub_ids_struct_out_description{ii}.poloidal_angle = ['from ' aa.data_fullpath];
- sub_ids_struct_out_description{ii}.toroidal_angle = 'set to 0 although BPOL_003 means sector 3';
+ sub_ids_struct_out_description{ii}.toroidal_angle = 'set to 0';
+ sub_ids_struct_out_description{ii}.area = 'from array in machine description file';
+ sub_ids_struct_out_description{ii}.length = 'from array in machine description file';
end
fixed_error = 0.009999999776483;
switch error_bar
case 'delta'
for ii=1:Nprobes
- sub_ids_struct_out{ii}.field.data_error_upper = fixed_error.*ones(size(sub_ids_struct_out{ii}.field.data));% TO BE FOUND;
+ sub_ids_struct_out{ii}.field.data_error_upper = fixed_error.*ones(size(sub_ids_struct_out{ii}.field.data));% TO BE FOUND;
sub_ids_struct_out_description{ii}.field_data_error_upper = ['from fixed error value in case ' error_bar];
sub_ids_struct_out_description{ii}.field_data_error_lower = ['not provided since symmetric'];
end
@@ -60,7 +79,7 @@ switch error_bar
sub_ids_struct_out_description{ii}.field_data_error_lower = ['from fixed error value in case ' error_bar];
end
case 'added'
- for ii=1:Nprobes
+ for ii=1:Nprobes
sub_ids_struct_out{ii}.field.data_error_upper = sub_ids_struct_out{ii}.field.data + fixed_error.*ones(size(sub_ids_struct_out{ii}.field.data));
sub_ids_struct_out{ii}.field.data_error_lower = sub_ids_struct_out{ii}.field.data - fixed_error.*ones(size(sub_ids_struct_out{ii}.field.data));
sub_ids_struct_out_description{ii}.field_data_error_upper = ['from data + fixed error value in case ' error_bar];
diff --git a/crpptbx/TCV_IMAS/tcv_ids_flux_loop.m b/crpptbx/TCV_IMAS/tcv_ids_flux_loop.m
index 786c6b7a..f5fe9e9d 100644
--- a/crpptbx/TCV_IMAS/tcv_ids_flux_loop.m
+++ b/crpptbx/TCV_IMAS/tcv_ids_flux_loop.m
@@ -1,8 +1,8 @@
function [sub_ids_struct_out,sub_ids_struct_out_description] = tcv_ids_flux_loop(shot, ids_structures, gdat_params, varargin)
%
% [sub_ids_struct_out,sub_ids_struct_out_description] = tcv_ids_flux_loop(shot, ids_structures, gdat_params, varargin)
-%
-% Get ids field magnetics.flux_loop
+%
+% Get ids field magnetics.flux_loop
%
% error_bar: from gdat_params.error_bar
% 'delta' (default): error_bar to be added inserted in "upper" only as mentioned in description
@@ -15,7 +15,7 @@ if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_par
error_bar = gdat_params.error_bar;
end
-% Get data
+% Get data
tmp = gdat_tcv(shot, 'tcv_idealloop("FLUX")');
names = tmp.dim{2};
time = tmp.dim{1};
@@ -33,6 +33,8 @@ for ii=1:Nprobes
sub_ids_struct_out_description{ii}.position_r = ['from ''STATIC("R_F" )[' sub_ids_struct_out{ii}.name ']'''];
sub_ids_struct_out{ii}.position{1}.z = mdsvalue('STATIC("Z_F" )[$1]',sub_ids_struct_out{ii}.name);
sub_ids_struct_out_description{ii}.position_z = ['from ''STATIC("Z_F" )[' sub_ids_struct_out{ii}.name ']'''];
+ sub_ids_struct_out{ii}.position{1}.phi = 0.;
+ sub_ids_struct_out_description{ii}.position_phi = ['assumed 0'];
sub_ids_struct_out{ii}.flux.data = data(:,ii);
sub_ids_struct_out_description{ii}.flux_data = 'from ''tcv_idealloop("FLUX")''';
sub_ids_struct_out{ii}.flux.time = time;
@@ -42,14 +44,14 @@ end
fixed_error = 0.001200000056997;
switch error_bar
case 'delta'
- for ii=1:Nprobes
+ for ii=1:Nprobes
sub_ids_struct_out{ii}.flux.data_error_upper = fixed_error.*ones(size(sub_ids_struct_out{ii}.flux.data));
sub_ids_struct_out_description{ii}.flux_data_error_upper = ['from fixed error value in case ' error_bar];
sub_ids_struct_out_description{ii}.flux_data_error_lower = ['not provided since symmetric'];
%(not filled if symmetric) sub_ids_struct_out{ii}.flux.data_error_lower = 0.0012;%.*ones(size(sub_ids_struct_out{ii}.flux.data));
end
case 'delta_with_lower'
- for ii=1:Nprobes
+ for ii=1:Nprobes
sub_ids_struct_out{ii}.flux.data_error_upper = fixed_error.*ones(size(sub_ids_struct_out{ii}.flux.data));
sub_ids_struct_out{ii}.flux.data_error_lower = sub_ids_struct_out{ii}.flux.data_error_upper;
sub_ids_struct_out_description{ii}.flux_data_error_upper = ['from fixed error value in case ' error_bar];
@@ -57,7 +59,7 @@ switch error_bar
%(not filled if symmetric) sub_ids_struct_out{ii}.flux.data_error_lower = 0.0012;%.*ones(size(sub_ids_struct_out{ii}.flux.data));
end
case 'added'
- for ii=1:Nprobes
+ for ii=1:Nprobes
sub_ids_struct_out{ii}.flux.data_error_upper = sub_ids_struct_out{ii}.flux.data + fixed_error.*ones(size(sub_ids_struct_out{ii}.flux.data));
sub_ids_struct_out{ii}.flux.data_error_lower = sub_ids_struct_out{ii}.flux.data - fixed_error.*ones(size(sub_ids_struct_out{ii}.flux.data));
sub_ids_struct_out_description{ii}.flux.data_error_upper = ['from data + fixed error value in case ' error_bar];
@@ -66,4 +68,3 @@ switch error_bar
otherwise
error(['tcv_ids_flux_loop: error_bar option not known: ' error_bar])
end
-
--
GitLab