Commit 35439908 authored by Olivier Sauter's avatar Olivier Sauter
Browse files

add IDS/IMAS to eqdsk directly

tested: gdat(131025,'eqdsk','run',43,'database_user','public','tokamak_name','ITER','machine','IMAS','time',208);
parent f5bc4eb7
Pipeline #94752 passed with stages
in 3 minutes and 2 seconds
% 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
......@@ -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
......
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