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