C C CROSSEC C Copyright (C) 1995 Don Cromer C C This code is distributed under the terms and conditions of the C CCP4 Program Suite Licence Agreement as a CCP4 Application. C A copy of the CCP4 licence can be obtained by writing to the C CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK. C C************************************************************* C C C************************************************************** C C PROGRAM TO INTERPOLATE CROSS SECTIONS C AND COMPUTE ANOMALOUS SCATTERING FACTORS C C Originally from the now-retired: C DON T. CROMER C C Some changes made at Daresbury/York and CCP4-ised. Beware that C illegal bizzareness in calling sequences with array/scalar C confusion has been eliminated, but the code may be fragile as a C consequence. C C May 1997 Input keyworded - M. Winn C C C INPUT CARDS C C COMPULSORY KEYWORDS: C ATOM - atomic symbol C either NWAV .... - list of wavelengths C or CWAV - wavelengths centred C on and separated by C END or - end input and run C C OPTIONAL KEYWORDS: C NORD - interpolation method (default 2) C VERB - verbose output (default only final table) C C C You can have multiple NWAV or CWAV cards. The total list C of wavelengths is sorted. The following variables are passed C from s/r RCARDS to the main program: C C ATMNAM=atom symbol C NW=TOTAL NUMBER OF XRAY WAVELENGTHS (maximum 1000) C XR=sorted list of wavelengths C NORD=INTERPOLATION TYPE(SEE BELOW) C LVERBOSE=verbosity flag C C NORD=0 INTERPOLATION IS BY FITTING THREE CLOSEST POINTS C TO A QUADRATIC C NORD=1,2,3 ETC. USE AITKEN'S INTERPOLATION METHOD AND NORD C IS THE INTERPOLATION ORDER. NORD=2 IS PROBABLY THE BEST C VALUE TO USE C C CROSS SECTION FILE IS READ FROM FILE IS. THIS COULD BE C THE SAME AS THE INPUT FILE. C C OUTPUT IS ON FILE IO C C *********NOTE******** C C THREE CONTRIBUTIONS TO FPRIME ARE LISTED SEPARATELY C IN THE OUTPUT. C 1. P.E. PHOTOELECTRIC CONTRIBUTION C 2. ETERM TERM DEPENDENT ON THE TOTAL ENERGY OF C THE ATOM. CROMER AND LIBERMAN J. CHEM. PHYS. C 53,1891-1898(1970). C 3. JENSEN TERM DEPENDENT ON THE XRAY ENERGY. C JENSEN, PHYSICS LETTERS,74A,41-44(1979) C C THE CROSS SECTION FILE WILL HAVE A NUMBER OF ORBITALS C FOR EACH ATOM FROM ATOMIC NUMBER 3-98 C THE FIRST MX CARDS (MX=5 FOR THIS XSECTION FILE AND IS SET IN PROG C FOR EACH ORBITAL WILL HAVE CROSS SECTIONS AT MX VALUES OF ENERGY C FROM ABOUT 1 TO 80 KEV APPROXIMATELY EQUALLY SPACED IN LOG(ENERGY) C THE NEXT FIVE RECORDS WILL BE CROSS SECTIONS AT ENERGIES SELECTED C BY GAUSS INTEGRATION SCHEME. IF THE FUNCTION TYPE IS 0 (IF=0) C ( FUNCTION TYPE GIVEN IN CROMER AND LIBERMAN C J. CHEM. PHYS. 53,1891-1898(1970) ) C A SIXTH VALUE IS READ IN FOR AN ENERGY=1.001*(BINDING ENERGY). C IF THE XRAY ENERGY IS LESS THAN BINDING ENERGY, FUNCTION SIGMA3 C WILL BE USED (CROMER AND LIBERMAN ACTA CRYST. A37,267-268(1981)) C C FILES IN IO AND IS ARE SET AT THE BEGINNING OF THE PROGRAM C THESE WILL DEPEND ON LOCAL CONVENTIONS C C***********WARNING************* C C IF AN XRAY ENERGY IS VERY CLOSE TO ONE OF THE ENERGIES C USED IN THE GAUSS INTEGRATION AN 'ANOMALOUS' ANOMALOUS SCATTERING C FACTOR MAY RESULT. THERE IS NO EASY WAY OUT OF THIS PROBLEM. A C SUGGESTED WAY IS TO COMPUTE SEVERAL VALUES AT NEARBY ENERGIES C AND DRAW A SMOOTH CURVE. THIS METHOD SHOULD WORK PROVIDED C THE POINTS DO NOT PASS THROUGH AN EDGE C C C .. Parameters .. INTEGER NWAV PARAMETER (NWAV=1000) C .. C .. Scalars in Common .. DOUBLE PRECISION BB, CX, RX, SEDGE INTEGER ICOUNT, NPLOTL CHARACTER*4 NAT LOGICAL LVERBOSE C .. C .. Arrays in Common .. DOUBLE PRECISION SIGG(5) INTEGER NSHEL(2, 24) C .. C .. Local Scalars .. DOUBLE PRECISION AU, BE, C, C1, CORR, EDG, PI, SIGEDG, XKEV, XW, + ZX INTEGER I, IF, IFAIL, IN, IO, IS, IZ, J, K, LDUM, M, + MF, MM, MX, N1, NJ, NNW, NO, NORD, NW, NX CHARACTER*4 ATMNAM C .. C .. Local Arrays .. DOUBLE PRECISION ABS(NWAV, 24), CXB(NWAV), EG(5), EL(14), + ET(NWAV), ETERM(98), EW(14), FP(NWAV, 24), + FPP(NWAV, 24), SIG(14), SL(14), SUMFP(NWAV), + SUMFPP(NWAV), T(32), WT(98), XJENSN(NWAV), + XK(NWAV), XR(NWAV) INTEGER MGAUSS(6), NORB(98) CHARACTER*2 ATMTAB(98) C .. C .. External Functions .. DOUBLE PRECISION AKNINT, GAUSS, SIGMA0, SIGMA1, SIGMA2, SIGMA3 EXTERNAL AKNINT, GAUSS, SIGMA0, SIGMA1, SIGMA2, SIGMA3 C .. C .. External Subroutines .. EXTERNAL CCPDPN, CCPERR, CCPFYP, CCPRCS, CCPUPC, + LGNDR, SORT, XSECT C .. C .. Intrinsic Functions .. INTRINSIC DEXP, DLOG, FLOAT C .. C .. Common blocks .. COMMON /EDGE/SEDGE COMMON /GAUS/CX, BB, SIGG, RX, ICOUNT COMMON /PLT/NPLOTL, NSHEL COMMON /PLTCH/NAT COMMON /CINPUT/ ATMNAM COMMON /IINPUT/ NORD,NW COMMON /DINPUT/ XR COMMON /LINPUT/ LVERBOSE C .. C .. Data statements .. DATA ATMTAB/' ', ' ', 'LI', 'BE', 'B ', 'C ', 'N ', 'O ', 'F ', + 'NE', 'NA', 'MG', 'AL', 'SI', 'P ', 'S ', 'CL', 'AR', 'K ', + 'CA', 'SC', 'TI', 'V ', 'CR', 'MN', 'FE', 'CO', 'NI', 'CU', + 'ZN', 'GA', 'GE', 'AS', 'SE', 'BR', 'KR', 'RB', 'SR', 'Y ', + 'ZR', 'NB', 'MO', 'TC', 'RU', 'RH', 'PD', 'AG', 'CD', 'IN', + 'SN', 'SB', 'TE', 'I ', 'XE', 'CS', 'BA', 'LA', 'CE', 'PR', + 'ND', 'PM', 'SM', 'EU', 'GD', 'TB', 'DY', 'HO', 'ER', 'TM', + 'YB', 'LU', 'HF', 'TA', 'W ', 'RE', 'OS', 'IR', 'PT', 'AU', + 'HG', 'TL', 'PB', 'BI', 'PO', 'AT', 'RN', 'FR', 'RA', 'AC', + 'TH', 'PA', 'U ', 'NP', 'PU', 'AM', 'CM', 'BK', 'CF'/ DATA (ETERM(I), I=1, 66)/.0, .0, .001, .001, .002, .003, .005, + .007, .009, .011, .014, .018, .021, .026, .030, .035, .041, + .047, .053, .060, .068, .075, .084, .093, .102, .113, .123, + .135, .146, .159, .172, .186, .200, .215, .231, .247, .264, + .282, .300, .319, .338, .359, .380, .401, .424, .447, .471, + .496, .521, .547, .575, .602, .631, .660, .690, .721, .753, + .786, .819, .854, .889, .925, .962, 1.0, 1.039, 1.079/ DATA (ETERM(I), I=67, 98)/1.119, 1.161, 1.204, 1.248, 1.293, + 1.338, 1.385, 1.433, 1.482, 1.532, 1.583, 1.636, 1.689, + 1.743, 1.799, 1.856, 1.914, 1.973, 2.033, 2.095, 2.157, + 2.221, 2.287, 2.353, 2.421, 2.490, 2.561, 2.633, 2.707, + 2.782, 2.858, 2.936/ DATA (WT(I), I=1, 54)/1.0079, 4.00260, 6.941, 9.01218, 10.81, + 12.011, 14.0067, 15.9994, 18.9984, 20.179, 22.98977, 24.305, + 26.98154, 28.0855, 30.97376, 32.06, 35.453, 39.948, 39.0983, + 40.08, 44.9559, 47.88, 50.9415, 51.996, 54.9380, 55.847, + 58.9332, 58.69, 63.546, 65.38, 69.72, 72.59, 74.9216, 78.96, + 79.904, 83.80, 85.4678, 87.62, 88.9059, 91.22, 92.9064, + 95.94, 98., 101.07, 102.9055, 106.42, 107.868, 112.41, + 114.82, 118.69, 121.75, 127.60, 126.9045, 131.29/ DATA (WT(I), I=55, 98)/132.9054, 137.33, 138.9055, 140.12, + 140.9077, 144.24, 145., 150.36, 151.96, 157.25, 158.9254, + 162.50, 164.9304, 167.26, 168.9342, 173.04, 174.67, 178.49, + 180.9479, 183.85, 186.207, 190.2, 192.22, 195.08, 196.9655, + 200.59, 204.383, 207.2, 208.9804, 209., 210., 222., 223., + 226.0254, 227.0278, 232.0381, 231.0359, 238.0482, 237.0482, + 239., 243., 247., 247., 251./ DATA (NORB(I), I=1, 80)/0, 0, 2, 2, 3, 3, 4, 4, 4, 4, 4, 4, 5, 6, + 7, 7, 7, 7, 7, 7, 7, 7, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, + 9, 9, 9, 12, 12, 13, 13, 14, 14, 14, 14, 14, 14, 14, 14, 14, + 14, 14, 14, 14, 17, 17, 17, 18, 18, 18, 18, 18, 18, 19, 19, + 19, 19, 19, 19, 19, 19, 20, 20, 20, 21, 21, 21, 21, 21, 22/ DATA (NORB(I), I=81, 98)/22, 23, 23, 24, 24, 24, 24, 24, 24, 24, + 24, 24, 24, 24, 24, 24, 24, 24/ DATA C/137.0367/ DATA PI/3.14159265/ DATA AU/2.80022E+7/ DATA C1/.02721/ DATA XKEV/12.397639/ DATA MGAUSS/5, 0, 0, 0, 0, 0/ C .. MX = 5 C CALL CCPFYP IFAIL = 0 CALL CCPDPN(5, 'DATA', 'READONLY', 'F', LDUM, IFAIL) CALL CCPDPN(6, 'PRINTER', 'PRINTER', 'F', LDUM, IFAIL) CALL CCPRCS(6, 'CROSSEC', '$Date: 2005/09/06 11:25:11 $') WRITE (6, FMT=9000) 9000 FORMAT (/' $TEXT:Reference1: $$ comment $$ ', + ' "Calculation of Anomalous Scattering Factors ', + 'at Arbitrary Wavelengths",', / + 'Don T. Cromer. J. Applied Cryst., Vol. 16 (1983) 437-8' + , //' $$', /' $SUMMARY :Reference1: $$ Crossec: $$', + /' :TEXT:Reference1: $$') C********** SET INPUT, OUTPUT AND X SECTION FILES***** IN = 5 IO = 6 IS = 1 WRITE(6,9200) ' COMPULSORY KEYWORDS:', + ' ATOM - atomic symbol', + ' either NWAV .... ', + ' - list of wavelengths', + ' or CWAV ', + ' - wavelengths centred on ', + ' and separated by ', + ' END or - end input and run' 9200 FORMAT(/,/,A,/,4X,A,/,4X,2A,/,4X,2A,/,30X,A,/,4X,A) WRITE(6,9210) ' OPTIONAL KEYWORDS:', + ' NORD - interpolation method (default 2)', + ' VERB - verbose output ', + ' (default only final table)' 9210 FORMAT(/,A,/,4X,A,/,4X,2A,/,/) C********** CALL CCPDPN(IS, 'CROSSECDATA', 'READONLY', 'F', 0, 0) C---Read keyworded input. Returns ATNAM,NORD,NW,XR. CALL RCARDS C*****READ ATOMIC SYMBOL AND ATOMIC NUMBER DO 10 IZ = 1, 98 IF (ATMNAM(1:2).EQ.ATMTAB(IZ)) GO TO 20 10 CONTINUE WRITE (IO, FMT='(3a)') ' Error - Atomic symbol not recognised', + ' Is this correct? ', ATMNAM 20 CONTINUE NO = NORB(IZ) WRITE (IO, FMT='(/,2a,3x,i4,/,a,i4,/,a,i4)') + ' Atom symbol and number ', ATMNAM, IZ, ' NORD value: ', NORD, + ' Number of wave lengths to analysis ', NW NNW = 0 DO 30 I = 1, NW IF (XR(I).LE.0.0) GO TO 30 NNW = NNW + 1 30 CONTINUE NW = NNW C*****READ NW VALUES OF XRAY(ANGSTROMS) AT WHICH INTERPOLATION C***** TO BE DONE C C*****CHANGE ANGSTROMS TO KEV DO 40 K = 1, NW XK(K) = XKEV/XR(K) CXB(K) = 0.0 40 CONTINUE REWIND IS DO 60 I = 1, NW DO 50 J = 1, NO FP(I, J) = 0.0 FPP(I, J) = 0.0 50 CONTINUE 60 CONTINUE DO 210 J = 1, NO DO 80 K = 1, MX 70 CONTINUE READ (IS, FMT=9031) NAT, NJ, NSHEL(1, J), NSHEL(2, J), XW, + EW(K), SIG(K) C NAT=ATOMIC SYMBOL C NJ=ORBITAL SEQUENCE NUMBER C NSHEL=ORBITAL TYPE. 1S1/2 ETC. C XW=ENERGY IN ANGSTROMS C EW=ENERGY IN KEV C SIG=CROSS SECTION IN BARNS C BE=BINDING ENERGY(KEV) C IF=0,1,2 FUNCTION TYPE C*****READ MX ENERGIES AND X SECTIONS FOR ORBITAL J IF (NAT.NE.ATMNAM) GO TO 70 IF (NJ.NE.J) GO TO 320 80 CONTINUE DO 90 K = 1, 5 READ (IS, FMT=9030) NAT, NJ, NSHEL(1, J), NSHEL(2, J), XW, + EG(K), SIGG(K), BE, IF C*****READ 5 ENERGIES AND X SECTIONS FOR ORBITAL J C***** FOR THE GAUSS INTEGRATION POINTS. BINDING ENERGY C AND FUNCTION TYPE ALSO IF (NAT.NE.ATMNAM) GO TO 320 IF (NJ.NE.J) GO TO 320 IF (BE.LE.0.) GO TO 320 EW(K+MX) = EG(K) SIG(K+MX) = SIGG(K) 90 CONTINUE NX = MX + 5 IF (IF.NE.0) GO TO 100 NX = NX + 1 C***** IF=0 SO READ X SEXTION AT ENERGY=1.001*BINDING ENERGY READ (IS, FMT=9031) NAT, NJ, NSHEL(1, J), NSHEL(2, J), XW, + EW(NX), SIG(NX) IF (NAT.NE.ATMNAM) GO TO 320 EDG = EW(NX) SIGEDG = SIG(NX) 100 CONTINUE BB = BE/C1 DO 110 I = 1, 5 SIGG(I) = SIGG(I)/AU 110 CONTINUE C*****SORT ALL X SECTIONS CALL SORT(NX, EW, SIG) C*****SORT THE FIVE CROSS SECTIONS AT INTEGRATION POINTS CALL SORT(5, EG, SIGG) C*****CHANGE TO LOGS FOR INTERPOLATION DO 120 K = 1, NX EL(K) = 0. SL(K) = 0. C EL(K) = DLOG(EW(K)) C IF (SIG(K).EQ.0.0) GO TO 120 C SL(K) = DLOG(SIG(K)) C 120 CONTINUE DO 200 K = 1, NW MF = 0 C ZX = DLOG(XK(K)) C C*****ZX=LOG OF XRAY(KEV) ENERGY CX = 0. IF (BE.GT.XK(K)) GO TO 170 IF (NORD.NE.0) GO TO 130 CALL XSECT(ZX, EL, SL, CX, NX) GO TO 160 130 CONTINUE DO 140 M = 1, NX N1 = M IF (SL(M).NE.0.0) GO TO 150 140 CONTINUE 150 CONTINUE MM = NX - N1 + 1 CX = AKNINT(ZX, MM, NORD, EL(N1), SL(N1), T) C CX = DEXP(CX) C C*****CX IS THE INTERPOLATED X SECTION IN BARNS 160 CONTINUE CXB(K) = CXB(K) + CX C *****CXB IS SUM TO GET MU/RHO C*****CHANGE CX TO ATOMIC UNITS CX = CX/AU 170 CONTINUE ICOUNT = 6 C*****RX=XRAY ENERGY IN ATOMIC UNITS RX = XK(K)/C1 IF (IF.NE.0) GO TO 180 IF (BE.LT.XK(K)) GO TO 180 C*****SEDGE IS XSECTION IN ATOMIC UNITS AT ENERGY=1.001*BE SEDGE = SIGEDG/AU CX = 0.0 FP(K, J) = GAUSS(1, SIGMA3, MGAUSS, LGNDR)*C/ (2.0*PI**2) MF = 3 GO TO 190 180 CONTINUE IF (IF.EQ.0) FP(K, J) = GAUSS(1, SIGMA0, MGAUSS, LGNDR)*C/ + (2.0*PI**2) IF (IF.EQ.1) FP(K, J) = GAUSS(1, SIGMA1, MGAUSS, LGNDR)*C/ + (2.0*PI**2) IF (IF.EQ.2) FP(K, J) = GAUSS(1, SIGMA2, MGAUSS, LGNDR)*C/ + (2.0*PI**2) 190 CONTINUE FPP(K, J) = 0. IF (CX.NE.0.0) FPP(K, J) = C*CX*RX/ (4.0*PI) CORR = 0.0 C IF (CX.NE.0.0) CORR = -CX*RX*0.5*DLOG((RX+BB)/ (RX-BB))*C/ + (2.0*PI**2) IF (MF.EQ.3) CORR = 0.5*SEDGE*BB**2*DLOG((-BB+RX)/ (-BB-RX))/ + RX*C/ (2.*PI**2) C FP(K, J) = FP(K, J) + CORR ABS(K, J) = CX*AU 200 CONTINUE C 210 CONTINUE IF (LVERBOSE) THEN WRITE (IO, FMT=9040) ATMNAM, IZ, NO, WT(IZ), ETERM(IZ) WRITE (IO, FMT=9010) NW, NORD WRITE (IO, FMT=9020) ATMNAM, (XR(I), I=1, NW) WRITE (IO, FMT=9050) ENDIF DO 220 I = 1, NW SUMFP(I) = 0.0 SUMFPP(I) = 0. 220 CONTINUE DO 240 J = 1, NO IF (LVERBOSE) WRITE (IO, FMT=9060) + NSHEL(1, J), NSHEL(2, J), (FP(K, J), K=1, NW) DO 230 K = 1, NW SUMFP(K) = SUMFP(K) + FP(K, J) 230 CONTINUE 240 CONTINUE IF (LVERBOSE) THEN WRITE (IO, FMT=9050) WRITE (IO, FMT=9090) (SUMFP(K), K=1, NW) ENDIF DO 250 K = 1, NW XJENSN(K) = -0.5*FLOAT(IZ)* (XK(K)/C1/137.0367**2)**2 250 CONTINUE DO 260 K = 1, NW ET(K) = -ETERM(IZ) 260 CONTINUE IF (LVERBOSE) THEN WRITE (IO, FMT=9100) (ET(K), K=1, NW) WRITE (IO, FMT=9080) (XJENSN(K), K=1, NW) ENDIF DO 270 K = 1, NW SUMFP(K) = SUMFP(K) + ET(K) + XJENSN(K) 270 CONTINUE IF (LVERBOSE) THEN WRITE (IO, FMT=9050) WRITE (IO, FMT=9070) (SUMFP(K), K=1, NW) WRITE (IO, FMT=9050) ENDIF DO 290 J = 1, NO IF (LVERBOSE) WRITE (IO, FMT=9060) + NSHEL(1, J), NSHEL(2, J), (FPP(K, J), K=1, NW) DO 280 K = 1, NW SUMFPP(K) = SUMFPP(K) + FPP(K, J) 280 CONTINUE 290 CONTINUE IF (LVERBOSE) THEN WRITE (IO, FMT=9050) WRITE (IO, FMT=9110) (SUMFPP(K), K=1, NW) ENDIF DO 300 K = 1, NW CXB(K) = CXB(K)*0.602472/WT(IZ) 300 CONTINUE IF (LVERBOSE) THEN WRITE (IO, FMT=9120) (CXB(K), K=1, NW) WRITE (IO, FMT=9050) WRITE (IO, FMT=9130) DO 310 J = 1, NO WRITE (IO, FMT=9140) + NSHEL(1, J), NSHEL(2, J), (ABS(K, J), K=1, NW) 310 CONTINUE ENDIF C Write xloggraph table: WRITE (IO, FMT='(/,3A,/,3A,/,A,/,A,/,A,/,A)') + ' $TABLE:Wave length v F'' and F"- ', ATMNAM, ' :', + ' $GRAPHS:Lambda v F'' and F" ', ATMNAM, ' :A:2,3,4: $$', + ' Atom_type Lambda F'' F" $$ ', ' Lambda F'' F" $$ ' WRITE (IO, FMT='(1000(A4,3x,3f11.4,/))') + (ATMNAM, XR(K), SUMFP(K), SUMFPP(K), K=1, NW) WRITE (IO, FMT='(A)') ' $$' C CLOSE(UNIT=IS,STATUS='KEEP') CALL CCPERR(0, 'Normal termination') 320 CONTINUE CLOSE(UNIT=IS,STATUS='KEEP') CALL CCPERR(1, 'Data file error') C 9010 FORMAT (' NW=', I5, ' NORD=', I5) 9020 FORMAT (' LAMBDA ', A4, 10F10.3, 9 (/12x, 10F10.4)) 9030 FORMAT (A4, I4, A4, A2, 1P, 3E15.8, E15.8, I2) 9031 FORMAT (A4, I4, A4, A2, 1P, 3E15.8) 9040 FORMAT (/,1X, A4, 2I5, 2F8.3) 9050 FORMAT (' ') 9060 FORMAT (2X, A4, A2, 10F10.3, 9 (/8x, 10F10.3)) 9070 FORMAT (' FP=', 10F10.3, 9 (/8x, 10F10.3)) 9080 FORMAT (' JENSEN=', 10F10.3, 9 (/8x, 10F10.3)) 9090 FORMAT (' P.E.', 10F10.3, 9 (/8x, 10F10.3)) 9100 FORMAT (' ETERM=', 10F10.3, 9 (/8x, 10F10.3)) 9110 FORMAT (' FPP=', 10F10.3, 9 (/8x, 10F10.3)) 9120 FORMAT (' MU/RHO=', 10F10.1, 9 (/8x, 10F10.1)) 9130 FORMAT (15X, 'CROSS SECTION(BARNS) AT XRAY ENERGY', /) 9140 FORMAT (2X, A4, A2, 1P, 10E10.3, 9 (/8x, 1P, 10E10.3)) END SUBROUTINE RCARDS INTEGER NPARM,NWAV PARAMETER (NPARM=200,NWAV=1000) C CHARACTER KEY*4,LINE*600,CVALUE(NPARM)*4,ATMNAM*4 REAL FVALUE(NPARM) LOGICAL LEND,LOGDEBUG,LATMNAM,LVERBOSE INTEGER ITOK,NTOK,IBEG(NPARM),IEND(NPARM),ITYPE(NPARM), + IDEC(NPARM),NORD,NW DOUBLE PRECISION XR(NWAV),DUMMY(NWAV) COMMON /CINPUT/ ATMNAM COMMON /IINPUT/ NORD,NW COMMON /DINPUT/ XR COMMON /LINPUT/ LVERBOSE C---Defaults NORD = 2 LVERBOSE = .FALSE. IWAV = 0 LATMNAM = .FALSE. LOGDEBUG = .TRUE. LEND = .FALSE. 10 CONTINUE NTOK=NPARM LINE=' ' CALL PARSER (KEY,LINE,IBEG,IEND,ITYPE,FVALUE,CVALUE,IDEC,NTOK, + LEND,LOGDEBUG) IF (LEND) GO TO 1000 C IF (KEY .EQ. 'END') THEN GO TO 1000 C ELSE IF (KEY .EQ. 'VERB') THEN LVERBOSE = .TRUE. C ELSE IF (KEY.EQ.'ATOM') THEN IF (NTOK.LT.2) + CALL CCPERR(1,' Argument expected for ATOM keyword.') ATMNAM = CVALUE(2) CALL CCPUPC(ATMNAM) LATMNAM = .TRUE. C ELSE IF (KEY.EQ.'NORD') THEN IF (NTOK.LT.2) THEN CALL CCPERR(3,' No argument for NORD keyword. Default used.') GO TO 10 ENDIF NORD = NINT(FVALUE(2)) C ELSE IF (KEY.EQ.'NWAV') THEN NWV = NINT(FVALUE(2)) IF (IWAV+NWV.GT.NWAV) CALL CCPERR(1, + ' Too many wavelengths requested: change NWAV in program.') IF (NTOK.LT.NWV+2) CALL CCPERR(1, + ' Too few arguments for NWAV. Some wavelengths missing?') DO ITOK = 3,NWV+2 IWAV = IWAV + 1 XR(IWAV) = FVALUE(ITOK) ENDDO C ELSE IF (KEY.EQ.'CWAV') THEN IF (NTOK.LT.4) + CALL CCPERR(1,' Not enough arguments for CWAV keyword.') NWV = NINT(FVALUE(2)) IF (IWAV+NWV.GT.NWAV) CALL CCPERR(1, + ' Too many wavelengths requested: change NWAV in program.') CENTRE = FVALUE(3) STEP = FVALUE(4) DO I = - NWV/2, (NWV - 1)/2 IWAV = IWAV + 1 XR(IWAV) = CENTRE + REAL(I)*STEP ENDDO ELSE C C----- ????? not understood C WRITE(6,*) ' Keyword ',KEY,' not recognised' END IF C GO TO 10 C 1000 CONTINUE C----Were compulsory keywords given? IF (.NOT.LATMNAM) CALL CCPERR(1,' ATOM must be given!') IF (IWAV.EQ.0) CALL CCPERR(1,' Wavelengths must be given!') NW = IWAV C----Wavelengths can come from several NWAV and CWAV cards, C so must be sorted. Use simple insertion routine included C here - OK for this small problem. CALL SORT(NW,XR,DUMMY) RETURN END C SUBROUTINE XSECT(ZX, EL, SL, CX, NX) C C C*****FIND EL(K) CLOSEST TO ZX C .. Scalar Arguments .. DOUBLE PRECISION CX, ZX INTEGER NX C .. C .. Array Arguments .. DOUBLE PRECISION EL(14), SL(14) C .. C .. Local Scalars .. DOUBLE PRECISION A0, A1, A2, DET, ER, P INTEGER L, LL C .. C .. Intrinsic Functions .. INTRINSIC DABS, DEXP C .. ER = 1000000. DO 10 L = 1, NX C P = DABS(ZX-EL(L)) C IF (P.GT.ER) GO TO 10 ER = P LL = L 10 CONTINUE LL = LL - 1 IF (LL.EQ.0) LL = 1 IF (LL.EQ.12) LL = 11 IF (SL(LL).EQ.0.) LL = LL + 1 DET = EL(LL+2)**2* (EL(LL+1)-EL(LL)) + + EL(LL+1)**2* (EL(LL)-EL(LL+2)) + + EL(LL)**2* (EL(LL+2)-EL(LL+1)) A0 = (EL(LL)**2* (SL(LL+1)*EL(LL+2)-SL(LL+2)*EL(LL+1))+ + EL(LL+1)**2* (SL(LL+2)*EL(LL)-SL(LL)*EL(LL+2))+ + EL(LL+2)**2* (SL(LL)*EL(LL+1)-SL(LL+1)*EL(LL)))/DET A1 = (EL(LL)**2* (SL(LL+2)-SL(LL+1))+ + EL(LL+1)**2* (SL(LL)-SL(LL+2))+ + EL(LL+2)**2* (SL(LL+1)-SL(LL)))/DET A2 = (SL(LL)* (EL(LL+2)-EL(LL+1))+SL(LL+1)* (EL(LL)-EL(LL+2))+ + SL(LL+2)* (EL(LL+1)-EL(LL)))/DET C CX = DEXP(A0+A1*ZX+A2*ZX**2) C END C*********************************************************************** C DOUBLE PRECISION FUNCTION SIGMA0(K, X) C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER K C .. C .. Scalars in Common .. DOUBLE PRECISION BB, CX, RX INTEGER ICOUNT C .. C .. Arrays in Common .. DOUBLE PRECISION SG(5) C .. C .. Common blocks .. COMMON /GAUS/CX, BB, SG, RX, ICOUNT C .. ICOUNT = ICOUNT - 1 SIGMA0 = SG(ICOUNT)*BB**3/X**2/ (RX**2*X**2-BB**2) - + BB*CX*RX**2/ (RX**2*X**2-BB**2) END C*********************************************************************** C DOUBLE PRECISION FUNCTION SIGMA1(K, X) C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER K C .. C .. Scalars in Common .. DOUBLE PRECISION BB, CX, RX INTEGER ICOUNT C .. C .. Arrays in Common .. DOUBLE PRECISION SG(5) C .. C .. Intrinsic Functions .. INTRINSIC DSQRT C .. C .. Common blocks .. COMMON /GAUS/CX, BB, SG, RX, ICOUNT C .. ICOUNT = ICOUNT - 1 C SIGMA1 = 0.5*BB**3*SG(ICOUNT)/ (DSQRT(X)* (RX**2*X**2-BB**2*X)) C END C*********************************************************************** C DOUBLE PRECISION FUNCTION SIGMA2(K, X) C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER K C .. C .. Scalars in Common .. DOUBLE PRECISION BB, CX, RX INTEGER ICOUNT C .. C .. Arrays in Common .. DOUBLE PRECISION SG(5) C .. C .. Local Scalars .. DOUBLE PRECISION DENOM C .. C .. Common blocks .. COMMON /GAUS/CX, BB, SG, RX, ICOUNT C .. ICOUNT = ICOUNT - 1 DENOM = X**3*RX**2 - BB**2/X SIGMA2 = 2.0*BB*SG(ICOUNT)*BB**2/X**4/DENOM - + 2.0*BB*CX*RX**2/DENOM END C***************************************************************** C DOUBLE PRECISION FUNCTION SIGMA3(K, X) C C .. Scalar Arguments .. DOUBLE PRECISION X INTEGER K C .. C .. Scalars in Common .. DOUBLE PRECISION BB, CX, RX, SEDGE INTEGER ICOUNT C .. C .. Arrays in Common .. DOUBLE PRECISION SG(5) C .. C .. Common blocks .. COMMON /EDGE/SEDGE COMMON /GAUS/CX, BB, SG, RX, ICOUNT C .. ICOUNT = ICOUNT - 1 SIGMA3 = BB**3* (SG(ICOUNT)-SEDGE*X**2)/ + (X**2* (X**2*RX**2-BB**2)) END C*********************************************************************** SUBROUTINE LGNDR(M, K, AA, Z) C C C .. Scalar Arguments .. DOUBLE PRECISION AA, Z INTEGER K, M C .. C .. Local Scalars .. DOUBLE PRECISION T INTEGER I4, IA, IH, IP, IS, KK C .. C .. Local Arrays .. DOUBLE PRECISION A(68), X(62) C .. C .. Intrinsic Functions .. INTRINSIC DSIGN, MOD C .. C .. Data statements .. DATA X(1)/.06943184420297/ DATA X(2)/.33000947820757/ DATA X(3)/.04691007703067/ DATA X(4)/.23076534494716/ DATA X(5)/.03376524289992/ DATA X(6)/.16939530676687/ DATA X(7)/.38069040695840/ DATA X(8)/.02544604382862/ DATA X(9)/.12923440720030/ DATA X(10)/.29707742431130/ DATA X(11)/.01985507175123/ DATA X(12)/.10166676129319/ DATA X(13)/.23723379504184/ DATA X(14)/.40828267875217/ DATA X(15)/.01591988024619/ DATA X(16)/.08198444633668/ DATA X(17)/.19331428364971/ DATA X(18)/.33787328829809/ DATA X(19)/.01304673574141/ DATA X(20)/.06746831665551/ DATA X(21)/.16029521585049/ DATA X(22)/.28330230293537/ DATA X(23)/.42556283050918/ DATA X(24)/.01088567092697/ DATA X(25)/.05646870011595/ DATA X(26)/.13492399721298/ DATA X(27)/.24045193539659/ DATA X(28)/.36522842202382/ DATA X(29)/.00921968287664/ DATA X(30)/.04794137181476/ DATA X(31)/.11504866290285/ DATA X(32)/.20634102285669/ DATA X(33)/.31608425050091/ DATA X(34)/.43738329574426/ DATA X(35)/.00790847264071/ DATA X(36)/.04120080038851/ DATA X(37)/.09921095463335/ DATA X(38)/.17882533027983/ DATA X(39)/.27575362448178/ DATA X(40)/.38477084202243/ DATA X(41)/.00685809565159/ DATA X(42)/.03578255816821/ DATA X(43)/.08639934246512/ DATA X(44)/.15635354759416/ DATA X(45)/.24237568182092/ DATA X(46)/.34044381553605/ DATA X(47)/.44597252564632/ DATA X(48)/.600374098758E-2/ DATA X(49)/.31363303799647E-1/ DATA X(50)/.75896708294787E-1/ DATA X(51)/.13779113431991/ DATA X(52)/.21451391369574/ DATA X(53)/.30292432646121/ DATA X(54)/.39940295300128/ DATA X(55)/.00529953250417/ DATA X(56)/.02771248846338/ DATA X(57)/.06718439880608/ DATA X(58)/.12229779582250/ DATA X(59)/.19106187779868/ DATA X(60)/.27099161117138/ DATA X(61)/.35919822461038/ DATA X(62)/.45249374508118/ DATA A(1)/.17392742256873/ DATA A(2)/.32607257743127/ DATA A(3)/.11846344252810/ DATA A(4)/.23931433524968/ DATA A(5)/.28444444444444/ DATA A(6)/.85662246189585E-1/ DATA A(7)/.18038078652407/ DATA A(8)/.23395696728635/ DATA A(9)/.06474248308443/ DATA A(10)/.13985269574464/ DATA A(11)/.19091502525256/ DATA A(12)/.20897959183674/ DATA A(13)/.05061426814519/ DATA A(14)/.11119051722669/ DATA A(15)/.15685332293894/ DATA A(16)/.18134189168918/ DATA A(17)/.04063719418079/ DATA A(18)/.09032408034743/ DATA A(19)/.13030534820147/ DATA A(20)/.15617353852000/ DATA A(21)/.16511967750063/ DATA A(23)/.07472567457529/ DATA A(24)/.10954318125799/ DATA A(25)/.13463335965500/ DATA A(26)/.14776211235738/ DATA A(27)/.02783428355809/ DATA A(28)/.06279018473245/ DATA A(29)/.09314510546387/ DATA A(30)/.11659688229599/ DATA A(31)/.13140227225512/ DATA A(32)/.13646254338895/ DATA A(33)/.02358766819326/ DATA A(34)/.05346966299766/ DATA A(35)/.08003916427167/ DATA A(36)/.10158371336153/ DATA A(37)/.11674626826918/ DATA A(38)/.12457352290670/ DATA A(39)/.02024200238266/ DATA A(40)/.04606074991886/ DATA A(41)/.06943675510989/ DATA A(42)/.08907299038097/ DATA A(43)/.10390802376845/ DATA A(44)/.11314159013145/ DATA A(45)/.11627577661544/ DATA A(46)/.01755973016588/ DATA A(47)/.04007904357988/ DATA A(48)/.06075928534395/ DATA A(49)/.07860158357910/ DATA A(50)/.09276919873897/ DATA A(51)/.10259923186065/ DATA A(52)/.10763192673158/ DATA A(53)/.01537662099806/ DATA A(54)/.03518302374405/ DATA A(55)/.05357961023359/ DATA A(56)/.06978533896308/ DATA A(57)/.08313460290850/ DATA A(58)/.09308050000778/ DATA A(59)/.09921574266356/ DATA A(60)/.10128912096278/ DATA A(61)/.01357622970588/ DATA A(62)/.03112676196932/ DATA A(63)/.04757925584125/ DATA A(64)/.06231448562777/ DATA A(65)/.07479799440829/ DATA A(66)/.08457825969750/ DATA A(67)/.09130170752246/ DATA A(68)/.09472530522754/ C .. KK = K IF ((M.GT.16) .OR. (M.LT.4)) KK = 4 IS = 0 IH = (M+1)/2 Z = .5 IF (MOD(M, 2).EQ.1) IS = -1 IP = KK T = 0. IF (IP.LE.IH) GO TO 10 IP = M + 1 - IP T = -1 10 CONTINUE I4 = M - 4 IA = (I4* (M+4)+IS)/4 + IP AA = A(IA) IF ((IP.EQ.IH) .AND. (IS.LT.0)) RETURN IA = IA - (I4+IS)/2 C Z = -T + DSIGN(X(IA), T) C END C*********************************************************************** C DOUBLE PRECISION FUNCTION GAUSS(N, Y, M, LTBL) C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. INTEGER M(6) C .. C .. Function Arguments .. DOUBLE PRECISION Y EXTERNAL Y C .. C .. Subroutine Arguments .. EXTERNAL LTBL C .. C .. Local Scalars .. INTEGER J, JJ, K, KK, L, LL, NN CHARACTER*50 CMT C .. C .. Local Arrays .. DOUBLE PRECISION A(6), G(6), Z(6) C .. C .. External Subroutines .. EXTERNAL LABRT C .. C .. Data statements .. DATA CMT/'GAUSS N NOT IN (1,6). RESULTS SET TO ZERO.'/ C .. C IF ((N.GT.6) .OR. (N.LT.1)) GO TO 140 NN = N GO TO (110, 90, 70, 50, 30, 10) NN 10 CONTINUE J = 1 G(6) = 0. 20 CONTINUE CALL LTBL(M(6), J, A(6), Z(6)) 30 CONTINUE K = 1 G(5) = 0. 40 CONTINUE CALL LTBL(M(5), K, A(5), Z(5)) 50 CONTINUE L = 1 G(4) = 0. 60 CONTINUE CALL LTBL(M(4), L, A(4), Z(4)) 70 CONTINUE JJ = 1 G(3) = 0. 80 CONTINUE CALL LTBL(M(3), JJ, A(3), Z(3)) 90 CONTINUE KK = 1 G(2) = 0. 100 CONTINUE CALL LTBL(M(2), KK, A(2), Z(2)) 110 CONTINUE LL = 1 G(1) = 0. 120 CONTINUE CALL LTBL(M(1), LL, A(1), Z(1)) G(1) = G(1) + A(1)*Y(1, Z(1)) LL = LL + 1 IF (LL.LE.M(1)) GO TO 120 IF (NN.EQ.1) GO TO 130 G(2) = G(2) + A(2)*Y(2, Z(1))*G(1) KK = KK + 1 IF (KK.LE.M(2)) GO TO 100 IF (NN.EQ.2) GO TO 130 G(3) = G(3) + A(3)*Y(3, Z(1))*G(2) JJ = JJ + 1 IF (JJ.LE.M(3)) GO TO 80 IF (NN.EQ.3) GO TO 130 G(4) = G(4) + A(4)*Y(4, Z(1))*G(3) L = L + 1 IF (L.LE.M(4)) GO TO 60 IF (NN.EQ.4) GO TO 130 G(5) = G(5) + A(5)*Y(5, Z(1))*G(4) K = K + 1 IF (K.LE.M(5)) GO TO 40 IF (NN.EQ.5) GO TO 130 G(6) = G(6) + A(6)*Y(6, Z(1))*G(5) J = J + 1 IF (J.LE.M(6)) GO TO 20 130 CONTINUE GAUSS = G(NN) RETURN 140 CONTINUE CALL LABRT(1, CMT, 1) GAUSS = 0. END C*********************************************************************** SUBROUTINE LABRT(ISW, LHOL, INX) C C C .. Scalar Arguments .. INTEGER INX, ISW CHARACTER*50 LHOL C .. C .. Local Scalars .. INTEGER NP LOGICAL PS, TS C .. C .. Data statements .. DATA NP/10/, PS/.TRUE./, TS/.FALSE./ C .. IF ((ISW.EQ.0) .OR. (ISW.GT.5)) RETURN GO TO (10, 20, 30, 40, 50) ISW 10 CONTINUE IF (PS .AND. (NP.GT.0)) PRINT 9000, LHOL, INX NP = NP - 1 CCC IF (TS) CALL EXIT RETURN 20 CONTINUE PS = .FALSE. RETURN 30 CONTINUE PS = .TRUE. NP = INX RETURN 40 CONTINUE TS = .TRUE. RETURN 50 CONTINUE TS = .FALSE. RETURN C 9000 FORMAT ('0', 9X, A50, 3X, I6) END C DOUBLE PRECISION FUNCTION AKNINT(XBAR, IN, IM, X, Y, T) C C AITKEN REPEATED INTERPOLATION C C XBAR = ABSCISSA AT WHICH INTERPOLATION IS DESIRED C IABS(IN) = NO. OF VALUES IN TABLE C IF IN.GT.0, CHK ORDERING OF X(I). C IF IN.LT.0, SKIP PRECEEDING TEST. C IM = DEGREE OF APPROXIMATING POLYNOMIAL C X = VECTOR OF IABS(IN) VALUES OF ABSCISSA C Y = VECTOR OF IABS(IN) VALUES OF ORDINATE C T = TEMPORARY STORAGE VECTOR OF 4*(M+1) LOCATIONS) C C C .. C .. Scalar Arguments .. DOUBLE PRECISION XBAR INTEGER IM, IN C .. C .. Array Arguments .. DOUBLE PRECISION T(80), X(14), Y(14) C .. C .. Local Scalars .. DOUBLE PRECISION DXBAR, S, Z INTEGER I, J, JJ, K, KK, M, MEND, N C .. C .. Intrinsic Functions .. INTRINSIC IABS, MAX0, MIN0 C .. DXBAR = XBAR N = IABS(IN) M = IM IF (M.GE.N) GO TO 130 10 CONTINUE K = N - 1 IF (N.LT.2) GO TO 120 S = X(2) - X(1) IF (IN.LT.0) GO TO 30 C CHK IF ORDER MONOTONIC IF (N.EQ.2) GO TO 30 DO 20 I = 3, N Z = (X(I)-X(I-1))*S IF (Z.LE.0.) GO TO 110 20 CONTINUE 30 CONTINUE IF (S.LT.0.) GO TO 50 C INCREASING ORDER DO 40 J = 1, N IF (XBAR.LE.X(J)) GO TO 70 40 CONTINUE J = N GO TO 70 C DECREASING ORDER 50 CONTINUE DO 60 J = 1, N IF (XBAR.GE.X(J)) GO TO 70 60 CONTINUE J = N 70 CONTINUE K = M M = M + 1 J = J - M/2 J = MAX0(J, 1) J = MIN0(J, N-K) MEND = J + K DO 80 I = J, MEND KK = I - J + 1 T(KK) = Y(I) T(KK+M) = X(I) - DXBAR 80 CONTINUE DO 100 I = 1, K KK = I + 1 DO 90 JJ = KK, M T(JJ) = (T(I)*T(JJ+M)-T(JJ)*T(I+M))/ (X(JJ+J-1)-X(I+J-1)) 90 CONTINUE 100 CONTINUE AKNINT = T(M) RETURN 110 CONTINUE PRINT 9000 9000 FORMAT (' AKNINT X(I) NOT SEQUENCED PROPERLY') 120 CONTINUE PRINT 9010 9010 FORMAT (' AKNINT N.LT.2 YBAR RETURNED AS Y(1)') AKNINT = Y(1) RETURN 130 CONTINUE PRINT 9020 9020 FORMAT (' AKNINT WARNING ORDER OF INTERPOLATION TOO LARGE') M = N - 1 GO TO 10 END C SUBROUTINE SORT(N, A, B) C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. DOUBLE PRECISION A(*), B(*) C .. C .. Local Scalars .. DOUBLE PRECISION X, Y INTEGER I, I1, J, M C .. M = N - 1 DO 20 I = 1, M I1 = I + 1 DO 10 J = I1, N IF (A(J).GT.A(I)) GO TO 10 X = A(J) Y = A(I) A(I) = X A(J) = Y X = B(J) Y = B(I) B(I) = X B(J) = Y 10 CONTINUE 20 CONTINUE END