      FUNCTION F1DIM(X)
      PARAMETER (NMAX=50)
      COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
      DIMENSION XT(NMAX)
      DO 11 J=1,NCOM
        XT(J)=PCOM(J)+X*XICOM(J)
11    CONTINUE
      F1DIM=FUNC(XT)
      RETURN
      END
      FUNCTION FACTLN(N)
      DIMENSION A(100)
      DATA A/100*-1./
      IF (N.LT.0) PAUSE 'negative factorial'
      IF (N.LE.99) THEN
        IF (A(N+1).LT.0.) A(N+1)=GAMMLN(N+1.)
        FACTLN=A(N+1)
      ELSE
        FACTLN=GAMMLN(N+1.)
      ENDIF
      RETURN
      END
      FUNCTION FACTRL(N)
      DIMENSION A(33)
      DATA NTOP,A(1)/0,1./
      IF (N.LT.0) THEN
        PAUSE 'negative factorial'
      ELSE IF (N.LE.NTOP) THEN
        FACTRL=A(N+1)
      ELSE IF (N.LE.32) THEN
        DO 11 J=NTOP+1,N
          A(J+1)=J*A(J)
11      CONTINUE
        NTOP=N
        FACTRL=A(N+1)
      ELSE
        FACTRL=EXP(GAMMLN(N+1.))
      ENDIF
      RETURN
      END
      SUBROUTINE FGAUSS(X,A,Y,DYDA,NA)
      DIMENSION A(NA),DYDA(NA)
      Y=0.
      DO 11 I=1,NA-1,3
        ARG=(X-A(I+1))/A(I+2)
        EX=EXP(-ARG**2)
        FAC=A(I)*EX*2.*ARG
        Y=Y+A(I)*EX
        DYDA(I)=EX
        DYDA(I+1)=FAC/A(I+2)
        DYDA(I+2)=FAC*ARG/A(I+2)
11    CONTINUE
      RETURN
      END
      SUBROUTINE FIT(X,Y,NDATA,SIG,MWT,A,B,SIGA,SIGB,CHI2,Q)
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA)
      SX=0.
      SY=0.
      ST2=0.
      B=0.
      IF(MWT.NE.0) THEN
        SS=0.
        DO 11 I=1,NDATA
          WT=1./(SIG(I)**2)
          SS=SS+WT
          SX=SX+X(I)*WT
          SY=SY+Y(I)*WT
11      CONTINUE
      ELSE
        DO 12 I=1,NDATA
          SX=SX+X(I)
          SY=SY+Y(I)
12      CONTINUE
        SS=FLOAT(NDATA)
      ENDIF
      SXOSS=SX/SS
      IF(MWT.NE.0) THEN
        DO 13 I=1,NDATA
          T=(X(I)-SXOSS)/SIG(I)
          ST2=ST2+T*T
          B=B+T*Y(I)/SIG(I)
13      CONTINUE
      ELSE
        DO 14 I=1,NDATA
          T=X(I)-SXOSS
          ST2=ST2+T*T
          B=B+T*Y(I)
14      CONTINUE
      ENDIF
      B=B/ST2
      A=(SY-SX*B)/SS
      SIGA=SQRT((1.+SX*SX/(SS*ST2))/SS)
      SIGB=SQRT(1./ST2)
      CHI2=0.
      IF(MWT.EQ.0) THEN
        DO 15 I=1,NDATA
          CHI2=CHI2+(Y(I)-A-B*X(I))**2
15      CONTINUE
        Q=1.
        SIGDAT=SQRT(CHI2/(NDATA-2))
        SIGA=SIGA*SIGDAT
        SIGB=SIGB*SIGDAT
      ELSE
        DO 16 I=1,NDATA
          CHI2=CHI2+((Y(I)-A-B*X(I))/SIG(I))**2
16      CONTINUE
        Q=GAMMQ(0.5*(NDATA-2),0.5*CHI2)
      ENDIF
      RETURN
      END
      SUBROUTINE FIXRTS(D,NPOLES)
      PARAMETER (NPMAX=100)
      DIMENSION D(NPOLES)
      LOGICAL POLISH
      COMPLEX A(NPMAX),ROOTS(NPMAX)
      A(NPOLES+1)=CMPLX(1.,0.)
      DO 11 J=NPOLES,1,-1
        A(J)=CMPLX(-D(NPOLES+1-J),0.)
11    CONTINUE
      POLISH=.TRUE.
      CALL ZROOTS(A,NPOLES,ROOTS,POLISH)
      DO 12 J=1,NPOLES
        IF(CABS(ROOTS(J)).GT.1.)THEN
          ROOTS(J)=1./CONJG(ROOTS(J))
        ENDIF
12    CONTINUE
      A(1)=-ROOTS(1)
      A(2)=CMPLX(1.,0.)
      DO 14 J=2,NPOLES
        A(J+1)=CMPLX(1.,0.)
        DO 13 I=J,2,-1
          A(I)=A(I-1)-ROOTS(J)*A(I)
13      CONTINUE
        A(1)=-ROOTS(J)*A(1)
14    CONTINUE
      DO 15 J=1,NPOLES
        D(NPOLES+1-J)=-REAL(A(J))
15    CONTINUE
      RETURN
      END
      SUBROUTINE FLEG(X,PL,NL)
      DIMENSION PL(NL)
      PL(1)=1.
      PL(2)=X
      IF(NL.GT.2) THEN
        TWOX=2.*X
        F2=X
        D=1.
        DO 11 J=3,NL
          F1=D
          F2=F2+TWOX
          D=D+1.
          PL(J)=(F2*PL(J-1)-F1*PL(J-2))/D
11      CONTINUE
      ENDIF
      RETURN
      END
      SUBROUTINE FLMOON(N,NPH,JD,FRAC)
      PARAMETER (RAD=0.017453293)
      C=N+NPH/4.
      T=C/1236.85
      T2=T**2
      AS=359.2242+29.105356*C
      AM=306.0253+385.816918*C+0.010730*T2
      JD=2415020+28*N+7*NPH
      XTRA=0.75933+1.53058868*C+(1.178E-4-1.55E-7*T)*T2
      IF(NPH.EQ.0.OR.NPH.EQ.2)THEN
        XTRA=XTRA+(0.1734-3.93E-4*T)*SIN(RAD*AS)-0.4068*SIN(RAD*AM)
      ELSE IF(NPH.EQ.1.OR.NPH.EQ.3)THEN
        XTRA=XTRA+(0.1721-4.E-4*T)*SIN(RAD*AS)-0.6280*SIN(RAD*AM)
      ELSE
        PAUSE 'NPH is unknown.'
      ENDIF
      IF(XTRA.GE.0.)THEN
        I=INT(XTRA)
      ELSE
        I=INT(XTRA-1.)
      ENDIF
      JD=JD+I
      FRAC=XTRA-I
      RETURN
      END
      SUBROUTINE FOUR1(DATA,NN,ISIGN)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION DATA(*)
      N=2*NN
      J=1
      DO 11 I=1,N,2
        IF(J.GT.I)THEN
          TEMPR=DATA(J)
          TEMPI=DATA(J+1)
          DATA(J)=DATA(I)
          DATA(J+1)=DATA(I+1)
          DATA(I)=TEMPR
          DATA(I+1)=TEMPI
        ENDIF
        M=N/2
1       IF ((M.GE.2).AND.(J.GT.M)) THEN
          J=J-M
          M=M/2
        GO TO 1
        ENDIF
        J=J+M
11    CONTINUE
      MMAX=2
