SUBROUTINE COMPT C VERSION 4.00 -- 26 JAN 1986/1900 C****************************************************************** C BUTCHER AND MESSEL'S CROSS SECTION EXPRESSION IS USED C (BUTCHER AND MESSEL, OP.CIT., P. 17-19,25), BUT THE C 1/EPSILON PART IS NOT SAMPLED IN THE WAY THAT THEY DO. C THIS ROUTINE CALLS THEIR 'EPSILON' VARIABLE BY THE NAME 'BR'. C BR=FINAL PHOTON ENERGY /INITIAL PHOTON ENERGY. C BR0 = MINIMUM BR = 1./(1.+2.*(E(NP)/RM)) C MAXIMUM BR IS 1. C BUTCHER AND MESSEL'S EXPRESSION FOR THE DIFFERENTIAL CROSS C SECTION IS PROPORTIONAL TO C (1./BR+BR)*(1.-BR*SINTHE**2/(1.+BR*BR)) C WE SHALL SAMPLE FROM THE FIRST FACTOR FROM THE INTERVAL (BR0,1) C AND USE THE SECOND FACTOR AS A REJECTION FUNCTION. C****************************************************************** DOUBLE PRECISION PEIG,PESG,PESE *KEEP,RANDPA. COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR DOUBLE PRECISION FAC,U1,U2 REAL RD(3000) INTEGER ISEED(103,10),NSEQ LOGICAL KNOR *KEEP,RUNPAR. COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB, * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN, * MONIOU,MDEBUG,NUCNUC, * CETAPE, * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL, * N1STTR,MDBASE, * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR, * FIX1I,FMUADD,FNKG,FPRINT,FDBASE * ,GHEISH,GHESIG COMMON /RUNPAC/ DSN,HOST,USER DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB REAL STEPFC INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC, * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL, * N1STTR,MDBASE INTEGER CETAPE CHARACTER*79 DSN CHARACTER*20 HOST,USER LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR, * FIX1I,FMUADD,FNKG,FPRINT,FDBASE * ,GHEISH,GHESIG *KEEP,STACKE. COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP DOUBLE PRECISION E(60),TIME(60) REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60) INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP *KEND. COMMON/THRESH/RMT2,RMSQ,ESCD2,AP,API,AE,UP,UE,TE,THMOLL COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2 DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL COMMON/ACLOCK/NCLOCK,JCLOCK C_____IF (NCLOCK.GT.JCLOCK) THEN C______WRITE(MDEBUG,* )' COMPT: NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP) C______CALL AUSGB2 C_____END IF PEIG=E(NP) EIG=PEIG EGP=EIG*RMI BR0I=1.+2.*EGP BR0=1./BR0I ALPH1=LOG(BR0I) ALPH2=EGP*(BR0I+1.)*BR0*BR0 SUMALP = ALPH1+ALPH2 371 CONTINUE CALL RMMAR(RNNO15,1,2) IF (ALPH1.GE.SUMALP*RNNO15) THEN CALL RMMAR(RNNO16,1,2) BR=EXP(ALPH1*RNNO16)*BR0 ELSE CALL RMMAR(RD,2,2) BRP=RD(1) RNNO18=RD(2) IF (EGP.GE.(EGP+1.)*RNNO18) THEN CALL RMMAR(RNNO19,1,2) BRP=MAX(BRP,RNNO19) END IF BR=((BR0I-1.)*BRP+1.)*BR0 END IF ESG=BR*EIG A1MIBR = 1.-BR TEMP=RM*A1MIBR/ESG SINTHE=MAX(0.0,TEMP*(2.0-TEMP)) CALL RMMAR(RNNO20,1,2) IF(((1.0-RNNO20)*(1.0+BR*BR).GE.BR*SINTHE))GO TO372 GO TO 371 372 CONTINUE SINTHE=SQRT(SINTHE) COSTHE=1.0-TEMP PESG=ESG PESE=PEIG-PESG+PRM ESE=PESE CALL UPHI(2,1) NP=NP+1 PSQ=ESE*ESE-RMSQ IF (PSQ.LE.0.0) THEN COSTHE=0. SINTHE=-1. ELSE COSTHE=(ESE+ESG)*A1MIBR/SQRT(PSQ) SINTHE=-SQRT(MAX(0.0,1.0-COSTHE*COSTHE)) END IF CALL UPHI(3,2) IF (ESE.LE.ESG) THEN IQ(NP)=3 E(NP)=PESE E(NP-1)=PESG ELSE IQ(NP)=1 IQ(NP-1)=3 E(NP)=PESG E(NP-1)=PESE T=U(NP) U(NP)=U(NP-1) U(NP-1)=T T=V(NP) V(NP)=V(NP-1) V(NP-1)=T T=W(NP) W(NP)=W(NP-1) W(NP-1)=T END IF RETURN END