PENGGUNAAN METODE RUNGE-KUTTA (RK4) PADA. SIMULASI LINTASAN BERKAS PROTON DALAM SIKLOTRON PET 13 MeV

Volume 15, Oktober 2013

ISSN 1411-1349

PENGGUNAAN METODE RUNGE-KUTTA (RK4) PADA SIMULASI LINTASAN BERKAS PROTON DALAM SIKLOTRON PET 13 MeV Pramudita Anggraita, Emy Mulyani Pusat Teknologi Akselerator dan Proses Bahan (PTAPB), Badan Tenaga Nuklir Nasional (BATAN) Jl. Babarsari, Kotak Pos 6101 ykbb Yogyakarta 55281 E-mail:[email protected]

ABSTRAK PENGGUNAAN METODE RUNGE-KUTTA (RK4) PADA SIMULASI LINTASAN BERKAS PROTON DALAM SIKLOTRON PET 13 MeV. Laju perubahan momentum proton dalam siklotron merupakan persamaan diferensial order satu dalam waktu sebagai fungsi gaya Lorentz yang dibangkitkan oleh medan listrik dan medan magnet dalam siklotron. Persamaan diferensial order satu ini dapat diselesaikan secara numerik dengan metode Runge-Kutta (RK4) dengan selang waktu dt yang tetap. Dibandingkan dengan penyelesaian secara numerik sebelumnya yaitu dengan selang jarak ds yang tetap, ternyata RK4 memberikan hasil yang lebih baik terutama pada linearitas kenaikan tenaga. Kata Kunci:Runge-Kutta (RK4), simulasi, lintasan proton, siklotron PET

ABSTRACT APPLICATION OF RUNGE-KUTTA METHOD (RK4) IN PROTON BEAM TRAJECTORY SIMULATION IN A 13 MeV PET CYCLOTRON. The rate of change of proton momentum in a cyclotron is an ordinary or firstorder differential equation in time as a function of Lorenz force generated by electric and magnetic fields in the cyclotron. This ordinary differential equation can be solved numerically by using Runge-Kutta (RK4) method with constant time increment dt. Comparison with previous solution by using constant spatial increment ds, RK4 method yields better solution especially in the linearity of the energy increase. Keywords:Runge-Kutta (RK4), simulation, proton trajectory, PET cyclotron

PENDAHULUAN

P

TAPB (Pusat Teknologi Akselerator dan Proses Bahan) – BATAN mempunyai tugas dan fungsi mengembangkan teknologi akselerator, salah satunya adalah pengembangan siklotron 13 MeV untuk PET (DECY-13, Development of Experimental Cyclotron in Yogyakarta) yang direncanakan selesai tahun 2019[1]. Hal yang penting dalam disain siklotron adalah mengetahui bagaimana lintasan berkas proton sejak dari sumber ion hingga mengenai target. Lintasan berkas proton tersebut tergantung pada distribusi medan listrik pemercepat dan medan magnet di dalam siklotron, yang dapat dihitung dengan perangkat lunak Opera3D (dalam 3 dimensi) dengan modul Tosca (untuk medan magnet dan listrik) dan modul Soprano atau RELAX3D (untuk medan listrik). Simulasi lintasan dilakukan dengan program iterasi menggunakan perangkat lunak Scilab, yang merupakan perangkat lunak open source serupa Mathlab (komersial). Hasil simulasi lintasan dapat digunakan untuk memperkirakan apakah disain

PENGGUNAAN METODE RUNGE-KUTTA (RK4) PADA SIMULASI LINTASAN BERKAS PROTON DALAM. SIKLOTRON PET 13 MeV Pramudita Anggraita, Emy Mulyani

siklotron yang dibuat dapat berfungsi seperti yang diinginkan, maupun untuk memperbaiki disain yang telah dibuat. Simulasi lintasan berkas elektron dalam akselerator di PTAPB pernah dilakukan untuk akselerator searah (dc) linear [2] dengan perhitungan distribusi medan listrik (statis) menggunakan iterasi dengan metode SOR (Simultaneous Over Relaxation [3] atau Succesive Over Relaxation[4]). Dengan perangkat lunak Opera3D dan modul Tosca/Soprano perhitungan distribusi medan listrik ini juga dapat dilakukan. Simulasi lintasan berkas proton sudah dilakukan dengan perangkat lunak Scilab 5.3.3 pada tahun 2012 untuk lintasan berkas H− hingga mencapai tenaga sekitar 13 MeV. Simulasi lintasan dilakukan dengan menggunakan persamaan gerak F = dp/dt = mγdv/dt, F gaya Lorentz, p momentum, m massa, γ koreksi relativistik, t waktu, kecepatan v = ds/dt = vo + dv, dan posisi x = xo + vdt. Selang jarak ds dibuat tetap sesuai ketelitian data simulasi medan yaitu 1 mm. Hasil simulasi ternyata bentuk lintasan dan kurva kenaikan tenaga belum teratur (Gambar 1)[5].

