function [nel_out,s_out,nes_out,varargout] = line_average(shot,time,Rdiag,Zdiag,ne_in,rhopol_in,Rmesh,Zmesh,psi_RZmesh,psi_axis,psi_edge);
%
% compute s and ne(s) in s_out and nes_out; where s is distance from first point into the plasma along Rdiag,Zdiag chord
% compute int[ne(s)ds]/int[ds] in nel_out
%
% can get Rmesh, Zmesh and psi_RZmesh from eqdsk at shot,time from gdat if not given, but for many times, should use equil output from gdat
%

nel_out = [];
s_out = [];
nes_out = [];

time_eff = time;

if ~exist('Rmesh') || isempty(Rmesh)
  % do only 1st time
  time_eff = time(1);
  ab=gdat(shot,'eqdsk','time',time_eff);
  eqdsk=ab.eqdsk;
  if isempty(Rdiag)
    % assume H-1 line
    par=sf2par('DCN',shot,'H1LFS_R','DCNgeo');
    Rlfs=par.value/1e3;
    par=sf2par('DCN',shot,'H1LFS_z','DCNgeo');
    Zlfs=par.value/1e3;
    par=sf2par('DCN',shot,'H1HFS_R','DCNgeo');
    Rhfs=par.value/1e3;
    par=sf2par('DCN',shot,'H1HFS_z','DCNgeo');
    Zhfs=par.value/1e3;
    Rdiag=[Rlfs Rhfs];
    Zdiag=[Zlfs Zhfs];
  end
  if length(Rdiag)==2
    nbpoints=50;
    Rout=linspace(Rdiag(1),Rdiag(2),nbpoints);
    Zout=linspace(Zdiag(1),Zdiag(2),nbpoints);
  elseif length(Rdiag)>10
    Rout=Rdiag;
    Zout=Zdiag;
  else
    disp(['unexpected nb of diag points: length(Rdiag) = ',length(Rdiag)]);
    return
  end
  Rmesh = eqdsk.rmesh;
  Zmesh = eqdsk.zmesh;
  psi_RZmesh(:,:,1) = eqdsk.psirz;
  psi_axis(1) = eqdsk.psiaxis;
  psi_edge(1) = eqdsk.psiedge;
end

if prod(size(Rmesh)) == max(size(Rmesh))
  Rmesh_eff = repmat(reshape(Rmesh,length(Rmesh),1),1,length(time_eff));
  Zmesh_eff = repmat(reshape(Zmesh,length(Zmesh),1),1,length(time_eff));
end

for it=1:length(time_eff)
  psi_at_routzout(:,it) = interpos2Dcartesian(Rmesh_eff(:,it),Zmesh(:,it),psi_RZmesh(:,:,it),Rout,Zout);
  ij_in=find((psi_at_routzout(:,it)-psi_axis(it)).*(psi_at_routzout(:,it)-psi_edge(it))<=0);
  if ~isempty(ij_in)
    if ij_in(1)==1
      Rstart=Rout(ij_in(1));
      Zstart=Zout(ij_in(1));
    else
      Rstart = (psi_edge(it)-psi_at_routzout(ij_in(1)-1),it)./diff(psi_at_routzout(ij_in(1)-1:ij_in(1),it)).*diff(Rout(ij_in(1)-1:ij_in(1))) + Rout(ij_in(1)-1);
      Zstart = (psi_edge(it)-psi_at_routzout(ij_in(1)-1),it)./diff(psi_at_routzout(ij_in(1)-1:ij_in(1),it)).*diff(Zout(ij_in(1)-1:ij_in(1))) + Zout(ij_in(1)-1);
    end
    if ij_in(end)==length(Rout)
      Rend=Rout(ij_in(end));
      Zend=Zout(ij_in(end));
    else
      Rend = (psi_edge(it)-psi_at_routzout(ij_in(end),it))./diff(psi_at_routzout(ij_in(end):ij_in(end)+1,it)).*diff(Rout(ij_in(end):ij_in(end)+1)) + Rout(ij_in(end));
      Zend = (psi_edge(it)-psi_at_routzout(ij_in(end),it))./diff(psi_at_routzout(ij_in(end):ij_in(end)+1,it)).*diff(Zout(ij_in(end):ij_in(end)+1)) + Zout(ij_in(end));
    end
    s_eff = sqrt((Rout(ij_in)-Rstart).^2 + (Zout(ij_in)-Zstart).^2);
    rhopol_eff=sqrt((psi_at_routzout(ij_in,it)-psi_axis(it))./(psi_edge(it)-psi_axis(it)));
    if ij_in(1)==1
      rhopol_eff = reshape(rhopol_eff,length(rhopol_eff),1);
      s_eff = reshape(s_eff,length(s_eff),1);
    else
      rhopol_eff=[1 ;reshape(rhopol_eff,length(rhopol_eff),1); 1];
      s_eff=[0 ;reshape(s_eff,length(s_eff),1); sqrt((Rend-Rstart).^2 + (Zend-Zstart).^2)];
    end
    [rho_in_eff,isort]=sort(rhopol_in(:,it));
    if rho_in_eff(1)<1e-3
      % assume 1st point is center
      nes = interpos(rho_in_eff,ne_in(isort,it),rhopol_eff,-0.1,[1 0],[0 0]);
    else
      rho_in_eff = [0 ; rho_in_eff];
      ne_in_eff = [ne_in(isort(1),it); ne_in(isort,it)];
      sigma_in = ones(size(ne_in_eff));
      sigma_in(1)=100;
      nes = interpos(rho_in_eff,ne_in_eff,rhopol_eff,-1,[1 0],[0 0],sigma_in);
    end
    nes_out(:,it) = nes;
    s_out(:,it) = s_eff;
    [dum1,~,~,nel_out_prof] = interpos(s_eff,nes,-0.1);
    nel_out(it) = nel_out_prof(end)./(s_eff(end)-s_eff(1));
  end
end