SUBROUTINE powellann(p,xi,n,np,ftol,fret,pmin,pmax,
     &   funkmin,f1dim,ITMAX)
! fret must be given on entry

      implicit none
      INTEGER iter,n,np,NMAX,ITMAX
      double precision fret,ftol,p(np),xi(np,np),TINY,
     &        pmin(np),pmax(np)

      PARAMETER (NMAX=50,TINY=1.0d-25)
CU    USES funkmin,annlinmin
      INTEGER i,ibig,j
      double precision del,fp,fptt,t,pt(NMAX),
     &    ptt(NMAX),xit(NMAX)
      external funkmin,f1dim
           
      do 11 j=1,n
        pt(j)=p(j)
11    continue
      iter=0
1     iter=iter+1
      fp=fret
      ibig=0
      del=0.0d0
      do 13 i=1,n
        do 12 j=1,n
          xit(j)=xi(j,i)
12      continue
        fptt=fret
        
        call annlinmin(p,pmin,pmax,xit,n,f1dim,fret)
        if(fptt-fret.gt.del)then
          del=fptt-fret
          ibig=i
        endif
13    continue
      if(2.0d0*(fp-fret).le.ftol*(dabs(fp)+dabs(fret))+TINY)return
      if(iter.eq.ITMAX)then
        write(*,*)'powell exceeding maximum iterations'
        return
      endif
      do 14 j=1,n
        ptt(j)=2.0d0*p(j)-pt(j)
        xit(j)=p(j)-pt(j)
        pt(j)=p(j)
14    continue
      call funkmin(n,ptt,fptt)
      if(fptt.ge.fp)goto 1
      t=2.0d0*(fp-2.0d0*fret+fptt)*(fp-fret-del)**2-
     &    del*(fp-fptt)**2
      if(t.ge.0.0d0)goto 1
      
      call annlinmin(p,pmin,pmax,xit,n,f1dim,fret)
            
      do 15 j=1,n
        xi(j,ibig)=xi(j,n)
        xi(j,n)=xit(j)
15    continue
      goto 1
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.

      SUBROUTINE annlinmin(p,pmin,pmax,xi,n,f1dim,fret)
      implicit none
      INTEGER n,NMAX
      double precision fret,p(n),xi(n),ftol,pmin(n),pmax(n)
      PARAMETER (NMAX=1000,ftol=1.0d-6)
CU    USES brent,f1dim,mnbrak
      INTEGER j,ncom,k,ierr
      double precision ax,bx,fa,fb,fx,xmin,xx,pcom(NMAX),
     &   xicom(NMAX),xxmin,xxmax,term,w(2),bph(2),
     &   paramnormsk,paramnormsb,ynorm,q(2),bend,terma,
     &   termb,termc,delta,reducer,root,f1dim,postannfunc
      parameter(reducer=0.25d0) 
     
      COMMON /f1com/ pcom,xicom,ncom
      
      integer maxm1dsamp,m1dsamp
      parameter(maxm1dsamp=1000)
      common /intannlinmin/m1dsamp
      double precision y1dsamp(maxm1dsamp),params1d(maxm1dsamp)
      common /dbleannlinmin/y1dsamp,params1d

      EXTERNAL f1dim

      ncom=n
      delta=0.5d0
      
      do j=1,n
        pcom(j)=p(j)
        xicom(j)=xi(j)
      enddo

      xxmax=1.0d+100
      xxmin=-1.0d+100
      do j=1,n
        if(xicom(j).gt.1.0d-100)then
!        if(xicom(j).gt.0.0d0)then
          xx=(pmax(j)-pcom(j))/xicom(j)
          ax=(pmin(j)-pcom(j))/xicom(j)
        else
          if(xicom(j).lt.(-1.0d-100))then
