-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbessel.pro
126 lines (103 loc) · 2.93 KB
/
bessel.pro
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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
nx=40L
ny=41L
nz=42L
x=findgen(nx)+0.0
y=findgen(ny)+0.0
z=findgen(nz)+0.0
xx=rebin(reform(x,nx,1,1),nx,ny,nz )/1.0/nx-0.5+1e-6
yy=rebin(reform(y,1,ny,1),nx,ny,nz)/1.0/ny-0.5+1e-6
zz=rebin(reform(z,1,1,nz),nx,ny,nz )/1.0/nz-0.5+1e-6
rr=sqrt(xx^2 + yy^2)
zz=zz
theta=atan(yy,xx )
ar=fltarr(nx,ny,nz)
atheta=fltarr(nx,ny,nz)
az=fltarr(nx,ny,nz)
help,ar
pitchangle=0.01
for i=0,nx-1 do begin
for j=0,ny-1 do begin
for k=0,nz-1 do begin
for n=1,5 do begin
a=2.0
scrh=n*pitchangle
r=rr[i,j,k]
z=zz[i,j,k]
thet=theta[i,j,k]
nka=scrh*a
nkr=scrh*rr[i,j,k]
phi=0.25
if (r lt a) then begin
KnkaInkrplus1=beselk(nka,n+1)*beseli(nkr,n+1)
ar[i,j,k] -= (beselk(nka,n+1)*beseli(nkr,n+1)-beselk(nka,n-1)*beseli(nkr,n-1))*sin ( n*( thet -phi - pitchangle*z) )
atheta[i,j,k] += (beselk(nka,n+1)*beseli(nkr,n+1)+beselk(nka,n-1)*beseli(nkr,n-1))*cos ( n*( thet -phi - pitchangle*z) )
az[i,j,k] += (beselk(nka,n)*beseli(nkr,n))*cos ( n*( thet -phi - pitchangle*z) )
endif else begin
ar[i,j,k] -= (beseli(nka,n+1)*beselk(nkr,n+1)-beseli(nka,n-1)*beselk(nkr,n-1))*sin ( n*( thet -phi - pitchangle*z) )
atheta[i,j,k] += (beseli(nka,n+1)*beselk(nkr,n+1)+beseli(nka,n-1)*beselk(nkr,n-1))*cos ( n*( thet -phi - pitchangle*z) )
az[i,j,k] += (beseli(nka,n)*beselk(nkr,n))*cos ( n*( thet -phi - pitchangle*z) )
endelse
;print, i,j,k
endfor
if (r lt a) then begin
ar[i,j,k] *= pitchangle*a
atheta[i,j,k] *= pitchangle*a
atheta[i,j,k] += pitchangle*rr[i,j,k]
az[i,j,k] -= alog (a)
endif else begin
ar[i,j,k] *= pitchangle*a
atheta[i,j,k] *= pitchangle*a
atheta[i,j,k] += pitchangle/rr[i,j,k]*a*a
az[i,j,k] -= alog (r)
endelse
endfor
endfor
endfor
ax=ar*cos(theta) -sin(theta) *atheta
ay=ar*sin(theta) +cos(theta) *atheta
curl, ax,ay,ax, bx,by,bz
plot9, ar,atheta,az
;plot9, bx,by,bz
j=1
;------------------------------------------------------------------------
;f = rf(j,datadir='./') ; change directory if needed
;------------------------------------------------------------------------
openw,2,'fields_electrons_'+STRTRIM(j,2)+'.vtk',/SWAP_ENDIAN
;B-field
grid0 = 0
grid1 = nx
grid3 = ny
grid5 = nz
grid_tot = (grid1)*(grid3)*(grid5)
spa1 = 1.
spa2 = 1.
spa3 = 1.
print,'Writing the binary data of snapshot',j
printf,2,'# vtk DataFile Version 2.0'
printf,2,'PIC'
printf,2,'BINARY'
printf,2,'DATASET STRUCTURED_POINTS'
printf,2,'DIMENSIONS '+STRTRIM(grid1)+STRTRIM(grid3)+STRTRIM(grid5)
printf,2,'ORIGIN 0 0 0'
printf,2,'SPACING '+STRTRIM(spa1)+STRTRIM(spa2)+STRTRIM(spa3)
printf,2,'POINT_DATA '+STRTRIM(grid_tot)
printf,2,'SCALARS Vx FLOAT'
printf,2,'LOOKUP_TABLE default'
writeu,2,ar
printf,2,'SCALARS Vy FLOAT'
printf,2,'LOOKUP_TABLE default'
writeu,2,atheta
printf,2,'SCALARS Vz FLOAT'
printf,2,'LOOKUP_TABLE default'
writeu,2,az
printf,2,'SCALARS Bx FLOAT'
printf,2,'LOOKUP_TABLE default'
writeu,2,bx
printf,2,'SCALARS By FLOAT'
printf,2,'LOOKUP_TABLE default'
writeu,2,by
printf,2,'SCALARS Bz FLOAT'
printf,2,'LOOKUP_TABLE default'
writeu,2,bz
close,2
end