get_grids_1d.m

function [gdat_data] = get_grids_1d(gdat_data,nbdim_x,nopt,nverbose);
%
% [gdat_data] = get_grids_1d(gdat_data,nbdim_x,nopt,nverbose);
%
% add various rhos in grids_1d
%
% assume gdat_data.x=rhopol, gdat_data.t, and gdat_data(x,t)
%
% nbdim_x = 1: same x(1:length(x)) (rhopol) for all times
%         = 2: x(:,t), rhopol depends on time like Thomson projection
%
% nopt = 0: do not fill in, just make the empty structure
%      = 1: do compute various fields
%
% compute psi, rhotor_edge, rhotornorm, volume_edge and rhovol and b0 used for rhotor
%
% use gdat calls to psi_axis, psi_edge, rhotor, etc with same basic parameters as data.gdat_params
%

gdat_data.grids_1d.rhopolnorm = gdat_data.x;
if (nopt == 0) || isempty(gdat_data.x) || isempty(gdat_data.t) || isempty(gdat_data.data) || ischar(gdat_data.data)
  gdat_data.grids_1d.rhotornorm = [];
  gdat_data.grids_1d.rhovolnorm = [];
  gdat_data.grids_1d.psi = [];
  gdat_data.grids_1d.rhotor_edge = [];
  gdat_data.grids_1d.volume_edge = [];
  gdat_data.grids_1d.b0 = [];
  return
end
params_eff = gdat_data.gdat_params;params_eff.doplot=0;
params_eff.data_request='rhotor_norm';
rhotor_norm = gdat(gdat_data.shot,params_eff);
ndim_x_rhotor = length(find(size(rhotor_norm.x)>1));
params_eff.data_request='rhovol';
rhovol = gdat(gdat_data.shot,params_eff);
ndim_x_rhovol = length(find(size(rhotor_norm.x)>1));
% check that rhotor and rhovol come from same equil psi
if sum(abs(size(rhotor_norm.x)-size(rhovol.x)))~=0 || sum(sum(abs(rhotor_norm.x-rhovol.x)))>1e-10 ...
      || sum(abs(rhotor_norm.t-rhovol.t))>1e-10;
  error(['get_grids_1d: x and t arrays for rhotor_norm and rhovol are different, although should refer to same equilibrium' char(10) ...
         'size(rhotor_norm.x) = ' num2str(size(rhotor_norm.x)) '; size(rhovol.x) = ' num2str(size(rhovol.x))]);
end
params_eff.data_request='psi_axis';
psi_axis = gdat(gdat_data.shot,params_eff);
params_eff.data_request='psi_edge';
psi_edge = gdat(gdat_data.shot,params_eff);
ij=~isnan(psi_axis.t);
psi0_edge = interpos(63,psi_axis.t(ij),psi_axis.data(ij)-psi_edge.data(ij),gdat_data.t,-0.01);
if (nbdim_x == 1)
  gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2*reshape(psi0_edge,1,length(psi0_edge));
elseif (nbdim_x == 2)
  gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2.*repmat(reshape(psi0_edge,1,length(psi0_edge)),size(gdat_data.grids_1d.rhopolnorm,1),1);
else
  if nverbose>=0; disp(['option: nbdim_x = ' numstr(nbdim_x) ' not implemented, check with O. Sauter']); end
  return
end
gdat_data.grids_1d.rhotornorm = NaN*ones(size(gdat_data.data));
gdat_data.grids_1d.rhovolnorm = NaN*ones(size(gdat_data.data));

if (isempty(rhotor_norm.x) ||isempty(rhotor_norm.t) || isempty(rhotor_norm.data)) ...
      || (isempty(rhovol.x) ||isempty(rhovol.t) || isempty(rhovol.data))
  return
end
it_rt = iround_os(rhotor_norm.t,gdat_data.t);
it_vol = iround_os(rhovol.t,gdat_data.t);
tens0=0.;
for it=1:length(gdat_data.t)
  % do an interpolation on closest point to avoid 2D interp
  it_rt_eff = it_rt(it);
  it_vol_eff = it_vol(it);
  if ndim_x_rhotor==1
    gdat_data.grids_1d.rhopolnorm_equil(:,it) = rhotor_norm.x; % to have same sizes for all 3 reference rhos from equil
  else
    gdat_data.grids_1d.rhopolnorm_equil(:,it) = rhotor_norm.x(:,it_rt_eff);
  end
  gdat_data.grids_1d.rhotornorm_equil(:,it) = rhotor_norm.data(:,it_rt_eff);
  % can assume same .x for rhotor and rhovol since called with same params thus equil (checked with sum before)
  gdat_data.grids_1d.rhovolnorm_equil(:,it) = rhovol.data(:,it_vol_eff);
  %
  if (nbdim_x == 1)
    ii=find(isfinite(gdat_data.grids_1d.rhopolnorm));
  else
    ii=find(isfinite(gdat_data.grids_1d.rhopolnorm(:,it)));
  end
  if (nbdim_x == 1)
    if length(ii) >0
      nb_ii = length(ii);
      if ~isempty(rhotor_norm.x) && ~isempty(rhotor_norm.data)
        if ndim_x_rhotor==1
          gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii),tens0);
        else
          gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii),tens0);
        end
      end
      if ~isempty(rhovol.x) && ~isempty(rhovol.data)
	if ndim_x_rhovol==1
	  gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii),tens0);
	else
	  gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii),tens0);
	end
      end
    end
  else
% $$$     if length(ii)==size(gdat_data.grids_1d.rhopolnorm,1)
    if length(ii) >0
      if ~isempty(rhotor_norm.x) && ~isempty(rhotor_norm.data)
        nb_ii = length(ii);
        if ndim_x_rhotor==1
          gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it),tens0);
        else
          gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it),tens0);
        end
      end
      if ~isempty(rhovol.x) && ~isempty(rhovol.data)
	if ndim_x_rhovol==1
	  gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it),tens0);
	else
	  gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it),tens0);
	end
      end
    end
  end
end
gdat_data.grids_1d.rhotor_edge=interpos(-63,rhotor_norm.t',rhotor_norm.rhotor_edge,gdat_data.t',-0.01);
gdat_data.grids_1d.volume_edge=interpos(-63,rhovol.t',rhovol.volume_edge,gdat_data.t',-0.01);
gdat_data.grids_1d.b0=interpos(-63,rhotor_norm.t',rhotor_norm.b0,gdat_data.t',-0.01);