1

Volume 15, Oktober 2013

(a)

ISSN 1411-1349

(b)

(c) −

Gambar 1. Hasil simulasi lintasan berkas H hingga mencapai tenaga sekitar 13 MeV[5]. (a) lintasan (hijau) dan posisi fase berkas (hitam), (b) kenaikan tenaga yang tidak linear terhadap putaran (biru), garis merah kenaikan tenaga maksimum, (c) simpangan vertikal (dalam mm). Simulasi lintasan juga dikembangkan untuk disain siklotron PET di Korea, yaitu dengan program pwheel [6] yang menggunakan penyelesaian numerik Runge-Kutta (RK4). Metode ini akan diterapkan pada program simulasi lintasan yang dikembangkan dengan menggunakan Scilab 5.4.1.

TEORI Perhitungan lintasan berkas ion muatan q kecepatan v didasarkan pada gaya Lorentz F yang diterima oleh partikel bermuatan dalam medan listrik E dan medan magnet B. Medan listrik E merupakan fungsi posisi (x,y,z) dan waktu (t), sedang medan magnet B merupakan fungsi posisi saja. Besar kedua medan ini dapat diperoleh dari simulasi perhitungan maupun dari pengukuran (mapping). Secara komponen persamaan gerak dapat dituliskan sebagai dpx/dt = Fx = q(Ex + vyBz − vzBy), dpy/dt = Fy = q(Ey + vzBx − vxBz), dpz/dt = Fz = q(Ez + vxBy − vyBx).

(1)

Penyelesaian masing-masing komponen persamaan (1) menggunakan metode Runge-Kutta orde-4 (RK4) adalah [7] (i = x,y,z) pi(to + dt) = pio + (ki1 + 2ki2 + 2ki3 + ki4)/6

(2)

dengan ki1 = Fi(to, po)dt, ki2 = Fi(to + dt/2, po + ki1/2)dt, ki3 = Fi(to + dt/2, po + ki2/2)dt, ki4 = Fi(to + dt, po + ki3)dt.

(3)

Setelah momentum baru pi diperoleh maka kecepatan vi = pi/γm dan posisi baru xi = xio + vidt dapat ditentukan. Prosiding Pertemuan dan Presentasi Ilmiah Teknologi Akselerator dan Aplikasinya Vol. 15 , Oktober 2013 : 1 - 8

TATA KERJA Program lintasan berkas H− dibuat dengan Scilab 5.4.1 seperti terlihat pada Lampiran III. Program tersebut memerlukan data kuat medan magnet yang dimasukkan dengan program loadB pada Lampiran I. Karena simetri, data kuat medan magnet ini cukup mencakup 1/8 dari seluruh kawasan medan yang berukuran 480 mm × 480 mm × 30 mm, yang diperoleh dari perhitungan dengan menggunakan perangkat lunak Opera3d dan modul Tosca [8]. Besar data binary sekitar 38 MB untuk masing-masing komponen Bx, By, dan Bz. Data potensial listrik dengan potensial dee 40 kV terhadap ground juga dihitung dengan perangkat lunak Opera3d dan modul Tosca dimasukkan dengan program loadV12 dan loadV34 pada Lampiran II. Karena tidak simetris data meliputi seluruh kawasan medan yang berukuran 480 mm × 480 mm × 30 mm. Supaya ukuran data tidak terlalu besar, data potensial dibagi menjadi 4 kuadran, yaitu V1, V2 (dimasukkan dengan loadV12), V3, dan V4 (dimasukkan dengan load V34). Besar data binary sekitar 56 MB untuk masing-masing V1, V2, V3, dan V4. Medan listrik Ex, Ey, dan Ez untuk masing-masing kuadran dihitung dari beda potensial pada ds = 1 mm ke arah masing-masing komponen medan. Berkas lintasan dimulai pada posisi lubang sumber ion x = 17,5 mm; y = 9 mm; z = 0. Kecepatan awal berkas diberikan berdasarkan kenaikan tenaga berkas setelah melintas 0,1 mm dari posisi awal, dengan arah gerak sesuai arah medan listrik pada posisi awal. Gerak, tenaga, dan posisi selanjutnya dihitung dengan menggunakan persamaan (1), (2), dan (3). Lintasan (hijau) diplot pada bidang XY, tenaga dan posisi vertikal diplot sebagai fungsi putaran. Karena siklotron bekerja pada harmonik 4, periode putaran sama dengan 4 periode frekuensi medan pemercepat. Posisi berkas pada tiap periode putaran 2

Volume 15, Oktober 2013

ISSN 1411-1349

