Module data

integer, public :: maxpart,npoints,min_part,max_part,nfiles, nsuffix
real, public, allocatable :: entropy(:,:)   ! contains logs of counts for a certain energy+number
                                    ! entropy (0,maxpart) is the starting (lowest) value of the energy
                                    ! entropy (-1,maxpart) is the highest active value of the energy
real, public :: width                ! width of energy distributions in each histogram file
real smix                    ! s field mixing parameter
real, public, allocatable :: beta(:,:),mu(:,:) ! temperature, chemical potential for each histogram file
real beta1,mu1            ! temperature and chem. pot. of interest for current calcs.
real, public :: avg_N    ! average N
real, allocatable :: p_l(:)            ! ordering operator distribution
integer ncrit                ! number of points for p_l

end module

program entropy3

! debug: gfortran -ffree-line-length-none -g -fdefault-real-8 -ffpe-trap=underflow,invalid,overflow  -o entropy3_dbg.out entropy3.f90
! compile: gfortran -ffree-line-length-none -O2 -fdefault-real-8 -o entropy3.out entropy3.f90
! version 3.0 - Apr. '24: now uses the Tsypin and Bloete function - 3D only
! version 2.4 - Dec. '23: automatic error bars, going over all suffix files separately
! version 2.39 - Jul '18: ported to gfortran
! version 2.38 - Oct '05: fixed problems for no observations at intermediate ipart
! version 2.37 - June '04: directory read from entr2_par.dat file
! version 2.36 - May '04: olddev, rewind(26), MATLAB compatibility
! version 2.35 - Jul '03: dmu_dT variable, new obj. function /w more weight at peaks
! version 2.34 - Jun '03: fixed bug /w 2D capability
! version 2.33 - Aug '02: 2D capability
! version 2.21 - Oct '99 (AZP)
! Diff from 2.2 - output formats
! skipped 2.1 to conform to entropy2 v. 2.2
! version 2.0 - February '97 (AZP)
! based on version 1.1 - July '96 (AZP)
! program takes as input observations of histograms of N and E
! and weights calculated by the
! Ferrenberg and Swedsen algorithm (PRL, 1989) (from program entropy2)
! and optimizes the Critical Temperature, Critical Chem. Pot., 
! field mixing parameter S and the number of bins.

use data

implicit none

integer, parameter :: nopt_dim = 2  ! number of optimizing directions: Mu_c, smix
real ptry(nopt_dim)
real, parameter :: min_ene = -999999.	! minimum energy
real, allocatable :: omega(:,:)   ! temporary holding variable, reads contr. to entropy from
								! each file being read
real energy_start               ! starting value of energy for each ipart in histogram file
real, allocatable :: counts(:)		! counts for every file being read per (E,N) point
real, allocatable :: counts_file(:,:)	! total counts for each file
integer ncount	! counts how many observations there are for a certain number of particles
integer ipart,i,idiff	! local counters
integer min_part1,max_part1		! min. and max. number of particles in file
real x						! ordering operator distribution x axis
real max_x, min_x, avg_x, dx	! limits of x, average x
integer index				! index for x
real variance				! variance of p_l
real, allocatable :: weight(:,:)		! f_n in Ferr. + Swed. paper
real num_P				! Numerator for P(S,K) in eq. (3) of Ferr. + Swed. paper
real probab				! probability of microstate at (E,N) for given T,mu
integer ifile,jfile,jsuffix	! how many files being read
integer iostat				! iostat becomes -1 on end-of-file
real Z					! partition function
real, allocatable :: dens(:) ! dens stores the probability of a certain number
integer iter
integer isuffix	! number of different run series
real width1, xdim1, ydim1, zdim1	! these are only read, not used anywhere

real p(nopt_dim+1,nopt_dim)		! input to optimization function
real y(nopt_dim+1)				! input to optimization function
real ftol ! tolerance for convergence of optimization
real trial_vec(nopt_dim)		! factors for initial search in Tcrit
real funk
external funk
real zfactor	! auxiliary variable in calculation of weight and Z
real Tstart,Tend,Tincr,dmu_dT	! Temperatures to be studied, change of mu vs. T at rho = rho_c
integer Ntemp,Itemp		! Number of temperatures, increment in temperature
real ene_liq,avgnum	! mean energy, average number of particles
real olddev
real mu0,smix0   ! initial values of parameters
real Tall,dall,Muall,Sall,Nall
integer ipass, suff_start, suff_end
real, allocatable :: Topt(:),dopt(:),Muopt(:),Sopt(:),Nopt(:) ! [nsuffix]
real Tstd,dstd,Mustd,Sstd,Nstd
character(20), allocatable :: filename(:)    ! histogram file name - as hisXX.dat; only the XX part is stored
character(2), allocatable :: suffix(:)    ! for dealing with series of runs

