c     given function f and h small we compute DR(f)h
c     f is given in X=(x^2-c)/r coordinates

      subroutine drenor( n, f, hh,c,r,drn)
      integer n
      double precision f(n), hh(n),drn(n)
      double precision c,r
      double precision workf(n), newf(n), workh(n), df(n),alphaf, betaf
      double precision workdf(n),workdrn(n), newftemp(n),aux, aux2
      double precision auxh, auxtemp,dalphafh, lambdaf, gammaf, dgammafh

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

      call evalu(n,hh,aux2,aux)
      auxh=(2.d0/r)*auxtemp*aux
      call deriv(n,f,df)
      call evalu(n,df,lambdaf,aux)
      dalphafh= aux*auxh 
      call evalu(n,hh,lambdaf,aux)
      dalphafh = dalphafh + aux
      dgammafh=2.d0*alphaf*dalphafh 
c -------------------------------------------------------------

c     we compute the funtion Gamma(f)(t)= (f^2(gammaf*t + betaf)-c)/r
c     which will be stored in workf

      call mataffine(n,f,newf,betaf,gammaf)
      call matmov(n,newf,newftemp)
      call prod(n,newf,newftemp,workf)
      workf(1) = workf(1) - c
      call matsca(n,workf,1.0d0/r)
c -------------------------------------------------------------

      call matzer(n,newf)
      call compos(n,workf,f, newf)
      call matsca(n,newf,dalphafh/gammaf)
      call compos(n,workf,hh,workh)
      call matsca(n,workh,1.d0/alphaf)
      call matsub(n,workh,newf,workdrn)
 
c --------------------------------------------------------------
    
c     computes f'(gammaf*t + betaf)*dgammafh*(t+c/r) + h(gammag*t + betaf)
c     and stores that in newf

      call mataffine(n,df,workdf,betaf,gammaf)
      call prodaffine(n,workdf,workdf,c/r,1.d0)
      call matsca(n,workdf,dgammafh)

      call matzer(n,workh)
      call matzer(n,newf)

      call mataffine(n,hh,workh,betaf,gammaf)
      call matadd(n,workdf,workh,newf)

c ---------------------------------------------------------------
      call matzer(n,workdf)
      call matzer(n,workh)
      call matzer(n,newftemp)
c     
      call mataffine(n,f,workdf,betaf,gammaf)
      call compos(n,workf,df,workh)

      call prod(n,workh,workdf,newftemp)

      call matzer(n,workdf)

      call prod(n,newftemp,newf,workdf)
      aux=r*alphaf
      aux=2.d0/aux
      call matsca(n,workdf,aux)
      call matadd(n, workdf,workdrn,drn)
      return
      end