C C C C BUNCH A nonrelativistic longitudinal optics program C R.C. Pardo 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. 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 of 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 C BVOLTAGE, PERIOD - Maximum buncher C voltage in kV, period of buncher in C nanoseconds. C C RELPHS - phase,in radians, of this C buncher relative to center of central C particle. Used only for TYPE =1. C 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 C TWAIST - time to waist, in nanoseconds, C used only for type = 3, C 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. No parameters. C Display a two dimensional phase space C plot of the beam ensemble at this point C along the beamline. 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 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. TMAX C - The maximum delta time, in nsec., to C 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 C Original Version : R. Pardo, ANL, written in BASIC, C received at UK Dec 1988. C C Modified : J. Vanhoy, Univ. of KY, changed to FORTRAN, C Jan '89 C C C***************************************************************************** 10 dimension tbins(50),ebins(50) 12 dimension tf(1000),ef(1000),vf(1000) real xpos(50),mass,m,l,lzero,hafmas 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, * tcenter,ecenter,tfwhm,efwhm C iseed=011159 14 p25 = .125 p50 = .25 p100 = .5 16 p200 = 1. p400 = 2. 18 manual = 0 io=1 write(*,19) 19 format(' enter input file name: ',$) read(*,20) fspec 20 format(a20) if(fspec.eq.' ') then io=5 go to 25 endif open(unit=1,file=fspec,type='old') 25 open(unit=2,file='out.dat',type='new') 28 write(io+1,29) 29 format(' ENTER 1 FOR GAUSSIAN OR 0 FOR UNIFORM DIST ',$) read(io,*) distr 30 continue ! cls max = 1000 c c initialize system c 1515 l = 0. tzero = 0. v = 0. p = 0. twaist = 0. 1520 prev = 0. s = 0. stab = 0. bl = 0. gflag = 1. c write(io+1,41) 41 format(' ENTER BEAM UZERO(keV);MASS;Q;DU(keV) ',$) read(io,*) u0,m,q,du write(io+1,50) 50 format(' ENTER BEAM BUNCH PERIOD (ns.) ',$) read(io,*) bprd write(2,57) fspec 57 format(' INPUT FILE NAME: ',a20) call getdate(cdate) call gettime(ctime) write(2,*) cdate,' ',ctime 60 phi = 0. csq = .09 65 c = .3 mass = 931507. * m 70 vzero = c * sqrt(2 * u0 / mass) beta0 = vzero / c write(2,58) m,q,u0,beta0,du 58 format(' ION M;Q;U0;BETA;DU: ',5F10.4) write(2,59) bprd 59 format(' BEAM BUNCH PERIOD: ',F10.3,'ns.') 75 go to 1400 ! RPHASE 85 continue write(io+1,*) ' bunching parameters ' 86 beta0 = vzero / c tarea = du * bprd ! pbrd write(io+1,88) ecenter,efwhm 88 format(' DELTA ENERGY CENTROID= ',f7.3, * ' keV FWHM= ',f7.3,' keV') write(io+1,90) tcenter,tfwhm 90 format(' DELTA TIME CENTROID= ',f7.3, * ' keV FWHM= ',f7.3,' keV') write(io+1,99) m,q,u0,beta0,du 99 format(' ION M;Q;U0(keV);BETA;DU(keV): ',5f10.3) write(io+1,101) v,p,twaist,bstab 101 format(' GRID VOLTAGE, PERIOD, TWAIST, STABILITY ',4f10.3) write(io+1,102) prev,s,stab 102 format(' ACCEL. VOLTAGE, LENGTH, STABILITY: ',3f10.3) write(io+1,107) l,lzero 107 format(' TOTAL DISTANCE,TRANSIT TIME: ',2f10.3) c c command section c 125 continue write(2,127) l,tzero 127 format(' TOTAL DISTANCE: ',f10.3,' TRANSIT TIME: ',f10.3) 130 continue write(io+1,131) 131 format(' COMMAND (a,b,d,g,h,l,q,p,r,s,i) : ',$) read(io,132) chk 132 format(a1) write(2,133) chk 133 format(' COMMAND: ',a1) 137 if(chk.eq.'c') go to 2100 !C Comment line in input file Echo to file 138 if(chk.eq.'e') go to 2040 !E End of input file data 139 if(chk.eq.'g') then call graph(tf,tmax,ef,emax,max) go to 85 !G Graph results endif 140 if(chk.eq.'l') go to 2200 !L PRINT LOG OF RESULT SO FAR 142 if(chk.eq.'p') go to 695 !P Print 145 if(chk.eq.'r') then close(1) close(2) go to 18 !R Restart endif 150 if(chk.eq.'s') go to 1730 !S Change graph scale 155 if(chk.eq.'q') go to 185 !Q Quit 160 if(chk.eq.'a') go to 470 !A Accelerator tube 165 if(chk.eq.'b') go to 775 !B Buncher 170 if(chk.eq.'d') go to 605 !D Drift region 175 if(chk.eq.'h') go to 1200 !H Histogram energy and time projections 176 if(chk.eq.'i') go to 230 !I Integrate phase space 177 if(chk.eq.'n') go to 28 !N New case, same file 180 go to 125 185 close(1) close(2) stop C C integrate area of phase space C 230 ara25 = 0. ara50 = 0. ara100 = 0. ara200 = 0. ara400 = 0. 231 do 265 j = 1,max 240 if(abs(tf(j)).le.p25) ara25 = ara25 + 1 245 if(abs(tf(j)).le.p50) ara50 = ara50 + 1 250 if(abs(tf(j)).le.p100) ara100 = ara100 + 1 255 if(abs(tf(j)).le.p200) ara200 = ara200 + 1 260 if(abs(tf(j)).le.p400) ara400 = ara400 + 1 265 continue 270 pc25 = 100. * ara25 / max pc50 = 100. * ara50 / max pc100 = 100. * ara100 / max 275 pc200 = 100. * ara200 / max pc400 = 100. * ara400 / max write(io+1,277) tarea,pc25,pc50,pc100,pc200,pc400 277 format(' AREA & % ',6f10.3) 285 go to 125 C ************************************************************************ C C GRAPH RESULTS C C330 SCALEY = -(1.3 * EMAX) / 75 C SCALEX = (1.3 * TMAX) / 150 C335 CALL ScreenTest(EM%, CR%, VL%, VR%, VT%, VB%): CLS : GFLAG = 1 C336 WINDOW SCREEN (0, 0)-(620, 200) C337 'VIEW (VL, VT)-(VR, VB) C340 LINE (470, 100)-(100, 100): LINE (320, 25)-(320, 175) C345 LINE (282, 97)-(282, 103): LINE (245, 97)-(245, 103): LINE (207, 97)-(207, 103) C350 LINE (316, 81)-(324, 81): LINE (316, 62)-(324, 62): LINE (316, 44)-(324, 44) C355 LINE (316, 119)-(324, 119): LINE (316, 138)-(324, 138): LINE (316, 156)-(324, 156) C360 LINE (358, 97)-(358, 103): LINE (395, 97)-(395, 103): LINE (433, 97)-(433, 103) C395 FOR INDX = 1 TO MAX C400 X% = TF(INDX) / SCALEX + 320 C405 Y% = EF(INDX) / SCALEY + 100 C410 IF X% < 1 OR X% > 639 OR Y% < 1 OR Y% > 199 THEN GOTO 420 C415 PSET (X%, Y%) C420 NEXT INDX C425 LOCATE 14, 61: PRINT CHR$(127); : PRINT "T ns"; C430 LOCATE 3, 42: PRINT CHR$(127); : PRINT "E kev"; C435 LOCATE 14, 30: PRINT ; -TMAX / 2; : LOCATE 14, 49: PRINT ; TMAX / 2; C440 LOCATE 8, 42: PRINT ; EMAX / 2; : LOCATE 18, 42: PRINT ; -EMAX / 2; C445 IF ITYP = 3 THEN LOCATE 2, 34: PRINT "LINEAR BUNCHING"; C450 IF ITYP = 2 THEN LOCATE 2, 33: PRINT "SAWTOOTH BUNCHING"; C455 IF ITYP = 1 THEN LOCATE 2, 33: PRINT "SINEWAVE BUNCHING"; C DO C LOOP WHILE INKEY$ = "" C SCREEN 0 C460 RETURN C ##################################################### C ************************************************************************ C C ACCELERATION REGION C 470 continue write(io+1,475) 475 format(' ENTER ACCELERATING VOLTAGE AND LENGTH ',$) 490 read(io,*) prev,s write(io+1,495) 495 format(' ENTER VOLTAGE STABILITY : ',$) 500 read(io,*) stab write(2,503) prev,s 503 format(/,' ACCELERATING VOLTAGE: ',f10.3, * 'kV LENGTH= ',f10.3,'m') write(2,506) stab 506 format(' VOLTAGE STABILITY= ',f10.3,' kV') 510 accel = .09 * q * prev / (s * mass) 515 u0 = u0 + q * prev 520 ac4 = 2. * accel * s ustab = stab * q * prev 525 term = (stab * q * prev) term = term * term 530 time = (-vzero + sqrt(vzero * vzero + ac4)) / accel 535 do 570 j = 1,max 540 t = (-vf(j) + sqrt(vf(j) * vf(j) + ac4)) / accel 545 tf(j) = tf(j) + time - t 550 theta = 6.283 * ran(iseed) 555 ef(j) = ef(j) + ustab * sin(theta) 560 vf(j) = c * sqrt(2 * (u0 + ef(j)) / mass) 570 continue 575 tzero = tzero + time 580 vzero = c * sqrt(2 * u0 / mass) 582 l = l + s write(2,584) m,q,u0,beta0,du 584 format(' ION M;Q;U0;BETA;DU: ',5f10.3) 586 call project(0,tmax,emax) ! project results onto e and t axes c 590 call graph (COMMENTED OUT) 595 go to 85 c ********************************************************************* c c drift region c 605 continue write(io+1,620) 620 format(' ENTER DRIFT LENGTH: ',$) read(io,*) drl write(2,625) drl 625 format(' DRIFT LENGTH ',f10.3,' m') 630 vzero = c * sqrt(2 * u0 / mass) 635 time = drl / vzero 640 do 655 j = 1,max 645 tf(j) = tf(j) + time - drl / vf(j) 655 continue 660 tzero = tzero + time 665 l = l + drl 666 call project(0,tmax,emax) ! project results onto e and t axes c 670 call graph 675 go to 85 c *********************************************************************** c c print results c 695 continue write(*,700) 700 format(///,' BUNCHING PARAMETERS') write(*,705) m,q,u0,du 705 format(' ION M,Q,U0,DU: ',4f10.3) write(*,710) l,tzero 710 format(' TOTAL DISTANCE & TRANSIT TIME: ',2f10.3) write(*,715) v,p,twaist 715 format(' GRID VOLTAGE, PERIOD,TWAIST: ',3f10.3) write(*,720) prev,s,stab 720 format(' ACCELERATION VOLTAGE, LENGTH, STABILITY: ', * 3f10.3,//) c 730'lprint " ef tf c ef tf ef tf" write(*,731) 731 format(1x,' CHANNEL TIME ENERGY ') write(*,732) 732 format(1x,' NUMBER PROJECTION PROJECTION') 733 do 745 j = 1,40 write(*,734) j,tbins(j),ebins(j) 734 format(1x,6x,i2,10x,f4.0,11x,f4.0) c 735 'for j=1 to max step 3 c 740 'lprint using "####.#### ";ef(j);tf(j);ef(j+1);tf(j+1);ef(j+2);tf(j+2) 745 continue write(*,746) tsize,esize 746 format(/,' TIME AND ENERGY CHANNEL CALIBRATON:', * f10.3,'ns/ch',f10.3,'kev/ch') write(*,747) tfwhm,efwhm 747 format(/,' TIME FWHM',f8.4,'ns ENERGY FWHM',f8.4,' keV',//) c 754 lprint chr$(12) 755 go to 85 c *********************************************************************** c c buncher c 775 continue write(io+1,781) 781 format(' buncher type 1=>sinewave 2=>sawtooth 3=>linear ',$) write(io+1,784) 784 format(' ENTER BUNCHER TYPE,VOLTAGE (kV),', * 'PERIOD (ns), STABILITY (dV/V)') read(io,*) ityp,v,p,bstab write(io+1,795) 795 format(' ENTER BUNCHER PHASE RELATIVE TO FIRST: ',$) read(io,*) phi 800 omega = 6.283 / p phi = .017453 * phi write(io+1,810) 810 format(' DISTANCE TO WAIST OR POINT OF INTEREST: (m.) ',$) read(io,*) bl write(io+1,820) 820 format(' ENTER TIME TO WAIST: (ns.) ',$) read(io,*) twaist if(ityp.eq.1) write(2,826) 826 format(' BUNCHER TYPE IS SINEWAVE') if(ityp.eq.2) write(2,827) 827 format(' BUNCHER TYPE IS SAWTOOTH') if(ityp.eq.3) write(2,828) 828 format(' BUNCHER TYPE IS LINEAR') write(2,829) v,p,phi 829 format(' BUNCHER VOLTAGE:',f10.3,'kV PERIOD:',f10.3, * 'ns. RELATIVE PHASE:',f10.3) write(2,830) bstab 830 format(' BUNCHER VOLTAGE STABILITY:',f10.3) write(2,831) bl,twaist 831 format(' DRIFT DISTANCE:',f10.3,' TIME TO WAIST:',f10.3) 832 vzero = c * sqrt(2 * u0 / mass) l = l + bl 835 tzero = bl / vzero 840 aterm = 2 * u0 * u0 * u0 / mass 845 if(bl.gt.0.) vexp = c * p / (2 * q * bl) * sqrt(aterm) 850 if(ityp.eq.1) go to 990 855 if(ityp.eq.2) go to 1045 860 if(ityp.eq.3) go to 940 write(io+1,861) write(*,861) 861 format(' INVALID BUNCHER TYPE') go to 125 865 do 890 j = 1,max 870 u = u0 + ef(j) 875 vf(j) = c * sqrt(2 * u / mass) 880 tf(j) = tf(j) + tzero - bl / vf(j) 890 continue 891 call project(0,tmax,emax) ! project results onto e and t axes c 892 call graph 893 go to 85 c *********************************************************************** c c linear bunching c 940 continue 945 if(twaist.eq.0.) twaist = tzero 950 hafmass = mass / 2 qv = u0 * (p/twaist + 3.*p*p/(4.*twaist*twaist)) 955 if(v.eq.0.) then v = qv / q u1 = u0 else u1 = u0 * v * q / qv endif 956 ustab = u1 * bstab do 975 j = 1,max 965 termf = tf(j) / twaist 966 theta = 6.283 * ran(iseed) 967 u2 = u1 + ustab * sin(theta) 970 ef(j) = ef(j)+ * u2*(-2*termf + 3*termf*termf - 4*termf*termf*termf) 975 continue 980 go to 865 c *********************************************************************** c c sinewave bunching c 990 continue 1000 if(v.eq.0.) v = .7136 * vexp 1005 qv = q * v 1006 ustab = qv * bstab 1010 do 1020 j = 1,max 1011 theta = 6.283 * ran(iseed) 1012 u2 = qv + ustab * sin(theta) 1015 ef(j) = ef(j) + u2 * sin(3.14159 + phi + omega * tf(j)) 1020 continue 1025 go to 865 c *********************************************************************** c c sawtooth bunching c 1045 if(v.eq.0.) v = 1.998 * vexp 1050 term = 2 * q * v / p 1051 ustab = term * bstab 1055 do 1065 j = 1,max 1056 theta = 6.283 * ran(iseed) 1060 ef(j) = ef(j) - tf(j) * (term + ustab * sin(theta)) 1065 continue 1070 go to 865 C *********************************************************************** C C PROJECT PARTICLES ONTO TIME AND ENERGY AXIS C C 1090 TMAX = 0 C EMAX = 0 C 1105 DO 1115 J = 1,MAX C 1110 IF(ABS(TF(J)).GT.TMAX) TMAX = ABS(TF(J)) C 1111 IF(ABS(EF(J)).GT.EMAX) EMAX = ABS(EF(J)) C 1115 CONTINUE C entry point inside project @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C 1120 TSIZE = TMAX / 19 C ESIZE = EMAX / 19 C 1121 DO 1123 J = 1,40 C 1122 TBINS(J) = 0 C EBINS(J) = 0 C 1123 CONTINUE C 1124 SUMSLO = 0 C N = 1 C SLORNG = .2 * TMAX C 1125 DO 1140 J = 1,MAX C 1126 IF(ABS(TF(J)).GT.SLORNG) GO TO 1130 C 1127 SUMSLO = SUMSLO + TF(J) / EF(J) C 1128 N = N + 1 C 1130 BINNO = TF(J) / TSIZE + 20.5 C 1135 IF(BINNO.GT.0.5) THEN C IF(BINNO.LT.40) TBINS(BINNO)=TBINS(BINNO)+1 C ENDIF C 1136 BINNO = EF(J) / ESIZE + 20.5 C 1137 IF(BINNO.GT.0.5) THEN C IF(BINNO.LT.40) EBINS(BINNO)=EBINS(BINNO)+1 C 1140 CONTINUE C 1145 TSUM = 0 C TSUMCHN = 0 C TSUMSQ = 0 C 1146 ESUM = 0 C ESUMCHN = 0 C ESUMSQ = 0 C 1150 SUMSLO = SUMSLO / (N - 1) C YTMAX = 0 C YEMAX = 0 C 1155 DO 1175 J = 1,40 C 1160 TSUM = TSUM + TBINS(J) C TSUMCHN = TSUMCHN + J * TBINS(J) C 1161 ESUM = ESUM + EBINS(J) C ESUMCHN = ESUMCHN + J * EBINS(J) C 1165 TSUMSQ = TSUMSQ + J * J * TBINS(J) C 1166 ESUMSQ = ESUMSQ + J * J * EBINS(J) C 1170 IF(TBINS(J).GT.YTMAX) YTMAX = TBINS(J) C 1171 IF(EBINS(J).GT.YEMAX) YEMAX = EBINS(J) C 1175 CONTINUE C 1176 TERM = TSUMCHN / TSUM C 1180 TCENTER = TSIZE * (TERM - 20.) C 1181 TSIGMA = TSIZE * SQR(TSUMSQ / TSUM - TERM * TERM) C 1182 TFWHM = 2.355 * TSIGMA C 1183 TERM = ESUMCHN / ESUM C 1184 ECENTER = ESIZE * (TERM - 20.) C 1185 ESIGMA = ESIZE * SQR(ESUMSQ / ESUM - TERM * TERM) C 1186 EFWHM = 2.355 * ESIGMA C WRITE(2,1190) C 1190 FORMAT(' PHASE SPACE PARAMETERS') C WRITE(2,1191) ECENTER,EFWHM C 1191 FORMAT(' DELTA ENERGY CENTROID=',F10.3,'keV FWHM=',F10.3,'keV') C WRITE(2,1192) EMAX C 1192 FORMAT(' ONE-HALF MAXIMUM ENERGY RANGE=',F10.3,'keV') C WRITE(2,1193) TCENTER,TFWHM C 1193 FORMAT(' DELTA TIME CENTROID=',F10.3,'ns. FWHM=',F10.3,'ns.') C WRITE(2,1194) TMAX C 1194 FORMAT(' ONE-HALF MAXIMUM TIME RANGE=',F10.3,'ns.') C WRITE(2,1195) SUMSLO C 1195 FORMAT(' BEAM SLOPE =',F10.3) C 1199 RETURN C END C************************************************************************* C C PLOT HISTOGRAMS OF ENERGY AND TIME C 1200 continue gflag=2. do 1201 j=1,50 xpos(j)=float(j) 1201 continue call pgbegin(0,'/re',1,1) c encode(7,1202,cc) tcenter encode(7,1202,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,1202,cc) ecenter encode(7,1202,cf) efwhm clabel='Centroid '//cc//' keV' call pgtext(0.68,0.35,clabel) clabel=' FWHM '//cf//' keV' call pgtext(0.68,0.3,clabel) 1202 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',' ',' ') c call pgline(40,xpos,tbins) 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',' ',' ') c call pgline(40,xpos,ebins) call pghist(1000,ef,-emax,emax,40,1) c c call pgupdt read(*,1210) ans 1210 format(a1) call pgadvance call pgend go to 85 C1200 TSCALEY = -YTMAX / 75: SCALEX = 5: ESCALEY = -YEMAX / 75 C1202 CALL ScreenTest(EM%, CR%, VL%, VR%, VT%, VB%): CLS : GFLAG = 2 C1203 WINDOW SCREEN (0, 0)-(620, 200) 'VIEW (VL, VT)-(VR, VB) C1205 LINE (100, 100)-(300, 100): LINE (100, 25)-(100, 100) C1210 LINE (175, 97)-(175, 100): LINE (150, 97)-(150, 100): LINE (125, 97)-(125, 100) C1215 LINE (100, 81)-(105, 81): LINE (100, 62)-(105, 62): LINE (100, 44)-(105, 44) C1220 LINE (200, 97)-(200, 103) C1225 LINE (225, 97)-(225, 100): LINE (250, 97)-(250, 100): LINE (275, 97)-(275, 100) C1230 XL0 = 100: YL0 = TBINS(1) / TSCALEY + 100 C1235 FOR J = 2 TO 40 C1240 XL = XL0 + SCALEX C1245 YL = TBINS(J) / TSCALEY + 100 C1250 LINE (XL0, YL0)-(XL, YL) C1255 XL0 = XL: YL0 = YL C1260 NEXT J C1265 LOCATE 15, 23: PRINT CHR$(127); : PRINT "T ns"; C1270 'SYMBOL (50,100),"EVENTS",1,1,7,3 C1275 LOCATE 8, 9: PRINT ; YTMAX / 2; C1280 LOCATE 14, 15: PRINT ; -10 * TSIZE; : LOCATE 14, 28: PRINT ; 10 * TSIZE; C1285 LOCATE 16, 12: PRINT USING "CENTROID ###.##"; TCENTER; C1290 PRINT USING " FWHM ###.#####"; TFWHM; C1315 LINE (420, 100)-(620, 100): LINE (420, 25)-(420, 100) C1320 LINE (495, 97)-(495, 100): LINE (470, 97)-(470, 100): LINE (445, 97)-(445, 100) C1325 LINE (420, 81)-(425, 81): LINE (420, 62)-(425, 62): LINE (420, 44)-(425, 44) C1330 LINE (520, 97)-(520, 103) C1335 LINE (545, 97)-(545, 100): LINE (570, 97)-(570, 100): LINE (595, 97)-(595, 100) C1340 XL0 = 420: YL0 = EBINS(1) / ESCALEY + 100 C1345 FOR J = 2 TO 40 C1350 XL = XL0 + SCALEX C1355 YL = EBINS(J) / ESCALEY + 100 C1360 LINE (XL0, YL0)-(XL, YL) C1365 XL0 = XL: YL0 = YL C1370 NEXT J C1371 IF ITYP = 3 THEN LOCATE 2, 34: PRINT "LINEAR BUNCHING"; C1372 IF ITYP = 2 THEN LOCATE 2, 33: PRINT "SAWTOOTH BUNCHING"; C1373 IF ITYP = 1 THEN LOCATE 2, 33: PRINT "SINEWAVE BUNCHING"; C1375 LOCATE 15, 62: PRINT CHR$(127); : PRINT "E kev"; C1380 'SYMBOL (375,100),"EVENTS",1,1,7,3 C1390 LOCATE 8, 49: PRINT ; YEMAX / 2; C1395 LOCATE 14, 56: PRINT ; -10 * ESIZE; : LOCATE 14, 68: PRINT ; 10 * ESIZE; C1396 LOCATE 16, 50: PRINT USING "CENTROID ###.##"; ECENTER; C1397 PRINT USING " FWHM ###.####"; EFWHM; C DO C LOOP WHILE INKEY$ = "" C SCREEN 0 C1398 GO TO 85 C *********************************************************************** C C C GENERATE RANDOMLY FILLED PHASE SPACE C 1400 continue 1410 tmax = bprd / 2 emax = du / 2 1411 a = du / 2.355 sign = 1 1412 scale = 1.414 * a 1415 sign = 1 1416 if(distr.eq.1.) go to 1450 1420 do 1440 j = 1,max 1425 tf(j) = -tmax + 2 * tmax * ran(iseed) 1430 ef(j) = -emax + 2 * emax * ran(iseed) 1435 vf(j) = c * sqrt(2 * (u0 + ef(j)) / mass) 1440 continue 1445 go to 85 1450 do 1480 j = 1,max 1451 term = 6.2832 * ran(iseed) term1 = -2. * log(ran(iseed)) 1452 asine = sin(term) 1453 x = asine * term1 1455 tf(j) = tmax * (2 * ran(iseed) - 1) 1460 ef(j) = emax * asine * term1 / 2.526 c 1461 'sign=-sign c 1465 'tf(j)=sign*tmax*y/2! c 1470 'ef(j)=sign*scale*sqrt(-term) 1475 vf(j) = c * sqrt(2 * (u0 + ef(j)) / mass) 1480 continue 1481 call project(0,tmax,emax) 1485 go to 85 C *********************************************************************** C C INITIALIZE SYSTEM C C 1515 L = 0 C TZERO = 0 C V = 0 C P = 0 C TWAIST = 0 C 1520 PREV = 0 C S = 0 C STAB = 0 C BL = 0 C GFLAG = 1 C 1525 RETURN C END C '************************************************** C ' C 'ARBITRARY ORIENTATION LINE C ' C'1600 DEL = 1: Z% = 1 C' IF ABS(XL = XL0) < ABS(YL - YL0) GOTO 1645 C' SLOPE = (YL - YL0) / (XL - XL0) C' IF XL < XL0 THEN DEL = -1 C' FOR X% = XL0 TO XL STEP DEL C' Y% = (X% - XL0) * SLOPE + YL0 C' CALL POINTH(X%, Y%, Z%, O%): NEXT X% C' RETURN C'1645 SLOPE = (XL - XL0) / (YL - YL0) C' IF YL < YL0 THEN DEL = -1 C' FOR Y% = YL0 TO YL STEP DEL C' X% = (Y% - YL0) * SLOPE + XL0 C' CALL POINTH(X%, Y%, Z%, O%): NEXT Y% C' RETURN C *********************************************************************** C C CHANGE GRAPH SCALE C 1730 continue write(io+1,1740) 1740 format(' ENTER ENERGY SCALE AND TIME SCALE : ',$) 1742 read(io,*) emax, tmax 1755 call project(1,tmax,emax) 1760 if(gflag.eq.1) call graph 1780 go to 85 C 2000 '********************************************** C 2010 'ERROR HANDLING ROUTINE C 2020 ' C LOCATE 17, 1: PRINT "ERROR #"; ERR; " LINE POSITION= "; ERL 2040 close(1) C 2050 OPEN "KYBD:" FOR INPUT ACCESS READ AS #1 io=5 2055 manual = 1 2060 go to 85 C **************************************************** C C PRINT COMMENT LINE TO FILE C 2100 read(io,2101) comment 2101 format(a40) 2105 write(2,*) comment 2110 go to 85 C 2190 '******************************************** C 2191 ' LIST LOG OF RESULTS TO PRINTER C 2192 ' 2200 continue C 2200 CLOSE #2: OPEN "OUT.DAT" FOR INPUT ACCESS READ AS #2 C 2210 WHILE NOT EOF(2) C 2220 INPUT #2, COMMNT$: LPRINT USING "&"; COMMNT$ C 2230 WEND C 2240 CLOSE #2: OPEN "OUT.DAT" FOR APPEND ACCESS READ WRITE AS 2 2250 go to 85 end C'************************************************** C'handle screen error interrupt C' CBadscreen: C EM = 0 C SCREEN 3 C VL = 10: VR = 600 C VT = 5: VB = 300 C CR = 3 ' Hercules graphics assumed C RESUME NEXT CDEFINT A-Z C'************************************************** C'Handle bad file specification C' CBadfile: C LOCATE 24, 1 C PRINT "Bad File Specification" C CLOSE #1 C CLOSE #2 C END C' ************************************************* C' initialize for keyboard input C' CKeyboard: C OPEN "KYBD:" FOR INPUT ACCESS READ AS #1 C manual = 1 C ON ERROR GOTO Badfile C RESUME 22 C C C' C' ======================== ScreenTest ======================== C' Tests to see if user has EGA hardware with SCREEN 8. C' If this causes an error, the EM flag is set to FALSE, C' and the screen is set with SCREEN 3. C' C' Also sets values for corners of viewport (VL = left, C' VR = right, VT = top, VB = bottom), scaled with the C' correct aspect ratio so viewport is a perfect square. C' ============================================================ C' CSUB ScreenTest (EM, CR, VL, VR, VT, VB) STATIC C EM = 1 'Assume that screen is EGA C VL = 110: VR = 529 C VT = 5: VB = 179 C CR = 15 ' 16 colors (0 - 15) C ON ERROR GOTO Badscreen C SCREEN 8 C ON ERROR GOTO 0 CEND SUB C C ************************************************************************ subroutine graph(tf,tmax,ef,emax,max) c dimension tf(1000),ef(1000) gflag=1 c call pgbegin(0,'/re',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(*,100) ans 100 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, * tcenter,ecenter,tfwhm,efwhm c c project particles onto time and energy axis c if(iq.eq.1) go to 1117 1090 tmax = 0. emax = 0. 1105 do 1115 j = 1,max 1110 if(abs(tf(j)).gt.tmax) tmax = abs(tf(j)) 1111 if(abs(ef(j)).gt.emax) emax = abs(ef(j)) 1115 continue 1117 continue ! start here if scale changed 1120 tsize = tmax / 19. esize = emax / 19. 1121 do 1123 j = 1,40 1122 tbins(j) = 0. ebins(j) = 0. 1123 continue 1124 sumslo = 0. n = 1 slorng = .2 * tmax 1125 do 1140 j = 1,max 1126 if(abs(tf(j)).gt.slorng) go to 1130 1127 sumslo = sumslo + tf(j) / ef(j) 1128 n = n + 1 1130 binno = tf(j) / tsize + 20.5 1135 if(binno.gt.0.5) then iloc=int(binno) if(binno.lt.40) tbins(iloc)=tbins(iloc)+1 endif 1136 binno = ef(j) / esize + 20.5 1137 if(binno.gt.0.5) then iloc=int(binno) if(binno.lt.40) ebins(iloc)=ebins(iloc)+1 endif 1140 continue 1145 tsum = 0. tsumchn = 0. tsumsq = 0. 1146 esum = 0. esumchn = 0. esumsq = 0. 1150 sumslo = sumslo / float(n - 1) ytmax = 0. yemax = 0. 1155 do 1175 j = 1,40 1160 tsum = tsum + tbins(j) tsumchn = tsumchn + j * tbins(j) 1161 esum = esum + ebins(j) esumchn = esumchn + j * ebins(j) 1165 tsumsq = tsumsq + j * j * tbins(j) 1166 esumsq = esumsq + j * j * ebins(j) 1170 if(tbins(j).gt.ytmax) ytmax = tbins(j) 1171 if(ebins(j).gt.yemax) yemax = ebins(j) 1175 continue 1176 term = tsumchn / tsum 1180 tcenter = tsize * (term - 20.) 1181 tsigma = tsize * sqrt(tsumsq / tsum - term * term) 1182 tfwhm = 2.355 * tsigma 1183 term = esumchn / esum 1184 ecenter = esize * (term - 20.) 1185 esigma = esize * sqrt(esumsq / esum - term * term) 1186 efwhm = 2.355 * esigma write(2,1190) 1190 format(' PHASE SPACE PARAMETERS') write(2,1191) ecenter,efwhm 1191 format(' DELTA ENERGY CENTROID=',f10.3,'keV FWHM=',f10.3,'keV') write(2,1192) emax 1192 format(' ONE-HALF MAXIMUM ENERGY RANGE=',f10.3,'keV') write(2,1193) tcenter,tfwhm 1193 format(' DELTA TIME CENTROID=',f10.3,'ns. FWHM=',f10.3,'ns.') write(2,1194) tmax 1194 format(' ONE-HALF MAXIMUM TIME RANGE=',f10.3,'ns.') write(2,1195) sumslo 1195 format(' BEAM SLOPE =',f10.3) 1199 return end subroutine getdate(cdate) character*9 cdate call date(cdate) return end subroutine gettime(ctime) character*8 ctime call time(ctime) return end