open (1,file='entr3_par.dat')
read (1,*) (trial_vec(i),i=1,nopt_dim)
read (1,*) ftol, dmu_dT
close (1)

open (1,file='input_hs2.dat')
read (1,*) nsuffix
allocate (suffix(nsuffix))
allocate (Topt(nsuffix),dopt(nsuffix),Muopt(nsuffix),Sopt(nsuffix),Nopt(nsuffix))
do isuffix = 1,nsuffix
	read (1,*) suffix(isuffix)
enddo
read (1,*) nfiles, npoints, maxpart
close(1)

allocate (entropy(-1:npoints,0:maxpart),omega(-1:npoints,0:maxpart))
allocate (counts(npoints),counts_file(nfiles,nsuffix))
allocate (beta(nfiles,nsuffix),mu(nfiles,nsuffix))
allocate (filename(nfiles),weight(nfiles,nsuffix))
allocate (dens(0:maxpart))
ncrit = min(maxpart+1,1000)
allocate (p_l(0:ncrit))

write (*,'(a)') ' -----------------------------------------------------------------'
write (*,'(a)') ' Optimization of critical parameters - version 3.0 (AZP, Apr. 24)'
write (*,'(a)') ' -----------------------------------------------------------------'

! -----------------------------------------------------------------
! Do Loop with no indent to deal with multiple passes over nsuffix
! -----------------------------------------------------------------
do ipass = 0,nsuffix ! ipass = 0 does all of them
counts = 0
entropy = 0.	! program looks at counts, so a zero value is not taken seriously
				! unless backed up by data
entropy(0,0:maxpart) = min_ene
entropy(-1,0:maxpart) = min_ene

open (1,file='input_hs2.dat')
read (1,*)
do isuffix = 1,nsuffix
    read (1,*)
enddo
read (1,*)
do ifile = 1,nfiles
    do isuffix = 1,nsuffix
        read (1,*) filename(ifile), counts_file(ifile,isuffix), &
                   weight(ifile,isuffix),beta(ifile,isuffix),mu(ifile,isuffix),width
        beta(ifile,isuffix) = 1./beta(ifile,isuffix)
        if (.not.(ifile==1 .and. isuffix==1)) weight(ifile,isuffix) = weight(ifile,isuffix) - weight(1,1)
    enddo
enddo
weight(1,1) = 0.
close (1)

olddev = 1e6
min_part = maxpart+1	! so that these are later set to correct values
max_part = -1	! so that these are later set to correct values

