From aa4e10dad1e30c6e0475eee3b26416bc909c0ff8 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Tue, 26 Mar 2019 11:15:04 +0000
Subject: [PATCH] add aug imas specific files at this stage in AUG_IMAS folder

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11644 d63d8f72-b253-0410-a779-e742ad2e26cf
---
 crpptbx/AUG/gdat_aug.m                     |   2 +-
 crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m | 341 +++++++++++++++++++++
 crpptbx/AUG_IMAS/aug_ids_headpart.m        |  55 ++++
 3 files changed, 397 insertions(+), 1 deletion(-)
 create mode 100644 crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m
 create mode 100644 crpptbx/AUG_IMAS/aug_ids_headpart.m

diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
index 446d8711..823f1f9d 100644
--- a/crpptbx/AUG/gdat_aug.m
+++ b/crpptbx/AUG/gdat_aug.m
@@ -1115,7 +1115,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
     end
     try
       if ~isempty(shot)
-        eval(['[ids_top,ids_top_description]=tcv_get_ids_' ids_top_name '(shot,equil_empty);'])
+        eval(['[ids_top,ids_top_description]=aug_get_ids_' ids_top_name '(shot,equil_empty);'])
         gdat_data.(ids_top_name) = ids_top;
         gdat_data.([ids_top_name '_description']) = ids_top_description;
       else
