function [signal_with_1Dgrids,varargout] = addpsi_fromRZ(main_signal,equil,varargin)
%
% [signal_with_1Dgrids,varargout] = addpsi_fromRZ(main_signal,equil,varargin);
%
% At this stage tries with cxrs which has cxrs.r(:,time) and cxrs.z(:,time) for the R,Z points
%

% $$$ shot=30382;
% $$$ time=3.;
% $$$ 
% $$$ main_signal=gdat(shot,'main_signal',0);
% $$$ equil=gdat(shot,'equil',0);

if ~exist('equil')
  disp('no equil')
  return
end
if ~exist('main_signal')
  disp('no main_signal')
  return
end

inb_chord_main_signal=size(main_signal.data,1);
inb_time_main_signal=size(main_signal.data,2);
psi_out = NaN*ones(inb_chord_main_signal,inb_time_main_signal);
rhopsinorm_out = NaN*ones(inb_chord_main_signal,inb_time_main_signal);
rhotornorm_out = NaN*ones(inb_chord_main_signal,inb_time_main_signal);
rhovolnorm_out = NaN*ones(inb_chord_main_signal,inb_time_main_signal);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
iok=find(main_signal.r(:,1)>0);
tic
for itequil=1:length(time_equil)-1
  rr=equil.Rmesh(:,itequil);
  zz=equil.Zmesh(:,itequil);
  psirz_in = equil.psi2D(:,:,itequil);
  it_main_signal_inequil = find(main_signal.t>=time_equil(itequil) & main_signal.t<=time_equil(itequil+1));
  if ~isempty(it_main_signal_inequil)
    rout=main_signal.r(iok,it_main_signal_inequil);
    zout=main_signal.z(iok,it_main_signal_inequil);
    psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
    psi_out(iok,it_main_signal_inequil) = reshape(psi_at_routzout,length(iok),length(it_main_signal_inequil));
    rhopsinorm_out(iok,it_main_signal_inequil) = sqrt((psi_out(iok,it_main_signal_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
    for it_cx=1:length(it_main_signal_inequil)
      rhotornorm_out(iok,it_main_signal_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(iok,it_main_signal_inequil(it_cx)),-3,[2 2],[0 1]);
      rhovolnorm_out(iok,it_main_signal_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(iok,it_main_signal_inequil(it_cx)),-3,[2 2],[0 1]);
    end
  end
end
toc

signal_with_1Dgrids = main_signal;
signal_with_1Dgrids.psi_on_rztime = psi_out;
signal_with_1Dgrids.rhopsinorm_on_rztime = rhopsinorm_out;
signal_with_1Dgrids.rhotornorm_on_rztime = rhotornorm_out;
signal_with_1Dgrids.rhovolnorm_on_rztime = rhovolnorm_out;

% $$$ iok=find(main_signal.r(:,1)>0);
% $$$ rout=main_signal.r(iok,1);
% $$$ zout=main_signal.z(iok,1);
% $$$ 
% $$$ % [farray_out,dfarrdx_out,dfarrdy_out,d2farrdx2_out,d2farrdy2_out,d2farrdxdy_out]= ...
% $$$ [farray_out]= interpos2Dcartesian(rr',zz',farray_in,rout,zout);
% $$$ 
% $$$ figure;
% $$$ contour(rr,zz,farray_in(:,:,1)',100);
% $$$ hold on
% $$$ plot(rout,zout,'*')
% $$$ axis equal
% $$$ colorbar
% $$$ 
% $$$ figure;
% $$$ contour(rr,zz,farray_in(:,:,30)',linspace(equil.psi_axis(30),0..*equil.psi_axis(30)+1.0.*equil.psi_lcfs(30),100));
% $$$ hold on
% $$$ plot(rout,zout,'*')
% $$$ axis equal
% $$$ colorbar
% $$$ 
% $$$ [zzz itime]=min(abs(equil.t-time-1.5));
% $$$ [zzz iz]=min(abs(zz-0.05));
% $$$ figure
% $$$ plot(rr,farray_in(:,iz,itime),'*')
% $$$ hold on
% $$$ plot(rout(iok),farray_out(iok,1,itime),'ro')
% $$$ 
% $$$ figure;
% $$$ contour(rr,zz,farray_in(:,:,itime)',farray_out(iok,1,itime));
% $$$ hold on
% $$$ plot(rout,zout,'*')
% $$$ axis equal
% $$$ colorbar