c     this routine finds g'(g(alpha x))) in terms of g(x)=f(X)
c     where  X=(x^2-c)/r 

      subroutine dfgamma(n,f,c,r,alphaf,betaf,gammaf,workf,dfatfalp)
      integer n,p
      double precision f(n),workf(n), fatalp(n), dfatfalp(n)
      double precision c,r
      double precision df(n),newftemp(n), newf(n)
      double precision aux, auxtemp
      double precision alphaf, betaf, gammaf

c     the terms alpha(f), gamma(f)=alpha^2(f), beta(f)=(c/r)*(gamma(f)-1),
c     are calculated
      
      aux= -c/r
      call evalu(n,f,aux,auxtemp)
      aux = (auxtemp*auxtemp - c)/r
      call evalu(n,f,aux,alphaf)
      gammaf=alphaf*alphaf
      betaf=c/r
      betaf = betaf*(gammaf -1.d0)
 
c -------------------------------------------------------------

c     we compute the funtion Gammaf(X)= (f^2(gammaf*X + betaf)-c)/r
c     and store it in workf.
c     function f(gammaf*X + betaf) is stored in fatalp

      call mataffine(n,f,fatalp,betaf,gammaf)
      call matmov(n,fatalp,newftemp)
      call prod(n,fatalp,newftemp,workf)
      workf(1) = workf(1) - c
      call matsca(n,workf,1.0d0/r)
c -------------------------------------------------------------
c we compute f'(Gammaf(X)) and store it in newf

      call deriv(n,f,df)
      call compos(n,workf,df,newf)

c     we compute -f'(Gammaf(X))*f(gamma*X+beta)*2/r
c     result is stored in dfatfalp

      call prod(n,newf,fatalp,dfatfalp)
      call matsca(n,dfatfalp,-2.d0/r)
      return
      end