fieldstep=1;

maxN=10e15;
Minwell=-1700;
tend=min(8000,length(M.t2d));
tend=length(M.t2d);
plot3d=false;

rAthetpos=25;

thefig=figure('Position',[0 0 1600 900]);
if plot3d
    filename=[M.name,'_well3D.avi'];
else
    filename=[M.name,'_well2D.avi'];
end
videowriterobj=VideoWriter([M.folder,'/',filename]);
videowriterobj.FrameRate=10;
open(videowriterobj);


plotaxes(1)=subplot(1,2,1,'Parent',thefig);
plotaxes(2)=subplot(1,2,2,'Parent',thefig);
if plot3d
    sf=surface(plotaxes(1),M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none');
    xlim(plotaxes(1),[M.zgrid(1) M.zgrid(end)])
    ylim(plotaxes(1),[M.rgrid(1) M.rgrid(sum(M.nnr(1:2))+5)])
    xlabel(plotaxes(1),'z [m]')
    ylabel(plotaxes(1),'r [m]')
    title(plotaxes(1),'Density')
    c = colorbar(plotaxes(1));
    c.Label.String= 'n[m^{-3}]';
    caxis(plotaxes(1),[-Inf maxN]);
    
    well=M.PotentialWell(fieldstep,false);
    sf=surface(plotaxes(2),M.zgrid(2:end),M.rgrid(2:end),well(2:end,2:end),'edgecolor','none');
    xlim(plotaxes(2),[M.zgrid(1) M.zgrid(end)])
    ylim(plotaxes(2),[M.rgrid(1) M.rgrid(sum(M.nnr(1:2))+5)])
    xlabel(plotaxes(2),'z [m]')
    ylabel(plotaxes(2),'r [m]')
    title(plotaxes(2),'Well')
    c = colorbar(plotaxes(2));
    c.Label.String= 'depth [V]';
    view(plotaxes(2),3)
    colormap(plotaxes(2),'jet');
    caxis(plotaxes(2),[Minwell 0]);
    zlim(plotaxes(2),[Minwell 0]);
    
    
    for i=fieldstep:5:tend
        sgtitle(thefig,sprintf('t= %1.2f [ns]',M.t2d(i)*1e9))
        dens=M.N(:,:,i);
        
        well=M.PotentialWell(i,false);
        well=well(2:end,2:end);
        
        plotaxes(1).Children.ZData=dens;
        plotaxes(1).Children.CData=dens;
        
        plotaxes(2).Children.ZData=well;
        plotaxes(2).Children.CData=well;
        writeVideo(videowriterobj,getframe(thefig));
    end
else
    dens=M.N(:,:,fieldstep);
    model=M.potentialwellmodel(fieldstep);
    z=model.z;
    r=model.r;
    pot=model.pot;
    rathet=model.rathet;
    [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,1));
    dens=griddata(Zmesh,M.rAthet,dens,Zmesh,Rmesh);
    plot(plotaxes(1),M.zgrid,dens(rAthetpos,1:end));
    xlim(plotaxes(1),[M.zgrid(1) M.zgrid(end)])
    ylim(plotaxes(1),[0 maxN]);
    xlabel(plotaxes(1),'z [m]')
    ylabel(plotaxes(1),'n [m^{-3}]')
    title(plotaxes(1),'Density')
    grid(plotaxes, 'on')
    
    pot=griddata(z,rathet,pot,Zmesh,Rmesh);
    plot(plotaxes(2),M.zgrid(1:end),pot(rAthetpos,1:end));
    xlim(plotaxes(2),[M.zgrid(1) M.zgrid(end)])
    ylim(plotaxes(2),[Minwell 0]);
    xlabel(plotaxes(2),'z [m]')
    ylabel(plotaxes(2),'depth [eV]')
    title(plotaxes(2),'Well')
    grid on
    
    for i=fieldstep:20:tend
        sgtitle(thefig,sprintf('rA_\\theta=%1.3g [Tm^2] t= %1.2f [ns]',Rmesh(rAthetpos,1),M.t2d(i)*1e9))
        dens=M.N(:,:,i);
        model=M.potentialwellmodel(i);
        z=model.z;
        r=model.r;
        pot=model.pot;
        rathet=model.rathet;
        [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,1));
        dens=griddata(Zmesh,M.rAthet,dens,Zmesh,Rmesh);
        pot=griddata(z,rathet,pot,Zmesh,Rmesh);
        
        plotaxes(1).Children.YData=dens(rAthetpos,1:end);
        plotaxes(2).Children.YData=pot(rAthetpos,1:end);
        writeVideo(videowriterobj,getframe(thefig));
    end
end
close(videowriterobj);