subroutine surffunc(nyvars,yvars,nxvars,
     &  xvars,ndim,beta,dydxp,idowhat)
      implicit none
!idowhat=0, value of the function only
!idowhat=1, derivative with respect to the independent variable x
!idowhat=2, derivative with respect to the parameters
      integer nyvars,nxvars,ndim,idowhat
      double precision yvars(nyvars),xvars(nxvars),
     &beta(5),dydxp(nyvars,(nxvars+5))
!-------------------------------------------------------
      double precision y0,a,b,c,x0,x,term,crit
      parameter(crit=300.0d0)
      a=beta(1)
      b=beta(2)
      c=beta(3)
      x=xvars(1)
      if(ndim.eq.3)then
        term=dexp(-b*x)
        if(idowhat.eq.0)yvars(1)=c+a*(1.0d0-term)
        if(idowhat.eq.1)then
          dydxp(1,1)=a*b*term
        endif
        if(idowhat.eq.2)then
          dydxp(1,1)=1.0d0-term
          dydxp(1,2)=a*x*term
          dydxp(1,3)=1.0d0
        endif
        return
      endif
      if(ndim.eq.4)then
        x0=beta(4)
        if(idowhat.eq.0)yvars(1)=a*(1.0d0-b*x)*(x-x0)/(1.0d0+c*x)
        if(idowhat.eq.1)then
          dydxp(1,1)=a*(1.0d0-2.0d0*b*x-b*c*x*x+(b+c)*x0)/
     &((1.0d0+c*x)*(1.0d0+c*x))
        endif
        if(idowhat.eq.2)then
          dydxp(1,1)=(1.0d0-b*x)*(x-x0)/(1.0d0+c*x)
          dydxp(1,2)=a*(1.0d0-x)*(x-x0)/(1.0d0+c*x)
          dydxp(1,3)=-a*(1.0d0-b*x)*(x-x0)*x/((1.0d0+c*x)*(1.0d0+c*x))
          dydxp(1,4)=-a*(1.0d0-b*x)/(1.0d0+c*x)
        endif
        return
      endif
!      if(ndim.eq.3)then
!        yvars(1)=(1.0d0+a*x)/(b+c*x)
!        if(idowhat.eq.0)return
!        if(idowhat.eq.1)then
!          dydxp(1,1)=(a-yvars(1)*c)/(b+c*x)
!          return
!        endif
!        if(idowhat.eq.2)then
!          dydxp(1,1)=x/(b+c*x)
!          dydxp(1,2)=-yvars(1)/(b+c*x)
!          dydxp(1,3)=-yvars(1)*dydxp(1,1)
!          return
!        endif
!      endif

!A/Ci or A/PAR curves
      x0=beta(4)
      y0=beta(5)
      if(idowhat.eq.0)then
        if((-(x-x0)/b).lt.crit)then
          term=dexp(-(x-x0)/b)
          yvars(1)=y0+a*(1.0d0/(1.0d0+term))**c
        else
          term=dexp((x-x0)/b)
          yvars(1)=y0+a*(term/(1.0d0+term))**c
        endif
      endif
      if(idowhat.eq.1)then
        if((-(x-x0)/b).lt.crit)then
          term=dexp(-(x-x0)/b)
          dydxp(1,1)=(a*c*term/b)*
     &      (1.0d0/(1.0d0+term))**(1.0d0+c)
        else
          term=(x-x0)/b
          dydxp(1,1)=(a*c/b)*(dexp(term*c/(c+1.0d0))/
     &      (1.0d0+dexp(term)))**(c+1.0d0)
        endif
      endif
      if(idowhat.eq.2)then
        dydxp(1,5)=1.0d0
        if((-(x-x0)/b).lt.crit)then
          term=dexp(-(x-x0)/b)
          dydxp(1,1)=(1.0d0/(1.0d0+term))**c
          dydxp(1,4)=-(a*c*term/b)*
     &      (1.0d0/(1.0d0+term))**(1.0d0+c)
          dydxp(1,2)=-(a*c*term*(x-x0)/(b*b))*
     &      (1.0d0/(1.0d0+term))**(1.0d0+c)
          dydxp(1,3)=-(a*dlog(1.0d0+term))*
     &      (1.0d0/(1.0d0+term))**c
        else
          term=(x-x0)/b
          dydxp(1,1)=(dexp(term)/(1.0d0+dexp(term)))**c
          dydxp(1,4)=-(a*c/b)*(dexp(term*c/(c+1))/
     &              (1.0d0+dexp(term)))**(c+1.0d0)
          dydxp(1,2)=-(a*c*(x-x0)/(b*b))*(dexp(term*c/
     &      (c+1.0d0))/(1.0d0+dexp(term)))**(1.0d0+c)
          dydxp(1,3)=-a*(dlog(1.0d0+dexp(term))-term)*
     &           (dexp(term)/(1.0d0+dexp(term)))**c
        endif
      endif
      return
      end
!==========================================================
      subroutine properties_surffunc(ndim,beta,root,der_root,fmax,
     &yinter,der_yinter,agivenx,der_agivenx,funcval_agivenx,
     &xmin,xmax,curvatmax,xcurvatmax)
      implicit none
      integer ndim
      double precision beta(5),root,der_root,fmax,yinter,der_yinter,
     &agivenx,der_agivenx,funcval_agivenx,xmin,xmax,curvatmax,xcurvatmax
      double precision a,b,c,x0,y0,term,term1,term2,term3,step,
     &deratx,der2atx,yvector(5),xvector(5),zvector(5)

      a=beta(1)
      b=beta(2)
      c=beta(3)

      if(ndim.eq.3)then