ini diplot pada lintasan berkas pada bidang XY, sehingga menunjukkan ketepatan fase berkas pada tiap putaran. Karena data medan magnet dan potensial listrik dalam ketelitian mm, maka nilai medan untuk tiap posisi dibulatkan ke nilai mm terdekat. Dalam perhitungan lintasan dengan pwheel digunakan interpolasi bspline 5-titik dan data potensial listrik mempunyai ketelitian 0,25 mm; 0,5 mm; 1 mm; 2 mm; dan 4 mm dalam ukuran 256×256×5 [6].

HASIL DAN PEMBAHASAN Dengan variasi frekuensi medan pemercepat dan fase awal berkas, lintasan dan capaian tenaga terbaik diperoleh pada frekuensi sekitar 77,64 MHz dan fase awal −50°. Selang waktu yang optimum untuk perhitungan dengan RK4 diperoleh dt = periode medan pemercepat/100. Hasil lintasan, kenaikan tenaga, dan simpangan vertikal dapat dilihat pada Gambar 2. Terlihat bahwa lintasan (hijau) lebih simetris, fase berkas (hitam) mendekati stabil untuk tiap putaran (a). Kenaikan tenaga linear dan 13 MeV dicapai pada sekitar putaran ke 89 (b) dan simpangan vertikal (z) mendekati stabil (c).

(c) simpangan vertikal (mm) Gambar 2.

Hasil simulasi lintasan berkas H− hingga mencapai tenaga sekitar 13 MeV dengan RK4.

Variasi fase awal memberikan rentang fase sebesar 16º (fase awal antara −37º dan −53º) yang memberikan lintasan dan kenaikan tenaga cukup baik hingga 13 MeV [9]. Sebagai perbandingan rentang fase pada siklotron KIRAMS-13 adalah sekitar 20º [10], sehingga desain DECY-13 masih perlu untuk dioptimasikan lebih lanjut.

KESIMPULAN Penggunaan RK4 dalam simulasi lintasan berkas dalam siklotron DECY-13 telah dapat memberikan lintasan, fase berkas, dan linearitas kenaikan tenaga yang lebih baik dibandingkan dengan simulasi sebelumnya. Rentang fase 16º yang diperoleh perlu ditingkatkan dengan optimasi desain sistem pemercepat yang lebih baik.

UCAPAN TERIMA KASIH (a) lintasan (hijau) dan posisi fase berkas (hitam)

Salah satu penulis (Emy Mulyani) mengucapkan terima kasih kepada Kementerian Riset dan Teknologi yang telah memberikan beasiswa program S2 (20112013) sehingga dapat melaksanakan sebagian tugas yang dilaporkan pada makalah ini.

DAFTAR PUSTAKA [1] ANONIM, Renstra BATAN 2010 – 2014, Badan Tenaga Atom Nasional 2010.

(b) kenaikan tenaga yang linear terhadap putaran (biru), garis merah kenaikan tenaga maksimum PENGGUNAAN METODE RUNGE-KUTTA (RK4) PADA SIMULASI LINTASAN BERKAS PROTON DALAM. SIKLOTRON PET 13 MeV Pramudita Anggraita, Emy Mulyani

[2] PRAMUDITA ANGGRAITA, dkk., Simulasi Lintasan Berkas Elektron dalam Tabung Pemercepat Mesin Berkas Elektron, Prosiding Pertemuan dan Presentasi Ilmiah PPNY-BATAN, Yogyakarta 26-27 Mei 1998, hal. 450 – 458.

3

Volume 15, Oktober 2013

[3] HERMANTO, A., dkk., Solusi Numerik Persamaan Poisson Tiga Dimensi dengan Komputer Pribadi, Laporan Penelitian FMIPAUGM, 1994. [4] BERNARD, M.Y., Particles and Fields: Fundamental Equations, dalam Focusing of Charged Particles, Vol. I, ed. A. Septier, Academic Press Inc., NY, 1967. [5] PRAMUDITA ANGGRAITA dkk., Optimasi Sistem Pemercepat Siklotron Pet 13 MeV dengan Simulasi Lintasan Berkas, Laporan Teknis, 2012. [6] YOON, M., Central Region Design and Beam Dynamics, BATAN Accelerator School, Yogyakarta, 2012. [7] LANDAU, R.H., PAEZ, M.J., BORDEIANU, C.C., A Survey Computational Physics, Princeton University Press, 2010. [8] PRAMUDITA ANGGRAITA dkk., Beam Tracking Simulation in the Central Region of a 13 MeV PET Cyclotron, AIP Conf. Proc. 1454, 178, 2012. [9] EMY MULYANI, Desain Central Region Siklotron Proton 13 MeV, tesis Program Studi S2 Ilmu Fisika Universitas Gadjah Mada, Agustus 2013. [10] AN, D.H., dkk., Beam Trajectory Simulations for KIRAMS-13 Cyclotron, Lab. Of Accelerator Development, KIRAMS, Seoul, Korea, 2004. http://epaper.kek.jp/c04/data/ CYC2004 _papers/20P05.pdf

