C C C C BUNCH A nonrelativistic longitudinal optics program C R.C. Pardo, ANL C C C*********************************************************************** C C C Original Version : R. Pardo, ANL, written in BASIC, C received at UK Dec 1988. C C Modified : J. Vanhoy, Univ. of KY, C 1) changed to FORTRAN, Jan '89 C 2) transit time effects in magnets, Feb '89 C C C*********************************************************************** C C BUNCH is a program which allows the calculation of longitudinal C optics and is especially suited for studies of beam bunching for C arbitrary paths with acceleration, drifts, and bunching modules C of various types. Transverse optics effects are totally ignored C in these calculations at present, so path length differences in C magnets are not treated (in original version < Feb 1989). C C The input is either from a file or from the keyboard. The first C input from the keyboard requests a file name. If a carriage C return is executed from the console, then console input is C assumed. If a file name is given, the file is opened, and the C commands are read from the file as if they had been manually C entered from the console. When the program reads an "E" command C in the file the file is closed and control is returned to the C console for further work or final exiting of program. C C Every command executed and overall beam parameters at that point C are echoed to an output file called "OUT.DAT". The file can be C printed at any time to get a record of the session. The activity C is appended to the existing file, therefore the user will want to C delete that file occasionally, or rename it to save an important C result. C C The first portion of input after the input file name is a C description of the beam parameters. The following sequence is C required: C I. 0 or 1 controls distribution of particles C 0 - produces a uniform distribution of C particles in both time and energy C 1 - produces a gaussian distribution in C energy but a uniform distribution in time C such as one would expect form a DC ion C source. C II. U0, M, Q, DU C U0 - the particle initial energy in keV, C M - the particle mass in AMU, C Q - the particles charge in proton charge C units, and C DU - the energy spread of the distribution. C III. BPRD - the beam bunch period in nanoseconds. C This parameters along with the energy spread C DU determines the initial phase space of the C particle ensemble. Now a thousand particles C in the phase space defined by DU and BPRD C will be produces as required by the C distribution function selected above. C C Now a sequence of commands can be entered which control output C results and the evolution of this bunch of particles. Each C command has certain control data which may be necessary for the C implementation of that command. The commands which are C understood by the program at the present and described below C along with the required input parameters for that command. C C A Acceleration region. C VACEL, LENGTH - total accelerating C voltage, in keV, and length, in meters, C of acceleration region. This is a C uniform acceleration region such as an C accelerator tube from an ion source C high voltage platform. C STAB - stability of voltage system in C deltaV/V. C C B Buncher. C TYPE - Buncher type. Presently three C are allowed. C TYPE = 1 sinewave buncher C TYPE = 2 Sawtooth waveform C TYPE = 3 linear bunching - velocity C variation is linear in time. Voltage is C quadratic. C BVOLTAGE, PERIOD - Maximum buncher C voltage in kV, period of buncher in C nanoseconds. C RELPHS - phase,in radians, of this C buncher relative to center of central C particle. Used only for TYPE =1. C DIST - Distance to waist, in meters, or C point of interest, can be left zero and C "D" command used after this step. C TWAIST - time to waist, in nanoseconds, C used only for type = 3. C BSTAB - This parameter allows one to C study the stability requirements for C the buncher amplitude. The value, in C kVolts, is range on total stability to C assume. The 'correct voltage' is C varied randomly, sine wieghted, by the C amount BSTAB. C C C Comment line C String$ - an ASCII string of your C chosing which is echoed to the output C file C C D Drift region. C DIST - Drift distance in meters. C C E End of input file data. No parameters. C This command closes the input file and C returns control to the console. C C G Graph results. C Display a two dimensional phase space C plot of the beam ensemble at this point C along the beamline. C DEVICE - pgplot output device, assumed C to be '/re' unless specified on line C with G command. C C H Histogram energy and time projections. C No parameters. Display projections of C time and energy spectra from the phase C plots. C C I Integrate phase space. No parameters. C Attempts to find area of phase space C occupied by particles. Doesn't work C very well. C C L Print log of results so far. No C parameters. C This command forces a listing of the C file OUT.DAT at this point. OUT.DAT is C essentially a journal file of C everything that has happened up to this C point. C C M Analyzing magnet with 90 deg bend. C Takes into account transit time differences C arising from finite beam size and C energy spread in beam. C C P Print. No parameters. C Print a statement of the beam C parameters at this point in the C process. The general parameters of the C ensemble are printed and the histograms C projected onto the energy and time axis C are also printed. C C R Restart. No parameters. C Start an entirely new case. All past C results are destroyed. C C S Change graph scale. C EMAX - The maximum delta energy, in C keV, to display in subsequent graph C displays. C TMAX - The maximum delta time, in nsec., C to display. C C T Input graph title. C STRING$ - ASCII string to be used as C title for graphs. C C Q Quit. No parameters. C End program. C C N New case, same file. Start new case C without closing the input file. Same as C "R" command from console. C C*********************************************************************** C dimension tbins(50),ebins(50) dimension tf(1000),ef(1000),vf(1000) real xpos(50),mass,m,l,hafmas character*3 device integer*4 iseed character*9 cdate character*8 ctime character*7 cc,cf character*20 fspec,clabel character*40 comment common/proj/tbins,ebins,ytmax,yemax,tf,ef,max, 1 tcenter,ecenter,tfwhm,efwhm C iseed=011159 p25 = .125 p50 = .25 p100 = .5 p200 = 1. p400 = 2. 5 manual = 0 in=1 iout=2 write(*,10) 10 format(' enter input file name: ',$) read(*,15) fspec 15 format(a20) if(fspec.eq.' ') then in=5 iout=6 go to 20 endif open(unit=1,file=fspec,type='old') 20 open(unit=2,file='out.dat',type='new') 25 write(iout,30) 30 format(' ENTER 1 FOR GAUSSIAN OR 0 FOR UNIFORM DIST ',$) read(in,*) distr continue ! cls max = 1000 c c initialize system c l = 0. tzero = 0. v = 0. p = 0. twaist = 0. prev = 0. s = 0. stab = 0. bl = 0. gflag = 1. c write(iout,35) 35 format(' ENTER BEAM UZERO(keV);MASS;Q;DU(keV) ',$) read(in,*) u0,m,q,du write(iout,40) 40 format(' ENTER BEAM BUNCH PERIOD (ns.) ',$) read(in,*) bprd write(2,45) fspec 45 format(' INPUT FILE NAME: ',a20) call getdate(cdate) call gettime(ctime) write(2,*) cdate,' ',ctime phi = 0. csq = .09 c = .3 mass = 931507. * m vzero = c * sqrt(2 * u0 / mass) beta0 = vzero / c write(2,50) m,q,u0,beta0,du 50 format(' ION M;Q;U0;BETA;DU: ',5F10.4) write(2,55) bprd 55 format(' BEAM BUNCH PERIOD: ',F10.3,'ns.') go to 385 ! generate randomized beam 60 continue write(iout,*) ' bunching parameters ' beta0 = vzero / c tarea = du * bprd ! pbrd write(iout,65) ecenter,efwhm 65 format(' DELTA ENERGY CENTROID= ',f7.3, 1 ' keV FWHM= ',f7.3,' keV') write(iout,70) tcenter,tfwhm 70 format(' DELTA TIME CENTROID= ',f7.3, 1 ' ns FWHM= ',f7.3,' ns') write(iout,75) m,q,u0,beta0,du 75 format(' ION M;Q;U0(keV);BETA;DU(keV): ',5f10.3) write(iout,80) v,p,twaist,bstab 80 format(' GRID VOLTAGE, PERIOD, TWAIST, STABILITY ',4f10.3) write(iout,85) prev,s,stab 85 format(' ACCEL. VOLTAGE, LENGTH, STABILITY: ',3f10.3) write(iout,90) l,tzero 90 format(' TOTAL DISTANCE,TRANSIT TIME: ',2f10.3) c c command section c 95 continue write(2,100) l,tzero 100 format(' TOTAL DISTANCE: ',f10.3,' TRANSIT TIME: ',f10.3) continue write(iout,105) 105 format(' COMMAND (a,b,d,g,h,l,m,p,q,r,s,i) : ',$) read(in,110) chk,device 110 format(a1,1x,a3) if(device.eq.' ') device='/re' write(2,115) chk 115 format(' COMMAND: ',a1) if(chk.eq.'c') go to 420 !C Comment line in input file Echo to file if(chk.eq.'e') go to 415 !E End of input file data if(chk.eq.'g') then call graph(tf,tmax,ef,emax,max,device) go to 60 !G Graph results endif if(chk.eq.'l') go to 430 !L PRINT LOG OF RESULT SO FAR if(chk.eq.'p') go to 195 !P Print if(chk.eq.'r') then close(1) close(2) go to 5 !R Restart endif if(chk.eq.'s') go to 405 !S Change graph scale if(chk.eq.'q') go to 120 !Q Quit if(chk.eq.'a') go to 140 !A Accelerator tube if(chk.eq.'b') go to 260 !B Buncher if(chk.eq.'d') go to 175 !D Drift region if(chk.eq.'h') go to 365 !H Histogram energy and time projection if(chk.eq.'i') go to 125 !I Integrate phase space if(chk.eq.'n') go to 25 !N New case, same file if(chk.eq.'m') go to 440 !M Magnet, 90deg analyzing go to 95 120 close(1) close(2) stop C C integrate area of phase space C 125 ara25 = 0. ara50 = 0. ara100 = 0. ara200 = 0. ara400 = 0. do 130 j = 1,max if(abs(tf(j)).le.p25) ara25 = ara25 + 1 if(abs(tf(j)).le.p50) ara50 = ara50 + 1 if(abs(tf(j)).le.p100) ara100 = ara100 + 1 if(abs(tf(j)).le.p200) ara200 = ara200 + 1 if(abs(tf(j)).le.p400) ara400 = ara400 + 1 130 continue pc25 = 100. * ara25 / max pc50 = 100. * ara50 / max pc100 = 100. * ara100 / max pc200 = 100. * ara200 / max pc400 = 100. * ara400 / max write(iout,135) tarea,pc25,pc50,pc100,pc200,pc400 135 format(' AREA & % ',6f10.3) go to 95 C ********************************************************************** C C ACCELERATION REGION C 140 continue write(iout,145) 145 format(' ENTER ACCELERATING VOLTAGE AND LENGTH : ',$) read(in,*) prev,s write(iout,150) 150 format(' ENTER VOLTAGE STABILITY dV [keV] : ',$) read(in,*) stab write(2,155) prev,s 155 format(/,' ACCELERATING VOLTAGE: ',f10.3, 1 'kV LENGTH= ',f10.3,'m') write(2,160) stab 160 format(' VOLTAGE STABILITY= ',f10.3,'keV') accel = .09 * q * prev / (s * mass) u0 = u0 + q * prev ac4 = 2. * accel * s c ustab = stab * q * prev c term = (stab * q * prev) ustab = stab * q term = (stab * q) term = term * term time = (-vzero + sqrt(vzero * vzero + ac4)) / accel do 165 j = 1,max t = (-vf(j) + sqrt(vf(j) * vf(j) + ac4)) / accel tf(j) = tf(j) + time - t theta = 6.283 * ran(iseed) ef(j) = ef(j) + ustab * sin(theta) vf(j) = c * sqrt(2 * (u0 + ef(j)) / mass) 165 continue tzero = tzero + time vzero = c * sqrt(2 * u0 / mass) l = l + s write(2,170) m,q,u0,beta0,du 170 format(' ION M;Q;U0;BETA;DU: ',5f10.3) call project(0,tmax,emax) ! project results onto e and t axes go to 60 c ********************************************************************* c c drift region c 175 continue write(iout,180) 180 format(' ENTER DRIFT LENGTH: ',$) read(in,*) drl write(2,185) drl 185 format(' DRIFT LENGTH ',f10.3,' m') vzero = c * sqrt(2 * u0 / mass) time = drl / vzero do 190 j = 1,max tf(j) = tf(j) + time - drl / vf(j) 190 continue tzero = tzero + time l = l + drl call project(0,tmax,emax) ! project results onto e and t axes go to 60 c ********************************************************************** c c print results c 195 continue write(*,200) 200 format(///,' BUNCHING PARAMETERS') write(*,205) m,q,u0,du 205 format(' ION M,Q,U0,DU: ',4f10.3) write(*,210) l,tzero 210 format(' TOTAL DISTANCE & TRANSIT TIME: ',2f10.3) write(*,215) v,p,twaist 215 format(' GRID VOLTAGE, PERIOD,TWAIST: ',3f10.3) write(*,220) prev,s,stab 220 format(' ACCELERATION VOLTAGE, LENGTH, STABILITY: ', 1 3f10.3,//) write(*,225) 225 format(1x,' CHANNEL TIME ENERGY ') write(*,230) 230 format(1x,' NUMBER PROJECTION PROJECTION') do 240 j = 1,40 write(*,235) j,tbins(j),ebins(j) 235 format(1x,6x,i2,10x,f4.0,11x,f4.0) 240 continue write(*,245) tsize,esize 245 format(/,' TIME AND ENERGY CHANNEL CALIBRATON:', 1 f10.3,'ns/ch',f10.3,'kev/ch') write(*,250) tfwhm,efwhm 250 format(/,' TIME FWHM',f8.4,'ns ENERGY FWHM',f8.4,' keV',//) 255 go to 60 c ********************************************************************** c c buncher c 260 continue write(iout,265) 265 format(' buncher type 1=>sinewave 2=>sawtooth 3=>linear ',$) write(iout,270) 270 format(' ENTER BUNCHER TYPE,VOLTAGE (kV),', 1 'PERIOD (ns), STABILITY (dV/V)') read(in,*) ityp,v,p,bstab write(iout,275) 275 format(' ENTER BUNCHER PHASE RELATIVE TO FIRST: ',$) read(in,*) phi omega = 6.283 / p phi = .017453 * phi write(iout,280) 280 format(' DISTANCE TO WAIST OR POINT OF INTEREST: (m.) ',$) read(in,*) bl write(iout,285) 285 format(' ENTER TIME TO WAIST: (ns.) ',$) read(in,*) twaist if(ityp.eq.1) write(2,290) 290 format(' BUNCHER TYPE IS SINEWAVE') if(ityp.eq.2) write(2,295) 295 format(' BUNCHER TYPE IS SAWTOOTH') if(ityp.eq.3) write(2,300) 300 format(' BUNCHER TYPE IS LINEAR') write(2,305) v,p,phi 305 format(' BUNCHER VOLTAGE:',f10.3,'kV PERIOD:',f10.3, 1 'ns. RELATIVE PHASE:',f10.3) write(2,310) bstab 310 format(' BUNCHER VOLTAGE STABILITY:',f10.3) write(2,315) bl,twaist 315 format(' DRIFT DISTANCE:',f10.3,' TIME TO WAIST:',f10.3) vzero = c * sqrt(2 * u0 / mass) l = l + bl tzero = bl / vzero aterm = 2 * u0 * u0 * u0 / mass if(bl.gt.0.) vexp = c * p / (2 * q * bl) * sqrt(aterm) if(ityp.eq.1) go to 345 if(ityp.eq.2) go to 355 if(ityp.eq.3) go to 335 write(iout,320) write(*,320) 320 format(' INVALID BUNCHER TYPE') go to 95 325 do 330 j = 1,max u = u0 + ef(j) vf(j) = c * sqrt(2 * u / mass) tf(j) = tf(j) + tzero - bl / vf(j) 330 continue call project(0,tmax,emax) ! project results onto e and t axes go to 60 c ********************************************************************* c c linear bunching c 335 continue if(twaist.eq.0.) twaist = tzero hafmass = mass / 2 qv = u0 * (p/twaist + 3.*p*p/(4.*twaist*twaist)) if(v.eq.0.) then v = qv / q u1 = u0 else u1 = u0 * v * q / qv endif ustab = u1 * bstab do 340 j = 1,max termf = tf(j) / twaist theta = 6.283 * ran(iseed) u2 = u1 + ustab * sin(theta) ef(j) = ef(j)+ 1 u2*(-2*termf + 3*termf*termf - 4*termf*termf*termf) 340 continue go to 325 c ********************************************************************** c c sinewave bunching c 345 continue if(v.eq.0.) v = .7136 * vexp qv = q * v ustab = qv * bstab do 350 j = 1,max theta = 6.283 * ran(iseed) u2 = qv + ustab * sin(theta) ef(j) = ef(j) + u2 * sin(3.14159 + phi + omega * tf(j)) 350 continue go to 325 c ********************************************************************** c c sawtooth bunching c 355 if(v.eq.0.) v = 1.998 * vexp term = 2 * q * v / p ustab = term * bstab do 360 j = 1,max theta = 6.283 * ran(iseed) ef(j) = ef(j) - tf(j) * (term + ustab * sin(theta)) 360 continue go to 325 C*********************************************************************** C C PLOT HISTOGRAMS OF ENERGY AND TIME C 365 continue gflag=2. do 370 j=1,50 xpos(j)=float(j) 370 continue call pgbegin(0,device,1,1) c encode(7,375,cc) tcenter encode(7,375,cf) tfwhm clabel='Centroid '//cc//' ns' call pgtext(0.2,0.35,clabel) clabel=' FWHM '//cf//' ns' call pgtext(0.2,0.3,clabel) c encode(7,375,cc) ecenter encode(7,375,cf) efwhm clabel='Centroid '//cc//' keV' call pgtext(0.68,0.35,clabel) clabel=' FWHM '//cf//' keV' call pgtext(0.68,0.3,clabel) 375 format(f7.3) c call pgvport(0.20,0.45,0.60,0.90) call pgwindow(-tmax,tmax,0.,1.2*ytmax) call pgbox('abnpst',0.,0,'abcnpstv',0.,0) call pglabel('T ns',' ',' ') call pghist(1000,tf,-tmax,tmax,40,1) c call pgvport(0.60,0.85,0.60,0.90) call pgwindow(-emax,emax,0.,1.2*yemax) call pgbox('abnpst',0.,0,'abcnpstv',0.,0) call pglabel('E keV',' ',' ') call pghist(1000,ef,-emax,emax,40,1) c c call pgupdt read(*,380) ans 380 format(a1) call pgadvance call pgend go to 60 C ********************************************************************** C C C GENERATE RANDOMLY FILLED PHASE SPACE C 385 continue tmax = bprd / 2 emax = du / 2 a = du / 2.355 sign = 1 scale = 1.414 * a sign = 1 if(distr.eq.1.) go to 395 do 390 j = 1,max tf(j) = -tmax + 2 * tmax * ran(iseed) ef(j) = -emax + 2 * emax * ran(iseed) vf(j) = c * sqrt(2 * (u0 + ef(j)) / mass) 390 continue go to 60 395 do 400 j = 1,max term = 6.2832 * ran(iseed) term1 = -2. * log(ran(iseed)) asine = sin(term) x = asine * term1 tf(j) = tmax * (2 * ran(iseed) - 1) ef(j) = emax * asine * term1 / 2.526 vf(j) = c * sqrt(2 * (u0 + ef(j)) / mass) 400 continue call project(0,tmax,emax) go to 60 C ********************************************************************** C C CHANGE GRAPH SCALE C 405 continue write(iout,410) 410 format(' ENTER ENERGY SCALE AND TIME SCALE : ',$) read(in,*) emax, tmax call project(1,tmax,emax) if(gflag.eq.1) call graph(tf,tmax,ef,emax,max,device) go to 60 C ********************************************** C ERROR HANDLING ROUTINE C 415 close(1) in=5 iout=6 manual = 1 go to 60 C **************************************************** C C PRINT COMMENT LINE TO FILE C 420 read(in,425) comment 425 format(a40) write(2,*) comment go to 60 430 continue 435 go to 60 C c ********************************************************************* c c analyzing magnet drift region c c models transit time difference effects c due to the size and divergence of the beam c entering the magnet and time dispersion c from energy spread not included. c 2/2/89 c 440 continue write(iout,445) 445 format(' ENTER RADIUS OF CURVATURE [m]: ',$) read(in,*) curv drl=3.14159*2.*curv/4. ! flight path for 90 deg bend write(iout,446) 446 format(' ENTER BEAM DIVERGENCE [mrad], BEAM RADIUS [cm]: ',$) read(in,*) div,radius write(2,450) curv,drl 450 format(' RAD CURVATURE ',f6.3, * 'm DRIFT LENGTH (90deg bend) ',f10.3,'m') write(2,452) div,radius 452 format(' BEAM DIVERGENCE, RADIUS [mrad,cm]: ',2f10.3) div=div/1000. radius=radius/100. sigd=curv*div sigr=curv*asin(radius/curv) dstab=sqrt(sigd**2.+sigr**2.) vzero = c * sqrt(2 * u0 / mass) time = drl / vzero do 455 j = 1,max theta=6.283 * ran(iseed) dist=drl + dstab * sin(theta) + drl * ef(j) /2. /u0 tf(j) = tf(j) + time - dist / vf(j) 455 continue tzero = tzero + time l = l + drl call project(0,tmax,emax) ! project results onto e and t axes go to 60 C end C C ********************************************************************** subroutine graph(tf,tmax,ef,emax,max,device) c dimension tf(1000),ef(1000) character*3 device gflag=1 c call pgbegin(0,device,1,1) call pgvport(0.25,0.75,0.2,0.8) call pgwindow(-tmax,tmax,-emax,emax) call pgbox('anpst',0.,0,'anpstv',0.,0) call pgpoint(max,tf,ef,1) call pglabel('t ns','Energy keV','PHASE PLOT') call pgupdt read(*,5) ans 5 format(a1) call pgask(.false.) call pgadvance call pgend return end C ********************************************************************** subroutine project(iq,tmax,emax) real tbins(50),ebins(50),tf(1000),ef(1000) common/proj/tbins,ebins,ytmax,yemax,tf,ef,max, 1 tcenter,ecenter,tfwhm,efwhm c c project particles onto time and energy axis c if(iq.eq.1) go to 10 tmax = 0. emax = 0. do 5 j = 1,max if(abs(tf(j)).gt.tmax) tmax = abs(tf(j)) if(abs(ef(j)).gt.emax) emax = abs(ef(j)) 5 continue 10 continue ! start here if scale changed tsize = tmax / 19. esize = emax / 19. do 15 j = 1,40 tbins(j) = 0. ebins(j) = 0. 15 continue sumslo = 0. n = 1 slorng = .2 * tmax do 25 j = 1,max if(abs(tf(j)).gt.slorng) go to 20 sumslo = sumslo + tf(j) / ef(j) n = n + 1 20 binno = tf(j) / tsize + 20.5 if(binno.gt.0.5) then iloc=int(binno) if(binno.lt.40) tbins(iloc)=tbins(iloc)+1 endif binno = ef(j) / esize + 20.5 if(binno.gt.0.5) then iloc=int(binno) if(binno.lt.40) ebins(iloc)=ebins(iloc)+1 endif 25 continue tsum = 0. tsumchn = 0. tsumsq = 0. esum = 0. esumchn = 0. esumsq = 0. sumslo = sumslo / float(n - 1) ytmax = 0. yemax = 0. do 30 j = 1,40 tsum = tsum + tbins(j) tsumchn = tsumchn + j * tbins(j) esum = esum + ebins(j) esumchn = esumchn + j * ebins(j) tsumsq = tsumsq + j * j * tbins(j) esumsq = esumsq + j * j * ebins(j) if(tbins(j).gt.ytmax) ytmax = tbins(j) if(ebins(j).gt.yemax) yemax = ebins(j) 30 continue term = tsumchn / tsum tcenter = tsize * (term - 20.) tsigma = tsize * sqrt(tsumsq / tsum - term * term) tfwhm = 2.355 * tsigma term = esumchn / esum ecenter = esize * (term - 20.) esigma = esize * sqrt(esumsq / esum - term * term) efwhm = 2.355 * esigma write(2,35) 35 format(' PHASE SPACE PARAMETERS') write(2,40) ecenter,efwhm 40 format(' DELTA ENERGY CENTROID=',f10.3,'keV FWHM=',f10.3,'keV') write(2,45) emax 45 format(' ONE-HALF MAXIMUM ENERGY RANGE=',f10.3,'keV') write(2,50) tcenter,tfwhm 50 format(' DELTA TIME CENTROID=',f10.3,'ns. FWHM=',f10.3,'ns.') write(2,55) tmax 55 format(' ONE-HALF MAXIMUM TIME RANGE=',f10.3,'ns.') write(2,60) sumslo 60 format(' BEAM SLOPE =',f10.3) return end C C C subroutine getdate(cdate) character*9 cdate call date(cdate) return end C C C subroutine gettime(ctime) character*8 ctime call time(ctime) return end