C ALGORITHM 629 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 2, C JUN., 1985, P. 85-96. 3 2 1 1 100 5 16 8 COEFFICIENTS.DECIMAL COEFFICIENTS.FORMATFREE C 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=9,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/7/,FILEDC/8/,FILEFF/9/ C C ****************************************** C C OBTAIN THE NAME OF THE INPUT FILE. 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) CALL EXIT 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 QUAN 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 QUADRATUR 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 UPPE 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 T 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, 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 F C DEGREES .LE. N:=DEGREE, AS INPUT. ON OUTPUT, THERE WILL BE N TABL 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 C 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 'R1MACH', WHICH IS ONE OF THE SUBPROGRAMS GIVEN IN THE C 'GNRT' PACKAGE. THE ROUTINE 'R1MACH' SHOULD BE RESET C APPROPRIATELY FOR THE USER'S COMPUTER. C INTEGER DEGREE,FILEDC,FILEFF,FILEIN CHARACTER NAMEDC*50,NAMEFF*50,NAMEIN*50 PARAMETER(MAXDEG=9,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/7/,FILEDC/8/,FILEFF/9/ C C ****************************************** C C OBTAIN THE NAME OF THE INPUT FILE. READ(1,'(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)=', * 1PE17.10,2E21.10,/,' SURFACE PARAMETER ALPHA=',E17.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) CALL EXIT 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 QUAN 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 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.0 IF(IFLAG .EQ. 0) GO TO 100 DXR = 0.0 DYR = 0.0 DZR = 0.0 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.0*UZ*UZ - 1.0 TEMP = SQRT(ALPHA + C2T*C2T) R = SQRT(C2T + TEMP) IF(IFLAG .EQ. 0) GO TO 100 DXR = 0.0 DYR = 0.0 DZR=2.0*UZ*(1.0+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.0+ALPHA*(UZ+1.0)**2 R=2.0-1.0/TEMP IF(IFLAG .EQ. 0) GO TO 100 DXR=0.0 DYR=0.0 DZR=2.0*ALPHA*(UZ+1.0)/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 QUADRATUR 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 UPPE 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 T C FOLLOWING SUBROUTINE 'GNRT2'. C 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*R1MACH(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 REAL 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.0,P5=0.5,ONE=1.0,TWO=2.0,FOUR=4.0) 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, 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(1PE20.10,5E20.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 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.0,ONE=1.0,TWO=2.0) 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 REAL 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.0) 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 PARAMETER(MAXDEG=20,LFORD=2*MAXDEG) DIMENSION VALUES(*),FACT(0:LFORD) PARAMETER(ONE=1.0,TWO=2.0,FOUR=4.0) 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 F C DEGREES .LE. N:=DEGREE, AS INPUT. ON OUTPUT, THERE WILL BE N TABL C FOR GALERKIN COEFFICIENTS OF DEGREES .LE. L, WITH 1 .LE. L .LE. N. C INTEGER DEGREE,DEGSQ,FILEFF REAL 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 PARAMETER(ONE=1.0,TWO=2.0) 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 INTEGER ARG(5) DIMENSION A1(4),A2(4),Q(68),Q0300(18),Q0603(16),Q0806(15), * Q1008(19) PARAMETER(C40D3=40.0/3.0) PARAMETER(ZERO=0.0,ONE=1.0,TWO=2.0,FOUR=4.0) PARAMETER(DIVSOR=7.0298584066096E-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.1307592530313668E-01, 0.8580903946529252E-02, * 0.1969565496511163E-02, -0.6785908047534174E-04, * 0.5157989084218512E-05, -0.3254283578286669E-06, * 0.3009165432294816E-07, -0.2853960792992456E-08, * 0.3125316663674964E-09, -0.369071915636499E-10, * 0.47226104267733E-11, -0.6452090127750E-12, * 0.935240772542E-13, -0.142880027692E-13, * 0.22888743282E-14, -0.3827976860E-15, * 0.665889416E-16, -0.120097537E-16, * 0.22395647E-17/ DATA Q0806/0.7530091699129213E-01, 0.2215896190278894E-01, * 0.1517662771681320E-02, -0.1374204580631322E-04, * 0.2501181491115358E-05, -0.2491206397236787E-07, * 0.5430948732810670E-08, -0.1406070367942514E-09, * 0.160397857131033E-10, -0.7706914621139E-12, * 0.656131517151E-13, -0.43257167630E-14, * 0.3459964512E-15, -0.264283638E-16, * 0.21607480E-17/ DATA Q0603/0.2285305746758029E+00, 0.5670149715650661E-01, * 0.3881611773394649E-02, 0.1483758828946990E-03, * 0.1986416462810431E-04, 0.1166284710859293E-05, * 0.1048168134503124E-06, 0.6572636994171403E-08, * 0.5344588684897204E-09, 0.335115172128537E-10, * 0.26225503158527E-11, 0.1606252080762E-12, * 0.124685606769E-13, 0.73251000E-15, * 0.579829277E-16, 0.31817034E-17/ DATA Q0300/0.5404466499651320E+00, 0.1034761954005156E+00, * 0.9090551216566596E-02, 0.9130257654553327E-03, * 0.1026685176623610E-03, 0.1035244509501565E-04, * 0.1034925154934048E-05, 0.1002050044392230E-06, * 0.9552837592432324E-08, 0.8941334072231475E-09, * 0.825727197051832E-10, 0.75255600668576E-11, * 0.6783483380205E-12, 0.605164239084E-13, * 0.53497536298E-14, 0.4689365198E-15, * 0.407909118E-16, 0.35230149E-17/ DATA A1/20.0,20.0,C40D3,C40D3/, * A2/18.0,14.0,6.0,2.0/, * ARG/1,20,35,51,69/ DATA P1/-.3/,P2/-.6/,P3/-.8/ C C SET THE MACHINE DEPENDENT CONSTANTS. UFLOW=2*R1MACH(1) EOFLOW=LOG(R1MACH(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 DIMENSION Z(109),W(109),ZZ(N),WW(N) DATA ONE/1.0/,TWO/2.0/ DATA Z(1),Z(2),Z(3),Z(4),Z(5),Z(6),Z(7),Z(8),Z(9),Z(10)/ *.577350269189626E0,.774596669241483E0,0.0E0,.861136311594053E0, *.339981043584856E0,.906179845938664E0,.538469310105683E0,0.0E0, *.932469514203152E0,.661209386466265E0/ DATA Z(11),Z(12),Z(13),Z(14),Z(15),Z(16),Z(17),Z(18),Z(19),Z(20)/ *.238619186083197E0,.949107912342759E0,.741531185599394E0, *.405845151377397E0,0.0E0,.960289856497536E0,.796666477413627E0, *.525532409916329E0,.183434642495650E0,.968160239507626E0/ DATA Z(21),Z(22),Z(23),Z(24),Z(25),Z(26),Z(27),Z(28),Z(29)/ *.836031107326636E0,.613371432700590E0,.324253423403809E0, *0.0E0,.973906528517172E0,.865063366688985E0,.679409568299024E0, *.433395394129247E0,.148874338981631E0/ DATA Z(30),Z(31),Z(32),Z(33),Z(34),Z(35),Z(36),Z(37),Z(38)/ *.978228658146057E0,.887062599768095E0,.730152005574049E0, *.519096129206812E0,.269543155952345E0,0.0E0, *.981560634246719E0,.904117256370475E0,.769902674194305E0/ DATA Z(39),Z(40),Z(41),Z(42),Z(43),Z(44),Z(45),Z(46),Z(47)/ *.587317954286617E0,.367831498998180E0,.125233408511469E0, *.984183054718588E0,.917598399222978E0,.801578090733310E0, *.642349339440340E0,.448492751036447E0,.230458315955135E0/ DATA Z(48),Z(49),Z(50),Z(51),Z(52),Z(53),Z(54),Z(55),Z(56)/ *0.0E0,.986283808696812E0,.928434883663574E0, *.827201315069765E0,.687292904811685E0,.515248636358154E0, *.319112368927890E0,.108054948707344E0,.987992518020485E0/ DATA Z(57),Z(58),Z(59),Z(60),Z(61),Z(62),Z(63),Z(64),Z(65)/ *.937273392400706E0,.848206583410427E0,.724417731360170E0, *.570972172608539E0,.394151347077563E0,.201194093997435E0, *0.0E0,.989400934991650E0,.944575023073233E0/ DATA Z(66),Z(67),Z(68),Z(69),Z(70),Z(71),Z(72),Z(73),Z(74)/ *.865631202387832E0,.755404408355003E0,.617876244402644E0, *.458016777657227E0,.281603550779259E0,.0950125098376374E0, *.990575475314417E0,.950675521768768E0,.880239153726986E0/ DATA Z(75),Z(76),Z(77),Z(78),Z(79),Z(80),Z(81),Z(82),Z(83)/ *.781514003896801E0,.657671159216691E0,.512690537086477E0, *.351231763453876E0,.178484181495848E0,0.0E0, *.991565168420931E0,.955823949571398E0,.892602466497556E0/ DATA Z(84),Z(85),Z(86),Z(87),Z(88),Z(89),Z(90),Z(91),Z(92)/ *.803704958972523E0,.691687043060353E0,.559770831073948E0, *.411751161462843E0,.251886225691506E0,.0847750130417353E0, *.992406843843584E0,.960208152134830E0,.903155903614818E0/ DATA Z(93),Z(94),Z(95),Z(96),Z(97),Z(98),Z(99),Z(100),Z(101)/ *.822714656537143E0,.720966177335229E0,.600545304661681E0, *.464570741375961E0,.316564099963630E0,.160358645640225E0, *0.0E0,.993128599185095E0,.963971927277914E0/ DATA Z(102),Z(103),Z(104),Z(105),Z(106),Z(107),Z(108),Z(109)/ *.912234428251326E0,.839116971822219E0,.746331906460151E0, *.636053680726515E0,.510867001950827E0,.373706088715420E0, *.227785851141645E0,.0765265211334973E0/ DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9),W(10)/ *1.0E0,.555555555555556E0,.888888888888889E0,.347854845137454E0, *.652145154862546E0,.236926885056189E0,.478628670499366E0, *.568888888888889E0,.171324492379170E0,.360761573048139E0/ DATA W(11),W(12),W(13),W(14),W(15),W(16),W(17),W(18),W(19),W(20)/ *.467913934572691E0,.129484966168870E0,.279705391489277E0, *.381830050505119E0,.417959183673469E0,.101228536290376E0, *.222381034453374E0,.313706645877887E0,.362683783378362E0, *.0812743883615744E0/ DATA W(21),W(22),W(23),W(24),W(25),W(26),W(27),W(28),W(29)/ *.180648160694857E0,.260610696402935E0,.312347077040003E0, *.330239355001260E0,.0666713443086881E0,.149451349150581E0, *.219086362515982E0,.269266719309996E0,.295524224714753E0/ DATA W(30),W(31),W(32),W(33),W(34),W(35),W(36),W(37),W(38)/ *.0556685671161737E0,.125580369464905E0,.186290210927734E0, *.233193764591990E0,.262804544510247E0,.272925086777901E0, *.0471753363865118E0,.106939325995318E0,.160078328543346E0/ DATA W(39),W(40),W(41),W(42),W(43),W(44),W(45),W(46),W(47)/ *.203167426723066E0,.233492536538355E0,.249147045813403E0, *.0404840047653159E0,.0921214998377284E0,.138873510219787E0, *.178145980761946E0,.207816047536889E0,.226283180262897E0/ DATA W(48),W(49),W(50),W(51),W(52),W(53),W(54),W(55),W(56)/ *.232551553230874E0,.0351194603317519E0,.0801580871597602E0, *.121518570687903E0,.157203167158194E0,.185538397477938E0, *.205198463721296E0,.215263853463158E0,.0307532419961173E0/ DATA W(57),W(58),W(59),W(60),W(61),W(62),W(63),W(64),W(65)/ *.0703660474881081E0,.107159220467172E0,.139570677926154E0, *.166269205816994E0,.186161000015562E0,.198431485327112E0, *.202578241925561E0,.0271524594117541E0,.0622535239386478E0/ DATA W(66),W(67),W(68),W(69),W(70),W(71),W(72),W(73),W(74)/ *.0951585116824928E0,.124628971255534E0,.149595988816577E0, *.169156519395003E0,.182603415044924E0,.189450610455068E0, *.0241483028685479E0,.0554595293739872E0,.0850361483171792E0/ DATA W(75),W(76),W(77),W(78),W(79),W(80),W(81),W(82),W(83)/ *.111883847193404E0,.135136368468525E0,.154045761076810E0, *.168004102156450E0,.176562705366993E0,.179446470356207E0, *.0216160135264833E0,.0497145488949698E0,.0764257302548891E0/ DATA W(84),W(85),W(86),W(87),W(88),W(89),W(90),W(91),W(92)/ *.100942044106287E0,.122555206711478E0,.140642914670651E0, *.154684675126265E0,.164276483745833E0,.169142382963144E0, *.0194617882297265E0,.0448142267656996E0,.0690445427376412E0/ DATA W(93),W(94),W(95),W(96),W(97),W(98),W(99),W(100),W(101)/ *.0914900216224500E0,.111566645547334E0,.128753962539336E0, *.142606702173607E0,.152766042065860E0,.158968843393954E0, *.161054449848784E0,.0176140071391521E0,.0406014298003869E0/ DATA W(102),W(103),W(104),W(105),W(106),W(107),W(108),W(109)/ *.0626720483341091E0,.0832767415767047E0,.101930119817240E0, *.118194531961518E0,.131688638449177E0,.142096109318382E0, *.149172986472604E0,.152753387130726E0/ 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 C 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=FILTRM) 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,FILTRM,FILALT,FILOUT,DEGREE 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/5/,FILTRM/1/,FILALT/6/ 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(1,'(A1)') T IF(T .EQ. 'Y') THEN FILOUT=FILTRM 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 CALL EXIT END IF PRINT *, ' GIVE THE DEGREE OF THE APPROXIMATION AND THE' PRINT *, ' INTEGRATION PARAMETER.' 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) CALL EXIT 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,5X,'AUXILIARY BOUNDARY', * ' FUNCTION PARAMETER=',1PD17.10,/,' DEGREE OF ', * 'APPROXIMATION=',I2,5X,'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 QUAN 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 ROUT 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 FOR 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 C 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=FILTRM) 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 'R1MACH', WHICH IS ONE OF THE SUBPROGRAMS GIVEN IN THE C 'LAPLAC' PACKAGE. THE ROUTINE 'R1MACH' SHOULD BE RESET C APPROPRIATELY FOR THE USER'S COMPUTER. C INTEGER FILEIN,FILTRM,FILALT,FILOUT,DEGREE 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/5/,FILTRM/1/,FILALT/6/ 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') THEN FILOUT=FILTRM 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,)=', * 1PE17.10,2E20.10,/,' AUXILIARY SURFACE PARAMETER ALPHA = ', * E17.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 CALL EXIT END IF PRINT *, ' GIVE THE DEGREE OF THE APPROXIMATION AND THE' PRINT *, ' INTEGRATION PARAMETER.' 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) CALL EXIT 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,5X,'AUXILIARY BOUNDARY', * ' FUNCTION PARAMETER=',1PE17.10,/,' DEGREE OF ', * 'APPROXIMATION=',I2,5X,'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,1PE21.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,1PE21.10,E14.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 COMMON/BDYPAR/NUMBF,A DIMENSION Q(3) PARAMETER(ONE=1.0,TWO=2.0,THREE=3.0,FOUR=4.0, * FIVE=5.0) 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 QUAN 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 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.0 IF(IFLAG .EQ. 0) GO TO 100 DXR = 0.0 DYR = 0.0 DZR = 0.0 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.0*UZ*UZ - 1.0 TEMP = SQRT(ALPHA + C2T*C2T) R = SQRT(C2T + TEMP) IF(IFLAG .EQ. 0) GO TO 100 DXR = 0.0 DYR = 0.0 DZR=2.0*UZ*(1.0+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.0+ALPHA*(UZ+1.0)**2 R=2.0-1.0/TEMP IF(IFLAG .EQ. 0) GO TO 100 DXR=0.0 DYR=0.0 DZR=2.0*ALPHA*(UZ+1.0)/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 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*R1MACH(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 ROUT C REAL 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.0,ONE=1.0,TWO=2.0,FOUR=4.0) 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 REAL KERNEL DIMENSION P(3),Q(3),QCROSS(3) COMMON/MACHCN/U100,U100SQ C U100 IS 100 TIMES THE MACHINE UNIT ROUND. PARAMETER(ZERO=0.0) 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 PARAMETER(MAXDEG=20,LFORD=2*MAXDEG) DIMENSION VALUES(*),FACT(0:LFORD) PARAMETER(ONE=1.0,TWO=2.0,FOUR=4.0) 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 DIMENSION Z(109),W(109),ZZ(N),WW(N) DATA ONE/1.0/,TWO/2.0/ DATA Z(1),Z(2),Z(3),Z(4),Z(5),Z(6),Z(7),Z(8),Z(9),Z(10)/ *.577350269189626E0,.774596669241483E0,0.0E0,.861136311594053E0, *.339981043584856E0,.906179845938664E0,.538469310105683E0,0.0E0, *.932469514203152E0,.661209386466265E0/ DATA Z(11),Z(12),Z(13),Z(14),Z(15),Z(16),Z(17),Z(18),Z(19),Z(20)/ *.238619186083197E0,.949107912342759E0,.741531185599394E0, *.405845151377397E0,0.0E0,.960289856497536E0,.796666477413627E0, *.525532409916329E0,.183434642495650E0,.968160239507626E0/ DATA Z(21),Z(22),Z(23),Z(24),Z(25),Z(26),Z(27),Z(28),Z(29)/ *.836031107326636E0,.613371432700590E0,.324253423403809E0, *0.0E0,.973906528517172E0,.865063366688985E0,.679409568299024E0, *.433395394129247E0,.148874338981631E0/ DATA Z(30),Z(31),Z(32),Z(33),Z(34),Z(35),Z(36),Z(37),Z(38)/ *.978228658146057E0,.887062599768095E0,.730152005574049E0, *.519096129206812E0,.269543155952345E0,0.0E0, *.981560634246719E0,.904117256370475E0,.769902674194305E0/ DATA Z(39),Z(40),Z(41),Z(42),Z(43),Z(44),Z(45),Z(46),Z(47)/ *.587317954286617E0,.367831498998180E0,.125233408511469E0, *.984183054718588E0,.917598399222978E0,.801578090733310E0, *.642349339440340E0,.448492751036447E0,.230458315955135E0/ DATA Z(48),Z(49),Z(50),Z(51),Z(52),Z(53),Z(54),Z(55),Z(56)/ *0.0E0,.986283808696812E0,.928434883663574E0, *.827201315069765E0,.687292904811685E0,.515248636358154E0, *.319112368927890E0,.108054948707344E0,.987992518020485E0/ DATA Z(57),Z(58),Z(59),Z(60),Z(61),Z(62),Z(63),Z(64),Z(65)/ *.937273392400706E0,.848206583410427E0,.724417731360170E0, *.570972172608539E0,.394151347077563E0,.201194093997435E0, *0.0E0,.989400934991650E0,.944575023073233E0/ DATA Z(66),Z(67),Z(68),Z(69),Z(70),Z(71),Z(72),Z(73),Z(74)/ *.865631202387832E0,.755404408355003E0,.617876244402644E0, *.458016777657227E0,.281603550779259E0,.0950125098376374E0, *.990575475314417E0,.950675521768768E0,.880239153726986E0/ DATA Z(75),Z(76),Z(77),Z(78),Z(79),Z(80),Z(81),Z(82),Z(83)/ *.781514003896801E0,.657671159216691E0,.512690537086477E0, *.351231763453876E0,.178484181495848E0,0.0E0, *.991565168420931E0,.955823949571398E0,.892602466497556E0/ DATA Z(84),Z(85),Z(86),Z(87),Z(88),Z(89),Z(90),Z(91),Z(92)/ *.803704958972523E0,.691687043060353E0,.559770831073948E0, *.411751161462843E0,.251886225691506E0,.0847750130417353E0, *.992406843843584E0,.960208152134830E0,.903155903614818E0/ DATA Z(93),Z(94),Z(95),Z(96),Z(97),Z(98),Z(99),Z(100),Z(101)/ *.822714656537143E0,.720966177335229E0,.600545304661681E0, *.464570741375961E0,.316564099963630E0,.160358645640225E0, *0.0E0,.993128599185095E0,.963971927277914E0/ DATA Z(102),Z(103),Z(104),Z(105),Z(106),Z(107),Z(108),Z(109)/ *.912234428251326E0,.839116971822219E0,.746331906460151E0, *.636053680726515E0,.510867001950827E0,.373706088715420E0, *.227785851141645E0,.0765265211334973E0/ DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9),W(10)/ *1.0E0,.555555555555556E0,.888888888888889E0,.347854845137454E0, *.652145154862546E0,.236926885056189E0,.478628670499366E0, *.568888888888889E0,.171324492379170E0,.360761573048139E0/ DATA W(11),W(12),W(13),W(14),W(15),W(16),W(17),W(18),W(19),W(20)/ *.467913934572691E0,.129484966168870E0,.279705391489277E0, *.381830050505119E0,.417959183673469E0,.101228536290376E0, *.222381034453374E0,.313706645877887E0,.362683783378362E0, *.0812743883615744E0/ DATA W(21),W(22),W(23),W(24),W(25),W(26),W(27),W(28),W(29)/ *.180648160694857E0,.260610696402935E0,.312347077040003E0, *.330239355001260E0,.0666713443086881E0,.149451349150581E0, *.219086362515982E0,.269266719309996E0,.295524224714753E0/ DATA W(30),W(31),W(32),W(33),W(34),W(35),W(36),W(37),W(38)/ *.0556685671161737E0,.125580369464905E0,.186290210927734E0, *.233193764591990E0,.262804544510247E0,.272925086777901E0, *.0471753363865118E0,.106939325995318E0,.160078328543346E0/ DATA W(39),W(40),W(41),W(42),W(43),W(44),W(45),W(46),W(47)/ *.203167426723066E0,.233492536538355E0,.249147045813403E0, *.0404840047653159E0,.0921214998377284E0,.138873510219787E0, *.178145980761946E0,.207816047536889E0,.226283180262897E0/ DATA W(48),W(49),W(50),W(51),W(52),W(53),W(54),W(55),W(56)/ *.232551553230874E0,.0351194603317519E0,.0801580871597602E0, *.121518570687903E0,.157203167158194E0,.185538397477938E0, *.205198463721296E0,.215263853463158E0,.0307532419961173E0/ DATA W(57),W(58),W(59),W(60),W(61),W(62),W(63),W(64),W(65)/ *.0703660474881081E0,.107159220467172E0,.139570677926154E0, *.166269205816994E0,.186161000015562E0,.198431485327112E0, *.202578241925561E0,.0271524594117541E0,.0622535239386478E0/ DATA W(66),W(67),W(68),W(69),W(70),W(71),W(72),W(73),W(74)/ *.0951585116824928E0,.124628971255534E0,.149595988816577E0, *.169156519395003E0,.182603415044924E0,.189450610455068E0, *.0241483028685479E0,.0554595293739872E0,.0850361483171792E0/ DATA W(75),W(76),W(77),W(78),W(79),W(80),W(81),W(82),W(83)/ *.111883847193404E0,.135136368468525E0,.154045761076810E0, *.168004102156450E0,.176562705366993E0,.179446470356207E0, *.0216160135264833E0,.0497145488949698E0,.0764257302548891E0/ DATA W(84),W(85),W(86),W(87),W(88),W(89),W(90),W(91),W(92)/ *.100942044106287E0,.122555206711478E0,.140642914670651E0, *.154684675126265E0,.164276483745833E0,.169142382963144E0, *.0194617882297265E0,.0448142267656996E0,.0690445427376412E0/ DATA W(93),W(94),W(95),W(96),W(97),W(98),W(99),W(100),W(101)/ *.0914900216224500E0,.111566645547334E0,.128753962539336E0, *.142606702173607E0,.152766042065860E0,.158968843393954E0, *.161054449848784E0,.0176140071391521E0,.0406014298003869E0/ DATA W(102),W(103),W(104),W(105),W(106),W(107),W(108),W(109)/ *.0626720483341091E0,.0832767415767047E0,.101930119817240E0, *.118194531961518E0,.131688638449177E0,.142096109318382E0, *.149172986472604E0,.152753387130726E0/ 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 FOR 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 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)/.577350269189626E0, * .861136311594053E0,.339981043584856E0,.960289856497536E0, * .796666477413627E0,.525532409916329E0,.183434642495650E0, * .989400934991650E0,.944575023073233E0,.865631202387832E0, * .755404408355003E0,.617876244402644E0,.458016777657227E0, * .281603550779259E0,.950125098376374E-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)/.997263861849482E0, * .985611511545268E0,.964762255587506E0,.934906075937740E0, * .896321155766052E0,.849367613732570E0,.794483795967942E0, * .732182118740290E0,.663044266930215E0,.587715757240762E0, * .506899908932229E0,.421351276130635E0,.331868602282128E0, * .239287362252137E0,.144471961582796E0,.483076656877383E-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)/.999305041735772E0, * .996340116771955E0,.991013371476744E0,.983336253884626E0, * .973326827789911E0,.961008799652054E0,.946411374858403E0, * .929569172131940E0,.910522137078503E0,.889315445995114E0, * .865999398154093E0,.840629296252580E0,.813265315122798E0, * .783972358943341E0,.752819907260532E0,.719881850171611E0/ 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)/.685236313054233E0, * .648965471254657E0,.611155355172393E0,.571895646202634E0, * .531279464019895E0,.489403145707053E0,.446366017253464E0, * .402270157963992E0,.357220158337668E0,.311322871990211E0, * .264687162208767E0,.217423643740007E0,.169644420423993E0, * .121462819296121E0,.729931217877990E-1,.243502926634244E-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)/.999824887947132E0, * .999077459977376E0,.997733248625514E0,.995792758534981E0, * .993257112900213E0,.990127818491734E0,.986406742724586E0, * .982096108435719E0,.977198491463907E0,.971716818747137E0, * .965654366431965E0,.959014757853700E0,.951801961341264E0, * .944020287830220E0,.935674388277916E0,.926769250878948E0/ 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)/ * .917310198080961E0,.907302883401757E0,.896753288049158E0, * .885667717345397E0,.874052796958032E0,.861915468939548E0, * .849262987577969E0,.836102915060907E0,.822443116955644E0, * .808291757507914E0,.793657294762193E0,.778548475506412E0, * .762974330044095E0,.746944166797062E0,.730467566741909E0, * .713554377683587E0,.696214708369514E0,.678458922447719E0/ 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)/ * .660297632272646E0,.641741692562308E0,.622802193910585E0, * .603490456158549E0,.583818021628763E0,.563796648226618E0, * .543438302412810E0,.522755152051175E0,.501759559136144E0, * .480464072404172E0,.458881419833552E0,.437024501037104E0, * .414906379552275E0,.392540275033267E0,.369939555349859E0, * .347117728597636E0,.324088435024413E0/ 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)/.300865438877677E0, * .277462620177904E0,.253893966422694E0,.230173564226660E0, * .206315590902079E0,.182334305985337E0,.158244042714225E0, * .134059199461188E0,.109794231127644E0, * .854636405045155E-1,.610819696041396E-1,.366637909687335E-1, * .122236989606158E-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)/ * .999956050018992E0,.999768437409263E0,.999430937466261E0, * .998943525843409E0,.998306266473006E0,.997519252756721E0, * .996582602023382E0,.995496454481096E0,.994260972922410E0, * .992876342608822E0,.991342771207583E0,.989660488745065E0, * .987829747564861E0,.985850822286126E0,.983724009760315E0/ 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)/ * .981449629025464E0,.979028021257622E0,.976459549719234E0, * .973744599704370E0,.970883578480743E0,.967876915228489E0, * .964725060975706E0,.961428488530732E0,.957987692411178E0, * .954403188769716E0,.950675515316628E0,.946805231239127E0, * .942792917117462E0,.938639174837814E0, .934344627502003E0/ 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)/ * .929909919334006E0,.925335715583316E0,.920622702425146E0, * .915771586857490E0,.910783096595065E0,.905657979960145E0, * .900397005770304E0,.895000963223085E0,.889470661777611E0, * .883806931033158E0,.878010620604707E0,.872082599995488E0, * .866023758466555E0,.859835004903376E0,.853517267679503E0/ 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)/ * .847071494517296E0,.840498652345763E0,.833799727155505E0, * .826975723850813E0,.820027666098917E0,.812956596176432E0, * .805763574812999E0,.798449681032171E0,.791016011989546E0, * .783463682808184E0,.775793826411326E0,.768007593352446E0, * .760106151642655E0,.752090686575492E0,.743962400549112E0/ 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)/ * .735722512885918E0,.727372259649652E0,.718912893459971E0, * .710345683304543E0,.701671914348685E0,.692892887742577E0, * .684009920426076E0,.675024344931163E0,.665937509182049E0, * .656750776292973E0,.647465524363725E0,.638083146272911E0, * .628605049469015E0,.619032655759261E0,.609367401096334E0/ 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)/ * .599610735362968E0,.589764122154454E0,.579829038559083E0, * .569806974936569E0,.559699434694481E0,.549507934062719E0, * .539234001866059E0,.528879179294822E0,.518445019673674E0, * .507933088228616E0,.497344961852181E0,.486682228866890E0, * .475946488786983E0,.465139352078479E0,.454262439917590E0/ 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)/ * .443317383947527E0,.432305826033741E0,.421229418017624E0, * .410089821468717E0,.398888707435459E0,.387627756194516E0, * .376308656998716E0,.364933107823654E0,.353502815112970E0, * .342019493522372E0,.330484865662417E0,.318900661840106E0, * .307268619799319E0,.295590484460136E0,.283868007657082E0/ 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)/ * .272102947876337E0,.260297069991943E0,.248452145001057E0, * .236569949758284E0,.224652266709132E0,.212700883622626E0, * .200717593323127E0,.188704193421389E0,.176662486044902E0, * .164594277567554E0,.152501378338656E0,.140385602411376E0, * .128248767270607E0,.116092693560333E0,.103919204810509E0/ DATA T(248),T(249),T(250),T(251),T(252),T(253),T(254),T(255)/ * .917301271635196E-1,.795272891002330E-1,.673125211657164E-1, * .550876556946340E-1,.428545265363791E-1,.306149687799790E-1, * .183708184788137E-1,.612391237518953E-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.0E0,.347854845137454E0, * .652145154862546E0,.101228536290376E0,.222381034453374E0, * .313706645877887E0,.362683783378362E0,.271524594117541E-1, * .622535239386479E-1,.951585116824928E-1,.124628971255534E0, * .149595988816577E0,.169156519395003E0,.182603415044924E0, * .189450610455068E0/ 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)/.701861000947010E-2, * .162743947309057E-1,.253920653092621E-1,.342738629130214E-1, * .428358980222267E-1,.509980592623762E-1,.586840934785355E-1, * .658222227763618E-1,.723457941088485E-1,.781938957870703E-1, * .833119242269468E-1,.876520930044038E-1,.911738786957639E-1, * .938443990808046E-1,.956387200792749E-1,.965400885147278E-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)/.178328072169643E-2, * .414703326056247E-2,.650445796897836E-2,.884675982636395E-2, * .111681394601311E-1,.134630478967186E-1,.157260304760247E-1, * .179517157756973E-1,.201348231535302E-1,.222701738083833E-1, * .243527025687109E-1,.263774697150547E-1,.283396726142595E-1, * .302346570724025E-1,.320579283548516E-1,.338051618371416E-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)/ .354722132568824E-1, * .370551285402400E-1,.385501531786156E-1,.399537411327203E-1, * .412625632426235E-1,.424735151236536E-1,.435837245293235E-1, * .445905581637566E-1,.454916279274181E-1, .462847965813144E-1, * .469681828162100E-1,.475401657148303E-1,.479993885964583E-1, * .483447622348030E-1,.485754674415034E-1,.486909570091397E-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)/ .449380960292090E-3, * .104581267934035E-2,.164250301866903E-2,.223828843096262E-2, * .283275147145799E-2,.342552604091022E-2,.401625498373864E-2, * .460458425670296E-2,.519016183267633E-2,.577263754286570E-2, * .635166316170719E-2,.692689256689881E-2,.749798192563473E-2, * .806458989048606E-2,.862637779861675E-2,.918300987166087E-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)/.973415341500681E-2, * .102794790158322E-1,.108186607395031E-1,.113513763240804E-1, * .118773073727403E-1,.123961395439509E-1,.129075627392673E-1, * .134112712886163E-1,.139069641329520E-1,.143943450041668E-1, * .148731226021473E-1,.153430107688651E-1,.158037286593993E-1, * .162550009097852E-1,.166965578015892E-1,.171281354231114E-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)/ * .175494758271177E-1,.179603271850087E-1,.183604439373313E-1, * .187495869405447E-1,.191275236099509E-1,.194940280587066E-1, * .198488812328309E-1,.201918710421300E-1,.205227924869601E-1, * .208414477807511E-1,.211476464682213E-1,.214412055392085E-1, * .217219495380521E-1,.219897106684605E-1,.222443288937998E-1, * .224856520327450E-1,.227135358502365E-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)/ * .229278441436868E-1,.231284488243870E-1,.233152299940628E-1, * .234880760165359E-1,.236468835844476E-1,.237915577810034E-1, * .239220121367035E-1,.240381686810241E-1,.241399579890193E-1, * .242273192228152E-1,.243002001679719E-1,.243585572646906E-1, * .244023556338496E-1,.244315690978500E-1,.244461801962625E-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)/ * .112789017822272E-3,.262534944296446E-3,.412463254426176E-3, * .562348954031410E-3,.712154163473321E-3,.861853701420089E-3, * .101142439320844E-2,.116084355756772E-2,.131008868190250E-2, * .145913733331073E-2,.160796713074933E-2,.175655573633073E-2, * .190488085349972E-2,.205292022796614E-2,.220065164983991E-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)/ * .234805295632731E-2,.249510203470371E-2,.264177682542749E-2, * .278805532532771E-2,.293391559082972E-2,.307933574119934E-2, * .322429396179420E-2,.336876850731555E-2,.351273770505631E-2, * .365617995814250E-2,.379907374876626E-2,.394139764140883E-2, * .408313028605267E-2,.422425042138154E-2,.436473687796806E-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)/ * .450456858144790E-2,.464372455568006E-2,.478218392589269E-2, * .491992592181387E-2,.505692988078684E-2,.519317525086928E-2, * .532864159391593E-2,.546330858864431E-2,.559715603368291E-2, * .573016385060144E-2,.586231208692265E-2,.599358091911534E-2, * .612395065556793E-2,.625340173954240E-2,.638191475210788E-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)/ * .650947041505366E-2,.663604959378107E-2,.676163330017380E-2, * .688620269544632E-2,.700973909296982E-2,.713222396107539E-2, * .725363892583391E-2,.737396577381235E-2,.749318645480588E-2, * .761128308454566E-2,.772823794738156E-2,.784403349893971E-2, * .795865236875435E-2,.807207736287350E-2,.818429146643827E-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)/ * .829527784623523E-2,.840501985322154E-2,.851350102502249E-2, * .862070508840101E-2,.872661596169881E-2,.883121775724875E-2, * .893449478375821E-2,.903643154866287E-2,.913701276045081E-2, * .923622333095630E-2,.933404837762327E-2,.943047322573775E-2, * .952548341062928E-2,.961906467984073E-2,.971120299526628E-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)/ * .980188453525733E-2,.989109569669583E-2,.997882309703491E-2, * .100650535763064E-1,.101497741990949E-1,.102329722564782E-1, * .103146352679340E-1,.103947509832117E-1,.104733073841704E-1, * .105502926865815E-1,.106256953418966E-1,.106995040389798E-1, * .107717077058046E-1,.108422955111148E-1,.109112568660490E-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)/ * .109785814257296E-1,.110442590908139E-1,.111082800090098E-1, * .111706345765534E-1,.112313134396497E-1,.112903074958755E-1, * .113476078955455E-1,.114032060430392E-1,.114570935980906E-1, * .115092624770395E-1,.115597048540436E-1,.116084131622531E-1, * .116553800949452E-1,.117005986066207E-1,.117440619140606E-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)/ * .117857634973434E-1,.118256971008240E-1,.118638567340711E-1, * .119002366727665E-1,.119348314595636E-1,.119676359049059E-1, $ .119986450878058E-1,.120278543565826E-1,.120552593295601E-1, * .120808558957245E-1,.121046402153405E-1,.121266087205273E-1, * .121467581157945E-1,.121650853785355E-1,.121815877594818E-1/ DATA W(248),W(249),W(250),W(251),W(252),W(253),W(254),W(255)/ * .121962627831147E-1,.122091082480372E-1,.122201222273040E-1, * .122293030687103E-1,.122366493950402E-1,.122421601042728E-1, * .122458343697479E-1,.122476716402898E-1/ SCALE=(B-A)/2.0 M=N/2 NPLACE=M-1 DO 1 I=1,M S=T(NPLACE+I) R=W(NPLACE+I)*SCALE TT(I)=(A*(1.0+S)+(1.0-S)*B)/2.0 TT( N+1-I)=(A*(1.0-S)+B*(1.0+S))/2.0 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 INTEGER PIVOT REAL 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.0) 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.0) 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.0 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 C PROGRAM PRINTCOF C C TITLE: PRINT GALERKIN COEFFICIENTS FILE. C ---------------------------------------- C C THIS PROGRAM WILL READ TABLES FROM A FORMATFREE FILE PRODUCED BY C SUBROUTINE 'GNRT'. THESE WILL THEN BE PRINTED IN A CONVENIENT FORM C IN DECIMAL. THE USER WILL BE ASKED FOR TWO QUANTITIES. C C INPUT C FNAME: THE NAME FOR THE FORMATFREE FILE OF COEFFICIENTS. C DEG: THE MAXIMUM DEGEE OF THE TABLES FROM FILE 'FNAME' THAT C ARE TO BE PRINTED. THEY WILL BE PRINTED IN DECIMAL, WITH C THE INDICES OF THE ASSOCIATED SPHERICAL HARMONICS. C *** NOTE *** THERE IS AN UPPER LIMIT ON 'DEG', C GIVEN AS 'MAXDEG' IN THE BELOW PARAMETER STATEMENT. C C OUTPUT C FNAMED: THE NAME OF THE OUTPUT FILE OF DECIMAL COEFFICIENTS. C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER DEG,ORDER,DEGREE,FILEIN,FILOUT DOUBLE PRECISION KERMAT CHARACTER FNAME*50,FNAMED*50 PARAMETER(MAXDEG=7,MAXORD=(MAXDEG+1)**2) DIMENSION KERMAT(MAXORD,MAXORD) C C ****** INPUT/OUTPUT **************************************************** C FOLLOWING ARE THE I/O FILE NUMBERS BEING USED. THESE MAY NEED TO C BE CHANGED FOR YOUR COMPUTER. DATA FILEIN/5/,FILOUT/6/ C ************************************************************************ C C OBTAIN INPUT INFORMATION. PRINT *, ' GIVE THE NAME OF YOUR INPUT FILE OF FORMATFREE' PRINT *, ' COEFFICIENTS.' READ(*,'(A)') FNAME PRINT *, ' GIVE A NAME TO BE USED TO DENOTE THE OUTPUT FILE' PRINT *, ' OF DECIMAL COEFFICIENTS.' READ(*,'(A)') FNAMED PRINT *, ' GIVE THE MAXIMUM DEGREE OF THE COEFFICIENTS THAT YOU' PRINT *, ' WANT TO HAVE OUTPUT.' READ *, DEG C L=LEN(FNAME) OPEN(FILEIN,FILE=FNAME(1:L),FORM=UNFORMATTED') L=LEN(FNAMED) OPEN(FILOUT,FILE=FNAMED(1:L)) C C MAIN LOOP. DO 30 K=1,DEG READ(FILEIN)DEGREE,NINNER,NOUTER WRITE(FILOUT,1000) DEGREE,NINNER,NOUTER 1000 FORMAT('1DEGREE=',I2,5X,'NINNER=',I2,5X,'NOUTER=',I2,//, * ' THE INDICES OF A SPHERICAL HARMONIC ARE (N,M,L),'/, * ' THE MEANING OF WHICH IS GIVEN IN FUNCTION LOCATN.',/, * ' THE SUFFIXES R AND C USED BELOW ARE TO INDICATE WHETHER',/, * ' THE ROW OR COLUMN IS BEING REFERENCED.',//) ORDER=(DEGREE+1)**2 READ(FILEIN)((KERMAT(I,J),J=1,ORDER),I=1,ORDER) DO 20 MR=0,DEGREE DO 20 NR=MR,DEGREE DO 20 LR=0,1 IF((LR .EQ. 1) .AND. (MR .EQ. 0)) GO TO 20 LOCR=LOCATN(NR,MR,LR,DEGREE) WRITE(FILOUT,1001) LOCR,NR,MR,LR 1001 FORMAT(/,' ROW ',I2,3X,'(NR,MR,LR)=',3I3,/, * 7X,'COL NC MC LC ELEMENT') DO 10 MC=0,DEGREE DO 10 NC=MC,DEGREE DO 10 LC=0,1 IF((LC .EQ. 1) .AND. (MC .EQ. 0)) GO TO 10 LOCC=LOCATN(NC,MC,LC,DEGREE) WRITE(FILOUT,1002) LOCC,NC,MC,LC,KERMAT(LOCR,LOCC) 1002 FORMAT(7X,I3,3X,I2,3X,I2,3X,I1,1PD20.10) 10 CONTINUE 20 CONTINUE 30 CONTINUE CLOSE(FILEIN) CLOSE(FILOUT) CALL EXIT 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 C PROGRAM PRINTCOF C C TITLE: PRINT GALERKIN COEFFICIENTS FILE. C ---------------------------------------- C C THIS PROGRAM WILL READ TABLES FROM A FORMATFREE FILE PRODUCED BY C SUBROUTINE 'GNRT'. THESE WILL THEN BE PRINTED IN A CONVENIENT FORM C IN DECIMAL. THE USER WILL BE ASKED FOR TWO QUANTITIES. C C INPUT C FNAME: THE NAME FOR THE FORMATFREE FILE OF COEFFICIENTS. C DEG: THE MAXIMUM DEGEE OF THE TABLES FROM FILE 'FNAME' THAT C ARE TO BE PRINTED. THEY WILL BE PRINTED IN DECIMAL, WITH C THE INDICES OF THE ASSOCIATED SPHERICAL HARMONICS. C *** NOTE *** THERE IS AN UPPER LIMIT ON 'DEG', C GIVEN AS 'MAXDEG' IN THE BELOW PARAMETER STATEMENT. C C OUTPUT C FNAMED: THE NAME OF THE OUTPUT FILE OF DECIMAL COEFFICIENTS. C INTEGER DEG,ORDER,DEGREE,FILEIN,FILOUT REAL KERMAT CHARACTER FNAME*50,FNAMED*50 PARAMETER(MAXDEG=7,MAXORD=(MAXDEG+1)**2) DIMENSION KERMAT(MAXORD,MAXORD) C C ****** INPUT/OUTPUT **************************************************** C FOLLOWING ARE THE I/O FILE NUMBERS BEING USED. THESE MAY NEED TO C BE CHANGED FOR YOUR COMPUTER. DATA FILEIN/5/,FILOUT/6/ C ************************************************************************ C C OBTAIN INPUT INFORMATION. PRINT *, ' GIVE THE NAME OF YOUR INPUT FILE OF FORMATFREE' PRINT *, ' COEFFICIENTS.' READ(*,'(A)') FNAME PRINT *, ' GIVE A NAME TO BE USED TO DENOTE THE OUTPUT FILE' PRINT *, ' OF DECIMAL COEFFICIENTS.' READ(*,'(A)') FNAMED PRINT *, ' GIVE THE MAXIMUM DEGREE OF THE COEFFICIENTS THAT YOU' PRINT *, ' WANT TO HAVE OUTPUT.' READ *, DEG C L=LEN(FNAME) OPEN(FILEIN,FILE=FNAME(1:L),FORM='UNFORMATTED') L=LEN(FNAMED) OPEN(FILOUT,FILE=FNAMED(1:L)) C C MAIN LOOP. DO 30 K=1,DEG READ(FILEIN)DEGREE,NINNER,NOUTER WRITE(FILOUT,1000) DEGREE,NINNER,NOUTER 1000 FORMAT('1DEGREE=',I2,5X,'NINNER=',I2,5X,'NOUTER=',I2,//, * ' THE INDICES OF A SPHERICAL HARMONIC ARE (N,M,L),'/, * ' THE MEANING OF WHICH IS GIVEN IN FUNCTION LOCATN.',/, * ' THE SUFFIXES R AND C USED BELOW ARE TO INDICATE WHETHER',/, * ' THE ROW OR COLUMN IS BEING REFERENCED.',//) ORDER=(DEGREE+1)**2 READ(FILEIN)((KERMAT(I,J),J=1,ORDER),I=1,ORDER) DO 20 MR=0,DEGREE DO 20 NR=MR,DEGREE DO 20 LR=0,1 IF((LR .EQ. 1) .AND. (MR .EQ. 0)) GO TO 20 LOCR=LOCATN(NR,MR,LR,DEGREE) WRITE(FILOUT,1001) LOCR,NR,MR,LR 1001 FORMAT(/,' ROW ',I2,3X,'(NR,MR,LR)=',3I3,/, * 7X,'COL NC MC LC ELEMENT') DO 10 MC=0,DEGREE DO 10 NC=MC,DEGREE DO 10 LC=0,1 IF((LC .EQ. 1) .AND. (MC .EQ. 0)) GO TO 10 LOCC=LOCATN(NC,MC,LC,DEGREE) WRITE(FILOUT,1002) LOCC,NC,MC,LC,KERMAT(LOCR,LOCC) 1002 FORMAT(7X,I3,3X,I2,3X,I2,3X,I1,1PE20.10) 10 CONTINUE 20 CONTINUE 30 CONTINUE CLOSE(FILEIN) CLOSE(FILOUT) CALL EXIT 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 DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 831014 (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 RETURNS 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 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 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 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 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 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 HP 9000 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 DENELCOR HEP C (THERE IS NO DOUBLE PRECISION FLOATING POINT) C C C***FIRST EXECUTABLE STATEMENT D1MACH IF (I .LT. 1 .OR. I .GT. 5)RETURN C D1MACH = DMACH(I) RETURN C END REAL FUNCTION R1MACH(I) C***BEGIN PROLOGUE R1MACH C***DATE WRITTEN 790101 (YYMMDD) C***REVISION DATE 831014 (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 RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS C***DESCRIPTION C C R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS C FOR THE LOCAL MACHINE ENVIRONMENT. IT IS A FUNCTION C SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED C AS FOLLOWS, FOR EXAMPLE C C A = R1MACH(I) C C WHERE I=1,...,5. THE (OUTPUT) VALUE OF A ABOVE IS C DETERMINED BY THE (INPUT) VALUE OF I. THE RESULTS FOR C VARIOUS VALUES OF I ARE DISCUSSED BELOW. C C SINGLE-PRECISION MACHINE CONSTANTS C R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. C R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. C R1MACH(5) = LOG10(B) C***REFERENCES FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR C A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE- C MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978, C PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE R1MACH C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) C REAL RMACH(5) C EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA RMACH(1) / Z400800000 / C DATA RMACH(2) / Z5FFFFFFFF / C DATA RMACH(3) / Z4E9800000 / C DATA RMACH(4) / Z4EA800000 / C DATA RMACH(5) / Z500E730E8 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS. C C DATA RMACH(1) / O1771000000000000 / C DATA RMACH(2) / O0777777777777777 / C DATA RMACH(3) / O1311000000000000 / C DATA RMACH(4) / O1301000000000000 / C DATA RMACH(5) / O1157163034761675 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C DATA RMACH(1) / 00564000000000000000B / C DATA RMACH(2) / 37767777777777777776B / C DATA RMACH(3) / 16414000000000000000B / C DATA RMACH(4) / 16424000000000000000B / C DATA RMACH(5) / 17164642023241175720B / C C MACHINE CONSTANTS FOR THE CRAY 1 C C DATA RMACH(1) / 200034000000000000000B / C DATA RMACH(2) / 577767777777777777776B / C DATA RMACH(3) / 377224000000000000000B / C DATA RMACH(4) / 377234000000000000000B / C DATA RMACH(5) / 377774642023241175720B / 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 RMACH(5) C C DATA SMALL/20K,0/,LARGE/77777K,177777K/ C DATA RIGHT/35420K,0/,DIVER/36020K,0/ C DATA LOG10/40423K,42023K/ C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '00000177 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000352 / C DATA DIVER(1),DIVER(2) / '20000000, '00000353 / C DATA LOG10(1),LOG10(2) / '23210115, '00000377 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C C DATA RMACH(1) / O402400000000 / C DATA RMACH(2) / O376777777777 / C DATA RMACH(3) / O714400000000 / C DATA RMACH(4) / O716400000000 / C DATA RMACH(5) / O776464202324 / C C MACHINE CONSTANTS FOR THE HP 2100 C C 3 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 1 / C DATA LARGE91), LARGE(2) / 77777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 325B / C DATA DIVER(1), DIVER(2) / 40000B, 327B / C DATA LOG10(1), LOG10(2) / 46420B, 46777B / 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 RMACH(1) / Z00100000 / C DATA RMACH(2) / Z7FFFFFFF / C DATA RMACH(3) / Z3B100000 / C DATA RMACH(4) / Z3C100000 / C DATA RMACH(5) / Z41134413 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR). C C DATA RMACH(1) / "000400000000 / C DATA RMACH(2) / "377777777777 / C DATA RMACH(3) / "146400000000 / C DATA RMACH(4) / "147400000000 / C DATA RMACH(5) / "177464202324 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1) / 8388608 / C DATA LARGE(1) / 2147483647 / C DATA RIGHT(1) / 880803840 / C DATA DIVER(1) / 889192448 / C DATA LOG10(1) / 1067065499 / C C DATA RMACH(1) / O00040000000 / C DATA RMACH(2) / O17777777777 / C DATA RMACH(3) / O06440000000 / C DATA RMACH(4) / O06500000000 / C DATA RMACH(5) / O07746420233 / 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 LARGE(1),LARGE(2) / 32767, -1 / C DATA RIGHT(1),RIGHT(2) / 13440, 0 / C DATA DIVER(1),DIVER(2) / 13568, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8347 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA RIGHT(1),RIGHT(2) / O032200, O000000 / C DATA DIVER(1),DIVER(2) / O032400, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020233 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA RMACH(1) / O000400000000 / C DATA RMACH(2) / O377777777777 / C DATA RMACH(3) / O146400000000 / C DATA RMACH(4) / O147400000000 / C DATA RMACH(5) / O177464202324 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS*** C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1) / 128 / C DATA LARGE(1) / -32769 / C DATA RIGHT(1) / 13440 / C DATA DIVER(1) / 13568 / C DATA LOG10(1) / 547045274 / C C DATA SMALL(1) / Z00000080 / C DATA LARGE(1) / ZFFFF7FFF / C DATA RIGHT(1) / Z00003480 / C DATA DIVER(1) / Z00003500 / C DATA LOG10(1) / Z209B3F9A / C C MACHINE CONSTANTS FOR THE HP 9000 C C DATA SMALL(1) / 00040000000B/ C DATA LARGE(1) / 17677777777B/ C DATA RIGHT(1) / 06340000000B/ C DATA DIVER(1) / 06400000000B/ C DATA LOG10(1) / 07646420233B/ C C MACHINE CONSTANTS FOR THE DENELCOR HEP C C DATA SMALL(1) / Z0010000000000000/ C DATA LARGE(1) / Z7FFFFFFFFFFFFFFF/ C DATA RIGHT(1) / Z3310000000000000/ C DATA DIVER(1) / Z3410000000000000/ C DATA LOG10(1) / Z41134413509F79FF/ C C MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR C C DATA SMALL(1),SMALL(2) / 0, 256/ C DATA LARGE(1),LARGE(2) / -1, -129/ C DATA RIGHT(1),RIGHT(2) / 0, 26880/ C DATA DIVER(1),DIVER(2) / 0, 27136/ C DATA LOG10(1),LOG10(2) / 8347, 32538/ C C C***FIRST EXECUTABLE STATEMENT R1MACH IF (I .LT. 1 .OR. I .GT. 5)RETURN C R1MACH = RMACH(I) RETURN C END