Skip to content
Snippets Groups Projects
Commit 0c4c8438 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 9569692a
No related branches found
No related tags found
1 merge request!117add IDS/IMAS to eqdsk directly
Pipeline #94749 passed
Stage: test
Stage: post-test
Hide whitespace changes
Inline Side-by-side
matlab/IMAS/ITERlimiter.data 0 → 100644
+ 346
− 0
View file @ 0c4c8438
% 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
matlab/IMAS/gdat_imas.m
+ 27
− 0
View file @ 0c4c8438
...@@ -518,6 +518,33 @@ elseif strcmp(mapping_for_imas.method,'switchcase') ...@@ -518,6 +518,33 @@ elseif strcmp(mapping_for_imas.method,'switchcase')
end end
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 otherwise
if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_imas']); end if (gdat_params.nverbose>=1); warning(['switchcase= ' data_request_eff ' not known in gdat_imas']); end
......
matlab/IMAS/ids2eqdsk.m 0 → 100644
+ 256
− 0
View file @ 0c4c8438
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment