c PARTICLE RAY TRACING CODE " ERAY " c c Written by Z.Q. Xie c 1 st Version, June 6, 1987 c c Note: Accuracy of this Code is the user's own risk. c c ERAY Code does the 3D particle ray tracing of beam of multiple ion c species extracted from the RT-ECR by using the iterative numerical solution c (Runge-Kutta method). The ray number limit of this version is 800, each ray c carries a portion of the total beam current partitioned by the user at their c desires. c c In this version, all the particles initially are assumed to follow the c magnetic field line and emitt from a plane surface with no radial velocity. c c The magnetic filed (including the source solenoid and the first focussing c solenoid coil magnetic fields) and the extraction electric field (including c the extraciotn and the puller) components are linearly interpolated based on c the input data files (which are generated with the help of POISSON Code, 4 c data files should be generated: 1. source solenoid magnetic field (fixed c field, usually a typical field); 2. end electrode electric field data file c in POISSON calculation, 10000 V is applied to the end electrode and 0 V is c applied to the puller); 3. puller electrode electric field data file (in c POISSON calculation -100 V is applied to the puller and 0 V is applied to the c end electrode); 4. focussing solenoid magnetic data file (in POISSON c calculation, 100 A should be applied to the calulation). Thus inside ERAY c routine, factors of 1/10000, -1/100 and 1/100 are scaled to the input data c accordingly to generate the real fields. The magnetic contribution from the c hexapole magnet can be analytically calculated from the "PERMMAG" Code which c is ingored in this first stage of development for the time being in order to c save CPU time. PERMAG will be incorporated in ERAY as a subrouting in the c future version. c c The source solenoid field input data file is formated in specific step size c of DZ=1 mm in the Z and DR=1 mm in the r direction when r<1 cm, and DR=5 mm c when r>1 cm. The electric field data file have the same Z step size but DR= c 1.27 mm in the region (Z=-.3 cm to EXEFL and R to 2.159 cm). The focussing c solenoid field data file is formated DZ = DR = 5mm. c The rest input data are: initial position in X and Y axises with Z=-3 mm is c assumed to be the starting point. Z=0 is at one of the extraction aperture c edges in the puller electrode side; c (In this version, the fields are tracted back to Z=-3 mm, the thikness c of the extraction aperture is assumed about 3 mm thick, in -3 mm to z=0. c region, considering the total current is contributed from the particles c under consideration which get through the extraction aperture, then in the c above mentioned region, the radial electric field is taken as zero.) c The extraction electric potentials (V); The input initial velocities in all c 3 dimensions in the unit of eV; Charge state Q and mass number A. c c The electric field inside the conductor should be zero but the POISSON c calculation does not give the correct values due to the deficiency of the c POISSON Code. When the particle get very close to the electrode within a c grid unit, instead of putting pseudo values into those ponits, the values c from POISSON just inside the conductor are used in the 2-D linearly c interpoltation of the field strength outside the conductor but within a c grid unit, and those particle ray will be notified in the output LOG data c file. c c The space charge are assumed to uniform distributed inside any circle of c radius r and has only the r component. Rays are allowed to cross each other. c The charge inside the circle is related to the total beam currnet by a c definition (total current)*(eclosed area)/(total area) which depends on z and c r. Details of space charge handling are refered to "PRELIMINARY REPORT ON A c 3D EXTRACTION CODE FOR ECR SOURCES", Int. Con. on ECR Ion Sources, NSCL c REPORT #MSUCP-47, DEC 1987, P420. c c In the output data file, all the notations are self-explained, all the c lengths are in mm unit, the magnetic field strengths are in Tesla, velocities c are in 10E6 m/s unit. After the ray tracing calculation, a routine name c UP1000Q can be called to read the ERAY output data file to generate a plotting c data file NN.PLT and a routine ERAYPLOT can be used to generate the output c plot graph on a computer (VAX) line printer. See 'explanation' for sample c input and more details. c c CHARACTER*13 optf,bscoil,eend,epull,bfcoil real*8 EX,EY,EZ,X(800),Y(800),X0(800),Y0(800),VX0(800),RM(800) real*8 VX1(800),VY1(800),VY0(800),VZ0(800),VZ1(800),B,BX,BY,BZ real*8 COMMON ERR,EK,CCK,BHX,BHY,BHZ,X01(800),Y01(800),QS(800) real*8 RCC(20),BCZ1(20),BCZ3(20),BCR1(20),BCR3(20),Q(800),A(800) REAL*8 BBCX1,BBCY1,BBCZ1,BBCX2,BBCY2,BBCZ2,BBCR1,BBCR2 REAL*8 BBCX,BBCY,BCZ,BCR,Z,DZ,Z0,VEXR1(20),VEXZ1(20),VEXR3(20) REAL*8 VEXZ3(20),EXVR1(20),EXVR3(20),EXVZ1(20),EXVZ3(20),RMX,RC REAL*8 PULVR1(20),PULVR3(20),PULVZ1(20),PULVZ3(20),AP(800) INTEGER NSTOP(800),NC(800) REAL ZSTOP(800),ZC(800),ZRC(800),TT(800),TL(800) REAL*8 BSR1(20),BSR3(20),BSZ1(20),BSZ3(20),E_TRANS(800) REAL*8 E_LONG(800) real*8 YI(800,6),CK(6),QM(800,6),ARK(4),BRK(4),CRK(4) INTEGER*4 IR c Runge-Kutta coefficients DATA ARK/.5D+0,.292893219D+0,1.707106781D+0,.1666666667D+0/, & BRK/2.D+0,1.D+0,1.D+0,2.D+0/, & CRK/-.5D+0,-.292893219D+0,-1.707106781D+0,-.5D+0/ CTAA Type *, ' Output data file name?' CTAA read(5,FMT='(A13)')optf Type *, ' Source magnetic data file name?' c read(5,FMT='(A13)')bscoil Type *, ' End electrode data file name?' read(5,FMT='(A13)')eend Type *, ' Puller electrode data file name?' read(5,FMT='(A13)')epull Type *, ' Focussing solenoid data file name?' read(5,FMT='(A13)')bfcoil OPEN(10,STATUS='NEW') OPEN(11,FILE=bscoil,STATUS='OLD') OPEN(12,FILE=eend,STATUS='OLD') OPEN(13,FILE=epull,STATUS='OLD') OPEN(14,FILE=bfcoil,STATUS='OLD') cc IR is the random generator for randomly generate the initial particle cc possitions, if IR=1, then the initial particle possitions are key board cc inputted. TYPE *, ' TYPE IR ' READ(5,*) IR TYPE *, ' HOW MANY RAYS, TYPE NP (NPmax=800)' READ (5,*) NP TYPE *, ' CHARGE STATES IN THE ORDER OF 1 TO NP' READ(5,*)(Q(I),I=1,NP) TYPE *, ' CURRENT ( UA ) PER RAY ' READ(5,*)(AP(I),I=1,NP) TYPE *,' PARTICLE MASS IN THE ORDER OF 1 TO NP' READ(5,*)(A(I),I=1,NP) IF(IR.EQ.1) GOTO 65 type *, ' INITIAL ENERGIES IN TRANSVERSE DIRECTION (eV) ?' READ(5,*) (E_trans(I),I=1,NP) TYPE *,' INITIAL ENERGIES IN LONGITUDINAL DIRECTION (eV) ?' read(5,*) (E_long(I),I=1,NP) GOTO 76 65 TYPE *, ' INITIAL POSITIONS IN THE ORDER X0 AND Y0 ( MM 1) OF 1 TO NP' READ(5,*)(X01(I),Y01(I),I=1,NP) TYPE *,' INITIAL TRANS. AND LONGI. ENERGIES IN eV IN THE ORDER OF 1 1 TO NP' READ(5,*)(TT(I),TL(I),I=1,NP) 76 Z0=-3. TYPE *,' Z STEP SIZE ( MM ) AND NUMBER OF POINTS WHICH WILL' TYPE *,' BE WRITTEN INTO THE OUTPUT DATA FILE. THE POINT ' TYPE *,' STEP SIZE IS NZ*(Z STEP SIZE ).' TYPE *,' TYPE DZ , NZ , M ' READ(5,*) DZ0,NZ,M TYPE *,' EXTRACTION APERTURE (RADIUS), PULLER APERTURE' TYPE *,' (RADIUS), EXTRACTION SPACING, RADIUS' TYPE *,' OF BEAM PIPE, EXTRACTION ELETRIC FIELD LENGTH,' TYPE *,' EXTRACTION MAGNETIC FIELD LENGTH AND FOCUSSING' TYPE *,' SOLENOID FIELD LENGTH, ALL IN METER UNIT (m).' READ(5,*) EXR,PUR,EXSP,RPIP,EXEFL,EXBFL,FSFL EXRS=EXR*1000. WRITE(10,10) 10 FORMAT(5X,' K ',10X,' A ',10X,' Q ',5X,'I (UA)/RAY',/) DO K=1,NP WRITE(10,11)K,A(K),Q(K),AP(K) 11 FORMAT(4X,I4,10X,F4.0,9X,F4.0,5X,F8.4) ENDDO AP1=0. DO K=1,NP AP1=AP1+AP(K)*.001 ENDDO c c neutralzation factor is used to reduce the space force strength. c the second neutralization factor (FNEU_NF) will oversee the first c neutralization FNEU1 starting from ZN1 till ZN2. c TYPE *, ' NEUTRALIZATION FACTOR ( < 1 ), IF NEUTRALIZATION ' TYPE *, ' OCCURS AT SOME REGION, TYPE 1 (ELSE 0), REGION SIZE' TYPE *, ' (ZN1 AND ZN2 (M)), NEUTRALIZATION FACTOR IN THAT REGION' READ(5,*)FNEU1,NF,ZN1,ZN2,FNEU_NF WRITE(10,13)DZ0,NZ,M,AP1,FNEU1 13 FORMAT(//,' DZ = ',F8.5,2X,' POINT STEP SIZE = ',I4, 12X,' MUNBER OF POINTS = ',I4,/,5X,' BEAM CURRENT = ', 1 F5.2,' mA ',5X,' FNEU1 =',F4.2,/) type *,' VEXTRACTION ( V ) AND PULLER VOLTAGE ( in V )' type *, ' Percent of the source B field intensity?' READ(5,*) V_ext,V_pul,bpercent c C bpercent is the scaling factor of the source magnetic field. c Which assumes all the 9 coils increase or decrease the same percentage c of its original input currents. c TYPE *, ' SOLENOID CURRENTS (A) AND LOCATION (M) ?' TYPE *, ' OBJECT SLIT WIDTH AND LOCATION (M) ?' READ(5,*) SOL_CURI,SOL_LOC,OSLIT,OSLT_LOC write(10,15)V_ext,V_pul,SOL_CURI,SOL_LOC 15 format(3x,'V_EXTR = ',F8.0,' V',3X,'V_PUL = ',F8.0,' V', 1' SOL_CURI = ',F8.3,' AMP',3X,'SOL_LOC = ',F5.3,' M',/) WRITE(10,16) 16 FORMAT(4X,'K',6X,'Z',10X,'X',10X,'Y',10X,'Vx',12x,'Vy',12X,'Vz', 19X,'Et ( eV )',5x,'El ( eV )',5X,'Bz0 ( T )',//,16x,' length 1dimensions are in mm and velocities are in 10E6 M/S ',//) IF(IR.EQ.1) GOTO 170 L=1 DO 169 I=1,10000 X01(L)=EXRS*(2.*RAN(IR)-1.) Y01(L)=EXRS*(2.*RAN(IR)-1.) R00=(X01(L)*X01(L)+Y01(L)*Y01(L))**.5 IF(R00.GT.EXRS*.999) GOTO 169 TT(L)=RAN(IR)*Q(L)*E_trans(L)/2. C ACCORDING TO THE POISSON CALCULATION, THE ELECTRIC POTENTIAL AT THE C STARTING POINT Z=-3 mm IS ABOUT .002 X (APPLIED VOLTAGE), SO AN ENERGY OF C AMOUNT Q*V_ext*.001 eV IS ADDED TO THE PARTICLE LONGITUDINAL ENERGY IN C ORDER TO MACTH THE REAL CASE CLOSER THAN ASSUMING THE PARTICLE AT THAT C POINT HAS ONLY THE THERMAL ENERGY. TL(L)=RAN(IR)*Q(L)*E_long(L)/2.+Q(L)*V_ext*.004 IF(L.EQ.NP) GOTO 170 L=L+1 169 CONTINUE 170 IF(IR.EQ.1) THEN DO I=1,NP TL(I)=TL(I)+Q(I)*V_ext*.004 ENDDO ENDIF do k=1,NP X0(k)=X01(k)*.001 Y0(k)=Y01(k)*.001 RM(K)=(X0(K)*X0(K)+Y0(K)*Y0(K))**.5 VT=(2.*96484570.*TT(K)/A(K))**.5 IF(RM(K).EQ.0.) THEN VX0(k)=0. VY0(k)=0. GOTO 17 ENDIF AN0=RAN(IR)*2*3.1415927 VX0(k)=VT*COS(AN0) VY0(k)=VT*SIN(AN0) 17 VZ0(k)=(2.*96484570.*TL(K)/A(K))**.5 VX1(K)=VX0(K)/1000000. VY1(K)=VY0(K)/1000000. VZ1(K)=VZ0(K)/1000000. enddo CC FNEU=1.-FNEU1 C ELECTRIC FIELD DATE FILE DUE TO EXTRACTION END PLATE IS CALCULATED WITH C 10 kV APPLIED VOLTAGE, ELECTRIC FIELD DATE FILE DUE TO PULLER IS C CALCULATED WITH 100 V APPLIED VOLTAGE, THE FOCUSSING SOLENOID MAGNETIC C DATE FIELD IS PRECALCULATED WITH 100 A. THE FOLLOWING EXF, EPUL AND BFSOL C ARE THE SCALE FIELD FACTORS WITH THE CORRESPONDING INPUTS. C Exf=V_ext/10000. Epul=-V_pul/100. BFSOL=SOL_CURI/100. C READ THE FIRST TWO INPUT DATA READ(11,*)ZC1,(BCR1(I),I=1,11) READ(11,*)ZC1,(BCR1(I),I=12,19) READ(11,*)ZC1,(BCZ1(I),I=1,11) READ(11,*)ZC1,(BCZ1(I),I=12,19) READ(11,*)ZC3,(BCR3(I),I=1,11) READ(11,*)ZC3,(BCR3(I),I=12,19) READ(11,*)ZC3,(BCZ3(I),I=1,11) READ(11,*)ZC3,(BCZ3(I),I=12,19) READ(12,*)ZZC1,(EXVR1(I),I=1,9) READ(12,*)ZZC1,(EXVR1(I),I=10,18) READ(12,*)ZZC1,(EXVZ1(I),I=1,9) READ(12,*)ZZC1,(EXVZ1(I),I=10,18) READ(13,*)ZZC1,(PULVR1(I),I=1,9) READ(13,*)ZZC1,(PULVR1(I),I=10,18) READ(13,*)ZZC1,(PULVZ1(I),I=1,9) READ(13,*)ZZC1,(PULVZ1(I),I=10,18) READ(12,*)ZZC3,(EXVR3(I),I=1,9) READ(12,*)ZZC3,(EXVR3(I),I=10,18) READ(12,*)ZZC3,(EXVZ3(I),I=1,9) READ(12,*)ZZC3,(EXVZ3(I),I=10,18) READ(13,*)ZZC3,(PULVR3(I),I=1,9) READ(13,*)ZZC3,(PULVR3(I),I=10,18) READ(13,*)ZZC3,(PULVZ3(I),I=1,9) READ(13,*)ZZC3,(PULVZ3(I),I=10,18) READ(14,*)ZSS1,(BSR1(I),I=1,8) READ(14,*)ZSS1,(BSR1(I),I=9,16) READ(14,*)ZSS1,(BSZ1(I),I=1,8) READ(14,*)ZSS1,(BSZ1(I),I=9,16) READ(14,*)ZSS3,(BSR3(I),I=1,8) READ(14,*)ZSS3,(BSR3(I),I=9,16) READ(14,*)ZSS3,(BSZ3(I),I=1,8) READ(14,*)ZSS3,(BSZ3(I),I=9,16) ZS1=SOL_LOC+ZSS1 ZS3=SOL_LOC+ZSS3 DO I=1,18 VEXR1(I)=EXF*EXVR1(I)+EPUL*PULVR1(I) VEXZ1(I)=EXF*EXVZ1(I)+EPUL*PULVZ1(I) VEXR3(I)=EXF*EXVR3(I)+EPUL*PULVR3(I) VEXZ3(I)=EXF*EXVZ3(I)+EPUL*PULVZ3(I) ENDDO Bz0=BCZ1(1)*bpercent do K=1,NP write(10,200)K,Z0,X01(K),Y01(K),VX1(K),VY1(K),VZ1(K),TT(K), 1TL(K),Bz0 enddo CCK=1.79836*(100000.**2.) DO K=1,NP YI(k,4)=X0(K) YI(K,5)=Y0(K) YI(K,1)=VX0(K) YI(K,2)=VY0(K) YI(K,3)=VZ0(K) ENDDO DZ=DZ0*.001 Z0=-0.003-DZ NT=0 MC=0 AH=1. AE=1. AGE=23.806*3.1415927/180. ZSOL=ZS1 do k=1,NP do i=1,5 QM(K,I)=0. enddo enddo NZZ=NZ*DZ0 5 DO 800 l=1,M IF(Z.GT.EXEFL) THEN NZ=NZZ DZ=.001 ENDIF DO 100 N=1,NZ Z=DZ*N+Z0 IF((NF.EQ.1).AND.((Z.GE.ZN1).AND.(Z.LE.ZN2))) THEN FNEU=1.-FNEU_NF ELSE FNEU=1.-FNEU1 ENDIF IF(Z.GT.SOL_LOC+FSFL) GOTO 301 IF(Z.GE.ZSOL) THEN IF(Z.GT.ZS3) THEN ZS1=ZS3 DO I=1,16 BSR1(I)=BSR3(I) BSZ1(I)=BSZ3(I) ENDDO READ(14,*)ZSS3,(BSR3(I),I=1,8) READ(14,*)ZSS3,(BSR3(I),I=9,16) READ(14,*)ZSS3,(BSZ3(I),I=1,8) READ(14,*)ZSS3,(BSZ3(I),I=9,16) ZS3=ZSS3+SOL_LOC ENDIF ENDIF IF(Z.GT.EXBFL) GOTO 301 IF(Z.GT.ZC3) THEN ZC1=ZC3 DO I=1,19 BCR1(I)=BCR3(I) BCZ1(I)=BCZ3(I) ENDDO DO I=1,18 VEXR1(I)=VEXR3(I) VEXZ1(I)=VEXZ3(I) ENDDO READ(11,*)ZC3,(BCR3(I),I=1,11) READ(11,*)ZC3,(BCR3(I),I=12,19) READ(11,*)ZC3,(BCZ3(I),I=1,11) READ(11,*)ZC3,(BCZ3(I),I=12,19) IF(Z.GT.EXEFL) GOTO 301 READ(12,*)ZZC3,(EXVR3(I),I=1,9) READ(12,*)ZZC3,(EXVR3(I),I=10,18) READ(12,*)ZZC3,(EXVZ3(I),I=1,9) READ(12,*)ZZC3,(EXVZ3(I),I=10,18) READ(13,*)ZZC3,(PULVR3(I),I=1,9) READ(13,*)ZZC3,(PULVR3(I),I=10,18) READ(13,*)ZZC3,(PULVZ3(I),I=1,9) READ(13,*)ZZC3,(PULVZ3(I),I=10,18) DO I=1,18 VEXR3(I)=EXF*EXVR3(I)+EPUL*PULVR3(I) VEXZ3(I)=EXF*EXVZ3(I)+EPUL*PULVZ3(I) ENDDO ENDIF CC C TAKE THE CENTER MASS VELOCITY OF THE M PARTICLES AS THE BEAM VELOCITY CC 301 DO 309 K=1,NP DO NS=1,NT IF(K.EQ.NSTOP(NS)) GOTO 309 ENDDO RC=RM(K) if((Z.LT..0).and.(RC.GT.EXR)) then nt=nt+1 nstop(nt)=k ZSTOP(NT)=Z QS(NT)=Q(K) goto 309 endif CC RAYS WILL BE STOPPED IF ITS RADIUS LARGER THAN THE PULLER RADIUS AT THE CC VICINITY OF THE PULLER SINCE THESE RAYS WILL HIT THE PULLER. LACKING OF CC THE EXACT SHAPE OF THE PULLER, THE 45 DEGREE TILT SHAPE OF THE PULLER OF CC THE MSU RT-ECR IS USED HERE AS THE CRITERION TO STOP THOSE RAYS CAN NOT GO CC THROUGH THE PULLER. USERS CAN MODIFY THIS PART TO MEET THEIR DESIRES. CC if((Z.ge.EXSP).and.(Z.le.EXSP+.01524)) AE=RC-PUR-(Z-EXSP) if(ABS(AE).le..0001) then nt=nt+1 nstop(nt)=k ZSTOP(NT)=Z QS(NT)=Q(K) AE=1. goto 309 endif if(RC.GT.RPIP) then nt=nt+1 nstop(nt)=k ZSTOP(NT)=Z QS(NT)=Q(K) goto 309 endif if((abs(Z-OSLT_LOC).le..001).and.(RC.gt.OSLIT)) then nt=nt+1 nstop(nt)=k ZSTOP(NT)=Z QS(NT)=Q(K) goto 309 endif IF(Z.GT.EXEFL) THEN ER=0. EZ=0. GOTO 35 ENDIF DO I=1,18 RCC(I)=.00127*((I-1)*1.) IF(RC.LT.RCC(I)) THEN CR1=RCC(I-1) CR2=RCC(I) if ((z.le.0.).and.(rc.gt..00381)) VEXZ3(I)=VEXZ3(I-1) EZ=(RC-CR1)*((VEXZ3(I)-VEXZ3(I-1)+VEXZ1(I-1)-VEXZ1(I)) 1 *(Z-ZC1)+(VEXZ1(I)-VEXZ1(I-1))*(ZC3-ZC1))/(ZC3-ZC1) 1 /(CR2-CR1)+(VEXZ3(I-1)-VEXZ1(I-1))*(Z-ZC1)/(ZC3-ZC1) 1 +VEXZ1(I-1) IF((RC.EQ.0.).or.(Z.le.0.)) THEN ER=0. GOTO 25 ENDIF ER=(Z-ZC1)*((VEXR3(I)-VEXR3(I-1)+VEXR1(I-1)-VEXR1(I)) 1 *(RC-CR1)+(VEXR3(I-1)-VEXR1(I-1))*(CR2-CR1))/(ZC3-ZC1) 1 /(CR2-CR1)+(VEXR1(I)-VEXR1(I-1))*(RC-CR1)/(CR2-CR1) 1 +VEXR1(I-1) 25 DO NS=1,NP IF(K.EQ.NC(NS)) GOTO 35 ENDDO CC KEEP TRACKING THE RAYS CLOSE TO THE END ELECTRODE SURFACE WITHIN ONE GRID CC UNIT. IF(Z.GT..0) THEN IF((((CR2-EXR-ZC1/TAN(AGE)).GT..0000001).AND. 1(ZC1.LT..009525)).OR.((CR2-PUR-(ZC1-EXSP)).GT..0000001) 1 .AND.((ZC1.LT..06223).AND.(ZC1.GE..03651))) THEN MC=MC+1 NC(MC)=K ZC(MC)=Z ZRC(MC)=RC ENDIF GOTO 35 ELSE IF(ABS(RC-EXR).LE..00127) THEN MC=MC+1 NC(MC)=K ZC(MC)=Z ZRC(MC)=RC ENDIF GOTO 35 ENDIF ENDIF ENDDO 35 if((Z.GT.EXBFL).or.(bpercent.eq.0.)) THEN BBCX1=0. BBCY1=0. BBCZ1=0. goto 39 ENDIF DO I=1,19 if(I.LE.11) then RCC(I)=.001*(I-1) else RCC(I)=.01+.005*(I-11) ENDIF IF(RC.LT.RCC(I)) THEN CR1=RCC(I-1) CR2=RCC(I) BBCZ1=((RC-CR1)*((BCZ3(I)-BCZ3(I-1)+BCZ1(I-1)-BCZ1(I)) 1 *(Z-ZC1)+(BCZ1(I)-BCZ1(I-1))*(ZC3-ZC1))/(ZC3-ZC1) 1 /(CR2-CR1)+(BCZ3(I-1)-BCZ1(I-1))*(Z-ZC1)/(ZC3- 1 ZC1)+BCZ1(I-1))*bpercent IF(RC.EQ.0.) THEN BBCX1=0. BBCY1=0. GOTO 39 ENDIF BBCR1=((Z-ZC1)*((BCR3(I)-BCR3(I-1)+BCR1(I-1)-BCR1(I)) 1 *(RC-CR1)+(BCR3(I-1)-BCR1(I-1))*(CR2-CR1))/(ZC3-ZC1) 1 /(CR2-CR1)+(BCR1(I)-BCR1(I-1))*(RC-CR1)/(CR2-CR1) 1 +BCR1(I-1))*bpercent BBCX1=-BBCR1*YI(k,4)/RC BBCY1=-BBCR1*YI(K,5)/RC GOTO 39 ENDIF ENDDO 39 IF((Z.GE.ZSOL).AND.(Z.LE.SOL_LOC+FSFL)) THEN DO I=1,16 RCC(I)=.005*(I-1) IF(RC.LT.RCC(I)) THEN CR1=RCC(I-1) CR2=RCC(I) BBCZ2=((RC-CR1)*((BSZ3(I)-BSZ3(I-1)+BSZ1(I-1)-BSZ1(I)) 1 *(Z-ZS1)+(BSZ1(I)-BSZ1(I-1))*(ZS3-ZS1))/(CR2-CR1) 1 /(ZS3-ZS1)+(BSZ3(I-1)-BSZ1(I-1))*(Z-ZS1)/(ZS3- 1 ZS1)+BSZ1(I-1))*bfsol IF(RC.EQ.0.) THEN BBCX2=0. BBCY2=0. GOTO 53 ENDIF BBCR2=((Z-ZS1)*((BSR3(I)-BSR3(I-1)+BSR1(I-1)-BSR1(I)) 1 *(RC-CR1)+(BSR3(I-1)-BSR1(I-1))*(CR2-CR1))/(ZS3-ZS1) 1 /(CR2-CR1)+(BSR1(I)-BSR1(I-1))*(RC-CR1)/(CR2-CR1) 1 +BSR1(I-1))*bfsol BBCX2=BBCR2*YI(k,4)/RC BBCY2=BBCR2*YI(K,5)/RC GOTO 53 ENDIF ENDDO ENDIF 53 IF(RC.EQ.0.) THEN EX=0. EY=0. GOTO 54 ENDIF TLD=0. DO 99 J=1,NP DO NQ=1,NT IF(J.EQ.NSTOP(NQ)) GOTO 99 ENDDO IF(RM(J).GT.RC) GOTO 99 TLD=TLD+AP(J)*.000001/YI(J,3) 99 CONTINUE ERR=FNEU*CCK*TLD/RC EX=(ER+ERR)*YI(k,4)/RC EY=(ER+ERR)*YI(K,5)/RC IF(Z.LT.0.) EX=0. IF(Z.LT.0.) EY=0. 54 BX=BBCX1+BBCX2 BY=BBCY1+BBCY2 BZ=BBCZ1+BBCZ2 IF(Z.GT.SOL_LOC+FSFL) THEN BX=0. BY=0. BZ=0. ENDIF DO 450 J=1,4 C----------------------------------------------------------------------- C EQUATIONS OF MOTION C----------------------------------------------------------------------- CK(1)=Q(K)/A(K)/YI(K,3)*(EX+YI(K,2)*BZ-YI(K,3)*BY)*96484570.*DZ CK(2)=Q(K)/A(K)/YI(K,3)*(EY+YI(K,3)*BX-YI(K,1)*BZ)*96484570.*DZ CK(3)=Q(K)/A(K)/YI(K,3)*(EZ+YI(K,1)*BY-YI(K,2)*BX)*96484570.*DZ CK(4)=YI(K,1)/YI(K,3)*DZ CK(5)=YI(K,2)/YI(K,3)*DZ DO 350 I=1,5 RRK=ARK(J)*(CK(I)-BRK(J)*QM(K,I)) YI(K,I)=YI(K,I)+RRK 350 QM(K,I)=QM(K,I)+3.*RRK+CRK(J)*CK(I) 450 CONTINUE 309 CONTINUE DO 104 I=1,NP RM(I)=(YI(I,4)*YI(I,4)+YI(I,5)*YI(I,5))**.5 104 CONTINUE 100 CONTINUE Z0=Z BZ0=BCZ1(1)*bpercent IF((Z.GE.ZSOL).AND.(Z.LE.EXBFL)) BZ0=BSZ1(1)*bfsol+ & BCZ1(1)*bpercent IF((Z.GT.EXBFL).AND.(Z.LE.SOL_LOC+FSFL)) BZ0=BSZ1(1)*bfsol IF(((Z.GT.EXBFL).AND.(Z.LT.ZSOL)).or.(Z.GT.SOL_LOC+FSFL)) BZ0=0. DO K=1,NP XF=1000.*YI(k,4) YF=1000.*YI(K,5) ZF=1000.*(Z+DZ) DO NS=1,NT IF(K.EQ.NSTOP(NS)) ZF=ZSTOP(NS)*1000. ENDDO VXF=YI(K,1)*.000001 VYF=YI(K,2)*.000001 VZF=YI(K,3)*.000001 TTE=1/2.*A(K)/96484570.*(YI(K,1)*YI(K,1)+YI(K,2)*YI(K,2)) TLE=1/2.*A(K)/96484570.*YI(K,3)*YI(K,3) WRITE(10,200)K,ZF,XF,YF,VXF,VYF,VZF,TTE,TLE,Bz0 200 FORMAT(I5,3(1X,F10.4),5(1X,E13.6),2x,F6.4) ENDDO 800 CONTINUE 999 DO I=1,NP IF(NSTOP(I).EQ.0) GOTO 1030 IF(ZSTOP(I).GT..1) THEN WRITE(6,1020)NSTOP(I) ELSE WRITE(6,1010)NSTOP(I) ENDIF 1010 FORMAT(/,2X,' PARTICLE',I4,' HIT THE PULLER.') 1020 FORMAT(/,2X,' PARTICLE',I4,' HIT THE BEAM PIPE.') ENDDO 1030 DO I=1,NP IF(NC(I).EQ.0) GOTO 1110 ZCL=ZC(I)*1000. ZRCL=ZRC(I)*1000. WRITE(6,1040) NC(I),ZCL,ZRCL 1040 FORMAT(/,2X,' PARTICLE',I4,' WAS CLOSE TO THE ELECTRODE WITHIN A 1 GRID UNIT AT Z =',F7.2,' AND R =',F7.2,' mm') ENDDO 1110 STOP END