PROGRAM DRIVER C C TITLE: TEST OF SUBROUTINE 'LAPLAC' C ---------------------------------- C C THIS IS AN INTERACTIVE DRIVER PROGRAM FOR TESTING THE SUBROUTINE C 'LAPLAC'. THE SUBROUTINE SOLVES THE INTERIOR DIRICHLET PROBLEM C FOR LAPLACE'S EQUATION ON A SIMPLY-CONNECTED REGION D WITH A C SMOOTH BOUNDARY S. THIS DRIVER PROGRAM IS SET UP TO SOLVE C PROBLEMS WITH KNOWN SOLUTIONS, GIVEN IN 'BDYFCN'; BUT IT CAN BE C USED WITH NO SIGNIFICANT CHANGE FOR THE CASE WHERE THE SOLUTION C IS UNKNOWN. C C THE PROGRAM WILL QUIZ THE USER FOR NEEDED INFORMATION ABOUT THE C PROBLEM. THE INFORMATION ASKED FOR IS AS FOLLOWS: C C NUMSUR,A,B,C,ALPHA: THESE VARIABLES DEFINE THE SURFACE S. THEY C WERE USED EARLIER IN CALCULATING THE GALERKIN COEFFICIENTS C THAT ARE NEEDED BY 'LAPLAC'. AS BEFORE, THE SURFACE IS C DEFINED BY SUBROUTINE 'SURFAC'. C FNAME: THIS IS THE NAME OF THE FILE CONTAINING THE GALERKIN C COEFFICIENTS, CALCULATED PREVIOUSLY WITH SUBROUTINE 'GNRT'. C THIS FILE MUST BE THE FORMATFREE OUTPUT FILE FROM 'GNRT'. C FILOUT: THIS IS THE NUMBER OF THE OUTPUT FILE, AND IT DEPENDS C ON WHETHER THE OUTPUT IS TO THE TERMINAL (FILOUT=TMLOUT) C OR TO A SEPARATE FILE (FILOUT=FILALT). IN THE LATTER CASE, C THE USER WILL BE ASKED TO SUPPLY A NAME FOR THE OUTPUT FILE. C THE ACTUAL UNIT NUMBERS USED ARE GIVEN IN THE DATA C STATEMENT BELOW. C NUMPTS: THE NUMBER OF POINTS IN THE REGION D AT WHICH THE C SOLUTION U (GIVEN IN USOLN) OF LAPLACE'S EQUATION IS BEING C SOUGHT. THERE IS AN UPPER LIMIT ON 'NUMPTS'. IT IS CALLED C 'MAXPTS', AND IT IS GIVEN BELOW IN THE PARAMETER STATEMENT. C POINTS: THIS ARRAY CONTAINS THE POINTS IN D AT WHICH THE C SOLUTION U IS TO BE FOUND. C NUMBF: THE INDEX OF THE BOUNDARY FUNCTION, GIVEN IN 'BDYFCN'. C AUXPAR: AN EXTRA PARAMETER THAT MAY BE USED IN DEFINING C THE BOUNDARY FUNCTION, IF DESIRED. C DEGREE: THE DEGREE OF THE SPHERICAL HARMONIC POLYNOMIAL TO BE C USED IN SOLVING THE DOUBLE LAYER INTEGRAL EQUATION, IN C SUBROUTINE 'LAPLAC'. THERE IS AN UPPER LIMIT FOR 'DEGREE'. C IT IS CALLED 'MAXDEG', AND IT IS GIVEN BELOW IN THE C PARAMETER STATEMENT. C NINTEG: THE INTEGRATION PARAMETER USED IN DEFINING THE C GAUSSIAN QUADRATURE OF THE DOUBLE LAYER POTENTIAL, IN C SUBROUTINE 'LAPLAC'. THERE IS AN UPPER LIMIT FOR 'NINTEG'. C IT IS CALLED 'MAXNIN', AND IT IS GIVEN BELOW IN THE C PARAMETER STATEMENT. THERE ARE TWO GAUSSIAN QUADRATURE C SUBROUTINES, 'ZEROLG' AND 'GAUSS2' BEING USED HERE. THE FIRST C COMPUTES THE GAUSS-LEGENDRE NODES AND WEIGHTS OF ORDER .LE. 20, C AND THE SECOND DOES SO FOR ORDERS .GE. 32, IN POWERS OF 2. C THUS IF 'NINTEG' IS CHOSEN LARGER THAN 20, IT SHOULD BE A C POWER OF 2. C C THE ARRAY 'WORK' GIVES THE WORKING STORAGE NEEDED BY 'LAPLAC'. C THE DIMENSION OF 'WORK' IS GIVEN IN THE PARAMETER STATEMENT C BELOW, BASED ON DEFAULT UPPER LIMITS FOR VARIOUS INPUT VARIABLES. C THERE IS ALSO A CHECK AFTER THE INPUT OF ALL VARIABLES THAT C THE DIMENSION OF 'WORK' IS ADEQUATE. C C MACHINE DEPENDENT CONSTANTS ARE SET THROUGH THE FUNCTION C 'D1MACH', WHICH IS ONE OF THE SUBPROGRAMS GIVEN IN THE C 'LAPLAC' PACKAGE. THE ROUTINE 'D1MACH' SHOULD BE RESET C APPROPRIATELY FOR THE USER'S COMPUTER. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER FILEIN, FILALT, FILOUT, DEGREE, TML_IN, TMLOUT CHARACTER FNAME*50, T, ANAME*50 PARAMETER(MAXDEG=9, MAXNIN=64 ,MAXPTS=20) PARAMETER(IWORK=(MAXDEG+1)**4+(MAXDEG+1)**2+2*MAXDEG * +MAXNIN*(4+(MAXDEG*(MAXDEG+3))/2)) DIMENSION USOLN(MAXPTS), POINTS(3,MAXPTS), WORK(IWORK) COMMON/BDYPAR/NUMBF, AUXPAR COMMON/BLKWRK/WORK COMMON/SURPAR/A, B, C, ALPHA, NUMSUR EXTERNAL SURFAC, BDYFCN C C ******* INPUT/OUTPUT ********************************** C THESE ARE THE I/O FILE NUMBERS, AND THEY MAY NEED C TO BE CHANGED FOR USE ON YOUR COMPUTER. DATA FILEIN/8/, TML_IN/5/, TMLOUT/6/, FILALT/9/ C ******************************************************* C C INPUT THE SURFACE PARAMETERS. PRINT *, ' GIVE THE SURFACE PARAMETERS (NUMSUR,A,B,C,ALPHA).' READ *, NUMSUR,A,B,C,ALPHA C C GET INFORMATION ABOUT INPUT AND OUTPUT FILES. PRINT *, ' GIVE THE NAME FOR YOUR FILE OF GALERKIN COEFFICIENTS.' READ(*,'(A)') FNAME L=LEN(FNAME) OPEN(FILEIN,FILE=FNAME(1:L),FORM='UNFORMATTED') PRINT *, ' DO YOU WANT THE ANSWERS WRITTEN TO THE TERMINAL (Y/N)?' READ(*,'(A1)') T IF((T .EQ. 'Y') .OR. (T .EQ. 'y')) THEN FILOUT=TMLOUT ELSE PRINT *,' GIVE THE NAME YOU WISH TO USE FOR YOUR OUTPUT FILE.' READ(*,'(A)') ANAME FILOUT=FILALT L=LEN(ANAME) OPEN(FILALT,FILE=ANAME(1:L)) END IF C C PRINT THE SURFACE AND GALERKIN COEFFICIENT PARAMETERS. REWIND FILEIN READ(FILEIN) NDEG,NINNER,NOUTER WRITE(FILOUT,1000) NUMSUR,A,B,C,ALPHA,NINNER,NOUTER 1000 FORMAT(//,' SURFACE PARAMETERS: SURFACE INDEX=',I1,/,' (A,B,C,)=', * 1PD17.10,2D20.10,/,' AUXILIARY SURFACE PARAMETER ALPHA = ', * D17.10,/,' GALERKIN COEFFICIENT PARAMETERS: NINNER,NOUTER=', * I3,I6,//) C C READ THE POINTS AT WHICH THE SOLUTION IS TO BE EVALUATED. PRINT *, ' AT HOW MANY POINTS INSIDE THE SURFACE IS THE SOLUTION' PRINT *, ' TO BE EVALUATED?' READ *, NUMPTS PRINT *, ' GIVE THE POINTS (X,Y,Z).' READ *, ((POINTS(I,J),I=1,3),J=1,NUMPTS) C C READ THE PARTICULAR PROBLEM PARAMETERS. THE PROGRAM WILL LOOP C BACK TO THIS POINT AS OFTEN AS DESIRED. 10 PRINT *, ' GIVE THE INDEX OF THE BOUNDARY FUNCTION, NUMBF, AND' PRINT *, ' THE BOUNDARY PARAMETER AUXPAR. TO STOP THE PROGRAM,' PRINT *, ' LET NUMBF=0.' READ *, NUMBF,AUXPAR IF(NUMBF .EQ. 0) THEN CLOSE(FILEIN) IF(FILOUT .EQ. FILALT) THEN CLOSE(FILALT) END IF STOP END IF PRINT *, ' GIVE THE DEGREE OF THE APPROXIMATION AND THE' PRINT *, ' INTEGRATION PARAMETER. ALSO, THE' PRINT *, ' INTEGRATION PARAMETER SHOULD BE A POWER OF TWO' PRINT *, ' AND LESS THAT OR EQUAL TO 256.' READ *, DEGREE,NINTEG C C CHECK ON THE NEEDED DIMENSION FOR 'WORK'. ICHECK=(DEGREE+1)**4+(DEGREE+1)**2+2*DEGREE * +NINTEG*(4+(DEGREE*(DEGREE+3))/2) IF(ICHECK .GT. IWORK) THEN C THERE IS INSUFFICIENT SPACE IN THE ARRAY 'WORK'. PRINT 1001, ICHECK 1001 FORMAT(' THE DIMENSION OF THE ARRAY WORK IS NOT SUFFICIENT.', * /,' THE DIMENSION SHOULD BE AT LEAST ',I5) STOP END IF C CALL LAPLAC(BDYFCN,SURFAC,DEGREE,NINTEG,FILEIN,NUMPTS,POINTS, * USOLN,0,WORK) C C PRINT RESULTS. WRITE(FILOUT,1002) NUMBF,AUXPAR,DEGREE,NINTEG 1002 FORMAT(///,' BOUNDARY FUNCTION INDEX=',I2,/,' AUXILIARY BOUNDARY', * ' FUNCTION PARAMETER=',1PD17.10,/,' DEGREE OF ', * 'APPROXIMATION=',I2,/,' INTEGRATION PARAMETER=',I3,/) WRITE(FILOUT,1003) IB=(DEGREE+1)**4+2*DEGREE+NINTEG*(4+(DEGREE*(DEGREE+3))/2) 1003 FORMAT(' LAPLACE EXPANSION COEFFICIENTS FOR DENSITY FUNCTION', * //,7X,'INDICES',11X,'COEFFICIENT',//,3X,'I',4X,'N',4X, * 'M',3X,'L') DO 15 I=1,(DEGREE+1)**2 CALL INDX(I,DEGREE,N,M,L) 15 WRITE(FILOUT,1004) I,N,M,L,WORK(IB+I) 1004 FORMAT(I4,I5,I5,I4,1PD21.10) WRITE(FILOUT,1005) 1005 FORMAT(//,' POTENTIAL SOLUTION U(X,Y,Z) AT DESIRED POINTS',//, * 2X,'I',7X,'X',9X,'Y',9X,'Z',12X,'U',18X,'ERROR') DO 20 I=1,NUMPTS TRUEU=BDYFCN(POINTS(1,I)) ERROR=TRUEU-USOLN(I) 20 WRITE(FILOUT,1006) I,(POINTS(J,I),J=1,3),USOLN(I),ERROR 1006 FORMAT(I3,3F10.4,1PD21.10,D14.3) GO TO 10 END FUNCTION BDYFCN(Q) C --------------- C C THIS DEFINES THE BOUNDARY VALUE FUNCTION OF THE UNKNOWN C SOLUTION OF LAPLACE'S EQUATION. IN THIS PARTICULAR INSTANCE, C IT ALSO GIVES THE TRUE SOLUTION OF LAPLACE'S EQUATION, FOR C DEMONSTRATION PURPOSES. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/BDYPAR/NUMBF,A DIMENSION Q(3) PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,FOUR=4.0D0, * FIVE=5.0D0) C X=Q(1) Y=Q(2) Z=Q(3) GO TO (10,20,30,40,50),NUMBF 10 BDYFCN=ONE RETURN 20 BDYFCN=X RETURN 30 BDYFCN=X*X+Y*Y-TWO*Z*Z RETURN 40 BDYFCN=EXP(A*X)*COS(A*Y)+EXP(A*Z)*SIN(A*X) RETURN 50 BDYFCN=ONE/SQRT((X-FIVE)**2+(Y-FOUR)**2+(Z-THREE)**2) RETURN 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/SURPAR/A,B,C,ALPHA,NUMSUR COMMON/MACHCN/U100,U100SQ C U100 IS 100 TIMES THE MACHINE UNIT ROUND. C C COMPUTE UQ. UX = ST*CP UY = ST*SP UZ = CT C C******************************************************************* 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 LAPLAC(BDYFCN,SURFAC,DEGREE,NINTEG,FILEIN,NUMPTS, * POINTS,USOLN,OPREPT,WORK) C ----------------- C C THIS CALCULATES THE SOLUTION OF THE INTERIOR DIRICHLET C PROBLEM FOR LAPLACE'S EQUATION ON A REGION D IN THREE C DIMENSIONS. THE REGION D IS ASSUMED TO BE SIMPLY-CONNECTED C WITH A SMOOTH BOUNDARY S. THE METHOD OF SOLUTION IS TO C REPRESENT THE UNKNOWN HARMONIC FUNCTION U(X,Y,Z) AS A DOUBLE C LAYER POTENTIAL. THIS IS FOUND BY NUMERICALLY SOLVING FOR THE C DOUBLE LAYER DENSITY FUNCTION RHO(X,Y,Z) ON S BY MEANS OF C GALERKIN'S METHOD WITH SPHERICAL HARMONICS AS BASIS FUNCTIONS. C THE METHOD IS DESCRIBED IN THE ARTICLE C K. ATKINSON, 'THE NUMERICAL SOLUTION OF LAPLACE'S EQUATION C IN THREE DIMENSIONS', SIAM J. NUM. ANAL. 19(1982),263-274. C THE DOUBLE LAYER POTENTIAL IS THEN EVALUATED BY USING GAUSSIAN C QUADRATURE OVER THE SPHERE, AS DESCRIBED IN C K. ATKINSON, 'NUMERICAL INTEGRATION ON THE SPHERE', J. C AUSTRALIAN MATH. SOC., (SERIES B) 23(1982),332-347. C C INPUT: C BDYFCN - THE BOUNDARY CONDITION FOR THE SOLUTION U. IT MUST BE C GIVEN AS A FUNCTION BDYFCN(Q), WITH Q A VECTOR OF LENGTH 3. C SURFAC - THIS IS THE NAME OF THE SUBROUTINE THAT DEFINES THE C SURFACE S. SEE THE EXAMPLE SUBROUTINE 'SURFAC' FOR DETAILS C OF THE DEFINITION AND THE CALLING SEQUENCE. C DEGREE - THE UPPER LIMIT ON THE DEGREE OF THE SPHERICAL HARMONICS C BEING USED IN THE APPROXIMATION OF RHO. C NINTEG - THE INTEGRATION PARAMETER BEING USED TO DETERMINE THE C GAUSSIAN QUADRATURE FORMULA FOR SPHERICAL INTEGRATION. AT C PRESENT, IF 'NINTEG' IS GREATER THAN 20, THEN IT SHOULD BE A C POWER OF 2. THIS IS DUE TO THE SUBROUTINE 'GAUSS2' BEING C USED TO PRODUCE THE GAUSS-LEGENDRE NODES AND WEIGHTS IN SUCH C A CASE. THE TOTAL NUMBER OF INTEGRATION NODES ON THE SURFACE C S WILL BE 2*NINTEG**2. C FILEIN - THE NUMBER OF THE INPUT FILE CONTAINING THE GALERKIN C COEFFICIENTS, AS GENERATED BY THE SUBROUTINE 'GNRT'. THIS C FILE MUST BE THE FORMATFREE FILE FROM 'GNRT'. C NUMPTS - THE NUMBER OF POINTS IN D AT WHICH THE SOLUTION U IS TO C BE EVALUATED. C POINTS - THE ARRAY CONTAINING THE POINTS FOR THE EVALUATION OF C U(X,Y,Z). THE COORDINATES (X,Y,Z) OF POINT #J ARE STORED C IN POINTS(1,J), POINTS(2,J), POINTS(3,J). C OPREPT - OPREPT=0 MEANS THAT THIS IS THE FIRST CALL ON 'LAPLAC' C WITH THIS SET OF PROBLEM PARAMETERS. C OPREPT=1 MEANS 'LAPLAC' HAS BEEN CALLED PREVIOUSLY, AND C IT IS DESIRED TO EVALUATE U AT A NEW SET OF POINTS IN D. C ONLY 'NUMPTS' AND 'POINTS' SHOULD HAVE BEEN CHANGED SINCE C THE PREVIOUS CALL ON 'LAPLAC'. C WORK - A TEMPORARY ARRAY OF WORKING STORAGE. ITS DIMENSION C SHOULD BE AT LEAST C (DEGREE+1)**4 + (DEGREE+1)**2 + 2*DEGREE C +NINTEG*(4+(DEGREE*(DEGREE+3))/2) C C OUTPUT: C USOLN - A ONE DIMENSIONAL ARRAY OF LENGTH AT LEAST 'NUMPTS'. C THE ELEMENT USOLN(I) CONTAINS THE APPROXIMATION TO U AT C POINT #I OF 'POINTS'. C GALERKIN COEFFICIENTS OF RHO - IF DESIRED, THESE ARE STORED IN C 'WORK', BEGINNING AT POSITION C (DEGREE+1)**4 + 2*DEGREE C + NINTEG*(4+(DEGREE*(DEGREE+3))/2) C THE SPHERICAL HARMONICS WITH WHICH THESE COEFFICIENTS ARE C ASSOCIATED CAN BE DETERMINED BY USING SUBROUTINE 'INDX' TO C CONVERT THE INDEX I OF THE COEFFICIENT #I. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER DEGREE,FILEIN,OPREPT,DEGSQ,DEG2 DIMENSION POINTS(3,*),USOLN(*),WORK(*) COMMON/MACHCN/U100,U100SQ C THE VARIABLE U100 IS 100 TIMES THE MACHINE UNIT ROUND. SAVE EXTERNAL BDYFCN,SURFAC C C SET MACHINE DEPENDENT CONSTANT. U100=100*D1MACH(4) U100SQ=U100**2 C IF(OPREPT .EQ. 1) GO TO 10 DEGSQ=(DEGREE+1)**2 DEG2=2*DEGREE+1 LGDORD=1+(DEGREE*(DEGREE+3))/2 C C SUBDIVIDE 'WORK' INTO SMALLER ONE AND TWO DIMENSIONAL ARRAYS FOR C USE IN SUBROUTINE 'LPLC2'. I1=1 I2=I1+DEGSQ**2 I3=I2+NINTEG I4=I3+NINTEG I5=I4+NINTEG I6=I5+DEGREE I7=I6+DEGREE I8=I7+LGDORD*NINTEG C 10 CALL LPLC2(BDYFCN,SURFAC,DEGREE,NINTEG,FILEIN,NUMPTS,POINTS,USOLN, * OPREPT,WORK(I1),WORK(I2),WORK(I3),WORK(I4),WORK(I5), * WORK(I6),WORK(I7),WORK(I8),DEGSQ,DEG2,LGDORD) RETURN END SUBROUTINE LPLC2(BDYFCN,SURFAC,DEGREE,NINTEG,FILEIN,NUMPTS, * POINTS,USOLN,OPREPT,KERMAT,CSTH,SNTH,W,CSPH,SNPH,VALUES, * RHS,DEGSQ,DEG2,LGDORD) C ---------------- C C THE ACTUAL SOLUTION OF LAPLACE'S EQUATION TAKES PLACE IN THIS ROUTINE. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION KERMAT,KERNEL INTEGER DEGREE,DEGSQ,DEG2,OPREPT,FILEIN DIMENSION KERMAT(DEGSQ,DEGSQ),CSTH(NINTEG),SNTH(NINTEG), * W(NINTEG),CSPH(DEGREE),SNPH(DEGREE),VALUES(LGDORD,NINTEG), * RHS(DEGSQ),POINTS(3,*),USOLN(*),Q(3),QCROSS(3) PARAMETER(ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0) SAVE C IF(OPREPT .EQ. 1) GO TO 200 C C MOVE POINTER OF INPUT FILE TO CORRECT SET OF TABLES. REWIND FILEIN DO 10 K=1,DEGREE-1 KT=(K+1)**4 READ(FILEIN)NDEG,NINNER,NOUTER 10 READ(FILEIN)(DUMMY,I=1,KT) C C READ MATRIX OF GALERKIN COEFFICIENTS. READ(FILEIN)NDEG,NINNER,NOUTER READ(FILEIN)((KERMAT(I,J),J=1,DEGSQ),I=1,DEGSQ) C C INITIALIZE. PI=FOUR*ATAN(ONE) TWOPI=TWO*PI HPHI=PI/NINTEG DO 20 I=1,DEGSQ 20 KERMAT(I,I)=KERMAT(I,I)+TWOPI C C INITIALIZE INTEGRATION NODES AND THE ASSOCIATED TRIG FUNCTIONS. IF(NINTEG .LE. 20) THEN CALL ZEROLG(NINTEG,CSTH,W,-ONE,ONE) ELSE CALL GAUSS2(-ONE,ONE,NINTEG,CSTH,W) END IF DO 30 I=1,NINTEG THETA=ACOS(CSTH(I)) SNTH(I)=SIN(THETA) 30 CALL LEGEND(CSTH(I),DEGREE,VALUES(1,I)) C C BEGIN SECTION TO CREATE RHS. C C CREATE RIGHT-HAND SIDES OF GALERKIN SYSTEM. DO 70 I=1,DEGSQ 70 RHS(I)=ZERO DO 100 J=1,2*NINTEG PHI=J*HPHI CP=COS(PHI) SP=SIN(PHI) SNPH(1)=SP CSPH(1)=CP DO 80 M=2,DEGREE CSPH(M)=CP*CSPH(M-1)-SP*SNPH(M-1) 80 SNPH(M)=SP*CSPH(M-1)+CP*SNPH(M-1) DO 100 I=1,NINTEG CALL SURFAC(Q,DUMMY,SNTH(I),CSTH(I),SP,CP,0) BF=BDYFCN(Q) DO 90 N=1,DEGREE+1 90 RHS(N)=RHS(N)+W(I)*VALUES(N,I)*BF IB=DEGREE+1 IBNM=DEGREE DO 100 M=1,DEGREE DO 100 N=M,DEGREE IB=IB+1 IBNM=IBNM+2 RHS(IBNM)=RHS(IBNM)+W(I)*BF*CSPH(M)*VALUES(IB,I) 100 RHS(IBNM+1)=RHS(IBNM+1)+W(I)*BF*SNPH(M)*VALUES(IB,I) DO 110 J=1,DEGSQ 110 RHS(J)=HPHI*RHS(J) C C END SECTION ON CREATING RHS. C C ********************************************************* C C SOLVE LINEAR SYSTEM FOR LAPLACE COEFFICIENTS CALL LINSYS(KERMAT,RHS,DEGSQ,DEGSQ,IER) C ********************************************************* C C INITIALIZE ALL VALUES OF U TO ZERO. 200 DO 210 L=1,NUMPTS 210 USOLN(L)=ZERO C C BEGIN LOOP TO CREATE VALUES OF USOLN AT POINTS (X,Y,Z) C GIVEN IN THE ARRAY 'POINTS'. C DO 250 J=1,2*NINTEG PHI=J*HPHI SP=SIN(PHI) CP=COS(PHI) SNPH(1)=SP CSPH(1)=CP DO 220 M=2,DEGREE SNPH(M)=SP*CSPH(M-1)+CP*SNPH(M-1) 220 CSPH(M)=CP*CSPH(M-1)-SP*SNPH(M-1) DO 250 I=1,NINTEG C CALCULATE RHO AT (THETA(I),PHI). RHO=ZERO DO 230 N=1,DEGREE+1 230 RHO=RHO+RHS(N)*VALUES(N,I) IB=DEGREE+1 IBNM=DEGREE DO 240 M=1,DEGREE DO 240 N=M,DEGREE IB=IB+1 IBNM=IBNM+2 240 RHO=RHO+VALUES(IB,I)*(CSPH(M)*RHS(IBNM) * +SNPH(M)*RHS(IBNM+1)) C CALCULATE VALUE OF U AT EACH (X,Y,Z) IN POINTS. DO 250 L=1,NUMPTS CALL SURFAC(Q,QCROSS,SNTH(I),CSTH(I),SP,CP,1) FM=KERNEL(POINTS(1,L),Q,QCROSS) USOLN(L)=USOLN(L)+W(I)*RHO*FM 250 CONTINUE DO 260 L=1,NUMPTS 260 USOLN(L)=HPHI*USOLN(L) 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 U100 IS 100 TIMES THE MACHINE UNIT ROUND. 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 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 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 SUBROUTINE GAUSS2(A,B,N,TT,WW) C ----------------- C C THE GAUSSIAN QUADRATURE NODES AND WEIGHTS ARE CALCULATED FOR THE FORMULA C C B I=N C INT F(T)*DT APPROX. = SUM WW(I)*F(TT(I)) C A I=1 C C THE NODES WILL BE STORED IN TT IN INCREASING ORDER, C A .LT. TT(1) .LT. ... .LT. TT(N) .LT. B. C THESE NODES ARE THE ZEROS OF THE LEGENDRE POLYNOMIAL ON (A,B) OF C DEGREE N. THE CORRESPONDING WEIGHT FOR TT(I) WILL BE STORED IN WW(I). C ***WARNING*** THE INTEGER N MUST BE A POWER OF TWO, C 2 .LE. N .LE. 256. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION WW(1),TT(1) DIMENSION W(255),T(255) DATA T(1),T(2),T(3),T(4),T(5),T(6),T(7),T(8),T(9),T(10),T(11), * T(12),T(13),T(14),T(15)/.577350269189626D0, * .861136311594053D0,.339981043584856D0,.960289856497536D0, * .796666477413627D0,.525532409916329D0,.183434642495650D0, * .989400934991650D0,.944575023073233D0,.865631202387832D0, * .755404408355003D0,.617876244402644D0,.458016777657227D0, * .281603550779259D0,.950125098376374D-1/ DATA T(16),T(17),T(18),T(19),T(20),T(21),T(22),T(23),T(24),T(25), * T(26),T(27),T(28),T(29),T(30),T(31)/.997263861849482D0, * .985611511545268D0,.964762255587506D0,.934906075937740D0, * .896321155766052D0,.849367613732570D0,.794483795967942D0, * .732182118740290D0,.663044266930215D0,.587715757240762D0, * .506899908932229D0,.421351276130635D0,.331868602282128D0, * .239287362252137D0,.144471961582796D0,.483076656877383D-1/ DATA T(32),T(33),T(34),T(35),T(36),T(37),T(38),T(39),T(40),T(41), * T(42),T(43),T(44),T(45),T(46),T(47)/.999305041735772D0, * .996340116771955D0,.991013371476744D0,.983336253884626D0, * .973326827789911D0,.961008799652054D0,.946411374858403D0, * .929569172131940D0,.910522137078503D0,.889315445995114D0, * .865999398154093D0,.840629296252580D0,.813265315122798D0, * .783972358943341D0,.752819907260532D0,.719881850171611D0/ DATA T(48),T(49),T(50),T(51),T(52),T(53),T(54),T(55),T(56),T(57), * T(58),T(59),T(60),T(61),T(62),T(63)/.685236313054233D0, * .648965471254657D0,.611155355172393D0,.571895646202634D0, * .531279464019895D0,.489403145707053D0,.446366017253464D0, * .402270157963992D0,.357220158337668D0,.311322871990211D0, * .264687162208767D0,.217423643740007D0,.169644420423993D0, * .121462819296121D0,.729931217877990D-1,.243502926634244D-1/ DATA T(64),T(65),T(66),T(67),T(68),T(69),T(70),T(71),T(72),T(73), * T(74),T(75),T(76),T(77),T(78),T(79)/.999824887947132D0, * .999077459977376D0,.997733248625514D0,.995792758534981D0, * .993257112900213D0,.990127818491734D0,.986406742724586D0, * .982096108435719D0,.977198491463907D0,.971716818747137D0, * .965654366431965D0,.959014757853700D0,.951801961341264D0, * .944020287830220D0,.935674388277916D0,.926769250878948D0/ DATA T(80),T(81),T(82),T(83),T(84),T(85),T(86),T(87),T(88),T(89), * T(90),T(91),T(92),T(93),T(94),T(95),T(96),T(97)/ * .917310198080961D0,.907302883401757D0,.896753288049158D0, * .885667717345397D0,.874052796958032D0,.861915468939548D0, * .849262987577969D0,.836102915060907D0,.822443116955644D0, * .808291757507914D0,.793657294762193D0,.778548475506412D0, * .762974330044095D0,.746944166797062D0,.730467566741909D0, * .713554377683587D0,.696214708369514D0,.678458922447719D0/ DATA T(98),T(99),T(100),T(101),T(102),T(103),T(104),T(105),T(106), * T(107),T(108),T(109),T(110),T(111),T(112),T(113),T(114)/ * .660297632272646D0,.641741692562308D0,.622802193910585D0, * .603490456158549D0,.583818021628763D0,.563796648226618D0, * .543438302412810D0,.522755152051175D0,.501759559136144D0, * .480464072404172D0,.458881419833552D0,.437024501037104D0, * .414906379552275D0,.392540275033267D0,.369939555349859D0, * .347117728597636D0,.324088435024413D0/ DATA T(115),T(116),T(117),T(118),T(119),T(120),T(121),T(122), * T(123),T(124),T(125),T(126),T(127)/.300865438877677D0, * .277462620177904D0,.253893966422694D0,.230173564226660D0, * .206315590902079D0,.182334305985337D0,.158244042714225D0, * .134059199461188D0,.109794231127644D0, * .854636405045155D-1,.610819696041396D-1,.366637909687335D-1, * .122236989606158D-1/ DATA T(128),T(129),T(130),T(131),T(132),T(133),T(134),T(135), * T(136),T(137),T(138),T(139),T(140),T(141),T(142)/ * .999956050018992D0,.999768437409263D0,.999430937466261D0, * .998943525843409D0,.998306266473006D0,.997519252756721D0, * .996582602023382D0,.995496454481096D0,.994260972922410D0, * .992876342608822D0,.991342771207583D0,.989660488745065D0, * .987829747564861D0,.985850822286126D0,.983724009760315D0/ DATA T(143),T(144),T(145),T(146),T(147),T(148),T(149),T(150), * T(151),T(152),T(153),T(154),T(155),T(156),T(157)/ * .981449629025464D0,.979028021257622D0,.976459549719234D0, * .973744599704370D0,.970883578480743D0,.967876915228489D0, * .964725060975706D0,.961428488530732D0,.957987692411178D0, * .954403188769716D0,.950675515316628D0,.946805231239127D0, * .942792917117462D0,.938639174837814D0, .934344627502003D0/ DATA T(158),T(159),T(160),T(161),T(162),T(163),T(164),T(165), * T(166),T(167),T(168),T(169),T(170),T(171),T(172)/ * .929909919334006D0,.925335715583316D0,.920622702425146D0, * .915771586857490D0,.910783096595065D0,.905657979960145D0, * .900397005770304D0,.895000963223085D0,.889470661777611D0, * .883806931033158D0,.878010620604707D0,.872082599995488D0, * .866023758466555D0,.859835004903376D0,.853517267679503D0/ DATA T(173),T(174),T(175),T(176),T(177),T(178),T(179),T(180), * T(181),T(182),T(183),T(184),T(185),T(186),T(187)/ * .847071494517296D0,.840498652345763D0,.833799727155505D0, * .826975723850813D0,.820027666098917D0,.812956596176432D0, * .805763574812999D0,.798449681032171D0,.791016011989546D0, * .783463682808184D0,.775793826411326D0,.768007593352446D0, * .760106151642655D0,.752090686575492D0,.743962400549112D0/ DATA T(188),T(189),T(190),T(191),T(192),T(193),T(194),T(195), * T(196),T(197),T(198),T(199),T(200),T(201),T(202)/ * .735722512885918D0,.727372259649652D0,.718912893459971D0, * .710345683304543D0,.701671914348685D0,.692892887742577D0, * .684009920426076D0,.675024344931163D0,.665937509182049D0, * .656750776292973D0,.647465524363725D0,.638083146272911D0, * .628605049469015D0,.619032655759261D0,.609367401096334D0/ DATA T(203),T(204),T(205),T(206),T(207),T(208),T(209),T(210), * T(211),T(212),T(213),T(214),T(215),T(216),T(217)/ * .599610735362968D0,.589764122154454D0,.579829038559083D0, * .569806974936569D0,.559699434694481D0,.549507934062719D0, * .539234001866059D0,.528879179294822D0,.518445019673674D0, * .507933088228616D0,.497344961852181D0,.486682228866890D0, * .475946488786983D0,.465139352078479D0,.454262439917590D0/ DATA T(218),T(219),T(220),T(221),T(222),T(223),T(224),T(225), * T(226),T(227),T(228),T(229),T(230),T(231),T(232)/ * .443317383947527D0,.432305826033741D0,.421229418017624D0, * .410089821468717D0,.398888707435459D0,.387627756194516D0, * .376308656998716D0,.364933107823654D0,.353502815112970D0, * .342019493522372D0,.330484865662417D0,.318900661840106D0, * .307268619799319D0,.295590484460136D0,.283868007657082D0/ DATA T(233),T(234),T(235),T(236),T(237),T(238),T(239),T(240), * T(241),T(242),T(243),T(244),T(245),T(246),T(247)/ * .272102947876337D0,.260297069991943D0,.248452145001057D0, * .236569949758284D0,.224652266709132D0,.212700883622626D0, * .200717593323127D0,.188704193421389D0,.176662486044902D0, * .164594277567554D0,.152501378338656D0,.140385602411376D0, * .128248767270607D0,.116092693560333D0,.103919204810509D0/ DATA T(248),T(249),T(250),T(251),T(252),T(253),T(254),T(255)/ * .917301271635196D-1,.795272891002330D-1,.673125211657164D-1, * .550876556946340D-1,.428545265363791D-1,.306149687799790D-1, * .183708184788137D-1,.612391237518953D-2/ DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9),W(10),W(11), * W(12),W(13),W(14),W(15)/1.0D0,.347854845137454D0, * .652145154862546D0,.101228536290376D0,.222381034453374D0, * .313706645877887D0,.362683783378362D0,.271524594117541D-1, * .622535239386479D-1,.951585116824928D-1,.124628971255534D0, * .149595988816577D0,.169156519395003D0,.182603415044924D0, * .189450610455068D0/ DATA W(16),W(17),W(18),W(19),W(20),W(21),W(22),W(23),W(24),W(25), * W(26),W(27),W(28),W(29),W(30),W(31)/.701861000947010D-2, * .162743947309057D-1,.253920653092621D-1,.342738629130214D-1, * .428358980222267D-1,.509980592623762D-1,.586840934785355D-1, * .658222227763618D-1,.723457941088485D-1,.781938957870703D-1, * .833119242269468D-1,.876520930044038D-1,.911738786957639D-1, * .938443990808046D-1,.956387200792749D-1,.965400885147278D-1/ DATA W(32),W(33),W(34),W(35),W(36),W(37),W(38),W(39),W(40),W(41), * W(42),W(43),W(44),W(45),W(46),W(47)/.178328072169643D-2, * .414703326056247D-2,.650445796897836D-2,.884675982636395D-2, * .111681394601311D-1,.134630478967186D-1,.157260304760247D-1, * .179517157756973D-1,.201348231535302D-1,.222701738083833D-1, * .243527025687109D-1,.263774697150547D-1,.283396726142595D-1, * .302346570724025D-1,.320579283548516D-1,.338051618371416D-1/ DATA W(48),W(49),W(50),W(51),W(52),W(53),W(54),W(55),W(56),W(57), * W(58),W(59),W(60),W(61),W(62),W(63)/ .354722132568824D-1, * .370551285402400D-1,.385501531786156D-1,.399537411327203D-1, * .412625632426235D-1,.424735151236536D-1,.435837245293235D-1, * .445905581637566D-1,.454916279274181D-1, .462847965813144D-1, * .469681828162100D-1,.475401657148303D-1,.479993885964583D-1, * .483447622348030D-1,.485754674415034D-1,.486909570091397D-1/ DATA W(64),W(65),W(66),W(67),W(68),W(69),W(70),W(71),W(72),W(73), * W(74),W(75),W(76),W(77),W(78),W(79)/ .449380960292090D-3, * .104581267934035D-2,.164250301866903D-2,.223828843096262D-2, * .283275147145799D-2,.342552604091022D-2,.401625498373864D-2, * .460458425670296D-2,.519016183267633D-2,.577263754286570D-2, * .635166316170719D-2,.692689256689881D-2,.749798192563473D-2, * .806458989048606D-2,.862637779861675D-2,.918300987166087D-2/ DATA W(80),W(81),W(82),W(83),W(84),W(85),W(86),W(87),W(88), * W(89),W(90),W(91),W(92),W(93),W(94),W(95)/.973415341500681D-2, * .102794790158322D-1,.108186607395031D-1,.113513763240804D-1, * .118773073727403D-1,.123961395439509D-1,.129075627392673D-1, * .134112712886163D-1,.139069641329520D-1,.143943450041668D-1, * .148731226021473D-1,.153430107688651D-1,.158037286593993D-1, * .162550009097852D-1,.166965578015892D-1,.171281354231114D-1/ DATA W(96),W(97),W(98),W(99),W(100),W(101),W(102),W(103),W(104), * W(105),W(106),W(107),W(108),W(109),W(110),W(111),W(112)/ * .175494758271177D-1,.179603271850087D-1,.183604439373313D-1, * .187495869405447D-1,.191275236099509D-1,.194940280587066D-1, * .198488812328309D-1,.201918710421300D-1,.205227924869601D-1, * .208414477807511D-1,.211476464682213D-1,.214412055392085D-1, * .217219495380521D-1,.219897106684605D-1,.222443288937998D-1, * .224856520327450D-1,.227135358502365D-1/ DATA W(113),W(114),W(115),W(116),W(117),W(118),W(119),W(120), * W(121),W(122),W(123),W(124),W(125),W(126),W(127)/ * .229278441436868D-1,.231284488243870D-1,.233152299940628D-1, * .234880760165359D-1,.236468835844476D-1,.237915577810034D-1, * .239220121367035D-1,.240381686810241D-1,.241399579890193D-1, * .242273192228152D-1,.243002001679719D-1,.243585572646906D-1, * .244023556338496D-1,.244315690978500D-1,.244461801962625D-1/ DATA W(128),W(129),W(130),W(131),W(132),W(133),W(134),W(135), * W(136),W(137),W(138),W(139),W(140),W(141),W(142)/ * .112789017822272D-3,.262534944296446D-3,.412463254426176D-3, * .562348954031410D-3,.712154163473321D-3,.861853701420089D-3, * .101142439320844D-2,.116084355756772D-2,.131008868190250D-2, * .145913733331073D-2,.160796713074933D-2,.175655573633073D-2, * .190488085349972D-2,.205292022796614D-2,.220065164983991D-2/ DATA W(143),W(144),W(145),W(146),W(147),W(148),W(149),W(150), * W(151),W(152),W(153),W(154),W(155),W(156),W(157)/ * .234805295632731D-2,.249510203470371D-2,.264177682542749D-2, * .278805532532771D-2,.293391559082972D-2,.307933574119934D-2, * .322429396179420D-2,.336876850731555D-2,.351273770505631D-2, * .365617995814250D-2,.379907374876626D-2,.394139764140883D-2, * .408313028605267D-2,.422425042138154D-2,.436473687796806D-2/ DATA W(158),W(159),W(160),W(161),W(162),W(163),W(164),W(165), * W(166),W(167),W(168),W(169),W(170),W(171),W(172)/ * .450456858144790D-2,.464372455568006D-2,.478218392589269D-2, * .491992592181387D-2,.505692988078684D-2,.519317525086928D-2, * .532864159391593D-2,.546330858864431D-2,.559715603368291D-2, * .573016385060144D-2,.586231208692265D-2,.599358091911534D-2, * .612395065556793D-2,.625340173954240D-2,.638191475210788D-2/ DATA W(173),W(174),W(175),W(176),W(177),W(178),W(179),W(180), * W(181),W(182),W(183),W(184),W(185),W(186),W(187)/ * .650947041505366D-2,.663604959378107D-2,.676163330017380D-2, * .688620269544632D-2,.700973909296982D-2,.713222396107539D-2, * .725363892583391D-2,.737396577381235D-2,.749318645480588D-2, * .761128308454566D-2,.772823794738156D-2,.784403349893971D-2, * .795865236875435D-2,.807207736287350D-2,.818429146643827D-2/ DATA W(188),W(189),W(190),W(191),W(192),W(193),W(194),W(195), * W(196),W(197),W(198),W(199),W(200),W(201),W(202)/ * .829527784623523D-2,.840501985322154D-2,.851350102502249D-2, * .862070508840101D-2,.872661596169881D-2,.883121775724875D-2, * .893449478375821D-2,.903643154866287D-2,.913701276045081D-2, * .923622333095630D-2,.933404837762327D-2,.943047322573775D-2, * .952548341062928D-2,.961906467984073D-2,.971120299526628D-2/ DATA W(203),W(204),W(205),W(206),W(207),W(208),W(209),W(210), * W(211),W(212),W(213),W(214),W(215),W(216),W(217)/ * .980188453525733D-2,.989109569669583D-2,.997882309703491D-2, * .100650535763064D-1,.101497741990949D-1,.102329722564782D-1, * .103146352679340D-1,.103947509832117D-1,.104733073841704D-1, * .105502926865815D-1,.106256953418966D-1,.106995040389798D-1, * .107717077058046D-1,.108422955111148D-1,.109112568660490D-1/ DATA W(218),W(219),W(220),W(221),W(222),W(223),W(224),W(225), * W(226),W(227),W(228),W(229),W(230),W(231),W(232)/ * .109785814257296D-1,.110442590908139D-1,.111082800090098D-1, * .111706345765534D-1,.112313134396497D-1,.112903074958755D-1, * .113476078955455D-1,.114032060430392D-1,.114570935980906D-1, * .115092624770395D-1,.115597048540436D-1,.116084131622531D-1, * .116553800949452D-1,.117005986066207D-1,.117440619140606D-1/ DATA W(233),W(234),W(235),W(236),W(237),W(238),W(239),W(240), * W(241),W(242),W(243),W(244),W(245),W(246),W(247)/ * .117857634973434D-1,.118256971008240D-1,.118638567340711D-1, * .119002366727665D-1,.119348314595636D-1,.119676359049059D-1, $ .119986450878058D-1,.120278543565826D-1,.120552593295601D-1, * .120808558957245D-1,.121046402153405D-1,.121266087205273D-1, * .121467581157945D-1,.121650853785355D-1,.121815877594818D-1/ DATA W(248),W(249),W(250),W(251),W(252),W(253),W(254),W(255)/ * .121962627831147D-1,.122091082480372D-1,.122201222273040D-1, * .122293030687103D-1,.122366493950402D-1,.122421601042728D-1, * .122458343697479D-1,.122476716402898D-1/ SCALE=(B-A)/2.0D0 M=N/2 NPLACE=M-1 DO 1 I=1,M S=T(NPLACE+I) R=W(NPLACE+I)*SCALE TT(I)=(A*(1.0D0+S)+(1.0D0-S)*B)/2.0D0 TT( N+1-I)=(A*(1.0D0-S)+B*(1.0D0+S))/2.0D0 WW(I)=R 1 WW(N+1-I)=R RETURN END SUBROUTINE LINSYS(MAT,B,N,MD,IER) C ----------------- C C THIS ROUTINE SOLVES A SYSTEM OF LINEAR EQUATIONS C A*X = B C THE METHOD USED IS GAUSSIAN ELIMINATION WITH C PARTIAL PIVOTING. C C INPUT: C THE COEFFICIENT MATRIX A IS STORED IN THE ARRAY MAT. C THE RIGHT SIDE CONSTANTS ARE IN THE ARRAY B. C THE ORDER OF THE LINEAR SYSTEM IS N. C THE VARIABLE MD IS THE NUMBER OF ROWS THAT MAT C IS DIMENSIONED AS HAVING IN THE CALLING PROGRAM. C C OUTPUT: C THE ARRAY B CONTAINS THE SOLUTION X. C MAT CONTAINS THE UPPER TRIANGULAR MATRIX U C OBTAINED BY ELIMINATION. THE ROW MULTIPLIERS C USED IN THE ELIMINATION ARE STORED IN THE C LOWER TRIANGULAR PART OF MAT. C IER=0 MEANS THE MATRIX A WAS COMPUTATIONALLY C NONSINGULAR, AND THE GAUSSIAN ELIMINATION C WAS COMPLETED SATISFACTORILY. C IER=1 MEANS THAT THE MATRIX A WAS C COMPUTATIONALLY SINGULAR. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER PIVOT DOUBLE PRECISION MULT,MAT DIMENSION MAT(MD,*),B(*) C C BEGIN ELIMINATION STEPS. DO 40 K=1,N-1 C CHOOSE PIVOT ROW. PIVOT=K AMAX=ABS(MAT(K,K)) DO 10 I=K+1,N ABSA=ABS(MAT(I,K)) IF(ABSA .GT. AMAX) THEN PIVOT=I AMAX=ABSA END IF 10 CONTINUE C IF(AMAX .EQ. 0.0D0) THEN C COEFFICIENT MATRIX IS SINGULAR. IER=1 RETURN END IF C IF(PIVOT .NE. K) THEN C SWITCH ROWS K AND PIVOT. TEMP=B(K) B(K)=B(PIVOT) B(PIVOT)=TEMP DO 20 J=K,N TEMP=MAT(K,J) MAT(K,J)=MAT(PIVOT,J) 20 MAT(PIVOT,J)=TEMP END IF C C PERFORM STEP #K OF ELIMINATION. DO 30 I=K+1,N MULT=MAT(I,K)/MAT(K,K) MAT(I,K)=MULT B(I)=B(I)-MULT*B(K) DO 30 J=K+1,N 30 MAT(I,J)=MAT(I,J)-MULT*MAT(K,J) 40 CONTINUE C IF(MAT(N,N) .EQ. 0.0D0) THEN C COEFFICIENT MATRIX IS SINGULAR. IER=1 RETURN END IF C C SOLVE FOR SOLUTION X USING BACK SUBSTITUTION. DO 60 I=N,1,-1 SUM=0.0D0 DO 50 J=I+1,N 50 SUM=SUM+MAT(I,J)*B(J) 60 B(I)=(B(I)-SUM)/MAT(I,I) IER=0 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