nc=200
np=1000

norm=nc/(np-3)

vamp=0.05
voff=0.
vsig=1.
vshift=0;

dt=0.05

//  ============================================
//  initialize particle positions and velocities
//  ============================================

for i=1:np-3
    voff=voff*vsig;
    x(i)=(i-0.5)*norm;
    v(i)=vamp*(rand()+rand()+rand()+rand()...
              +rand()+rand()+rand()+rand()...
              +rand()+rand()+rand()+rand()-6.)+voff;
end

/// =============
/// electron beam
/// =============

bunchpos=100;
for i=np-2:np
    x(i)=bunchpos;
    v(i)=5.;
end

plot(x,v,,'linest','none','marker','.','Markers',5)

for k=1:nc
    cell(k)=k;
end

loop=0
while loop<1000
    loop=loop+1;
    printf('%d\n',loop)


//  ======================
//  compute charge density
//  ======================

    for k=1:nc
        rho(k)=0.0;
    end

    for i=1:np
        k=int(x(i));
        k1=k+1;
        k2=k+2;
        if k2>nc then
           k2=1;
        end

        dx=x(i)-k;
        rho(k1)=rho(k1)-(1.-dx);
        rho(k2)=rho(k2)-dx;
    end   
 
    for k=1:nc
        rho(k)=norm*rho(k)+1.;
    end

//  ======================
//  compute electric field
//  ======================

    esum=0.;
    for k=1:nc
        for k1=1:k
            k2=k1+1;
            if k2>nc then
               k2=1;
            end

            esum=esum+(rho(k1)+rho(k2))/2.;
        end
    end

    e(1)=-esum/nc;

    for k=2:nc
        e(k)=e(k-1)+(rho(k)+rho(k-1))/2.;
    end

//  ==================
//  move the particles
//  ==================

    for i=1:np
        k=int(x(i));
        k1=k+1;
        k2=k+2;
        if k2>nc then
           k2=1;
        end

        dx=x(i)-k;
        force=-e(k1)*(1.-dx)-e(k2)*dx;
       
        v(i)=v(i)+force*dt;
        x(i)=x(i)+v(i)*dt;

        while x(i)<0., x(i)=x(i)+nc; end
        while x(i)>nc, x(i)=x(i)-nc; end      
    end
    
    if modulo(loop-1,20) == 0 then
       plot(x,v+vshift,,'linest','none','marker','.','Markers',5)
       vshift=vshift+10.;
    end
end



