From f0eb3ddcb9ec550ba42ca5cdf648500126cf1610 Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Mon, 20 Jan 2020 22:50:16 +0100
Subject: [PATCH 1/6] Make labels more readable when plotting connection
matrix.
---
matlab/TCV_IMAS/plot_connection_matrices.m | 4 +++-
1 file changed, 3 insertions(+), 1 deletion(-)
diff --git a/matlab/TCV_IMAS/plot_connection_matrices.m b/matlab/TCV_IMAS/plot_connection_matrices.m
index 0d452f90..d575d048 100644
--- a/matlab/TCV_IMAS/plot_connection_matrices.m
+++ b/matlab/TCV_IMAS/plot_connection_matrices.m
@@ -44,6 +44,8 @@ for kk=1:info_pf_active.ntotcircuits
xlab{index} = [addcolor info_pf_active.coil{ii}.name 'out'];
end
+ xlab = strrep(xlab,'_','\_'); % since we need the tex interpreter for colors
+
shg
axis ij
ax = gca;
@@ -51,6 +53,6 @@ for kk=1:info_pf_active.ntotcircuits
xlabel('Element name. (red) Elemement belonging to circuit. (green) power supplies. (black) coils');
ylabel('Node');
title(['Circuit ' num2str(kk) ': ' info_pf_active.circuit{kk}.name])
- set(ax,'Xtick', xti, 'Ytick', yti, 'XTickLabel', xlab, 'YTickLabel', ylab')
+ set(ax,'Xtick', xti, 'Ytick', yti, 'XTickLabel', xlab, 'YTickLabel', ylab', 'XTickLabelRotation', 90)
end
--
GitLab
From 14319f15c54bc9dbfe91205188a121e200c72d9d Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Tue, 21 Jan 2020 10:36:02 +0100
Subject: [PATCH 2/6] Get probes area from MDS. Get all probes location in one
go and ensure order is preserved.
---
matlab/TCV_IMAS/tcv_ids_bpol_probe.m | 52 +++++++++++-----------------
1 file changed, 21 insertions(+), 31 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_ids_bpol_probe.m b/matlab/TCV_IMAS/tcv_ids_bpol_probe.m
index 2cb28fae..453c94b9 100644
--- a/matlab/TCV_IMAS/tcv_ids_bpol_probe.m
+++ b/matlab/TCV_IMAS/tcv_ids_bpol_probe.m
@@ -25,24 +25,17 @@ tmp = gdat_tcv( shot,params_eff);
names = tmp.dim{2};
time = tmp.dim{1};
data = tmp.data;
-params_eff.data_request = 'static("ANG_M")';
-% ang_m = gdat_tcv(shot,params_eff);
-% ang_m_desc = params_eff.data_request;
+ang_m = mdsvalue('static("ANG_M")[$1]',names);
+ang_m_desc = 'from ''static("ANG_M")''';
+r_m = mdsvalue('static("R_M")[$1]',names);
+r_m_desc = 'from ''static("R_M")''';
+z_m = mdsvalue('static("Z_M")[$1]',names);
+z_m_desc = 'from ''static("Z_M")''';
+area = mdsvalue('\magnetics::bpol_003:area');
+area_desc = 'from ''\magnetics::bpol_003_area''';
-%TODO find this either in the static or the magnetics tree.
% 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];
+probes_length = 2.400000D-02;
% Preallocate output structure
Nprobes = size(tmp.data,2);
@@ -50,26 +43,23 @@ sub_ids_struct_out(1:Nprobes) = ids_structures;
sub_ids_struct_out_description = cell(1,Nprobes);
for ii=1:Nprobes
- 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_description{ii}.position_z = ['from ''STATIC("Z_M" )[' sub_ids_struct_out{ii}.name ']'''];
+ sub_ids_struct_out{ii}.name = names{ii};
+ sub_ids_struct_out_description{ii}.name = names{ii};
+ sub_ids_struct_out{ii}.position.r = r_m(ii);
+ sub_ids_struct_out_description{ii}.position_r = r_m_desc;
+ sub_ids_struct_out{ii}.position.z = z_m(ii);
+ sub_ids_struct_out_description{ii}.position_z = z_m_desc;
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}.toroidal_angle = 0.; % to see if should match sector 3 (bpol003)
- params_eff.data_request = ['static("ANG_M")[$1]'',''' sub_ids_struct_out{ii}.name ''];
- aa = gdat_tcv(shot,params_eff);
- sub_ids_struct_out{ii}.poloidal_angle = -aa.data;
- ij=strcmp(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{ii}.toroidal_angle = 0.;
+ sub_ids_struct_out{ii}.poloidal_angle = -ang_m(ii);
+ sub_ids_struct_out{ii}.area = area(ii);
+ sub_ids_struct_out{ii}.length = probes_length;
+ sub_ids_struct_out_description{ii}.poloidal_angle = ang_m_desc;
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}.area = area_desc;
sub_ids_struct_out_description{ii}.length = 'from array in machine description file';
end
--
GitLab
From 4c609bcb1a861b5fc3bdeb3ec769370bbe3afbbb Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Wed, 22 Jan 2020 15:59:29 +0100
Subject: [PATCH 3/6] Fix typo.
---
matlab/TCV_IMAS/tcv_get_ids_tf.m | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_tf.m b/matlab/TCV_IMAS/tcv_get_ids_tf.m
index 9bb06703..c8605e17 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_tf.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_tf.m
@@ -26,7 +26,7 @@ if ~ischar(tmp.data)
ids_tf_description.r0 = ' b0.r0 from gdat_tcv(shot, ''b0'')';
ids_tf.b_field_tor_vacuum_r.data = tmp.data .* ids_tf.r0;
ids_tf_description.b_field_tor_vacuum_r.data = ' gdat_tcv(shot, ''b0'') * .r0';
- ids_tf.b_field_tor_vaccum_r.data_error_upper = 19.2e-6*200; %TODO: Check the 200A value (might be more since TOR goes up to 70kA)
+ ids_tf.b_field_tor_vacuum_r.data_error_upper = 19.2e-6*200; %TODO: Check the 200A value (might be more since TOR goes up to 70kA)
ids_tf_description.b_field_tor_vacuum_r.data_error_upper = ' Error on coil measurement is 200A, b0.r0 is 19.2e-6*iphi';
ids_tf.b_field_tor_vacuum_r.time = tmp.t;
else
--
GitLab
From e4b796368ee2a66d0cde4248d9ce8051b258a333 Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Wed, 22 Jan 2020 16:35:03 +0100
Subject: [PATCH 4/6] Fix incomplete second profiles_1d/ion substructure.
---
matlab/TCV_IMAS/tcv_get_ids_core_profiles.m | 13 ++++++++-----
1 file changed, 8 insertions(+), 5 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
index c6db74e2..9f37ddb0 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_core_profiles.m
@@ -1,4 +1,4 @@
-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);
%
@@ -13,7 +13,7 @@ function [ids_core_profiles,ids_core_profiles_description,varargout] = tcv_get_i
%
error_bar = 'delta';
-if exist('gdat_params') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar)
+if exist('gdat_params','var') && isfield(gdat_params,'error_bar') && ~isempty(gdat_params.error_bar)
error_bar = gdat_params.error_bar;
end
@@ -21,7 +21,7 @@ machine = 'tcv';
tens_time = -1;
tens_rho = -0.1;
-if exist('gdat_params')
+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{:});
else
[ids_core_profiles, params_cores_profiles] = tcv_ids_headpart(shot,ids_equil_empty,'cores_profiles',varargin{:});
@@ -29,7 +29,7 @@ else
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)
+try params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def)
ids_core_profiles_description = [];
@@ -262,6 +262,9 @@ it_ti = iround_os(temp_1d.ti.fit.t,ids_core_profiles.time);
% assumed 1 impurity with Zp=6
Zp = 6.;
for it=1:length(ids_core_profiles.time)
+ % Duplicate ion substructure
+ ids_core_profiles.profiles_1d{it}.ion(1:2) = ids_core_profiles.profiles_1d{it}.ion(1);
+
ids_core_profiles.profiles_1d{it}.ion{1}.element{1}.a = 2.;
ids_core_profiles.profiles_1d{it}.ion{1}.element{1}.z_n = 1;
ids_core_profiles.profiles_1d{it}.ion{1}.element{1}.atoms_n = 1;
@@ -351,7 +354,7 @@ for it=1:length(ids_core_profiles.time)
end
% cocos automatic transform
-if exist('ids_generic_cocos_nodes_transformation_symbolic') == 2
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
[ids_core_profiles,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_core_profiles,'core_profiles',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);
--
GitLab
From 1532d2a8ed4827943a394297b92065e331e1a822 Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Fri, 24 Jan 2020 12:45:29 +0100
Subject: [PATCH 5/6] Fix bug when computing G-coil geometry.
---
matlab/TCV_IMAS/tcv_ids_coil.m | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/matlab/TCV_IMAS/tcv_ids_coil.m b/matlab/TCV_IMAS/tcv_ids_coil.m
index 8fc6ae81..9f2bda74 100644
--- a/matlab/TCV_IMAS/tcv_ids_coil.m
+++ b/matlab/TCV_IMAS/tcv_ids_coil.m
@@ -60,7 +60,7 @@ iT = strcmp('T_003',namepfc); % Return current
nt_c(iT)=1;
% Approximate circular G-coil conductor with a square of the same crosssection
-iG=strcmp('G_00',namepfc);
+iG=strncmp('G_00',namepfc,4);
w_c(iG)=sqrt(xsect_c(iG));
h_c(iG)=sqrt(xsect_c(iG));
--
GitLab
From 4dffbd0fd036b0c6b7081c87062c72476d091a6e Mon Sep 17 00:00:00 2001
From: Antoine Merle <antoine.merle@epfl.ch>
Date: Fri, 24 Jan 2020 12:46:34 +0100
Subject: [PATCH 6/6] Add equilibrium constraints when fetching LIUQE data.
---
matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 157 ++++++++++++++++++----
matlab/TCV_IMAS/tcv_get_ids_magnetics.m | 2 +-
2 files changed, 129 insertions(+), 30 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index 391d8e44..bb275385 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -1,4 +1,4 @@
-function [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_equilibrium(shot,ids_equil_empty, gdat_params,varargin);
+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);
%
@@ -6,7 +6,7 @@ function [ids_equilibrium,ids_equilibrium_description,varargout] = tcv_get_ids_e
% gdat_params: gdat_data.gdat_params to get all params passed from original call, in particular error_bar and cocos_out options
%
-if exist('gdat_params')
+if exist('gdat_params','var')
[ids_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium','gdat_params',gdat_params,varargin{:});
else
[ids_equilibrium, params_equilibrium] = tcv_ids_headpart(shot,ids_equil_empty,'equilibrium',varargin{:});
@@ -22,14 +22,14 @@ end
%
%% liuqe.m at this stage is missing correction 0.996, so use std source by time of default liuqe to make sure
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)
+try params_eff_ref=rmfield(params_eff_ref,'source');catch;end % make sure no source (from ids def)
params_eff = params_eff_ref;
params_eff.data_request='b0'; params_eff.source='liuqe'; % to get liuqe time array
bb = gdat(params_equilibrium.shot,params_eff);
params_eff = rmfield(params_eff,'source'); % to get original magnetics data
vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,params_eff);
-ij_ok = [isfinite(vacuum_toroidal_field.b0.data)];
+ij_ok = isfinite(vacuum_toroidal_field.b0.data);
vacuum_toroidal_field.b0.data = interpos(21,vacuum_toroidal_field.b0.t(ij_ok),vacuum_toroidal_field.b0.data(ij_ok),bb.t);
vacuum_toroidal_field.b0.t = bb.t;
vacuum_toroidal_field.b0.dim = {vacuum_toroidal_field.b0.t};
@@ -41,15 +41,16 @@ ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
ids_equilibrium.time = vacuum_toroidal_field.b0.t;
ids_equilibrium_description.time = '.t subfield from: [''b0'',''source'',''liuqe'']';
+ntime = numel(ids_equilibrium.time);
-ids_equilibrium.time_slice(1:numel(ids_equilibrium.time)) = ids_equilibrium.time_slice(1);
+ids_equilibrium.time_slice(1:ntime) = ids_equilibrium.time_slice(1);
% load time array data to copy to time_slices
% 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
% brute force solution load all eqdsks
-% $$$ for it=1:numel(ids_equilibrium.time)
+% $$$ for it=1:ntime
% $$$ ids_equilibrium.time(it)
% $$$ temp.eqdsks{it}=gdat(params_equilibrium.shot,'eqdsk','time',ids_equilibrium.time(it),'write',0,'machine',gdat_params.machine);
% $$$ end
@@ -109,7 +110,7 @@ global_quantities_desc.w_mhd = params_eff.data_request;
global_quantities_fieldnames = fieldnames(global_quantities);
special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
for i=1:numel(global_quantities_fieldnames)
if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}))
@@ -123,7 +124,7 @@ for it=1:numel(ids_equilibrium.time)
end
% special case
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.r = temp.r_magnetic_axis.data(it);
ids_equilibrium.time_slice{it}.global_quantities.magnetic_axis.z = temp.z_magnetic_axis.data(it);
ids_equilibrium.time_slice{it}.global_quantities.psi_axis = temp.psi_axis.data(it) + ...
@@ -181,7 +182,7 @@ temp_desc.z_lcfs = params_eff.data_request;
boundary_fieldnames = fieldnames(boundary);
special_fields = {'lcfs', 'geometric_axis', 'x_point'}; % fields needing non-automatic treatments
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
for i=1:numel(boundary_fieldnames)
if ~any(strcmp(boundary_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}))
@@ -195,7 +196,7 @@ for it=1:numel(ids_equilibrium.time)
end
% special cases
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
ids_equilibrium.time_slice{it}.boundary.outline.r = temp.r_lcfs.data(:,it);
ids_equilibrium.time_slice{it}.boundary.outline.z = temp.z_lcfs.data(:,it);
ids_equilibrium.time_slice{it}.boundary.lcfs.r = ids_equilibrium.time_slice{it}.boundary.outline.r;
@@ -214,6 +215,115 @@ for it=1:numel(ids_equilibrium.time)
end
end
+%% constraints
+% TODO: Add description
+
+% Measured values
+liuqe_time = ids_equilibrium.time; % Not necessarily magnetics time so far
+mag_time = mdsvalue('\magnetics::bpol_003:dim0');
+itime = iround_os(mag_time, liuqe_time);
+mag_time_req = mdscvt(mag_time(itime),'f');
+bpol = mdsvalue('\magnetics::bpol_003[$1,*]',mag_time_req);
+flux = mdsvalue('\magnetics::flux[$1,*]',mag_time_req);
+diam = mdsvalue('\results::dmlcor[$1]',mag_time_req);
+ip = mdsvalue('\magnetics::iplasma:trapeze[$1]',mag_time_req);
+% Coil currents since dim of constraints pf_current is IDS:pf_active/coil
+dim_pol = {'OH_001','OH_002','OH_002','OH_002','OH_002','OH_002','OH_002',...
+ 'E_001','E_002','E_003','E_004','E_005','E_006','E_007','E_008',...
+ 'F_001','F_002','F_003','F_004','F_005','F_006','F_007','F_008',...
+ 'G_001','G_001','G_001','G_001','G_001','G_001'};
+ipol = mdsvalue('\magnetics::ipol[$1,$2]',mag_time_req,dim_pol);
+dim_pol(30:32) = {'TOR_001'};
+ipol(:,30:32) = 0; % TOR_001 is not used in LIUQE
+
+% Reconstructed values
+ipol_liuqe_order = [18,19*ones(1,6),1:16,17*ones(1,6)]; % LIUQE order is E F G OH
+bpol_liuqe = mdsvalue('tcv_eq("b_probe","liuqe.m")');
+flux_liuqe = mdsvalue('tcv_eq("psi_loop","liuqe.m")');
+diam_liuqe = mdsvalue('tcv_eq("tor_flux_dml","liuqe.m")');
+ip_liuqe = mdsvalue('tcv_eq("i_pl","liuqe.m")');
+ipol_liuqe = mdsvalue('tcv_eq("i_pol","liuqe.m")');
+ipol_liuqe = ipol_liuqe(ipol_liuqe_order,:);
+ipol_liuqe(20,:) = 0; % ... TOR
+
+% Weights (using old parameters tree for now)
+bpol_err = mdsvalue('\results::parameters:berr')./mdsvalue('\results::parameters:vvv[0:37]');
+flux_err = mdsvalue('\results::parameters:ferr')./mdsvalue('\results::parameters:www[0:37]')*2*pi;
+diam_err = 0.13e-3./mdsvalue('\results::parameters:idml');
+ip_err = mdsvalue('\results::parameters:plcerr')*1e3;
+ipol_err = mdsvalue('\results::parameters:cerr')./mdsvalue('\results::parameters:uuu[0:18]')*1e3;
+ipol_err = ipol_err(ipol_liuqe_order);
+ipol_err(20) = NaN;
+
+if ntime > 0
+ constraints_orig = ids_equilibrium.time_slice{1}.constraints;
+ % Remove unused arrays
+ ununsed_constraints = {'faraday_angle','mse_polarisation_angle','iron_core_segment',...
+ 'n_e','n_e_line','pressure','q','x_point'};
+ for name = ununsed_constraints, constraints_orig.(name{1})={}; end
+end
+for it = 1:ntime
+ constraints = constraints_orig;
+ % bpol_probe
+ nbpol = size(bpol,2);
+ bpol_probe(1:nbpol) = constraints.bpol_probe(1);
+ for ib = 1:nbpol
+ bpol_probe{ib}.measured = bpol(it,ib);
+ bpol_probe{ib}.source = sprintf('IDS:magnetics/bpol_probe[%02d]/field',ib);
+ bpol_probe{ib}.time_measurement = mag_time(itime(it));
+ bpol_probe{ib}.exact = 0;
+ bpol_probe{ib}.weight = 1/(bpol_err(ib)).^2;
+ bpol_probe{ib}.reconstructed = bpol_liuqe(ib,it);
+ end
+ constraints.bpol_probe = bpol_probe;
+ % flux_loop
+ nflux = size(flux,2);
+ flux_loop(1:nflux) = constraints.flux_loop(1);
+ for il = 1:nflux
+ flux_loop{il}.measured = flux(it,il);
+ flux_loop{il}.source = sprintf('IDS:magnetics/flux_loop[%02d]/flux',il);
+ flux_loop{il}.time_measurement = mag_time(itime(it));
+ flux_loop{il}.exact = 0;
+ flux_loop{il}.weight = 1/(flux_err(il)).^2;
+ flux_loop{il}.reconstructed = flux_liuqe(il,it);
+ end
+ constraints.flux_loop = flux_loop;
+ % ip
+ constraints.ip.measured = ip(it);
+ constraints.ip.source = 'IDS:magnetics/method[1]/ip';
+ constraints.ip.time_measurement = mag_time(itime(it));
+ constraints.ip.exact = 0;
+ constraints.ip.weight = 1/(ip_err).^2;
+ constraints.ip.reconstructed = ip_liuqe(it);
+ % diamagnetic_flux
+ constraints.diamagnetic_flux.measured = diam(it);
+ constraints.diamagnetic_flux.source = 'IDS:magnetics/method[1]/diamagnetic_flux';
+ constraints.diamagnetic_flux.time_measurement = mag_time(itime(it));
+ constraints.diamagnetic_flux.exact = 0;
+ constraints.diamagnetic_flux.weight = 1/(diam_err).^2;
+ constraints.diamagnetic_flux.reconstructed = diam_liuqe(it);
+ % pf_current
+ nipol = size(ipol,2);
+ pf_current(1:nipol) = constraints.pf_current(1);
+ for ic = 1:nipol
+ pf_current{ic}.measured = ipol(it,ic);
+ pf_current{ic}.source = sprintf('IDS:pf_active/coil[%02d]/current',ic);
+ pf_current{ic}.time_measurement = mag_time(itime(it));
+ if strcmp(dim_pol{ic},'TOR_001')
+ pf_current{ic}.source = [pf_current{ic}.source,' replaced with 0'];
+ pf_current{ic}.exact = 1;
+ else
+ pf_current{ic}.exact = 0;
+ pf_current{ic}.weight = 1/(ipol_err(ic)).^2;
+ pf_current{ic}.reconstructed = ipol_liuqe(ic,it);
+ end
+ end
+ constraints.pf_current = pf_current;
+
+ ids_equilibrium.time_slice{it}.constraints = constraints;
+end
+
+
%
%% profiles_1d (cannot use eqdsk since not same radial mesh)
%
@@ -314,7 +424,9 @@ profiles_1d.gm9.data = tmp_gm.x;
profiles_1d_desc.gm9 = 'FS_av(1./R.^1)';
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)
+nrho = numel(profiles_1d.rho_tor.x);
+gradpsi_over_R_sq = NaN(nrho,ntime);
+for it=1:ntime
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
@@ -335,7 +447,7 @@ profiles_1d.j_parallel.data = j_par./repmat(ids_equilibrium.vacuum_toroidal_fiel
profiles_1d_fieldnames = fieldnames(profiles_1d);
special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
for i=1:numel(profiles_1d_fieldnames)
if ~any(strcmp(profiles_1d_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}))
@@ -352,9 +464,7 @@ for it=1:numel(ids_equilibrium.time)
end
% special cases
-nrho = numel(profiles_1d.rho_tor.x);
-ntime = numel(temp.psi_axis.data);
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
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;
@@ -390,7 +500,7 @@ profiles_2d_desc.psi = [params_eff.data_request ' adding psi_bound in a 2nd step
profiles_2d_fieldnames = fieldnames(profiles_2d);
special_fields = {'grid', 'grid_type'}; % fields needing non-automatic treatments
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
for i=1:numel(profiles_2d_fieldnames)
if ~any(strcmp(profiles_2d_fieldnames{i},special_fields))
if ~isstruct(ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}))
@@ -407,7 +517,7 @@ for it=1:numel(ids_equilibrium.time)
end
% special cases
-for it=1:numel(ids_equilibrium.time)
+for it=1:ntime
ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.name = profiles_2d.grid_type.name;
ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.index = profiles_2d.grid_type.index;
ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.description = profiles_2d.grid_type.description;
@@ -424,17 +534,6 @@ for it=1:numel(ids_equilibrium.time_slice)
ids_equilibrium.time_slice{it}.boundary.strike_point = {};
ids_equilibrium.time_slice{it}.boundary_separatrix.x_point = {};
ids_equilibrium.time_slice{it}.boundary_separatrix.strike_point = {};
- ids_equilibrium.time_slice{it}.constraints.bpol_probe = {};
- ids_equilibrium.time_slice{it}.constraints.faraday_angle = {};
- ids_equilibrium.time_slice{it}.constraints.mse_polarisation_angle = {};
- ids_equilibrium.time_slice{it}.constraints.flux_loop = {};
- ids_equilibrium.time_slice{it}.constraints.iron_core_segment = {};
- ids_equilibrium.time_slice{it}.constraints.n_e = {};
- ids_equilibrium.time_slice{it}.constraints.n_e_line = {};
- ids_equilibrium.time_slice{it}.constraints.pf_current = {};
- ids_equilibrium.time_slice{it}.constraints.pressure = {};
- ids_equilibrium.time_slice{it}.constraints.q = {};
- ids_equilibrium.time_slice{it}.constraints.x_point = {};
end
% special test matrix cocos transform
@@ -465,7 +564,7 @@ end
% $$$ end
% cocos automatic transform
-if exist('ids_generic_cocos_nodes_transformation_symbolic') == 2
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
[ids_equilibrium,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_equilibrium,'equilibrium',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);
diff --git a/matlab/TCV_IMAS/tcv_get_ids_magnetics.m b/matlab/TCV_IMAS/tcv_get_ids_magnetics.m
index 8f405901..b20f74f1 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_magnetics.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_magnetics.m
@@ -24,7 +24,7 @@ ids_magnetics.b_field_tor_probe = {};
ids_magnetics.rogowski_coil = {};
% cocos automatic transform
-if exist('ids_generic_cocos_nodes_transformation_symbolic','file')
+if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic'))
[ids_magnetics,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_magnetics,'magnetics',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);
--
GitLab