c     this program calculated  the analytic hyperbilic fixed point 
c     of the RG operator using a Newton-type method. 
c     the external subroutines series.f, renor.f and drenor.f are needed
c     the critical logistic map is taken as initial point.
c     the method uses coordinates of the form g(x)=G((x^2-c)/r)

      integer n 
      parameter (n =128)
      double precision cmag,rmag
      common/param/cmag,rmag
      double precision magic
      parameter (magic =  0.14011551922975566714d1)
      double precision g(n)
      double precision grenorm(n)
      double precision drenop(n,n),gfix(n),gfixsc(n), sol(n), dgfix(n)
      double precision aux, aux1,aux2, aux3, auxtemp, workfix(n)
      double precision errdim(n), error, a, b
      integer i,j, iter
      open(1,file='errornum.dat')	
      cmag = .1d1
      rmag = .25d1

      call matzer(n,g)

      g(1) = -( magic*cmag - 1.d0)
      g(2) = - magic*rmag
      call matmov(n,g,gfix)
      iter=8
c   we perform a newton method to find fixed point of R(f)=f
c   with initial value gfix = g iter times

      do 10 i=1,iter
      call renor( n,gfix,cmag, rmag, grenorm)  
      call drenoprtor(n,gfix,cmag,rmag, drenop)
      call renfixed(n,drenop,gfix,grenorm,sol)
      call matsub(n,gfix,sol,workfix)
      call normerr(n,workfix,error)
      errdim(i)=error
      call matmov(n,sol,gfix)
 10   continue
   
      do 20 i=1,iter-1
      aux1=errdim(i)
      aux1=dabs(aux1)
      aux2=errdim(i+1)
      aux2=dabs(aux2)
      write(1,*)dlog(aux1),dlog(aux2)
 20   continue

c     here we estimate the norm |R(g)-g| where g=gfix is out 
c     estimate for the fixed point

      call renor(n,gfix,cmag, rmag, grenorm)
      call matsub(n,gfix,grenorm,workfix)
      call normerr(n,workfix,error)
      write(*,*)'residual=|R(g)-g|=',error
 

c     we scale back to the basis G(X)=h(r*X +c), where h(x^2)=g(x)
      call rnscale(n,gfix,gfixsc,cmag,rmag)
      call coefprint(n,gfixsc,gfix)
c we compute a few numbers that we need to compute
c among them are g(1)=lambda, g(lambda), g''(0), g'(lambda)
c the constants A anb B
      call deriv(n,gfix,dgfix)
	aux=(1.d0-cmag)/rmag
	call evalu(n,gfix,aux,aux1)
	aux=(aux1*aux1-cmag)/rmag
	call evalu(n,gfix,aux,aux2)
	write(*,*)'lambda= ',aux1
	write(*,*)'g(lambda)= ' , aux2
	aux=-cmag/rmag
	call evalu(n,dgfix,aux,auxtemp)
        aux3=2*auxtemp/rmag
	write(*,*) 'd2g(0)=', aux3
        aux=(aux1*aux1-cmag)/rmag
        call evalu(n,dgfix,aux,auxtemp)
        aux=2.d0*auxtemp*aux1/rmag
        write(*,*) 'Dg(lambda)=', aux
        a=aux*aux*aux2*aux2 / (aux1*aux1*aux3)
	b=aux1*aux3*aux3/ (aux*aux2)
	write(*,*)'A=', a
	write(*,*)'B=', b
	call deriv(n,gfix,dgfix)
      call kurto(n,gfix,dgfix,cmag,rmag)
      stop
      end


c     this subroutine calculates the derivative of the RG operator
c     at the point f in the standard basis

      subroutine drenoprtor(n,f,c, r,drenopf)
      integer n
      integer i
      double precision f(n), drenopf(n,n), c, r
      double precision h(n), drenh(n)
c     we write the derivative of renorm operator in standar basis. we
c     call the derivative at g drenop.
      do 20 i=1,n
         call matzer(n,h)
         h(i)=1.d0
         call drenor(n, f, h,c,r, drenh)
         call matrmov(n,drenh,i,drenopf)
 20   continue 
      return
      end



c     this subroutine performs a newton-type algorithm to compute 
c     the fixed point of the renormalization operator

      subroutine renfixed(n,oprtor,f, rf,sol)
      integer n
      integer i, info, lwork,rank
      double precision rcond
      double precision oprtor(n,n),f(n),rf(n),hf(n),sol(n)
      double precision s(n), work(5*n), hfsvd(n,1)
      call matsub(n, f, rf, hf) 

c  in this part we solve (DR(f)-Id)*h = R(f)-f

      do 10 i=1,n
         oprtor(i,i)=oprtor(i,i)-1.d0
 10   continue
      rcond=-1.d0
      do 20 i=1,n
         hfsvd(i,1)=hf(i)
 20   continue

c     this calls lapack subroutine dgless.
c     SUBROUTINE DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
c    $                   WORK, LWORK, INFO )
c     the solution to the problem inf| Ax-b|_2 is stored in b in the output
      
      call dgelss(n,n,1,oprtor,n,hfsvd,n,s,rcond,rank,
     $            work, 5*n, info)
      do 30 i=1,n
       hf(i)=hfsvd(i,1)
 30    continue
      call matadd(n,f,hf,sol)
      return
      end
    
      subroutine coefprint(n,f,ff)
      integer n
      double precision f(n),ff(n)
      integer i
      open(2,file='coef.dat')
      open(3,file='gfix.dat')	
      do 10 i=1,n
      write(2,'(d25.15)') f(i)
      write(3,'(d25.15)') ff(i)	
 10   continue
      return
      end

	subroutine kurto(n,f,df,c,r)
      integer n, nfin
      parameter ( nfin=2**14)
      real *8 f(n), df(n)
      real *8 c, r
      real *8 aux, aux2, daux, x, dx, dx2, dx4, delta4, delta2, kurt
      real *8 dlogkurt
      integer i
      open(1,file='kurt.dat')
      x=1.0d0
      delta2=1.d0
      delta4=1.d0
      kurt=1.d0
      i=1
      write(1,*)dlog(dble(i))/dlog(2.d0),dlog(kurt)/dlog(2.d0)

c	Here we do the iterations to find derivatives
	
      do i=2,nfin
         aux=(x*x - c)/r
      call evalu(n,f,aux,x)
         aux2=(x*x - c)/r
      call evalu(n,df,aux2,daux)
	dx=2.d0*x*daux/r
	dx2=dx*dx
	dx4=dx2*dx2
	delta2= 1.d0+dx2*delta2
	delta4= 1.d0+ dx4*delta4
	dlogkurt=dlog(delta4)-2.d0*dlog(delta2)
	write(1,*)dlog(dble(i))/dlog(2.d0), dlogkurt/dlog(2.d0)
	enddo
	return
	end