ISSN 1411-1349

TANYA JAWAB Silakhuddin − Apakah program komputer yang disajikan sama dengan program pwheel? − Jika tidak sama, apakah kelebihan program ini dibandingkan dengan pwheel? Pramudita Anggraita − Sama-sama memakai RK4, pwheel juga memakai bspline, program Scilab belum − pwheel memakai resolusi medan ¼, ½, 1, 2 mm, program Scilab hanya memakai resolusi 1 mm. Bambang Siswanto − Dengan menggunakan RK4 linearitas energinya sudah diperoleh, namun tidak mungkin berimpit dengan kondisi maksimum. Jadi berapa derajat maksimum yang diperbolehkan? − Dalam penyimpangan tersebut masih dikatakan memenuhi syarat? Pramudita Anggraita − Secara teori diberikan oleh faktor waktu transit sin(g/r)/(g/r), g = celah antara tepi dee dan tepi liner, r = ruji lintasan. Jika g ≈ rθ, faktor waktu transit ~ 0,99 pada sudut celah 30º (JJ Livingstone, “Principle of Cyclic Particle Accelerators,” hal. 124125).

LAMPIRAN I: Program loadB tic(); clear all stacksize('max'); //medan magnet (tesla) load('F:\simulasi\Bx133040.dat','B1'); load('F:\simulasi\By133040.dat','B2'); load('F:\simulasi\Bz133040.dat','B3'); funcprot(0);tl=toc;//8s

LAMPIRAN II: Program loadV12 dan loadV34 clear all stacksize('max') //medan listrik kuadran 3,4 (volt) load('F:\simulasi\V1gap30.dat','V1'); load('F:\simulasi\V2gap30.dat','V2'); funcprot(0)

Prosiding Pertemuan dan Presentasi Ilmiah Teknologi Akselerator dan Aplikasinya Vol. 15 , Oktober 2013 : 1 - 8

4

Volume 15, Oktober 2013

ISSN 1411-1349

clear all stacksize('max') //medan listrik kuadran 3,4 (volt) load('F:\simulasi\V3gap30.dat','V3'); load('F:\simulasi\V4gap30.dat','V4');

LAMPIRAN III: Program lintasan berkas dengan menggunakan RK4 clear x1;clear y1;clear z1;clear t2;clear T1;clear X;clear Y;clear F; stacksize('max') e=-1.60217933*10^-19;mp=1.672631*10^-27;me=9.1093897*10^-31; c=2.997925*10^8;a=1000;ds=0.001; t=0;bx=0;by=0;//waktu awal, pusat B f=-50;f1=f*%pi/180;//fase awal (derajat) n=40000;//banyaknya titik lintasan, 84s ff=77.64;omf=2*%pi*ff*10^6;//frekuensi RF harmonik 4 (MHz),frek sudut (Hz) mp=mp+2*me;//massa proton -> massa Hx=0.0175;y=0.009;z=-0.000;//posisi awal (m), kuadran 1 x1(1)=x*a;y1(1)=y*a;z1(1)=z*a;//posisi awal (mm) xr=round(x*a);yr=round(y*a);zr=round(z*a);//indeks awal Ex=-(V1(xr+2,yr+1,zr+16)-V1(xr+1,yr+1,zr+16))/ds;//medan listrik Ey=-(V1(xr+1,yr+2,zr+16)-V1(xr+1,yr+1,zr+16))/ds;//pada posisi Ez=-(V1(xr+1,yr+1,zr+17)-V1(xr+1,yr+1,zr+16))/ds;//sumber ion E=(Ex^2+Ey^2+Ez^2)^(0.5);T=E*ds/10;T1(1)=T/10^6;//T awal (MeV) v=abs((2*e*T/(mp))^(0.5))//perkiraan kecepatan awal berkas gam=mp/(1-v^2/c^2)^0.5;p=gam*v;//pusa relativistik px=-p*Ex/E;py=-p*Ey/E;pz=-p*Ez/E;//pusa awal searah E (negatif) function [fx, fy, fz]=fp(Ex, Ey, Ez, Bx, By, Bz, px1, py1, pz1, t1, dt) s=sin(omf*t1+f1); Ex=Ex*s; Ey=Ey*s; Ez=Ez*s;//medan listrik fx=e*(Ex+(py1*Bz-pz1*By)/gam); fy=e*(Ey+(pz1*Bx-px1*Bz)/gam); fz=e*(Ez+(px1*By-py1*Bx)/gam);//F(t,p) endfunction N=0;T1(2)=0; for i=2:n; xr=round(x*a);yr=round(y*a);zr=round(z*a);//indeks awal medan listrik if (x>=0&y>=0) then Ex=-(V1(xr+2,yr+1,zr+16)-V1(xr+1,yr+1,zr+16))/ds; Ey=-(V1(xr+1,yr+2,zr+16)-V1(xr+1,yr+1,zr+16))/ds; Ez=-(V1(xr+1,yr+1,zr+17)-V1(xr+1,yr+1,zr+16))/ds;end if (x<=0&y>=0) then Ex=-(V2(xr+481,yr+1,zr+16)-V2(xr+480,yr+1,zr+16))/ds; Ey=-(V2(xr+481,yr+2,zr+16)-V2(xr+481,yr+1,zr+16))/ds; Ez=-(V2(xr+481,yr+1,zr+17)-V2(xr+481,yr+1,zr+16))/ds;end if (x<=0&y<=0) then Ex=-(V3(xr+481,yr+481,zr+16)-V3(xr+480,yr+481,zr+16))/ds; Ey=-(V3(xr+481,yr+481,zr+16)-V3(xr+481,yr+480,zr+16))/ds; Ez=-(V3(xr+481,yr+481,zr+17)-V3(xr+481,yr+481,zr+16))/ds;end if (x>=0&y<=0) then Ex=-(V4(xr+2,yr+481,zr+16)-V4(xr+1,yr+481,zr+16))/ds; Ey=-(V4(xr+1,yr+481,zr+16)-V4(xr+1,yr+480,zr+16))/ds; Ez=-(V4(xr+1,yr+481,zr+17)-V4(xr+1,yr+481,zr+16))/ds;end //medan magnet Bx=B1(abs(xr)+1+bx,abs(yr)+1+by,abs(zr)+1)*sign(x)*sign(z); PENGGUNAAN METODE RUNGE-KUTTA (RK4) PADA SIMULASI LINTASAN BERKAS PROTON DALAM. SIKLOTRON PET 13 MeV Pramudita Anggraita, Emy Mulyani

