program gaussfit
      implicit none
      integer nka
      parameter (nka=40)
      integer nkmax
      double precision kx(nka,nka,nka),ky(nka,nka,nka),kz(nka,nka,nka),
     +     c(nka,nka,nka),d(nka,nka,nka)
      common /arrayinfo/nkmax,kx,ky,kz,c,d
      integer ix,iy,iz
      double precision rx,ry,rz,lambda,chisquare
      double precision XOUT(4),ACC,WK(33),FOUT,F
      integer NEV,N
      external SIMPLX,f
      N=4
      rx=5.d0
      ry=5.d0
      rz=5.d0
      lambda=1.d0
      write(6,*) 'what are guesses rx,ry,rz and lambda?'
      read(5,*) rx,ry,rz,lambda
      XOUT(1)=rx
      XOUT(2)=ry
      XOUT(3)=rz
      XOUT(4)=lambda
      NEV=1000
      WK(1)=0.2
      write(6,*) 'What is the accuracy parameter, ACC?'
      read(5,*) ACC
      open(67,form='formatted',status='unknown',name='results3d.txt')
      read(67,*) nkmax
      do ix=1,nkmax
         do iy=1,nkmax
            do iz=1,nkmax
               read(67,*) kx(ix,iy,iz),ky(ix,iy,iz),kz(ix,iy,iz),
     +              c(ix,iy,iz),d(ix,iy,iz)
            enddo
         enddo
      enddo
      chisquare=f(XOUT)
      write(6,*) 'initial guess yields chisquare=',chisquare
      call SIMPLX(F,N,XOUT,FOUT,NEV,ACC,WK)
      rx=XOUT(1)
      ry=XOUT(2)
      rz=XOUT(3)
      lambda=XOUT(4)
      chisquare=FOUT
      write(6,*) 'final chisquare=',chisquare
      write(6,*) rx,ry,rz,lambda
      stop
      end

      double precision function f(XOUT)
      implicit none
      integer nka
      parameter(nka=40)
      integer nkmax
      double precision kx(nka,nka,nka),ky(nka,nka,nka),kz(nka,nka,nka),
     +     c(nka,nka,nka),d(nka,nka,nka)
      common /arrayinfo/nkmax,kx,ky,kz,c,d
      integer ix,iy,iz
      double precision XOUT(4)
      double precision arg,answer
      double precision rx,ry,rz,lambda
      answer=0.d0
      rx=XOUT(1)
      ry=XOUT(2)
      rz=XOUT(3)
      lambda=XOUT(4)
      do ix=1,nkmax
         do iy=1,nkmax
            do iz=1,nkmax
               arg=(kx(ix,iy,iz)*rx)**2 +(ky(ix,iy,iz)*ry)**2
     +              +(kz(ix,iy,iz)*rz)**2
               arg=arg/(197.323*197.323)
               answer=answer+d(ix,iy,iz)*
     +              (c(ix,iy,iz)-1.0-lambda*exp(-arg))**2
            enddo
         enddo
      enddo
      answer=(6.0/3.141592654d0)*answer/(nkmax*nkmax*nkmax)
