clear all
global N Rtop Rbottom Rleft Rright Vtb Vlr q
m1=120e-26
N=1000;
M=1000;
epsilon=8.85e-12;
q=1.602e-19;
Vtbold=400;
Vlrold=-Vtbold;
r=10;
tmax=0.0001;
dt=0.000000001
nstep=55;
Q=0.4;
dd=2*r+1;
omega=sqrt(4*1.602e-19*Vtbold/m1/2e-6/Q);
U=0.04*m1*omega^2*2e-6/8/1.602e-19;
m2=8*U*1.602e-19/omega^2/2e-6/0.5;
Rtop=[N/2+r+1;8*N/10+1+r];
Rbottom=[N/2+1+r;N/10+1+r];
Rleft=[8*N/10+1+r;N/2+1+r];
Rright=[N/10+1+r;N/2+1+r];
Y1=zeros(6,1);
Y1(1)=-500;
Y1(2)=Rright(1)+(Rleft(1)-Rright(1))/2;
Y1(3)=(Rtop(2)-Rbottom(2))/2+Rbottom(2)+100;
Y1(4)=100000000;
V=m1*(Y1(4)/10^5)^2/2/1.602e-19;
Y2=Y1;
Y2(2)=N/2-5;
xpos1(1)=Y1(1);
ypos1(1)=Y1(2);
zpos1(1)=Y1(3);
xpos2(1)=Y2(1);
ypos2(1)=Y2(2);
zpos2(1)=Y2(3);
f=2;
g=2;
p=1;
actop=zeros(N+1,M+1);
V1=zeros(N+1,M+1);
V2=V1;
btop=actop;
dtop=btop;
bbottom=btop;
dbottom=dtop;
acbottom=actop;
for j=1:dd
    add=r/2*sin((acos(abs(r+1-j)/r)));
    actop(8*N/10+j,:)=N/2+1+add;
    actop(N/10+j,:)=N/2+1+add;
    btop(N/2-r/2+j,:)=8*N/10+1+add;
    acbottom(8*N/10+j,:)=N/2+1-add;
    acbottom(N/10+j,:)=N/2+1-add;
    bbottom(N/2-r/2+j,:)=8*N/10+1-add;
    dtop(N/2-r/2+j,:)=N/10+1+add;
    dbottom(N/2-r/2+j,:)=N/10+1-add;
end
V1old=actop;
V2old=btop+dtop;
x=1:M+1;
for i=1:10
    Vtb=U+sin(50*i*dt*0.5*omega-500*dt*0.5*omega)*Vtbold;
    Vlr=-U+sin(50*i*dt*0.5*omega+pi-500*dt*0.5*omega)*Vtbold;
    V1=Vtb*ones(N,M);
    V2=Vlr*ones(N,M);
    Y1=Y1+[50;0;0;0;0;0];
    Y2=Y2+[50;0;0;0;0;0];
     h0=mesh(x,8*N/10+1:8*N/10+dd,actop(8*N/10+1:8*N/10+dd,:),V1(8*N/10+1:8*N/10+dd,:))
