fieldstep=1;

maxN=3e17;
tend=min(8000,length(M.t2d));
tend=length(M.t2d);

filename=[M.name,'_dens.avi'];
videowriterobj=VideoWriter([M.folder,'/',filename]);
videowriterobj.FrameRate=10;
videowriterobj.Quality=90;
open(videowriterobj);

step=60;


%% Fields
f=figure('Name', sprintf('%s fields',M.name),'Position',[0 0 1600 900]);
ax1=gca;

dens=M.N(:,:,fieldstep);
dispdens=dens;
geomw=M.geomweight(:,:,1);
geomw(geomw>0)=1;
geomw(geomw<=0)=0;
dispdens(geomw<=0)=0;
[h,curve]=contourf(ax1,M.zgrid*1000,M.rgrid*1000,dispdens,'Displayname','n_e [m^{-3}]');
hold on;
pot=M.pot(:,:,fieldstep);
pot(geomw<0)=0;
Blines=zeros(size(M.Athet));
for i=1:length(M.zgrid)
    Blines(i,:)=M.Athet(i,:).*M.rgrid';
end
zindex=floor(length(M.zgrid/2));
levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),10);
[~,h1]=contour(ax1,M.zgrid*1000,M.rgrid*1000,Blines',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines');

[c1,h2]=contour(ax1,M.zgrid*1000,M.rgrid*1000,pot./1e3,'m--','ShowText','on','linewidth',1.2,'Displayname','Equipotentials [kV]');
clabel(c1,h2,'Color','white')

%contour(ax1,M.zgrid*1000,M.rgrid*1000,M.geomweight(:,:,1),[0 0],'w-','linewidth',1.5);

[c1,hContour]=contourf(ax1,M.zgrid*1000,M.rgrid*1000,geomw);
drawnow;
xlim(ax1,[M.zgrid(1)*1000 M.zgrid(end)*1000])
if(M.conformgeom)
    ylim(ax1,[M.rgrid(1)*1000 M.rgrid(rgridend)*1000])
else
    ylim(ax1,[M.rgrid(1)*1000 M.rgrid(end)*1000])
end
legend([h1,h2],{'Magnetic field lines','Equipotentials [kV]'},'location','southwest')
xlabel(ax1,'z [mm]')
ylabel(ax1,'r [mm]')
%title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens)))
c = colorbar(ax1);
c.Label.String= 'n[m^{-3}]';
view(ax1,2)
%set(h,'edgecolor','none');
grid on;

hFills=hContour.FacePrims;
[hFills.ColorType] = deal('truecoloralpha');  % default = 'truecolor'
try
    hFills(1).ColorData = uint8([150;150;150;255]);
    for idx = 2 : numel(hFills)
        hFills(idx).ColorData(4) = 0;   % default=255
    end
catch
end
rectangle('Position',[M.zgrid(1) M.rgrid(1)-2e-4 M.zgrid(end)-M.zgrid(1) 2e-4]*1000,'FaceColor',[150 150 150]/255)
rectangle('Position',[M.zgrid(1) M.rgrid(end) M.zgrid(end)-M.zgrid(1) 2e-4]*1000,'FaceColor',[150 150 150]/255)
sgtitle(f,sprintf('t= %1.2f [ns]',M.t2d(i)*1e9))
if isfield(M.maxwellsrce,'time_start') && M.maxwellsrce.time_start<=M.t2d(1)
  rectangle('Position',[M.maxwellsrce.zlim(1) M.maxwellsrce.rlim(1) diff(M.maxwellsrce.zlim) diff(M.maxwellsrce.rlim)]*1000,'linestyle','--')
end  

%set(ax1,'ColorScale','log')
ylim([M.rgrid(1)-2e-4 M.rgrid(end)+2e-4]*1000)

caxis([-inf maxN])
for i=fieldstep:step:tend
    sgtitle(f,sprintf('t= %1.2f [ns]',M.t2d(i)*1e9))
    dens=M.N(:,:,i);
    dispdens=dens;
    dispdens(geomw<=0)=0;
    pot=M.pot(:,:,i)./1e3;
    pot(geomw<0)=0;
    
    curve.ZData=dispdens;
    h2.ZData=pot;
    %Curve2.YData=pot(rAthetpos,1:end);
    drawnow;
    hFills=hContour.FacePrims;
[hFills.ColorType] = deal('truecoloralpha');  % default = 'truecolor'
try
    hFills(1).ColorData = uint8([150;150;150;255]);
    for idx = 2 : numel(hFills)
        hFills(idx).ColorData(4) = 0;   % default=255
    end
catch
end
    writeVideo(videowriterobj,getframe(f));
end
close(videowriterobj);