PROGRAM PERC ************************************************************************ * * * Percolation Model of Nuclear Fragmentation * * ------------------------------------------ * * Nuclei are represented as three-dimensional simple cubic * * lattices. Nucleons are at the lattice sites, and their * * interaction is represented by nearest neighbor bonds. The * * bonds are broken with a probability, PBREAK. The resulting * * connected clusters are identified with the fragments. * * * * Publications: * * W. Bauer et al., Phys. Lett. 150B, 53 (1985) * * W. Bauer et al., Nucl. Phys. A452, 699 (1986) * * W. Bauer, Phys. Rev. C 38, 1927 (1988) * * * ************************************************************************ * * * Author: Wolfgang Bauer * * Last modified: Nov-10-1994 * * * ************************************************************************ * Logical unit assignments: * * INPUT : Input data and parameters * * ISUM : Output data set for summary * * IMOM : Output data set for event-by-event moments * * ISAV : Temporary safety storage for summary * ************************************************************************ PARAMETER (INPUT=5, IMOM=15, ISUM=16, ISAV=17) PARAMETER (PI=3.141592654) * LOGICAL * 1 CLSTMB(0:11,0:11,0:11), CON(10,10,10,3), $ NONUCL(10,10,10) logical fixed * INTEGER NEWPTS(1000,3), CLSTSZ(1000), $ MULTPL(0:1000), NUCNUM(10,10,10), $ A1, A2, $ AF * INTEGER MEVENT(1000) *----------------------------------------------------------------------- * COMMON-BLOCKS * COMMON CLSTMB, NEWPTS, CON, NONUCL COMMON /C1/ A1, A2, AF, MASS, $ NX, NY, NZ, NORUNS COMMON /C2/ R1, R2, B, PBREAK, ALPHA, IRAND0 *----------------------------------------------------------------------- * INPUT * READ(INPUT,*) A1 C.. = Mass number of target READ(INPUT,*) A2 C.. = Mass number of projectile READ(INPUT,*) PBR0 C.. = Breaking proability for individual bonds READ(INPUT,*) IRAND0 IRANKP = IRAND0 C.. = Seed for random number generator READ(INPUT,*) IMPPAR C.. = 0 => Impact parameter independent PBREAK C.. = 1 => Woods-Saxon b-dependence for PBREAK C.. = 2 => PBREAK b-dependence obtained from Glauber model READ(INPUT,*) STOSS C.. = impact parameter (only if IMPPAR=0) READ(INPUT,*) BUNTEN, BOBEN C.. = Lower and upper bounds for impact parameter interval READ(INPUT,*) ALPHA ALPHA0 = ALPHA C.. = Angle of reaction plane relative to x-z plane C.. (only if IMPPAR=0) READ(INPUT,*) NORUNS C.. = Number of generated events NRN10 = NORUNS / 10 READ(INPUT,*) MCUT C.. = Upper cutoff mass for moments analysis READ(INPUT,*) MLOW, MHIGH C.. = Cutoff masses for multiplicity distributions read(input,*) ifix C.. = 0, if fixed fraction of bonds are broken, else = 1 if (ifix .eq. 0) then fixed = .true. else fixed = .false. end if read(input,*) io_mom C.. = 1 => write out lowest moments, else not read(input,*) io_sav C.. = 1 => write out intermediate histograms, else not *----------------------------------------------------------------------- * Geometry: * BMAX = 1.124 * (A1**(1./3.) + A2**(1./3.)) IF (A2 .EQ. 1) $ BMAX = 1.124 * A1**(1./3.) + 0.7 RTAR = 1.124 * A1**(1./3.) DELTA = 1.0 SURF = 0.55 BHIGH = BMAX + 1.0 NZINT = 25 ZLOW = - BHIGH ZHIGH = BHIGH DELTAZ = (ZHIGH-ZLOW) / FLOAT(NZINT) TOTAL = 0. DO 4 IZ = -NZINT,NZINT Z = FLOAT(IZ) * DELTAZ RHOZ = 1. / (1. + EXP((SQRT(Z**2) - RTAR) / SURF)) TOTAL = TOTAL + RHOZ 4 CONTINUE PNORM = TOTAL * DELTAZ * *----------------------------------------------------------------------- * Initialization of output-fields: * MULTPL(0) = 0 DO 5 I = 1,A1 CLSTSZ(I) = 0 MULTPL(I) = 0 5 CONTINUE * ====================================================================== * == == * == Main loop over events: == * == == * ====================================================================== DO 2000 NRUN = 1,NORUNS 10 CONTINUE *----------------------------------------------------------------------- * Random choice of impact parameter (if desired): * IF (IMPPAR .EQ. 0) THEN B = STOSS PBREAK = PBR0 ELSE 50 CONTINUE XXX = 1.0 - 2.0 * RAN(IRAND0) YYY = 1.0 - 2.0 * RAN(IRAND0) RRR = XXX**2 + YYY**2 IF (RRR .GT. 1.0) GOTO 50 B = SQRT(RRR) * BMAX ALPHA = 360.0 * RAN(IRAND0) IF (IMPPAR .EQ. 1) THEN PBREAK = PBR0 / ( 1.0 + EXP ( (B-RTAR) / DELTA ) ) ELSE * This is only correct for A2 = 1. For A2 not equal 1, the full * Glauber approximation has to be used. TOTAL = 0.0 DO 55 IZ = -NZINT,NZINT Z = FLOAT(IZ) * DELTAZ RHOZ = 1. / (1. + EXP((SQRT(B**2+Z**2) - RTAR) / SURF)) TOTAL = TOTAL + RHOZ 55 CONTINUE PBREAK = TOTAL * DELTAZ * PBR0 / PNORM END IF END IF IF (BUNTEN .GT. B .OR. BOBEN .LT. B) GOTO 2000 *----------------------------------------------------------------------- * Determination of the nucleon distribution on the lattice: * CALL SHAPE *----------------------------------------------------------------------- * NXP1 = NX + 1 NYP1 = NY + 1 NZP1 = NZ + 1 *----------------------------------------------------------------------- * Surface is assigned .TRUE. so that surface nucleons can be * treated in the same manner as bulk nucleons: * DO 200 IX = 1,NX DO 100 IY = 1,NY CLSTMB(IX,IY, 0 ) = .TRUE. CLSTMB(IX,IY,NZP1) = .TRUE. 100 CONTINUE 200 CONTINUE * DO 400 IX = 1,NX DO 300 IZ = 1,NZ CLSTMB(IX, 0 ,IZ) = .TRUE. CLSTMB(IX,NYP1,IZ) = .TRUE. 300 CONTINUE 400 CONTINUE * DO 600 IY = 1,NY DO 500 IZ = 1,NZ CLSTMB( 0 ,IY,IZ) = .TRUE. CLSTMB(NXP1,IY,IZ) = .TRUE. 500 CONTINUE 600 CONTINUE *----------------------------------------------------------------------- * Determination of broken bonds: * CALL BREAK(CON,NONUCL,IRAND0,PBREAK,NBONDS,NX,NY,NZ,fixed) * *----------------------------------------------------------------------- * Cluster recognition algorithm: * * 1) Initialization of array CLSTMB: * DO 1300 IZ = 1,NZ DO 1200 IY = 1,NY DO 1100 IX = 1,NX CLSTMB(IX,IY,IZ) = NONUCL(IX,IY,IZ) 1100 CONTINUE 1200 CONTINUE 1300 CONTINUE * MULTIP = 0 *----------------------------------------------------------------------- * 2) Cluster recognition and size determination: * DO 1800 IZ = 1,NZ DO 1700 IY = 1,NY DO 1600 IX = 1,NX IF(.NOT.CLSTMB(IX,IY,IZ)) THEN * * New cluster, -- start -- * CLSTMB(IX,IY,IZ) = .TRUE. NEWPTS(1,1) = IX NEWPTS(1,2) = IY NEWPTS(1,3) = IZ IMINO = 1 IMAXO = 1 * 1400 CONTINUE IMIN = IMAXO + 1 IMAX = IMIN DO 1500 I = IMINO,IMAXO IIX = NEWPTS(I,1) IIY = NEWPTS(I,2) IIZ = NEWPTS(I,3) CALL NEWTST(IIX-1,IIY ,IIZ ,IIX-1,IIY ,IIZ ,1, $ IMAX) CALL NEWTST(IIX+1,IIY ,IIZ ,IIX ,IIY ,IIZ ,1, $ IMAX) CALL NEWTST(IIX ,IIY-1,IIZ ,IIX ,IIY-1,IIZ ,2, $ IMAX) CALL NEWTST(IIX ,IIY+1,IIZ ,IIX ,IIY ,IIZ ,2, $ IMAX) CALL NEWTST(IIX ,IIY ,IIZ-1,IIX ,IIY ,IIZ-1,3, $ IMAX) CALL NEWTST(IIX ,IIY ,IIZ+1,IIX ,IIY ,IIZ ,3, $ IMAX) 1500 CONTINUE IF(IMAX .NE. IMIN) THEN IMINO = IMIN IMAXO = IMAX - 1 GOTO 1400 ENDIF * ; Increment total multiplicity IF (IMAXO .GE. MLOW .AND. IMAXO .LE. MHIGH) & MULTIP = MULTIP + 1 * ; No fragments => repeat event IF (IMAXO .EQ. A1) GOTO 10 * ; Increment mass yield CLSTSZ(IMAXO) = CLSTSZ(IMAXO)+1 * ; Store mass for moments analysis MEVENT(MULTIP) = IMAXO * * New cluster, -- end -- * END IF 1600 CONTINUE 1700 CONTINUE 1800 CONTINUE *----------------------------------------------------------------------- * Store total multiplicity of event: * MULTPL(MULTIP) = MULTPL(MULTIP) + 1 *----------------------------------------------------------------------- * Analysis of moments of mass distribution: * RMOM1 = 0.0 RMOM2 = 0.0 RMOM3 = 0.0 RMOM4 = 0.0 RMOM5 = 0.0 DO 1810 IFRG = 1,MULTIP * Calculate moments: IF (MEVENT(IFRG) .LE. MCUT) THEN RMOM1 = RMOM1 + REAL(MEVENT(IFRG)) RMOM2 = RMOM2 + REAL(MEVENT(IFRG))**2 RMOM3 = RMOM3 + REAL(MEVENT(IFRG))**3 RMOM4 = RMOM4 + REAL(MEVENT(IFRG))**4 RMOM5 = RMOM5 + REAL(MEVENT(IFRG))**5 END IF 1810 CONTINUE * Write out logarithms of moments: IF (RMOM1 .GT. 0.0 .and. io_mom .eq. 1) THEN RLOG1 = ALOG ( RMOM1 ) RLOG2 = ALOG ( RMOM2 ) RLOG3 = ALOG ( RMOM3 ) RLOG4 = ALOG ( RMOM4 ) RLOG5 = ALOG ( RMOM5 ) WRITE(IMOM,'(F10.5,I6,4X,5F10.5)') [ PBREAK, MULTIP, RLOG1, RLOG2, RLOG3, RLOG4, RLOG5 END IF *----------------------------------------------------------------------- IF (NRUN/NRN10*NRN10 .EQ. NRUN) THEN if (io_sav .eq. 1) then * Temporary safety storage of summary output (in case job crashes!) WRITE(ISAV,'(/'' ***> RUNS COMPLETED '',I10)') NRUN ABSNRM = 10.0 * PI * BMAX**2 / FLOAT(NRUN) WRITE(ISAV,'(I10,2I15)') 0, 0, MULTPL(0) DO 1950 I = 1,A1 WRITE(ISAV,'(I10,F15.6,I15)') I,ABSNRM*CLSTSZ(I),MULTPL(I) 1950 CONTINUE CLOSE(ISAV) end if write(*,'('' ==> '',i7,'' events completed'')') nrun END IF *----------------------------------------------------------------------- 2000 CONTINUE * ====================================================================== * == == * == End of events loop == * == == * ====================================================================== * ************************************************************************ * Final output * ************************************************************************ if (fixed) then write(isum,'(''( Fixed number of broken bonds'')') else write(isum,'(''( Fluctuating number of broken bonds'')') end if WRITE(ISUM,'(''( TARGET - MASS : '',I5)') A1 WRITE(ISUM,'(''( PROJECTILE - MASS: '',I5)') A2 WRITE(ISUM,'(''( PBR0: '',F13.6)') PBR0 WRITE(ISUM,'(''( IRAND0: '',I13)') IRANKP WRITE(ISUM,'(''( IMPPAR: '',I5)') IMPPAR WRITE(ISUM,'(''( B-UNTEN: '',F13.6)') BUNTEN WRITE(ISUM,'(''( B-OBEN: '',F13.6)') BOBEN WRITE(ISUM,'(''( STOSS: '',F13.6)') STOSS WRITE(ISUM,'(''( ALPHA: '',F13.6)') ALPHA0 WRITE(ISUM,'(''( MLOW; '',I5)') MLOW WRITE(ISUM,'(''( MHIGH: '',I5)') MHIGH WRITE(ISUM,'(''( NORUNS: '',I8)') NORUNS WRITE(ISUM,'(''( I MASS-YIELD[I] MULTIPLICITY[I]'')') *----------------------------------------------------------------------- * Normalization of cross section: * ABSNRM = 10.0 * PI * BMAX**2 / FLOAT(NORUNS) WRITE(ISUM,'(I10,2I15)') 0, 0, MULTPL(0) DO 5000 I = 1,A1 WRITE(ISUM,'(I10,F15.6,I15)') I, ABSNRM * CLSTSZ(I), MULTPL(I) 5000 CONTINUE C**********************************************************************C STOP END ************************************************************************ *SUB1 * SUBROUTINE NEWTST(IPX,IPY,IPZ,ICX,ICY,ICZ,IDIR,IMAX) * * * tests if the point (IPX,IPY,IPZ) has already been identified * * to belong to a cluster. If no, we test if the bond in the * * direction IDIR with the indices (ICX,ICY,ICZ) is already * * present. If yes, the point is identified to belong to a * * cluster by setting the corresponding array element of CLSTMB * * to the value .TRUE.. In addition, the counter for the number * * of points in the present cluster is incremented. The point * * is appended to the list of members of the present cluster. * * * ************************************************************************ LOGICAL * 1 CLSTMB(0:11,0:11,0:11), CON(10,10,10,3), $ NONUCL(10,10,10) * INTEGER NEWPTS(1000,3) * COMMON CLSTMB, NEWPTS, CON, NONUCL C======================================================================= IF(CLSTMB(IPX,IPY,IPZ)) THEN RETURN ELSE IF(.NOT.CON(ICX,ICY,ICZ,IDIR)) RETURN NEWPTS(IMAX,1) = IPX NEWPTS(IMAX,2) = IPY NEWPTS(IMAX,3) = IPZ IMAX = IMAX + 1 CLSTMB(IPX,IPY,IPZ) = .TRUE. RETURN END IF C======================================================================= END ************************************************************************ *SUB2 * SUBROUTINE SHAPE * * * chooses occupied lattice sites. * * * * A1, R1 : Mass and radius of target nucleus * * A2, R2 : Mass and radius of projectile nucleus * * AF : Mass of fireball trace * * B : Impact parameter * * X0, Y0, Z0 : Location of center of sphere * * XTRACE, YTRACE : Location of fireball trace * * RTRACE : Radius of trace * * MASS : Total number of spectator - nucleons * ************************************************************************ PARAMETER (PI=3.141592654) PARAMETER (XS=0.05, YS=0.07, ZS=0.13) PARAMETER (INPUT=5) *----------------------------------------------------------------------- * Declarations: * LOGICAL * 1 CLSTMB(0:11,0:11,0:11), CON(10,10,10,3), $ NONUCL(10,10,10) * INTEGER NEWPTS(1000,3), A1, $ A2, AF * REAL D(10,10,10) *----------------------------------------------------------------------- * Common-blocks: * COMMON CLSTMB, NEWPTS, CON, NONUCL COMMON /C1/ A1, A2, AF, MASS, $ NX, NY, NZ, NORUNS COMMON /C2/ R1, R2, B, PBREAK, ALPHA, IRAND0 *----------------------------------------------------------------------- * Deletion of all nucleons: * DO 300 IZ = 1,10 DO 200 IY = 1,10 DO 100 IX = 1,10 NO NUCL(IX,IY,IZ) = .TRUE. 100 CONTINUE 200 CONTINUE 300 CONTINUE *----------------------------------------------------------------------- * Treatment of the target: * R1 = (0.75*REAL(A1)/PI) ** (1./3.) NX = 2*(INT(R1)+1) NY = 2*(INT(R1)+1) NZ = 2*(INT(R1)+1) X0 = (REAL(NX)+1.)/2. + XS Y0 = (REAL(NY)+1.)/2. + YS Z0 = (REAL(NZ)+1.)/2. + ZS * DO 600 IX = 1,NX DO 500 IY = 1,NY DO 400 IZ = 1,NZ D(IX,IY,IZ) = $ SQRT((REAL(IX)-X0)**2+(REAL(IY)-Y0)**2+(REAL(IZ)-Z0)**2) 400 CONTINUE 500 CONTINUE 600 CONTINUE * DO 1000 N1 = 1,A1 DMIN = 999. DO 900 IX = 1,NX DO 800 IY = 1,NY DO 700 IZ = 1,NZ DXYZ = D(IX,IY,IZ) IF(DXYZ .LT. DMIN) THEN IXDMIN = IX IYDMIN = IY IZDMIN = IZ DMIN = DXYZ ENDIF 700 CONTINUE 800 CONTINUE 900 CONTINUE NONUCL(IXDMIN,IYDMIN,IZDMIN) = .FALSE. D(IXDMIN,IYDMIN,IZDMIN) = 9999. 1000 CONTINUE * MASS = A1 *----------------------------------------------------------------------- * Treatment of the projectile: * R2 = (0.75*REAL(A2)/PI) ** (1./3.) XTRACE = X0 + B * COS( ALPHA*PI/180. ) / 1.81 YTRACE = Y0 + B * SIN( ALPHA*PI/180. ) / 1.81 RTRACE = R2 AF = 0 DO 1300 IX = 1,NX DO 1200 IY = 1,NY DO 1100 IZ = 1,NZ DIS = SQRT((REAL(IX)-XTRACE)**2 + (REAL(IY)-YTRACE)**2) IF(DIS .LT. RTRACE .AND. .NOT. NONUCL(IX,IY,IZ)) THEN NONUCL(IX,IY,IZ) = .TRUE. AF = AF + 1 ENDIF 1100 CONTINUE 1200 CONTINUE 1300 CONTINUE * MASS = MASS - AF *----------------------------------------------------------------------- RETURN END ************************************************************************ *SUB4 * SUBROUTINE BREAK(CON,NONUCL,IRAND0,PBREAK,NBONDS,NX,NY,NZ,fixed) * * * Purpose: breaking the connections between nucleons located * * next to each other on the lattice with a breaking * * probability pbreak * * * * VERSION 1: 01.10.1985 * * VERSION 2: 21.04.1994 * * (fixed number of bonds are broken, if desired) * * Variables: * * CON - Logical-array (output) * * .TRUE. ==> Nucleons are connected * * .FALSE. ==> Bond is broken * * NONUCL - Logical-array (input) * * .TRUE. ==> No nucleon at this lattice site * * .FALSE. ==> Nucleon at this lattice site * * IRAND0 - Seed for random number generator (integer) * * PBREAK - Real number (breaking probability, input) * * NBONDS - Number of unbroken bonds (output, integer) * * NX,NY,NZ - Size of lattice (input, integer) * * fixed - Fixed number of bonds or not (input,logical) * * .true. => number of bonds fixed * * * * EXTERNALS: RAN - Random number generator (is called * * with IRAND0 and modifies it, too) * * * ************************************************************************ LOGICAL * 1 CON(10,10,10,3), NONUCL(10,10,10) logical fixed * REAL PBREAK * INTEGER IXNEXT, IYNEXT, $ IZNEXT, NBONDS, $ IRAND0, NX, $ NY, NZ, $ mix(2000), miy(2000), miz(2000), mid(2000) *----------------------------------------------------------------------- * find total number of bonds: NBONDS = 0 * DO 1000 IDIR = 1,3 IXNEXT = 0 IYNEXT = 0 IZNEXT = 0 IF(IDIR.EQ.1) IXNEXT = 1 IF(IDIR.EQ.2) IYNEXT = 1 IF(IDIR.EQ.3) IZNEXT = 1 DO 900 IZ = 1,NZ DO 800 IY = 1,NY DO 700 IX = 1,NX IF(NONUCL(IX,IY,IZ)) THEN CON(IX,IY,IZ,IDIR) = .FALSE. GOTO 700 END IF IF(NONUCL(IX+IXNEXT,IY+IYNEXT,IZ+IZNEXT)) THEN CON(IX,IY,IZ,IDIR) = .FALSE. ELSE if (fixed) then * => just connect the bond and count the total, later we * will break a fixed number of bonds randomly. CON(IX,IY,IZ,IDIR) = .TRUE. NBONDS = NBONDS + 1 mix(nbonds) = ix miy(nbonds) = iy miz(nbonds) = iz mid(nbonds) = idir else * => since the total number of bonds is not fixed now, we can * break each bond independently IF(RAN(IRAND0).LT.PBREAK) THEN CON(IX,IY,IZ,IDIR) = .FALSE. ELSE CON(IX,IY,IZ,IDIR) = .TRUE. NBONDS = NBONDS + 1 END IF end if END IF 700 CONTINUE 800 CONTINUE 900 CONTINUE 1000 CONTINUE * nbreak = nint(pbreak*nbonds) if (fixed) then do n = 1,nbreak 1100 i = nint(0.5+ran(iseed)*nbonds) if (i.le.0) goto 1100 if (i.gt.nbonds) goto 1100 1200 continue if (con(mix(i),miy(i),miz(i),mid(i))) then con(mix(i),miy(i),miz(i),mid(i)) = .false. else i = i+1 if (i.gt.nbonds) i=1 goto 1200 end if end do end if * control: count remaining bonds c Nleft = 0 c DO IDIR = 1,3 c DO IZ = 1,NZ c DO IY = 1,NY c DO IX = 1,NX c if (CON(IX,IY,IZ,IDIR)) Nleft = Nleft + 1 c end do c end do c end do c end do c write(*,'('' nbonds, nbreak, nleft: '',3i6)') c & nbonds, nbreak, nleft RETURN END