SUBROUTINE GOODWN(INX,INY,N,F) C C THIS ROUTINE COMPUTES INTEGRALS J(X,Y;N) AS DEFINED BY C CROUCH & SPIEGELMAN, JASA, 1990, USING GOODWN'S METHOD C C INPUT PARAMETERS: INX PARAMETER X OF J(X,Y;N) C INY PARAMETER Y OF J(X,Y;N) C N PARAMETER N OF J(X,Y;N) C OUTPUT F F=J(X,Y;N) C IMPLICIT REAL*8 (A-Z) INTEGER N,NCOM PARAMETER(EPS=1.D-16) EXTERNAL LGFT COMMON /COMFOF/X,Y,C,NCOM ZERO=0.D0 ONE=1.D0 TWO=2.D0 PI=3.14159265358979324D0 RTPI=dSQRT(PI) C C DEFINE X TO GET PARAMETERIZATION FOR J(X,Y;N) C Y=dABS(INY) C=DBLE(N)*INX+DBLE(N)**2*Y**2/4.D0 X=INX+Y**2/TWO*DBLE(N) XN=N+1 NCOM=N+1 C C SINCE I(X,Y;N)=I(X,-Y;N), SET Y TO dABS(Y) C C C FIRST, FIND H C K=PI/TWO/dSQRT(dLOG(TWO/EPS)) IF(Y.GT.K)THEN H=4.D0*Y*K**2/(Y**2+K**2) ELSE H=TWO*K ENDIF T0=MAXPNT(X,Y,XN) C C NOW, WE'RE READY TO SUM UP THE INTEGRAL C 15 CONTINUE B=X+Y*T0 A=-T0**2+C IF(B*XN.GT.30.D0)THEN F=dEXP(A-XN*B)/(ONE+dEXP(-B))**(N+1) ELSE F=dEXP(-T0**2+C)/(ONE+dEXP(X+Y*T0))**(N+1) ENDIF IF(F.EQ.ZERO)THEN RETURN ENDIF C C SUM FROM MIDDLE UP C K=ZERO 20 CONTINUE K=K+ONE C C THIS FORM OF THE ARITHMETIC OPERATION REDUCES CHANCES OF OVERFLOW C FOR CASES WHERE K GETS LARGE C A=C-(T0+K*H)**2 B=X+Y*(T0+K*H) iF(b*xn.gt.30.d0)then TERM=dEXP(A-(N+1)*B)/(ONE+dEXP(-B))**(N+1) else term=dexp(a)/(one+dexp(b))**(n+1) endif F=F+TERM IF(dABS(TERM/F).GT.EPS)GO TO 20 K=ZERO 30 CONTINUE C C SUM FROM MIDDLE DOWN C K=K-ONE A=C-(T0+K*H)**2 B=X+Y*(T0+K*H) IF(B*XN.GT.30.D0)THEN TERM=dEXP(A-XN*B)/(ONE+dEXP(-B))**(N+1) ELSE TERM=dEXP(A)/(ONE+dEXP(B))**(N+1) ENDIF F=F+TERM IF(dABS(TERM/F).GT.EPS)GO TO 30 F=H*F/RTPI RETURN END C C C C C C REAL FUNCTION MAXPNT*8(X,Y,N) IMPLICIT REAL*8 (A-Z) INTEGER COMN EXTERNAL LGFT COMMON /COMFOF/COMX,COMY,COMC,COMN ONE=1.D0 TWO=2.D0 XGUESS=-DBLE(COMN)*COMY/TWO c c erset is an imsl routine - delete if you find the minimum of function LGFT c some other way c CALL ERSET(5,-1,0) CALL ERSET(4,-1,0) IF(COMY.LT.0.D0)THEN LB=0.D0 UB=10.D0*XGUESS ELSE LB=10.D0*XGUESS UB=0.D0 ENDIF c c duvmif is an imsl routine which finds the minimum of a function. If imsl is c not available to you, you must find another equivalent function, or write c one yourself. Numerical recipes has some good choices at very low cost. c CALL DUVMIF(LGFT,XGUESS,1.D0,1.d6,1.D-6,100,T0) CALL ERSET(5,1,2) CALL ERSET(4,1,2) MAXPNT=T0 RETURN END C C C REAL FUNCTION LGFT*8(T) COMMON /COMFOF/X,Y,C,N REAL*8 T,X,Y,C,a,b INTEGER N ONE=1.D0 A=C-T**2 B=X+Y*T IF(B.GT.30.D0)THEN LGFT=A-DBLE(N)*B ELSE LGFT=A-DBLE(N)*dLOG(ONE+dEXP(B)) ENDIF LGFT=-LGFT RETURN END