PROGRAM XCOM1 C C XCOM1 (Version 1.0); Martin J. Berger, 16 Dec 1987. C Similar to XCOM, but with direct-accessd unformatted C database file UDAT. PARAMETER (ME=550,MEA=400,NPOL=1) DIMENSION E(108),AFIT(108),BFIT(108),CFIT(108),DFIT(108), 1 SCATCO(108),SCATIN(108),PHOT(108),PAIRAT(108),PAIREL(108), 2 SATCO(94),SATIN(94),POT(94),PDIF(94),PAIRT(94),PAIRL(94), 3 PR1(55),PR2(51),EPHOT(40),PHC(35,14),WEIGHT(100),NZ(100), 4 JM(100),JZ(100),KM(ME),KZ(ME),LM(14),LZ(14),JDG(ME),LEN(100), 5 SCTCO(ME),SCTIN(ME),PHT(ME),PHDIF(ME),ENB(80),ATWTS(100), 6 PRAT(ME),PREL(ME),AT(ME),ATNC(ME),EN(ME),ENL(ME),EAD(MEA), 7 IDG(14),IDGR(14),MAX(14),EDGEN(14),X(ME,8) CHARACTER XFIL1*30,XFIL2*30,SUB*72,ADG(14)*2,ALAB(ME)*2, 1 EDGE(94)*2,FF*1 DATA EPAIR1/1.022007E+06/,EPAIR2/2.044014E+06/ DATA NENG/80/,AVOG/0.60221367/ DATA ENB/1.0E+03,1.5E+03,2.0E+03,3.0E+03,4.0E+03,5.0E+03, 1 6.0E+03,8.0E+03,1.0E+04,1.5E+04,2.0E+04,3.0E+04,4.0E+04, 2 5.0E+04,6.0E+04,8.0E+04,1.0E+05,1.5E+05,2.0E+05,3.0E+05, 3 4.0E+05,5.0E+05,6.0E+05,8.0E+05,1.0E+06,1.022E+06,1.25E+06, 4 1.5E+06,2.0E+06,2.044E+06,3.0E+06,4.0E+06,5.0E+06,6.0E+06, 5 7.0E+06,8.0E+06,9.0E+06,1.0E+07,1.1E+07,1.2E+07,1.3E+07, 6 1.4E+07,1.5E+07,1.6E+07,1.8E+07,2.0E+07,2.2E+07,2.4E+07, 7 2.6E+07,2.8E+07,3.0E+07,4.0E+07,5.0E+07,6.0E+07,8.0E+07, 8 1.0E+08,1.5E+08,2.0E+08,3.0E+08,4.0E+08,5.0E+08,6.0E+08, 9 8.0E+08,1.0E+09,1.5E+09,2.0E+09,3.0E+09,4.0E+09,5.0E+09, 1 6.0E+09,8.0E+09,1.0E+10,1.5E+10,2.0E+10,3.0E+10,4.0E+10, 2 5.0E+10,6.0E+10,8.0E+10,1.0E+11/ DATA (ATWTS(KA),KA=1,72)/ 1 1.00794, 4.002602, 6.941, 9.012182, 2 10.811, 12.011, 14.00674, 15.9994, 3 18.9984032, 20.1797, 22.989768, 24.3050, 4 26.981539, 28.0855, 30.973762, 32.066, 5 35.4527, 39.948, 39.0983, 40.078, 6 44.955910, 47.88, 50.9415, 51.9961, 7 54.93805, 55.847, 58.93320, 58.69, 8 63.546, 65.39, 69.723, 72.61, 9 74.92159, 78.96, 79.904, 83.80, 1 85.4678, 87.62, 88.90585, 91.224, 2 92.90638, 95.94, 97.9072, 101.07, 3 102.9055, 106.42, 107.8682, 112.411, 4 114.82, 118.710, 121.75, 127.60, 5 126.90447, 131.29, 132.90543, 137.327, 6 138.9055, 140.115, 140.90765, 144.24, 7 144.9127, 150.36, 151.965, 157.25, 8 158.92534, 162.50, 164.93032, 167.26, 9 168.93421, 173.04, 174.967, 178.49/ DATA (ATWTS(KA),KA=73,100)/ 1 180.9479, 183.85, 186.207, 190.2, 2 192.22, 195.08, 196.96654, 200.59, 3 204.3833, 207.2, 208.98037, 208.9824, 4 209.9871, 222.0176, 223.0197, 226.0254, 5 227.0278, 232.0381, 231.03588, 238.0289, 6 237.0482, 239.0522, 243.0614, 247.0703, 7 247.0703, 251.0796, 252.083, 257.0951/ 5 FORMAT(A) 10 FORMAT(1H ) FF=CHAR(12) PRINT *,' Program XCOM1 (Version 1.1)' PRINT *,' M. J. Berger, 16 Dec 1987' PRINT 10 OPEN (UNIT=10,FILE='UDAT',FORM='UNFORMATTED', 1 ACCESS='DIRECT',RECL=4740) CALL SPEC1(SUB,NF,KMAX,NZ,WEIGHT,JENG,EAD,NEGO) GO TO (60,30,15),NEGO 15 NENG=JENG DO 20 N=1,NENG EN(N)=EAD(N) KZ(N)=-1 20 KM(N)=0 MMAX=1 LEN(1)=NENG GO TO 220 30 DO 40 N=1,NENG EN(N)=ENB(N) KZ(N)=0 40 KM(N)=N DO 50 J=1,JENG JZ(J)=-1 50 JM(J)=0 CALL MERGE(EN,KZ,KM,NENG,EAD,JZ,JM,JENG) GO TO 75 60 DO 70 N=1,NENG EN(N)=ENB(N) KZ(N)=0 70 KM(N)=N 75 DO 110 K=1,KMAX KR=NZ(K) READ (10,REC=KR) IZ,ATWT,MAXEDG,MAXE,IDG,ADG,EDGEN IF(MAXEDG)110,110,80 80 DO 90 I=1,MAXEDG LZ(I)=IZ 90 LM(I)=I+80 CALL MERGE(EN,KZ,KM,NENG,EDGEN,LZ,LM,MAXEDG) 110 CONTINUE IF(KZ(2))130,130,120 120 KOUNT=KOUNT+1 EAD(KOUNT)=SQRT(EN(1)*EN(2)) 130 DO 150 N=2,NENG IF(KZ(N))150,150,140 140 IF(KZ(N-1))150,150,141 141 KOUNT=KOUNT+1 IF(EN(N)-EN(N-1))143,142,143 142 EAD(KOUNT)=EN(N) EN(N-1)=EN(N-1)*0.99995 EN(N)=EN(N)*1.00005 GO TO 150 143 EAD(KOUNT)=SQRT(EN(N)*EN(N-1)) 150 CONTINUE IF(KOUNT)180,180,160 160 DO 170 J=1,KOUNT JZ(J)=-2 170 JM(J)=0 CALL MERGE(EN,KZ,KM,NENG,EAD,JZ,JM,KOUNT) 180 MX=0 DO 200 N=1,NENG IF(KZ(N))200,200,190 190 MX=MX+1 JDG(MX)=N 200 CONTINUE MMAX=MX+1 IF(MMAX-1)201,201,205 201 LEN(1)=NENG GO TO 220 205 MXED=MX LEN(1)=JDG(1) DO 210 M=2,MX 210 LEN(M)=JDG(M)-JDG(M-1)+1 LEN(MMAX)=NENG-JDG(MXED)+1 220 DO 230 N=1,NENG ENL(N)=LOG(EN(N)) SCTCO(N)=0.0 SCTIN(N)=0.0 PHT(N)=0.0 PRAT(N)=0.0 230 PREL(N)=0.0 DO 610 K=1,KMAX KR=NZ(K) READ (10,REC=KR) IZ,ATWT1,MAXEDG,MAXE,IDG,ADG,EDGEN,E, 1 SCATCO,SCATIN,PHOT,PAIRAT,PAIREL,MAXEP,PHC ATWT=ATWTS(KR) MAXK=0 IF(MAXEDG)250,250,240 240 MAXK=MAXE-IDG(MAXEDG)+1 GO TO 310 250 DO 300 M=1,MAXE SATCO(M)=SCATCO(M) SATIN(M)=SCATIN(M) POT(M)=PHOT(M) PDIF(M)=0.0 PAIRT(M)=PAIRAT(M) 300 PAIRL(M)=PAIREL(M) GO TO 395 310 MAX(1)=MAXEP DO 320 I=2,MAXEDG 320 MAX(I)=MAXEP-IDG(I-1)+I IRV=MAXEDG DO 325 I=1,MAXEDG IG=IDG(I) IP=I+80 SATCO(IP)=SCATCO(IG) SATIN(IP)=SCATIN(IG) POT(IP)=PHOT(IG) PDIF(IP)=PHOT(IG)-PHOT(IG-1) PAIRT(IP)=PAIRAT(IG) PAIRL(IP)=PAIREL(IG) EDGE(IP)=ADG(IRV) 325 IRV=IRV-1 MB=0 DO 340 M=1,MAXE DO 335 I=1,MAXEDG IF(M-IDG(I)+1)330,340,330 330 IF(M-IDG(I))335,340,335 335 CONTINUE MB=MB+1 SATCO(MB)=SCATCO(M) SATIN(MB)=SCATIN(M) POT(MB)=PHOT(M) PDIF(MB)=0.0 PAIRT(MB)=PAIRAT(M) PAIRL(MB)=PAIREL(M) 340 CONTINUE MS=0 DO 360 M=1,MAXE DO 350 I=1,MAXEDG IF(M-IDG(I)+1)350,360,350 350 CONTINUE MS=MS+1 E(MS)=E(M) SCATCO(MS)=SCATCO(M) SCATIN(MS)=SCATIN(M) PHOT(MS)=PHOT(M) PAIRAT(MS)=PAIRAT(M) PAIREL(MS)=PAIREL(M) 360 CONTINUE MAXE=MS 395 CALL REV(MAXE,E) CALL REV(MAXE,SCATCO) CALL REV(MAXE,SCATIN) CALL REV(MAXE,PHOT) CALL REV(MAXE,PAIRAT) CALL REV(MAXE,PAIREL) IF(MAXEDG)410,410,400 400 DO 405 M=1,MAXEP M1=M+MAXE-MAXEP 405 EPHOT(M)=E(M1) 410 E(51)=EPAIR2 E(55)=EPAIR1 DO 420 M=1,54 TERM=(E(M)-EPAIR1)/E(M) 420 PR1(M)=LOG(PAIRAT(M)/(TERM**3)) PR1(55)=3.006275*PR1(54)-2.577757*PR1(53)+0.571482*PR1(52) DO 425 M=1,50 TERM=(E(M)-EPAIR2)/E(M) 425 PR2(M)=LOG(PAIREL(M)/(TERM**3)) PR2(51)=3.006275*PR2(50)-2.577757*PR2(49)+0.571482*PR2(48) DO 430 M=1,MAXE E(M)=LOG(E(M)) SCATCO(M)=LOG(SCATCO(M)) SCATIN(M)=LOG(SCATIN(M)) 430 PHOT(M)=LOG(PHOT(M)) DO 435 M=1,MAXEP 435 EPHOT(M)=LOG(EPHOT(M)) GO TO (436,436,437), NF 436 FRAC=1.0 GO TO 438 437 FRAC=AVOG*WEIGHT(K)/ATWT 438 ECUT=EDGEN(MAXEDG) IF(NEGO.EQ.3) GO TO 460 IF(NENG-80)460,440,460 440 DO 450 N=1,NENG SCTCO(N)=SCTCO(N)+FRAC*SATCO(N) SCTIN(N)=SCTIN(N)+FRAC*SATIN(N) PHT(N)=PHT(N)+FRAC*POT(N) PRAT(N)=PRAT(N)+FRAC*PAIRT(N) 450 PREL(N)=PREL(N)+FRAC*PAIRL(N) GO TO 610 460 IMP=1 DO 600 N=1,NENG IF(KZ(N))500,490,470 470 IF(KZ(N)-IZ)500,480,500 480 NN=KM(N) PHDIF(N)=FRAC*PDIF(NN) ALAB(N)=EDGE(NN) 490 NN=KM(N) SCTCO(N)=SCTCO(N)+FRAC*SATCO(NN) SCTIN(N)=SCTIN(N)+FRAC*SATIN(NN) PHT(N)=PHT(N)+FRAC*POT(NN) PRAT(N)=PRAT(N)+FRAC*PAIRT(NN) PREL(N)=PREL(N)+FRAC*PAIRL(NN) GO TO 600 500 T=EN(N) TL=ENL(N) CALL SCOF(E,SCATCO,MAXE,AFIT,BFIT,CFIT,DFIT) CALL BSPOL(TL,E,AFIT,BFIT,CFIT,DFIT,MAXE,RES) SCTCO(N)=SCTCO(N)+FRAC*EXP(RES) CALL SCOF(E,SCATIN,MAXE,AFIT,BFIT,CFIT,DFIT) CALL BSPOL(TL,E,AFIT,BFIT,CFIT,DFIT,MAXE,RES) SCTIN(N)=SCTIN(N)+FRAC*EXP(RES) IF(T-EPAIR1)530,530,510 510 TERM=(T-EPAIR1)/T CALL SCOF(E,PR1,55,AFIT,BFIT,CFIT,DFIT) CALL BSPOL(TL,E,AFIT,BFIT,CFIT,DFIT,55,RES) PRAT(N)=PRAT(N)+FRAC*(TERM**3)*EXP(RES) IF(T-EPAIR2)530,530,520 520 TERM=(T-EPAIR2)/T CALL SCOF(E,PR2,51,AFIT,BFIT,CFIT,DFIT) CALL BSPOL(TL,E,AFIT,BFIT,CFIT,DFIT,51,RES) PREL(N)=PREL(N)+FRAC*(TERM**3)*EXP(RES) 530 IF(MAXEDG)540,540,550 540 CALL SCOF(E,PHOT,MAXE,AFIT,BFIT,CFIT,DFIT) CALL BSPOL(TL,E,AFIT,BFIT,CFIT,DFIT,MAXE,RES) PHT(N)=PHT(N)+FRAC*EXP(RES) GO TO 600 550 IF(T-ECUT)570,560,560 560 CALL SCOF(E,PHOT,MAXK,AFIT,BFIT,CFIT,DFIT) CALL BSPOL(TL,E,AFIT,BFIT,CFIT,DFIT,MAXK,RES) PHT(N)=PHT(N)+FRAC*EXP(RES) GO TO 600 570 IF(T-EDGEN(IMP))590,580,580 580 IMP=IMP+1 GO TO 570 590 MAXX=MAX(IMP) GO TO (592,594), NPOL 592 CALL BLIN(TL,EPHOT,PHC(1,IMP),MAXX,RES) GO TO 596 594 CALL SCOF(EPHOT,PHC(1,IMP),MAXX,AFIT,BFIT,CFIT,DFIT) CALL BSPOL(TL,EPHOT,AFIT,BFIT,CFIT,DFIT,MAXX,RES) 596 PHT(N)=PHT(N)+FRAC*EXP(RES) 600 CONTINUE 610 CONTINUE DO 622 N=1,NENG ATNC(N)=SCTIN(N)+PHT(N)+PRAT(N)+PREL(N) AT(N)=ATNC(N)+SCTCO(N) GO TO (622,621,622), NF 621 AT(N)=AVOG*AT(N)/ATWT ATNC(N)=AVOG*ATNC(N)/ATWT 622 CONTINUE PRINT *,' Specify file on which output (cross section table)' PRINT *,' is to be stored: ' READ 5, XFIL1 OPEN (UNIT=8,FILE=XFIL1) NN=1 JTAL=0 DO 880 N=1,NENG IF(JTAL)640,640,805 640 IF(N-1)650,650,645 645 WRITE (8,*) FF 650 WRITE (8,660) SUB 660 FORMAT(8X,A) WRITE (8,10) WRITE (8,662) 662 FORMAT(8X,'Constituents (Atomic Number:Fraction by Weight)') WRITE (8,665) (NZ(K),WEIGHT(K),K=1,KMAX) 665 FORMAT((8X,6(2X,I2,':',F7.5))) WRITE (8,10) GO TO (670,675,690), NF 670 WRITE (8,672) 672 FORMAT(8X,'Cross Sections') WRITE (8,685) GO TO 710 675 WRITE (8,680) 680 FORMAT(8X,'Cross Sections and Attenuation Coefficients') WRITE (8,685) 685 FORMAT(8X,'(Note that 1 b(arn) = 10**(-24) cm2)') GO TO 710 690 WRITE (8,700) 700 FORMAT(8X,'Partial Interaction Coefficients', 1' and Total Attenuation Coefficients') 710 WRITE (8,10) GO TO (711,715,715), NF 711 WRITE (8,712) 712 FORMAT(8X,'PHOTON',6X,'SCATTERING', 16X,'PHOTO-',4X,'PAIR PRODUCTION', 22X,'TOT.CROSS SECTION') GO TO 724 715 WRITE (8,720) 720 FORMAT(8X,'PHOTON',6X,'SCATTERING', 16X,'PHOTO-',4X,'PAIR PRODUCTION', 22X,'TOTAL ATTENUATION') 724 WRITE (8,730) 730 FORMAT(8X,'ENERGY',3X,'COHERENT INCOHER. ELECTRIC', 14X,'IN',7X,'IN',7X,'WITH',3X,'WITHOUT') WRITE (8,740) 740 FORMAT(34X,'ABSORPTION NUCLEAR ELECTRON COHERENT COHERENT') WRITE (8,750) 750 FORMAT(46X,'FIELD',4X,'FIELD',4X,'SCATT.',3X,'SCATT.') GO TO (791,780,760), NF 760 WRITE (8,770) 770 FORMAT(9X,'(MeV)',2X,7(2X,'(cm2/g)')) GO TO 800 780 WRITE (8,790) 790 FORMAT(9X,'(MeV)',3X,5('(b/atom) '),2X,'(cm2/g)',2X,'(cm2/g)') GO TO 800 791 WRITE (8,792) 792 FORMAT(9X,'(MeV)',3x,7('(b/atom) ')) 800 WRITE (8,10) 805 ENN=EN(N)*(1.0E-06) IF(KZ(N))870,870,810 810 PHT1=PHT(N)-PHDIF(N) GO TO (812,811,812), NF 811 AT1=AT(N)-0.6022045*PHDIF(N)/ATWT ATNC1=ATNC(N)-AVOG*PHDIF(N)/ATWT GO TO 815 812 AT1=AT(N)-PHDIF(N) ATNC1=ATNC(N)-PHDIF(N) 815 WRITE(8,820) ENN,SCTCO(N),SCTIN(N),PHT1, 1 PRAT(N),PREL(N),AT1,ATNC1 820 FORMAT(7X,1PE10.3,1P7E9.2) WRITE (8,10) X(NN,1)=ENN X(NN,2)=SCTCO(N) X(NN,3)=SCTIN(N) X(NN,4)=PHT1 X(NN,5)=PRAT(N) X(NN,6)=PREL(N) X(NN,7)=AT1 X(NN,8)=ATNC1 NN=NN+1 WRITE(8,840) KZ(N),ALAB(N),ENN,SCTCO(N),SCTIN(N),PHT(N), 1 PRAT(N),PREL(N),AT(N),ATNC(N) 840 FORMAT(1X,I2,1X,A2,1X,1PE10.3,1P7E9.2) X(NN,1)=ENN X(NN,2)=SCTCO(N) X(NN,3)=SCTIN(N) X(NN,4)=PHT(N) X(NN,5)=PRAT(N) X(NN,6)=PREL(N) X(NN,7)=AT(N) X(NN,8)=ATNC(N) NN=NN+1 JTAL=JTAL+3 850 IF(JTAL-43)880,880,860 860 JTAL=0 GO TO 880 870 WRITE (8,820) ENN,SCTCO(N),SCTIN(N),PHT(N),PRAT(N),PREL(N), 1 AT(N),ATNC(N) X(NN,1)=ENN X(NN,2)=SCTCO(N) X(NN,3)=SCTIN(N) X(NN,4)=PHT(N) X(NN,5)=PRAT(N) X(NN,6)=PREL(N) X(NN,7)=AT(N) X(NN,8)=ATNC(N) NN=NN+1 JTAL=JTAL+1 GO TO 850 880 CONTINUE CLOSE (8) PRINT 890, XFIL1 890 FORMAT(' Cross-section table with headings has been stored', 1' on file ',A) PRINT *,' ' PRINT *,' Options for further output:' PRINT *,' 1. No more output' PRINT *,' 2. Selected arrays stored on disk' PRINT *,' ' PRINT *,' Select option 1 or 2 --> ' READ *, NOUT GO TO (1020,920), NOUT 920 PRINT *,' ' PRINT *,' Specify file on which selected arrays are to be stored.' READ 5, XFIL2 OPEN (UNIT=9,FILE=XFIL2) WRITE (9,5) SUB WRITE (9,940) KMAX 940 FORMAT(12I6) WRITE (9,940) (NZ(K),K=1,KMAX) WRITE (9,950) (WEIGHT(K),K=1,KMAX) 950 FORMAT(6F9.6) WRITE (9,940) MMAX WRITE (9,940) (LEN(M),M=1,MMAX) 960 PRINT *,' Options for array to be stored on on disk:' PRINT *,' 0. Quit' PRINT *,' 1. Energy list' PRINT *,' 2. Coherent scattering cross section' PRINT *,' 3. Incoherent scattering cross section' PRINT *,' 4. Photoelectric absorption cross section' PRINT *,' 5. Pair prod. cross section (atomic nucleus)' PRINT *,' 6. Pair prod. cross section (atomic electrons)' PRINT *,' 7. Total attenuation coeff.(or cross section)' PRINT *,' 8. Attenuation coeff.(or cross section)' PRINT *,' without contribution from coherent scattering' PRINT *,' ' PRINT *,' Select option --> ' READ *, MAR IF(MAR)1000,1000,970 970 L1=1 GO TO (971,973,975,977,979,981,983,988), MAR 971 WRITE (9,972) 972 FORMAT(' ENERGY LIST') GO TO 993 973 WRITE (9,974) 974 FORMAT(' COHERENT SCATTERING CROSS SECTION') GO TO 993 975 WRITE (9,976) 976 FORMAT(' INCOHERENT SCATTERING CROSS SECTION') GO TO 993 977 WRITE (9,978) 978 FORMAT(' PHOTOELECTRIC ABSORPTION CROSS SECTION') GO TO 993 979 WRITE (9,980) 980 FORMAT(' PAIR PROD. CROSS SECTION (ATOMIC NUCLEUS)') GO TO 993 981 WRITE (9,982) 982 FORMAT(' PAIR PROD. CROSS SECTION (ATOMIC ELECTRONS)') GO TO 993 983 GO TO (984,986,986), NF 984 WRITE (9,985) 985 FORMAT(' TOTAL CROSS SECTION') GO TO 993 986 WRITE (9,987) 987 FORMAT(' TOTAL ATTENUATION COEFFICIENT') GO TO 993 988 GO TO (989,991,991), NF 989 WRITE (9,990) 990 FORMAT(' TOTAL CROSS SECTION WITHOUT COH. SCATTERING') GO TO 993 991 WRITE (9,992) 992 FORMAT(' TOTAL ATTENUATION COEFF. WITHOUT COH. SCATTERING') 993 DO 995 M=1,MMAX L2=L1+LEN(M)-1 WRITE (9,994) (X(L,MAR),L=L1,L2) 994 FORMAT(1P6E12.4) 995 L1=L2+1 PRINT *,' Make next selection or quit.' GO TO 960 1000 PRINT 1010, XFIL2 1010 FORMAT(' Selected cross section arrays have been stored', 1' on file ',A) 1020 PRINT *,' Calculation is finished.' STOP END SUBROUTINE SPEC1(SUBST,NFORM,MMAX,JZ,WT,JING,EADD,NELL) C SPEC1, 26 NOV 87. Revised version of SPEC, C to be used with XCOM1. PARAMETER (MEA=400) DIMENSION JZ(100),WT(100),JZ1(100),WT1(100),LH(100),WATE(100), 1 FRAC(100),EADD(MEA) CHARACTER*72 FORMLA,FRM(100),SUBST CHARACTER*30 ENGIN NFORM=3 10 FORMAT(1H ) 12 FORMAT(A) PRINT *,' Enter name of substance > ' READ 12, SUBST PRINT 10 PRINT *,' Options for characterization of substance:' PRINT *,' 1. Elemental substance, specified by atomic number' PRINT *,' 2. Elemental substance, specified by chemical symbol' PRINT *,' 3. Compound, specified by chemical formula' PRINT *,' 4. Mixture of elements and/or compounds' PRINT *,' ' PRINT *,' Enter desired option (1,2,3 or 4) --> ' READ *,NSUB PRINT 10 GO TO(15,20,50,55),NSUB 15 PRINT *,' Enter atomic number of element --> ' READ *, JZ(1) MMAX=1 WT(1)=1.0 GO TO 22 20 PRINT *,' Enter chemical symbol for element --> ' READ 12,FORMLA PRINT 10 CALL FORM(FORMLA,MMAX,JZ,WT) 22 PRINT *,' Options for output quantities:' PRINT *,' 1. Cross sections in barns/atom' PRINT *,' 2. Cross sections in barns/atom, and' PRINT *,' attenuation coefficients in cm2/g' PRINT *,' 3. Partial interaction coefficients and' PRINT *,' attenuation coefficients in cm2/g' PRINT *,' ' PRINT *,' Enter desired option (1, 2 or 3) --> ' READ *,NFORM PRINT 10 25 MAX=MMAX-1 DO 40 M=1,MAX MP1=M+1 DO 40 N=MP1,MMAX IF(JZ(M)-JZ(N))40,40,30 30 JZTEMP=JZ(M) WTTEMP=WT(M) JZ(M)=JZ(N) WT(M)=WT(N) JZ(N)=JZTEMP WT(N)=WTTEMP 40 CONTINUE GO TO 150 50 PRINT *,'Enter chemical formula for compound --> ' READ 12,FORMLA PRINT 10 CALL FORM(FORMLA,MMAX,JZ,WT) GO TO 25 55 PRINT *,' How many components in mixture? Enter number --> ' READ *,NCOMP PRINT 10 DO 65 N=1,NCOMP PRINT 60,N 60 FORMAT(' Enter chemical symbol or formula for component', 1 I3,'--> ') READ 12, FRM(N) PRINT 61,N 61 FORMAT(' Enter fraction by weight for component',I3,'--> ') READ *, FRAC(N) 65 CONTINUE PRINT 10 SUMF=0.0 DO 70 N=1,NCOMP 70 SUMF=SUMF+FRAC(N) PRINT 72 72 FORMAT(' Component Fraction') PRINT 73 73 FORMAT(' by Weight') PRINT 10 DO 75 N=1,NCOMP PRINT 74,N,FRAC(N),FRM(N) 74 FORMAT(I12,F12.6,6X,A) 75 CONTINUE PRINT 10 PRINT 80, SUMF 80 FORMAT(6X,'Sum = ',F12.6) PRINT 10 PRINT *,' Options for accepting or rejecting composition data:' PRINT *,' 1. Accept, but let program normalize fractions' PRINT *,' by weight so that their sum is unity' PRINT *,' 2. Reject, and enter different set of fractions' PRINT *,' ' PRINT *,' Enter desired option (1 or 2) --> ' READ *,MSUMGO PRINT 10 GO TO (85,55), MSUMGO 85 DO 90 N=1,NCOMP 90 FRAC(N)=FRAC(N)/SUMF DO 95 L=1,100 95 LH(L)=0 DO 120 N=1,NCOMP CALL FORM (FRM(N),MAX,JZ1,WT1) DO 120 M=1,MAX IN=JZ1(M) IF(LH(IN))100,100,110 100 LH(IN)=1 WATE(IN)=FRAC(N)*WT1(M) GO TO 120 110 WATE(IN)=WATE(IN)+FRAC(N)*WT1(M) 120 CONTINUE LL=0 DO 140 L=1,100 IF(LH(L))140,140,130 130 LL=LL+1 JZ(LL)=L WT(LL)=WATE(L) 140 CONTINUE MMAX=LL 150 PRINT *,' Options for energy list for output data:' PRINT *,' 1. Standard energy grid only' PRINT *,' 2. Standard grid plus additional energies' PRINT *,' 3. Additional energies only' PRINT *,' ' PRINT *,' Enter option desired (1,2 or 3) --> ' READ *,NELL PRINT 10 GO TO (240,190,190), NELL 190 PRINT *,' Modes of entering additional energies:' PRINT *,' 1. Entry from keyboard' PRINT *,' 2. Entry from prepared input file' PRINT *,' ' PRINT *,' Enter option desired (1 or 2) --> ' READ *, INEN GO TO (200,210), INEN 200 PRINT *,' How many additional energies are wanted? --> ' READ *,JING PRINT 10 DO 205 J=1,JING IF(J-1)201,201,202 201 PRINT *,' Enter first energy (in MeV) -->' GO TO 203 202 PRINT *,' Enter next energy (in MeV) --> ' 203 READ *, EADD(J) 205 CONTINUE PRINT 10 GO TO 220 210 PRINT *,' Specify file that contains input energy list.' PRINT *,' Specification can include drive and path.' PRINT *,' If necessary insert floppy disk with input data' PRINT *,' into selected drive before pressing "enter" --> ' READ 12, ENGIN PRINT 10 OPEN (UNIT=7,FILE=ENGIN) READ (7,*) JING,(EADD(J),J=1,JING) CLOSE (7) 220 DO 230 J=1,JING 230 EADD(J)=EADD(J)*1.0E+06 240 RETURN END SUBROUTINE FORM (W,MMAX,JZ,WT) C FORM, 26 Nov 87 DIMENSION MASH1(26),MASH2(418),IC(72),K(72),JZ(100),NZ(100), 1 MS(100),ATWTS(100),WT(100) CHARACTER*72,W DATA MASH1/0,5,6,0,0,9,0,1,53,0,19,0,0,7,8,15,0,0,16, 1 0,92,23,74,0,39,0/ DATA (MASH2(I),I=1,111)/70,36*0,54,3*0,73,8*0,65,8*0,43,51,88, 1 7*0,21,37,6*0,52,3*0,91,5*0,34,2*0,82,6*0,75,30,2*0,11,2*0, 2 90,3*0,46,0,41/ DATA (MASH2(I),I=112,205)/2*0,22,7*0,57,0,14,45,3*0,60,5*0,40, 1 2*0,10,2*0,81,8*0,69,9*0,62,5*0,12,2*0,50,2*0,31,0,28,4*0, 2 86,10*0,61,3*0,3,3*0,2,64,5*0,38/ DATA (MASH2(I),I=206,288)/0,72,84,3*0,20,3*0,80,0,26,3*0,56,6*0, 1 25,5*0,59,0,93,42,48,2*0,5 44,5*0,58,0,89,2*0,78,76,2*0,98, 2 4,3*0,94,6*0,49,15*0,36,47,0,67/ DATA (MASH2(I),I=289,418)/0,100,3*0,83,7*0,71,2*0,77,5*0,17,97, 1 7*0,96,10*0,13,3*0,87,2*0,27,0,95,4*0,68,8*0,99,10*0,24, 2 6*0,63,0,55,35,9*0,18,6*0,29,0,33,8*0,85,8*0,79,5*0,66/ DATA (ATWTS(KA),KA=1,72)/ 1 1.00794, 4.002602, 6.941, 9.012182, 2 10.811, 12.011, 14.00674, 15.9994, 3 18.9984032, 20.1797, 22.989768, 24.3050, 4 26.981539, 28.0855, 30.973762, 32.066, 5 35.4527, 39.948, 39.0983, 40.078, 6 44.955910, 47.88, 50.9415, 51.9961, 7 54.93805, 55.847, 58.93320, 58.69, 8 63.546, 65.39, 69.723, 72.61, 9 74.92159, 78.96, 79.904, 83.80, 1 85.4678, 87.62, 88.90585, 91.224, 2 92.90638, 95.94, 97.9072, 101.07, 3 102.9055, 106.42, 107.8682, 112.411, 4 114.82, 118.710, 121.75, 127.60, 5 126.90447, 131.29, 132.90543, 137.327, 6 138.9055, 140.115, 140.90765, 144.24, 7 144.9127, 150.36, 151.965, 157.25, 8 158.92534, 162.50, 164.93032, 167.26, 9 168.93421, 173.04, 174.967, 178.49/ DATA (ATWTS(KA),KA=73,100)/ 1 180.9479, 183.85, 186.207, 190.2, 2 192.22, 195.08, 196.96654, 200.59, 3 204.3833, 207.2, 208.98037, 208.9824, 4 209.9871, 222.0176, 223.0197, 226.0254, 5 227.0278, 232.0381, 231.03588, 238.0289, 6 237.0482, 239.0522, 243.0614, 247.0703, 7 247.0703, 251.0796, 252.083, 257.0951/ DO 116 L=1,72 IC(L)=ICHAR(W(L:L)) IF(IC(L)-32)101,102,103 101 K(L)=1 GO TO 116 102 K(L)=2 GO TO 116 103 IF(IC(L)-48)104,105,105 104 K(L)=1 GO TO 116 105 IF(IC(L)-58)106,107,107 106 K(L)=3 GO TO 116 107 IF(IC(L)-65)108,109,109 108 K(L)=1 GO TO 116 109 IF(IC(L)-91)110,111,111 110 K(L)=4 GO TO 116 111 IF(IC(L)-97)112,113,113 112 K(L)=1 GO TO 116 113 IF(IC(L)-123)114,115,115 114 K(L)=5 GO TO 116 115 K(L)=1 116 CONTINUE L=1 M=0 117 IF(K(L)-2)118,118,119 118 L=L+1 GO TO 117 119 LMIN=L 120 KG=K(L) IF(L-LMIN)130,130,140 130 GO TO (150,150,150,160,150), KG 140 GO TO (150,470,150,160,150), KG 150 STOP 1 160 KG1=K(L+1) GO TO (170,180,180,180,240), KG1 170 STOP 2 180 ICC=IC(L)-64 JT=MASH1(ICC) IF(JT)190,190,200 190 STOP 3 200 M=M+1 JZ(M)=JT GO TO (170,210,230,220,240), KG1 210 NZ(M)=1 GO TO 470 220 NZ(M)=1 L=L+1 GO TO 120 230 IN=L+1 GO TO 390 240 ICC=9*IC(L+1)-10*IC(L)+9 IF(ICC-1)310,250,250 250 IF(ICC-418)260,260,310 260 IF(ICC-208)300,270,300 270 M=M+1 IF(IC(L)-71)290,280,290 280 JZ(M)=32 GO TO 330 290 JZ(M)=84 GO TO 330 300 JT=MASH2(ICC) IF(JT)310,310,320 310 STOP 4 320 M=M+1 JZ(M)=JT 330 KG2=K(L+2) GO TO (340,350,380,360,370), KG2 340 STOP 5 350 NZ(M)=1 GO TO 470 360 NZ(M)=1 L=L+2 GO TO 120 370 STOP 6 380 IN=L+2 390 INN=IN IS=0 NZ(M)=0 400 IF(K(INN)-3)420,410,420 410 IS=IS+1 MS(IS)=IC(INN)-48 INN=INN+1 GO TO 400 420 ISM=IS KFAC=1 430 NZ(M)=NZ(M)+KFAC*MS(IS) KFAC=10*KFAC IS=IS-1 IF(IS)440,440,430 440 IF(NZ(M))450,450,460 450 STOP 7 460 L=IN+ISM GO TO 120 470 MMAX=M ASUM=0.0 DO 480 M=1,MMAX JM=JZ(M) 480 ASUM=ASUM+ATWTS(JM)*REAL(NZ(M)) DO 490 M=1,MMAX JM=JZ(M) 490 WT(M)=ATWTS(JM)*REAL(NZ(M))/ASUM RETURN END SUBROUTINE MERGE(E1,K1,L1,MMAX,E2,K2,L2,NMAX) DIMENSION E1(1),K1(1),L1(1),E2(1),K2(1),L2(1) DATA MLIM/200/ DO 50 N=1,NMAX M=2 10 IF(E2(N)-E1(M))30,30,20 20 M=M+1 GO TO 10 30 MC=M MF=M+1 MMAX=MMAX+1 IF(MMAX-MLIM)35,35,32 32 PRINT 33,MMAX,MLIM 33 FORMAT(6H MMAX=,I3,3X,6H MLIM=,I3) STOP 35 DO 40 M=MMAX,MF,-1 E1(M)=E1(M-1) K1(M)=K1(M-1) 40 L1(M)=L1(M-1) E1(MC)=E2(N) K1(MC)=K2(N) 50 L1(MC)=L2(N) RETURN END SUBROUTINE REV(NMAX,X) DIMENSION X(300) NH=NMAX/2 DO 10 N=1,NH N1=NMAX-N+1 T=X(N1) X(N1)=X(N) 10 X(N)=T RETURN END SUBROUTINE SCOF(X,F,NMAX,A,B,C,D) C SUBRROUTINE SCOF(X,F,NMAX,A,B,C,D), 22 FEB 83 C IF S LIES BETWEEN X(M) AND X(M+1), THEN C F(S)=((D(M)*S+C(M))*S+B(M))*S+A(M) DIMENSION X(1000),F(1000),A(1000),B(1000),C(1000),D(1000) M1=2 M2=NMAX-1 S=0.0 DO 10 M=1,M2 D(M)=X(M+1)-X(M) R=(F(M+1)-F(M))/D(M) C(M)=R-S 10 S=R S=0.0 R=0.0 C(1)=0.0 C(NMAX)=0.0 DO 20 M=M1,M2 C(M)=C(M)+R*C(M-1) B(M)=(X(M-1)-X(M+1))*2.0-R*S S=D(M) 20 R=S/B(M) MR=M2 DO 30 M=M1,M2 C(MR)=(D(MR)*C(MR+1)-C(MR))/B(MR) 30 MR=MR-1 DO 40 M=1,M2 S=D(M) R=C(M+1)-C(M) D(M)=R/S C(M)=C(M)*3.0 B(M)=(F(M+1)-F(M))/S-(C(M)+R)*S 40 A(M)=F(M) RETURN END SUBROUTINE BSPOL(S,X,A,B,C,D,N,G) C SUBROUTINE BSPOL(S,X,A,B,C,D,N,G), 22 FEB 83 DIMENSION X(1000),A(1000),B(1000),C(1000),D(1000) IF (X(1).GT.X(N)) GO TO 10 IDIR=0 MLB=0 MUB=N GO TO 20 10 IDIR=1 MLB=N MUB=0 20 IF (S.GE.X(MUB+IDIR)) GO TO 60 IF (S.LE.X(MLB+1-IDIR)) GO TO 70 ML=MLB MU=MUB GO TO 40 30 IF (IABS(MU-ML).LE.1) GO TO 80 40 MAV=(ML+MU)/2 IF (S.LT.X(MAV)) GO TO 50 ML=MAV GO TO 30 50 MU=MAV GO TO 30 60 MU=MUB+2*IDIR-1 GO TO 90 70 MU=MLB-2*IDIR+1 GO TO 90 80 MU=MU+IDIR-1 90 Q=S-X(MU) G=((D(MU)*Q+C(MU))*Q+B(MU))*Q+A(MU) RETURN END SUBROUTINE BLIN(S,X,Y,N,T) C SUBROUTINE BLIN, 12 APR 87 DIMENSION X(1000),Y(1000) IF (X(1).GT.X(N)) GO TO 10 IDIR=0 MLB=0 MUB=N GO TO 20 10 IDIR=1 MLB=N MUB=0 20 IF (S.GE.X(MUB+IDIR)) GO TO 60 IF (S.LE.X(MLB+1-IDIR)) GO TO 70 ML=MLB MU=MUB GO TO 40 30 IF (IABS(MU-ML).LE.1) GO TO 80 40 MAV=(ML+MU)/2 IF (S.LT.X(MAV)) GO TO 50 ML=MAV GO TO 30 50 MU=MAV GO TO 30 60 MU=MUB+2*IDIR-1 GO TO 90 70 MU=MLB-2*IDIR+1 GO TO 90 80 MU=MU+IDIR-1 90 Q=S-X(MU) T=Y(MU)+Q*(Y(MU+1)-Y(MU))/(X(MU+1)-X(MU)) RETURN END