!INDENT M3; !INDENT F2; "******************************************************************" "********************************************* KYOTO UNIVERSITY *" "*********** DEMO_PARALLEL **************************************" "******************************** EGS4 USER CODE -- APR 01 2006 *" "******************************************************************" " " " PROGRAMMER: YOSHITOMO ISHIHARA " " " "******************************************************************" " " " LAST UPDATED: JUL 17, 2007 " " " "******************************************************************" " " " F E A T U R E S " " " " - USES ENERGY CONSERVATION PROGRAM CALLED ECNSV1 " " - USES 'COUNTER' ROUTINE CALLED NTALLY " " " "******************************************************************" " " " THE FOLLOWING 'STEPS' REFER TO THE STEPS OUTLINED " " IN THE EGS3 USER MANUAL (SLAC-210). " " VARIOUS EGS USER NOTES (EUN'S) HAVE BEEN CREATED " " TO SUPPLEMENT SLAC-210 FOR THE CORRECTIONS, CHANGES " " AND ADDITIONS THAT ARE IN EGS4. " " " "******************************************************************" "******** STEP 1. USER-OVER-RIDE-OF-EGS-MACROS *******************" "******************************************************************" %C80 !NEWCONDITIONAL; "------------------------------------------------------------------" "Select random number generater: 0=RAN6 1=RANMAR " "RANMAR is a Lagged-Fibonacci Method pseudo random number generator" "devised by George Marsaglia and Arif Zaman. " "------------------------------------------------------------------" REPLACE {$RNGEN} WITH {2} "STEP 1. USER-OVER-RIDE-OF-EGS-MACROS" REPLACE {;COMIN/RANDOM/;} WITH { {SETR B=$RNGEN} [IF] {COPY B}=0 [ ;COMMON/RANDOM/URNDM(97),IXX,IXXST; ] [IF] {COPY B}=1 [ "This is ranmar.correlations (SID 1.8 last edited 18 Dec 1996)" " by Alex F Bielajew " "RANDOM VARIABLE COMMON" "RANDM0, RANNDM1, RANDM2 ARE SHADOW AREAS USED FOR CORRELATIONS" ;COMMON/RANDOMM/URNDM(97),CRNDM,CDRNDM,CMRNDM,IXX,JXX,IDUM2 ; COMMON/RANDM0/UDM0(97), CDM0, CDDM0, CMDM0, IXXDM0,JXXDM0; COMMON/RANDM1/UDM1(97), CDM1, CDDM1, CMDM1, IXXDM1,JXXDM1; COMMON/RANDM2/UDM2(97), CDM2, CDDM2, CMDM2, IXXDM2,JXXDM2; REAL URNDM,CRNDM,CDRNDM,CMRNDM,UDM0,CDM0,CDDM0,CMDM0,UDM1,CDM1, CDDM1,CMDM1,UDM2,CDM2,CDDM2,CMDM2, r4opt; INTEGER IXX,JXX,IDUM2,IXXDM0,JXXDM0,IXXDM1,JXXDM1,IXXDM2,JXXDM2, IXXIN,JXXIN; ] [IF] {COPY B}=2 [ ;COMMON/RANDOM/URNDM(97),stream; ] } REPLACE {$RANDOMSET#;} WITH { {SETR B=$RNGEN} [IF] {COPY B}=0 [ IXX=IXX*663608941;{P1}=IXX*0.23283064E-09;IF(IXX.LT.0){P1}={P1}+1.0; IF(IXX.EQ.IXXST) [OUTPUT;(' WARNING !'/ ' Same random number will be produced.'/ ' It is better to use RANNMAR random number generator.')] ] [IF] {COPY B}=1 [ {P1}=URNDM(IXX)-URNDM(JXX); IF({P1}.LT.0.) {P1}={P1}+1.; URNDM(IXX) = {P1}; IXX=IXX-1; IF(IXX.EQ.0) IXX=97; JXX=JXX-1; IF(JXX.EQ.0) JXX=97; CRNDM=CRNDM-CDRNDM; IF(CRNDM.LT.0.) CRNDM=CRNDM+CMRNDM; {P1}={P1}-CRNDM; IF({P1}.LT.0.) {P1}={P1}+1.; ] [IF] {COPY B}=2 [ {P1}=SPRNGCALL(stream); {P1}={P1}*(1-1.E-3)+1.E-3;] } "This should be called somewhere near the beginning of the main routine" "before any random numbers are asked for"; REPLACE {$RNG-INITIALIZATION;} WITH {; {SETR B=$RNGEN} [IF] {COPY B}=0 [;] [IF] {COPY B}=1 [ IXX=0; JXX=0; CALL RMARIN; DO II=1,20005[ $RANDOMSET XRANM; IF(II.GT.20000) OUTPUT (MOD(INT(XRANM*16.**JJ),16),JJ=1,7); (8X,7I3); ] ] [IF] {COPY B}=2 [ IXX=985456376; stream=init_sprng(4,MYRANK,NPROCS,IXX,SPRNG_DEFAULT); ] } "-------------ranmar.correlations---------------" "This is ranmar.correlations (SID 1.8 last edited 18 Dec 1996)" %C80 ; "**********************************************************************" " " " ranmar.correlations " " ******************* " " " " Macro set for doing restarts and correlations " " " " Coding for the EGS4 system by: " " " " Alex F Bielajew 89/12/21 Version 1.0 " " Institute for National Measurement Standards " " National Research Council of Canada " " Ottawa, CANADA " " KIA 0R6 " " " " 95/11/15 changed STORE-RNG so -ve argument does not cause bomb DR " " also eplicitly included types in RANDOM " " " " 96/08/09 added r4opt to random def'n and in RANDOMSET DR " " " " 96/12/18 Removed redundant declaration of " " $RANDOMSET " " $COMMON-RANDOM-DECLARATION-IN-BLOCK-DAT " " and " " $RNG-INITIALIZATION " " Only home for these is in ranmar.macros - BLIF " " " "**********************************************************************" ; REPLACE {$STORE-RNG(#);} WITH { ;IDUM2 = {P1}; ;IF(IDUM2.EQ.0)[ DO IDUM=1,97[UDM0(IDUM)=URNDM(IDUM);] CDM0=CRNDM;CDDM0=CDRNDM;CMDM0=CMRNDM;IXXDM0=IXX;JXXDM0=JXX; ] ELSEIF(IDUM2.EQ.-1)[ DO IDUM=1,97[UDM1(IDUM)=URNDM(IDUM);] CDM1=CRNDM;CDDM1=CDRNDM;CMDM1=CMRNDM;IXXDM1=IXX;JXXDM1=JXX; ] ELSEIF(IDUM2.EQ.-2)[ DO IDUM=1,97[UDM2(IDUM)=URNDM(IDUM);] CDM2=CRNDM;CDDM2=CDRNDM;CMDM2=CMRNDM;IXXDM2=IXX;JXXDM2=JXX; ] ELSE[ WRITE(IDUM2,*)URNDM,CRNDM,CDRNDM,CMRNDM,IXX,JXX; ] } ; REPLACE {$RESET-RNG(#);} WITH { ;IDUM2 = {P1}; ;IF(IDUM2.EQ.0)[ DO IDUM=1,97[URNDM(IDUM)=UDM0(IDUM);] CRNDM=CDM0;CDRNDM=CDDM0;CMRNDM=CMDM0;IXX=IXXDM0;JXX=JXXDM0; ] ELSEIF(IDUM2.EQ.-1)[ DO IDUM=1,97[URNDM(IDUM)=UDM1(IDUM);] CRNDM=CDM1;CDRNDM=CDDM1;CMRNDM=CMDM1;IXX=IXXDM1;JXX=JXXDM1; ] ELSEIF(IDUM2.EQ.-2)[ DO IDUM=1,97[URNDM(IDUM)=UDM2(IDUM);] CRNDM=CDM2;CDRNDM=CDDM2;CMRNDM=CMDM2;IXX=IXXDM2;JXX=JXXDM2; ] ELSE[ READ(IDUM2,*)URNDM,CRNDM,CDRNDM,CMRNDM,IXX,JXX; ] } ; "-------------end of ranmar.correlations---------------" "------------------------------------------------------------------" " PLACE CO MPILER DEPENDENT SUBROUTINE CALL HERE " "------------------------------------------------------------------" " Select Fortran Compiler. " " Lahey Fortran = 1 " " Microsoft Fortran = 2 " " Other PC = 3 (default) " " g77 on Windows = 4 " " Sun UNIX Workstation and Linux = 5 " " Other UNIX Workstation =6 " "------------------------------------------------------------------" REPLACE {$COMPILER} WITH {5} "------------------------------------------------------------------" "------------------------------------------------------------------" "Macro to select the timer to be used (compiler dependent). " "------------------------------------------------------------------" REPLACE {$TIMERSET;} WITH {;} REPLACE {$TIME-NOW#;} WITH { {SETR B=$COMPILER} [IF] {COPY B}=1 [CALL TIMER({P1});] [IF] {COPY B}=2 [CALL GETTIM(IHR,IMIN,ISEC,I100); {P1}=(IHR*3600+IMIN*60+ISEC)*100+I100;] [IF] {COPY B}=3 [{P1}=0.0;] [IF] {COPY B}=4 [TT=ETIME(TARRAY); {P1}=TARRAY(1)*100.0;] [IF] {COPY B}=5 [TT=ETIME(TARRAY); {P1}=TARRAY(1)*100.0;] [IF] {COPY B}=6 [{P1}=0.0;] } "------------------------------------------------------------------" "Macro to define DIMENSION related to the timer to be used " " (compiler dependent). " "------------------------------------------------------------------" REPLACE {$TIME-DIM;} WITH { {SETR B=$COMPILER} [IF] {COPY B}=1 [;] [IF] {COPY B}=2 [;] [IF] {COPY B}=3 [;] [IF] {COPY B}=4 [REAL TARRAY(2);] [IF] {COPY B}=5 [REAL TARRAY(2);] [IF] {COPY B}=6 [;] } "------------------------------------------------------------------" "Macro to define OPEN STATEMENT for Cross-section data. " " (UNIX or PC). " "------------------------------------------------------------------" REPLACE {$OPEN;} WITH { {SETR B=$COMPILER} [IF] {COPY B}<5 [OPEN(12,FILE='mortjob.xse',status='old'); OPEN(6,FILE='mortjob.out',status='new'); OPEN(8,FILE='mortjob.dum',status='new');] [ELSE] [OPEN(12,FILE='mortjob.xsec',status='old'); OPEN(6,FILE='mortjob.output',status='new'); OPEN(8,FILE='mortjob.dummy',status='new');] } REPLACE{$CALL-HOWNEAR(#);} WITH {$CALL-HOWNEAR-FOR-SLAC-RECT-PLANE-GEOMETRY({P1});} " **** " "THIS IS THE MACRO THAT SHOULD RETURN THE CLOSEST PERPENDICULAR " "DISTANCE TO ANY SURFACE WHIC FORMS A BOUNDARY FOR THE CURRENT " "REGION. IN THIS APPLICATION T IS REPLACED BY THE MACRO FOLLOW- " "ING WHICH IS SPECIALIZED FOR THE SLAC RECTANGULAR PLANE " "GEOMETRY. IT IS THE USER'S RESPONSIBILITY TO PROVIDE THIS MACRO" "FOR HIS OWN GEOMETRY. " "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ " ; "BUFFER FLUSH" REPLACE{$CALL-HOWNEAR-FOR-SLAC-RECT-PLANE-GEOMETRY(#);} WITH " ==================****========================" {;ZL=Z(NP);XL=X(NP);YL=Y(NP);IRL=IR(NP); IP=(IRL-2)/IXYP+1;JP=(IRL-2-IXYP*(IP-1))/$NXBIN+1; KP=IRL-1-IXYP*(IP-1)-$NXBIN*(JP-1); JP=$NZBIN+1+JP;KP=$NZBIN+$NYBIN+2+KP; ZLEFT=ZL-PCOORD(3,IP);ZRIGHT=PCOORD(3,IP+1)-ZL; {P1}=MIN(ZLEFT,ZRIGHT); YLEFT=YL-PCOORD(2,JP);YRIGHT=PCOORD(2,JP+1)-YL; {P1}=MIN({P1},YLEFT,YRIGHT); XLEFT=XL-PCOORD(1,KP);XRIGHT=PCOORD(1,KP+1)-XL; {P1}=MIN({P1},XLEFT,XRIGHT);} "THIS ROUTINE IS INTENDED TO BE USED TO CALCULATE THE MINIMUM " "PERPENDICLAR TO THE NEAREST BOUNDING SURFACE. THIS VERSION IS " "SPECIALLY DESIGNED FOR THE SLAC PLANE GEOMETRY PACKAGE. " "A DIFFERENT VERSION IS NEEDED FOR OTHER GEOMETRY PACKAGES. " "GEOMETRICAL INFORMATION" "SLAC DEFINITION OF /GEOM/ AND RE-DEFINITIONS FOR /PLADTA/ " "**** " REPLACE {;COMIN/GEOM/;} WITH {;COMIN/PLADTA/;} REPLACE {;COMIN/PLADTA/;} WITH " ================" {;COMMON/PLADTA/PCOORD(3,$MXPLNS),PNORM(3,$MXPLNS),IXYP;} "IXYP must be defined at MAIN. " "COMMON to define variables to score at AUSGAB" "DEPE:deposited energy inside the detector" "DELTAE:energy bin width in MeV" "SPG:Gamma spectrum, SPE:Electron spectrum, SPP:Positron spectrum" REPLACE {;COMIN/TOTALS/;} WITH {;COMMON/TOTALS/PROFILED(27,5,$NXBIN),DEPTHD(27,$NZBIN),IPRIM;} "COMMON of print-out parameter" REPLACE {;COMIN/LINES/;} WITH {;COMMON/LINES/NLINES,NWRITE,NCOUNT,ILINES;} "COMMON of geometry related parameter" REPLACE {;COMIN/PASSIT/;} WITH {;COMMON/PASSIT/NREG,IXY,NPLAN,NXPLAN,NYPLAN,NZPLAN;} REPLACE {$USER-CONTROLS-NEGATIVE-USTEP;} WITH " ============================" {;IF(USTEP<-1.E-2)[ IERUST=IERUST+1;OUTPUT IERUST,USTEP,IR(NP),IRNEW, IROLD,X(NP),Y(NP),Z(NP),SQRT(X(NP)**2+Y(NP)**2); (I6,' NEGATIVE USTEP=',E12.6,' IR,IRNEW,IROLD=',3I6,'X,Y,Z,R=', 4E10.3); IF(IERUST>1000)[OUTPUT;(///'0CALLED EXIT, TOO MANY USTEP ERRORS'///); STOP; "CALL EXIT(0); changed to STOP to avoid MS-Fortran error 5SEP95/YN"]] USTEP=0.0;} PARAMETER $VOXELXIN=1.0; "VOXEL DIMENSION IN X-DIRECTION (IN-FIELD)" PARAMETER $VOXELXOUT=0.25; "VOXEL DIMENSION IN X-DIRECTION (OUT-FIELD)" PARAMETER $VOXELY=0.5; "VOXEL DIMENSION IN Y-DIRECTION" PARAMETER $VOXELZ=0.5; "VOXEL DIMENSION IN Z-DIRECTION" PARAMETER $FIELDX=7; "SIZE TO MEASURE DOSE PROFILES IN $VOXELXIN SPACING" "THE NUMBER OF $FIELDX/$VOXELXIN MUST BE ODD" PARAMETER $SSD=100; "SOURCE-TO-SURFACE DISTANCE" PARAMETER $DEPTH1=3; "SCORING PROFILES" PARAMETER $DEPTH2=5; "SCORING PROFILES" PARAMETER $DEPTH3=10; "SCORING PROFILES" PARAMETER $DEPTH4=15; "SCORING PROFILES" PARAMETER $DEPTH5=20; "SCORING PROFILES" PARAMETER $PHANTOMX=40; "PHANTOM DIMENSION IN X-DIRECTION" PARAMETER $PHANTOMY=40; "PHANTOM DIMENSION IN Y-DIRECTION" PARAMETER $PHANTOMZ=40; "PHANTOM DIMENSION IN Z-DIRECTION" PARAMETER $NXBIN=(($PHANTOMX-$FIELDX)/$VOXELXOUT+$FIELDX/$VOXELXIN); "NUMBER OF BINS IN X-DIRECTION" PARAMETER $NXBININ=($FIELDX/$VOXELXIN); "NUMBER OF BINS IN X-DIRECTION (COARSE)" PARAMETER $NXBINOUT=(($PHANTOMX-$FIELDX)/(2*$VOXELXOUT)); "NUMBER OF BINS IN X-DIRECTION (TOUGH)" PARAMETER $NYBIN=3; "NUMBER OF BINS IN Y-DIRECTION" PARAMETER $NZBIN=($PHANTOMZ/$VOXELZ); "NUMBER OF BINS IN Z-DIRECTION" PARAMETER $IE=8; PARAMETER $NOFRECY=40; PARAMETER $MXREG=50000; PARAMETER $MXPLNS=500; "******************************************************************" "******************* ADDITIONAL (NON-EGS) MACROS ******************" "******************************************************************" REPLACE {;COMIN/MPI/;} WITH {;COMMON/MPI/NPROCS,MYRANK;} "******************************************************************" "************************** DECLARATIONS **************************" "******************************************************************" ;COMIN/DEBUG,BOUNDS,BREMPR,EDGE,ELECIN,ETALY1,GEOM,LINES,MEDIA,MISC, NTALY1,PASSIT,RANDOM,STACK,THRESH,TOTALS,UPHIOT,USEFUL,USER,MPI/; DIMENSION PHASE($IE),DEPTHDD(27,$NZBIN),DEPTHDD2(27,$NZBIN), DEPTHDDSUM($NZBIN),DEPTHDD2SUM($NZBIN),SUDEPTHD($NZBIN),RELASUD($NZBIN), PROFILEDD(27,5,$NXBIN),PROFILEDD2(27,5,$NXBIN),PROFILEDDSUM(5,$NXBIN), PROFILEDD2SUM(5,$NXBIN),SUPROFILED(5,$NXBIN),RELASUP(5,$NXBIN),NEVENT(27), NINCELECT(27); REAL*8 TOTKE,AVAILE; "NEEDED FOR ENERGY CONSERVATION TABULATION" CHARACTER*15 AAA,BBB,CCC,DDD,EEE,FFF,GGG,HHH,III,JJJ; $TYPE MEDARR(24,2); DATA MEDARR/$S'WATER',19*' ', $S'AIR',21*' '/; "******************************************************************" "*************** START OF EXECUTABLE CODE *************************" "******************************************************************" INCLUDE 'mpif.h'; #include ; SPRNG_POINTER stream; CALL MPI_INIT(IERR); CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NPROCS,IERR); CALL MPI_COMM_RANK(MPI_COMM_WORLD,MYRANK,IERR); IF(MYRANK.EQ.0)[ OPEN(8,file='/home/egs4-2/mortjob0.dummy'); OPEN(13,file='/home/egs4-2/mortjob0.output');] ELSEIF(MYRANK.EQ.1) [OPEN(8,file='/home/egs4-2/mortjob01.dummy');] ELSEIF(MYRANK.EQ.2) [OPEN(8,file='/home/egs4-2/mortjob02.dummy');] ELSEIF(MYRANK.EQ.3) [OPEN(8,file='/home/egs4-2/mortjob03.dummy');] ELSEIF(MYRANK.EQ.4) [OPEN(8,file='/home/egs4-2/mortjob04.dummy');] ELSEIF(MYRANK.EQ.5) [OPEN(8,file='/home/egs4-2/mortjob05.dummy');] ELSEIF(MYRANK.EQ.6) [OPEN(8,file='/home/egs4-2/mortjob06.dummy');] ELSEIF(MYRANK.EQ.7) [OPEN(8,file='/home/egs4-2/mortjob07.dummy');] ELSEIF(MYRANK.EQ.8) [OPEN(8,file='/home/egs4-2/mortjob08.dummy');] ELSEIF(MYRANK.EQ.9) [OPEN(8,file='/home/egs4-2/mortjob09.dummy');] ELSEIF(MYRANK.EQ.10) [OPEN(8,file='/home/egs4-2/mortjob10.dummy');] ELSEIF(MYRANK.EQ.11) [OPEN(8,file='/home/egs4-2/mortjob11.dummy');] ELSEIF(MYRANK.EQ.12) [OPEN(8,file='/home/egs4-2/mortjob12.dummy');] ELSEIF(MYRANK.EQ.13) [OPEN(8,file='/home/egs4-2/mortjob13.dummy');] ELSEIF(MYRANK.EQ.14) [OPEN(8,file='/home/egs4-2/mortjob14.dummy');] ELSEIF(MYRANK.EQ.15) [OPEN(8,file='/home/egs4-2/mortjob15.dummy');] ELSEIF(MYRANK.EQ.16) [OPEN(8,file='/home/egs4-2/mortjob16.dummy');] ELSEIF(MYRANK.EQ.17) [OPEN(8,file='/home/egs4-2/mortjob17.dummy');] ELSEIF(MYRANK.EQ.18) [OPEN(8,file='/home/egs4-2/mortjob18.dummy');] ELSEIF(MYRANK.EQ.19) [OPEN(8,file='/home/egs4-2/mortjob19.dummy');] ELSEIF(MYRANK.EQ.20) [OPEN(8,file='/home/egs4-2/mortjob20.dummy');] ELSEIF(MYRANK.EQ.21) [OPEN(8,file='/home/egs4-2/mortjob21.dummy');] ELSEIF(MYRANK.EQ.22) [OPEN(8,file='/home/egs4-2/mortjob22.dummy');] ELSEIF(MYRANK.EQ.23) [OPEN(8,file='/home/egs4-2/mortjob23.dummy');] ELSEIF(MYRANK.EQ.24) [OPEN(8,file='/home/egs4-2/mortjob24.dummy');] ELSEIF(MYRANK.EQ.25) [OPEN(8,file='/home/egs4-2/mortjob25.dummy');] ELSEIF(MYRANK.EQ.26) [OPEN(8,file='/home/egs4-2/mortjob26.dummy');] ELSEIF(MYRANK.EQ.27) [OPEN(8,file='/home/egs4-2/mortjob27.dummy');] OPEN(12,file='waterphantom.data',status='old'); CALL MPI_BARRIER(MPI_COMM_WORLD,IERR); TS=MPI_WTIME(); OUTPUT NPROCS; (I10); "******************************************************************" "******** STEP 2. PRE-HATCH-CALL-INITIALIZATION COMES NEXT *******" "******************************************************************" NMED=2; "NUMBER OF MEDIA" DO J=1,NMED [ DO I=1,24 [MEDIA(I,J)=MEDARR(I,J);]] NPLAN=$NXBIN+$NYBIN+($NZBIN+1)+3; "ACTUAL NUMBER OF PLANES" IXY=$NXBIN*$NYBIN; IXYZ=IXY*$NZBIN; IXYP=IXY; "PRESTA" NREG=IXYZ+7; "ACTUAL NUMBER OF REGION USED" "SET MEDIUM INDEX FOR EACH REGION" DO I=1,NREG [MED(I)=2;] "Set all region to the air at firtst" DO I=3,NREG-5 [MED(I)=1;] "PHANTOM MATERIAL" "Change the cut-off energy" DO I=1,NREG [ECUT(I)=0.7; PCUT(I)=0.01;] "******************************************************************" "***************** STEP 3. HATCH-CALL COMES NEXT *****************" "******************************************************************" CALL HATCH; "OUTPUT VARIOUS QUANTITIES ASSOCIATED WITH THE MEDIA" "OUTPUT; ('1QUANTITIES ASSOCIATED WITH EACH MEDIA:',//); DO J=1,NMED [ OUTPUT (MEDIA(I,J),I=1,24); (/,1X,24A1); OUTPUT RHO(J),RLC(J); (5X,' RHO=',G15.7,' G/CM**3 RLC=', G15.7,' CM'); OUTPUT AE(J),UE(J); (5X,' AE=',G15.7,' MEV UE=',G15.7,' MEV'); OUTPUT AP(J),UP(J); (5X,' AP=',G15.7,' MEV UP=',G15.7,' MEV'); ] OUTPUT;(/' INFORMATION OF MEDIUM AND CUT-OFFFOR EACH REGION'//); DO I=1,NREG [ IF(MED(I).EQ.0) [OUTPUT I,ECUT(I),PCUT(I); (' MEDIUM(',I3,')=VACUUM',18X,'ECUT=',G10.5,' MEV, PCUT=',G10.5,' MEV'); ] ELSE [OUTPUT I,(MEDIA(II,MED(I)),II=1,24),ECUT(I),PCUT(I); (' MEDIUM(',I3,')=',24A1,'ECUT=',G10.5,' MEV, PCUT=',G10.5,' MEV');] ]" "******************************************************************" "*********** STEP 4. HOWFAR-INITIALIZATION COMES NEXT ************" "******************************************************************" "DEFINITION OF PLANES" "SET ALL COORDINATES AND NORMALS TO ZERO TO BEGIN WITH" DO J=1,NPLAN [ PCOORD(1,J)=0.0; PCOORD(2,J)=0.0; PCOORD(3,J)=0.0; PNORM(1,J)=0.0; PNORM(2,J)=0.0; PNORM(3,J)=0.0; ] ID1=2; "First Z plane #" ID2=ID1+$NZBIN; "Last Z plane #" ID3=ID2+1; "First Y plane #" ID4=ID3+$NYBIN; "Last Y plane #" ID5=ID4+1; "First X plane #" ID6=ID5+$NXBIN; "Last X plane #" NYPLAN=ID2+1; NXPLAN=ID4+1; DO I=1,ID2 [PNORM(3,I)=1.0;] DO I=ID3,ID4 [PNORM(2,I)=1.0;] DO I=ID5,ID6 [PNORM(1,I)=1.0;] "NOW PUT IN THE EXCEPTIONS" "Z-direction" PCOORD(3,1)=-($SSD-53.8718); PCOORD(3,ID1)=0; PCOORD(3,ID1+1)=PCOORD(3,ID1)+$VOXELZ*0.5; DO I=ID1+2,ID2 [PCOORD(3,I)=PCOORD(3,I-1)+$VOXELZ;] "Y-direction" PCOORD(2,ID3)=-$PHANTOMY*0.5; PCOORD(2,ID3+1)=-$VOXELY*0.5; PCOORD(2,ID3+2)=$VOXELY*0.5; PCOORD(2,ID4)=$PHANTOMY*0.5; "X-direction" PCOORD(1,ID5)=-$PHANTOMX*0.5; DO I=ID5+1,ID5+$NXBINOUT [PCOORD(1,I)=PCOORD(1,I-1)+$VOXELXOUT;] DO I=ID5+$NXBINOUT+1,ID5+$NXBINOUT+$NXBININ [ PCOORD(1,I)=PCOORD(1,I-1)+$VOXELXIN;] DO I=ID5+$NXBINOUT+$NXBININ+1,ID6 [PCOORD(1,I)=PCOORD(1,I-1)+$VOXELXOUT;] "OUTPUT; ('1PCOORD AND PNORM VALUES FOR EACH J-PLANE (I=1,3):',//); DO J=1,NPLAN [ OUTPUT J,(PCOORD(I,J),I=1,3),(PNORM(I,J),I=1,3); (I5,6G15.7);]" "******************************************************************" "********** STEP 5. INITIALIZATION FOR AUSGAB COMES NEXT *********" "******************************************************************" CALL ECNSV1(0,NREG,TOTKE);" INITIALIZE ESUM ARRAY FOR ENERGY" " CONSERVATION CALCULATION." " NREG=NUMBER OF REGIONS" " TOTKE=TOTAL KE (DUMMY VARIABLE HERE)" " (MUST BE REAL*8)" CALL NTALLY(0,NREG); NCOUNT=0; "PARTICLE HISTORY COUNTER" ILINES=0; "INITIALIZE LINE-OUTPUT COUNTER" NEVENT(MYRANK)=0.0; DO I=1,$NZBIN [ DEPTHD(MYRANK,I)=0.0; DEPTHDD(MYRANK,I)=0.0; DEPTHDD2(MYRANK,I)=0.0;] DO J=1,5 [ DO K=1,$NXBIN [ PROFILED(MYRANK,J,K)=0.0; PROFILEDD(MYRANK,J,K)=0.0; PROFILEDD2(MYRANK,J,K)=0.0;]] "******************************************************************" "****** STEP 6. DETERMINATION OF INCIDENT PARTICLE PROPERTIES ****" "******************************************************************" ECUTMN=ECUT(2); EK0=10; "*PRESTA*" $PRESTA-INPUTS; "INPUT THE *PRESTA* VARIABLES" IF(MYRANK.EQ.0) [OPEN(7,file='/home/egs4-2/mlc0.output');] ELSEIF(MYRANK.EQ.1) [OPEN(7,file='/home/egs4-2/mlc01.output');] ELSEIF(MYRANK.EQ.2) [OPEN(7,file='/home/egs4-2/mlc02.output');] ELSEIF(MYRANK.EQ.3) [OPEN(7,file='/home/egs4-2/mlc03.output');] ELSEIF(MYRANK.EQ.4) [OPEN(7,file='/home/egs4-2/mlc04.output');] ELSEIF(MYRANK.EQ.5) [OPEN(7,file='/home/egs4-2/mlc05.output');] ELSEIF(MYRANK.EQ.6) [OPEN(7,file='/home/egs4-2/mlc06.output');] ELSEIF(MYRANK.EQ.7) [OPEN(7,file='/home/egs4-2/mlc07.output');] ELSEIF(MYRANK.EQ.8) [OPEN(7,file='/home/egs4-2/mlc08.output');] ELSEIF(MYRANK.EQ.9) [OPEN(7,file='/home/egs4-2/mlc09.output');] ELSEIF(MYRANK.EQ.10) [OPEN(7,file='/home/egs4-2/mlc10.output');] ELSEIF(MYRANK.EQ.11) [OPEN(7,file='/home/egs4-2/mlc11.output');] ELSEIF(MYRANK.EQ.12) [OPEN(7,file='/home/egs4-2/mlc12.output');] ELSEIF(MYRANK.EQ.13) [OPEN(7,file='/home/egs4-2/mlc13.output');] ELSEIF(MYRANK.EQ.14) [OPEN(7,file='/home/egs4-2/mlc14.output');] ELSEIF(MYRANK.EQ.15) [OPEN(7,file='/home/egs4-2/mlc15.output');] ELSEIF(MYRANK.EQ.16) [OPEN(7,file='/home/egs4-2/mlc16.output');] ELSEIF(MYRANK.EQ.17) [OPEN(7,file='/home/egs4-2/mlc17.output');] ELSEIF(MYRANK.EQ.18) [OPEN(7,file='/home/egs4-2/mlc18.output');] ELSEIF(MYRANK.EQ.19) [OPEN(7,file='/home/egs4-2/mlc19.output');] ELSEIF(MYRANK.EQ.20) [OPEN(7,file='/home/egs4-2/mlc20.output');] ELSEIF(MYRANK.EQ.21) [OPEN(7,file='/home/egs4-2/mlc21.output');] ELSEIF(MYRANK.EQ.22) [OPEN(7,file='/home/egs4-2/mlc22.output');] ELSEIF(MYRANK.EQ.23) [OPEN(7,file='/home/egs4-2/mlc23.output');] ELSEIF(MYRANK.EQ.24) [OPEN(7,file='/home/egs4-2/mlc24.output');] ELSEIF(MYRANK.EQ.25) [OPEN(7,file='/home/egs4-2/mlc25.output');] ELSEIF(MYRANK.EQ.26) [OPEN(7,file='/home/egs4-2/mlc26.output');] ELSEIF(MYRANK.EQ.27) [OPEN(7,file='/home/egs4-2/mlc27.output');] READ(7,*) NPROCJAW,NPS1USED,NPS2,NINCELECT(MYRANK); "NPRIM: Number of particles in clinac.ps used to calculate mlc.ps" "NCASES: Number of particles in mlc.ps" IF(MYRANK.EQ.0) [OPEN(10,file='/home/egs4-2/mlc0.ps');] ELSEIF(MYRANK.EQ.1) [OPEN(10,file='/home/egs4-2/mlc01.ps');] ELSEIF(MYRANK.EQ.2) [OPEN(10,file='/home/egs4-2/mlc02.ps');] ELSEIF(MYRANK.EQ.3) [OPEN(10,file='/home/egs4-2/mlc03.ps');] ELSEIF(MYRANK.EQ.4) [OPEN(10,file='/home/egs4-2/mlc04.ps');] ELSEIF(MYRANK.EQ.5) [OPEN(10,file='/home/egs4-2/mlc05.ps');] ELSEIF(MYRANK.EQ.6) [OPEN(10,file='/home/egs4-2/mlc06.ps');] ELSEIF(MYRANK.EQ.7) [OPEN(10,file='/home/egs4-2/mlc07.ps');] ELSEIF(MYRANK.EQ.8) [OPEN(10,file='/home/egs4-2/mlc08.ps');] ELSEIF(MYRANK.EQ.9) [OPEN(10,file='/home/egs4-2/mlc09.ps');] ELSEIF(MYRANK.EQ.10) [OPEN(10,file='/home/egs4-2/mlc10.ps');] ELSEIF(MYRANK.EQ.11) [OPEN(10,file='/home/egs4-2/mlc11.ps');] ELSEIF(MYRANK.EQ.12) [OPEN(10,file='/home/egs4-2/mlc12.ps');] ELSEIF(MYRANK.EQ.13) [OPEN(10,file='/home/egs4-2/mlc13.ps');] ELSEIF(MYRANK.EQ.14) [OPEN(10,file='/home/egs4-2/mlc14.ps');] ELSEIF(MYRANK.EQ.15) [OPEN(10,file='/home/egs4-2/mlc15.ps');] ELSEIF(MYRANK.EQ.16) [OPEN(10,file='/home/egs4-2/mlc16.ps');] ELSEIF(MYRANK.EQ.17) [OPEN(10,file='/home/egs4-2/mlc17.ps');] ELSEIF(MYRANK.EQ.18) [OPEN(10,file='/home/egs4-2/mlc18.ps');] ELSEIF(MYRANK.EQ.19) [OPEN(10,file='/home/egs4-2/mlc19.ps');] ELSEIF(MYRANK.EQ.20) [OPEN(10,file='/home/egs4-2/mlc20.ps');] ELSEIF(MYRANK.EQ.21) [OPEN(10,file='/home/egs4-2/mlc21.ps');] ELSEIF(MYRANK.EQ.22) [OPEN(10,file='/home/egs4-2/mlc22.ps');] ELSEIF(MYRANK.EQ.23) [OPEN(10,file='/home/egs4-2/mlc23.ps');] ELSEIF(MYRANK.EQ.24) [OPEN(10,file='/home/egs4-2/mlc24.ps');] ELSEIF(MYRANK.EQ.25) [OPEN(10,file='/home/egs4-2/mlc25.ps');] ELSEIF(MYRANK.EQ.26) [OPEN(10,file='/home/egs4-2/mlc26.ps');] ELSEIF(MYRANK.EQ.27) [OPEN(10,file='/home/egs4-2/mlc27.ps');] READ(10,*) IPRIMOLD; "Primary history number of the first incident particle" CLOSE(10); IDINC=-1; "AN IDENTIFIER (LIKE IARG) TO MARK INCIDENT PARTICLES" IXXST=17847465; IXX=IXXST; "INITIALIZED RANDOM NUMBER WITH STARTING SEED" $RNG-INITIALIZATION; NOFRECY=$NOFRECY; "NUMBER OF PHASE SPACE RECYCLES" NWRITE=10; "NUMBER OF INCIDENT CASES TO PRINT OUT" NLINES=15; "NUMBER OF LINES TO PRINT OUT" "******************************************************************" "******************** STEP 7. SHOWER-CALL---NEXT ******************" "******************************************************************" IF(MYRANK.EQ.0) [OPEN(10,file='/home/egs4-2/mlc0.ps');] ELSEIF(MYRANK.EQ.1) [OPEN(10,file='/home/egs4-2/mlc01.ps');] ELSEIF(MYRANK.EQ.2) [OPEN(10,file='/home/egs4-2/mlc02.ps');] ELSEIF(MYRANK.EQ.3) [OPEN(10,file='/home/egs4-2/mlc03.ps');] ELSEIF(MYRANK.EQ.4) [OPEN(10,file='/home/egs4-2/mlc04.ps');] ELSEIF(MYRANK.EQ.5) [OPEN(10,file='/home/egs4-2/mlc05.ps');] ELSEIF(MYRANK.EQ.6) [OPEN(10,file='/home/egs4-2/mlc06.ps');] ELSEIF(MYRANK.EQ.7) [OPEN(10,file='/home/egs4-2/mlc07.ps');] ELSEIF(MYRANK.EQ.8) [OPEN(10,file='/home/egs4-2/mlc08.ps');] ELSEIF(MYRANK.EQ.9) [OPEN(10,file='/home/egs4-2/mlc09.ps');] ELSEIF(MYRANK.EQ.10) [OPEN(10,file='/home/egs4-2/mlc10.ps');] ELSEIF(MYRANK.EQ.11) [OPEN(10,file='/home/egs4-2/mlc11.ps');] ELSEIF(MYRANK.EQ.12) [OPEN(10,file='/home/egs4-2/mlc12.ps');] ELSEIF(MYRANK.EQ.13) [OPEN(10,file='/home/egs4-2/mlc13.ps');] ELSEIF(MYRANK.EQ.14) [OPEN(10,file='/home/egs4-2/mlc14.ps');] ELSEIF(MYRANK.EQ.15) [OPEN(10,file='/home/egs4-2/mlc15.ps');] ELSEIF(MYRANK.EQ.16) [OPEN(10,file='/home/egs4-2/mlc16.ps');] ELSEIF(MYRANK.EQ.17) [OPEN(10,file='/home/egs4-2/mlc17.ps');] ELSEIF(MYRANK.EQ.18) [OPEN(10,file='/home/egs4-2/mlc18.ps');] ELSEIF(MYRANK.EQ.19) [OPEN(10,file='/home/egs4-2/mlc19.ps');] ELSEIF(MYRANK.EQ.20) [OPEN(10,file='/home/egs4-2/mlc20.ps');] ELSEIF(MYRANK.EQ.21) [OPEN(10,file='/home/egs4-2/mlc21.ps');] ELSEIF(MYRANK.EQ.22) [OPEN(10,file='/home/egs4-2/mlc22.ps');] ELSEIF(MYRANK.EQ.23) [OPEN(10,file='/home/egs4-2/mlc23.ps');] ELSEIF(MYRANK.EQ.24) [OPEN(10,file='/home/egs4-2/mlc24.ps');] ELSEIF(MYRANK.EQ.25) [OPEN(10,file='/home/egs4-2/mlc25.ps');] ELSEIF(MYRANK.EQ.26) [OPEN(10,file='/home/egs4-2/mlc26.ps');] ELSEIF(MYRANK.EQ.27) [OPEN(10,file='/home/egs4-2/mlc27.ps');] DO NOFPART=1,NPS2 [ "START OF SHOWER CALL LOOP" "OUTPUT NOFBAT,I,IXX;(' NOFBAT, I=',2I10,' IXX=',I12);" "Read the PS file" READ(10,*) IPRIM,(PHASE(IE),IE=1,8); IQI=PHASE(1); EI=PHASE(2); XI=PHASE(3); YI=PHASE(4); ZI=PCOORD(3,1); UI=PHASE(5); VI=PHASE(6); WI=PHASE(7); WTI=PHASE(8); "OUTPUT IQI,EI,XI,YI,ZI,UI,VI,WI,IRI,WTI; (10G15.5);" IF(NOFPART.LT.NPS2) [ IF(IPRIM.NE.IPRIMOLD) [ "Calculations for statistical estimator" NEVENT(MYRANK)=NEVENT(MYRANK)+1; "Counter of primary histories" DO I=1,$NZBIN [ "DEPTH-DOSE" DEPTHDD(MYRANK,I)=DEPTHDD(MYRANK,I)+DEPTHD(MYRANK,I); DEPTHDD2(MYRANK,I)=DEPTHDD2(MYRANK,I)+DEPTHD(MYRANK,I)**2; DEPTHD(MYRANK,I)=0; "Initialization"] DO J=1,5 [ DO K=1,$NXBIN [ "DOSE-PROFILE" PROFILEDD(MYRANK,J,K)=PROFILEDD(MYRANK,J,K)+PROFILED(MYRANK,J,K); PROFILEDD2(MYRANK,J,K)= PROFILEDD2(MYRANK,J,K)+PROFILED(MYRANK,J,K)**2; PROFILED(MYRANK,J,K)=0; "Initialization"]]] IPRIMOLD=IPRIM; "Previous primary history number"] DO I=1,NOFRECY [ "START OF PHASE SPACE RECYCLE LOOP" IRI=2; "IF(NCOUNT.LE.NWRITE.AND.ILINES.LE.NLINES) [ OUTPUT EI,XI,YI,ZI,UI,VI,WI, IQI,IRI,IDINC; (7G15.7,3I5); ILINES=ILINES+1;]" CALL SHOWER(IQI,EI,XI,YI,ZI,UI,VI,WI,IRI,WTI); NCOUNT=NCOUNT + 1; IXXEND=IXX; "LAST RANDOM NUMBER USED" ] "END OF PHASE SPACE RECYCLE LOOP" ] "END OF SHOWER CALL LOOP" IF(IPRIMOLD.EQ.IPRIM) [ DO I=1,$NZBIN [ "DEPTH-DOSE" DEPTHDD(MYRANK,I)=DEPTHDD(MYRANK,I)+DEPTHD(MYRANK,I); DEPTHDD2(MYRANK,I)=DEPTHDD2(MYRANK,I)+DEPTHD(MYRANK,I)**2;] DO J=1,5 [ DO K=1,$NXBIN [ "DOSE-PROFILE" PROFILEDD(MYRANK,J,K)=PROFILEDD(MYRANK,J,K)+PROFILED(MYRANK,J,K); PROFILEDD2(MYRANK,J,K)= PROFILEDD2(MYRANK,J,K)+PROFILED(MYRANK,J,K)**2;]]] ELSE [ NEVENT(MYRANK)=NEVENT(MYRANK)+1; "Counter of primary histories" DO I=1,$NZBIN [ "DEPTH-DOSE" DEPTHDD(MYRANK,I)=DEPTHDD(MYRANK,I)+DEPTHD(MYRANK,I); DEPTHDD2(MYRANK,I)=DEPTHDD2(MYRANK,I)+DEPTHD(MYRANK,I)**2;] DO J=1,5 [ DO K=1,$NXBIN [ "DOSE-PROFILE" PROFILEDD(MYRANK,J,K)=PROFILEDD(MYRANK,J,K)+PROFILED(MYRANK,J,K); PROFILEDD2(MYRANK,J,K)= PROFILEDD2(MYRANK,J,K)+PROFILED(MYRANK,J,K)**2;]]] OUTPUT MYRANK; ('Process# ',I3,' Finished'); "******************************************************************" "******************** STEP 8. OUTPUT OF RESULTS *******************" "******************************************************************" DO I=1,$NZBIN [ CALL MPI_REDUCE(DEPTHDD(MYRANK,I),DEPTHDDSUM(I),1,MPI_REAL,MPI_SUM, 0,MPI_COMM_WORLD,IERR); CALL MPI_REDUCE(DEPTHDD2(MYRANK,I),DEPTHDD2SUM(I),1,MPI_REAL,MPI_SUM, 0,MPI_COMM_WORLD,IERR);] DO J=1,5 [ DO K=1,$NXBIN [ CALL MPI_REDUCE(PROFILEDD(MYRANK,J,K),PROFILEDDSUM(J,K),1,MPI_REAL, MPI_SUM,0,MPI_COMM_WORLD,IERR); CALL MPI_REDUCE(PROFILEDD2(MYRANK,J,K),PROFILEDD2SUM(J,K),1,MPI_REAL, MPI_SUM,0,MPI_COMM_WORLD,IERR);]] CALL MPI_REDUCE(NEVENT(MYRANK),NEVENTSUM,1,MPI_INTEGER,MPI_SUM,0, MPI_COMM_WORLD,IERR); CALL MPI_REDUCE(NINCELECT(MYRANK),NINCELECTSUM,1,MPI_INTEGER,MPI_SUM,0, MPI_COMM_WORLD,IERR); "CALCULATE STATISTICAL UNCERTAINTY FOR EACH VOXEL" "DEPTH-DOSE" IF(MYRANK.EQ.0) [ WRITE(13,:FMT10:); :FMT10: FORMAT('Depth dose curve',/,'Depth(cm)',2X, 'Dose (x10^-8 cGy/particle)',2X,'Relative uncertainty (%)'); DO I=1,$NZBIN [ "statistical uncertainty" SUDEPTHD(I)=SQRT((1.0/(NEVENTSUM-1))* ((DEPTHDD2SUM(I)/NEVENTSUM)-(DEPTHDDSUM(I)/NEVENTSUM)**2)); "relative uncertainty" RELASUD(I)=100*SUDEPTHD(I)/(DEPTHDDSUM(I)/NEVENTSUM); IF(I.EQ.1) [DEPTH=0;] ELSE [DEPTH=$VOXELZ*(I-1);] WRITE(13,:FMT11:)DEPTH,DEPTHDDSUM(I)/NINCELECTSUM/$NOFRECY,RELASUD(I); :FMT11: FORMAT(G10.3,7X,G12.5,8X,'+/-',3X,G11.5);] "DOSE-PROFILE" IDK1=1; IDK2=($NXBINOUT); IDK3=($NXBINOUT+1); IDK4=($NXBINOUT+$NXBININ); IDK5=($NXBINOUT+$NXBININ+1); IDK6=($NXBIN); DO J=1,5 [ IF(J.EQ.1) [ WRITE(13,:FMT12:)$DEPTH1; :FMT12: FORMAT(//,'Dose profile curve at',I3,'cm depth',/, 'Distance(cm)',1X,'Dose (x10^-8 cGy/particle)',2X, 'Relative uncertainty (%)');] ELSEIF(J.EQ.2) [ WRITE(13,:FMT13:)$DEPTH2; :FMT13: FORMAT(//,'Dose profile curve at',I3,'cm depth',/, 'Distance(cm)',1X,'Dose (x10^-8 cGy/particle)',2X, 'Relative uncertainty (%)');] ELSEIF(J.EQ.3) [ WRITE(13,:FMT14:)$DEPTH3; :FMT14: FORMAT(//,'Dose profile curve at',I3,'cm depth',/, 'Distance(cm)',1X,'Dose (x10^-8 cGy/particle)',2X, 'Relative uncertainty (%)');] ELSEIF(J.EQ.4) [ WRITE(13,:FMT15:)$DEPTH4; :FMT15: FORMAT(//,'Dose profile curve at',I3,'cm depth',/, 'Distance(cm)',1X,'Dose (x10^-8 cGy/particle)',2X, 'Relative uncertainty (%)');] ELSE [ WRITE(13,:FMT16:)$DEPTH5; :FMT16: FORMAT(//,'Dose profile curve at',I3,'cm depth',/, 'Distance(cm)',1X,'Dose (x10^-8 cGy/particle)',2X, 'Relative uncertainty (%)');] DO K=1,$NXBIN [ "statistical uncertainty" SUPROFILED(J,K)=SQRT((1.0/(NEVENTSUM-1))* ((PROFILEDD2SUM(J,K)/NEVENTSUM)-(PROFILEDDSUM(J,K)/NEVENTSUM)**2)); "relative uncertainty" RELASUP(J,K)=100*SUPROFILED(J,K)/(PROFILEDDSUM(J,K)/NEVENTSUM); IF(K.GE.IDK1.AND.K.LE.IDK2) [ OFFAXIS=(-$FIELDX*0.5+$VOXELXOUT*0.5)+$VOXELXOUT*(K-$NXBINOUT-1);] ELSEIF(K.GE.IDK3.AND.K.LE.IDK4) [ OFFAXIS=$VOXELXIN*(K-INT($NXBIN/2)-1);] ELSE [OFFAXIS=$FIELDX*0.5-$VOXELXOUT*0.5+$VOXELXOUT*(K-IDK5+1);] WRITE(13,:FMT17:)OFFAXIS,PROFILEDDSUM(J,K)/NINCELECTSUM/$NOFRECY, RELASUP(J,K); :FMT17: FORMAT(G11.5,6X,G12.5,8X,'+/-',3X,G11.5);]] ] "NEXT, CALL THE SUBROUTINE ECNSV1 TO WRITE-OUT THE ENERGY DEPOSITION" "TOTALS---TO CHECK ENERGY CONSERVATION FOR ONE THING" "CALL ECNSV1(1,NREG,TOTKE);" "CALL NTALLY(1,NREG);" CALL MPI_FINALIZE(IERR); STOP; END; "END OF MAIN PROGRAM" %E "******************************************************************" REAL*8 FUNCTION SPRNGCALL(stream); "******************************************************************" "SPRNG RANDOM NUMBER GENERATING SUBROUTINE" #include ; SPRNG_POINTER stream; [ SPRNGCALL=sprng(stream); ] RETURN; END; "END OF FUNCTION SPRNGCALL" %E "******************************************************************" " STANFORD LINEAR ACCELERATOR CENTER" SUBROUTINE AUSGAB(IARG); " EGS4 SUBPROGRAM - 8 MAY 1983/1730" "******************************************************************" ;COMIN/DEBUG,EPCONT,ETALY1,LINES,MISC,NTALY1,PASSIT,STACK,TOTALS,USEFUL, MPI/; REAL*8 DPWT; IRL=IR(NP); "SET LOCAL VARIABLE" DPWT=WT(NP); "KEEP TRACK OF THE ENERGY DEPOSITION---FOR CONSERVATION PURPOSES" ESUM(IQ(NP)+2,IRL,IARG+1)=ESUM(IQ(NP)+2,IRL,IARG+1)+EDEP*DPWT; NSUM(IQ(NP)+2,IRL,IARG+1)=NSUM(IQ(NP)+2,IRL,IARG+1) + 1; I=(IRL-3)/IXY+1; "SLAB NUMBER" J=(IRL-3-IXY*(I-1))/$NXBIN+1; "COLUMN NUMBER" K=IRL-2-IXY*(I-1)-$NXBIN*(J-1); "ROW NUMBER" "DEPTH-DOSE" IF(J.EQ.2.AND.K.EQ.INT($NXBIN/2+1)) [ IF(I.EQ.1) [DEPTHD(MYRANK,I)=DEPTHD(MYRANK,I)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXIN*$VOXELY*$VOXELZ*0.5);] ELSE [DEPTHD(MYRANK,I)=DEPTHD(MYRANK,I)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXIN*$VOXELY*$VOXELZ);]] "DOSE-PROFILE" IDK1=($NXBINOUT+1); IDK2=($NXBINOUT+$NXBININ); IF(I.EQ.INT($DEPTH1/$VOXELZ)+1.AND.J.EQ.2) [ IF(K.GE.IDK1.AND.K.LE.IDK2) [ PROFILED(MYRANK,1,K)=PROFILED(MYRANK,1,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXIN*$VOXELY*$VOXELZ);] ELSE [ PROFILED(MYRANK,1,K)=PROFILED(MYRANK,1,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXOUT*$VOXELY*$VOXELZ);]] ELSEIF(I.EQ.INT($DEPTH2/$VOXELZ)+1.AND.J.EQ.2) [ IF(K.GE.IDK1.AND.K.LE.IDK2) [ PROFILED(MYRANK,2,K)=PROFILED(MYRANK,2,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXIN*$VOXELY*$VOXELZ);] ELSE [ PROFILED(MYRANK,2,K)=PROFILED(MYRANK,2,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXOUT*$VOXELY*$VOXELZ);]] ELSEIF(I.EQ.INT($DEPTH3/$VOXELZ)+1.AND.J.EQ.2) [ IF(K.GE.IDK1.AND.K.LE.IDK2) [ PROFILED(MYRANK,3,K)=PROFILED(MYRANK,3,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXIN*$VOXELY*$VOXELZ);] ELSE [ PROFILED(MYRANK,3,K)=PROFILED(MYRANK,3,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXOUT*$VOXELY*$VOXELZ);]] ELSEIF(I.EQ.INT($DEPTH4/$VOXELZ)+1.AND.J.EQ.2) [ IF(K.GE.IDK1.AND.K.LE.IDK2) [ PROFILED(MYRANK,4,K)=PROFILED(MYRANK,4,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXIN*$VOXELY*$VOXELZ);] ELSE [ PROFILED(MYRANK,4,K)=PROFILED(MYRANK,4,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXOUT*$VOXELY*$VOXELZ);]] ELSEIF(I.EQ.INT($DEPTH5/$VOXELZ)+1.AND.J.EQ.2) [ IF(K.GE.IDK1.AND.K.LE.IDK2) [ PROFILED(MYRANK,5,K)=PROFILED(MYRANK,5,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXIN*$VOXELY*$VOXELZ);] ELSE [ PROFILED(MYRANK,5,K)=PROFILED(MYRANK,5,K)+EDEP*1.6*WT(NP)/ (RHOR(IRL)*$VOXELXOUT*$VOXELY*$VOXELZ);]] "IF(NCOUNT.LE.NWRITE.AND.ILINES.LE.NLINES) [ OUTPUT E(NP),X(NP),Y(NP),Z(NP),U(NP),V(NP),W(NP), IQ(NP),IRL,IARG; (7G15.7,3I5); ILINES=ILINES+1;]" RETURN; END; "END OF SUBROUTINE AUSGAB" %E "******************************************************************" " STANFORD LINEAR ACCELERATOR CENTER" SUBROUTINE HOWFAR; " EGS4 SUBPROGRAM - 8 MAY 1983/1730" "******************************************************************" ;COMIN/DEBUG,EPCONT,GEOM,LINES,PASSIT,STACK,THRESH/; IRL=IR(NP); "SET LOCAL VARIABLE" IF(IRL.EQ.1.OR.IRL.GE.NREG-4) [IDISC=1; RETURN;] ELSEIF(IRL.EQ.2) [ "Z-DIRECTION" $PLANE1(1,-1,IHIT,TPLN); IF(IHIT.EQ.1) [$CHGTR(TPLN,1);] $PLANE1(2,1,IHIT,TVAL); IF(IHIT.EQ.1) [ $FINVAL(TVAL,XX,YY,ZZ); XXX=XX+$PHANTOMX/2; YYY=YY+$PHANTOMY/2; IF(0.LE.XXX.AND.XXX.LT.$NXBINOUT*$VOXELXOUT) [IDX=XXX/$VOXELXOUT+1;] ELSEIF($NXBINOUT*$VOXELXOUT.LE.XXX.AND.XXX.LT.$NXBINOUT*$VOXELXOUT+ $NXBININ*$VOXELXIN) [ IDX=$NXBINOUT+((XXX-$NXBINOUT*$VOXELXOUT)/$VOXELXIN+1);] ELSE [ IDX=$NXBINOUT+$NXBININ+ ((XXX-$NXBINOUT*$VOXELXOUT-$NXBININ*$VOXELXIN)/$VOXELXOUT+1);] IF(0.LE.YYY.AND.YYY.LT.$PHANTOMY*0.5-$VOXELY*0.5) [IDY=1;] ELSEIF($PHANTOMY*0.5+$VOXELY*0.5.LE.YYY) [IDY=3;] ELSE [IDY=2;] $CHGTR(TVAL,2+(IDY-1)*$NXBIN+IDX);] "Y-DIRECTION" $PLAN2P(NXPLAN-1,NREG-3,1,NYPLAN,NREG-4,-1); "X-DIRECTION" $PLAN2P(NPLAN,NREG-1,1,NXPLAN,NREG-2,-1);] ELSE [ I=(IRL-3)/IXY+1; "SLAB NUMBER" J=(IRL-3-IXY*(I-1))/$NXBIN+1; "COLUMN NUMBER" K=IRL-2-IXY*(I-1)-$NXBIN*(J-1); "ROW NUMBER" "Z-DIRECTION" NPL1=I+2; NPL2=I+1; IF(I.LT.$NZBIN) [NRG1=IRL+IXY;] ELSE [NRG1=NREG;] IF(I.GT.1) [NRG2=IRL-IXY;] ELSE [NRG2=2;] $PLAN2P(NPL1,NRG1,1,NPL2,NRG2,-1); "Y-DIRECTION" NPL1=NYPLAN+J; NPL2=NPL1-1; IF(J.LT.$NYBIN) [NRG1=IRL+$NXBIN;] ELSE [NRG1=NREG-3;] IF(J.GT.1) [NRG2=IRL-$NXBIN;] ELSE [NRG2=NREG-4;] $PLAN2P(NPL1,NRG1,1,NPL2,NRG2,-1); "X-DIRECTION" NPL1=NXPLAN+K; NPL2=NPL1-1; IF(K.LT.$NXBIN) [NRG1=IRL+1;] ELSE [NRG1=NREG-1;] IF(K.GT.1) [NRG2=IRL-1;] ELSE [NRG2=NREG-2;] $PLAN2P(NPL1,NRG1,1,NPL2,NRG2,-1);] RETURN; END; "END OF SUBROUTINE HOWFAR" %E "******************************************************************" " STANFORD LINEAR ACCELERATOR CENTER" SUBROUTINE ECNSV1(NTREE,NREG,TOTKE); " EGS4 SUBPROGRAM - 8 MAY 1983/1730" "******************************************************************" " SUBROUTINE FOR KEEPING TRACK OF ENERGY CONSERVATION---TO BE " " USED WITH EGS USER CODES. WHEN NTREE=0, THE PROGRAM IS " " ENTERED IN ORDER TO INITIALIZE THE ESUM ARRAY TO ZERO. " " OTHERWISE, IT IS ENTERED FOR TOTALING AND OUTPUTTING THE " " RESULTS. THE ESUM ARRAY IS NEEDED IN SUBROUTINE AUSGAB, " " WHERE EDEP (ENERGY DEPOSITION) IS ADDED TO THE ELEMENT OF " " THE ARRAY CORRESPONDING TO THE VALUE OF IQ, IR, AND IARG. " "******************************************************************" COMIN/DEBUG,ETALY1/; "INSERT IN ALL SUBPROGRAMS THAT USE ESUM" REAL*8 ROWSUM(4,$MXREG),COLSUM(4,5),SUMSUM(4),GSUM,TOTKE; " CHECK WHETHER NREG IS GE $MXREG. IF IT IS, STOP AND OUTPUT. " IF(NREG.GT.$MXREG) [ MDUMMY=$MXREG; OUTPUT NREG,MDUMMY; (///,' ***** NOTE: STOPPED IN SUBROUTINE ECNSV1 BECAUSE NREG= ', I5,' IS LARGER THAN $MXREG= ',I5,' *****'); STOP;] IF(NTREE.EQ.0) [ "INITIALIZE ESUM TO ZERO AND RETURN" DO I=1,4 [DO J=1,NREG [DO K=1,5 [ESUM(I,J,K)=0.D0;]]] RETURN;] " REACH THIS POINT WHEN FINAL TALLY IS TO BE MADE." " FIRST, INITIALIZE SUMS" GSUM=0.D0; DO IQ=1,4 [ SUMSUM(IQ)=0.D0; DO IR=1,NREG [ROWSUM(IQ,IR)=0.D0;] DO ICODE=1,5 [COLSUM(IQ,ICODE)=0.D0;] " END OF IQ-LOOP"] " SUM IQ=1,2,3 INTO IQ=4 OF ESUM FOR ALL IR- AND ICODE-VALUES" DO IR=1,NREG [ DO ICODE=1,5 [ DO IQ=1,3 [ ESUM(4,IR,ICODE)=ESUM(4,IR,ICODE) + ESUM(IQ,IR,ICODE); ]]] " NORMALIZE DATA TO TOTKE" DO IQ=1,4 [ DO IR=1,NREG [ DO ICODE=1,5 [ ESUM(IQ,IR,ICODE)=ESUM(IQ,IR,ICODE)/TOTKE; ]]] " SUM-UP COLUMNS AND ROWS" DO IQ=1,4 [ DO IR=1,NREG [ DO ICODE=1,5 [ ROWSUM(IQ,IR)=ROWSUM(IQ,IR)+ESUM(IQ,IR,ICODE); ]] DO ICODE=1,5 [ DO IR=1,NREG [ COLSUM(IQ,ICODE)=COLSUM(IQ,ICODE)+ESUM(IQ,IR,ICODE); ]] " END OF IQ-LOOP"] " NOW GET TOTAL FOR IQ AND GRAND TOTAL" DO IQ=1,4 [ DO IR=1,NREG [ DO ICODE=1,5 [ SUMSUM(IQ)=SUMSUM(IQ)+ESUM(IQ,IR,ICODE); IF(IQ.LE.3) [GSUM=GSUM+ESUM(IQ,IR,ICODE);] ]]] " NOW WRITE-OUT THE RESULTS OF THE ENERGY DEPOSITION SUMMARY" DO IQ=1,4 [ IF(IQ.LE.3) [ IQNOW=IQ-2; OUTPUT IQNOW; ('1ENERGY DEPOSITION SUMMARY FOR PARTICLES WITH IQ=',I2,///, 55X,'IARG',/,19X,'0',15X,'1',13X,'2',14X,'3',14X,'4',16X,'ROW SUM', /,3X,'REGION',/); ] ELSE [ OUTPUT; ('1ENERGY DEPOSITION SUMMARY FOR ALL PARTICLES:',///, 55X,'IARG',/,19X,'0',15X,'1',13X,'2',14X,'3',14X,'4',16X,'ROW SUM', /,3X,'REGION',/); ] DO IR=1,NREG [ OUTPUT IR,(ESUM(IQ,IR,ICODE),ICODE=1,5),ROWSUM(IQ,IR); (I7,5X,5G15.7,5X,G15.7); " END OF IR-LOOP"] OUTPUT (COLSUM(IQ,ICODE),ICODE=1,5),SUMSUM(IQ); (/,3X,'COL SUM',2X,5G15.7,5X,G15.7); " END OF IQ-LOOP"] OUTPUT GSUM; (/////,' TOTAL FRACTION=',G15.7, ' NOTE: THIS NUMBER SHOULD BE VERY CLOSE TO UNITY'); RETURN; END; "END OF SUBROUTINE ECNSV1" %E "******************************************************************" " STANFORD LINEAR ACCELERATOR CENTER" SUBROUTINE NTALLY(NTREE,NREG); " EGS4 SUBPROGRAM - 8 MAY 1983/1730" "******************************************************************" " THIS PROGRAM KEEPS 'TALLY' OF THE NUMBER OF TIMES AUSGAB IS " " ENTERED FOR VARIOUS REASONS, ETC. IT GIVES THE USER A 'ROUGH' " " IDEA OF WHETHER CERTAIN EVENTS ARE RARE OR NOT. " " CAUTION *** DO NOT USE THESE NUMBERS IN ANY STATISTICAL SENSE] " "******************************************************************" COMIN/DEBUG,NTALY1/; INTEGER ROWSUM(4,$MXREG),COLSUM(4,5),SUMSUM(4),GSUM; " CHECK WHETHER NREG IS GE $MXREG. IF IT IS, STOP AND OUTPUT." IF(NREG.GT.$MXREG) [ MDUMMY=$MXREG; OUTPUT NREG,MDUMMY; (///,' ***** NOTE: STOPPED IN SUBROUTINE NTALLY BECAUSE NREG= ', I5,' IS LARGER THAN $MXREG= ',I5,' *****'); STOP;] IF(NTREE.EQ.0) [ "INITIALIZE NSUM TO ZERO AND RETURN" DO I=1,4 [DO J=1,NREG [DO K=1,5 [NSUM(I,J,K)=0;]]] RETURN;] " REACH THIS POINT WHEN FINAL TALLY IS TO BE MADE." " FIRST, INITIALIZE SUMS" GSUM=0; DO IQ=1,4 [ SUMSUM(IQ)=0; DO IR=1,NREG [ROWSUM(IQ,IR)=0;] DO ICODE=1,5 [COLSUM(IQ,ICODE)=0;] " END OF IQ-LOOP"] " SUM IQ=1,2,3 INTO IQ=4 OF NSUM FOR ALL IR- AND ICODE-VALUES" DO IR=1,NREG [ DO ICODE=1,5 [ DO IQ=1,3 [ NSUM(4,IR,ICODE)=NSUM(4,IR,ICODE) + NSUM(IQ,IR,ICODE); ]]] " SUM-UP COLUMNS AND ROWS" DO IQ=1,4 [ DO IR=1,NREG [ DO ICODE=1,5 [ ROWSUM(IQ,IR)=ROWSUM(IQ,IR)+NSUM(IQ,IR,ICODE); ]] DO ICODE=1,5 [ DO IR=1,NREG [ COLSUM(IQ,ICODE)=COLSUM(IQ,ICODE)+NSUM(IQ,IR,ICODE); ]] " END OF IQ-LOOP"] " NOW GET TOTAL FOR IQ AND GRAND TOTAL" DO IQ=1,4 [ DO IR=1,NREG [ DO ICODE=1,5 [ SUMSUM(IQ)=SUMSUM(IQ)+NSUM(IQ,IR,ICODE); IF(IQ.LE.3) [GSUM=GSUM+NSUM(IQ,IR,ICODE);] ]]] " NOW WRITE-OUT THE RESULTS" DO IQ=1,4 [ IF(IQ.LE.3) [ IQNOW=IQ-2; OUTPUT IQNOW; ('1SUMMARY OF EVENT COUNT FOR PARTICLES WITH IQ=',I2,///, 55X,'IARG',/,19X,'0',15X,'1',13X,'2',14X,'3',14X,'4',16X,'ROW SUM', /,3X,'REGION',/); ] ELSE [ OUTPUT; ('1SUMMARY OF EVENT COUNT FOR ALL PARTICLES:',///, 55X,'IARG',/,19X,'0',15X,'1',13X,'2',14X,'3',14X,'4',16X,'ROW SUM', /,3X,'REGION',/); ] DO IR=1,NREG [ OUTPUT IR,(NSUM(IQ,IR,ICODE),ICODE=1,5),ROWSUM(IQ,IR); (I7,5X,5I15,5X,I15); " END OF IR-LOOP"] OUTPUT (COLSUM(IQ,ICODE),ICODE=1,5),SUMSUM(IQ); (/,3X,'COL SUM',2X,5I15,5X,I15); " END OF IQ-LOOP"] OUTPUT GSUM; (/////,' TOTAL NUMBER OF EVENTS=',I15); RETURN; END; "END OF SUBROUTINE NTALLY" %E %C80 "-----------------------------------------------------------------" " Start of PRSTAAUX MORTRAN - Auxiliary codes required by PRESTA. " "-------------------------- Taken from NRCCAUXP MORTRAN. " " 24 August 1989 WRN " "-----------------------------------------------------------------" " " "*********************************************************************" " * * " " * FIXTMX * " " * * " " ********** " SUBROUTINE FIXTMX(ESTEP,MEDIUM); " " " THIS ROUTINE CHANGES THE STEP SIZE ALGORITHM USED IN EGS SO THAT " " THE STEP SIZE ARRAYS FOR TMXS CORRESPOND TO AN ARBITRARY,BUT " " FIXED FRACTIONAL ENERGY LOSS ESTEPE. " " IT IS ONLY NECESSARY FOR LOW ENERGY ELECTRON PROBLEMS SINCE " " TYPICALLY THE 200*TEFF0 RESTRICTION ON TMXS IS MORE STRINGENT " " FOR ELECTRONS WITH ENERGIES ABOVE A FEW MEV " " " " NOTE THAT THE $TMXS-OVER-RIDE MACRO MAY STILL BE IN FORCE IN EGS. " " " " THE ROUTINE CHANGES THE VALUES ONLY FOR THE MEDIUM 'MEDIUM' " " AND IT SHOULD PROBABLY BE USED FOR ALL MEDIA IN A PROBLEM. " ;" " " THE ROUTINE MUST BE CALLED AFTER HATCH HAS BEEN CALLED AND BEFORE " " THE SIMULATION IS BEGUN. " " " " THE ROUTINE IS INDEPENDENT OF WHAT UNITS ARE BEING USED, AS LONG " " AS THEY ARE CONSISTENT( E.G. CM, RL OR G/CM**2 ) " " " " IF CALLED WITH IOLDTM=0 (PASSSED IN COMIN USER) THE TMXS ARRAYS ARE" " ADJUSTED TO GIVE A FIXED ESTEPE AND ARE SUBJECTED TO THE TMIN AND " " CONSTRAINTS. " " IF CALLED WITH IOLDTM=1 THE CURRENT EGS ALGORITHM IS USED. " " IF CALLED WITH IOLDTM=0 AND ESTEPE=0 THE CURRENT EGS ALGORITHM IS " " USED. " " IF CALLED WITH IOLDTM=1 AND ESTEPE=0 THEN ESTEPE=1.0 IS USED. " " " " FOR A DETAILED DISCUSSION OF THE USE OF THIS ROUTINE, SEE " " 'Low Energy Electron Transport with EGS' in Nuclear Instr. and " " Methods A227 (1984)535-548. D.W.O. Rogers " " " " FOR A DISCUSSION OF THE NEW FEATURES (V03+) OF THIS ROUTINE, " " ESPECIALLY WITH REGARD TO THE NEW UPPER AND LOWER LIMITS, SEE " " 'PRESTA-the Parameter Reduced Electron-Step Transport Algorithm- " " for Electron Monte Carlo Transport' by A.F.Bielajew & D.W.O.Rogers," " NRCC Internal Report PIRS-042 obtainable by contacting the above. " " " " V01 DEC 10,1981 DAVE ROGERS NRCC " " V02 DEC 1984 EGS4 VERSION " " V03 JAN 1986 ALEX BIELAJEW NRCC REVISED FOR PRESTA " "*********************************************************************" ;COMIN/ELECIN,MEDIA,USER/; ESTEPE=ESTEP; IF(MEDIUM > $MXMED)["ERROR" OUTPUT MEDIUM; (///'0********* MEDIUM=',I4,' IN FIXTMX IS TOO LARGE');RETURN;] IF((ESTEPE = 0.) & (IOLDTM = 1)) RETURN; "USE THE CURRENT ALGORITHM " IF(ESTEPE = 0.) ESTEPE=1.; "NEW VERSION DEFAULTS TO TOTAL ENERGY LOSS" IF(IOLDTM = 0)[BLCCC=BLCC(MEDIUM);XCC2=XCC(MEDIUM)**2; "NEEDED BY ROOTMX"] "SET UP SOME VARIABLES FOR FIRST PASS THROUGH LOOP" EI =EXP( (1.-EKE0(MEDIUM))/EKE1(MEDIUM));"ENERGY OF FIRST TABLE ENTRY" EIL = ALOG(EI); LEIL=1; "THIS IS EQUIVALENT TO $SETINTERVAL EIL,EKE; BUT AVOIDS ROUNDOFF" $EVALUATE EDEDX USING EDEDX(EIL);"GET THE ELECTRON STOPPPING AT EI" "NOW CALCULATE STEP REQUIRED TO CAUSE AN ESTEPE REDUCTION IN ENERGY" IF(IOLDTM = 1)[SI=ESTEPE*EI/EDEDX;]ELSE[SI=ROOTMX(EI,ESTEPE);] "TABULATED ENERGIES ARE IN A FIXED RATIO - CALC LOG OF THE RATIO" ERATIO=-1./EKE1(MEDIUM); NEKE=MEKE(MEDIUM);"NUMBER OF ELEMENTS IN STORAGE ARRAY" DO I=1,NEKE-1[ EIP1=EXP((FLOAT(I+1)-EKE0(MEDIUM))/EKE1(MEDIUM));"ENERGY AT I+1" EIP1L=ALOG(EIP1);LEIP1L=I+1;"DESIGNED THIS WAY=$SETINTERVAL" $EVALUATE EDEDX USING EDEDX(EIP1L); IF(IOLDTM = 1)[SIP1=ESTEPE*EIP1/EDEDX;]ELSE[SIP1=ROOTMX(EIP1,ESTEPE);] "NOW SOLVE THESE EQUATIONS " " SI = TMXS1 * EIL + TMXS0 " " SIP1 = TMXS1 * EIP1L + TMXS0 " " " TMXS1(I,MEDIUM)=(SI-SIP1)/ERATIO;TMXS0(I,MEDIUM)=SI-TMXS1(I,MEDIUM)*EIL; "TRANSFER VALUES FOR NEXT LOOP" EIL=EIP1L;SI=SIP1;] "NOW PICK UP LAST TABLE ENTRY WHICH APPLIES ONLY TO LAST ENERGY" TMXS0(NEKE,MEDIUM)=TMXS0(NEKE-1,MEDIUM); TMXS1(NEKE,MEDIUM)=TMXS1(NEKE-1,MEDIUM); RETURN;END; %E "******************************************************************************" " * * " " * ROOTMX * " " * * " " ********** " " " FUNCTION ROOTMX(EI,ESTEP); " " " THIS ROUTINE RETURNS MAX(TMIN,MIN(TMAX,ESTEPE*EI/DEDX)) WHERE " " TMAX IS THE MAXIMUM STEP ALLOWED BY THE MOLIERE MULTIPLE SCATTERING " " THEORY, TMIN IS THE THE MINIMUN STEP AND ESTEPE*EI/DEDX IS THE GREATEST " " STEP ALLOWED DUE TO CONTINUOUS ENERGY LOSS PROCESSES. " " " " NOTE THE USE OF ITS AUXILLIARY FUNCTION FTMX APPENDED TO ROOTMX. " " BECAUSE THE TMAX FUNCTION IS STRONGLY ENERGY DEPENDENT, IT WAS FOUND " " NECESSARY TO INCLUDE A CORRECTION FOR ENERGY LOSS IN IT. OTHERWISE THE " " UPPER LIMIT COULD BE GREATLY EXCEEDED - BY AS MUCH AS 50% IN SOME CASES. " " CORRECTING FOR ENERGY LOSS NECESSITATES USING A ROOT FINDING METHOD TO " " OBTAIN TMAX (HENCE THE NAME ROOTMX). TMIN IS ALSO STRONGLY ENERGY " " DEPENDENT BUT IT DOES NOT MATTER WITHIN THE LOGIC OF THE CODE IF THIS " " QUANTITY IS AS MUCH AS 50% HIGH SINCE NO PHYSICS CONSTRAINTS WILL BE " " VIOLATED. " " " " THE ZERO-FINDING ROUTINE IS A CRUDE ONE BASED ON THE ASSUMPTION THAT " " THE FUNCTION FTMX IS MONOTONIC AND THAT THE FUNCTION EVALUATED AT THE TWO " " STARTING POINTS RETURNS DIFFERENT SIGNS. IF THE SIGNS ARE THE SAME THEN " " EITHER THE ENERGY-LOSS STEP-SIZE IS MORE RESTRICTIVE OR THE STEP-SIZE IS " " BELOW TMIN. " " " " ALTHOUGH THIS ROUTINE COMES WITH THE PRESTA PACKAGE IT IS REALLY " " INDEPENDENT OF IT AND IT IS AN IMPROVEMENT OVER THE PREVIOUS TMXS METHODS." " THE OLD TMXS ROUTINE ALLOWED BOTH THE TMAX AND TMIN BOUNDS TO BE VIOLATED." " EXCEEDING TMAX TAKES ONE OUT OF THE REGION OF VALIDITY OF THE MOLIERE " " THEORY STILL ALLOWING A MULTIPLE SCATTERING SELECTION BUT OF UNPREDICTABLE" " WORTH. GOING LOWER THAN TMIN CAUSES THE MULTIPLE SCATTERING TO GET " " SWITCHED OFF (STARTING WITH THE LOWER ENERGIES). THIS CAN SOMETIMES LEAD " " TO CALCULATIONAL ARTIFACTS. ONE WORD OF CAUTION] USING THIS ROUTINE AT " " VERY LOW ELECTRON ENERGIES .LE.10 keV CAUSES NEGATIVE USTEP ERRORS IF THE " " OLD EGS PATHLENGTH CORRECTION ALGORITHM (BASED ON FERMI-EYGES THEORY) IS " " USED. THE OLD EGS LESSENED THIS PROBLEM BY REDUCING THE UPPER LIMIT TO " " 0.8 THE VALUE USED IN THIS ROUTINE. THE PRESTA PATHLENGTH CORRECTION DOES " " NOT GIVE NEGATIVE USTEPS IN ANY OF THE CASES WE HAVE TESTED. " " " " VERSION 1 ALEX BIELAJEW JAN. 86 " " VERSION 1.1 ALEX BIELAJEW OCT. 87 " " Lower limit ESTEPE violation fixed " " " "******************************************************************************" ; COMIN/USEFUL,USER/; ESTEPE=ESTEP; TMIN=2.718282*EI*(EI+2.*RM)/(BLCCC*(EI+RM)**2); "LOWER LIMIT, eq.(2-8)" X1=TMIN; "INITIAL LOWER STARTING POINT OF THE SEARCH" X2=ESTEPE*EI/EDEDX; "INITIAL UPPER STARTING POINT OF THE SEARCH" "THIS IS THE FIX-UP FOR THE MINIMUM STEP-SIZE" IF( X2 <= X1 ) [ROOTMX=X1;RETURN;] F1=FTMX(X1,EI);F2=FTMX(X2,EI); AF1=ABS(F1);AF2=ABS(F2); SF1=SIGN(1.,F1);SF2=SIGN(1.,F2); "FIRST CHECK TO SEE IF EITHER OF THE STARTING POINTS IS ALREADY GOOD ENOUGH." IF((AF1 <= $ROOTMX_PRECISION) | (AF2 <= $ROOTMX_PRECISION))[ IF(AF1 <= AF2)[ROOTMX=X1;]ELSE[ROOTMX=X2;]] "NOW CHECK TO SEE IF EITHER THE ENERGY LOSS IS MORE RESTRICTIVE THAN THE " "UPPER LIMIT TMAX (TRUE FOR HIGH ENERGIES) OR IF IT MORE RESTRICTIVE THAN " "TMIN (TRUE FOR LOW ENERGIES WITH A SMALL ENOUGH ESTEPE). " ELSEIF(SF1 = SF2)[ROOTMX=X2;] "OTHERWISE A SEARCH FOR TMAX MUST BE UNDERTAKEN." ELSE[ "ITERATE" ITI=0; "NUMBER OF ITERATIONS COUNTER" XL=X1; "LAST X FOUND" :SEARCH-ROOT:LOOP[ ITI=ITI+1; IF(ITI > 1000)[ "QUIT IF THIS HAPPENS" OUTPUT;(' SEARCH FOR TMXS ABORTED. TOO MANY ITERATIONS');STOP;] XT=(X1*F2-X2*F1)/(F2-F1); IF(XT = XL)[ROOTMX=XT;EXIT:SEARCH-ROOT:;"CONVERGENCE OBTAINED"] FT=FTMX(XT,EI);AFT=ABS(FT); IF(AFT <= $ROOTMX_PRECISION)[ROOTMX=XT;EXIT:SEARCH-ROOT:;"CONVERGENCE OBTAINED"] ELSE[ "RE-ITERATE" SFT=SIGN(1.,FT); IF(SFT = SF1)[X1=XT;F1=FT;AF1=AFT;SF1=SFT;]ELSE[X2=XT;F2=FT;AF2=AFT;SF2=SFT;] XL=XT; "UPDATE LAST X FOUND" ] ] "END OF SEARCH FOR ROOT LOOP" ] "END OF ITERATE ELSE" RETURN;END; FUNCTION FTMX(T,EI); "When t=tmax as defined in eq.(2-10) this function returns 0. It is used by " "FUNCTION ROOTMX in the search for tmax. " COMIN/USEFUL,USER/; "Energy dependent quantities are evaluated at the energy mid-point of the step." "See section IV of the report PIRS-042. " EK=AMAX1(0.0001,EI-0.5*EDEDX*T);E=EK+RM;BETA2=EK*(E+RM)/E**2; A=BLCCC/BETA2;G=XCC2/(E*BETA2)**2; FTMX=1./ALOG(A/G)-G*T; RETURN;END; %E "******************************************************************" SUBROUTINE RMARIN; "******************************************************************" COMIN/RANDOM/; IF((IXX.LE.0).OR.(IXX.GT.31328)) IXX=1802; "SETS MARSAGLIA DEFAULT" " BUG. In the following line the assignment previous to 90/09/18 " " was to IXX. This DID NOT upset the randomness of the sequence, " " just the initial starting point. BLIF 90/09/18. " IF((JXX.LE.0).OR.(JXX.GT.30081)) JXX=9373; "SETS MARSAGLIA DEFAULT" I = MOD(IXX/177,177) + 2; J = MOD(IXX, 177) + 2; K = MOD(JXX/169,178) + 1; L = MOD(JXX, 169) ; DO II=1,97[ S=0.0;T=0.5; DO JJ=1,24[ M=MOD(MOD(I*J,179)*K,179); I=J;J=K;K=M;L=MOD(53*L+1,169); IF(MOD(L*M,64).GE.32) S=S+T; T=0.5*T; ] URNDM(II)=S; ] CRNDM = 362436./16777216.; CDRNDM = 7654321./16777216.; CMRNDM = 16777213./16777216.; IXX = 97; JXX = 33; RETURN;END; "******************************************************************" "************************** UCNAI4P END ***************************" "******************************************************************" %N