Commit 4a340b16 authored by Thomas Planche's avatar Thomas Planche

finished the TOSCA section of the tutorial, now working on the MAP2D-E section

parent f2f64fc7
set term postscript eps color enhanced "Times-Roman" 18 lw 2
set xlabel 'Z [cm]'
set ylabel 'E [MV/cm??]'
set grid
titletext="Testing EMULT vs MAP2D-E"
labeltext="(c) `whoami`, "
today="`date +%Y/%b/%d`"
set label labeltext.today at screen .01, screen .02
set title titletext
#set xrange[-22:22]
set output "plots/Ez_vs_z.eps"
plot "zgoubi.plt" u 22:37 w lp lc 1
set output "plots/Ex_vs_z.eps"
plot "zgoubi.plt" u 22:38 w lp lc 1
set output "plots/Ey_vs_z.eps"
plot "zgoubi.plt" u 22:39 w lp lc 1
** TEST MULTIPOLE
'OBJET' 1
546.198 !60 keV U1+ = 546.198
1 !1:multi particles on a grid
1 1 1 1 1 1
0.3 1.0 0.3 0.0 0.0 1.
0.3 0.0 0.3 0.0 0.0 1.
'PARTICUL' 2
22172.3 1.602176487E-19 0.0 0.0 0.0
! 'OPTIONS'
! 1 1
! CONSTY ON
'MAP2D-E' 5
1 2 !print the map (no print = 0); output along partile(s) trajectory(ies) = 2
0.00001 0.1 0.1!Normalization factor for E (to be put in MV/cm), X (cm), and Y(cm)
HEADER_0 !this actually only works if you have n=0 header lines. It crashes if n>0 and I could not figure out why
45 27 !Number of longitudinal and horizontal-transverse nodes of the mesh (the Z elevation is arbitrary)
test.map !file name
0 0. 0. 0.
4
1.
1 0. 0. 0.
'DRIFT'
20.0
'FAISTORE' 3
zgoubi_map.fai
1
'END'
** TEST MULTIPOLE
'OBJET' 1
546.198 !60 keV U1+ = 546.198
1 !1:multi particles on a grid
11 1 11 1 1 1
0.3 1.0 0.3 0.0 0.0 1.
0.0 0.0 0.0 0.0 0.0 1.
'PARTICUL' 2
22172.3 1.602176487E-19 0.0 0.0 0.0
'ELMULT' 10
2 !IL=7 for output in zgoubi.impdev.out
10. 3. 0. 1000. 20000. 0. 0. 0. 0. 0. 0. 0. ! Length[cm]; tip radius[cm]; dipole tip field[V/m]; quad[V/m]; sextu[V/m]; etc...
20. 10. 1. 1. 1. 1. 1. 1. 1. 1. 1. ! ENTRANCE FACE: Integration zone[cm]; Dipole fringe field extend (FFE)[cm]; quad FFE relative to dipole FFE [unitless]; sextu FFE/dipole FFE [unitless]; etc...
4 0. 1.8 0. 0. 0. 0. !unused; C0...C5
20. 10. 1. 1. 1. 1. 1. 1. 1. 1. 1. ! EXIT FACE: Integration zone[cm]; Dipole fringe field extend (FFE)[cm]; quad FFE relative to dipole FFE [unitless]; sextu FFE/dipole FFE [unitless]; etc...
4 0. 1.8 0. 0. 0. 0. !unused; C0...C5
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. !skew angles[rad] of each component
1.0 !Integration step size[cm]
1 0. 0. 0. !KPOS (1=no misalagnment ); XCE[cm]; YCE[cm]; ALE[rad]
'DRIFT'
40.0
'FAISTORE' 3
zgoubi_emult.fai
1
'END'
set term postscript eps color enhanced "Times-Roman" 18 lw 2
set xlabel 'X [cm]'
set ylabel 'E [MV/cm??]'
set grid
titletext="Testing EMULT vs MAP2D-E"
labeltext="(c) `whoami`, "
today="`date +%Y/%b/%d`"
set label labeltext.today at screen .01, screen .02
set title titletext
set key bot left
set output "plots/Ex_vs_x.eps"
plot "ELMULT.plt" u 22:37 w p lc 1 ,\
"MAP2D.plt" u 22:37 w p lc 3
set key top right
set output "plots/Ey_vs_x.eps"
plot "ELMULT.plt" u 22:38 w p lc 1 ,\
"MAP2D.plt" u 22:38 w p lc 3
set output "plots/Ez_vs_x.eps"
plot "ELMULT.plt" u 22:39 w p lc 1 ,\
"MAP2D.plt" u 22:39 w p lc 3
** TEST MULTIPOLE
'OBJET' 1
546.198 !60 keV U1+ = 546.198
1 !1:multi particles on a grid
11 1 11 1 1 1
0.3 1.0 0.3 0.0 0.0 1.
0.0 0.0 0.0 0.0 0.0 1.
'PARTICUL' 2
22172.3 1.602176487E-19 0.0 0.0 0.0
! 'OPTIONS'
! 1 1
! CONSTY ON
'MAP2D-E' 5
1 2 !print the map (no print = 0); output along partile(s) trajectory(ies) = 2
1. 1. 1.
HEADER_0
51 15
impdev2EMap.out
0 0. 0. 0.
2
1.
1 0. 0. 0.
'DRIFT'
20.0
'FAISTORE' 3
zgoubi_map.fai
1
'END'
** TEST MULTIPOLE
'OBJET' 1
546.198 !60 keV U1+ = 546.198
1 !1:multi particles on a grid
11 1 1 1 1 1
0.3 1.0 0.3 0.0 0.0 1.
0.0 0.0 0.0 0.0 0.0 1.
'PARTICUL' 2
22172.3 1.602176487E-19 0.0 0.0 0.0
'ELMULT' 10
2 !IL=7 for output in zgoubi.impdev.out
10. 3. 0. 1000. 20000. 0. 0. 0. 0. 0. 0. 0. ! Length[cm]; tip radius[cm]; dipole tip field[V/m]; quad[V/m]; sextu[V/m]; etc...
20. 10. 1. 1. 1. 1. 1. 1. 1. 1. 1. ! ENTRANCE FACE: Integration zone[cm]; Dipole fringe field extend (FFE)[cm]; quad FFE relative to dipole FFE [unitless]; sextu FFE/dipole FFE [unitless]; etc...
4 0. 1.8 0. 0. 0. 0. !unused; C0...C5
20. 10. 1. 1. 1. 1. 1. 1. 1. 1. 1. ! EXIT FACE: Integration zone[cm]; Dipole fringe field extend (FFE)[cm]; quad FFE relative to dipole FFE [unitless]; sextu FFE/dipole FFE [unitless]; etc...
4 0. 1.8 0. 0. 0. 0. !unused; C0...C5
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. !skew angles[rad] of each component
1.0 !Integration step size[cm]
1 0. 0. 0. !KPOS (1=no misalagnment ); XCE[cm]; YCE[cm]; ALE[rad]
'DRIFT'
40.0
! 'FAISTORE' 3
! forFrancois.fai
! 1
'END'
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
set term postscript eps color enhanced "Times-Roman" 18 lw 2
set xlabel 'X [cm]'
set ylabel 'E [MV/cm??]'
set grid
titletext="Testing EMULT with fringe field"
labeltext="(c) `whoami`, "
today="`date +%Y/%b/%d`"
set label labeltext.today at screen .01, screen .02
set title titletext
set key bot left
set output "Ex_vs_x.eps"
plot "zgoubi.plt" u 22:37 w p lc 1
set key top right
set output "Ey_vs_x.eps"
plot "zgoubi.plt" u 22:38 w p lc 1
set output "Ez_vs_x.eps"
plot "zgoubi.plt" u 22:39 w p lc 1
** HSR no quad using DIPOLE
'OBJET' 1
546.198 !60 keV U1+ = 546.198
1 !1:multi particles on a grid
15 1 1 1 1 1
0.3 1.0 0.0 0.0 0.0 1.
0.0 0.0 0.0 0.0 0.0 1.
'PARTICUL' 2 2
22172.3 1.602176487E-19 0.0 0.0 0.0
'OPTIONS'
1 1
CONSTY ON
'ELMULT' 10
7 !IL=7 for output in zgoubi.impdev.out
10. 3. 0. 1000. 20000. 0. 0. 0. 0. 0. 0. 0. ! Length[cm]; tip radius[cm]; dipole tip field[V/m]; quad[V/m]; sextu[V/m]; etc...
20. 10. 1. 1. 1. 1. 1. 1. 1. 1. 1. ! ENTRANCE FACE: Integration zone[cm]; Dipole fringe field extend (FFE)[cm]; quad FFE relative to dipole FFE [unitless]; sextu FFE/dipole FFE [unitless]; etc...
4 0. 1.8 0. 0. 0. 0. !unused; C0...C5
20. 10. 1. 1. 1. 1. 1. 1. 1. 1. 1. ! EXIT FACE: Integration zone[cm]; Dipole fringe field extend (FFE)[cm]; quad FFE relative to dipole FFE [unitless]; sextu FFE/dipole FFE [unitless]; etc...
4 0. 1.8 0. 0. 0. 0. !unused; C0...C5
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. !skew angles[rad] of each component
1.0 !Integration step size[cm]
1 0. 0. 0. !KPOS (1=no misalagnment ); XCE[cm]; YCE[cm]; ALE[rad]
! 'FAISTORE' 3
! zgoubi.fai
! 1
'END'
implicit double precision (a-h,o-z)
character(200) txt200
parameter (mxs = 600, mxy = 400, mxz = 100)
! dimension ss(mxs), yy(mxy), zz(mxz)
dimension ex(mxs,mxy,mxz),ey(mxs,mxy,mxz),ez(mxs,mxy,mxz)
dimension ss(mxs,mxy,mxz),yy(mxs,mxy,mxz),zz(mxs,mxy,mxz)
parameter (rad = (4.d0 * atan(1.d0)) /180.d0)
INTEGER NMAX
PARAMETER(NMAX=411) !number of value per line in zgoubi.impdev.out
REAL TTT(NMAX)
open(unit=1,file='zgoubi.impdev.out')
open(unit=2,file='impdev2EMap.out')
C read 2-line header
read(1,fmt='(a)')
read(1,fmt='(a)')
is = 51 ! number of mesh nodes, longitudinal i.e. nb of steps
jy = 15 ! MUST BE ODD: number of mesh nodes, transverse i.e nb of test particles
kz = 1 ! MUST BE ODD: number of mesh nodes, vertical =1 (2D)
write(*,*) 'Number of mesh nodes in X/Y/Z : ',is,'/',jy,'/',kz
! write(*,*) 'type enter to continue'
! read(*,*)
write(*,*) 'Ok, now busy, making map ... wait'
write(*,*) 'Map will be in impdev2FieldMap.out'
if(is .gt. mxs) stop ' Too many angles'
if(jy .gt. mxy) stop ' Too many radii'
if(kz .gt. mxz) stop ' Too many Zs'
exmi = 1d10
exma = -1d10
eymi = 1d10
eyma = -1d10
ezmi = 1d10
ezma = -1d10
do jj = 1, jy
do kk = 1, kz
do i = 1, is
read(1,*,err=10,end=10) (TTT(III),III=1,NMAX)
if(i.eq.is) read(1,*,err=10,end=10) ! because integr writes twice at the last integration step
yy(i,jj,kk) = TTT(2)
zz(i,jj,kk) = TTT(4)
ss(i,jj,kk) = TTT(6)
ex(i,jj,kk) = TTT(371)
ey(i,jj,kk) = TTT(372)
ez(i,jj,kk) = TTT(373)
if(exmi .gt. ex(i,jj,kk)) exmi = ex(i,jj,kk)
if(exma .lt. ex(i,jj,kk)) exma = ex(i,jj,kk)
if(eymi .gt. ey(i,jj,kk)) eymi = ey(i,jj,kk)
if(eyma .lt. ey(i,jj,kk)) eyma = ey(i,jj,kk)
if(ezmi .gt. ez(i,jj,kk)) ezmi = ez(i,jj,kk)
if(ezma .lt. ez(i,jj,kk)) ezma = ez(i,jj,kk)
enddo
enddo
enddo
10 continue
write(*,*) ' Done reading ! '
write(*,*) ' Min - max field components : '
write(*,*) ' ex : ',exmi,'/',exma
write(*,*) ' ey : ',eymi,'/',eyma
write(*,*) ' ez : ',ezmi,'/',ezma
do j = 1, jy
if(j .le. (jy+1)/2) then
jj = jy-2*(j-1)
else
jj=2*(j-(jy+1)/2)
end if
! PRINT *,jj,ex(is/2,jj,1)
do k = 1, kz
if(k .le. (kz+1)/2) then !this part not really tested yet (only tried with kz=1, i.e 2D)
kk = kz-2*(k-1)
else
kk=2*(k-(kz+1)/2)
end if
do i = 1, is
write(2,fmt='(1p,f11.5,1x,
> 0p,2(f11.5,1x),
> 1p,3(e18.10,1x),3(i5,1x))') yy(i,jj,kk)/10.,zz(i,jj,kk)/10., !/10 because y and z come in mm, and we want them in cm
> ss(i,jj,kk),
> ey(i,jj,kk), ez(i,jj,kk), ex(i,jj,kk),j,k,i
enddo
enddo
enddo
stop
end
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
set term postscript eps color enhanced "Times-Roman" 18 lw 2
set output "plots/thinBeam.eps"
set xlabel 'Y [cm]'
set ylabel 'Y` [mrad]'
set grid
titletext="Testing EMULT"
labeltext="(c) `whoami`, "
today="`date +%Y/%b/%d`"
set label labeltext.today at screen .01, screen .02
set title titletext
set key top left
cm2mm=10.0
#set xrange[-80:80]
#set yrange[-60:60]
plot "zgoubi_emult.fai" u 10:11 with points,\
"zgoubi_map.fai" u 10:11 with points lc 3
set output "plots/thinBeam_z.eps"
set xlabel 'Z [cm]'
set ylabel 'Z` [mrad]'
plot "zgoubi_emult.fai" u 12:13 with points,\
"zgoubi_map.fai" u 12:13 with points lc 3
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
set term postscript eps color enhanced "Times-Roman" 18 lw 2
set output "plots/thinBeam.eps"
set xlabel 'Y [cm]'
set ylabel 'Y` [mrad]'
set grid
titletext="Testing EMULT"
labeltext="(c) `whoami`, "
today="`date +%Y/%b/%d`"
set label labeltext.today at screen .01, screen .02
set title titletext
set key top left
cm2mm=10.0
#set xrange[-80:80]
#set yrange[-60:60]
plot "zgoubi_emult.fai" u 10:11 with points,\
"zgoubi_map.fai" u 10:11 with points lc 3
set output "plots/thinBeam_z.eps"
set xlabel 'Z [cm]'
set ylabel 'Z` [mrad]'
plot "zgoubi_emult.fai" u 12:13 with points,\
"zgoubi_map.fai" u 12:13 with points lc 3
This diff is collapsed.
......@@ -3,9 +3,7 @@ $const #z_step 10.
$const #x_width 130.*2
$const #z_width 220*2
GRID X0=-#x_width/2. Y0=0 Z0=-#z_width/2. DXG=#x_step DYG=1 DZG=#z_step NXG=INT(#x_width/#x_step)+1 NYG=1 NZG=INT(#z_width/#z_step)+1 FILE='/home/tplanche/text/designs/HRS/opera/multipole/from_Carla/test.table' BINARY=NO FORMAT=2 F1=X UNIT1G=LENGU F2=Y UNIT2G=LENGU F3=Z UNIT3G=LENGU F4=Ex UNIT4G=1 F5=Ey UNIT5G=1 F6=Ez UNIT6G=1 F7=INT((X+#x_width/2.)/#x_step) UNIT7G=1 F8=1 UNIT8G=1 F9=INT((Z+#z_width/2.)/#z_step) UNIT9G=1 UNIT10G=1 UNIT11G=1 UNIT12G=1
y
$const #a INT(#z_width/#z_step)+1
$const #a INT(#x_width/#x_step)+1
$const #print_nb_step_z INT(#z_width/#z_step)+1
$const #print_nb_step_x INT(#x_width/#x_step)+1
This diff is collapsed.
set term postscript eps color enhanced "Times-Roman" 18 lw 2
set output "plots/BzXplot.eps"
set xlabel 'Theta [deg.]'
set ylabel 'Bz [T]'
set grid
titletext="HRS TOSCA reference trajectory"
labeltext="(c) `whoami`, "
today="`date +%Y/%b/%d`"
set label labeltext.today at screen .01, screen .02
set title titletext
set key top left
cm2mm=10.0
plot "zgoubi.plt" using ($22/pi*180.):($25/10.) with p pt 0
This diff is collapsed.
This diff is collapsed.
** HSR using TOSCA map
'OBJET' 1
546.198 !60 keV U1+ = 546.198
1
1 121 1 1 1 1
0.0 1.0 0.0 0.0 0.0 1.
0.0 0.0 0.0 0.0 0.0 1.
! 'OBJET' 1
! 546.198 !60 keV U1+ = 546.198
! 5 !5: generate 11 particles used to calculate transfer matrix ('MATRIX')
! 0.01 0.01 0.01 0.01 .1 .001 !step size in Y; T; Z; P; S; D
! 0.0 0. 0. 0. 0. 1.
'PARTICUL' 2
22172.3 1.602176487E-19 0.0 0.0 0.0
'DRIFT' 3
25.500
'TOSCA' 4
0 0 !print the map (no print = 0); output along partile(s) trajectory(ies) = 2
1.0365 1. 1. 1. !Magnetic field scaling; X coordinate scaling; Y; Z
HEADER_2 ! Title
281 65 1 22. ! Nb of nodes in Theta; R; Z (=1 for 2D map); MOD(see manual, =0 for Cartesian, with mid-plane symmetry)
Dipole2D.map !File name
0 0 0 0. ! ID (see manual);
4 !IORDRE(=2, 4 or 25 if 2D, unused if 3D)
0.1 !Integration step size[cm]
2 !KPOS, normally=2
132.5586547 -0.436332313 132.5586547 0.436332313 !RE has been fitted to be =RS; TE=-25deg(in rad) whic is the distance between the effective edge and the edge of th efield map; RS=RE; TS=-TE
'DRIFT' 5
2.0430 !24.0430 - 22.0 (which is the half length of the multipole field map)
'MAP2D-E' 6
0 0 !print the map (no print = 0); output along partile(s) trajectory(ies) = 2
0.000001 0.1 0.1!Normalization factor for E (to be put in MV/cm), X (cm), and Y(cm)
HEADER_0 !this actually only works if you have n=0 header lines. It crashes if n>0 and I could not figure out why
45 27 !Number of longitudinal and horizontal-transverse nodes of the mesh (the Z elevation is arbitrary)
Multipole2D.map !file name
0 0. 0. 0.
4
1.
1 0. 0. 0.
'DRIFT' 7
2.0430 !24.0430 - 22.0 (which is the half length of the multipole field map)
'TOSCA' 8
0 0 !print the map (no print = 0); output along partile(s) trajectory(ies) = 2
1.0365 1. 1. 1. !Magnetic field scaling; X coordinate scaling; Y; Z
HEADER_2 ! Title
281 65 1 22. ! Nb of nodes in Theta; R; Z (=1 for 2D map); MOD(see manual, =0 for Cartesian, with mid-plane symmetry)
Dipole2D.map !File name
0 0 0 0. ! ID (see manual);
4 !IORDRE(=2, 4 or 25 if 2D, unused if 3D)
0.1 !Integration step size[cm]
2 !KPOS, normally=2
132.5586547 -0.436332313 132.5586547 0.436332313 !RE has been fitted to be =RS; TE=-25deg(in rad) whic is the distance between the effective edge and the edge of th efield map; RS=RE; TS=-TE
'DRIFT' 9
25.500
'FAISTORE' 10
zgoubi.fai
1
! 'MATRIX' 9
! 1 0 !order of the matrix/map; 0:means calculate the matrix here
! 'FIT' 10
! 1 !Number of physical parameters to be varied
! 2 1 7.001 0.5
! 2
! 1 1 2 9 0.0 1. 0 !transfer matrix;
! 1 2 1 9 0.0 1. 0 !transfer matrix;
'END' 11
THREED TYPE=SURFACE VECTOR=NO XORIGIN=0 YORIGIN=1200 ZORIGIN=106 ROTX=1.0E-04 ROTY=179.999 ROTZ=1.0E-04 XASPECT=1 YASPECT=1 ZASPECT=1 SIZE=704.188791082898 FACETANGLE=10 PERSPECTIVE=NO LINECOLOUR=YES OPTION=SETVIEW
/++ GRID X0=0 Y0=500 Z0=0 DXG=20 DYG=20 DZG=1 NXG=101 NYG=51 NZG=1 FILE='forZgoubi' BINARY=NO FORMAT=2 F1=Y UNIT1G=LENGU F2=Z UNIT2G=LENGU F3=X UNIT3G=LENGU F4=0 UNIT4G=FLUXU F5=Bz UNIT5G=FLUXU F6=0 UNIT6G=FLUXU UNIT7G=1 UNIT8G=1 UNIT9G=1 UNIT10G=1 UNIT11G=1 UNIT12G=1
/++ yes
//global WCS
SET XLOCAL=0 YLOCAL=0 ZLOCAL=0 TLOCAL=0,PLOCAL=0,SLOCAL=0 LOOK=ANYWHERE ABORT=YES
$const #mm2cm 0.1
$const #T2kG 10.
$const #pi ATAN(1)*4
$const #rmin 104. |//cm
$const #rmax 136. |//cm good field region : 120+/-16 cm
$const #Dr 1.0 |//cm
$const #Dtheta 0.5 |//deg
$const #Dz 0 |//cm
//open output file
$open 1 forZgoubi_MOD22.dat overwrite
$format 1 fixed 16 8
$format 2 character 16
/++ R0,DR,DTTA,DZ
$assign 1 1 1 1
$write 1 #rmin #Dr #Dtheta #Dz
$assign 2 2 2 2 2 2
$write 1 'Y[cm]' 'Z[cm]' 'X[cm]' 'Br[kG]' 'Bth[kG]' 'Bz[kG]'
$assign 1 1 1 1 1 1
/++ $do #r 900 1500 20
$do #r #rmin/#mm2cm #rmax/#mm2cm #Dr/#mm2cm
$do #th -25 115 #Dtheta
$const #x #r*SIND(-#th+45.0)
$const #y #r*COSD(-#th+45.0)
$const #z 0
POINT XP=#x YP=#y ZP=0 COMP=bz
$write 1 #y*#mm2cm #z*#mm2cm #x*#mm2cm 0.0 BZ*#T2kG 0.0
/++ $write 1 #r*#mm2cm #th*#pi/180. #z*#mm2cm 0.0 0.0 BZ*#T2kG
$end do
$end do