c     write(6,*) 'inside f, ',rx,ry,rz,lambda,answer
      f=answer
      return
      end

      SUBROUTINE SIMPLX(F,N,XOUT,FOUT,NEV,ACC,WK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     
C     ROUTINE SEARCHES FOR A MINIMUM IN N-DIM SPACE OF PARAMETERS
C     
C     --INPUT--
C     F   - FUNCTION, IN A CALLING PROGRAM F(X), WHERE X IS AN ARRAY OF N
C     PARAMETERS, MUST APPEAR IN AN EXTERNAL
C     N   - NUMBER OF PARAMETERS
C     XOUT- GUESS FOR A MINIMUM
C     NEV - LIMIT ON FUNCTION EVALUATIONS
C     ACC - REQUIRED ACCURACY, .GT.0. ABSOLUTE, .LT.0. RELATIVE
C     WK  - WORK ARRAY, MUST BE DIMENSIONED WK(N*N+4*N+1), AT AN ENTRY
C     WK(1) IS A SCALE PARAMETER, IF WK(1).LT.0. THEN WK(I)
C     I=1,...,N  CONTAINS SCALE PARAMETERS, ATTENTION - CANNOT BE TO
C     SMALL
C     
C     --OUTPUT--
C     XOUT- LOCATION OF MINIMUM
C     FOUT- VALUE OF MINIMUM
C     NEV - ACTUAL NUMBER OF EVALUTIONS
C     ACC - ESTIMATED ACCURACY
C     
      EXTERNAL F
      DIMENSION XOUT(*),WK(*)
C     
      CALL SIMPLE(F,N,N+1,XOUT,FOUT,NEV,ACC,WK(N+2),WK,WK(N*N+2*N+2),
     S     WK(N*N+3*N+2))
C     
      RETURN
      END


      SUBROUTINE SIMPLE(F,N,M,XOUT,FOUT,NEV,ACC,X,FX,XR,XRR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      EXTERNAL F
      LOGICAL L
      DIMENSION X(N,M)
      DIMENSION XOUT(N),FX(M),XR(N),XRR(N)
C
      L=ACC.LT.0.D0
      PN=1.D0/N
      N1=N+1
      PN1=1.D0/N1
C
C  INITIAL POINTS
      IF(FX(1).EQ.0.D0)FX(1)=1.D0
C
      IF(FX(1).GT.0.D0)GOTO 1
C
      FX(1)=-FX(1)
      GOTO 3
C
 1    DO 2 I=1,N
 2    FX(I)=FX(1)
C
 3    DO 4 I=1,N
      FR=XOUT(I) !-PN1*FX(I)
      DO 4 J=1,N1
 4    X(I,J)=FR
C
      DO 5 J=1,N
 5    X(J,J+1)=X(J,J+1)+FX(J)
C
C
C  VALUES AT INITIAL POINTS
      DO 10 J=1,N1
 10   FX(J)=F(X(1,J))
      IEV=N1+1
C
C  THE SEARCH SETS IN
 20   JL=1
      JH=N1
      DO 30 J=1,N1
      IF(FX(J).GT.FX(JH))JH=J
      IF(FX(J).LT.FX(JL))JL=J
 30   CONTINUE
C
C  TESTING FOR EXIT
      FRR=0.D0
      DO 32 J=1,N1
 32   FRR=FRR+(FX(J)-FX(JL))**2
      FRR=SQRT(PN*FRR)
      FR=ACC
      IF(L)FR=ABS(FX(JL)*ACC)
      IF(FRR.LE.FR.OR.IEV.GE.NEV)GOTO 190
C
C  EVALUATING A MEAN
 38   DO 50 I=1,N
      FR=0.D0
      DO 40 J=1,N1
 40   FR=FR+X(I,J)
 50   XOUT(I)=(FR-X(I,JH))*PN
C
C  REFLECTED POINT
      DO 60 I=1,N
 60   XR(I)=2.D0*XOUT(I)-X(I,JH)
      FR=F(XR)
      IEV=IEV+1
C
C  TESTING FOR EXPANSION
      IF(FR.LT.FX(JL))GOTO 130
C
C  TESTING FOR REFLECTED POINT
      DO 70 J=1,N1
      IF(FR.LE.FX(J).AND.J.NE.JH)GOTO 170
 70   CONTINUE
C
C  TESTING FOR A DIRECT REDUCED ATTEMPT
      IF(FR.GT.FX(JH))GOTO 90
C
C  INSERTION OF THE REFLECTED POINT
      FX(JH)=FR
      DO 80 I=1,N
 80   X(I,JH)=XR(I)
C
C  REDUCED ATTEMPT
 90   DO 100 I=1,N
 100  XRR(I)=.5D0*(X(I,JH)+XOUT(I))
      FRR=F(XRR)
      IEV=IEV+1
C
C  TESTING FOR THE EXPANDED POINT
      IF(FRR.LE.FX(JH))GOTO 150
C
C  REDUCING SIMPLEX
      DO 120 J=1,N1
      IF(J.EQ.JL)GOTO 120
      DO 110 I=1,N
 110  X(I,J)=.5D0*(X(I,J)+X(I,JL))
      FX(J)=F(X(1,J))
 120  CONTINUE
      IEV=IEV+N
      GOTO 20
C
C  EXPANDING IN THE DIRECTION
 130  DO 140 I=1,N
 140  XRR(I)=4.D0*XOUT(I)-3.D0*X(I,JH)
      FRR=F(XRR)
      IEV=IEV+1
C
C  TESTING FOR THE REFLECTED POINT
      IF(FRR.GE.FX(JL))GOTO 170
C
C  ACCEPTING EXPANDED POINT
 150  FX(JH)=FRR
      DO 160 I=1,N
 160  X(I,JH)=XRR(I)
      GOTO 20
C
C  ACCEPTING REFLECTED POINT
 170  FX(JH)=FR
      DO 180 I=1,N
 180  X(I,JH)=XR(I)
      GOTO 20
C
C  EXIT
 190  NEV=IEV
C
      DO 210 I=1,N
      FR=0.D0
      DO 200 J=1,N1
 200  FR=FR+X(I,J)
 210  XOUT(I)=FR*PN1
      FOUT=F(XOUT)
C
      IF(FOUT.LE.FX(JL))GOTO 230
C
      FOUT=FX(JL)
      DO 220 I=1,N
 220  XOUT(I)=X(I,JL)
C
 230  ACC=FRR
      IF(L.AND.FOUT.NE.0.)ACC=ACC/ABS(FOUT)
C
      RETURN
      END