diff --git a/crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m b/crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m
new file mode 100644
index 00000000..710a5c3d
--- /dev/null
+++ b/crpptbx/AUG_IMAS/aug_get_ids_equilibrium.m
@@ -0,0 +1,341 @@
+function [ids_equilibrium,ids_equilibrium_description,varargout] = aug_get_ids_equilibrium(shot,ids_equil_empty,varargin);
+%
+%  [ids_equilibrium,ids_equilibrium_description,varargout] = aug_get_ids_equilibrium(shot,ids_equil_empty,varargin);
+% 
+%
+
+machine = 'aug';
+
+[ids_equilibrium, params_equilibrium] = aug_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:
+% "global_quantities" with subfields being the relevant data to get and a local structure:
+% "global_quantities_desc" which contains the same subfields themselves containing the gdat string aftre shot used
+%
+% vacuum_toroidal_field and time, using homogeneous and q_rho time base
+%
+% time may be "too long" in shotfile compared to effective equilibrium reconstructed times
+temp.q_rho = gdat(params_equilibrium.shot,'q_rho','machine',machine);
+temp_desc.q_rho = 'q_rho';
+ids_equilibrium.time = temp.q_rho.t;
+ids_equilibrium_description.time = '.t subfield from: [''q_rho'']';
+%
+tens_time = -1;
+
+vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,'b0','machine',machine);
+vacuum_toroidal_field_desc.b0 = ['''b0'' with interpos using tens_time = ' num2str(tens_time)];
+vacuum_toroidal_field_desc.r0 = '.r0 subfield from: [''b0'']';
+ids_equilibrium.vacuum_toroidal_field.r0 = vacuum_toroidal_field.b0.r0;
+ids_equilibrium.vacuum_toroidal_field.b0 = interpos(vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data, ids_equilibrium.time,tens_time);
+ids_equilibrium_description.vacuum_toroidal_field = vacuum_toroidal_field_desc;
+
+ids_equilibrium.time_slice(1:length(ids_equilibrium.time)) = 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:length(ids_equilibrium.time)
+% $$$   ids_equilibrium.time(it)
+% $$$   temp.eqdsks{it}=gdat(params_equilibrium.shot,'eqdsk','time',ids_equilibrium.time(it),'write',0,'machine',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_desc.area = 'area_edge';
+global_quantities.beta_normal = gdat(params_equilibrium.shot,'betan','machine',machine);
+global_quantities_desc.beta_normal = 'betan';
+global_quantities.beta_pol = gdat(params_equilibrium.shot,'betap','machine',machine);
+global_quantities_desc.beta_pol = 'betap';
+global_quantities.beta_tor = gdat(params_equilibrium.shot,'beta','machine',machine);
+global_quantities_desc.beta_tor = 'beta';
+global_quantities.energy_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',machine);
+global_quantities_desc.energy_mhd = 'w_mhd';
+global_quantities.ip = gdat(params_equilibrium.shot,'ip','machine',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);
+global_quantities_desc.li_3 = 'li';
+temp.r_magnetic_axis = gdat(params_equilibrium.shot,'r_axis','machine',machine);
+temp_desc.r_magnetic_axis = 'r_axis';
+temp.z_magnetic_axis = gdat(params_equilibrium.shot,'z_axis','machine',machine);
+temp_desc.z_magnetic_axis = 'z_axis';
+temp.psi_axis = gdat(params_equilibrium.shot,'psi_axis','machine',machine); % needs to add psi_edge since psi_axis ??? assuming 0 dege value
+temp_desc.psi_axis = 'psi_axis';
+global_quantities.psi_boundary = gdat(params_equilibrium.shot,'psi_edge','machine',machine);
+global_quantities_desc.psi_boundary = 'psi_edge';
+global_quantities.q_95 = gdat(params_equilibrium.shot,'q95','machine',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_desc.q_axis = 'q0';
+% surface = gdat(params_equilibrium.shot,'surface','machine',machine); % to be added
+global_quantities.volume = gdat(params_equilibrium.shot,'volume','machine',machine);
+global_quantities_desc.volume = 'volume';
+global_quantities.w_mhd = gdat(params_equilibrium.shot,'w_mhd','machine',machine);
+global_quantities_desc.w_mhd = 'w_mhd';
+
+global_quantities_fieldnames = fieldnames(global_quantities);
+special_fields = {'magnetic_axis', 'psi_axis', 'q_min'}; % fields needing non-automatic treatments
+for i=1:length(global_quantities_fieldnames)
+  try
+    if ~any(strcmp(global_quantities_fieldnames{i},special_fields))
+      if ~isstruct(ids_equilibrium.time_slice{1}.global_quantities.(global_quantities_fieldnames{i}))
+	if length(global_quantities.(global_quantities_fieldnames{i}).data) > length(ids_equilibrium_description.time)
+	  % assume need interpolation with constant extrapolation
+	  aa = interpos(global_quantities.(global_quantities_fieldnames{i}).t, ...
+	  global_quantities.(global_quantities_fieldnames{i}).data,ids_equilibrium.time,tens_time);
+	  for it=1:length(ids_equilibrium.time)
+	    ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = aa(it);
+	  end
+	elseif length(global_quantities.(global_quantities_fieldnames{i}).data) == length(ids_equilibrium_description.time)
+	  % assume same time
+	  if sum(abs(global_quantities.(global_quantities_fieldnames{i}).t-ids_equilibrium_description.time))>0
+	    warning(['not same time for ' global_quantities_fieldnames{i}])
+	  else
+	    for it=1:length(ids_equilibrium.time)
+	      ids_equilibrium.time_slice{it}.global_quantities.(global_quantities_fieldnames{i}) = aa(it);
+	    end
+	  end
+	elseif isempty(global_quantities.(global_quantities_fieldnames{i}).data)
+	  % do not do anything, no data
+	else
+	  warning(['less points than expected, cannot treat for ' global_quantities_fieldnames{i}])
+	end
+      else
+	warning([global_quantities_fieldnames{i} ' is a structure should be in special cases?']);
+      end
+    else
+      special_fields{end+1} = global_quantities_fieldnames{i};
+    end
+  catch ME_global
+    throw(ME_global)
+    keyboard
+  end
+end
+
+% special case
+for it=1:length(ids_equilibrium.time)
+  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) + ...
+      ids_equilibrium.time_slice{it}.global_quantities.psi_boundary;
+  [zz,izz]=min(temp.q_rho.data(:,it));
+  ids_equilibrium.time_slice{it}.global_quantities.q_min.value = zz;
+  ids_equilibrium.time_slice{it}.global_quantities.q_min.rho_tor_norm = temp.q_rho.grids_1d.rhotornorm(izz);
+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);
+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);
+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);
+boundary_desc.triangularity = 'delta';
+boundary.triangularity_lower = gdat(params_equilibrium.shot,'delta_bottom','machine',machine);
+boundary_desc.triangularity_lower = 'delta_bottom';
+boundary.triangularity_upper = gdat(params_equilibrium.shot,'delta_top','machine',machine);
+boundary_desc.triangularity_upper = 'delta_top';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TCV LIUQE specific, to deal later...
+temp.n_x_point = gdat(params_equilibrium.shot,'tcv_eq(''''n_xpts'''',''''liuqe.m'''')','machine',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_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_desc.z_x_point = '''tcv_eq(''''z_xpts'''',''''liuqe.m'''')''';
+temp.rgeom = gdat(params_equilibrium.shot,'rgeom','machine',machine);
+temp_desc.rgeom = 'rgeom';
+temp.zgeom = gdat(params_equilibrium.shot,'zgeom','machine',machine);
+temp_desc.zgeom = 'zgeom';
+temp.r_lcfs = gdat(params_equilibrium.shot,'r_contour_edge','machine',machine);
+temp_desc.r_lcfs = 'r_contour_edge';
+temp.z_lcfs = gdat(params_equilibrium.shot,'z_contour_edge','machine',machine);
+temp_desc.z_lcfs = 'z_contour_edge';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+temp_2d.equil = gdat(params_equilibrium.shot,'equil','machine',machine);
+% at this stage use geteqdsk to get R,Z boundary since already included
+for it=1:length(temp_2d.equil.t)
+  [eqdskAUG{it}, dummy, equil_t_index]=geteqdskAUG(temp_2d.equil,temp_2d.equil.t(it),[],[],[],[],[],'fignb',999);
+end
+
+boundary_fieldnames = fieldnames(boundary);
+special_fields = {'lcfs', 'geometric_axis', 'x_point','outline'}; % fields needing non-automatic treatments
+for it=1:length(ids_equilibrium.time)
+  for i=1:length(boundary_fieldnames)
+    if ~any(strcmp(boundary_fieldnames{i},special_fields))
+      if ~isstruct(ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}))
+        ids_equilibrium.time_slice{it}.boundary.(boundary_fieldnames{i}) = ...
+            boundary.(boundary_fieldnames{i}).data(it);
+      else
+        special_fields{end+1} = boundary_fieldnames{i};
+      end
+    end
+  end
+end
+
+% special cases
+for it=1:length(ids_equilibrium.time)
+  if ~isempty(eqdskAUG{it}.rplas)
+    ids_equilibrium.time_slice{it}.boundary.outline.r = eqdskAUG{it}.rplas;
+    ids_equilibrium.time_slice{it}.boundary.outline.z = eqdskAUG{it}.zplas;
+    ids_equilibrium.time_slice{it}.boundary.lcfs.r = ids_equilibrium.time_slice{it}.boundary.outline.r;
+    ids_equilibrium.time_slice{it}.boundary.lcfs.z = ids_equilibrium.time_slice{it}.boundary.outline.z;
+  end
+  if ~isempty(temp.rgeom.data)
+    ids_equilibrium.time_slice{it}.boundary.geometric_axis.r = temp.rgeom.data(it);
+  end
+  if ~isempty(temp.zgeom.data)
+    ids_equilibrium.time_slice{it}.boundary.geometric_axis.z = temp.zgeom.data(it);
+  end
+  % could use temp_2d.equil.Xpoints but not sure which points to take
+% $$$   if ~isempty(temp.n_x_point.data) && temp.n_x_point.data(it) > 0
+% $$$     ids_equilibrium.time_slice{it}.boundary.x_point(1:temp.n_x_point.data(it)) = ids_equilibrium.time_slice{it}.boundary.x_point(1);
+% $$$     for i=1:length(temp.n_x_point.data(it))
+% $$$       ids_equilibrium.time_slice{it}.boundary.x_point{i}.r = temp.r_x_point.data(i,it);
+% $$$       ids_equilibrium.time_slice{it}.boundary.x_point{i}.z = temp.z_x_point.data(i,it);
+% $$$     end
+% $$$   end
+
+end
+
+%
+%% profiles_1d (cannot use eqdsk since not same radial mesh), but can interpolate and use base from equil
+%
+% area = gdat(params_equilibrium.shot,'area','machine',machine);
+area.data = temp_2d.equil.area;
+% b_average = gdat(params_equilibrium.shot,'b_average','machine',machine);
+% obsolescent: b_average.data = temp_2d.equil.bave;
+% (profile) beta_pol = gdat(params_equilibrium.shot,'','machine',machine);
+b_field_average = temp_2d.equil.bave;
+% 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.data = temp_2d.equil.dpressuredpsi;
+% dpsi_drho_tor = gdat(params_equilibrium.shot,'dpsi_drho_tor','machine',machine);
+dvolume_dpsi = temp_2d.equil.dvoldpsi;
+% 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.data = temp_2d.equil.ffprime;
+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);
+
+profiles_1d_fieldnames = fieldnames(profiles_1d);
+special_fields = {'geometric_axis', 'rho_tor_norm', 'psi'}; % fields needing non-automatic treatments
+for it=1:length(ids_equilibrium.time)
+  for i=1:length(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}))
+        if ~ischar(profiles_1d.(profiles_1d_fieldnames{i}).data) && ~isempty(profiles_1d.(profiles_1d_fieldnames{i}).data) ...
+              && size(profiles_1d.(profiles_1d_fieldnames{i}).data,2)>=it
+          ids_equilibrium.time_slice{it}.profiles_1d.(profiles_1d_fieldnames{i}) = ...
+              profiles_1d.(profiles_1d_fieldnames{i}).data(:,it);
+        end
+      else
+        special_fields{end+1} = profiles_1d_fieldnames{i};
+      end
+    end
+  end
+end
+
+% special cases
+nrho = length(profiles_1d.rho_tor.x);
+ntime = length(temp.psi_axis.data);
+for it=1:length(ids_equilibrium.time)
+  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 = (1-profiles_1d.rho_tor.x.^2).*temp.psi_axis.data(it) + ...
+      global_quantities.psi_boundary.data(it);
+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
+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
+
+profiles_2d_fieldnames = fieldnames(profiles_2d);
+special_fields = {'grid', 'grid_type', 'psi'}; % fields needing non-automatic treatments
+for it=1:length(ids_equilibrium.time)
+  for i=1:length(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}))
+        if ~ischar(profiles_2d.(profiles_2d_fieldnames{i}).data) && ~isempty(profiles_2d.(profiles_2d_fieldnames{i}).data) ...
+              && size(profiles_2d.(profiles_2d_fieldnames{i}).data,3)>=it
+          ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}) = ...
+              profiles_2d.(profiles_2d_fieldnames{i}).data(:,:,it);
+        end
+      else
+        special_fields{end+1} = profiles_2d_fieldnames{i};
+      end
+    end
+  end
+end
+
+% special cases
+for it=1:length(ids_equilibrium.time)
+  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;
+  ids_equilibrium.time_slice{it}.profiles_2d{1}.grid.dim1 = temp_2d.equil.Rmesh(:,it);
+  ids_equilibrium.time_slice{it}.profiles_2d{1}.grid.dim2 = temp_2d.equil.Zmesh(:,it);
+  ids_equilibrium.time_slice{it}.profiles_2d{1}.psi(:,:) =  temp_2d.equil.psi2D(:,:,it);
+end
diff --git a/crpptbx/AUG_IMAS/aug_ids_headpart.m b/crpptbx/AUG_IMAS/aug_ids_headpart.m
new file mode 100644
index 00000000..da7351ef
--- /dev/null
+++ b/crpptbx/AUG_IMAS/aug_ids_headpart.m
@@ -0,0 +1,55 @@
+function [ids_generic, params_ids_generic] = aug_ids_headpart(shot,ids_in,ids_name,varargin);
+%
+% [ids_generic, params_ids_generic] = aug_ids_headpart(shot,ids_in,ids_name,varargin);
+%
+% parses inputs and fill in ids_properties
+%
+%
+% varargin options:
+%
+% 'comment': comment to include in ids_properties
+% 'homogeneous_time': homogeneous_time in ids_properties: 1 (default) if the whole ids has same time, 0 otherwise
+%
+%
+% example:
+%   [ids_equilibrium, params_ids_equilibrium] = aug_ids_headpart(shot,ids_equil_empty,'equilibrium','comment','your comment');
+%
+%
+
+% initialize input parser
+p = inputParser;
+% no required inputs here so one can call with empty args and get defaults parameters
+% effectively required inputs should have defaults as empty
+p.addOptional('shot', [], @(x) (isnumeric(x) && isscalar(x) && (x == round(x)))); % integer
+p.addOptional('ids_in', struct([]), @(x) (isstruct(x)));
+p.addOptional('ids_name', '', @(x) (ischar(x))); % char
+p.addOptional('comment', 'default comment', @(x) isempty(x) || ischar(x)); % char
+p.addOptional('homogeneous_time', 1, @(x) (isnumeric(x) && isscalar(x) && (x == round(x))));
+
+p.parse(shot,ids_in,ids_name);
+defaults_ids_generic = p.Results; % to keep track of defaults
+params = p.Results;
+if nargin >= 4
+  p.parse(shot,'ids_in',ids_in,'ids_name',ids_name,varargin{:});
+  params = p.Results;
+end
+
+% replace empty inputs with defaults
+names = fieldnames(params);
+mask = structfun(@isempty,params);
+if any(mask),
+  params = rmfield(params,unique([names(mask); p.UsingDefaults.']));
+end
+
+params_ids_generic = params;
+
+ids_generic = ids_in;
+
+%
+% ids_properties
+%
+ids_generic.ids_properties.comment = params_ids_generic.comment;
+ids_generic.ids_properties.homogeneous_time = params_ids_generic.homogeneous_time;
+ids_generic.ids_properties.source = ['AUG shot #' num2str(params_ids_generic.shot)];
+ids_generic.ids_properties.provider = ['aug_get_ids_' ids_name];
+ids_generic.ids_properties.creation_date = date;
-- 
GitLab