2     IF (N.GT.MMAX) THEN
        ISTEP=2*MMAX
        THETA=6.28318530717959D0/(ISIGN*MMAX)
        WPR=-2.D0*DSIN(0.5D0*THETA)**2
        WPI=DSIN(THETA)
        WR=1.D0
        WI=0.D0
        DO 13 M=1,MMAX,2
          DO 12 I=M,N,ISTEP
            J=I+MMAX
            TEMPR=SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+1)
            TEMPI=SNGL(WR)*DATA(J+1)+SNGL(WI)*DATA(J)
            DATA(J)=DATA(I)-TEMPR
            DATA(J+1)=DATA(I+1)-TEMPI
            DATA(I)=DATA(I)+TEMPR
            DATA(I+1)=DATA(I+1)+TEMPI
12        CONTINUE
          WTEMP=WR
          WR=WR*WPR-WI*WPI+WR
          WI=WI*WPR+WTEMP*WPI+WI
13      CONTINUE
        MMAX=ISTEP
      GO TO 2
      ENDIF
      RETURN
      END
      SUBROUTINE FOURN(DATA,NN,NDIM,ISIGN)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION NN(NDIM),DATA(*)
      NTOT=1
      DO 11 IDIM=1,NDIM
        NTOT=NTOT*NN(IDIM)
11    CONTINUE
      NPREV=1
      DO 18 IDIM=1,NDIM
        N=NN(IDIM)
        NREM=NTOT/(N*NPREV)
        IP1=2*NPREV
        IP2=IP1*N
        IP3=IP2*NREM
        I2REV=1
        DO 14 I2=1,IP2,IP1
          IF(I2.LT.I2REV)THEN
            DO 13 I1=I2,I2+IP1-2,2
              DO 12 I3=I1,IP3,IP2
                I3REV=I2REV+I3-I2
                TEMPR=DATA(I3)
                TEMPI=DATA(I3+1)
                DATA(I3)=DATA(I3REV)
                DATA(I3+1)=DATA(I3REV+1)
                DATA(I3REV)=TEMPR
                DATA(I3REV+1)=TEMPI
12            CONTINUE
13          CONTINUE
          ENDIF
          IBIT=IP2/2
1         IF ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) THEN
            I2REV=I2REV-IBIT
            IBIT=IBIT/2
          GO TO 1
          ENDIF
          I2REV=I2REV+IBIT
14      CONTINUE
        IFP1=IP1
2       IF(IFP1.LT.IP2)THEN
          IFP2=2*IFP1
          THETA=ISIGN*6.28318530717959D0/(IFP2/IP1)
          WPR=-2.D0*DSIN(0.5D0*THETA)**2
          WPI=DSIN(THETA)
          WR=1.D0
          WI=0.D0
          DO 17 I3=1,IFP1,IP1
            DO 16 I1=I3,I3+IP1-2,2
              DO 15 I2=I1,IP3,IFP2
                K1=I2
                K2=K1+IFP1
                TEMPR=SNGL(WR)*DATA(K2)-SNGL(WI)*DATA(K2+1)
                TEMPI=SNGL(WR)*DATA(K2+1)+SNGL(WI)*DATA(K2)
                DATA(K2)=DATA(K1)-TEMPR
                DATA(K2+1)=DATA(K1+1)-TEMPI
                DATA(K1)=DATA(K1)+TEMPR
                DATA(K1+1)=DATA(K1+1)+TEMPI
15            CONTINUE
16          CONTINUE
            WTEMP=WR
            WR=WR*WPR-WI*WPI+WR
            WI=WI*WPR+WTEMP*WPI+WI
17        CONTINUE
          IFP1=IFP2
        GO TO 2
        ENDIF
        NPREV=N*NPREV
18    CONTINUE
      RETURN
      END
	PROGRAM D12R13
C	Driver for routine FOURN
	PARAMETER(NDIM=3,NDAT=1024)
	DIMENSION NN(NDIM),DATA(NDAT)
	DO 11 I=1,NDIM
		NN(I)=2*(2**I)
