From 0c4c8438ae7528716d1706968544ef5035cc1d7d Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Mon, 31 Jan 2022 15:20:23 +0100
Subject: [PATCH] add IDS/IMAS to eqdsk directly
tested: gdat(131025,'eqdsk','run',43,'database_user','public','tokamak_name','ITER','machine','IMAS','time',208);
---
matlab/IMAS/ITERlimiter.data | 346 +++++++++++++++++++++++++++++++++++
matlab/IMAS/gdat_imas.m | 27 +++
matlab/IMAS/ids2eqdsk.m | 256 ++++++++++++++++++++++++++
3 files changed, 629 insertions(+)
create mode 100644 matlab/IMAS/ITERlimiter.data
create mode 100644 matlab/IMAS/ids2eqdsk.m
diff --git a/matlab/IMAS/ITERlimiter.data b/matlab/IMAS/ITERlimiter.data
new file mode 100644
index 00000000..b4637a5d
--- /dev/null
+++ b/matlab/IMAS/ITERlimiter.data
@@ -0,0 +1,346 @@
+% R Z_wall_131025run43t208ITERpublic
+4.105300 -2.503700
+4.105300 -1.491400
+4.105300 -0.476100
+4.105300 0.540200
+4.105300 1.556600
+4.105300 2.571900
+4.125300 3.585200
+4.324600 4.337100
+4.935400 4.719600
+5.753400 4.534400
+6.552400 3.925600
+7.401500 3.179600
+7.906200 2.464700
+8.269700 1.684700
+8.393800 0.635400
+8.305700 -0.420000
+7.900200 -1.337200
+7.282400 -2.254400
+6.266100 -3.043400
+6.363200 -3.244600
+6.265100 -3.228600
+6.154900 -3.229600
+6.027800 -3.253600
+5.907600 -3.302700
+5.799500 -3.374800
+5.708400 -3.466900
+5.637300 -3.575100
+5.589200 -3.696200
+5.569200 -3.793300
+5.564200 -3.923500
+5.564200 -4.574300
+5.556600 -4.545000
+5.272300 -4.260800
+5.252600 -3.982400
+5.192900 -3.883600
+5.133200 -3.820900
+5.062800 -3.770200
+4.984300 -3.733500
+4.900400 -3.712000
+4.813900 -3.706400
+4.727900 -3.716900
+4.645300 -3.743200
+4.491500 -3.906500
+4.188200 -3.882600
+4.162400 -3.917500
+4.474500 -3.277500
+4.510400 -3.179600
+4.526000 -3.076600
+4.520500 -2.972500
+4.494300 -2.871600
+4.448300 -2.778100
+4.384400 -2.695800
+4.305300 -2.627900
+4.214100 -2.577400
+4.114600 -2.546200
+3.942100 -2.535700
+4.105300 -2.503700
+3.539600 -2.532100
+3.539600 -2.334900
+3.539600 -2.137800
+3.539600 -1.940700
+3.539600 -1.743500
+3.539600 -1.546400
+3.539600 -1.349300
+3.539600 -1.152100
+3.539600 -0.955000
+3.539600 -0.757900
+3.539600 -0.560700
+3.539600 -0.363600
+3.539600 -0.166500
+3.539600 0.030700
+3.539600 0.227800
+3.539600 0.425000
+3.539600 0.622100
+3.539600 0.819200
+3.539600 1.016400
+3.539600 1.213500
+3.539600 1.410600
+3.539600 1.607800
+3.539600 1.804900
+3.539600 2.002000
+3.539600 2.199200
+3.539600 2.396300
+3.539600 2.593500
+3.539600 2.790600
+3.539600 2.987700
+3.539600 3.184900
+3.539600 3.382000
+3.539600 3.579100
+3.551400 3.775300
+3.586900 3.968600
+3.645300 4.156200
+3.726000 4.335400
+3.827600 4.503500
+3.948900 4.658200
+4.087800 4.797100
+4.242500 4.918300
+4.410700 5.020000
+4.589900 5.100600
+4.777500 5.159000
+4.970800 5.194400
+5.167000 5.206200
+5.363100 5.194300
+5.556400 5.158900
+5.744000 5.100400
+5.923200 5.019700
+6.091300 4.918000
+6.247400 4.810200
+6.403400 4.702400
+6.559500 4.594600
+6.715600 4.486800
+6.871600 4.379000
+7.031500 4.264000
+7.186700 4.142900
+7.337100 4.015800
+7.482500 3.883000
+7.622600 3.744600
+7.757200 3.600900
+7.886100 3.452100
+8.009200 3.298400
+8.126200 3.140000
+8.237000 2.977200
+8.341300 2.810200
+8.439100 2.639300
+8.530200 2.464700
+8.614500 2.286700
+8.691800 2.105600
+8.761900 1.921600
+8.824900 1.735100
+8.876400 1.552900
+8.916600 1.367900
+8.945300 1.180900
+8.963800 1.007300
+8.976800 0.833300
+8.984400 0.659000
+8.986500 0.484500
+8.983200 0.310000
+8.974400 0.135700
+8.960200 -0.038200
+8.936400 -0.208800
+8.898500 -0.376800
+8.846800 -0.541100
+8.781700 -0.700600
+8.702300 -0.874300
+8.622800 -1.048000
+8.543300 -1.221700
+8.463800 -1.395400
+8.384400 -1.569100
+8.304900 -1.742800
+8.225400 -1.916600
+8.145900 -2.090300
+8.066500 -2.264000
+7.987000 -2.437700
+7.907500 -2.611400
+7.828000 -2.785100
+7.748600 -2.958800
+7.669100 -3.132500
+7.589600 -3.306200
+7.510200 -3.479900
+7.430700 -3.653600
+7.345800 -3.822200
+7.248600 -3.984000
+7.139600 -4.138000
+7.019500 -4.283500
+6.888800 -4.419600
+6.748300 -4.545600
+6.598800 -4.660800
+6.441100 -4.764500
+6.276200 -4.856100
+6.104800 -4.935200
+5.928100 -5.001300
+5.746900 -5.054000
+5.562300 -5.093100
+5.375300 -5.118300
+5.186900 -5.129500
+4.989700 -5.121700
+4.794900 -5.089900
+4.605400 -5.034500
+4.424200 -4.956400
+4.253800 -4.856700
+4.096900 -4.736900
+3.955800 -4.598900
+3.832700 -4.444600
+3.729300 -4.276500
+3.647200 -4.097000
+3.587700 -3.908800
+3.551600 -3.714800
+3.539600 -3.517800
+3.539600 -3.320600
+3.539600 -3.123500
+3.539600 -2.926400
+3.539600 -2.729200
+3.539600 -2.532100
+3.262200 -2.533300
+3.262200 -2.335900
+3.262200 -2.138600
+3.262200 -1.941200
+3.262200 -1.743900
+3.262200 -1.546600
+3.262200 -1.349200
+3.262200 -1.151900
+3.262200 -0.954500
+3.262200 -0.757200
+3.262200 -0.559900
+3.262200 -0.362500
+3.262200 -0.165200
+3.262200 0.032200
+3.262200 0.229500
+3.262200 0.426900
+3.262200 0.624200
+3.262200 0.821500
+3.262200 1.018900
+3.262200 1.216200
+3.262200 1.413600
+3.262200 1.610900
+3.262200 1.808300
+3.262200 2.005600
+3.262200 2.202900
+3.262200 2.400300
+3.262200 2.597600
+3.262200 2.795000
+3.262200 2.992300
+3.262200 3.189700
+3.262200 3.387000
+3.262200 3.584300
+3.262200 3.781700
+3.271700 3.969700
+3.300200 4.155900
+3.347300 4.338200
+3.412600 4.514800
+3.495300 4.683900
+3.594800 4.843800
+3.709800 4.992900
+3.839400 5.129500
+3.982000 5.252400
+4.136400 5.360300
+4.300800 5.452000
+4.473700 5.526600
+4.653200 5.583400
+4.827400 5.622300
+5.004000 5.648300
+5.182000 5.661400
+5.360500 5.661400
+5.538500 5.648400
+5.715100 5.622400
+5.889400 5.583600
+6.060300 5.532200
+6.227000 5.468400
+6.409000 5.387300
+6.588000 5.300000
+6.763900 5.206400
+6.936500 5.106700
+7.105400 5.001100
+7.270600 4.889700
+7.431800 4.772600
+7.588800 4.649900
+7.741400 4.521800
+7.889500 4.388500
+8.032800 4.250100
+8.171200 4.106800
+8.304600 3.958800
+8.432700 3.806200
+8.555400 3.649200
+8.672600 3.488100
+8.784000 3.322900
+8.889700 3.154000
+8.989400 2.981500
+9.083000 2.805600
+9.170500 2.626600
+9.251600 2.444700
+9.326400 2.260000
+9.394600 2.072800
+9.456300 1.883300
+9.511400 1.691900
+9.559800 1.498600
+9.601400 1.303700
+9.636100 1.107600
+9.661600 0.915300
+9.677100 0.722000
+9.682500 0.528200
+9.677700 0.334300
+9.662900 0.141000
+9.638000 -0.051300
+9.603100 -0.242100
+9.558300 -0.430800
+9.503800 -0.616900
+9.439600 -0.799900
+9.366000 -0.979300
+9.290000 -1.151400
+9.214100 -1.323500
+9.138200 -1.495600
+9.062300 -1.667700
+8.986300 -1.839800
+8.910400 -2.011900
+8.834500 -2.184000
+8.758500 -2.356100
+8.682600 -2.528200
+8.606700 -2.700300
+8.530800 -2.872400
+8.454800 -3.044500
+8.378900 -3.216600
+8.303000 -3.388700
+8.227000 -3.560800
+8.151100 -3.732900
+8.068100 -3.906600
+7.974800 -4.074900
+7.871400 -4.237300
+7.758400 -4.393100
+7.636200 -4.541800
+7.505200 -4.682800
+7.365900 -4.815700
+7.218800 -4.939800
+7.064500 -5.054900
+6.903500 -5.160400
+6.736400 -5.255900
+6.563800 -5.341200
+6.386400 -5.415800
+6.204800 -5.479600
+6.019600 -5.532300
+5.831700 -5.573700
+5.641500 -5.603600
+5.449900 -5.622000
+5.257500 -5.628700
+5.065100 -5.623800
+4.873400 -5.603400
+4.684800 -5.563700
+4.501100 -5.505200
+4.324300 -5.428500
+4.156100 -5.334300
+3.998200 -5.223700
+3.852300 -5.097700
+3.719900 -4.957600
+3.602200 -4.804900
+3.500500 -4.641100
+3.415900 -4.467900
+3.349200 -4.287100
+3.301000 -4.100400
+3.271900 -3.909900
+3.262200 -3.717300
+3.262200 -3.520000
+3.262200 -3.322600
+3.262200 -3.125300
+3.262200 -2.928000
+3.262200 -2.730600
+3.262200 -2.533300
diff --git a/matlab/IMAS/gdat_imas.m b/matlab/IMAS/gdat_imas.m
index 5978eb9f..24c3e6db 100644
--- a/matlab/IMAS/gdat_imas.m
+++ b/matlab/IMAS/gdat_imas.m
@@ -518,6 +518,33 @@ elseif strcmp(mapping_for_imas.method,'switchcase')
end
end
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'eqdsk'}
+ % add defaults
+ if ~isfield(gdat_data.gdat_params,'time')
+ gdat_data.gdat_params.time = 1.; % default time
+ end
+ if ~isfield(gdat_data.gdat_params,'write')
+ gdat_data.gdat_params.write = 1; % default write eqdsk to /tmp/username
+ end
+ params_eff = gdat_data.gdat_params;
+ params_eff = rmfield(params_eff,'time');
+ params_eff.data_request = 'ids';
+ params_eff.source = 'equilibrium';
+ gdat_equil = gdat_imas(shot,params_eff);
+ params_eff.source = 'wall';
+ gdat_wall = gdat_imas(shot,params_eff);
+ %
+ eqdsk_filename_suffix = [num2str(gdat_data.shot) '_' num2str(gdat_data.gdat_params.run) 't' num2str(gdat_data.gdat_params.time) ...
+ gdat_data.gdat_params.tokamak_name '_' gdat_data.gdat_params.database_user];
+ [eqdsk_out] = ids2eqdsk(gdat_equil.equilibrium,gdat_wall.wall,gdat_data.gdat_params.time,eqdsk_filename_suffix,gdat_data.gdat_params.write);
+ gdat_data.equilibrium = gdat_equil.equilibrium;
+ gdat_data.wall = gdat_wall.wall;
+ gdat_data.eqdsk = eqdsk_out;
+ if isfield(gdat_data.gdat_params,'doplot') && gdat_data.gdat_params.doplot == 1
+ plot_eqdsk(gdat_data.eqdsk)
+ end
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
otherwise
if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_imas']); end
diff --git a/matlab/IMAS/ids2eqdsk.m b/matlab/IMAS/ids2eqdsk.m
new file mode 100644
index 00000000..5f31a442
--- /dev/null
+++ b/matlab/IMAS/ids2eqdsk.m
@@ -0,0 +1,256 @@
+function [eqdsk_out,varargout] = ids2eqdsk(ids_equilibrium,ids_wall,time_ref,eqdsk_filename_suffix,write_opt)
+%
+% [eqdsk_out, varout] = ids2eqdsk(ids_equilibrium, ids_wall, time_ref, eqdsk_filename_suffix);
+%
+% ids_equilibrium: ids equilibrium structure to extract eqdsk related quantities
+% ids_wall: ids wall structure to extract limiter and vessel points if provided
+% time_ref: time to extract data
+% eqdsk_filename_suffix: string to be used for end of eqdsk filename, e.g. shot_run_text
+% write_opt: 1 (default) do write eqdsk to /tmp/username, 0: no writing
+%
+% OUTPUT
+% eqdsk_out: eqdsk structure
+%
+if ~exist('time_ref','var') || isempty(time_ref)
+ % take 1st time point
+ if ids_equilibrium.ids_properties.homogeneous_time == 1
+ time_ref = ids_equilibrium.time(1);
+ else
+ time_ref = ids_equilibrium.time_slice{1}.time;
+ end
+end
+if ~exist('eqdsk_filename','var') || isempty(eqdsk_filename)
+ eqdsk_filename = ['t' num2str(time_ref)];
+end
+if ~exist('write_opt','var') || isempty(write_opt)
+ write_opt = 1;
+end
+
+% find time_slice for equilibrium ids
+if ids_equilibrium.ids_properties.homogeneous_time == 1
+ i_equil = iround_os(ids_equilibrium.time,time_ref);
+else
+ for j=1:numel(ids_equilibrium.time_slice{j})
+ time_array(j) = ids_equilibrium.time_slice{j}.time;
+ end
+ i_equil = iround_os(time_array,time_ref);
+end
+% find time_slice for wall ids
+if ids_wall.ids_properties.homogeneous_time == 1
+ i_wall = iround_os(ids_wall.time,time_ref);
+else
+ for j=1:numel(ids_wall.time_slice{j})
+ time_array(j) = ids_wall.time_slice{j}.time;
+ end
+ i_wall = iround_os(time_array,time_ref);
+end
+
+eqdsk_out.cocos = 11; % assume 11 within IMAS at this stage, should come from ids_equilibrium metadata
+
+% check if needed fields are present and non empty
+if ids_equilibrium.time_slice{i_equil}.global_quantities.ip ~= -9.0000e+40
+ eqdsk_out.ip = ids_equilibrium.time_slice{i_equil}.global_quantities.ip;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.global_quantities.ip not provided']);
+end
+if ids_equilibrium.vacuum_toroidal_field.b0 ~= -9.0000e+40
+ eqdsk_out.b0 = ids_equilibrium.vacuum_toroidal_field.b0;
+else
+ error(['ids_equilibrium.vacuum_toroidal_field.b0 not provided']);
+end
+if ids_equilibrium.vacuum_toroidal_field.r0 ~= -9.0000e+40
+ eqdsk_out.r0 = ids_equilibrium.vacuum_toroidal_field.r0;
+else
+ error(['ids_equilibrium.vacuum_toroidal_field.r0 not provided']);
+end
+if ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.r ~= -9.0000e+40
+ eqdsk_out.raxis = ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.r;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.r not provided']);
+end
+if ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.z ~= -9.0000e+40
+ eqdsk_out.zaxis = ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.z;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.global_quantities.magnetic_axis.z not provided']);
+end
+eqdsk_out.zshift = 0.; % can add an option to shift eqdsk but can use eqdsk_transform
+
+if ids_equilibrium.time_slice{i_equil}.profiles_1d.dpressure_dpsi ~= -9.0000e+40
+ dpressure_dpsi = ids_equilibrium.time_slice{i_equil}.profiles_1d.dpressure_dpsi;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.dpressure_dpsi not provided']);
+end
+if ids_equilibrium.time_slice{i_equil}.profiles_1d.pressure ~= -9.0000e+40
+ pressure = ids_equilibrium.time_slice{i_equil}.profiles_1d.pressure;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.pressure not provided']);
+end
+if ids_equilibrium.time_slice{i_equil}.profiles_1d.f ~= -9.0000e+40
+ f_rbphi = ids_equilibrium.time_slice{i_equil}.profiles_1d.f;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.f not provided']);
+end
+if ids_equilibrium.time_slice{i_equil}.profiles_1d.f_df_dpsi ~= -9.0000e+40
+ f_df_dpsi = ids_equilibrium.time_slice{i_equil}.profiles_1d.f_df_dpsi;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.f_df_dpsi not provided']);
+end
+if ids_equilibrium.time_slice{i_equil}.global_quantities.psi_axis ~= -9.0000e+40
+ psi_axis = ids_equilibrium.time_slice{i_equil}.global_quantities.psi_axis;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.global_quantities.psi_axis not provided']);
+end
+if ids_equilibrium.time_slice{i_equil}.global_quantities.psi_boundary ~= -9.0000e+40
+ psi_boundary = ids_equilibrium.time_slice{i_equil}.global_quantities.psi_boundary;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.global_quantities.psi_boundary not provided']);
+end
+if ids_equilibrium.time_slice{i_equil}.profiles_1d.psi ~= -9.0000e+40
+ psi = ids_equilibrium.time_slice{i_equil}.profiles_1d.psi;
+ rhopol_norm = sqrt((psi-psi(1))./(psi(end)-psi(1)));
+ rhopol_norm(1) = 0.; rhopol_norm(end) = 1.;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.psi not provided, could use rho_tor_norm, not implemented yet']);
+end
+if ids_equilibrium.time_slice{i_equil}.profiles_1d.q ~= -9.0000e+40
+ q_profile = ids_equilibrium.time_slice{i_equil}.profiles_1d.q;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.profiles_1d.q not provided']);
+end
+if numel(ids_equilibrium.time_slice{i_equil}.profiles_2d) >= 1
+ i_2d = 0;
+ % find eqdsk type grid
+ for j=1:numel(ids_equilibrium.time_slice{i_equil}.profiles_2d)
+ if ids_equilibrium.time_slice{i_equil}.profiles_2d{j}.grid_type.index == 1
+ i_2d = j;
+ end
+ end
+ if i_2d > 0
+ % grid ala eqdsk so can use psi, dim1, dim2
+ rmesh = ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.grid.dim1;
+ zmesh = ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.grid.dim2;
+ if ~isempty(ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.psi)
+ psi_2d = ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.psi;
+ else
+ error('ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.psi not provided');
+ end
+ else
+ error([' cannot find ids_equilibrium.time_slice{i_equil}.profiles_2d{i_2d}.grid_type.index == 1']);
+ end
+end
+if ~isempty(ids_equilibrium.time_slice{i_equil}.boundary.outline.r)
+ r_lcfs = ids_equilibrium.time_slice{i_equil}.boundary.outline.r;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.boundary.outline.r not provided']);
+end
+if ~isempty(ids_equilibrium.time_slice{i_equil}.boundary.outline.z)
+ z_lcfs = ids_equilibrium.time_slice{i_equil}.boundary.outline.z;
+else
+ error(['ids_equilibrium.time_slice{i_equil}.boundary.outline.z not provided']);
+end
+
+eqdsk_out.cocos = 11;
+eqdsk_out.nr = numel(rmesh);
+eqdsk_out.nz = numel(zmesh);
+eqdsk_out.rmesh = rmesh;
+eqdsk_out.zmesh = zmesh;
+eqdsk_out.zshift = 0.;
+eqdsk_out.psi = psi_2d;
+eqdsk_out.psirz = psi_2d;
+eqdsk_out.rboxlen = eqdsk_out.rmesh(end) - eqdsk_out.rmesh(1);
+eqdsk_out.rboxleft=eqdsk_out.rmesh(1);
+eqdsk_out.zboxlen = eqdsk_out.zmesh(end) - eqdsk_out.zmesh(1);
+eqdsk_out.zmid = 0.5*(eqdsk_out.zmesh(end) + eqdsk_out.zmesh(1)) ;
+eqdsk_out.psiaxis = psi_axis
+eqdsk_out.psiedge = psi_boundary;
+if ids_equilibrium.time_slice{i_equil}.boundary_separatrix.type == 1
+ % diverted
+ % could check psi value with psi boundary
+end
+
+% need psimesh ot be equidistant and same nb of points as nr
+eqdsk_out.psimesh=linspace(0,1,eqdsk_out.nr)';
+eqdsk_out.rhopsi = sqrt(eqdsk_out.psimesh);
+% psimesh assumed from 0 on axis (so psi_edge is delta_psi) to have correct d/dpsi but to have sign(dpsi) easily
+% (psimesh from 0 to 1) eqdsk_out.psimesh = (eqdsk_out.psiedge-eqdsk_out.psiaxis) .* eqdsk_out.psimesh;
+nrho = size(pressure,1);
+eqdsk_out.p = interpos(rhopol_norm,pressure,eqdsk_out.rhopsi,-0.01,[1 2],[0 pressure(end)]);
+eqdsk_out.pprime = interpos(rhopol_norm,dpressure_dpsi,eqdsk_out.rhopsi,-0.01,[1 2],[0 dpressure_dpsi(end)]);
+eqdsk_out.F = interpos(rhopol_norm,f_rbphi,eqdsk_out.rhopsi,-0.01,[1 2],[0 f_rbphi(end)]);
+eqdsk_out.FFprime = interpos(rhopol_norm,f_df_dpsi,eqdsk_out.rhopsi,-0.01,[1 2],[0 f_df_dpsi(end)]);
+eqdsk_out.q = interpos(rhopol_norm,q_profile,eqdsk_out.rhopsi,-0.01,[1 2],[0 q_profile(end)]);
+
+eqdsk_out.stitle = [eqdsk_filename_suffix ' tref=' num2str(time_ref)];
+
+eqdsk_out.nbbound = numel(ids_equilibrium.time_slice{i_equil}.boundary.outline.r);
+if ~isempty(ids_equilibrium.time_slice{i_equil}.boundary.outline.r) && eqdsk_out.nbbound > 1 && ...
+ eqdsk_out.nbbound == numel(ids_equilibrium.time_slice{i_equil}.boundary.outline.z)
+ eqdsk_out.rplas = ids_equilibrium.time_slice{i_equil}.boundary.outline.r;
+ eqdsk_out.zplas = ids_equilibrium.time_slice{i_equil}.boundary.outline.z;
+else
+ error('Plasma boundary not provided in time_slice{i_equil}.boundary.outline.r');
+end
+
+% extract limiter + double ids_wall if available
+if numel(ids_wall.description_2d)>0 && numel(ids_wall.description_2d{1}.limiter.unit)>0 && ...
+ ~isempty(ids_wall.description_2d{1}.limiter.unit{1}.outline.r) && ...
+ numel(ids_wall.description_2d{1}.vessel.unit)>=2 && ...
+ ~isempty(ids_wall.description_2d{1}.vessel.unit{1}.annular.centreline.r) && ...
+ ~isempty(ids_wall.description_2d{1}.vessel.unit{2}.annular.centreline.r)
+ Rlimiter = ids_wall.description_2d{1}.limiter.unit{1}.outline.r;
+ Zlimiter = ids_wall.description_2d{1}.limiter.unit{1}.outline.z;
+ if numel(ids_wall.description_2d{1}.limiter.unit)>=2 && ~isempty(ids_wall.description_2d{1}.limiter.unit{2}.outline.r)
+ nbpts2 = numel(ids_wall.description_2d{1}.limiter.unit{2}.outline.r);
+ Rlimiter(end+1:end+nbpts2) = ids_wall.description_2d{1}.limiter.unit{2}.outline.r(nbpts2:-1:1);
+ Zlimiter(end+1:end+nbpts2) = ids_wall.description_2d{1}.limiter.unit{2}.outline.z(nbpts2:-1:1);
+ end
+ Rlimiter(end+1) = Rlimiter(1);
+ Zlimiter(end+1) = Zlimiter(1);
+
+ % vessel
+ Rvessel_1 = ids_wall.description_2d{1}.vessel.unit{1}.annular.centreline.r;
+ Zvessel_1 = ids_wall.description_2d{1}.vessel.unit{1}.annular.centreline.z;
+ Rvessel_2 = ids_wall.description_2d{1}.vessel.unit{2}.annular.centreline.r;
+ Zvessel_2 = ids_wall.description_2d{1}.vessel.unit{2}.annular.centreline.z;
+
+ aa = [Rlimiter Zlimiter];
+ nlim = size(aa,1);
+ nvessel_1 = numel(Rvessel_1);
+ nvessel_2 = numel(Rvessel_2);
+
+ d_lim_vessel_1 = sqrt((Rvessel_1-Rlimiter(end)).^2 + (Zvessel_1-Zlimiter(end)).^2);
+ [zzz,idmin] = min(d_lim_vessel_1);
+ aa(nlim+1:nlim+(nvessel_1-idmin+1),1) = Rvessel_1(idmin:end);
+ aa(nlim+(nvessel_1-idmin+1)+1:nlim+nvessel_1+1,1) = Rvessel_1(1:idmin); % end return point
+ aa(nlim+1:nlim+(nvessel_1-idmin+1),2) = Zvessel_1(idmin:end);
+ aa(nlim+(nvessel_1-idmin+1)+1:nlim+nvessel_1+1,2) = Zvessel_1(1:idmin);
+
+ d_lim_vessel_2 = sqrt((Rvessel_2-Rlimiter(end)).^2 + (Zvessel_2-Zlimiter(end)).^2);
+ [zzz,idmin] = min(d_lim_vessel_2);
+ aa(nlim+nvessel_1+1+1:nlim+nvessel_1+1+(nvessel_2-idmin+1),1) = Rvessel_2(idmin:end);
+ aa(nlim+nvessel_1+1+(nvessel_2-idmin+1)+1:nlim+nvessel_1+1+nvessel_2+1,1) = Rvessel_2(1:idmin);
+ aa(nlim+nvessel_1+1+1:nlim+nvessel_1+1+(nvessel_2-idmin+1),2) = Zvessel_2(idmin:end);
+ aa(nlim+nvessel_1+1+(nvessel_2-idmin+1)+1:nlim+nvessel_1+1+nvessel_2+1,2) = Zvessel_2(1:idmin);
+
+else
+ % read from previously defined boundary
+ try
+ aa = load('ITERlimiter.data');
+ catch
+ aa = [];
+ end
+end
+
+if ~isempty(aa)
+ eqdsk_out.nblim = size(aa,1);
+ eqdsk_out.rlim = aa(:,1);
+ eqdsk_out.zlim = aa(:,2);
+else
+ eqdsk_out.nblim = 0;
+ eqdsk_out.rlim = [];
+ eqdsk_out.zlim = [];
+end
+
+if write_opt == 1
+ write_eqdsk(['EQDSK_COCOS' num2str(eqdsk_out.cocos) '_' eqdsk_filename_suffix],eqdsk_out,[eqdsk_out.cocos eqdsk_out.cocos]);
+end
--
GitLab