t2dlength=size(M.t2d,1);
tstep=t2dlength;
rpos=96;
zpos=320;
rAthet_pos=M.rAthet(rpos,zpos);
vpar0=sqrt(M.kb*22000/M.msim*M.weight);


model=M.potentialwellmodel(tstep);
N=M.N(:,:,tstep);
Psieval=rAthet_pos*ones(M.nz+1,1);
Zeval=M.zgrid;

%rdisp=sort(unique([M.rAthet(:,end);M.rAthet(end,:)]));
[Zmesh,Rmesh]=meshgrid(M.zgrid,M.rgrid);
pot=griddata(model.z,model.rathet,model.pot,Zeval,Psieval);
[Zinit,~]=meshgrid(M.zgrid,M.rAthet(:,1));
N0=griddata(Zinit,M.rAthet,N,Zeval,Psieval);
pot=pot+min(pot);

% Magnetic field mirror ratio at each grid position
% compared to local position
R=griddata(Zmesh,M.rAthet,M.B',Zeval,Psieval,'natural')/M.B(zpos,rpos);

% Electrostatic potential on magnetic field line
% coordinates
phis=griddata(Zmesh,M.rAthet,M.pot(:,:,tstep),Zeval,Psieval,'natural')-M.pot(rpos,zpos,tstep);

rposline=griddata(Zmesh,M.rAthet,Rmesh,Zeval,Psieval,'natural');

%%
figure
subplot(3,1,1)
contourf(M.zgrid,M.rgrid,N,'edgecolor','none')
hold on
plot(M.zgrid(zpos),M.rgrid(rpos),'rx')
plot(Zeval,rposline,'r--')
contour(M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'w-.','linewidth',1.5,'Displayname','Vessel Boundaries');
ylim([1e-3 12e-3])

subplot(3,1,2)
plot(M.zgrid,pot)
ylabel('-e\Delta\phi [eV]')
yyaxis right
plot(M.zgrid,N0)
ylabel('n [m^{-3}]')
xlabel('z [m]')

subplot(3,1,3)
EparB=0.5*M.msim/M.weight*((M.Ez(rpos,zpos,tstep)*M.Br(zpos,rpos)-M.Er(rpos,zpos,tstep)*M.Bz(zpos,rpos))/M.B(zpos,rpos)^2)^2*(1-R)/M.qe;
Eparphi=-M.qsim/M.weight*phis/M.qe;
Epar=EparB+Eparphi+0.5*M.msim/M.weight/M.qe*vpar0^2;
plot(M.zgrid,Epar,'displayname','E_{par,tot}')
hold on
plot(M.zgrid,EparB,'displayname','E_{par,mag}')
plot(M.zgrid,Eparphi,'displayname','E_{par,elec}')
plot(M.zgrid,0.5*M.msim/M.weight*(vpar0^2)/M.qe*ones(size(M.zgrid)),'displayname','E_{par,0}')
grid on
xlabel('z [m]')
ylabel('E_{par} [eV]')
legend('location','south')
%ylim(1000*[-1 1])