C THIS IS THE KENT STATE NEUTRON DETECTOR EFFICIENCY CODE C REFERENCE: R.CECIL, B.ANDERSON AND R.MADEY, NUCL. INSTRUM. C AND METH. 161, 439 (1979). C CONTACT PERSON: B.ANDERSON, PHYSICS DEPT., KENT STATE C UNIV., KENT, OHIO 44242 C EMAIL: ANDERSON@KSUVXD.KENT.EDU C C--------------------------------------------------------------------- C C MONTE CARLO NEUTRON COUNTER EFFICIENCY PROGRAM C MAIN ROUTINE C IN THIS ROUTINE WE READ IN ALL INPUT DATA NECESSARY TO C CALCULATE A NEUTRON EFFICIENCY FOR A PARTICULAR C SCINTILLATOR GEOMETRY C C WE THEN CALL FOLNUT TO PROPAGATE A NEUTRON THROUGH THE C SCINT AND PRODUCE THE LIGHT DEPOSITED BY ANY CHARGED C PARTICLES PRODUCED IN THE PROPAGATION PROCESS C C THE MAIN ROUTINE THEN BINS THE LIGHT OUTPUT OF FOLNUT C AS REQUESTED BY INPUT PARAMETERS PREVIOUSLY READ C C ------------------------------------------------------------------- C DIMENSION X(3),C(3),EL(9),S(10),XI(3),CI(3),COMMNT(20) C LT, LCON, IN, AND NEV ARE MADE FLOATING C BECAUSE OF THE LIMITED INTEGER RANGE ON PDP-11'S DIMENSION LHI(9),ELOPE(3) REAL*4 LCON(9),LT,EFF,IN,NEV VIRTUAL LT(1000,9),EFF(1000) COMMON/SCINP/PA,PB,PC,PR,PS,PT COMMON/SCINA/AG,AH,AO,AP,AQ,AU COMMON /CRBANG/ED(32),AD(32,8) COMMON /SIGDAT/ EDAT(128),SDAT(128,9) REAL*4 NDUMP,MDUMP COMMON/SWBLK/NSW(6),NDUMP,MDUMP COMMON /FIRSC/ NFS(6,5,2) COMMON/SCIN/DHYD,DCARB,IGEO,XB,YB,ZB COMMON /TIMING/ TBIAS, DTIM, DPOS, NTIM(100),NPOS(100) C TO REDUCE THE OVERHEAD IN CALLING UNIRND - ALL C REFERENCES HAVE BEEN CHANGED TO RAN(KS1,KS2) THE C RANDOM NUMBER ROUTINE N THE PDP-11 FORTRAN LIBRARY C THE FOLLOWING ARE SEEDS FOR THE RANDOM NUMBER GENERATOR. C C VAX-11 USES ONLY KS1 ! C ALL REFERENCES TO RAN(KS1,KS2) HAVE BEEN CHANGED C TO RAN(KS1). C C COMMON /RANGEN/ KS1,KS2 C CHARACTER*32 INPUT,OUTPUT,CFILE CFILE='crosec.in' C C VAX-11 DEVICE ASSIGNMENT C 77777 FORMAT(A) C 55 FORMAT (20A4) 56 FORMAT (/1X,20A4) 9001 FORMAT (6I1,F4.0) 9002 FORMAT (1X,6I1,F4.0) C READ IN CROSS SECTION DATA C EDAT(I) CONTAINS THE ENERGY (IN MEV) FOR THE ITH DATA POINT C SDAT(I,J) CONTAINS TOTAL CROSS-SECTION (IN BARNS) FOR REACTION C CHANNEL J FOR ITH DATA POINT C J=1 NP ELASTIC SCATTERING C J=2 NC NONDIFFRACTIVE ELASTIC SCATTERING C J=3 N+C+N+C+GAMMA C J=4 N+C=ALPHA+9BE C J=5 N+C=N+3ALPHAS C J=6 N+C=(N+P+11B) AND (2N+11C) OR (P+12B) AND (2N+11C) C J=7 NC DIFFRACTIVE ELASTIC SCATTERING C J=8 N2N CHANNEL (PREVIOUSLY INCLUDED IN CHANNEL 6) open(unit=3,name=cfile,status="old") DO 2 I=1,128 2 EDAT(I)=1.0E7 READ(3,55) COMMNT READ(3,1) N 1 FORMAT (I3) READ(3,55) COMMNT C (N CROSS-SECTION CARDS FOLLOW) DO 1331 I=1,N 1331 READ(3,3) EDAT(I), (SDAT(I,J),J=1,9) 603 FORMAT(16F5.3) 605 FORMAT(1X,16F8.3) 3 FORMAT (10F5.3) 1334 FORMAT (F8.2,5X,7F7.3) C READ IN COEFFICIENT ARRAY FOR NC NON DIFFRACTIVE ELASTIC SCATTERING C ED(I) CONTAINS ENERGY FOR ITH DATA POINT C AD(I,J),J=1 TO 5, CONTAINS JTH SHAPE PARAMETER FOR ITH DATA POINT READ(3,55) COMMNT READ(3,1) M READ(3,55) COMMNT C (M DATA CARDS FOLLOW) DO 1349 J=1,32 1349 ED(J)=1.0E7 DO 1350 I=1,M READ(3,4) ED(I),(AD(I,J),J=1,5) C AD(I,J),J=6 TO 8, CONTAINS JTH AREA PARAMETER FOR ITH DATA POINT AD(I,6)=2.*AD(I,1) AD(I,7)=(1.-AD(I,3))*(AD(I,2)-AD(I,1))*.5 IF (AD(I,7) .LT. 0.) AD(I,7)=0. AD(I,8)=(1.+AD(I,5))*(AD(I,4)-AD(I,1))*.5 IF (AD(I,8) .LT. 0.) AD(I,8)=0. 1350 CONTINUE C C INPUT SENSE SWITCH SETTINGS AND MAXIMUM NO. OF EVENTS TO DUMP C NSW(I), THE "SENSE SWITCHES" CONTROL HOW MUCH OUTPUT IS DUMPED READ(5,55) COMMNT WRITE(6,56) COMMNT READ(5,55) COMMNT WRITE(6,56) COMMNT READ(5,9001)NSW,MDUMP C C READ IN BIAS AND BIN WIDTH INFORMATION FOR TIMING AND POSITION STUDY READ(5,55) COMMNT READ(5,4) TBIAS, DTIM, DPOS READ(5,55) COMMNT READ(5,5) (ELOPE(I),I=2,3) ELOPE(1)=0. C READ SCINTILLATOR BOUNDARIES AND 1 PHOTOELECTRON LEVELS READ(5,55) COMMNT WRITE(6,56)COMMNT READ(5,5006)KS1,KS2 5006 FORMAT(2I5) WRITE(6,5007) KS1,KS2 5007 FORMAT(1X,2I5) READ(5,55) COMMNT WRITE(6,56) COMMNT READ(5,4) RHO, COMP WRITE(6,7) RHO, COMP READ(5,55)COMMNT WRITE(6,56)COMMNT READ(5,104)PA,PB,PC,PR,PS,PT C READ COEFFICIENTS FOR LIGHT RESPONSE FUNCTIONS FOR PROTONS WRITE(6,107)PA,PB,PC,PR,PS,PT READ(5,55)COMMNT WRITE(6,56)COMMNT READ(5,104)AG,AH,AO,AP,AQ,AU C READ COEFFICIENTS FOR LIGHT RESPONSE FUNCTIONS FOR ALPHA WRITE(6,107)AG,AH,AO,AP,AQ,AU 104 FORMAT(6F12.8) 107 FORMAT(1X,6F12.8) 4 FORMAT(10F8.3) 7 FORMAT(1X,10F8.3) DCARB=.6025*2.54*RHO/(12.+COMP) DHYD=DCARB*COMP READ(5,55) COMMNT WRITE(6,56) COMMNT READ(5,5) XB,YB,ZB,IGEO,IRANP WRITE(6,5) XB,YB,ZB,IGEO,IRANP C IGEO DEFINES SHAPE OF SCINTILLATOR C IGEO=0 GIVES RECTANGULAR C+ IGEO=1 CYLINDER NEUTRONS INCIDENT ON END C IGEO=2 CYLINDER NEUTRONS INCIDENT ON A CURVED SIDE 5 FORMAT(3F6.2,2I2) C INITIALIZE HISTOGRAM ARRAYS AND COUNTERS READ(5,55) COMMNT READ(5,55) COMMNT READ(5,55) COMMNT 1000 DO 10 I=1,9 DO 9 J=1,1000 9 LT(J,I)=0. LCON(I)=0. 10 LHI(I)=0 NDET=0 NDUMP=0. IN=0. C ZERO TIMING AND POSITION ARRAYS DO 11 I=1,100 NTIM(I)=0 11 NPOS(I)=0 DO 12 I=1,5 DO 12 J=1,6 NFS(J,I,1)=0 12 NFS(J,I,2)=0 4100 READ(5,15,END=9999)(XI(I),I=1,3),(CI(I),I=1,3),E1,DE,NEV,BINW 15 FORMAT(3F5.2,3F7.4,2F8.3,F8.0,F8.3) C PDP-11/VAX-11 PROCESSOR CLOCK CALL CWMZ PROST=SECNDS(0.) C C DE IS THE DESIRED BINWIDTH CENTERED ABOUT E1 IF(DE .LE. 0.) DE=0. 2000 IN=IN+1. DO 30 I=1,3 X(I)=XI(I) 30 C(I)=CI(I) IF (IRANP .NE. 1) GO TO 2006 IF (IGEO-1) 2002,2004,2005 2002 X(2)=YB*(1.-2.*RAN(KS1)) X(1)=XB*(1.-2.*RAN(KS1)) GO TO 2006 C IF CYLINDRICAL GEOMETRY, CHOOSE RANDOM RADIUS FROM LINEAR DISTRIBUTI 2004 R=XB*SQRT(RAN(KS1)) PHI=6.2832*RAN(KS1) X(1)=R*COS(PHI) X(2)=R*SIN(PHI) GOTO 2006 2005 X(3)=ZB*RAN(KS1) X(1)=XB*(1.-2.*RAN(KS1)) X(2)=SQRT(YB*YB-X(1)*X(1)) 2006 CONTINUE E=E1+DE*(.5-RAN(KS1)) DO 40 I=1,9 40 EL(I)=0. C CALL SUBROUTINE FOLNUT, WHICH FOLLOWS THE NEUTRON UNTIL IT EITHER C ESCAPES FROM THE SCINTILLATOR OR DROPS BELOW 0.1 MEV IN ENERGY CALL FOLNUT(X,C,E,.1,EL,MS) IF(EL(7) .LE. 0.) GO TO 46 DO 44 J=2,3 K=J+6 SIG=SQRT(EL(7)/ELOPE(J)+.5)*ELOPE(J) C GAUSS SUBROUTINE REPLACED BELOW C APPROXIMATION TO NORMAL DISTRIBUTION USING CENTRAL C LIMIT THEORM AND UNIFORM RANDOM NUMBER GENERATOR GAUSSB=-6. DO 9998 IGAUSS=1,12 9998 GAUSSB=GAUSSB+RAN(KS1) 44 EL(K)=EL(7)+SIG*GAUSSB C GAUSS SUBROUTINE REPLACED ABOVE 46 DO 50 I=1,9 IF (EL(I) .LE. 0.) GO TO 50 L=IFIX(EL(I)/BINW)+1 IF(L .GE. 1000) L=1000 LT(L,I)=LT(L,I)+1. LCON(I)=LCON(I)+1. LHI(I)=MAX0(LHI(I),L) 50 CONTINUE IF (EL(7) .GT. 0.) NDET=NDET+1 NSW(3)=NSW(3)+1 51 FORMAT (2I8) IF(IN .GE. NEV) GO TO 3000 GO TO 2000 3000 CONTINUE XBB=2.*XB YBB=2.*YB CWMZ PROST=SECNDS(PROST)/60. CWMZ IF(PROST .LT. 0.0) PROST=PROST+1440. GOTO (3002,3004,3004),IGEO+1 3002 WRITE(6,3003) ZB,XBB,YBB GO TO 3006 3004 WRITE(6,3005) ZB,XB 3003 FORMAT(///' RECTANGULAR SCINTILLATOR (INCHES) ' , 1 F8.3,' DEEP BY ',F8.3,' BY ',F8.3) 3005 FORMAT(' CYLINDRICAL SCINTILLATOR (INCHES) ' , 1 F8.3,' DEEP BY ',F8.3,' IN RADIUS ' ) 3006 CONTINUE WRITE(6,56) COMMNT WRITE(6,57) RHO, COMP 57 FORMAT (' DENSITY= ', F7.3,' GM/CM**3 ',5X,' H/C RATIO= ',F7.3) WRITE(6,300) IN,E1,DE,(XI(I),I=1,3),(CI(I),I=1,3) 300 FORMAT(1X,F8.0,' NEUTRONS OF ENERGY', F6.2,' BINW' ,F6.3,' MEV' / 1 ' INCIDENT AT X,Y,Z=' ,3F7.2, ' INCHES',3X 1 ,/,' DIRECTION COSINES= ',3F7.4/) IF (IRANP .EQ. 1) WRITE(6,301) 301 FORMAT (' ** RANDOM INITIAL POSITION ** ' ) CWMZ WRITE(6,3102)PROST 3102 FORMAT(' PROCESSSOR TIME = ',F8.2,' MINUTES.') IF (NSW(1).EQ.1) GOTO 3360 WRITE(6,3007) WRITE(6,3008) 3007 FORMAT (' DIFFERENTIAL PULSE HEIGHT SPECTRA ' /) 3008 FORMAT (' LIGHT OUTPUT IN MEV ELECTRON EQUIVALENT ' ) DO 350 I=1,6 NC=LHI(I) WRITE(6,302) I,BINW,LCON(I) IF (LCON(I) .EQ. 0.) GO TO 350 302 FORMAT( ' CHANNEL ' ,I3,' ONLY ' ,5X,' BIN WIDTH ' ,F5.3,6X,F8.0, 1 ,' EVENTS ' ) EN=LCON(I) DO 351 J=1,NC EFF(J)=EN/IN 351 EN=EN-LT(J,I) WRITE(6,304)(LT(L,I),L=1,NC) WRITE(6,3009) WRITE(6,310)(EFF(J),J=1,NC) 350 CONTINUE 304 FORMAT (1X,10F5.0) DO 360 I=1,3 K=I+6 WRITE(6,306) ELOPE(I),BINW,LCON(K) NC=LHI(K) 306 FORMAT (' LIGHT OUTPUT FROM ALL CHANNELS' ,5X, 1 ' ONE PHOTOELECTRON LEVEL AT ' , F7.3, ' MEV ' / 1 1X,' BIN WIDTH ' ,F5.3,' MEV ' ,6X,F7.0, ' EVENTS' ) 360 WRITE(6,304) (LT(L,K),L=1,NC) WRITE(6,3009) 3009 FORMAT (//1X,' INTEGRAL PULSE HEIGHT SPECTRA ' /) WRITE(6,3008) 3360 K=8 NC=LHI(K) WRITE(6,308) ELOPE(2),BINW 308 FORMAT (1X/' INTEGRAL EFFICIENCY,ONE P.E.= ' ,F7.3,' MEV ' ,3X, 1 ' BIN WIDTH ' ,F5.3,' MEV ' ) WRITE(6,309) 309 FORMAT (' NOTE FIRST BIN CORRESPONDS TO ZERO BIAS. ') EN=LCON(K) DO 370 J=1,NC EFF(J)=EN/IN 370 EN=EN-LT(J,K) WRITE (6,3333) E1,BINW,NC 3333 FORMAT (1X ,F7.2,F7.3,I7) 380 WRITE(6,310) (EFF(J),J=1,NC ) 310 FORMAT ((1X ,10(2PF7.3))) IF (NSW(1).EQ.1) GOTO 1000 WRITE(6,400) 400 FORMAT (' FIRST SCATTER DATA ' ) WRITE(6,3010) 3010 FORMAT(' BIAS ',6X,'CH.1 CH.2 CH.3 CH.4 CH.5 CH.6 ') 420 FORMAT (1X,F6.2,6X,6I6) BIAS=.125 DO 410 I=1,5 BIAS=BIAS*2. 410 WRITE(6,420) BIAS,(NFS(J,I,1),J=1,6) C OUTPUT SECOND SCATTER INFORAMATION WRITE(6,500) 500 FORMAT (' SECOND SCATTER DATA ') WRITE(6,3010) BIAS=.125 DO 510 I=1,5 BIAS=BIAS*2. 510 WRITE(6,420) BIAS,(NFS(J,I,2),J=1,6) C OUTPUT TIMING AND POSITION INFORMATION WRITE(6,520) TBIAS 520 FORMAT (' HISTOGRAMS OF TIME DELAY AND X-COORD.', 1 ' DRIFT BEFORE BIAS OF ',F5.3,' MEV ATTAINED ' ) WRITE(6,521) DTIM 521 FORMAT (' TIMING ' ,F5.3,' NSEC/BIN ' ) WRITE(6,522) NTIM 522 FORMAT ( 2(10I5,5X)) WRITE(6,523) DPOS 523 FORMAT (' POSITION ' ,F5.3,' INCH/BIN ' ) WRITE(6,522) NPOS GO TO 1000 9999 STOP END C C----------------------------------------------------------------- C FUNCTION EEQUIV(E,X,C,COSLP,L) C THIS SUBROUTINE COMPUTES THE ELECTRON EQUIVALENTS OF C LIGHT DEPOSITED IN THE SCINTILLATER, WHEN A PARTICLE C OF TYPE "L" (1=PROTON, 4=ALPHA) WITH ENERGY E C SCATTERS FROM COORDINATES X(I) WITH DIRECTION COSINES C(I) C AND SCATTERING ANGLE COSLP C C THE SCINT RESPONSE FUNC IS ELM(ENERGY,PAR1,PAR2,PAR3,PAR4) C C ------------------------------------------------------------------ C DIMENSION C(3),CN(3),X(3) COMMON/SCIN/DHYD,DCARB,IGEO,X0,Y0,Z0 COMMON/SCINP/A,B,R,S,T,V COMMON/SCINA/G,H,O,P,Q,W COMMON /RANGEN/ KS1,KS2 ELM(Y,A,B,R,S)=A*(1.-EXP(B*Y**R))+S*Y ELT=0 IF (E.LT..02) GOTO 300 IF (E.LE.0.2) GOTO 100 U=6.2832*RAN(KS1) CALL SCATTR(C,COSLP,U,CN) C FIND DISTANCE DSM TO WALL IF(CN(3) .GT. 0.0) DSM=(Z0-X(3))/CN(3) IF(CN(3) .EQ. 0.0) DSM=1.0E+6 IF(CN(3) .LT. 0.0) DSM=-X(3)/CN(3) IF (IGEO-1) 10,20,20 C FOR RECTANGULAR SCINTILLATOR DISTANCE TO WALL IS SIMPLE 10 D=1.0E+6 IF(CN(1) .NE. 0.0) D=ABS((X0-(ABS(CN(1))/CN(1))*X(1))/CN(1)) IF(D .LT. DSM) DSM=D IF(CN(2) .NE. 0.0) D=ABS((Y0-(ABS(CN(2))/CN(2))*X(2))/CN(2)) IF(D .LT. DSM) DSM=D GOTO 30 C FOR CYLINDER WE FIND DISTANCE A HARDER WAY 20 AQ=CN(1)*CN(1)+CN(2)*CN(2) BQ=2*(X(1)*CN(1)+X(2)*CN(2)) CQ=X(1)*X(1)+X(2)*X(2)-X0*X0 IF(AQ .NE. 0.0) D=(-BQ+SQRT(BQ*BQ-4.*AQ*CQ))/(2.*AQ) 30 IF(D.LE.DSM) DSM=D D=E IF (L.EQ.4) D=D*.25 D=ALOG(D) C NEXT TWO LINES FIND RANGE OF PARTICLE IN SCINTILLATOR F=(((3.1471E-3-2.321E-4*D)*D-.020364)*D+.0819283)*D D=EXP((F+1.61711)*D-3.8103)/25.4 C IF PARTICLE WONT ESCAPE SCINT I.E. D.LE.DSM CALCULATION C IS SIMPLIFIED TO JUST PLUGING ENERGY INTO ELM(E, , , , , ) IF(D.LE.DSM) GOTO 100 C RANGE IS TO BIG SO WE MUST CALCULATE ENERGY LOST C IN TRAVELING A DISTANCE D-DSM OUTSIDE OF SCINT D=ALOG(25.4*(D-DSM)) F=(((D*2.742E-5-1.8209E-4)*D-8.88849E-5)*D+1.00554E-3)*D EEXT=EXP((F+.561839)*D+2.1964) C EEXT=ENERGY FOR PRO OR ALPHA OUTSIDE SCIN IF (EEXT.LE..02) GOTO 100 IF (L.EQ.4) GOTO 150 ELT=ELM(EEXT,A,B,R,S) 100 IF (L.EQ.4) GOTO 400 ELT=ELM(E,A,B,R,S)-ELT GOTO 300 150 ELT=ELM(4.0*EEXT,G,H,O,P) 400 ELT=ELM(E,G,H,O,P)-ELT 300 IF(ELT.LE.0.) ELT=0. EEQUIV=ELT RETURN END C C--------------------------------------------------------------- C C FOLNUT IS CALLED ONCE PER INCIDENT NEUTRON TO FOLLOW THE C NEUTRON THROUGH THE SCINTILLATOR AND VARIOUS REACTIONS C ESSENTIALLY IT RETURNS ONLY ( EL ) THE LIGHT PRODUCED C C MANY SUBROUTINES THAT WERE PREVIOUSELY CODED AS SUBROUTINES C AND WERE CALLED FROM WITHIN FOLNUT HAVE BEEN RECODED AS C PART OF FOLNUT , THIS ELIMINATES MANY JUMPS TO SUBROUTINE C AND MUCH MEMORY SWAPPING AND OVER LAYING ON SOME COMPUTERS C THE RESULTING CODE EXECUTES FASTER C C----------------------------------------------------------------- C SUBROUTINE FOLNUT(X,C,E ,EMIN,EL,NS) DIMENSION X1S(3),C1S(3) DIMENSION X(3),C(3), EL(7),S(10),XG(3),CG(3),NF(5) COMMON/FIRSC/NFS(6,5,2) COMMON/SCIN/DHYD,DCARB,IGEO,XB,YB,ZB REAL*4 NDUMP,MDUMP COMMON/SWBLK/NSW(6),NDUMP,MDUMP COMMON/TIMING/TBIAS,DTIM, DPOS,NTIM(100),NPOS(100) COMMON/CRBANG/ED(32),AD(32,8) COMMON/MASS/RM1,RM2,RM3,RM4,RM2A,QQ,RM5,RM6 COMMON /RANGEN/KS1,KS2 RM1=939.55 RM2=11178. RM3=3728.3 RM4=8394.6 RM2A=11182.4 QQ=-7.26 RM5=1878.5 RM6=10255. L=7 NSMAX=10 NS=0 TIME=0. ICOUNT=0 X0=X(1) DO 4 I=1,L 4 EL(I)=0. C PREPARE TO PRINT INDIVIDUAL HISTORIES, IF DESIRED. NDUMP=NDUMP+1. IF (NDUMP .GT. MDUMP) GO TO 2000 IF (NDUMP .EQ. 1.) WRITE(6,9) 9 FORMAT( ' DUMP OF INDIVIDUAL HISTORIES' ,/, 1 1X,' THE FOLLOWING VARIABLE ARE DUMPED AT EACH INTERACTION ' ,/, 1 1X,' (SEE FOLNUT LISITNG FOR DEFINITION) ' ,/, 1 1X,' E,EL(7),MS,ELT,(X(I),I=1,3),(C(I),I=1,3),COSL,PHI,TIME ') WRITE(6,11) NDUMP,E 11 FORMAT(' EVENT ',F5.0,5X,' INITIAL ENERGY ' ,F8.3,' MEV ' ) C FIND TOTAL CROSS SECTIONS 2000 CONTINUE IF (E .LT. EMIN) RETURN CALL SIGTOT(E,S) ELT=0. C FIND DISTANCE TO NEXT SCATTER. C SHH AND SCC ARE THE TOTAL CROSS-SECTIONS IN HYDROGEN AND CARBON. C DHYD AND DCARB ARE THE NUMBER OF H AND C NUCLEI PER BARN-INCH SHH=S(1) SCC=S(2)+S(10)+S(7) PATH=1./(DHYD*SHH+DCARB*SCC) D=-PATH*ALOG(RAN(KS1)) C PROPAGATE NEUTRON A DISTANCE D. DO 20 I=1,3 20 X(I)=X(I)+C(I)*D IF(INBNDS(X) .LE. 0) RETURN NS=NS+1 TIME=TIME+D*1.84/SQRT(E) IF(NS .GT. NSMAX) RETURN C DECIDE WHETHER CARBON INTERACTION OCCURRED. IF(DCARB*SCC*PATH .GT. RAN(KS1)) GO TO 200 C ELASTIC SCATTER FROM HYDROGEN OCCURRED (REACTION CHANNEL 1) 100 MS=1 CALL NPNP(E,W,COSL,X,C,ELT) GO TO 1000 C CARBON INTERACTION OCCURED. DECIDE WHETHER ELASTIC OR INELASTIC. 200 IF (S(10)/SCC.GT.RAN(KS1)) GO TO 300 C ELASTIC SCATTER FROM CARBON OCCURRED (REACTION CHANNEL 2) MS=2 CALL NCEL(E,W,COSL,S(2),S(10),S(7),ELT) C W,COSL, ELT ARE DIFINED UNDER NPNP ABOVE GO TO 1000 C INELASTIC CARBON INTERACTION OCCURRED. CHOOSE WHICH. 300 IF (S(3)/S(10) .LE. RAN(KS1)) GO TO 400 C N+C--N+C+GAMMA OCCURRED (REACION CHANNEL 3) MS=3 CALL NCNCGM(E,W,COSL,X,C,ELT) GO TO 1000 400 SR=S(4)+S(5)+S(6)+S(8)+S(9) U=RAN(KS1) IF (S(4)/SR .LE. U) GO TO 500 C N+C--ALPHA+9BE OCCURRED (REACTION CHANNEL 4) MS=4 C START OF NALPHA-------------------------------------------------- C SIMULATE THE RACTION N+C--ALPHA+9BE C ISOTROPIC TWO-BODY REACTION ASSUMED. C DATA STATEMENT ABOVE FOR RM1-RM4 CNAL=-1.+2.*RAN(KS1) CALL KINNR(RM1,RM2,RM3,RM4,E,CNAL,VB ,WCM ,PCM ,COSL,T ) ELT=EEQUIV(T,X,C,COSL,4) C END OF NALPHA---------------------------------------------------- GO TO 950 500 IF((S(4)+S(5))/SR .LE. U) GO TO 600 C N+C--N+THREE ALPHAS OCCURRED (REACTION CHANNEL 5) MS=5 CALL NN3AL(E,W,COSL,X,C,ELT) GO TO 1000 C N+C--P+12B OCCURRED (REACTION CHANNEL 6) 600 MS=6 IF((S(4)+S(5)+S(6))/SR.LE.U) GOTO 700 CALL NPB(E,W,COSL,X,C,ELT) GOTO 1000 700 MS=6 IF ((S(4)+S(5)+S(6)+S(8))/SR.LE.U) GOTO 800 C (N,2N) , START OF------------------------------------------------- C TREAT N2N SAME AS NPB KINEMATICALLY EXCEPT FOR THRESHOLD 10222 SQUNI=SQRT(RAN(KS1)) CALL KINNR(RM1,RM2,1882.8,RM6,E,SQUNI,VB,WCM,PCM,COSL,E) CALL SCATTR(C,COSL,6.2832*RAN(KS1),C1S) DO 9223 I=1,3 9223 X1S(I)=X(I) W=E*RAN(KS1) E=E-W CALL SECNUT(X1S,C1S,E,.1,ELT,NS1S) C FOLLOW ONE NEUTRON WITH SECNUT AND THE OTHER STAYS IN FOLNUT C (N,2N) , END OF--------------------------------------------------- GOTO 1000 800 MS=6 C START OF NPEVAP------------------------------ 950 W=0. COSL=1. C INCREMENT THE ELEMENT MS=1 TO 6 OF THE LIGHT OUTPUT ARRAY ELT C CORRESPONDING TO THE CHANNEL OF INTERACTION, AND ELT(7),THE TOTAL C LIGHT OUTPUT 1000 EL(MS)=EL(MS)+ELT EL(7)=EL(7)+ELT E=W PHI=6.2832*RAN(KS1) CALL SCATTR(C,COSL,PHI,C) IF(NS-2) 810,820,1100 810 MC=MS DO 802 I=1,5 802 NF(I)=1 820 BIAS=.125 DO 910 I=1,5 BIAS=BIAS*2. IF(NF(I))910,910,902 902 IF (ELT-BIAS) 910,905,905 905 NFS(MC,I,NS)=NFS(MC,I,NS)+1 NF(I)=0 910 CONTINUE 1100 IF (ICOUNT .EQ. 1) GO TO 1110 C CHECK WHETHER TIMING BIAS IS EXCEEDED. IF SO, SET FLAG AND INCREMENT C ARRAYS. IF(ELT-TBIAS) 1110,1105,1105 1105 ICOUNT=1 JT=TIME/DTIM+1. IF(JT .GT. 100) JT=100 IF (JT .LT. 1) JT=1 JP=(X(1)-X0)/DPOS+50. IF (JP .GT. 100) JP=100 IF(JP .LT. 1) JP=1 NTIM(JT)=NTIM(JT)+1 NPOS(JP)=NPOS(JP)+1 C PRINT DETAILED INFORMATION FOR THIS EVENT 1110 IF (NDUMP .GE. MDUMP) GO TO 2000 WRITE(6,1) E,EL(L),MS,ELT,(X(I),I=1,3), (C(I),I=1,3),COSL,PHI 1,TIME 1 FORMAT(F6.1,3X,F6.2,I3,F6.2,3X,3F6.2,3X,3F6.2,3X,3F7.4) GO TO 2000 END C C------------------------------------------------------------------ C C THIS SUBROUTINE IS IDENTICAL TO FOLNUT, EXCEPT FOR A C MODIFICATION TO THE NPB ROUTINE. IT IS CALLED ONLY FROM C THE N2N ROUTINE WITHIN FOLNUT TO FOLLOW THE PATH OF THE C SECOND NEUTRON THROUGH THE SCINTILLATOR. C C ON SOME MACHINES THAT SUPPORT RECURSIVE CALLS TO SUBROUTINES C THIS SUBROUTINE CAN BE REPLACED BY A CALL TO FOLNUT WITHIN C FOLNUT ITSELF, HOWEVER CODING IN THIS FASHION ALLOWS C THE PROGRAM TO BE RUN ON ANY MACHINE EASILY C C----------------------------------------------------------------- C SUBROUTINE SECNUT(X,C,E ,EMIN,EL,NS) DIMENSION X(3),C(3),S(10),XG(3),CG(3) COMMON/SCIN/DHYD,DCARB,IGEO,XB,YB,ZB COMMON/TIMING/TBIAS,DTIM, DPOS,NTIM(100),NPOS(100) COMMON/CRBANG/ED(32),AD(32,8) COMMON/MASS/RM1,RM2,RM3,RM4,RM2A,QQ,RM5,RM6 COMMON /RANGEN/ KS1,KS2 L=7 NSMAX=10 NS=0 TIME=0. ICOUNT=0 X0=X(1) C FIND TOTAL CROSS SECTIONS 2000 CONTINUE IF (E .LT. EMIN) RETURN CALL SIGTOT(E,S) ELT=0. EL=0. C FIND DISTANCE TO NEXT SCATTER. C SHH AND SCC ARE THE TOTAL CROSS-SECTIONS IN HYDROGEN AND CARBON. C DHYD AND DCARB ARE THE NUMBER OF H AND C NUCLEI PER BARN-INCH SHH=S(1) SCC=S(2)+S(10)+S(7) PATH=1./(DHYD*SHH+DCARB*SCC) D=-PATH*ALOG(RAN(KS1)) C PROPAGATE NEUTRON A DISTANCE D. DO 20 I=1,3 20 X(I)=X(I)+C(I)*D IF(INBNDS(X) .LE. 0) RETURN NS=NS+1 TIME=TIME+D*1.84/SQRT(E) IF(NS .GT. NSMAX) RETURN C DECIDE WHETHER CARBON INTERACTION OCCURRED. IF(DCARB*SCC*PATH .GT. RAN(KS1)) GO TO 200 C ELASTIC SCATTER FROM HYDROGEN OCCURRED (REACTION CHANNEL 1) 100 MS=6 CALL NPNP(E,W,COSL,X,C,ELT) GO TO 1000 C CARBON INTERACTION OCCURED. DECIDE WHETHER ELASTIC OR INELASTIC. 200 IF (S(10)/SCC.GT.RAN(KS1)) GO TO 300 C ELASTIC SCATTER FROM CARBON OCCURRED (REACTION CHANNEL 2) CALL NCEL(E,W,COSL,S(2),S(10),S(7),ELT) C W,COSL, ELT ARE DIFINED UNDER NPNP ABOVE GO TO 1000 C INELASTIC CARBON INTERACTION OCCURRED. CHOOSE WHICH. 300 IF (S(3)/S(10) .LE. RAN(KS1)) GO TO 400 C N+C--N+C+GAMMA OCCURRED (REACION CHANNEL 3) CALL NCNCGM(E,W,COSL,X,C,ELT) GO TO 1000 400 SR=S(4)+S(5)+S(6)+S(8)+S(9) U=RAN(KS1) IF (S(4)/SR .LE. U) GO TO 500 C N+C--ALPHA+9BE OCCURRED (REACTION CHANNEL 4) C START OF NALPHA-------------------------------------------------- C SIMULATE THE RACTION N+C--ALPHA+9BE C ISOTROPIC TWO-BODY REACTION ASSUMED. C DATA STATEMENT ABOVE FOR RM1-RM4 CNAL=-1.+2.*RAN(KS1) CALL KINNR(RM1,RM2,RM3,RM4,E,CNAL,VB ,WCM ,PCM ,COSL,T ) ELT=EEQUIV(T,X,C,COSL,4) C END OF NALPHA---------------------------------------------------- GO TO 950 500 IF((S(4)+S(5))/SR .LE. U) GO TO 600 C N+C--N+THREE ALPHAS OCCURRED (REACTION CHANNEL 5) CALL NN3AL(E,W,COSL,X,C,ELT) GO TO 1000 C N+C--P+12B OCCURRED (REACTION CHANNEL 6) 600 CALL NPB(E,W,COSL,X,C,ELT) IF((S(4)+S(5)+S(6))/SR.LE.U) GOTO 950 GOTO 1000 950 W=0. COSL=1. C INCREMENT THE ELEMENT MS=1 TO 6 OF THE LIGHT OUTPUT ARRAY ELT C CORRESPONDING TO THE CHANNEL OF INTERACTION, AND ELT(7),THE TOTAL C LIGHT OUTPUT 1000 EL=EL+ELT E=W PHI=6.2832*RAN(KS1) CALL SCATTR(C,COSL,PHI,C) 1100 IF (ICOUNT .EQ. 1) GO TO 2000 C CHECK WHETHER TIMING BIAS IS EXCEEDED. IF SO, SET FLAG AND INCREMENT C ARRAYS. IF(ELT-TBIAS) 2000,1105,1105 1105 ICOUNT=1 JT=TIME/DTIM+1. IF(JT .GT. 100) JT=100 IF (JT .LT. 1) JT=1 JP=(X(1)-X0)/DPOS+50. IF(JP .GT. 100) JP=100 IF(JP .LT. 1) JP=1 NTIM(JT)=NTIM(JT)+1 NPOS(JP)=NPOS(JP)+1 END SUBROUTINE CMPSCT(NTRY,EGAM,SIGT,EEL,EGAMP,COSGAM) COMMON /RANGEN/ KS1,KS2 C COMPUTE TOTAL CROSS SECTION EPS=.001 G=EGAM/.5110 G4=(G*G)**2 A=1.+2.*G TM=2.*G*G/A B=(1.+G)/G**2 S=B*(2.*(1.+G)/A-ALOG(A)/G)+ALOG(A)*.5/G-(1.+3.*G)/A**2 SIGT=0.4989*S IF(NTRY .EQ. 0) RETURN U=RAN(KS1) C=0. D=TM 6 D=.5*D E=C+D F=G/(G-E) FL=ALOG(F)/G P=(B*(B*.5*E-FL)+.5*FL+.5*E*(F-E*G*.5)/G4)/S IF(ABS(U-P)-EPS) 20,20,10 10 IF(U-P) 6,20,12 12 C=E GO TO 6 20 EEL=0.5110*E EGAMP=EGAM-EEL COSGAM=1.-(EGAM/EGAMP-1.)/G IF (COSGAM .GT. 1.0) COSGAM=1. RETURN END FUNCTION IBNSH(A,B,N) C BINARY SEARCH FOR AN ELEMENT L OF ARRAY B SUCH THAT C A .GE. B(I) AND A .LT. B(I+1). C NOTE-DIMENSION N OF B MUST BE A POWER OF 2. DIMENSION B(2) 4 J=0 K=N 6 K=K/2 L=K+J IF(A-B(L))6,20,10 10 IF(A-B(L+1)) 20,12,12 12 J=L GO TO 6 20 IBNSH=L RETURN END SUBROUTINE FOLGAM(X,C,E,EMIN,ELIGHT) C MONTE CARLO GAMMA RAY TRACER C GIVEN A GAMMA-RAY OF ENERGY E AT POSITION X(I), MOVING ALONG DIRECTI C N COSINES C(I), FOLGAM FOLLOWS THE GAMMA THRU SUCCISSIVE COMPTON C COLLISIONS UNTIL IT ESCAPES FROM THE SCINTILLATOR OR ITS ENERGY FALLS C BELOW EMIN. FOLGAM RETURNS ELIGHT, SUM OF ELECTRON REAL ENERGIES DIMENSION X(3),C(3) COMMON/SCIN/DHYD,DCARB,IGEO,XB,YB,ZB COMMON /RANGEN/ KS1,KS2 ELIGHT=0. NSMAX=10 NS=0 C CMPSCT(0,E,SIGT,Y,Y,Y) RETURNS THE TOTAL COMPTON SCATTERING CROSS C SECTION SIGT (BARNS) PER ELECTRON FOR A GAMMA OF ENERGY E. Y IS A C DUMMY. DHYD AND 6.*DCARB ARE THE NUMBER OF HYDROGEN AND CARBON ELEC- C TRONS PER BARN-INCH 1000 CALL CMPSCT(0,E,SIGT,Y,Y,Y) PATH=1./(SIGT*(DHYD+6.*DCARB)) C PROPAGATE GAMMA A DISTANCE D TO NEXT SCATTER D=-PATH*ALOG(RAN(KS1)) DO 20 I=1,3 20 X(I)=X(I)+C(I)*D C IS GAMMA STILL IN SCINTILLATOR IF(INBNDS(X) .LE. 0) RETURN NS=NS+1 IF(NS .GT. NSMAX) RETURN C SIMULATE A COMPTON COLLISION. CALL CMPSCT(1,E,SIGT,EEL,E,COSGAM) ELIGHT=ELIGHT+EEL IF(E .LE. EMIN) RETURN PHI=6.2832*RAN(KS1) C FIND NEW DIRECTION COSINES CALL SCATTR(C,COSGAM,PHI,C) GO TO 1000 END SUBROUTINE ANGDIS(X,F,DF) C F+V IS THE INTEGRAL DISTRIBUTION, DF IS ITS DERIVATIVE. COMMON/COEF/A,B,C,V XS=X**2 F=X*(A+XS*(B+C*XS))-V DF=A+XS*(3.*B+XS*5.*C) IF (DF .LE. 0.0001) DF=.0001 RETURN END SUBROUTINE ENGDIS(X,F,DF) C PHASE SPACE ENERGY DISTRIBUTION FOR USE WITH THE NEWTONS C METHOD SUBROUTINE RTINI COMMON/ENG/U IF(X .LE. .01) GO TO 2 XR=SQRT(X)*6.562 XS=X*X DF=XR*(1.-X)**2 F=XR*X*(.6667-.8*X+.2857*XS)-U RETURN 2 DF=.6431 F=.01+DF*(X-.01)-U RETURN END SUBROUTINE SCATTR(C,COSTH,PHI, CNEW) DIMENSION C(3), CNEW(3) SI(X)=SQRT(1.-X**2) C PROTECT THE ANGLES CT=C(3) IF(CT-1.)6,4,2 2 CT=1. 4 ST=0. SP=0. CP=1. GO TO 30 6 IF(-CT-1.)20,4,8 8 CT=-1. GO TO 4 20 ST=SI(CT) CP=C(1)/ST IF(CP .GT. 1.) CP=1. IF(CP .LT. -1.) CP=-1. SP=C(2)/ST IF(SP .GT. 1.) SP=1. IF(SP .LT. -1.) SP=-1. 30 IF(COSTH .GT. 1.) COSTH=1. IF(COSTH .LT. -1.) COSTH=-1. S=SI(COSTH) CXP=S*COS(PHI) CYP=S*SIN(PHI) CZP=COSTH CNEW(1)=CXP*CT*CP-CYP*SP+CZP*ST*CP CNEW(2)=CXP*CT*SP+CYP*CP+CZP*ST*SP CNEW(3)=-CXP*ST+CZP*CT RETURN END FUNCTION INBNDS(X) C RETURNS 1 OR 0 ACCORDING TO WHETHER POSITION VECTOR X IS C INSIDE OR OUTSIDE COUNTER BOUNDARIES DIMENSION X(3) COMMON/SCIN/DHYD,DCARB,IGEO,X0,Y0,Z0 INBNDS=1 IF(X(3))9,2,2 2 IF(X(3)-Z0) 3,3,9 3 IF (IGEO-1) 32,34,34 C BOX 32 IF(ABS(X(1))-X0) 33,33,9 33 IF(ABS(X(2))-Y0) 10,10,9 C CYLINDER 34 R=SQRT(X(1)**2+X(2)**2) IF(R-X0) 10,10,9 9 INBNDS=0 10 RETURN END C C------------------------------------------------------------------ C C RTNI SOLVES NONLINEAR EQUATIONS OF THE FORM F(X)=0 C FOR THEIR ROOTS BY NEWTONS METHOD C C FCT IS THE EXTERNAL FUNCTION TO BE SOLVED C C X=ROOT, F=FUNCTION VALUE AT X C DERF=VAL OF DERIVATIVE AT X C XST=INITIAL GUESS OF X EPS=ERROR ALLOWED IN ROOT C IEND=MAX ITERATIONS ALLOWED C IER=ERROR PARAM, IER=0 THEN NO ERROR C C----------------------------------------------------------------- C SUBROUTINE RTNI(X,F,DERF,FCT,XST,EPS,IEND,IER) C PREPARE ITERATION IER=0 X=XST TOL=X CALL FCT(TOL,F,DERF) TOLF=100.*EPS C START ITERATION LOOP DO 6 I=1,IEND IF(F)1,7,1 C EQUATION IS NOT SATISFIED BY X 1 IF(DERF)2,8,2 C ITERATION IS POSSIBLE 2 DX=F/DERF X=X-DX TOL=X CALL FCT(TOL,F,DERF) C TEST ON SATISFACTORY ACCURACY TOL=EPS A=ABS(X) IF(A-1.)4,4,3 3 TOL=TOL*A 4 IF(ABS(DX)-TOL)5,5,6 5 IF(ABS(F)-TOLF)7,7,6 6 CONTINUE C END OF ITERATION LOOP C NO CONVERGENCE AFTER IEND ITERATION STEPS. ERROR RETURN. IER=1 7 RETURN C ERROR RETURN IN CASE OF ZERO DIVISOR 8 IER=2 RETURN END SUBROUTINE NCEL(E,W,COSL,S2,S7,S8,ELT) DIMENSION CG(3) COMMON/CRBANG/ED(32),AD(32,8) COMMON/MASS/RM1,RM2,RM3,RM4,RM2A,QQ,RM5,RM6 COMMON /RANGEN/ KS1,KS2 C START OF NCEL---------------------------------------------------- C COMMON CRBANG ABOVE C DECIDE WHETHER DIFFRACTION SCATTER OCCURRED IF(S8/(S2+S8) .GE. RAN(KS1)) GO TO 10350 C INTERPOLATE IN TABLE LC=IBNSH(E,ED,32) DEL =(E-ED(LC))/(ED(LC+1)-ED(LC)) DO 10302 I =1,3 CG(I)=AD(LC,I +5)+(AD(LC+1,I +5)-AD(LC,I +5))*DEL 10302 IF(CG(I) .LT. 0.) CG(I)=0. AT=CG(1)+CG(2)+CG(3) C CHOOSE PORTION OF ANGULAR DISTRIBUTION IF(CG(1)/AT .GE. RAN(KS1)) GO TO 9310 IF(CG(2)/(CG(2)+CG(3)).GE. RAN(KS1)) GOTO 10320 GOTO 10330 C RECTANGLE (ISOTROPIC) 9310 COSCM=-1.0 + 2.0*RAN(KS1) GOTO 10311 C FORWARD TRIANGLE 10320 X1 =AD(LC,3)+(AD(LC+1,3)-AD(LC,3))*DEL COSCM =X1 +(1.-X1 )*SQRT(RAN(KS1)) GOTO 10311 C BACKWARD TRIANGLE 10330 X1 =AD(LC,5)+(AD(LC+1,5)-AD(LC,5))*DEL COSCM =X1 -(X1 +1.)*SQRT(RAN(KS1)) C CALL THE KINEMATICS SUBROUTINE KINNR. 10311 IF(COSCM .GT. 1.) COSCM =1. IF(COSCM .LT.-1.) COSCM =-1. CALL KINNR(RM1,RM2,RM1,RM2,E,COSCM,VB,WCM,PCM,COSL,W ) C LIGHT OUTPUT FROM CARBON ELT=0.017*(E-W ) GOTO 9333 C DIFFRACTION ELASTIC SCATTERING. EXPONENTIAL FORWARD PEAK, WITH WIDTH C INVERSELY PROPORTIONAL TO THE TOTAL INELASTIC CROSS SECTION, AND C TO THE NEUTRON ENERGY E. 10350 COSCM =1.+ALOG(1.-RAN(KS1))/1.17/S7/E GOTO 10311 9333 RETURN END SUBROUTINE NCNCGM(E,W,COSL,X,C,ELT) COMMON/MASS/RM1,RM2,RM3,RM4,RM2A,QQ,RM5,RM6 COMMON/COEF/C1,C2,C3,V COMMON /RANGEN/ KS1,KS2 DIMENSION X(3),XG(3),C(3),CG(3) EXTERNAL ANGDIS C START OF NCNCGM-------------------------------------------------- C SIMULATE THE REACTION N+C--N+C+GAMMA C DIMENSION X(3),XG(3),C(3),CG(3) ABOVE C DATA EMU-EMJ ABOVE C FIND THE GAMMA ANGLE. THE ANGULAR DISTRIBUTION (TAKEN FROM 14 MEV C DATA) IS GIVEN BY ANGDIS, WHICH IS INVERTED BY RTNI. C C1,C2,C3, AND V ARE COEFFICIENTS FOR ANGDIS. V=RAN(KS1) C1=.7444 C2=.4342 C3=-.1786 CALL RTNI(CSG,F,D,ANGDIS,.5,.001,10,IER) IF(IER .NE. 0) WRITE(6,9111) IER,CSG,F,D 9111 FORMAT (' *** ',I3,3E12.4) IF(RAN(KS1) .GT. .5) CSG=-CSG C FOLLOW THE GAMMA C SCATTR FINDS THD DIRECTION COSINES CG(I) OF THE GAMMA FROM THE C INITIAL NEUTRON DIRECTION COSINES C(I) AND THE LAB COSINE CSG OF THE C GAMMA EMISSION ANGLE AND A UNIFORMLY-DISTRIBUTIED AZIMUTH ANGLE PHI DO 9150 IM=1,3 9150 XG(IM)=X(IM) EMM=0.1 PHI=6.2832*RAN(KS1) EGM=4.43 CALL SCATTR(C,CSG,PHI,CG) CALL FOLGAM(XG,CG,EGM,EMM,ELT) C DETERMINE THE NEUTRON ANGLE. THE ANGULAR DISTRIBUTION IS TAKEN FROM C 10 MEV PP DATA. C C1,C2,C3, AND V ARE COEFFICIENTS FOR ANGDIS. C1=.5357 C2=.8929 C3=-.4826 V=RAN(KS1) CALL RTNI (CSN,F,D,ANGDIS,.5,.001,10,IERR) IF(IERR .NE. 0) WRITE(6,9111) IERR,CSN,F,D IF(RAN(KS1) .GT. .5) CSN=-CSN C TRANSFORM NEUTRON TO LAB COORDINATES C KINNR IS THE KINEMATICS SUBROUTINE CALL KINNR(RM1,RM2,RM1,RM2A,E,CSN,VB,WCM,PCM,COSL,W ) C END OF NCNCGM---------------------------------------------------- RETURN END SUBROUTINE NN3AL(E,W,COSL,X,C,ELT) DIMENSION X(3),C(3) COMMON/MASS/RM1,RM2,RM3,RM4,RM2A,QQ,RM5,RM6 COMMON/ENG/U COMMON /RANGEN/ KS1,KS2 EXTERNAL ENGDIS C START OF NN3AL-------------------------------------------------- C SIMULATES THE REACTION N+C--N+THREE ALPHAS C PHASE SPECE ENERGY DISTRIBUTION OF OUTGOING NEUTRON ASSUMED. C DATA EMN, EMC,QQ/940.,11280.,-7.26/ WCM=E*.923+QQ IF(WCM) 2,2,15004 2 W =E COSL=1. GOTO 9500 C MAXIMUM NEUTRON ENERGY IN CM 15004 TM=.923*WCM C SIMULATE PHASE SPECE DISTRIBUTION FOR OUTGOING NEUTRON U=RAN(KS1) X3=.01 IF(U .LE. .01) GO TO 3 CALL RTNI(X3,FN,DV,ENGDIS,.5,.001,20,IER) IF(IER .NE. 0) WRITE(6,15001) IER,X3,FN,DV 15001 FORMAT ( ' ERRNN3AL',I3,3E12.4) 3 T=TM*X3 IF(T .LE. 0.) T=0. C FIND THE LAB ENERGY EN AND COSINE COSL OF THE LAB ASCATTERING ANGLE C FIND THE LAB ENERGY EN AND COSINE COSL OF THE LAB ASCATTERING ANGLE C FOR THE OUTGOING NEUTRON PCM=SQRT(2.*RM1*T) VB=SQRT(2.*RM1*E)/(RM1+RM2) C ISOTROPIC CM ANGULAR DISTRIBUTION OF THE NEUTRON IS ASSUMED. COSCM=-1.+2.*RAN(KS1) W =T+PCM*VB*COSCM+.5*RM1*VB**2 IF(W .LE. 0.) W =0. PN=SQRT(2.*RM1*W ) COSL=(PCM*COSCM+RM1*VB)/PN TR=E-W +QQ DO 15 I =1,2 TA=TR*RAN(KS1) ELT=EEQUIV(TA,X,C, RAN(KS1),4)+ELT 15 TR=TR-TA ELT=EEQUIV(TR,X,C, RAN(KS1),4)+ELT 9500 RETURN C END OF NN3AL----------------------------------------------------- END SUBROUTINE NPNP(E,W,COSL,X,C,ELT) DIMENSION X(3),C(3) COMMON/COEF/C1,C2,C3,V COMMON/MASS/RM1,RM2,RM3,RM4,RM2A,QQ,RM5,RM6 COMMON /RANGEN/ KS1,KS2 EXTERNAL ANGDIS C THE SCATTERING IS ASSUMED ISOTROPIC IN THE CM BELOW 30 7EV. ABO E C 30 MEV THE ANGULAR DISTRIBUTION BECOMES INCFEASINGLY PARABOLIC. C COMMON COEF ADDED IF(E.GT.30.) GOTO 10110 COSCM =-1.+2.*RAN(KS1) 9105 CALL KINNR (RM1,RM1,RM1,RM1,E,COSCM,VB,WCM,PCM,COSL,W) C NOTE: ARCOS(X) = 1.5703 - ARCTAN(X/SQRT(1-X^2)) COSLP=COS(ATAN2(COSL,SQRT(1-COSL*COSL))) C ELT=EEQUIV(E-W ,X,C,COSLP,1) RETURN C ABOVE 30 MEV THE ANGULAR DISTRIBUTION IS GIVEN BY ANGDIS, WITH C THE COEFFICIENTS C,D,F AND U (COMMON/COEF/). RTNI INVERTS ANGDIS. 10110 R =E/30. C1 =3./(3.+R ) C2 =R *C1 /3. C3 =0. V =RAN(KS1) CALL RTNI(COSCM ,FN ,DV ,ANGDIS,.5,.001,20,IER ) IF(RAN(KS1).GT..5) COSCM =-COSCM GOTO 9105 C END OF NPNP ROUTINE--------------------------------------------- C NPNP RETURNS THE NEW NEUTRON ENERGY W, THE COSINE OF THE LAB C SCATTERING ANGLE COSL, AND THE LIGHT OUTPUT ELT. END SUBROUTINE SIGTOT(E,S) DIMENSION S(10) COMMON/SIGDAT/EDAT(128),SDAT(128,9) C START OF SIGTOT--------------------------------------------------- C COMMON SIGDAT ADDED FOR SIGTOT ROUTINE C PURPOSE--TO RETURN THE ARRAY S(J) CONTAINING THE TOTAL CROSS SECTION C (IN BARNS) AT ENERGY E IN REACTION CHANNEL J. LSIG=IBNSH(E,EDAT,128) DEL=(E-EDAT(LSIG))/(EDAT(LSIG+1)-EDAT(LSIG)) DO 10010 I=1,9 10010 S(I)=SDAT(LSIG,I)+(SDAT(LSIG+1,I)-SDAT(LSIG,I))*DEL S(10)=S(3)+S(4)+S(5)+S(6)+S(8)+S(9) C END OF SIGTOT----------------------------------------------------- RETURN END SUBROUTINE NPB(E,W,COSL,X,C,ELT) DIMENSION X(3),C(3) COMMON/MASS/RM1,RM2,RM3,RM4,RM2A,QQ,RM5,RM6 COMMON /RANGEN/ KS1,KS2 C START OF NPB------------------------------------------------ RM9=10239.05+15.95 9224 CALL KINNR(RM1,RM2,RM5,RM9,E,SQRT(RAN(KS1)),VB,WCM,PCM,COSL,E) W=E*RAN(KS1) C PROTON ENERGY UNIFORMLY DISTRIBUTED BETWEEN 0 AND MAX VALUE C IF RESCATTERED NEUTRON ASSUME SAME SCATTERING ANGLE AS PROTON ELT=EEQUIV(W,X,C,COSL,1) C APPROX PCT N+C=N+P+B FOR RESCATTERED NEUTRON IF (RAN(KS1).LE.0.1) GO TO 9220 W=E-W C GIVE NEUTRON ENERGY NOT GIVEN TO PROTON RETURN C EEQUIV RETURNS LIGHT OUTPUT FOR PROTON OF ENERGY EN C (NO NEUTRON RELEASED) 9220 W=0. COSL=1. RETURN C END OF NPB------------------------------------------------------- END C C---------------------------------------------------------------- C C RELATIVISTIC KINEMATICS FOR TWO BODY SCATTERING PROBLEM C C ADAPTED FROM J. B. MARION AND J. L. FOWLER C "FAST NEUTRON PHYSICS, PART 1", PAGES 51-55, VOL 4 OF C INTERSCIENCE MONOGRAPHS AND TEXTS IN PHY.+ASTRO., 1960 NY.NY. C C------------------------------------------------------------------- C SUBROUTINE KINNR(EM1,EM2,EM3,EM4,T1,CSCM,VB,WCM,PCM,CSL,T3) P1=SQRT(T1*(T1+2.*EM1)) W1=T1+EM1 WCM=.5*(EM1*EM1+EM2*EM2+EM3*EM3-EM4*EM4+2.*EM2*W1) IF(WCM.LE.0.0) GOTO 10 WCM=WCM/SQRT(EM1*EM1+EM2*EM2+2.*EM2*W1) IF (WCM.LE.EM3) GOTO 10 PCM=SQRT(WCM*WCM-EM3*EM3) B0=P1/(W1+EM2) G0=1./SQRT(1.-B0*B0) T3=G0*(WCM+B0*PCM*CSCM)-EM3 IF (T3 .LT. .001) T3=.001 CSL=G0*(PCM*CSCM+B0*WCM)/SQRT(T3*(T3+2.*EM3)) IF (CSL .GT. 1.0) CSL=1.0 VB=PCM/(G0*EM3) RETURN 10 WCM=0. PCM=0. T3=0. VB=0. CSL=1. RETURN END c random number generator added by Steeve Wood--7/13/2005----Ramesh real*4 function ran(seed) integer seed integer init save init data init/0/ double precision grnd double precision result if(init.eq.0) then call sgrnd(seed) init = 1 endif result = grnd() ran = result return end