c     this program computes the m-th order cummulant of feigenbaum's
c     fixed point g. here g(x)=G((x^2-c)/r) and we use the standard basis 
c     in the coordinate system h(x)=H((x^2-c)/r)
c     this programs depends on series.f,  kum.f and dfgamma.f

      integer n, lwork, pw
      parameter (n =128, lwork=2*n, pw=12)
      double precision cmag,rmag
      parameter (cmag=.1d1, rmag=.25d1)
      double precision alpha, beta, gamma
      double precision g(n), ggam(n), dgofgal(n), dggalpw(n)
      double precision h(n),kummgh(n)
      double precision kummmat(n,n), matop(n,n)
      double precision reeig(n),imeig(n)
      integer  ilo, ihi, iwork(2*n-2)
      double precision scale,  abnrm, eignrm, eigmax
      double precision rconde(n),rcondv(n), vl(n,n), vr(n,n)
      double precision work(lwork)
      integer info
      integer i,j, mcount

      open(1,file='gfix.dat')
      open(2,file='absradspec.dat')

      do 10 i=1,n
         read(1,'(d25.15)')g(i)
 10   continue
      close(1)

c     we compute  -g'(Gamma(x))*g(gamma*X+beta)*2/r --> dgatgalp
c     Gamma(X)=(g^2(gamma *X + beta)-c)/r --> ggam
c     gamma=alpha^2 beta=(c/r)*(gamma - 1)
      call dfgamma(n,g,cmag,rmag,alpha,beta,gamma,ggam,dgofgal)
c     we compute  [-g'(Gamma(x))*g(gamma*X+beta)*2/r]^m --> dggalpw
	write(*,'(a,6x,a)') 'order', 'spec_rad'
	write(*,*)'---------------------------------'
      do mcount=1, pw
      call fm(n,mcount,dgofgal,dggalpw)
c     we use the standard basis h_k(x)=((x^2-c)/r)^k to get a representation
c     of the operator Kg_m
      do 20 i=1,n
         call matzer(n,h)
         h(i)=1.d0
         call kumabs(n, mcount, h,cmag,rmag, alpha, beta, gamma,
     *        dggalpw, ggam, kummgh)
         call matrmov(n,kummgh,i,kummmat)
 20   continue
      do 30 i=1,n
         do 40 j=1,n
            matop(i,j)=kummmat(i,j)
 40      continue
 30   continue


c     SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
c     $                VL,LDVL,VR,LDVR,ILO,IHI,SCALE,ABNRM,
c     $                RCONDE,RCONDV,WORK, LWORK, IWORK, INFO)
c     this LAPACK routine computes the eigenvalues of operator matop

      call dgeevx('b','n','n','n',n,matop, n,reeig,imeig,
     $        vl,n,vr,n,ilo,ihi,scale,abnrm,
     $        rconde,rcondv,work,lwork,iwork,info)


c      do 45 i = 1, n
c         write(2,'(d25.15, d25.15)')  reeig(i), imeig(i)
c 45   continue
      eigmax=reeig(1)*reeig(1)+ imeig(1)*imeig(1)
      eigmax=dsqrt(eigmax)
      do 50 i=2,n
         eignrm=reeig(i)*reeig(i) + imeig(i)*imeig(i)
         eignrm=dsqrt(eignrm)
         eigmax=dmax1(eignrm,eigmax)
 50   continue
	write(2,'(i4,2x,d25.15)') mcount, eigmax
      write(*,'(i4,3x,d25.15)')mcount, eigmax
      enddo

      close(2)
 60   stop 
      end