view(114,8);
caxis manual
caxis([-Vtbold Vtbold]);
colormap(winter)
colorbar
axis manual
axis([-500 M+500 0 N 0 N -Vtbold Vtbold]);
hold on
axis off
        h1=mesh(x,8*N/10+1:8*N/10+dd,acbottom(8*N/10+1:8*N/10+dd,:),V1(8*N/10+1:8*N/10+dd,:))
    h2=mesh(x,N/2-r/2+1:N/2-r/2+dd,btop(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h3=mesh(x,N/10+1:N/10+dd,acbottom(N/10+1:N/10+dd,:),V1(N/10+1:N/10+dd,:))
    h4=mesh(x,N/10+1:N/10+dd,actop(N/10+1:N/10+dd,:),V1(N/10+1:N/10+dd,:))
    h5=mesh(x,N/2-r/2+1:N/2-r/2+dd,bbottom(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h6=mesh(x,N/2-r/2+1:N/2-r/2+dd,dtop(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h7=mesh(x,N/2-r/2+1:N/2-r/2+dd,dbottom(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    refreshdata
        pause(0.1)
        hold on
 
            xpos1(f)=Y1(1);
            ypos1(f)=Y1(2);
            zpos1(f)=Y1(3);
            f=f+1;
            Plot1=plot3(xpos1,ypos1,zpos1);
            hold on
            if f>6
                xpos11=xpos1;
                ypos11=ypos1;
                zpos11=zpos1;
                for j=1:5
                    xpos1(j)=xpos11(j+1);
                    ypos1(j)=ypos11(j+1);
                    zpos1(j)=zpos11(j+1);                    
                end
                f=f-1;
        end
            xpos2(g)=Y2(1);
            ypos2(g)=Y2(2);
            zpos2(g)=Y2(3);
            g=g+1;
            plot3(xpos2,ypos2,zpos2,'-r')
            if g>6
                xpos21=xpos2;
                ypos21=ypos2;
                zpos21=zpos2;
                for j=1:5
                    xpos2(j)=xpos21(j+1);
                    ypos2(j)=ypos21(j+1);
                    zpos2(j)=zpos21(j+1);                    
                end
                g=g-1;
            end
                   saveas(gcf, ['1rfma',sprintf('%04d',p),'.jpg'])
        p=p+1;
        close all
end
    

for i=1:tmax/dt
    i
    Vtb=U+sin(i*dt*0.5*omega)*Vtbold;
    Vlr=-U+sin(i*dt*0.5*omega+pi)*Vtbold;
    V1=Vtb*ones(N,M);
    V2=Vlr*ones(N,M);
    
    f11 = quad(i*dt,Y1,m1);
    f21 = quad(i*dt+dt/2,Y1+f11.*dt/2,m1);
    f31 = quad(i*dt+dt/2,Y1+f21 .*  (dt/2),m1  );
    f41 = quad(i*dt+dt ,Y1+f31 .*   dt,m1     );
    Y1 = Y1 + (f11 + f21.*2 + f31.*2 + f41) .* (dt/6);
    if mod(i,nstep)==0
        if ypos1(f-1)<1.5*Rleft(1) & ypos1(f-1)>-Rright(1) & zpos1(f-1)>-Rbottom(2) &zpos1(f-1)<1.5*Rtop(2)
            xpos1(f)=Y1(1);
            ypos1(f)=Y1(2);
            zpos1(f)=Y1(3);
            f=f+1;
            Plot1=plot3(xpos1,ypos1,zpos1);
            hold on
            if f>6
                xpos11=xpos1;
                ypos11=ypos1;
                zpos11=zpos1;
                for j=1:5
                    xpos1(j)=xpos11(j+1);
                    ypos1(j)=ypos11(j+1);
                    zpos1(j)=zpos11(j+1);                    
                end
                f=f-1;
            end
        end
    end
    f12 = quad(i*dt,Y2,m1);
    f22 = quad(i*dt+dt/2,Y2+f11.*dt/2,m1);
    f32 = quad(i*dt+dt/2,Y2+f21 .*  (dt/2),m1  );
    f42 = quad(i*dt+dt ,Y2+f31 .*   dt,m1     );
    Y2 = Y2 + (f12 + f22.*2 + f32.*2 + f42) .*(dt/6);
    if mod(i,nstep)==0
        if ypos2(g-1)<1.5*Rleft(1) & ypos2(g-1)>-Rright(1) & zpos2(g-1)<1.5*Rtop(2) & zpos2(g-1)>-Rbottom(2)
            xpos2(g)=Y2(1);
            ypos2(g)=Y2(2);
            zpos2(g)=Y2(3);
            g=g+1;
            plot3(xpos2,ypos2,zpos2,'-r')
            if g>6
                xpos21=xpos2;
                ypos21=ypos2;
                zpos21=zpos2;
                for j=1:5
                    xpos2(j)=xpos21(j+1);
                    ypos2(j)=ypos21(j+1);
                    zpos2(j)=zpos21(j+1);                    
                end
                g=g-1;
            end
        end
    end
    if mod(i,nstep)==0
        h0=mesh(x,8*N/10+1:8*N/10+dd,actop(8*N/10+1:8*N/10+dd,:),V1(8*N/10+1:8*N/10+dd,:))
view(114,8);
caxis manual
caxis([-Vtbold Vtbold]);
colormap(winter)
colorbar
axis manual
axis([-500 M+500 0 N 0 N -Vtbold Vtbold]);
hold on
axis off
        h1=mesh(x,8*N/10+1:8*N/10+dd,acbottom(8*N/10+1:8*N/10+dd,:),V1(8*N/10+1:8*N/10+dd,:))
    h2=mesh(x,N/2-r/2+1:N/2-r/2+dd,btop(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h3=mesh(x,N/10+1:N/10+dd,acbottom(N/10+1:N/10+dd,:),V1(N/10+1:N/10+dd,:))
    h4=mesh(x,N/10+1:N/10+dd,actop(N/10+1:N/10+dd,:),V1(N/10+1:N/10+dd,:))
    h5=mesh(x,N/2-r/2+1:N/2-r/2+dd,bbottom(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h6=mesh(x,N/2-r/2+1:N/2-r/2+dd,dtop(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h7=mesh(x,N/2-r/2+1:N/2-r/2+dd,dbottom(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    %set(h0,'Cdata',V1(8*N/10+1:8*N/10+dd,:))
        %set(h1,'Cdata',V1(8*N/10+1:8*N/10+dd,:))
        %set(h2,'Cdata',V2(N/2-r/2+1:N/2-r/2+dd,:))
        %set(h3,'Cdata',V1(N/10+1:N/10+dd,:))
        %set(h4,'Cdata',V1(N/10+1:N/10+dd,:))
        %set(h5,'Cdata',V2(N/2-r/2+1:N/2-r/2+dd,:))
        %set(h6,'Cdata',V2(N/2-r/2+1:N/2-r/2+dd,:))
        %set(h7,'Cdata',V2(N/2-r/2+1:N/2-r/2+dd,:))
        refreshdata
        pause(0.1)
        hold on
        saveas(gcf, ['1rfma',sprintf('%04d',p),'.jpg'])
        p=p+1;
    end
    if max(xpos1)>M | max(xpos2)>M
        last=i
        break
    end
    close all
end
for i=1:8
    Vtb=U+sin(50*i*dt*0.5*omega+last*dt*0.5*omega)*Vtbold;
    Vlr=-U+sin(50*i*dt*0.5*omega+pi+last*dt*0.5*omega)*Vtbold;
    V1=Vtb*ones(N,M);
    V2=Vlr*ones(N,M);
    Y1=Y1-[0;0;0;0;Y1(5);Y1(6)];
    Y1=Y1+[50;0;0;0;0;0];
     h0=mesh(x,8*N/10+1:8*N/10+dd,actop(8*N/10+1:8*N/10+dd,:),V1(8*N/10+1:8*N/10+dd,:))
view(114,8);
caxis manual
caxis([-Vtbold Vtbold]);
colormap(winter)
colorbar
axis manual
axis([-500 M+500 0 N 0 N -Vtbold Vtbold]);
hold on
axis off
        h1=mesh(x,8*N/10+1:8*N/10+dd,acbottom(8*N/10+1:8*N/10+dd,:),V1(8*N/10+1:8*N/10+dd,:))
    h2=mesh(x,N/2-r/2+1:N/2-r/2+dd,btop(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h3=mesh(x,N/10+1:N/10+dd,acbottom(N/10+1:N/10+dd,:),V1(N/10+1:N/10+dd,:))
    h4=mesh(x,N/10+1:N/10+dd,actop(N/10+1:N/10+dd,:),V1(N/10+1:N/10+dd,:))
    h5=mesh(x,N/2-r/2+1:N/2-r/2+dd,bbottom(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h6=mesh(x,N/2-r/2+1:N/2-r/2+dd,dtop(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    h7=mesh(x,N/2-r/2+1:N/2-r/2+dd,dbottom(N/2-r/2+1:N/2-r/2+dd,:),V2(N/2-r/2+1:N/2-r/2+dd,:))
    refreshdata
        pause(0.1)
        hold on
 
            xpos1(f)=Y1(1);
            ypos1(f)=Y1(2);
            zpos1(f)=Y1(3);
            f=f+1;
            Plot1=plot3(xpos1,ypos1,zpos1);
            hold on
            if f>6
                xpos11=xpos1;
                ypos11=ypos1;
                zpos11=zpos1;
                for j=1:5
                    xpos1(j)=xpos11(j+1);
                    ypos1(j)=ypos11(j+1);
                    zpos1(j)=zpos11(j+1);                    
                end
                f=f-1;
        end
           
                   saveas(gcf, ['1rfma',sprintf('%04d',p),'.jpg'])
        p=p+1;
        close all
end