!y=c+a(1-exp(-bx))
        root=-dlog(1.0d0+c/a)/b
        der_root=a*b*dexp(-b*root)
        fmax=c+a
        yinter=c
        der_yinter=a*b
        funcval_agivenx=c+a*(1.0d0-dexp(-b*agivenx))
        der_agivenx=a*b*dexp(-b*agivenx)
        xcurvatmax=dlog(2*a*a*b*b)/(2.0d0*b)
!        curvatmax=-a*b*b*dexp(-b*xcurvatmax)/
!     &((1.0d0+a*a*b*b*dexp(-2.0d0*b*xcurvatmax))**1.5d0)
        curvatmax=-b*0.3849d0
        curvatmax=dabs(curvatmax)*1000.0d0
        return
      endif

      if(ndim.eq.4)then
!y=a*(1-bx)*(x-x0)/(1+c*x)
!we ignore the other root
!dydxp(1,1)=a*(1.0d0-2.0d0*b*x-b*c*x*x+(b+c)*x0)/((1.0d0+c*x)*(1.0d0+c*x))
        x0=beta(4)
        root=x0
        der_root=a*(1.0d0-2.0d0*b*root-b*c*root*root+(b+c)*x0)/
     &((1.0d0+c*root)*(1.0d0+c*root))
        term=(dsqrt((b+c)*(1.0d0+c*x0)/b)-1.0d0)/c
        fmax=a*(1.0d0-b*term)*(term-x0)/(1.0d0+c*term)
        yinter=-a*x0
        der_yinter=a*(1.0d0+(b+c)*x0)
        funcval_agivenx=a*(1.0d0-b*agivenx)*(agivenx-x0)/
     &(1.0d0+c*agivenx)
        der_agivenx=
     &a*(1.0d0-2.0d0*b*agivenx-b*c*agivenx*agivenx+(b+c)*x0)/
     &((1.0d0+c*agivenx)*(1.0d0+c*agivenx))
        xcurvatmax=-9999.0d0
        curvatmax=-9999.0d0
        return
      endif
!      if(ndim.eq.3)then
!y=(1+a*x)/(b+c*x)
!        root=-1.0d0/a
!        der_root=a/(b-c/a)
!        fmax=a/c
!        yinter=1.0d0/b
!        der_yinter=(a*b-c)/(b*b)
!        return
!      endif
      x0=beta(4)
      y0=beta(5)
      if((-a/y0).gt.0.0d0)then
        term=(-a/y0)**(1.0d0/c)-1.0d0
        root=x0-b*dlog(term)
        term=dexp(-(root-x0)/b)
        der_root=(a*c*term/b)*(1.0d0/(1.0d0+term))**(1.0d0+c)
      else
        root=-9999.0d0
        der_root=-9999.0d0
      endif
      fmax=a+y0
      xvector(1)=0.0d0
      call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),0)
      yinter=yvector(1)
      call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),1)
      der_yinter=zvector(1)
      xvector(1)=agivenx
      call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),1)
      der_agivenx=zvector(1)
      call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),0)
      funcval_agivenx=yvector(1)

      curvatmax=-9999.0d0
      xcurvatmax=-9999.0d0
      step=(xmax-xmin)/1000.0d0
      if(step.le.0.0d0)return
!      do term=xmin,xmax,step
!        call surffunc(1,term1,1,term,ndim,beta,deratx,1)
!        term2=dexp(-(term-x0)/b)
!        der2atx=-deratx/b+
!     &(1.0d0+c)*deratx*deratx*((1.0d0+term2)**c)/(a*c)
!        term3=dabs(der2atx/((1.0d0+deratx*deratx)**1.5d0))
!        if(term3.gt.curvatmax)then
!          curvatmax=term3
!          xcurvatmax=term
!        endif
!      enddo

      term=xmin
      do while (term.le.xmax)
        xvector(1)=term
        call surffunc(1,yvector(1),1,xvector(1),ndim,beta,zvector(1),1)
        term1=yvector(1)
        deratx=zvector(1)
        term2=dexp(-(term-x0)/b)
        der2atx=-deratx/b+
     &(1.0d0+c)*deratx*deratx*((1.0d0+term2)**c)/(a*c)
        term3=dabs(der2atx/((1.0d0+deratx*deratx)**1.5d0))
        if(term3.gt.curvatmax)then
          curvatmax=term3
          xcurvatmax=term
        endif
        term=term+step
      enddo

      if(dabs(xcurvatmax-xmin).le.step.or.
     &dabs(xcurvatmax-xmax).le.step)then
        curvatmax=-9999.0d0
        xcurvatmax=-9999.0d0
      else
        curvatmax=dabs(curvatmax)*1000.0d0
      endif
      return
      end
!==========================================================
      subroutine indices_surffunc(ndim,beta,root,
     &  der_root,fmax)
      implicit none
      integer ndim
      double precision beta(ndim),root,der_root,fmax
      double precision a,b,c,x0,y0,term
      a=beta(1)
      b=beta(2)
      c=beta(3)
      x0=beta(4)
      y0=beta(5)
      term=(-a/y0)**(1.0d0/c)-1.0d0
      root=x0-b*dlog(term)
      term=dexp(-(root-x0)/b)
      der_root=(a*c*term/b)*(1.0d0/(1.0d0+term))**(1.0d0+c)
      fmax=a+y0
      return
      end