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

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

inb_chord_cxrs=size(cxrs.r_time,1);
inb_time_cxrs=size(cxrs.r_time,2);
psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhopsinorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
% 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(cxrs.r_time(:,1)>0);
tic
for itequil=1:length(time_equil)-1
  rr=equil.Rmesh(:,itequil);
  zz=equil.Zmesh(:,itequil);
  psirz_in = equil.psi2D(:,:,itequil);
  it_cxrs_inequil = find(cxrs.t>=time_equil(itequil) & cxrs.t<=time_equil(itequil+1));
  if ~isempty(it_cxrs_inequil)
    rout=cxrs.r_time(iok,it_cxrs_inequil);
    zout=cxrs.z_time(iok,it_cxrs_inequil);
    psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout);
    psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil));
    rhopsinorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)));
    for it_cx=1:length(it_cxrs_inequil)
      rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
      rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopsinorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]);
    end
  end
end
toc

cxrs_all = cxrs;
cxrs_all.psi_on_rztime = psi_out;
cxrs_all.rhopsinorm_on_rztime = rhopsinorm_out;
cxrs_all.rhotornorm_on_rztime = rhotornorm_out;
cxrs_all.rhovolnorm_on_rztime = rhovolnorm_out;

% $$$ iok=find(cxrs.r_time(:,1)>0);
% $$$ rout=cxrs.r_time(iok,1);
% $$$ zout=cxrs.z_time(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