-
Notifications
You must be signed in to change notification settings - Fork 0
/
Transitorio.m
90 lines (74 loc) · 1.64 KB
/
Transitorio.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
clear all
close all
clc
colordef white
N=1e5;
R=1;
v=0.2665;
x=-1+2*rand(N,1);
y=-1+2*rand(N,1);
z=-1+2*rand(N,1);
vx=v*randn(N,1);
vy=v*randn(N,1);
vz=v*randn(N,1);
dentro=x.^2+y.^2+z.^2<=R^2;
N0=sum(dentro);
x=x(dentro);
y=y(dentro);
z=z(dentro);
vx=vx(dentro);
vy=vy(dentro);
vz=vz(dentro);
dt=0.1;
t=[0:dt:50];
% fig=figure;
% set(fig,'color','k');
% mov = avifile('en_distr.avi','fps',25,'quality',100,'compression','Indeo5');
for it=1:length(t)
t(it);
r2=x.^2+y.^2+z.^2;
[r2,ord]=sort(r2);
x=x(ord);
y=y(ord);
z=z(ord);
vx=vx(ord);
vy=vy(ord);
vz=vz(ord);
r=sqrt(r2);
N=((0:1:N0-1)+0.5)';
E=(min(1,r).^3-N/N0)./r2;
vx=vx-dt*E.*x./r;
vy=vy-dt*E.*y./r;
vz=vz-dt*E.*z./r;
x=x+vx*dt;
y=y+vy*dt;
z=z+vz*dt;
phi(N0)=0;
for j=N0-1:-1:1
phi(j)=phi(j+1)+0.5*(E(j)+E(j+1))*(r(j+1)-r(j));
end
Energia=(vx.^2+vy.^2+vz.^2)/2-phi';
plot(r,Energia,'.r',r,-phi,'Linewidth',1.5,'MarkerSize',1)
xlabel('x/R','fontsize',16)
grid on
axis([0 12 -.8 0.1])
drawnow
% F=getframe(gcf);
% mov=addframe(mov,F);
interno=x.^2+y.^2+z.^2<=R^2;
int(it)=sum(interno);
end
figure (2)
plot(r,-phi,'Linewidth',3)
axis([0 12 -0.5 0.1])
grid on
xlabel('x/R','fontsize',16)
ylabel('Potential','fontsize',16)
figure (3)
plot(t,int/N0*100,'Linewidth',2)
axis([0 50 0 100])
hold on
plot([0 t(length(t))],[100 100],'r','Linewidth',2)
xlabel('t\omega','fontsize',16)
ylabel('N/N_0 (%)','fontsize',16)
legend('Electrons','Protons')