5

Volume 15, Oktober 2013

ISSN 1411-1349

By=B2(abs(xr)+1+bx,abs(yr)+1+by,abs(zr)+1)*sign(y)*sign(z); Bz=B3(abs(xr)+1+bx,abs(yr)+1+by,abs(zr)+1); dt=1/(100*ff*10^6);//interval waktu t1=t;px1=px;py1=py;pz1=pz; [fx,fy,fz]=fp(Ex,Ey,Ez,Bx,By,Bz,px1,py1,pz1,t1,dt) k1x=dt*fx;k1y=dt*fy;k1z=dt*fz;//k1 t1=t+dt/2;px1=px+k1x/2;py1=py+k1y/2;pz1=pz+k1z/2; [fx,fy,fz]=fp(Ex,Ey,Ez,Bx,By,Bz,px1,py1,pz1,t1,dt) k2x=dt*fx;k2y=dt*fy;k2z=dt*fz;//k2 t1=t+dt/2;px1=px+k2x/2;py1=py+k2y/2;pz1=pz+k2z/2; [fx,fy,fz]=fp(Ex,Ey,Ez,Bx,By,Bz,px1,py1,pz1,t1,dt) k3x=dt*fx;k3y=dt*fy;k3z=dt*fz;//k3 t1=t+dt;px1=px+k3x;py1=py+k3y;pz1=pz+k3z; [fx,fy,fz]=fp(Ex,Ey,Ez,Bx,By,Bz,px1,py1,pz1,t1,dt) k4x=dt*fx;k4y=dt*fy;k4z=dt*fz;//k4 px=px+(k1x+2*k2x+2*k3x+k4x)/6;vx=px/gam;x=x+vx*dt;//pusa & posisi baru py=py+(k1y+2*k2y+2*k3y+k4y)/6;vy=py/gam;y=y+vy*dt; pz=pz+(k1z+2*k2z+2*k3z+k4z)/6;vz=pz/gam;z=z+vz*dt;t=t+dt; x1(i)=1000*x; y1(i)=1000*y;z1(i)=1000*z;//x,y,z in mm v=(vx^2+vy^2+vz^2)^0.5;ga=(1-v^2/c^2)^-0.5;gam=ga*mp;//faktor & massa relativistik T1(i)=(ga-1)*mp*c^2/(-e*10^6);t2(i)=t*10^6*ff/4;//T(MeV), turn number E(i)=(Ex^2+Ey^2+Ez^2)^0.5;F(i)=E(i)*sin(omf*t+f1);//E dalam V/m tn=(omf*t+f1-%pi/2)/(8*%pi);Z=abs(round(tn)-tn); if (Z<0.001) then N=N+1;X(N)=x1(i);Y(N)=y1(i);end end subplot(221) a1=100; isoview(-a1,a1,-a1,a1); plot(x1,y1,'g');//lintasan berkas //dees xd=[238.19 49.49 49.49 46.57 40.38 4.06 -25.63 -79.97 -92.48 -98.99 -98.99 -440.63 -238.19 -49.49 -49.49 -48.39 41.06 0 5.89 6.19 6.5 6.98 7.33 7.88 8.35 8.98 9.7 10.42 11.08 12.01 12.62 12.7 13.14 13.61 14.21 14.6 15.07 15.9 16.34 16.81 17.4 17.95 18.52 18.91 19.36 20.65 21.07 21.71 22.08 22.41 22.93 23.45 23.97 24.46 24.88 25.35 25.93 25.54 27.26 79.34 90.18 98.99 98.99 440.63 238.19]; yd=[440.63 98.99 98.99 92.02 81.25 15.73 -13.96 -41.65 -46.12 -49.49 -49.49 -238.19 -440.63 -98.99 -98.99 -90.19 -80.57 0 5.89 8.14 9.36 10.68 11.45 12.63 13.23 14.09 14.49 15.67 16.26 16.98 17.43 17.4 17.73 17.97 18.27 18.47 18.64 18.94 19.06 19.21 19.31 19.43 19.56 19.6 19.6 19.68 19.66 19.61 19.6 19.53 19.4 19.33 19.26 19.11 18.96 18.82 18.57 18.32 17.99 42.28 48.42 49.49 49.49 238.19 440.63]; //CR1 (liner1) xl1=[37.48 85.91 88.08 89.4 89.9 89.55 87.33 83.38 37.48]; yl1=[12.02 26.83 18.71 10.02 1.18 -10.16 -22.53 -33.88 12.02]; //CR2 (liner2) xl2=[2.82 56.08 49.8 41.59 33.09 21.16 7.9 -2.84 -11.03 -19.45 -27.73 2.82]; yl2=[-16.97 -70.39 -74.96 -79.8 -83.69 -87.47 -89.65 -89.96 -89.32 -87.87 -85.62 -16.97]; //CR3 (liner3) xl3=[-85.9 -29.86 -0.17 27.69 16.69 10 -4.66 -18.78 -32.29 -46.69 -63.64 -81.33 -86.67 -89.8 -89.31 -85.9]; yl3=[-26.85 -9.72 19.98 85.63 88.44 89.44 89.88 80.03 84.0 79.94 63.64 38.54 24.26 5.99 -11.10 -26.85]; //puller1 xp1=[17.2015 19.6259 21.0322 20.5516 20.1726 19.6259 16.7646 12.3956 11.6674 9.6750 9.6331 9.8082 10.000 12.7369 15.0974 15.3635 16.2817 16.5194 17.2015]; yp1=[14.8436 15.4869 18.6743 19.6768 20.4040 21.3728 25.3966 29.5743 29.4466 25.8522 25.5377 25.2447 25.0867 22.4119 19.2011 18.7572 17.0297 16.5194 14.8436]; //puller2 xp2=[7.9374 8.1022 8.3605 14.4747 14.8161 14.8116 13.7709 13.5216 12.4517 7.4622 7.0391 6.7400 4.5023 4.5864 7.4960 7.9374]; yp2=[12.5772 12.4982 12.4977 14.1209 14.4329 14.8116 16.2941 16.5418 17.2823 20.7358 20.8051 20.5664 16.5150 15.9197 13.0100 12.5772]; xbg=[0.6244 1.2104 1.6306 3.0962 3.1629 2.6709 2.3812 1.6306 1.2160 -1.0714 -4.0202 -5.7096 -8.9661 -12.4021 14.8518 -16.3992 -19.2422 -21.2550 -22.3784 -23.2210 -23.3658 -23.3191 -23.2874 -22.9997 -22.2926 -21.6288 Prosiding Pertemuan dan Presentasi Ilmiah Teknologi Akselerator dan Aplikasinya Vol. 15 , Oktober 2013 : 1 - 8

