subroutine simplex(x,n,nmax,f,epsf,epsx,ff,x0,xr,xe,xc)
external f ; dimension x(nmax,nmax+1),ff(n+1),x0(n),xr(n),xe(n),xc(n)
logical converged

do k=1,n+1 ; ff(k)=f(x(1,k),n) ; end do
call hilo(ff,n,khi,klo)
call centroid(x,n,nmax,khi,x0)

do while(.not. converged(x,n,nmax,ff,khi,klo,x0,epsf,epsx))
call reflection(x,n,nmax,khi,x0,xr)
fr=f(xr,n)
if (fr<ff(klo)) then                        !very good: expand and accept xe
	call expansion(x0,xr,xe,n)
	ff(khi)=f(xe,n)
	do i=1,n ; x(i,khi)=xe(i) ; end do
	print *,'expansion'
else if (fr<ff(khi)) then                   !just good: accept xr
	ff(khi)=f(xr,n)
	do i=1,n ; x(i,khi)=xr(i) ; end do
	print *,'reflection'
else 
	call contraction(x,n,nmax,khi,x0,xc)
	fc=f(xc,n)
	if (fc<ff(khi)) then                 !accept xc
		ff(khi)=fc
		do i=1,n ; x(i,khi)=xc(i) ; end do
		print *,'contraction'
	else
		call reduction(x,n,nmax,klo)
		do k=1,n+1 ; ff(k)=f(x(1,k),n) ; end do
		print *,'reduction'
	end if
end if
call hilo(ff,n,khi,klo)
call centroid(x,n,nmax,khi,x0)
end do
print *,'x0=',x0
print *,'fhi,flo=',ff(khi),ff(klo)
return; end

subroutine reflection(x,n,nmax,khi,x0,xr)
dimension x(nmax,nmax+1),x0(n),xr(n) ; parameter(alpha=1)
do i=1,n
	xr(i)=(1+alpha)*x0(i)-alpha*x(i,khi)
end do
return ; end

subroutine expansion(x0,xr,xe,n)
dimension x0(n),xr(n),xe(n) ; parameter(beta=2)
do i=1,n; xe(i)=(1-beta)*x0(i)+beta*xr(i); end do
return ; end

subroutine contraction(x,n,nmax,khi,x0,xc)
dimension x(nmax,nmax+1),x0(n),xc(n) ; parameter(gamma=.5)
do i=1,n
	xc(i)=x0(i)+gamma*(x(i,khi)-x0(i))
end do
return ; end

subroutine reduction(x,n,nmax,klo)
dimension x(nmax,nmax+1)
do k=1,n+1
do i=1,n ; x(i,k)=(x(i,k)+x(i,klo))/2 ; end do
end do
return ; end

subroutine hilo(ff,n,khi,klo)
dimension ff(n+1)                  ! function at the vertices
khi=1 ; klo=1
fhi=ff(khi) ; flo=ff(klo)
do k=2,n+1
	if (ff(k)>fhi) then
		fhi=ff(k); khi=k
	end if
	if (ff(k)<flo) then
		flo=ff(k); klo=k
	end if
end do
return ; end

subroutine centroid(x,n,nmax,khi,x0)
dimension x(nmax,nmax+1),x0(n)         ! simplex and centroid
do i=1,n
	sum=0
	do k=1,n+1
		sum=sum+x(i,k)    
	end do
	x0(i)=(sum-x(i,khi))/n
end do
return ; end

logical function converged(x,n,nmax,ff,khi,klo,x0,epsf,epsx)
dimension x(nmax,nmax+1),ff(n)

size=distance(x(1,1),x0,n)
do k=2,n+1
	dx=distance(x(1,k),x0,n)
	if (dx>size) size=dx
end do
converged=.false.
if(ff(khi)-ff(klo) < epsf) then
	converged=.true.
	print *,'f_high-f_low < ',epsf
endif
if(size < epsx) then
	converged=.true.
	print *,'simplex size < ',epsx
endif
return ; end

function distance(x,y,n)
dimension x(n),y(n)
d=0
do i=1,n
	d=d+(x(i)-y(i))**2
end do
distance=sqrt(d)
return;end