!          if(xicom(j).lt.0.0d0)then
            ax=(pmax(j)-pcom(j))/xicom(j)
            xx=(pmin(j)-pcom(j))/xicom(j)
          else
            xx=1.0d+100
            ax=-1.0d+100
          endif
        endif
        if(xxmax.gt.xx)then
          xxmax=xx
        endif
        if(xxmin.lt.ax)then
          xxmin=ax
        endif
      enddo
      ax=0.0d0
      if(dabs(xxmax).gt.dabs(xxmin))then
        xx=0.25d0*xxmax
      else
        xx=0.25d0*xxmin  
      endif
      
      m1dsamp=0
      call mnbrak(ax,xx,bx,fa,fx,fb,xxmin,xxmax,ierr,f1dim)
      if(ierr.ne.0)then
! bracketing failed
        fret=fx
        do j=1,n
          xi(j)=xx*xi(j)
          p(j)=p(j)+xi(j)
        enddo
        return
      endif
      
      if(ax.gt.xx)then
        xxmax=ax
        xxmin=bx
      else
        xxmax=bx
        xxmin=ax
      endif
      
      write(*,*)m1dsamp
                 
      if(m1dsamp.lt.6)then
        term=xx+(xxmax-xx)*(delta+1.0d0)/2.0d0       
        fret=f1dim(term)
        if(fret.lt.fx)then
          xxmin=xx
          xx=term
          fx=fret
        else
          xxmax=term
        endif
        term=xx+(xxmax-xx)*(1.0d0-delta)/2.0d0       
        fret=f1dim(term)
        if(fret.lt.fx)then
          xxmin=xx
          xx=term
          fx=fret
        else
          xxmax=term
        endif
        term=xx-(xx-xxmin)*(1.0d0+delta)/2.0d0         
        fret=f1dim(term)
        if(fret.lt.fx)then
          xxmax=xx
          xx=term
          fx=fret
        else
          xxmin=term
        endif
        term=xx-(xx-xxmin)*(1.0d0-delta)/2.0d0         
        fret=f1dim(term)
        if(fret.lt.fx)then
          xxmax=xx
          xx=term
          fx=fret
        else
          xxmin=term
        endif
      endif
            
      w(1)=0.75334d0
      w(2)=0.13425d0
      bph(1)=0.01d0
      bph(2)=-0.07d0
      q(1)=1.2d0
      q(2)=-2.0d0
      
100   call annfitting1d(m1dsamp,params1d,y1dsamp,paramnormsk,
     &  paramnormsb,ynorm,w,bph,q,bend)
     
      
      term=q(1)*w(1)
      fret=q(2)*w(2)
      terma=term*w(2)*w(2)+fret*w(1)*w(1)
      termb=2.0d0*(term*w(2)*bph(2)+fret*w(1)*bph(1))
      termc=term*(1.0d0+bph(2)*bph(2))+
     &      fret*(1.0d0+bph(1)*bph(1))
      
      if(terma.eq.0.0d0)then
        if(termb.eq.0.0d0)goto 200
        root=-termc/termb
        goto 110
      endif
      fret=termb*termb-4.0d0*terma*termc      
      if(fret.lt.0.0d0)goto 200
      term=-0.5d0*(termb+dsign(1.0d0,termb)*dsqrt(fret))
      root=termc/term
      root=(root-paramnormsb)/paramnormsk
      if(root.lt.xxmax.and.root.gt.xxmin)then
        if(dabs(root-xx).lt.ftol)goto 1000
        fret=f1dim(root)
        if(fret.lt.fx)then          
          if(root.lt.xx)then
            xxmax=xx
          else
            xxmin=xx
          endif
          xx=root
          fx=fret
          goto 100
        else
          if(root.lt.xx)then
            xxmin=root
          else
            xxmax=root
          endif
        endif
      endif 
      root=term/terma 
110   root=(root-paramnormsb)/paramnormsk   
      if(root.lt.xxmax.and.root.gt.xxmin)then
        if(dabs(root-xx).lt.ftol)goto 1000
        fret=f1dim(root)        
        if(fret.lt.fx)then
          if(root.lt.xx)then
            xxmax=xx
          else
            xxmin=xx
          endif
          xx=root
          fx=fret
          goto 100
        else
          if(root.lt.xx)then
            xxmin=root
          else
            xxmax=root
          endif
        endif
      endif 