11	CONTINUE
	DO 14 I=1,NN(3)
		DO 13 J=1,NN(2)
			DO 12 K=1,NN(1)
				L=K+(J-1)*NN(1)+(I-1)*NN(2)*NN(1)
 				LL=2*L-1
				DATA(LL)=FLOAT(LL)
				DATA(LL+1)=FLOAT(LL+1)
12			CONTINUE
13		CONTINUE
14	CONTINUE
	ISIGN=+1
	CALL FOURN(DATA,NN,NDIM,ISIGN)
	ISIGN=-1
	WRITE(*,'(1X,A)') 'Double 3-dimensional Transform'
	WRITE(*,'(/1X,T10,A,T35,A,T63,A)') 'Double Transf.',
     *		'Original Data','Ratio'
	WRITE(*,'(1X,T8,A,T20,A,T33,A,T45,A,T57,A,T69,A/)') 
     *		'Real','Imag.','Real','Imag.','Real','Imag.'
	CALL FOURN(DATA,NN,NDIM,ISIGN)
	DO 15 I=1,4
		J=2*I
		K=2*J
		L=K+(J-1)*NN(1)+(I-1)*NN(2)*NN(1)
		LL=2*L-1
		WRITE(*,'(1X,6F12.2)') DATA(LL),DATA(LL+1),FLOAT(LL),FLOAT(LL+1),
     *			DATA(LL)/LL,DATA(LL+1)/(LL+1)
15	CONTINUE
	WRITE(*,'(/1X,A,I4)') 'The product of transform lengths is:',
     *		NN(1)*NN(2)*NN(3)
	END
      SUBROUTINE FPOLY(X,P,NP)
      DIMENSION P(NP)
      P(1)=1.
      DO 11 J=2,NP
        P(J)=P(J-1)*X
11    CONTINUE
      RETURN
      END
      SUBROUTINE FRPRMN(P,N,FTOL,ITER,FRET)
      PARAMETER (NMAX=50,ITMAX=200,EPS=1.E-10)
      DIMENSION P(N),G(NMAX),H(NMAX),XI(NMAX)
      FP=FUNC(P)
      CALL DFUNC(P,XI)
      DO 11 J=1,N
        G(J)=-XI(J)
        H(J)=G(J)
        XI(J)=H(J)
11    CONTINUE
      DO 14 ITS=1,ITMAX
        ITER=ITS
        CALL LINMIN(P,XI,N,FRET)
        IF(2.*ABS(FRET-FP).LE.FTOL*(ABS(FRET)+ABS(FP)+EPS))RETURN
        FP=FUNC(P)
        CALL DFUNC(P,XI)
        GG=0.
        DGG=0.
        DO 12 J=1,N
          GG=GG+G(J)**2
C         DGG=DGG+XI(J)**2
          DGG=DGG+(XI(J)+G(J))*XI(J)
12      CONTINUE
        IF(GG.EQ.0.)RETURN
        GAM=DGG/GG
        DO 13 J=1,N
          G(J)=-XI(J)
          H(J)=G(J)+GAM*H(J)
          XI(J)=H(J)
13      CONTINUE
14    CONTINUE
      PAUSE 'FRPR maximum iterations exceeded'
      RETURN
      END
      SUBROUTINE FTEST(DATA1,N1,DATA2,N2,F,PROB)
      DIMENSION DATA1(N1),DATA2(N2)
      CALL AVEVAR(DATA1,N1,AVE1,VAR1)
      CALL AVEVAR(DATA2,N2,AVE2,VAR2)
      IF(VAR1.GT.VAR2)THEN
        F=VAR1/VAR2
        DF1=N1-1
        DF2=N2-1
      ELSE
        F=VAR2/VAR1
        DF1=N2-1
        DF2=N1-1
      ENDIF
      PROB = BETAI(0.5*DF2,0.5*DF1,DF2/(DF2+DF1*F))
     *    +(1.-BETAI(0.5*DF1,0.5*DF2,DF1/(DF1+DF2/F)))
      RETURN
      END
