PROGRAM DRIVER C C TITLE: CALCULATE GALERKIN COEFFICIENTS USING GNRT. C -------------------------------------------------- C C THIS IS A DRIVER PROGRAM FOR SUBROUTINE 'GNRT'. THE USER MUST C SUPPLY THE INPUT DATA IN A FILE, WHICH WILL BE DENOTED BY C 'NAMEIN' IN THIS DRIVER PROGRAM. THE ACTUAL NAME OF THE C INPUT FILE WILL BE READ BY THIS PROGRAM FROM THE STANDARD C INPUT UNIT, AND THE FILE NAME MUST BE A CHARACTER STRING. C THE INPUT PARAMETERS WILL THEN BE READ FROM FILE 'NAMEIN'. C C THE PARAMETERS THAT MUST BE SUPPLIED IN THE FILE 'NAMEIN' ARE C GIVEN BELOW, IN THE ORDER IN WHICH THEY MUST BE GIVEN IN THE FILE. C C NUMSUR THE INDEX OF THE SURFACE S IN SUBROUTINE 'SURFAC'. C A,B,C THE SURFACE PARAMETERS THAT REFER TO THE ELLIPSOID C PART IN THE DEFINITION OF S. SEE SUBROUTINE C 'SURFAC' FOR MORE DETAILS. C ALPHA AN AUXILIARY PARAMETER THAT MAY BE NEEDED IN C DEFINING THE SURFACE S. C DEGREE THE MAXIMUM DEGREE OF THE SPHERICAL HARMONICS TO C BE USED IN PRODUCING THE GALERKIN COEFFICIENTS. C THERE IS AN UPPER LIMIT ON 'DEGREE'. IT IS CALLED C 'MAXDEG', AND IT IS GIVEN BELOW IN THE PARAMETER C STATEMENT. C NINNER THE INTEGRATION PARAMETER USED IN APPROXIMATING THE C INNER SURFACE INTEGRAL FOR THE GALERKIN COEFFICIENTS. C NOUTER THE INTEGRATION PARAMETER USED IN EVALUATING THE C OUTER SURFACE INTEGRAL. THERE IS AN UPPER LIMIT C OF 20 FOR 'NOUTER', DUE TO THE PARTICULAR GAUSSIAN C QUADRATURE ROUTINE 'ZEROLG' BEING USED BY 'GNRT'. C NAMEDC THE NAME OF THE OUTPUT FILE THAT WILL CONTAIN THE C DECIMAL VERSION OF THE GALERKIN COEFFICIENTS. C NAMEFF THE NAME OF THE OUTPUT FILE THAT WILL CONTAIN THE C FORMATFREE VERSION OF THE GALERKIN COEFFICIENTS. C THEY ARE OUTPUT IN THIS FORM FOR MORE RAPID INPUT C TO THE PROGRAM 'LAPLAC' FOR SOLVING LAPLACE'S C EQUATION. C C THE FILE NUMBERS FOR THE INPUT AND OUTPUT FILES HAVE BEEN SET TO C DEFAULT VALUES, GIVEN IN A DATA STATEMENT BELOW. C C THE SUBROUTINE 'GNRT' REQUIRES A WORKING STORAGE ARRAY 'WORK' C TO BE SUPPLIED TO IT. THE DIMENSION OF THIS ARRAY MUST BE AT C LEAST THE MAXIMUM OF THE FOLLOWING QUANTITIES: C IT1=(DEGREE+1)**4 + 2*DEGREE**2 + 7*DEGREE +3 C + (3+4*DEGREE)*NOUTER + 7*NINNER C IT2=2*(DEGREE+1)**4 C THE DIMENSION OF WORK HAS BEEN DEFAULTED TO IT2, USING A MAXIMUM C DEGREE SPECIFIED IN THE PARAMETER STATEMENT BELOW. THIS WAS C CHOSEN BECAUSE IT2 IS USUALLY LARGER THAN IT1. IN ADDITION, C THERE IS A CHECK AT THE TIME OF INPUT AS TO WHETHER THE C DIMENSION OF 'WORK' IS LARGE ENOUGH. C C MACHINE DEPENDENT CONSTANTS ARE SET THROUGH THE FUNCTION C 'D1MACH', WHICH IS ONE OF THE SUBPROGRAMS GIVEN IN THE C 'GNRT' PACKAGE. THE ROUTINE 'D1MACH' SHOULD BE RESET C APPROPRIATELY FOR THE USER'S COMPUTER. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER DEGREE,FILEDC,FILEFF,FILEIN CHARACTER NAMEDC*50,NAMEFF*50,NAMEIN*50 PARAMETER(MAXDEG=10,IWORK=2*(MAXDEG+1)**4) DIMENSION WORK(IWORK) COMMON/SURPAR/A,B,C,ALPHA,NUMSUR COMMON/BLKWRK/WORK EXTERNAL SURFAC C C ****************************************** C THESE ARE THE INPUT/OUTPUT UNIT NUMBERS. C DATA FILEIN/8/,FILEDC/9/,FILEFF/10/ C C ****************************************** C C OBTAIN THE NAME OF THE INPUT FILE. PRINT *, ' GIVE THE NAME OF THE FILE CONTAINING THE' PRINT *, ' INPUT DATA.' READ(*,'(A)') NAMEIN C C OBTAIN INPUT DATA. L=LEN(NAMEIN) OPEN(FILEIN,FILE=NAMEIN(1:L)) READ(FILEIN,*)NUMSUR,A,B,C,ALPHA,DEGREE,NINNER,NOUTER READ(FILEIN,'(A)') NAMEDC,NAMEFF CLOSE(FILEIN) C C OPEN OUTPUT FILES AND PRINT INPUT PARAMETERS. L=LEN(NAMEDC) OPEN(FILEDC,FILE=NAMEDC(1:L)) L=LEN(NAMEFF) OPEN(FILEFF,FILE=NAMEFF(1:L),FORM='UNFORMATTED') WRITE(FILEDC,1000)NUMSUR,A,B,C,ALPHA,DEGREE,NINNER,NOUTER 1000 FORMAT('1SURFACE NUMBER ',I1,/,' SURFACE PARAMETERS (A,B,C)=', * 1PD17.10,2D21.10,/,' SURFACE PARAMETER ALPHA=',D17.10,/, * ' UPPER LIMIT ON DEGREE=',I2,/,' INTEGRATION PARAMETERS ', * 'NINNER, NOUTER = ',I2,I5,///) C C PRODUCE NEEDED SIZE FOR WORKING ARRAY 'WORK'. IT1=(DEGREE+1)**4+7*NINNER+(3+4*DEGREE)*NOUTER+2*DEGREE**2 * +7*DEGREE+3 IT2=2*(DEGREE+1)**4 NCHECK=MAX(IT1,IT2) C IF(NCHECK .LE. IWORK) THEN C THERE IS ENOUGH TEMPORARY STORAGE IN ARRAY 'WORK'. C CALL GNRT TO GENERATE THE GALERKIN COEFFICIENTS. CALL GNRT(SURFAC,DEGREE,NINNER,NOUTER,FILEDC,FILEFF,WORK) ELSE C THERE IS NOT ENOUGH WORKING STORAGE IN ARRAY 'WORK'. WRITE(FILEDC,1001) NCHECK 1001 FORMAT(' THERE IS NOT ENOUGH STORAGE IN THE ARRAY WORK.',/, * ' YOU NEED AT LEAST ',I5,' POSITIONS.') END IF CLOSE(FILEDC) CLOSE(FILEFF) STOP END SUBROUTINE SURFAC(Q,QCROSS,ST,CT,SP,CP,IFLAG) C ----------------- C C THIS CALCULATES A POINT Q ON THE SURFACE S, BASED ON S BEING THE C CONTINUOUS ONE-TO-ONE IMAGE OF THE UNIT SPHERE U. LET Q=(X,Y,Z) C DENOTE THE DESIRED POINT ON S. THEN IT IS TO CORRESPOND TO C THE POINT C UQ = (UX,UY,UZ) C ON U, WITH C UX = SIN(THETA)*COS(PHI) = ST*CP C UY = SIN(THETA)*SIN(PHI) = ST*SP C UZ = COS(THETA) = CT C C INPUT: C ST,CT,SP,CP THESE ARE THE GIVEN TRIGONOMETRIC FUNCTION VALUES C FOR THE POINT UQ. C IFLAG = 0 COMPUTE Q. C = 1 COMPUTE Q AND THE CROSS-PRODUCT C QCROSS = DPHI(Q) X DTHETA(Q)/SIN(THETA) C THE NOTATION DPHI(Q) MEANS THE DERIVATIVE OF Q REGARDED AS C AS A FUNCTION OF PHI, AND SIMILARLY FOR DTHETA(Q). THE QUANTITY C QCROSS IS USED IN COMPUTING THE NORMAL TO S AT Q; AND ITS C MAGNITUDE IS THE JACOBIAN IN THE CHANGE OF THE SURFACE C DIFFERENTIAL FOR CHANGING AN INTEGRAL OVER S TO AN INTEGRAL C OVER U. THE DIVISION BY SIN(THETA) REMOVES A SINGULARITY C IN THE MAPPING DUE TO THE USE OF SPHERICAL COORDINATES C IN REPRESENTING POINTS OF U. C OUTPUT: C Q THE POINT ON S CORRESPONDING TO UQ ON THE SPHERE U. C QCROSS THE CROSS-PRODUCT REFERRED TO ABOVE. C C THIS PARTICULAR SUBROUTINE ASSUMES THAT THE SURFACE S IS GIVEN BY C Q = (X,Y,Z) C = R*(A*UX,B*UY,C*UZ) C WITH R = R(UQ) A POSITIVE SMOOTH CONTINUOUS FUNCTION ON U. THE C CONSTANTS A,B,C ARE ASSUMED POSITIVE; AND FOR R=1, THE SURFACE C WILL BE AN ELLIPSOID WITH SEMI-AXES A,B,C. THIS GENERAL FORM C IS EQUIVALENT TO ASSUMING THAT THE REGION ENCLOSED BY S IS C STARLIKE WITH RESPECT TO THE ORIGIN. IT IS ALSO EQUIVALENT TO C THE STANDARD PARAMETRIC WAY OF DEFINING A SURFACE S IN TERMS C OF VARIABLES PHI AND THETA VARYING OVER (O,2*PI) AND (0,PI), C RESPECTIVELY. C C IT IS ASSUMED THAT R CAN BE EXTENDED AS A SMOOTH FUNCTION TO A C NEIGHBORHOOD OF U, SO THAT THE FIRST DERIVATIVES OF R(UX,UY,UZ) C CAN BE COMPUTED AT ALL POINTS OF U. THE USER MUST SUPPLY R AND C ITS DERIVATIVES WITH RESPECT TO UX,UY, AND UZ, DENOTED BY DXR, C DYR, AND DZR, RESPECTIVELY. EXAMPLES ARE GIVEN BELOW. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION Q(3),QCROSS(3) COMMON/MACHCN/U100,U100SQ C THE VARIABLE U100 IS 100 TIMES THE MACHINE UNIT ROUND. COMMON/SURPAR/A,B,C,ALPHA,NUMSUR C C COMPUTE UQ. UX = ST*CP UY = ST*SP UZ = CT C C ***** NOTE TO USER ***************************************************** C THIS BEGINS THE SECTION THAT THE USER MUST SPECIFY. C GO TO (10,20,30), NUMSUR C C SURFACE #1. C S IS AN ELLIPSOID. 10 R = 1.0D0 IF(IFLAG .EQ. 0) GO TO 100 DXR = 0.0D0 DYR = 0.0D0 DZR = 0.0D0 GO TO 100 C C SURFACE #2. C S IS A PEANUT-SHAPED REGION, BASED ON THE OVALS OF CASSINI. C THE PARAMETER ALPHA IS A NON-ZERO POSITIVE CONSTANT. AS IT C APPROACHES ZERO, THE SURFACE S BECOMES MORE CONSTRICTED IN C THE XY-PLANE. 20 C2T = 2.0D0*UZ*UZ - 1.0D0 TEMP = SQRT(ALPHA + C2T*C2T) R = SQRT(C2T + TEMP) IF(IFLAG .EQ. 0) GO TO 100 DXR = 0.0D0 DYR = 0.0D0 DZR=2.0D0*UZ*(1.0D0+C2T/TEMP)/R GO TO 100 C C SURFACE #3. C S IS A HEART-SHAPED REGION, INCREASINGLY SO AS THE POSITIVE C CONSTANT ALPHA INCREASES. 30 TEMP=1.0D0+ALPHA*(UZ+1.0D0)**2 R=2.0D0-1.0D0/TEMP IF(IFLAG .EQ. 0) GO TO 100 DXR=0.0D0 DYR=0.0D0 DZR=2.0D0*ALPHA*(UZ+1.0D0)/TEMP**2 GO TO 100 C C THIS ENDS THE SECTION THAT THE USER MUST SUPPLY. C ************************************************************************ C C COMPUTE THE POINT Q. 100 Q(1) = A*R*UX Q(2) = B*R*UY Q(3) = C*R*UZ IF(IFLAG .EQ. 0) RETURN C C COMPUTE CROSS-PRODUCT. C IF(ST .LE. U100) GO TO 200 C BEGIN BY COMPUTING DERIVATIVES OF Q WITH RESPECT TO PHI AND THETA. DPHIR = -UY*DXR + UX*DYR DTHER = CT*(CP*DXR+SP*DYR) - ST*DZR DPHIQX = A*ST*(-SP*R + CP*DPHIR) DTHEQX = A*CP*(CT*R + ST*DTHER) DPHIQY = B*ST*(CP*R + SP*DPHIR) DTHEQY = B*SP*(CT*R + ST*DTHER) DPHIQZ = C*CT*DPHIR DTHEQZ = C*(-ST*R + CT*DTHER) QCROSS(1) = (DPHIQY*DTHEQZ-DTHEQY*DPHIQZ)/ST QCROSS(2) = (DPHIQZ*DTHEQX-DTHEQZ*DPHIQX)/ST QCROSS(3) = (DPHIQX*DTHEQY-DTHEQX*DPHIQY)/ST RETURN C C COMPUTE CROSS-PRODUCT FOR Q AT A POLE. 200 QCROSS(1) = B*C*R*DXR QCROSS(2) = A*C*R*DYR QCROSS(3) = -A*B*UZ*R*R RETURN END SUBROUTINE GNRT(SURFAC,DEGREE,NINNER,NOUTER,FILEDC,FILEFF,WORK) C --------------- C C THIS ROUTINE CALCULATES THE GALERKIN COEFFICIENTS FOR SOLVING C THE DOUBLE LAYER POTENTIAL INTEGRAL EQUATION C (2*PI + K)*RHO = F . C IN THIS FORMULA, IT IS ASSUMED THAT THE VARIABLES HAVE BEEN C CHANGED SO THAT ALL FUNCTIONS ARE DEFINED ON THE UNIT SPHERE U. C THIS IS DONE AUTOMATICALLY IN SUBROUTINE 'GNRT', USING THE C SURFACE S AS SPECIFIED IN THE USER SUPPLIED SUBROUTINE 'SURFAC'. C C LET H(I) DENOTE THE SPHERICAL HARMONIC WITH INDEX I, BASED C ON THE ORDERING OF THEM AS DEFINED IN FUNCTION 'LOCATN' AND C SUBROUTINE 'INDX'. THEN 'GNRT' PRODUCES THE GALERKIN COEFFICIENTS C (H(I),K*H(J)) , C FOR 0 .LE. DEG(H(I)),DEG(H(J)) .LE. N, N='DEGREE' AN INPUT C VARIABLE. THESE ARE OUTPUT IN N TABLES. TABLE #L CONTAINS C THE GALERKIN COEFFICIENTS FOR ALL HARMONICS OF DEGREE .LE. L; C AND 1 .LE. L .LE. N. C C INPUT/OUTPUT C SURFAC: A USER SUPPLIED SUBROUTINE SPECIFYING THE SURFACE S. C SEE THE EXAMPLE DRIVER PROGRAM FOR THE SPECIFICATION OF C 'SURFAC'. C DEGREE: THE MAXIMUM DEGREE OF THE SPHERICAL HARMONICS TO BE C USED IN CALCULATING THE GALERKIN COEFFICIENTS. C NINNER: THE INNER INTEGRAL K*H(J) IS TO BE INTEGRATED BY THE C RULE SPECIFIED IN THE ACCOMPANYING ARTICLE, WITH THE C INTEGRATION PARAMETER 'NINNER'. C NOUTER: THE OUTER INTEGRAL, THE INNER PRODUCT (H(I),K*H(J)), C IS TO BE INTEGRATED WITH THE SPHERICAL GAUSSIAN QUADRATURE C FORMULA WITH INTEGRATION PARAMETER 'NOUTER'. THERE IS A C LIMIT OF 20 FOR 'NOUTER', DUE TO THE SUBROUTINE 'ZEROLG' C THAT IS BEING USED TO CALCULATE THE GAUSS-LEGENDRE NODES C AND WEIGHTS. USE ANOTHER SUCH ROUTINE WITH A LARGER UPPER C LIMIT TO REMOVE THE LIMIT ON 'NOUTER'. C FILEDC: THE FILE NUMBER FOR THE OUTPUT FILE CONTAINING THE C DECIMAL VERSION OF THE GALERKIN COEFFICIENTS. C FILEFF: THE FILE NUMBER FOR THE OUTPUT FILE CONTAINING THE C FORMATFREE VERSION OF THE GALERKIN COEFFICIENTS. C WORK: A TEMPORARY WORKING ARRAY FOR USE IN GNRT. ITS ORDER C MUST BE AT LEAST EQUAL TO THE MAXIMUM OF THE FOLLOWING C TWO QUANTITIES: C IT1 = (DEGREE+1)**4 + 2*DEGREE**2 + 7*DEGREE + 3 C + (3+4*DEGREE)*NOUTER + 7*NINNER C IT2 = 2*(DEGREE+1)**4 C IN THE PRESENT ROUTINE, THE ARRAY 'WORK' IS BROKEN APART C INTO SMALLER ONE AND TWO DIMENSIONAL ARRAYS, FOR USE IN THE C FOLLOWING SUBROUTINE 'GNRT2'. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER DEGREE,DEG2,DEGSQ,FILEDC,FILEFF DIMENSION WORK(*) COMMON/MACHCN/U100,U100SQ C THE VARIABLE U100 IS 100 TIMES THE MACHINE UNIT ROUND. EXTERNAL SURFAC C C SET MACHINE DEPENDENT CONSTANT. U100=100*D1MACH(4) U100SQ=U100**2 C C INITIALIZE VARIABLES TO BE USED IN GNRT2. DEG2=2*DEGREE+1 DEGSQ=(DEGREE+1)**2 NINNR2=2*NINNER NOUTR2=2*NOUTER LGDORD=1+(DEGREE*(DEGREE+3))/2 C C SET UP INITIAL ADDRESSES FOR BREAKING WORK INTO SMALLER C WORKING ARRAYS. I1=1 I2=I1+DEGSQ**2 I3=I2+NINNR2 I4=I3+NINNR2 I5=I4+DEGREE I6=I5+DEGREE I7=I6+NINNER I8=I7+NINNER I9=I8+NINNER I10=I9+LGDORD I11=I10+DEGREE*NOUTR2 I12=I11+DEGREE*NOUTR2 I13=I12+NOUTER I14=I13+NOUTER I15=I14+NOUTER I16=I15+LGDORD C C GO TO GNRT2 FOR ACTUAL COMPUTATION OF GALERKIN COEFFICIENTS. CALL GNRT2(SURFAC,DEGREE,NINNER,NOUTER,FILEDC,WORK(I1),WORK(I2), * WORK(I3),WORK(I4),WORK(I5),WORK(I6),WORK(I7),WORK(I8), * WORK(I9),WORK(I10),WORK(I11),WORK(I12),WORK(I13),WORK(I14), * WORK(I15),WORK(I16),NINNR2,NOUTR2,LGDORD,DEG2,DEGSQ) C C EXPAND AND WRITE THE FILE OF GALERKIN COEFFICIENTS. CALL EXPAND(WORK(I1),WORK(I2),DEGREE,DEGSQ,NINNER,NOUTER, * FILEFF) RETURN END SUBROUTINE GNRT2(SURFAC,DEGREE,NINNER,NOUTER,FILEDC, * KERMAT,SNPHI,CSPHI,SNPHI2,CSPHI2,SNTHI,CSTHI,WI,VALGDI, * SNPHE,CSPHE,SNTHE,CSTHE,WE,VALGDE,TEMPKH,NINNR2,NOUTR2, * LGDORD,DEG2,DEGSQ) C ---------------- C C THIS ROUTINE IS WHERE THE GALERKIN COEFFICIENTS ARE C ACTUALLY CALCULATED. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION KERMAT,KERNEL INTEGER DEGREE,DEG2,DEGSQ,FILEDC COMMON/BLOCKR/V,STI,CTI,SPI,CPI,UPX,UPY,UPZ DIMENSION KERMAT(DEGSQ,DEGSQ),SNPHI(NINNR2),CSPHI(NINNR2), * SNPHI2(DEGREE),CSPHI2(DEGREE),SNTHI(NINNER),CSTHI(NINNER), * WI(NINNER),VALGDI(LGDORD),SNPHE(DEGREE,NOUTR2), * CSPHE(DEGREE,NOUTR2),SNTHE(NOUTER),CSTHE(NOUTER),WE(NOUTER), * VALGDE(LGDORD),TEMPKH(DEGSQ) DIMENSION P(3),Q(3),QCROSS(3),V(3) PARAMETER(ZERO=0.0D0,P5=0.5D0,ONE=1.0D0,TWO=2.0D0, * FOUR=4.0D0) C C CALCULATE NODES,WEIGHTS, AND TRIG FUNCTIONS FOR INTERIOR AND C EXTERIOR INTEGRATIONS. PI=FOUR*ATAN(ONE) TWOPI=TWO*PI HI=PI/NINNER HE=PI/NOUTER C CALCULATE TRIG VALUES FOR INNER INTEGRATION. CALL INTMI(-ONE,ONE,NINNER,CSTHI,WI) DO 10 J=1,NINNER 10 SNTHI(J)=SIN(ACOS(CSTHI(J))) DO 20 L=1,NINNR2 PHI=(L-P5)*HI SNPHI(L)=SIN(PHI) 20 CSPHI(L)=COS(PHI) C C CALCULATE TRIG FUNCTIONS FOR EXTERIOR INTEGRATION. CALL ZEROLG(NOUTER,CSTHE,WE,-ONE,ONE) DO 30 J=1,NOUTER THE=ACOS(CSTHE(J)) 30 SNTHE(J)=SIN(THE) DO 40 L=1,NOUTR2 PHE=L*HE SNP=SIN(PHE) CSP=COS(PHE) SNPHE(1,L)=SNP CSPHE(1,L)=CSP DO 40 M=2,DEGREE SNPHE(M,L)=SNPHE(M-1,L)*CSP+CSPHE(M-1,L)*SNP 40 CSPHE(M,L)=CSPHE(M-1,L)*CSP-SNPHE(M-1,L)*SNP C C INITIALIZE FOR BEGINNING MAIN LOOP TO COMPUTE ALL (TNM,K(TRS)). C SET KERMAT TO ZERO. DO 50 I=1,DEGSQ DO 50 J=1,DEGSQ 50 KERMAT(I,J)=ZERO C C ********* C MAIN LOOP C ********* C DO 150 I=1,NOUTER CALL LEGEND(CSTHE(I),DEGREE,VALGDE) DO 150 J=1,NOUTR2 C PRODUCE SURFACE COORDINATES OF EXTERNAL VARIABLE POINT P=(X,Y,Z). ST=SNTHE(I) CT=CSTHE(I) SP=SNPHE(1,J) CP=CSPHE(1,J) CALL SURFAC(P,DUMMY,ST,CT,SP,CP,0) C C PRODUCE NEEDED DATA FOR ROTATION OF COORDINATES. UPX=CP*ST UPY=SP*ST UPZ=CT ALPHA=SIGN(ONE,-UPZ) V(3)=SQRT((ONE-UPZ/ALPHA)/TWO) PV=-TWO*ALPHA*V(3) V(2)=UPY/PV V(1)=UPX/PV C C INITIALIZE TEMPKH TO ZERO. DO 60 II=1,DEGSQ 60 TEMPKH(II)=ZERO C C LOOP TO INTEGRATE K(TRS) EVALUATED AT (THETA(I),PHI(J)). DO 100 II=1,NINNER DO 100 JJ=1,NINNR2 UPX=CSPHI(JJ)*SNTHI(II) UPY=SNPHI(JJ)*SNTHI(II) UPZ=CSTHI(II) C PRODUCE ROTATED POINT CALL ROTATE C C PRODUCE NEEDED TRIG AND LEGENDRE FUNCTIONS. CALL LEGEND(CTI,DEGREE,VALGDI) CSPHI2(1)=CPI SNPHI2(1)=SPI DO 70 M=2,DEGREE CSPHI2(M)=CPI*CSPHI2(M-1)-SPI*SNPHI2(M-1) 70 SNPHI2(M)=SPI*CSPHI2(M-1)+CPI*SNPHI2(M-1) C C CALCULATE SURFACE POSITION CALL SURFAC(Q,QCROSS,STI,CTI,SPI,CPI,1) TEMPK=KERNEL(P,Q,QCROSS) C DO 80 NR=1,DEGREE+1 80 TEMPKH(NR)=TEMPKH(NR)+WI(II)*TEMPK*(VALGDI(NR)-VALGDE(NR)) IBRS=DEGREE IB=DEGREE+1 DO 90 MS=1,DEGREE DO 90 NR=MS,DEGREE IB=IB+1 IBRS=IBRS+2 TEMPKH(IBRS)=TEMPKH(IBRS)+WI(II)*TEMPK*(VALGDI(IB) * *CSPHI2(MS)-VALGDE(IB)*CSPHE(MS,J)) 90 TEMPKH(IBRS+1)=TEMPKH(IBRS+1)+WI(II)*TEMPK* * (VALGDI(IB)*SNPHI2(MS)-VALGDE(IB)*SNPHE(MS,J)) 100 CONTINUE C C ADD BACK IN THE CONTRIBUTION SUBTRACTED EARLIER. DO 110 NR=1,DEGREE+1 110 TEMPKH(NR)=HI*TEMPKH(NR)+TWOPI*VALGDE(NR) IB=DEGREE+1 IBRS=DEGREE DO 120 MS=1,DEGREE DO 120 NR=MS,DEGREE IB=IB+1 IBRS=IBRS+2 TEMPKH(IBRS)=HI*TEMPKH(IBRS)+TWOPI*VALGDE(IB)*CSPHE(MS,J) 120 TEMPKH(IBRS+1)=HI*TEMPKH(IBRS+1) * +TWOPI*VALGDE(IB)*SNPHE(MS,J) C C PRODUCE THE NEW CONTRIBUTION, AT (TH(I),PH(J)), TO THE C ELEMENTS (TNM,K(TRS)) STORED IN KERMAT. IBNM=0 IB=0 TEMP1=WE(I)*HE DO 130 N=1,DEGREE+1 IBNM=IBNM+1 IB=IB+1 TEMP=TEMP1*VALGDE(N) DO 130 IBRS=1,DEGSQ 130 KERMAT(IBNM,IBRS)=KERMAT(IBNM,IBRS)+TEMP*TEMPKH(IBRS) IBNM=IBNM-1 DO 140 M=1,DEGREE DO 140 N=M,DEGREE IBNM=IBNM+2 IB=IB+1 TEMP=TEMP1*VALGDE(IB) DO 140 IBRS=1,DEGSQ KERMAT(IBNM,IBRS)=KERMAT(IBNM,IBRS) * +TEMP*CSPHE(M,J)*TEMPKH(IBRS) 140 KERMAT(IBNM+1,IBRS)=KERMAT(IBNM+1,IBRS) * +TEMP*SNPHE(M,J)*TEMPKH(IBRS) 150 CONTINUE C C OUTPUT VALUES OF KERMAT, AND ALSO STORE THEM IN COEFF. DO 200 I=1,DEGSQ CALL INDX(I,DEGREE,N,M,L) WRITE(FILEDC,1003)I,N,M,L 1003 FORMAT(' ROW=',I2,5X,'N,M,L=',3I3) WRITE(FILEDC,1002)(KERMAT(I,J),J=1,DEGSQ) 1002 FORMAT(1PD20.10,5D20.10) 200 CONTINUE RETURN END SUBROUTINE ROTATE C ----------------- C C CALCULATE (TX,TY,TZ)=(I-2*V*VT)*(X,Y,Z), THE C HOUSEHOLDER TRANSFORMATION OF (X,Y,Z). C IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/BLOCKR/V,ST,CT,SP,CP,X,Y,Z DIMENSION V(3) COMMON/MACHCN/U100,U100SQ C THE VARIABLE U100 IS 100 TIMES THE MACHINE EPSILON. PARAMETER(ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0) C C CALCULATE ROTATED POINT T. VM=TWO*(V(1)*X+V(2)*Y+V(3)*Z) TX=X-VM*V(1) TY=Y-VM*V(2) TZ=Z-VM*V(3) C C CALCULATE TRIG FUNCTIONS FOR ANGLES DETERMINED BY T. CT=TZ ST=SQRT(TX*TX+TY*TY) IF(ST .LE. U100) THEN CP=ONE SP=ZERO ELSE SP=TY/ST CP=TX/ST END IF RETURN END FUNCTION KERNEL(P,Q,QCROSS) C --------------- C C EVALUATE THE NORMAL DERIVATIVE OF 1/ABS(P-Q) WITH RESPECT TO THE C POINT Q ON THE SURFACE S. ALSO INCLUDE THE EFFECT OF THE CHANGE C OF SURFACE DIFFERENTIAL IN TRANSFORMING TO AN INTEGRAL OVER THE C UNIT SPHERE U. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION KERNEL DIMENSION P(3),Q(3),QCROSS(3) COMMON/MACHCN/U100,U100SQ C THE VARIABLE U100 IS 100 TIMES THE MACHINE EPSILON. PARAMETER(ZERO=0.0D0) C RS=(P(1)-Q(1))**2+(P(2)-Q(2))**2+(P(3)-Q(3))**2 IF(RS .LE. U100SQ) THEN KERNEL=ZERO ELSE DENOM=RS*SQRT(RS) TEMP=ZERO DO 10 J=1,3 10 TEMP=TEMP+(P(J)-Q(J))*QCROSS(J) KERNEL=TEMP/DENOM END IF RETURN END SUBROUTINE LEGEND(X,NMAX,VALUES) C ----------------- C C THIS ROUTINE WILL EVALUATE THE NORMALIZED ASSOCIATED LEGENDRE C FUNCTIONS OF THE FIRST KIND AT X, C J C Y(N,J)=P (X) C N C C THE DEFINITION IS TAKEN FROM CHAPTER 8 OF ABRAMOWITZ AND STEGUN, C HANDBOOK OF MATHEMATICAL FUNCTIONS. THESE WILL BE EVALUATED AND C STORED IN THE VECTOR 'VALUES' FOR 0 .LE. J .LE. N .LE. NMAX, WITH C NMAX .GE. 1. THE VALUE Y(N,J) IS STORED IN JB+N-J, WITH C JB=(2+(2*NMAX+3)*J-J**2)/2 C JB IS THE INDEX FOR Y(J,J). THE NORMALIZATION IS DONE SO THAT C THE SUBSEQUENT CONSTRUCTION OF SPHERICAL HARMONICS WILL RESULT C IN FUNCTIONS OF NORM 1 OVER THE UNIT SPHERE. C THE EVALUATION IS DONE USING THE TRIPLE RECURSION RELATION C (8.5.3) ON PAGE 334 OF THE ABOVE REFERENCE. THIS IS PRECEDED BY C A SPECIAL LOOP TO SET Y(N,N) AND Y(N,N-1). FOLLOWING THAT, THE C FUNCTIONS ARE NORMALIZED APPROPRIATELY. C C NOTE: THERE IS AN UPPER LIMIT ON NMAX, SPECIFIED BY THE C VARIABLE 'MAXDEG' IN THE PARAMETER STATEMENT BELOW. TO C INCREASE THE UPPER LIMIT, MERELY INCREASE 'MAXDEG'. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(MAXDEG=20,LFORD=2*MAXDEG) DIMENSION VALUES(*),FACT(0:LFORD) PARAMETER(ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0) C C SET INITIAL VALUES Y(N,N) AND Y(N,N-1), 0 .LE. N .LE. NMAX. RT=SQRT(ONE-X*X) VALUES(1)=ONE VALUES(2)=X IF(NMAX .EQ. 1) THEN VALUES(3)=-RT GO TO 30 END IF C JB=1 DO 10 J=1,NMAX-1 JBN=JB+NMAX-J+2 VALUES(JBN)=-(2*J-1)*RT*VALUES(JB) VALUES(JBN+1)=-(2*J+1)*RT*VALUES(JB+1) 10 JB=JBN VALUES(JB+2)=-(2*NMAX-1)*RT*VALUES(JB) C C PRODUCE THE REMAINDER OF THE TABLE VALUES Y(N,J). JB=1 DO 20 J=0,NMAX-2 JB=JB+2 DO 20 N=J+2,NMAX VALUES(JB)=((2*N-1)*X*VALUES(JB-1)-(N+J-1)*VALUES(JB-2))/(N-J) 20 JB=JB+1 C C NORMALIZE THE LEGENDRE FUNCTIONS. 30 FACT(0)=ONE DO 40 J=1,2*NMAX 40 FACT(J)=J*FACT(J-1) IB=0 PI=4*ATAN(ONE) TWOPI=TWO*PI FOURPI=FOUR*PI DO 50 N=0,NMAX 50 VALUES(N+1)=VALUES(N+1)*SQRT((2*N+1)/FOURPI) JB=NMAX+1 DO 60 M=1,NMAX DO 60 N=M,NMAX JB=JB+1 COEFF=SQRT(((2*N+1)/TWOPI)*FACT(N-M)/FACT(N+M)) 60 VALUES(JB)=COEFF*VALUES(JB) RETURN END FUNCTION LOCATN(N,M,L,NDEG) C --------------- C C GIVES THE INDEX OF THE SPHERICAL HARMONIC BASIS FUNCTION C WITH INDICES (N,M,L). FOR Q=(X,Y,Z) ON THE UNIT SPHERE, WRITE C Q = (SIN(THETA)*COS(PHI),SIN(THETA)*SIN(PHI),COS(THETA)) C FOR M=0, THE SPHERICAL HARMONIC CORRESPONDING TO (N,M,L) IS C C H(Q) = P (Z). C N C C FOR M>0 AND L=0, THE SPHERICAL HARMONIC IS C C M C H(Q) = P (Z)*COS(M*PHI) C N C C FOR M>0 AND L=1, THE SPHERICAL HARMONIC IS C C M C H(Q) = P (Z)*SIN(M*PHI) C N C C M C THE FUNCTIONS P (Z) AND P (Z) ARE THE LEGENDRE POLYNOMIALS AND C N N C THE ASSOCIATED LEGENDRE FUNCTIONS, DEFINED IN SUBROUTINE 'LEGEND'. C IF(M .EQ. 0) THEN LOCATN=N+1 ELSE LOCATN=NDEG*(2*M-1)-M*(M-3)+2*(N-M)+L END IF RETURN END SUBROUTINE INDX(I,NDEG,N,M,L) C --------------- C C THIS PROGRAM TAKES THE INDEX I OF A SPHERICAL HARMONIC, C 1 .LE. I .LE. (NDEG+1)**2, C AND IT FINDS THE INDICES (N,M,L) TO WHICH I CORRESPONDS. C SEE FUNCTION 'INDX' FOR THE DEFINITION OF (N,M,L). C KOL(MM)=NDEG*(MM-1)-(MM*(MM-3))/2-1 C IF(I .LE. NDEG+1) THEN N=I-1 M=0 L=0 RETURN END IF C J=I-NDEG-2 IF(J .EQ. 2*(J/2)) THEN L=0 ELSE L=1 END IF J=(J-L)/2 DO 40 M=1,NDEG IF(J .LT. KOL(M+1)) THEN N=J-KOL(M)+M RETURN END IF 40 CONTINUE END SUBROUTINE EXPAND(KERMAT,DUMMY,DEGREE,DEGSQ,NINNER,NOUTER, * FILEFF) C ----------------- C C THIS PROGRAM EXPANDS A TABLE OF GALERKIN COEFFICIENTS FROM THOSE FOR C DEGREES .LE. N:=DEGREE, AS INPUT. ON OUTPUT, THERE WILL BE N TABLES, C FOR GALERKIN COEFFICIENTS OF DEGREES .LE. L, WITH 1 .LE. L .LE. N. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER DEGREE,DEGSQ,FILEFF DOUBLE PRECISION KERMAT DIMENSION KERMAT(DEGSQ,DEGSQ),DUMMY(DEGSQ,DEGSQ) C REWIND FILEFF C C LOOP TO CONSTRUCT TABLES FOR SMALLER VALUES OF DEGREE. DO 20 ND=1,DEGREE NB=(ND+1)**2 DO 10 MR=0,ND DO 10 NR=MR,ND DO 10 LR=0,1 INEW=LOCATN(NR,MR,LR,ND) IOLD=LOCATN(NR,MR,LR,DEGREE) DO 10 MC=0,ND DO 10 NC=MC,ND DO 10 LC=0,1 JNEW=LOCATN(NC,MC,LC,ND) JOLD=LOCATN(NC,MC,LC,DEGREE) DUMMY(INEW,JNEW)=KERMAT(IOLD,JOLD) 10 CONTINUE WRITE(FILEFF)ND,NINNER,NOUTER 20 WRITE(FILEFF)((DUMMY(I,J),J=1,NB),I=1,NB) RETURN END SUBROUTINE INTMI(A,B,N,X,W) C ---------------- C C CALCULATE THE NODES AND WEIGHTS FOR THE INTEGRATION C C B N C INT F(X)*DX APPROX= SUM W(I)*F(X(I)) C A I=1 C C THE METHOD IS DUE TO IRI, MORIGUTI, AND TAKASAWA. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(ONE=1.0D0,TWO=2.0D0) DIMENSION X(*),W(*) C C INITIALIZE H=TWO/(N+1) IF(N .EQ. 2*(N/2)) THEN M=N/2 ELSE M=(N+1)/2 END IF FAC=B-A FAC2=H*(B-A) C C CALCULATE NODES AND WEIGHTS DO 10 I=1,M CALL TRNSFM(I*H-ONE,T,S) X(I)=A+T*FAC X(N+1-I)=B-T*FAC W(I)=FAC2*S W(N+1-I)=W(I) 10 CONTINUE RETURN END SUBROUTINE TRNSFM(Z,FI,SI) C ----------------- C C.......................................................................... C C FOR GIVEN Z, CALCULATE C C SI(Z) = CON*EXP(-4/(1-Z*Z)) C C AND C Z C FI(Z) = INT SI(T)*DT C -1 C C WHERE CON IS CHOSEN TO NORMALIZE FI(1)=1. C FI(Z) IS COMPUTED USING A TRUNCATED CHEBYSHEV SERIES. C C THIS ROUTINE IS TAKEN FROM THE PAPER: C IAN ROBINSON AND ELISE DE DONCKER, ALGORITHM 45: AUTOMATIC C COMPUTATION OF IMPROPER INTEGRALS OVER A BOUNDED OR C UNBOUNDED PLANAR REGION, COMPUTING 27(1981), PP. 253-284. C C.......................................................................... C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER ARG(5) DIMENSION A1(4),A2(4),Q(68),Q0300(18),Q0603(16),Q0806(15), * Q1008(19) PARAMETER(C40D3=40.0D0/3.0D0) PARAMETER(ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0) PARAMETER(DIVSOR=7.0298584066096D-3,DV1=TWO*DIVSOR, * DV2=FOUR*DIVSOR) EQUIVALENCE(Q(1),Q1008(1)),(Q(20),Q0806(1)),(Q(35),Q0603(1)), * (Q(51),Q0300(1)) DATA Q1008/0.1307592530313668D-01, 0.8580903946529252D-02, * 0.1969565496511163D-02, -0.6785908047534174D-04, * 0.5157989084218512D-05, -0.3254283578286669D-06, * 0.3009165432294816D-07, -0.2853960792992456D-08, * 0.3125316663674964D-09, -0.369071915636499D-10, * 0.47226104267733D-11, -0.6452090127750D-12, * 0.935240772542D-13, -0.142880027692D-13, * 0.22888743282D-14, -0.3827976860D-15, * 0.665889416D-16, -0.120097537D-16, * 0.22395647D-17/ DATA Q0806/0.7530091699129213D-01, 0.2215896190278894D-01, * 0.1517662771681320D-02, -0.1374204580631322D-04, * 0.2501181491115358D-05, -0.2491206397236787D-07, * 0.5430948732810670D-08, -0.1406070367942514D-09, * 0.160397857131033D-10, -0.7706914621139D-12, * 0.656131517151D-13, -0.43257167630D-14, * 0.3459964512D-15, -0.264283638D-16, * 0.21607480D-17/ DATA Q0603/0.2285305746758029D+00, 0.5670149715650661D-01, * 0.3881611773394649D-02, 0.1483758828946990D-03, * 0.1986416462810431D-04, 0.1166284710859293D-05, * 0.1048168134503124D-06, 0.6572636994171403D-08, * 0.5344588684897204D-09, 0.335115172128537D-10, * 0.26225503158527D-11, 0.1606252080762D-12, * 0.124685606769D-13, 0.73251000D-15, * 0.579829277D-16, 0.31817034D-17/ DATA Q0300/0.5404466499651320D+00, 0.1034761954005156D+00, * 0.9090551216566596D-02, 0.9130257654553327D-03, * 0.1026685176623610D-03, 0.1035244509501565D-04, * 0.1034925154934048D-05, 0.1002050044392230D-06, * 0.9552837592432324D-08, 0.8941334072231475D-09, * 0.825727197051832D-10, 0.75255600668576D-11, * 0.6783483380205D-12, 0.605164239084D-13, * 0.53497536298D-14, 0.4689365198D-15, * 0.407909118D-16, 0.35230149D-17/ DATA A1/20.0D0,20.0D0,C40D3,C40D3/, * A2/18.0D0,14.0D0,6.0D0,2.0D0/, * ARG/1,20,35,51,69/ DATA P1/-.3D0/,P2/-.6D0/,P3/-.8D0/ C C SET THE MACHINE DEPENDENT CONSTANTS. UFLOW=2*D1MACH(1) EOFLOW=LOG(D1MACH(2)/2) C C 'UFLOW' IS ABOUT EQUAL TO, BUT MUST NOT BE LESS THAN THE C SMALLEST POSITIVE REAL NUMBER REPRESENTABLE BY C THE MACHINE. C 'EOFLOW' IS ABOUT EQUAL TO, BUT MUST NOT EXCEED THE NAPERIAN C LOGARITHM OF THE MACHINE'S LARGEST REAL NUMBER. C C Q CONTAINS THE COEFFICIENTS OF THE CHEBYSHEV EXPANSION OF FI. C TO ACHIEVE GREATER THAN 16-DIGIT ACCURACY, IT IS NECESSARY TO C INCREASE THE NUMBER OF DIGITS IN THESE COEFFICIENTS AND ALSO C THE NUMBER OF COEFFICIENTS USED IN THE CHEBYSHEV SERIES. THIS C REQUIRES CHANGES TO THE DIMENSION OF Q, Q0300, Q0603, Q0806, C AND Q1008, TO THE VALUES OF ARG, AND TO THE CORRESPONDING C SUBSCRIPTS OF Q IN THE EQUIVALENCE STATEMENT. REFER TO C APPENDIX 3 OF THE FOLLOWING FOR DETAILS: C A HAEGEMANS, ALGORITHM 34: AN ALGORITHM FOR THE AUTOMATIC C INTEGRATION OVER A TRIANGLE, COMPUTING 19(1977),PP. 179-187. C C SET THE ARGUMENT X FOR EVALUATION. C IF(Z .LE. ZERO) THEN X=Z ELSE X=-Z END IF C C COMPUTE SI. C IF(X .LE. -ONE) THEN FI=ZERO SI=ZERO GO TO 20 ELSE EXPONT=FOUR/(ONE-X*X) IF(EXPONT .LE. EOFLOW) THEN SI=EXP(-EXPONT) ELSE FI=ZERO SI=ZERO GO TO 20 END IF END IF C C DETERMINE COEFFICIENTS TO BE USED IN CHEBYSHEV SERIES C IF(X .GE. P1) THEN INT=4 ELSE IF(X .GE. P2) THEN INT=3 ELSE IF(X .GE. P3) THEN INT=2 ELSE INT=1 END IF END IF END IF Y=A1(INT)*X+A2(INT) KFIRST=ARG(INT) KLAST=ARG(INT+1)-2 C C COMPUTE FI. C T2=ZERO T3=Q(KLAST+1) KSUM=KFIRST+KLAST DO 10 K=KFIRST,KLAST KREV=KSUM-K T1=T2 T2=T3 T3=Y*T2-T1+Q(KREV) 10 CONTINUE FAC=(T3-T1)/DV2 IF(FAC .GT. UFLOW/SI) THEN FI=FAC*SI ELSE FI=ZERO END IF 20 CONTINUE SI=SI/DV1 IF(Z .GT. ZERO) FI=ONE-FI RETURN END SUBROUTINE ZEROLG(N,ZZ,WW,A,B) C ----------------- C C PRODUCE THE NODES AND WEIGHTS FOR GAUSS-LEGENDRE QUADRATURE C ON (A,B), OF ORDER N. THE NUMBER OF NODES N IS RESTRICTED TO C BE .LE. 20. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION Z(109),W(109),ZZ(N),WW(N) DATA ONE/1.0D0/,TWO/2.0D0/ DATA Z(1),Z(2),Z(3),Z(4),Z(5),Z(6),Z(7),Z(8),Z(9),Z(10)/ *.577350269189626D0,.774596669241483D0,0.0D0,.861136311594053D0, *.339981043584856D0,.906179845938664D0,.538469310105683D0,0.0D0, *.932469514203152D0,.661209386466265D0/ DATA Z(11),Z(12),Z(13),Z(14),Z(15),Z(16),Z(17),Z(18),Z(19),Z(20)/ *.238619186083197D0,.949107912342759D0,.741531185599394D0, *.405845151377397D0,0.0D0,.960289856497536D0,.796666477413627D0, *.525532409916329D0,.183434642495650D0,.968160239507626D0/ DATA Z(21),Z(22),Z(23),Z(24),Z(25),Z(26),Z(27),Z(28),Z(29)/ *.836031107326636D0,.613371432700590D0,.324253423403809D0, *0.0D0,.973906528517172D0,.865063366688985D0,.679409568299024D0, *.433395394129247D0,.148874338981631D0/ DATA Z(30),Z(31),Z(32),Z(33),Z(34),Z(35),Z(36),Z(37),Z(38)/ *.978228658146057D0,.887062599768095D0,.730152005574049D0, *.519096129206812D0,.269543155952345D0,0.0D0, *.981560634246719D0,.904117256370475D0,.769902674194305D0/ DATA Z(39),Z(40),Z(41),Z(42),Z(43),Z(44),Z(45),Z(46),Z(47)/ *.587317954286617D0,.367831498998180D0,.125233408511469D0, *.984183054718588D0,.917598399222978D0,.801578090733310D0, *.642349339440340D0,.448492751036447D0,.230458315955135D0/ DATA Z(48),Z(49),Z(50),Z(51),Z(52),Z(53),Z(54),Z(55),Z(56)/ *0.0D0,.986283808696812D0,.928434883663574D0, *.827201315069765D0,.687292904811685D0,.515248636358154D0, *.319112368927890D0,.108054948707344D0,.987992518020485D0/ DATA Z(57),Z(58),Z(59),Z(60),Z(61),Z(62),Z(63),Z(64),Z(65)/ *.937273392400706D0,.848206583410427D0,.724417731360170D0, *.570972172608539D0,.394151347077563D0,.201194093997435D0, *0.0D0,.989400934991650D0,.944575023073233D0/ DATA Z(66),Z(67),Z(68),Z(69),Z(70),Z(71),Z(72),Z(73),Z(74)/ *.865631202387832D0,.755404408355003D0,.617876244402644D0, *.458016777657227D0,.281603550779259D0,.0950125098376374D0, *.990575475314417D0,.950675521768768D0,.880239153726986D0/ DATA Z(75),Z(76),Z(77),Z(78),Z(79),Z(80),Z(81),Z(82),Z(83)/ *.781514003896801D0,.657671159216691D0,.512690537086477D0, *.351231763453876D0,.178484181495848D0,0.0D0, *.991565168420931D0,.955823949571398D0,.892602466497556D0/ DATA Z(84),Z(85),Z(86),Z(87),Z(88),Z(89),Z(90),Z(91),Z(92)/ *.803704958972523D0,.691687043060353D0,.559770831073948D0, *.411751161462843D0,.251886225691506D0,.0847750130417353D0, *.992406843843584D0,.960208152134830D0,.903155903614818D0/ DATA Z(93),Z(94),Z(95),Z(96),Z(97),Z(98),Z(99),Z(100),Z(101)/ *.822714656537143D0,.720966177335229D0,.600545304661681D0, *.464570741375961D0,.316564099963630D0,.160358645640225D0, *0.0D0,.993128599185095D0,.963971927277914D0/ DATA Z(102),Z(103),Z(104),Z(105),Z(106),Z(107),Z(108),Z(109)/ *.912234428251326D0,.839116971822219D0,.746331906460151D0, *.636053680726515D0,.510867001950827D0,.373706088715420D0, *.227785851141645D0,.0765265211334973D0/ DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9),W(10)/ *1.0D0,.555555555555556D0,.888888888888889D0,.347854845137454D0, *.652145154862546D0,.236926885056189D0,.478628670499366D0, *.568888888888889D0,.171324492379170D0,.360761573048139D0/ DATA W(11),W(12),W(13),W(14),W(15),W(16),W(17),W(18),W(19),W(20)/ *.467913934572691D0,.129484966168870D0,.279705391489277D0, *.381830050505119D0,.417959183673469D0,.101228536290376D0, *.222381034453374D0,.313706645877887D0,.362683783378362D0, *.0812743883615744D0/ DATA W(21),W(22),W(23),W(24),W(25),W(26),W(27),W(28),W(29)/ *.180648160694857D0,.260610696402935D0,.312347077040003D0, *.330239355001260D0,.0666713443086881D0,.149451349150581D0, *.219086362515982D0,.269266719309996D0,.295524224714753D0/ DATA W(30),W(31),W(32),W(33),W(34),W(35),W(36),W(37),W(38)/ *.0556685671161737D0,.125580369464905D0,.186290210927734D0, *.233193764591990D0,.262804544510247D0,.272925086777901D0, *.0471753363865118D0,.106939325995318D0,.160078328543346D0/ DATA W(39),W(40),W(41),W(42),W(43),W(44),W(45),W(46),W(47)/ *.203167426723066D0,.233492536538355D0,.249147045813403D0, *.0404840047653159D0,.0921214998377284D0,.138873510219787D0, *.178145980761946D0,.207816047536889D0,.226283180262897D0/ DATA W(48),W(49),W(50),W(51),W(52),W(53),W(54),W(55),W(56)/ *.232551553230874D0,.0351194603317519D0,.0801580871597602D0, *.121518570687903D0,.157203167158194D0,.185538397477938D0, *.205198463721296D0,.215263853463158D0,.0307532419961173D0/ DATA W(57),W(58),W(59),W(60),W(61),W(62),W(63),W(64),W(65)/ *.0703660474881081D0,.107159220467172D0,.139570677926154D0, *.166269205816994D0,.186161000015562D0,.198431485327112D0, *.202578241925561D0,.0271524594117541D0,.0622535239386478D0/ DATA W(66),W(67),W(68),W(69),W(70),W(71),W(72),W(73),W(74)/ *.0951585116824928D0,.124628971255534D0,.149595988816577D0, *.169156519395003D0,.182603415044924D0,.189450610455068D0, *.0241483028685479D0,.0554595293739872D0,.0850361483171792D0/ DATA W(75),W(76),W(77),W(78),W(79),W(80),W(81),W(82),W(83)/ *.111883847193404D0,.135136368468525D0,.154045761076810D0, *.168004102156450D0,.176562705366993D0,.179446470356207D0, *.0216160135264833D0,.0497145488949698D0,.0764257302548891D0/ DATA W(84),W(85),W(86),W(87),W(88),W(89),W(90),W(91),W(92)/ *.100942044106287D0,.122555206711478D0,.140642914670651D0, *.154684675126265D0,.164276483745833D0,.169142382963144D0, *.0194617882297265D0,.0448142267656996D0,.0690445427376412D0/ DATA W(93),W(94),W(95),W(96),W(97),W(98),W(99),W(100),W(101)/ *.0914900216224500D0,.111566645547334D0,.128753962539336D0, *.142606702173607D0,.152766042065860D0,.158968843393954D0, *.161054449848784D0,.0176140071391521D0,.0406014298003869D0/ DATA W(102),W(103),W(104),W(105),W(106),W(107),W(108),W(109)/ *.0626720483341091D0,.0832767415767047D0,.101930119817240D0, *.118194531961518D0,.131688638449177D0,.142096109318382D0, *.149172986472604D0,.152753387130726D0/ C C CALCULATE THE WEIGHTS AND NODES. C SCALE=(B-A)/TWO IF((N/2)*2 .EQ. N) THEN M=N/2 IBASE=M*M ELSE IBASE=(N*N-1)/4 M=(N-1)/2 ZZ(M+1)=(A+B)/TWO WW(M+1)=W(IBASE+M)*SCALE END IF DO 100 I=1,M T=Z(IBASE+I-1) ZZ(I)=(A*(ONE+T)+(ONE-T)*B)/TWO ZZ(N+1-I)=(A*(ONE-T)+(ONE+T)*B)/TWO WW(I)=W(IBASE+I-1)*SCALE 100 WW(N+1-I)=WW(I) RETURN END DOUBLE PRECISION FUNCTION D1MACH (I) C -------------------------------- C C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 860501 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE RETURN DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS. C***DESCRIPTION C C D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C D = D1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF D ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C DOUBLE-PRECISION MACHINE CONSTANTS C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C D1MACH( 5) = LOG10(B) C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C FOR FTN4 C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777777B / C C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C FOR FTN5 C C DATA SMALL(1) / O"00564000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C C DATA LARGE(1) / O"37757777777777777777" / C DATA LARGE(2) / O"37157777777777777777" / C C DATA RIGHT(1) / O"15624000000000000000" / C DATA RIGHT(2) / O"00000000000000000000" / C C DATA DIVER(1) / O"15634000000000000000" / C DATA DIVER(2) / O"00000000000000000000" / C C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777774B / C C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/ C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1), SMALL(2) / '20000000, '00000201 / C DATA LARGE(1), LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1), DIVER(2) / '20000000, '00000334 / C DATA LOG10(1), LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 2100 C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C MACHINE CONSTANTS FOR THE HP 9000 C C D1MACH(1) = 2.8480954D-306 C D1MACH(2) = 1.40444776D+306 C D1MACH(3) = 2.22044605D-16 C D1MACH(4) = 4.44089210D-16 C D1MACH(5) = 3.01029996D-1 C C DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B / C DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B / C DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B / C DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B / C DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 8388608, 0 / C DATA LARGE(1), LARGE(2) / 2147483647, -1 / C DATA RIGHT(1), RIGHT(2) / 612368384, 0 / C DATA DIVER(1), DIVER(2) / 620756992, 0 / C DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA SMALL(3), SMALL(4) / 0, 0 / C C DATA LARGE(1), LARGE(2) / 32767, -1 / C DATA LARGE(3), LARGE(4) / -1, -1 / C C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 0 / C C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA DIVER(3), DIVER(4) / 0, 0 / C C DATA LOG10(1), LOG10(2) / 16282, 8346 / C DATA LOG10(3), LOG10(4) / -31493, -12296 / C C DATA SMALL(1), SMALL(2) / O000200, O000000 / C DATA SMALL(3), SMALL(4) / O000000, O000000 / C C DATA LARGE(1), LARGE(2) / O077777, O177777 / C DATA LARGE(3), LARGE(4) / O177777, O177777 / C C DATA RIGHT(1), RIGHT(2) / O022200, O000000 / C DATA RIGHT(3), RIGHT(4) / O000000, O000000 / C C DATA DIVER(1), DIVER(2) / O022400, O000000 / C DATA DIVER(3), DIVER(4) / O000000, O000000 / C C DATA LOG10(1), LOG10(2) / O037632, O020232 / C DATA LOG10(3), LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 / C C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / C C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / C C MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING) C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION) C C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / C DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X / C DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X / C DATA LOG10(1), DIVER(2) / '3FD34413'X,'509F79FF'X / C C MACHINE CONSTANTS FOR THE IBM PC - MICROSOFT FORTRAN C C DATA SMALL(1), SMALL(2) / #00000000, #00100000 / C DATA LARGE(1), LARGE(2) / #FFFFFFFF, #7FEFFFFF / C DATA RIGHT(1), RIGHT(2) / #00000000, #3CA00000 / C DATA DIVER(1), DIVER(2) / #00000000, #3CB00000 / C DATA LOG10(1), LOG10(2) / #509F79FF, #3FD34413 / C C MACHINE CONSTANTS FOR THE IBM PC - PROFESSIONAL FORTRAN C AND LAHEY FORTRAN C C DATA SMALL(1), SMALL(2) / Z'00000000', Z'00100000' / C DATA LARGE(1), LARGE(2) / Z'FFFFFFFF', Z'7FEFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'00000000', Z'3CA00000' / C DATA DIVER(1), DIVER(2) / Z'00000000', Z'3CB00000' / C DATA LOG10(1), LOG10(2) / Z'509F79FF', Z'3FD34413' / C C MACHINE CONSTANTS FOR THE PRIME 750/850. {9/23/86 TSENG} C C DATA SMALL(1), SMALL(2) / :10000000000, :00000100000 / C DATA LARGE(1), LARGE(2) / :17777777777, :37777677602 / C DATA RIGHT(1), RIGHT(2) / :10000000000, :00000000122 / C DATA DIVER(1), DIVER(2) / :10000000000, :00000000123 / C DATA LOG10(1), LOG10(2) / :11504046502, :17572000177 / C C MACHINE CONSTANTS FOR THE APOLLO. {9/23/86 TSENG} C DATA SMALL(1), SMALL(2) / 8#00004000000, 8#00000000000 / DATA LARGE(1), LARGE(2) / 8#17773777777, 8#37777777777 / DATA RIGHT(1), RIGHT(2) / 8#07450000000, 8#00000000000 / DATA DIVER(1), DIVER(2) / 8#07454000000, 8#00000000000 / DATA LOG10(1), LOG10(2) / 8#07764642023, 8#12047674777 / C C C***FIRST EXECUTABLE STATEMENT D1MACH C IF (I .LT. 1 .OR. I .GT. 5) C 1 CALL XERROR( 'D1MACH -- I OUT OF BOUNDS',25,1,2) C D1MACH = DMACH(I) RETURN C END