6

Volume 15, Oktober 2013

ISSN 1411-1349

18.6787 -18.4703 -18.3867 -18.3011 -18.0923 -17.7306 -16.7467 -15.2209 -12.1216 -10.0508 -7.3575 -4.1200 1.7426 0.2481 0.6244]; ybg=[18.2922 18.4007 18.8233 21.6388 22.4013 22.9876 23.0856 23.2033 23.2570 23.4109 23.257 22.9876 22.0722 20.4659 18.8233 17.5059 14.3070 10.8465 7.8536 3.9722 2.5974 0.1012 -0.8303 -1.4354 -1.7283 -1.4763 1.4713 1.7762 2.1359 3.2631 4.6438 6.1439 8.7307 11.3149 14.6078 16.0489 17.3361 18.1947 18.4100 18.3331 18.2922]; //beam guide 1 xbg1=[4.90 5.68 6.23 7.52 7.57 7.01 3.54 -0.05 -1.49 -2.42 -3.19 -3.64 -3.74 -3.85 -3.21 -1.91 -0.73 0.61 1.86 3.20 4.14 4.90]; ybg1=[27.26 27.23 27.73 30.44 31.22 31.79 33.08 34.02 33.96 33.51 32.65 31.42 30.49 29.57 29.49 29.31 29.03 28.73 28.36 27.92 27.57 27.26]; //beam guide 2 xbg2=[-9.02 -6.75 -8.71 -10.35 -12.08 -13.49 -14.70 -15.83 -19.54 -18.23 -16.78 -15.03 -13.25 -11.17 -9.02]; ybg2=[-17.69 -13.24 -12.16 -11.00 -9.49 -7.86 -6.17 -4.16 -7.87 -9.94 -11.83 -13.59 -15.11 -16.51 -17.69]; //head sumber ion xhs=[13.79 15.18 17.73 19.34 21.61 22.84 24.48 25.18 24.04 22.84 21.13 18.56 16.56 14.92 13.54 12.69 12.74 13.14 13.79]; yhs=[7.99 8.35 8.98 9.39 9.97 9.42 7.25 4.46 0.41 -0.67 -1.64 -2.04 -1.53 -0.49 0.98 3.38 5.83 7.03 7.99]; //batang sumber ion xts=[23.41 14.43 360.64 369.63 23.41]; yts=[8.74 -0.24 -346.46 -337.48 8.74]; //liner xl1=[84.76 443.64 443.64 188.26 188.26 30.23 30.23 84.76 84.76]; yl1=[30.26 188.43 188.43 443.71 443.71 84.77 84.77 30.26 30.26]; xl2=[-30.23 -30.23 -84.76 -84.76 -443.64 -443.64 -188.27 -188.27 -30.23]; yl2=[-84.77 -84.77 -30.26 -30.26 -188.43 -188.43 -443.71 -443.71 -84.77]; xl3=[-85.90 -86.80 -88.44 -89.26 -89.67 -90.08 -89.26 -88.85 -87.21 -85.16 -83.12 -81.07 -79.84 -78.61 -74.51 72.47 -69.19 -63.05 -56.50 -51.17 -46.67 -42.58 -38.07 -32.33 -26.61 -20.05 -12.68 -3.67 4.52 6.97 11.88 20.48 23.83 27.66 25.75 8.16 5.22 3.15 -1.86 -18.38 -29.88 -85.90]; yl3=[-26.85 -23.49 -18.59 -12.45 -5.91 -1.81 8.82 12.92 22.33 30.51 36.64 38.28 42.37 44.83 50.15 53.01 58.74 64.88 70.19 73.47 76.74 79.60 81.24 83.69 86.15 87.38 88.61 89.42 90.24 89.43 88.61 87.79 86.68 85.4 80.93 39.41 32.6 27.8 18.36 1.82 -9.73 -26.85]; //hills xh1=[-27.06 -69.16 -95.57 -106.07 -117.03 -128.4 -141.83 -145.25 -164.54 -183.67 -145.07 -94.18 -51.93 -1.04 51.93 94.18 145.07 183.67 164.54 145.25 141.83 128.4 117.03 106.07 95.57 69.16 27.06 20.12 11.92 0 -11.92 20.12 -27.06]; yh1=[85.83 219.35 289.83 312.49 335.15 357.65 379.36 383.41 397.27 443.47 455.67 470.08 476.8 479.68 476.8 470.08 455.67 443.47 397.27 383.41 379.36 357.65 335.15 312.49 289.83 219.35 85.83 87.72 89.21 90 89.21 87.72 85.83] xh2=[85.83 219.35 289.83 312.49 335.15 357.65 379.36 383.41 397.27 443.47 455.67 470.08 476.8 479.68 476.8 470.08 455.67 443.47 397.27 383.41 379.36 357.65 335.15 312.49 289.83 219.35 85.83 87.72 89.21 90 89.21 87.72 85.83]; yh2=[-27.06 -69.16 -95.57 -106.07 -117.03 -128.4 -141.83 -145.25 -164.54 -183.67 -145.07 -94.18 -51.93 -1.04 51.93 94.18 145.07 183.67 164.54 145.25 141.83 128.4 117.03 106.07 95.57 69.16 27.06 20.12 11.92 0 -11.92 20.12 -27.06]; xh3=[-27.06 -69.16 -95.57 -106.07 -117.03 -128.4 -141.83 -145.25 -164.54 -183.67 -145.07 -94.18 -51.93 -1.04 51.93 94.18 145.07 183.67 164.54 145.25 141.83 128.4 117.03 106.07 95.57 69.16 27.06 20.12 11.92 0 -11.92 20.12 -27.06]; yh3=[-85.83 -219.35 -289.83 -312.49 -335.15 -357.65 -379.36 -383.41 -397.27 -443.47 -455.67 -470.08 -476.8 479.68 -476.8 -470.08 -455.67 -443.47 -397.27 -383.41 -379.36 -357.65 -335.15 -312.49 -289.83 -219.35 -85.83 87.72 -89.21 -90 -89.21 -87.72 -85.83]; xh4=[-85.83 -219.35 -289.83 -312.49 -335.15 -357.65 -379.36 -383.41 -397.27 -443.47 -455.67 -470.08 -476.8 479.68 -476.8 -470.08 -455.67 -443.47 -397.27 -383.41 -379.36 -357.65 -335.15 -312.49 -289.83 -219.35 -85.83 87.72 -89.21 -90 -89.21 -87.72 -85.83]; yh4=[-27.06 -69.16 -95.57 -106.07 -117.03 -128.4 -141.83 -145.25 -164.54 -183.67 -145.07 -94.18 -51.93 -1.04 51.93 94.18 145.07 183.67 164.54 145.25 141.83 128.4 117.03 106.07 95.57 69.16 27.06 20.12 11.92 0 -11.92 20.12 -27.06]; //Central pole/bump PENGGUNAAN METODE RUNGE-KUTTA (RK4) PADA SIMULASI LINTASAN BERKAS PROTON DALAM. SIKLOTRON PET 13 MeV Pramudita Anggraita, Emy Mulyani

