function [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
%
% 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 = [];
  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);
params_eff.data_request='rhovol';
rhovol = gdat(gdat_data.shot,params_eff);
params_eff.data_request='psi_axis';
psi_axis = gdat(gdat_data.shot,params_eff);
psi0 = interpos(psi_axis.t,psi_axis.data,gdat_data.t,-0.01);

if (nbdim_x == 1)
  gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2*reshape(psi0,1,length(psi0));
elseif (nbdim_x == 2)
  gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2.*repmat(reshape(psi0,1,length(psi0)),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));

it_rt = iround_os(rhotor_norm.t,gdat_data.t);
it_vol = iround_os(rhovol.t,gdat_data.t);
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 (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)==length(gdat_data.grids_1d.rhopolnorm)
      gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm);
      if ~isempty(rhovol.x)
        gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm);
      end
    end
  else
    if length(ii)==size(gdat_data.grids_1d.rhopolnorm,1)
      gdat_data.grids_1d.rhotornorm(:,it)=interpos(rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
      if ~isempty(rhovol.x)
        gdat_data.grids_1d.rhovolnorm(:,it)=interpos(rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it));
      end
    end
  end
end
gdat_data.grids_1d.rhotor_edge=interpos(rhotor_norm.t',rhotor_norm.rhotor_edge,gdat_data.t',-0.01);
gdat_data.grids_1d.volume_edge=interpos(rhovol.t',rhovol.volume_edge,gdat_data.t',-0.01);