do ifile = 1,nfiles

    if (ipass == 0) then
        suff_start = 1; suff_end = nsuffix
    else
        suff_start = ipass ; suff_end = ipass
    endif
	do isuffix = suff_start,suff_end
		open (23,file='his'//trim(filename(ifile))//trim(suffix(isuffix))//'.dat')
		read (23,*)
		read (23,*) beta1,mu1,width1,xdim1,ydim1,zdim1	! although not actually used, these need to 
														! be read so that the correct number of lines are
														! skipped in the input file.
		min_part1 = maxpart + 1
		max_part1 = - 1
		iostat = 0
		omega = 0.
		omega(0,0:maxpart) = min_ene
		omega(-1,0:maxpart) = min_ene
		counts = 0
		
		do while (iostat/=-1)
			read (23,*,iostat=iostat) ipart,ncount,energy_start
			if (iostat/=-1) then
				if (ipart>maxpart.or.ncount>npoints) then
					write (*,*) 'Maxpart or Npoints exceeded, file, ipart, ncount = ',trim(filename(ifile))//trim(suffix(isuffix)),ipart,ncount
					stop
				endif
				omega(0,ipart) = energy_start
				omega(-1,ipart) = energy_start + ncount*width
				read (23,*) (counts(i),i=1,ncount)
				do i = 1,ncount
					if (counts(i)>0) then
						num_P = 0.
						do jfile = 1,nfiles
							do jsuffix = 1,nsuffix
								! to keep exponentials within reason, all are referenced to beta(1,1) and
								! mu(1,1)
								zfactor =  + beta(1,1)*(omega(0,ipart)+i*width) - &
								   beta(1,1)*mu(1,1)*ipart &
								   - beta(jfile,jsuffix)*(omega(0,ipart)+i*width) + &
								   beta(jfile,jsuffix)*mu(jfile,jsuffix)*ipart - weight(jfile,jsuffix)
								if (counts_file(jfile,jsuffix)>0) num_P = num_P + exp( log(counts_file(jfile,jsuffix)) + zfactor )
							enddo
						enddo
						omega(i,ipart) = counts(i)/num_P
					endif
				enddo
				min_part = min(min_part,ipart)
				max_part = max(max_part,ipart)
				min_part1 = min(min_part1,ipart)
				max_part1 = max(max_part1,ipart)
			endif
		enddo
		close(23)
		
		if (ifile==1.and.isuffix==suff_start) then
			entropy = omega
		else ! patch current data with previous data 
			do ipart = min_part,max_part
				idiff = 0
				if (omega(0,ipart)/=min_ene.and.entropy(0,ipart)/=min_ene) idiff = nint (( entropy(0,ipart) - omega(0,ipart)  ) / width)
				if ( idiff > 0) then	! shift entropy
					if ((max(entropy(-1,ipart),omega(-1,ipart))-omega(0,ipart))/width <= npoints) then
						entropy(0,ipart) = omega(0,ipart)
						entropy(-1,ipart) = max(entropy(-1,ipart),omega(-1,ipart))
						do i = npoints,1+idiff,-1
							entropy(i,ipart) = entropy(i-idiff,ipart)
						enddo
						do i = 1,1+idiff-1
							entropy(i,ipart) = 0.
						enddo
					else
						write (*,*) ' Width exceeded (low), ipart,file = ',ipart,trim(filename(ifile))//trim(suffix(isuffix))
						stop
					endif
				else if (idiff < 0) then
					if ((max(entropy(-1,ipart),omega(-1,ipart))-entropy(0,ipart))/width <= npoints) then
						entropy(-1,ipart) = max(entropy(-1,ipart),omega(-1,ipart))
						omega(0,ipart) = entropy(0,ipart)
						do i = npoints,1-idiff,-1
							omega(i,ipart) = omega(i+idiff,ipart)
						enddo
						do i = 1,1-idiff-1
							omega(i,ipart) = 0.
						enddo
					else
						write (*,*) ' Width exceeded (high), ipart,file = ',ipart,trim(filename(ifile))//trim(suffix(isuffix))
						stop
					endif
				endif
			enddo	
		
			do ipart = min_part,max_part
				do i = -1,0
					if (entropy(i,ipart) == min_ene) then
						entropy(i,ipart) = omega(i,ipart)
					endif
				enddo
				do i = 1,npoints
					entropy(i,ipart) = entropy(i,ipart) + omega(i,ipart)
				enddo
			enddo						

		endif

	enddo	! over isuffix
enddo	! over ifile

if (ipass == 0) then
    write (*,*) ' Enter Tstart, Tend, Tincr, Mu, smix'
    read (*,*) Tstart,Tend,Tincr,mu0,smix0
    write (*,*) ' Perturbations:      Mu           Smix      <N>       Dev'
endif

mu1 = mu0 ; smix = smix0
beta1 = 1./Tstart
p(1,1) = mu1
p(1,2) = smix
ptry(1:nopt_dim) = p(1,1:nopt_dim)
y(1) = funk(ptry)
if (ipass == 0) write(*,'(18x,g13.6,f9.4,f10.2,f10.3)') ptry(1),ptry(2),avg_N,y(1)
p(2,1) = mu1+trial_vec(1)
p(2,2) = smix
ptry(1:nopt_dim) = p(2,1:nopt_dim)
y(2) = funk(ptry)
if (ipass == 0) write(*,'(18x,g13.6,f9.4,f10.2,f10.3)') ptry(1),ptry(2),avg_N,y(2)
p(3,1) = mu1
p(3,2) = smix + trial_vec(2)
ptry(1:nopt_dim) = p(3,1:nopt_dim)
y(3) = funk(ptry)
if (ipass == 0) then
    write(*,'(18x,g13.6,f9.4,f10.2,f10.3)') ptry(1),ptry(2),avg_N,y(3)
    write (*,*) '-----------------------------------------------------'
    write (*,*) '    T         Dev        Mu      Smix        <N>    '
    write (*,*) '-----------------------------------------------------'
endif
if (ipass>0) write (*,'(2A)') '   Suffix = ', suffix(ipass)

Ntemp = int((Tend - Tstart)/Tincr) + 1
if (ipass == 0) then
    open (26,file='crit0.dat')
else
    open (26,file='crit_'//trim(suffix(ipass))//'.dat')
endif
do itemp = 1,Ntemp

	if (itemp /= 1) then	! Set initial amoeba for new temperature
		beta1 = 1./(Tstart + (itemp-1)*Tincr)
		p(1:3,1) = p(1:3,1) + dmu_dT*Tincr	! set chemical potential to a more reasonable value
		ptry(1:nopt_dim) = p(1,1:nopt_dim)
		y(1) = funk(ptry)
		p(2,1) = p(2,1) + trial_vec(1)/3. ! divide by 3 since already near optimum
		ptry(1:nopt_dim) = p(2,1:nopt_dim)
		y(2) = funk(ptry)
		p(3,2) = p(3,2) + trial_vec(2)/3.
		ptry(1:nopt_dim) = p(3,1:nopt_dim)
		y(3) = funk(ptry)
	endif

	call amoeba(p,y,nopt_dim+1,nopt_dim,nopt_dim,ftol,funk,iter)

	mu1 = p(1,1)
	smix = p(1,2)
	ptry(1:nopt_dim) = p(1,1:nopt_dim)
	y(1) = funk(ptry)
	Z = 0.
	avgnum = 0.
	ene_liq = 0.	
	dens = 0.
	do ipart = min_part,max_part
		do i = 1,npoints
			zfactor =    + beta(1,1)*(entropy(0,ipart)+i*width) - &
					   beta(1,1)*mu(1,1)*ipart &
			 		 - beta1*(entropy(0,ipart)+i*width)+beta1*mu1*ipart
			if (entropy(i,ipart)>0) Z = Z + exp( log(entropy(i,ipart)) + zfactor )
		enddo
	enddo
	do ipart = min_part,max_part
		do i = 1,npoints
			zfactor = + beta(1,1)*(entropy(0,ipart)+i*width) - &
					   beta(1,1)*mu(1,1)*ipart &
			 	  - beta1*(entropy(0,ipart)+i*width)+beta1*mu1*ipart
			probab = 0.
			if (entropy(i,ipart)>0) probab = exp(log(entropy(i,ipart)) + zfactor)
			dens(ipart) = dens(ipart) + probab/Z
			ene_liq = ene_liq + (entropy(0,ipart)+i*width)*probab/Z
		enddo
		avgnum = avgnum + ipart * dens(ipart)
	enddo

	write (*,'(g13.5,f6.4,f11.4,g12.3,f9.2)') 1/beta1,y(1),mu1,smix,avgnum
	if (y(1)<olddev) then
        if (ipass == 0) then
            Tall = 1/beta1 ; dall = y(1) ; muall = mu1; Sall = smix; Nall = avgnum
        else
            Topt(ipass)=1/beta1 ; dopt(ipass) = y(1) ; muopt(ipass) = mu1
            Sopt(ipass) = smix; Nopt(ipass) = avgnum
        endif
		olddev=y(1)
		rewind(26)
		write (26,'(A,f8.4,A,f8.4,A,20A5,3h *//3h /*,50A,3h *//3h /*,50A,3h */)') '/* T=',1/beta1,'; mu=',mu1,'; Suffixes and Files = '&
			,(suffix(isuffix),isuffix=1,nsuffix),(trim(filename(ifile)),ifile=1,nfiles),' */'
		p_l = 0.
		avg_x = 0.
		min_x = min_part - smix*entropy(-1,min_part)	! we are looking for extreme values
		max_x = min_x
		min_x = min(min_x, min_part - smix*entropy(0,min_part))
		max_x = max(max_x, min_part - smix*entropy(0,min_part))
		min_x = min(min_x, min_part - smix*(entropy(0,min_part)+npoints*width))
		max_x = max(max_x, min_part - smix*(entropy(0,min_part)+npoints*width))
		min_x = min(min_x, max_part - smix*entropy(0,max_part))
		max_x = max(max_x, max_part - smix*entropy(0,max_part))
		min_x = min(min_x, max_part - smix*(entropy(0,min_part)+npoints*width))
		max_x = max(max_x, max_part - smix*(entropy(0,min_part)+npoints*width))
        dx = (max_x - min_x) / ncrit
		do ipart = min_part,max_part
			if (entropy(0,ipart) /= min_ene) then	! to avoid problems if no observations at given ipart
				do i = 1,npoints
					x = ipart - smix*(entropy(0,ipart)+i*width)
					zfactor = + beta(1,1)*(entropy(0,ipart)+i*width) - &
							   beta(1,1)*mu(1,1)*ipart &
			 			  - beta1*(entropy(0,ipart)+i*width)+beta1*mu1*ipart
					probab = 0.
					if (entropy(i,ipart)>0) probab = exp(log(entropy(i,ipart))+zfactor)
					index = nint(ncrit*(x-min_x)/(max_x-min_x))
					p_l(index) = p_l(index) + probab/Z
					avg_x = avg_x + x * probab/Z
				enddo
			endif
		enddo
		variance = 0
		do i = 0,ncrit
			variance = variance + p_l(i)*(min_x+i*(max_x-min_x)/ncrit-avg_x)**2
		enddo

		do i = 0,ncrit
            write (26,'(2g12.5)') (min_x+i*(max_x-min_x)/ncrit-avg_x)/sqrt(variance), &
							   p_l(i)*sqrt(variance)/(max_x-min_x)*ncrit
		enddo
		write (26,'(A,2g14.3,A)') '/* <x> =',avg_x,variance,' */'
	endif
enddo

if (ipass>0) then
    write (*,*) '-----------------------------------------------------------------'
    write (*,'(g13.5,f6.4,f11.4,g12.3,f9.2)') &
        Topt(ipass),dopt(ipass),Muopt(ipass),Sopt(ipass),Nopt(ipass)
endif

enddo ! end of non-indented loop over ipass

if (nsuffix>1) then
    call stats(nsuffix,Topt,Tstd); call stats(nsuffix,Muopt,Mustd);call stats(nsuffix,Sopt,Sstd)
    call stats(nsuffix,Nopt,Nstd); call stats(nsuffix,dopt,dstd)

    write (*,*) '------------------------------------------------------------------------'
    write (*,*) '   T      ±        μ      ±     dev     ±      s      ±      Nc     ±'
    write (*,'(  f8.4,   f7.4,   f9.4,   f7.4   f6.3,  f6.3,  f8.4,  f7.4,   f8.2,  f7.2)') &
            Tall,Tstd,Muall,Mustd,dall,dstd,Sall,Sstd,Nall,Nstd
endif

stop

end

real function funk(ptry)
use data
integer, parameter :: nopt_dim = 2  ! number of optimizing directions: Mu_c, smix
real ptry(nopt_dim)
real Z					! partition function
integer ipart,i	! local counters
real x						! ordering operator distribution x axis
real max_x, min_x, avg_x	! limits of x, average x
real xpoint, ypoint, yvalue
integer inonzero	! number of actual points in distribution
real zfactor
real probab,variance,fvalue
integer index

mu1 = ptry(1)
smix = ptry(2)
Z = 0.
do ipart = min_part,max_part
	do i = 1,npoints
		zfactor = + beta(1,1)*(entropy(0,ipart)+i*width) - &
					beta(1,1)*mu(1,1)*ipart &
 				  - beta1*(entropy(0,ipart)+i*width)+beta1*mu1*ipart
		if (entropy(i,ipart) > 0) Z = Z +  exp( log(entropy(i,ipart)) + zfactor )
	enddo
enddo
p_l = 0.
avg_x = 0.
avg_N = 0.
min_x = min_part - smix*entropy(-1,min_part)    ! we are looking for extreme values
max_x = min_x
min_x = min(min_x, min_part - smix*entropy(0,min_part))
max_x = max(max_x, min_part - smix*entropy(0,min_part))
min_x = min(min_x, min_part - smix*(entropy(0,min_part)+npoints*width))
max_x = max(max_x, min_part - smix*(entropy(0,min_part)+npoints*width))
min_x = min(min_x, max_part - smix*entropy(0,max_part))
max_x = max(max_x, max_part - smix*entropy(0,max_part))
min_x = min(min_x, max_part - smix*(entropy(0,max_part)+npoints*width))
max_x = max(max_x, max_part - smix*(entropy(0,max_part)+npoints*width))
dx = (max_x - min_x) / ncrit
do ipart = min_part,max_part
	do i = 1,npoints
		x = ipart - smix*(entropy(0,ipart)+i*width)
		zfactor = + beta(1,1)*(entropy(0,ipart)+i*width) - &
					beta(1,1)*mu(1,1)*ipart &
				  - beta1*(entropy(0,ipart)+i*width)+beta1*mu1*ipart
		probab = 0.
		if (entropy(i,ipart)>0) probab = exp(log(entropy(i,ipart))+zfactor)
		index = nint((x-min_x)/dx)
        if (index<=ncrit .and. index >= 0) then
            p_l(index) = p_l(index) + probab/Z
            avg_N = avg_N + ipart * probab/Z
        endif
	enddo
enddo

do i = 0,ncrit
    avg_x = avg_x + p_l(i)*i
enddo

variance = 0
do i = 0,ncrit
    variance = variance + p_l(i)*(i-avg_x)**2
enddo

fvalue = 0.
inonzero = 0
do i = 0,ncrit
    xpoint = (i-avg_x)/sqrt(variance)
	ypoint = p_l(i)*sqrt(variance)
    yvalue = exp(-(0.7774*xpoint*xpoint-1)**2*(0.158*0.7774*xpoint*xpoint+0.776))/2.3306
    fvalue = fvalue + yvalue*abs(ypoint - yvalue)
!    write(*,*) i,xpoint,ypoint,p_l(i),val
	if (ypoint > 1e-6 .and. yvalue > 1e-6) then
		inonzero = inonzero + 1
	endif
enddo
fvalue = fvalue/inonzero
funk = 100*fvalue
!write (*,'(2f12.5,f10.3,f13.7,1x,f10.3)') 1./beta1,mu1,avg_N,smix,100.*fvalue

return
end

      FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac)
      INTEGER ihi,mp,ndim,np,NMAX
      REAL amotry,fac,p(mp,np),psum(np),y(mp),funk
      PARAMETER (NMAX=2)
      EXTERNAL funk  !    USES funk
      INTEGER j
      REAL fac1,fac2,ytry,ptry(NMAX)

      fac1=(1.-fac)/ndim
      fac2=fac1-fac
      do 11 j=1,ndim
        ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11    continue
      ytry=funk(ptry)
      if (ytry.lt.y(ihi)) then
        y(ihi)=ytry
        do 12 j=1,ndim
          psum(j)=psum(j)-p(ihi,j)+ptry(j)
          p(ihi,j)=ptry(j)

12      continue
      endif
      amotry=ytry
!	  write(*,'(2g13.5,g13.6,f10.2)') ytry,ptry(1),ptry(2),ptry(3)
      return
      END

SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter)
INTEGER iter,mp,ndim,np,NMAX,ITMAX
REAL ftol,p(mp,np),y(mp),funk
PARAMETER (NMAX=2,ITMAX=5000)
EXTERNAL funk