7

Volume 15, Oktober 2013

ISSN 1411-1349

xcp=[27.69 16.69 10 -4.66 -18.78 -32.29 -46.69 -63.64 -81.33 -86.67 -89.8 -89.31 -85.9 -82.43 -76.77 -69.62 -59.51 -49.2 -41.05 -32.74 -27.73 -11.03 -2.84 7.9 21.16 33.09 41.59 49.8 56.08 62.97 70.05 76.76 83.38 87.33 89.55 89.9 89.4 88.08 85.91 82.53 79.34 72.56 63.64 52.64 40.38 32.29 27.69]; ycp=[85.63 88.44 89.44 89.88 80.03 84.0 79.94 63.64 38.54 24.26 5.99 -11.10 -26.85 -36.14 -46.97 -57.03 -67.52 75.36 -80.57 -83.84 -85.62 -89.32 -89.96 -89.65 -87.47 -83.69 -79.8 -74.96 -70.39 -64.29 -55.98 -46.98 -33.88 22.53 -10.16 1.18 10.02 18.71 26.83 35.90 42.28 53.24 63.64 73 81.24 84 85.63]; xb=[90 87.74 81.09 70.36 56.11 39.05 20.03 0 -20.03 -39.05 -56.11 -70.36 -81.09 -87.74 -90 -87.74 -81.09 -70.36 56.11 -39.05 -20.03 0 20.03 39.05 56.11 70.36 81.09 87.74 90]; yb=[0 -20.03 -39.05 -56.11 -70.36 -81.09 -87.74 -90 -87.74 -81.09 -70.36 -56.11 -39.05 -20.03 0 20.03 39.05 56.11 70.36 81.09 87.74 90 87.74 81.09 70.36 56.11 39.05 20.03 0]; plot(xbg2,ybg2,'b'); plot(xd,yd,'b',xl1,yl1,'r',xl2,yl2,'r',xl3,yl3,'r');plot(xbg,ybg,'b') plot(xb,yb,'m');plot(xts,yts,'b');plot(xhs,yhs,'b'); plot(xh1,yh1,'m',xh2,yh2,'m',xh3,yh3,'m',xh4,yh4,'m'); plot(xl1,yl1,'b',xl2,yl2,'b',xl3,yl3,'b');plot(xbg1,ybg1,'b'); plot(xp1,yp1,'b',xp2,yp2,'b');plot(X,Y,'black'); title("H- trajectory (mm) RK4"); a=gca(); // Handle on axes entity a.x_location = "origin"; a.y_location = "origin"; subplot(224) plot(t2,z1,'b'); xtitle('vertical movement', 'turn number', 'z (mm)'); subplot(223) plot(t2,T1,'b',t2,t2*0.16,'r');//T vs putaran, T maks xtitle( 'Energy (MeV)', 'turn number', 'E (MeV)'); subplot(222) plot(t2,F/10^6,'b'); xtitle('Electric Field Strength', 'turn number', 'E (MV/m)'); funcprot(0)

Prosiding Pertemuan dan Presentasi Ilmiah Teknologi Akselerator dan Aplikasinya Vol. 15 , Oktober 2013 : 1 - 8

8