200   term=xx+delta*(xxmax-xx)
      
      write(*,*)'ann failed'
      
      fret=f1dim(term)
      if(fret.lt.fx)then
        xxmin=xx
        fx=fret
        xx=term
        goto 100
      endif
      xxmax=term
      term=xx-delta*(xx-xxmin)
      fret=f1dim(term)
      if(fret.lt.fx)then
        xxmax=xx
        fx=fret
        xx=term
        goto 100
      endif
      xxmin=term
      delta=delta*reducer
      if(delta.gt.ftol)goto 200
      
1000  fret=fx
      do j=1,n
        xi(j)=xx*xi(j)
        p(j)=p(j)+xi(j)
      enddo

      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.

      SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,xxmin,xxmax,
     &                  ierr,func)
      implicit none
      double precision ax,bx,cx,fa,fb,fc,
     &   func,GOLD,GLIMIT,TINY
      EXTERNAL func
      PARAMETER(GOLD=1.618034d0,GLIMIT=100.0d0,TINY=1.0d-20)
      double precision dum,fu,q,r,u,ulim,xxmin,xxmax
      integer ierr
      ierr=0
      fa=func(ax)
      fb=func(bx)
      if(fb.gt.fa)then
        dum=ax
        ax=bx
        bx=dum
        dum=fb
        fb=fa
        fa=dum
      endif
      if(fa.eq.fb)then
        cx=(bx+ax)/2.0d0
        fc=func(cx)
        if(fc.le.fa)return
      endif
      cx=bx+GOLD*(bx-ax)
      if(cx.le.xxmin)then
        cx=0.5d0*(dmin1(ax,bx)+xxmin)
      endif
      if(cx.ge.xxmax)then
        cx=0.5d0*(dmax1(ax,bx)+xxmax)
      endif      
      fc=func(cx)      
1     if(fb.ge.fc)then
        r=(bx-ax)*(fb-fc)
        q=(bx-cx)*(fb-fa)
        u=bx-((bx-cx)*q-(bx-ax)*r)/
     &    (2.0d0*dsign(dmax1(dabs(q-r),TINY),q-r))
        ulim=bx+GLIMIT*(cx-bx)
        if(ulim.ge.xxmax)then
          ulim=xxmax-tiny
        endif
        if(ulim.le.xxmin)then
          ulim=xxmin+tiny
        endif
        if((bx-u)*(u-cx).gt.0.0d0)then
          fu=func(u)
          if(fu.lt.fc)then
            ax=bx
            fa=fb
            bx=u
            fb=fu
            return
          elseif(fu.gt.fb)then
            cx=u
            fc=fu
            return
          endif
          u=cx+GOLD*(cx-bx)
          if(u.gt.xxmax)then
            u=cx+0.5d0*(xxmax-cx)
          endif
          if(u.lt.xxmin)then
            u=cx+0.5d0*(xxmin-cx)
          endif
          fu=func(u)
        elseif((cx-u)*(u-ulim).gt.0.0d0)then
          fu=func(u)
          if(fu.lt.fc)then
            bx=cx
            cx=u
            u=cx+GOLD*(cx-bx)
            if(u.gt.xxmax)then
              u=cx+0.5d0*(xxmax-cx)
            endif
            if(u.lt.xxmin)then
              u=cx+0.5d0*(xxmin-cx)
            endif
            fb=fc
            fc=fu
            fu=func(u)
          endif
        else if((u-ulim)*(ulim-cx).ge.0.0d0)then
          u=ulim
          fu=func(u)
        else
          u=cx+GOLD*(cx-bx)
          if(u.gt.xxmax)then
            u=cx+0.5d0*(xxmax-cx)
          endif
          if(u.lt.xxmin)then
            u=cx+0.5d0*(xxmin-cx)
          endif
          fu=func(u)
        endif
        ax=bx
        bx=cx
        cx=u
        fa=fb
        fb=fc
        fc=fu
        r=dmin1(dabs(ax-bx),dabs(ax-cx))
        r=dmin1(r,dabs(bx-cx))
        if(r.lt.tiny)then
! bracketing failed
          ierr=1
          return
        endif
        goto 1
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.