!     USES amotry,funk
      INTEGER i,ihi,ilo,inhi,j,m,n
      REAL rtol,sum,swap,ysave,ytry,psum(NMAX),amotry
      iter=0
1     do 12 n=1,ndim
        sum=0.
        do 11 m=1,ndim+1
          sum=sum+p(m,n)
11      continue
        psum(n)=sum
12    continue
2     ilo=1
      if (y(1).gt.y(2)) then
        ihi=1
        inhi=2
      else
        ihi=2

        inhi=1
      endif
      do 13 i=1,ndim+1
        if(y(i).le.y(ilo)) ilo=i
        if(y(i).gt.y(ihi)) then
          inhi=ihi
          ihi=i
        else if(y(i).gt.y(inhi)) then
          if(i.ne.ihi) inhi=i
        endif
13    continue
      rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
      if (rtol.lt.ftol) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue

        return
      endif
      if (iter.ge.ITMAX) write (*,*) 'ITMAX exceeded in amoeba'
      iter=iter+2
      ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0)
      if (ytry.le.y(ilo)) then
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0)
      else if (ytry.ge.y(inhi)) then
        ysave=y(ihi)
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5)
        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
              y(i)=funk(psum)
            endif
16        continue
          iter=iter+ndim
          goto 1

        endif
      else
        iter=iter-1
      endif
      goto 2
      END

subroutine stats (N,X,Std)

! This routine calculates the standard deviation of a sample

IMPLICIT NONE
INTEGER                 :: N   ! Number of points in sample (input)
REAL*8, DIMENSION(N)   :: X    ! vector containing sample values (input)
REAL*8                 :: Xavg ! mean value    (output)
REAL*8                 :: Std    ! standard deviation (output)
INTEGER              :: I    ! index (auxiliary variable)

Std = 0 ; Xavg = 0
DO I = 1, N
  Xavg = Xavg + X(I)
END DO

Xavg = Xavg / N

DO I = 1, N
  Std = Std + (X(I) - Xavg) ** 2
END DO

Std = SQRT( Std / (N - 1) )

return
end
