subroutine function_generic(ndim,beta,nxvars,
     &  xvars,nyvars,ymod)
      integer ndim,nxvars,nyvars
      double precision beta(ndim),xvars(nxvars),
     &  ymod(nyvars)
      double precision y0,a,b,c,x0,x,term,crit
      parameter(crit=300.0d0)
      a=beta(1)
      b=beta(2)
      c=beta(3)
      x0=beta(4)
      y0=beta(5)
      x=xvars(1)
      if((-(x-x0)/b).lt.crit)then
        term=dexp(-(x-x0)/b)
        ymod(1)=y0+a*(1.0d0/(1.0d0+term))**c
      else
        term=dexp((x-x0)/b)
        ymod(1)=y0+a*(term/(1.0d0+term))**c
      endif
      return
      end
      
      subroutine der_function_generic(np,beta,m,
     &    xvars,nq,der_beta)
      implicit none
      integer np,m,nq
      double precision beta(np),xvars(m),der_beta(np,nq)
      double precision y0,a,b,c,x0,x,term,crit
      parameter(crit=300.0d0)
      a=beta(1)
      b=beta(2)
      c=beta(3)
      x0=beta(4)
      y0=beta(5)
      x=xvars(1)
      der_beta(5,1)=1.0d0
      if((-(x-x0)/b).lt.crit)then
        term=dexp(-(x-x0)/b)
        der_beta(1,1)=(1.0d0/(1.0d0+term))**c
!        der_x=(a*c*term/b)*(1.0d0/(1.0d0+term))**(1.0d0+c)
        der_beta(4,1)=-(a*c*term/b)*
     &    (1.0d0/(1.0d0+term))**(1.0d0+c)
        der_beta(2,1)=-(a*c*term*(x-x0)/(b*b))*
     &    (1.0d0/(1.0d0+term))**(1.0d0+c)
        der_beta(3,1)=-(a*dlog(1.0d0+term))*
     &    (1.0d0/(1.0d0+term))**c
      else
        term=(x-x0)/b
        der_beta(1,1)=(dexp(term)/(1.0d0+dexp(term)))**c
!        der_x=(a*c/b)*(dexp(term*c/(c+1.0d0))/
!     &              (1.0d0+dexp(term)))**(c+1.0d0)
        der_beta(4,1)=-(a*c/b)*(dexp(term*c/(c+1))/
     &              (1.0d0+dexp(term)))**(c+1.0d0)
        der_beta(2,1)=-(a*c*(x-x0)/(b*b))*(dexp(term*c/
     &    (c+1.0d0))/(1.0d0+dexp(term)))**(1.0d0+c)
        der_beta(3,1)=-a*(dlog(1.0d0+dexp(term))-term)*
     &           (dexp(term)/(1.0d0+dexp(term)))**c
      endif
      return
      end

      subroutine indices_function_generic(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