c	simulation of the feigenbaum fixed point with small noise
c	depends on file gfix.dat created by newtonren 

	integer n, ndata
	parameter (n=128, ndata=5000) 
	real *8 gfix(n)
	real *8 cmag, rmag
	real*8 sigma,x(ndata),xm,xs,ave,d,dkstat,sdev,epsi,prob,var,uu, aux
	real *8 lgstep
	cmag= 1.d0
	rmag= 2.5d0

	xs=0.d0
	nitir=1024
c  The coefficients of the fixed point of the RG are read
	open(1,file='gfix.dat')
	do k=1,n
	   read(1,'(d25.15)') gfix(k)
	enddo

c initializes the seed for our random generator
	write(*,*)'give seed'
	read(*,*)idum

	open(2,file='noise.dat')
	sigma=.1d-9
c	open(3,file='answer.dat') 
	do ksigma=1,4
	   do j=1,ndata
	      xm=0.d0
	      do m=1,nitir	
		 u=ran1(idum)
		 uu=dble(u)
		 uu=uu-.5d0
		 aux= (xm*xm-cmag)/rmag
		 call evalu(n,gfix,aux,xm)
		 xm = xm + uu*sigma
	      enddo
	      x(j)=xm
	   enddo

c	begining of statistical kolmogorov goodness test
	   call moment(x,ndata,ave,var,sdev)
	   write(*,'(1x,a,d25.15,1x,a,d25.15)') 'mean=',ave, 'deviation=',sdev
	   write(*,'(1x,a,d25.15)') 'sigma/sdev=',sdev/sigma
	   do j=1,ndata
	      x(j)=(x(j)-ave)/sdev
c	      write(3,'(d25.15)') x(j)
	      enddo
	   call sort(ndata,x)
	   call ksone(x,ndata,ave,sdev,d,dkstat,prob)
	   write(*,'(a,d25.15,1x,a,d25.15,/)') 'ksstat=',dkstat,'ls=',prob
	lgstep=dlog(dble(nitir))
	   write(2,'(d25.15,1x,d25.15)')lgstep, dlog(1/sigma) + dlog(sdev)
	   sigma=sigma*.25d0
	   nitir=nitir*2
	enddo
 	stop
       	end		

	double precision function normal(y,mean,sd)
	real*8 y,yp,mean,sd
	yp=(y-mean)/(dsqrt(2.d0)*sd)
	normal= .5d0*derfc(-yp)
	return
	end


    	function ran1(idum)
	integer idum,ia,im,iq,ir,ntab,ndiv
	real*8 ran1,am,eps,rnmx
	parameter (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836,
     *		ntab=32,ndiv=1+(im-1)/ntab,eps=1.2e-7,rnmx=1.-eps)
	
	integer j,k,iv(ntab),iy
	save iv,iy
	data iv/ntab*0/, iy/0/
	if (idum.le.0.or.iy.eq.0) then
	    idum=max(-idum,1)
	    do  j=ntab+8,1,-1
	    	k=idum/iq
	    	idum=ia*(idum-k*iq)-ir*k
	    	if (idum.lt.0) idum =idum + im
	    	if (j.le.ntab) iv(j)=idum
	    enddo 
	    iy=iv(1)
	endif
	k=idum/iq
	idum=ia*(idum-k*iq)-ir*k
	if (idum.lt.0) idum=idum+im
	j=1+iy/ndiv
	iy=iv(j)
	iv(j)=idum
	ran1=min(am*iy,rnmx)
	return
	end 		

      subroutine ksone(data,n,mu,sig,d,dkstat,prob)
      integer n
      real*8 d,data(n),prob, dkstat,mu,sig
      integer j
      real*8 dt,en,ff,fn,fo,normal,probks
      en=n
      d=0.d0
      fo=0.d0
      do j=1,n
         fn=j/en
         ff=normal(data(j),0.d0,1.d0)
         dt=max(abs(fo-ff),abs(fn-ff))
         if (dt.gt.d) d=dt
         fo = fn
      enddo
      en=sqrt(en)
      dkstat=en*d
      prob=probks((en+0.12d0 + 0.11d0/en)*d)
      return
      end

      function probks(alam)
      real*8 probks, alam, eps1, eps2
      parameter (eps1=0.001d0, eps2=1.d-8)
c     kolmogorov-smirnov function
      integer j
      real*8 a2, fac, term, termbf
      a2=-2.d0*alam**2
      fac=2.d0
      probks=0.d0
c     below previous term in sum
      termbf=0.d0
      do j=1,100
         term=fac*exp(a2*j**2)
         probks=probks+term
         if(dabs(term).le.eps1*termbf.or.dabs(term).le.eps2*probks) then            
      return
      else
      fac=-fac
         termbf=dabs(term)
         endif
      enddo
c      get here only by failing to converge
      probks=1.d0  
      return   
      end
 
      subroutine sort(n,data)
      integer n,m,nstack
      real*8 data(n)
      parameter (m=7,nstack=50)
c     sorts array(1:n) into ascending numerical order using
c     the quicksort algorithm. n is input; arr is replaced
c     on output by its sorted arrangement.
      integer i,ir,j,jstack,k,l,istack(nstack)
      real*8 a,temp
      jstack=0
      l=1
      ir=n
 1    if (ir-l.lt.m) then
         do j=l+1,ir
            a=data(j)
            do i=j-1,1,-1
               if(data(i).le.a) goto 2
               data(i+1)=data(i)
            enddo
            i=0
 2          data(i+1)=a
         enddo
         if(jstack.eq.0) return
         ir=istack(jstack)
         l=istack(jstack-1)
         jstack=jstack-2
      else
         k=(l+ir)/2
         temp=data(k)
         data(k)=data(l+1)
         data(l+1)=temp
         if(data(l+1).gt.data(ir)) then
            temp=data(l+1)
            data(l+1)=data(ir)
            data(ir)=temp
         endif
         if(data(l).gt.data(ir)) then
            temp=data(l)
            data(l)=data(ir)
            data(ir)=temp
         endif
         if(data(l+1).gt.data(l)) then
            temp=data(l+1)
            data(l+1)=data(l)
            data(l)=temp
         endif
         i=l+1
         j=ir
         a=data(l)
 3       continue
         i=i+1
         if(data(i).lt.a) goto 3
 4       continue
         j=j-1
         if(data(j).gt.a) goto 4
         if(j.lt.i)goto 5
         temp=data(i)
         data(i)=data(j)
         data(j)=temp
         goto 3
 5       data(l)=data(j)
         data(j)=a
         jstack=jstack+2
         if(jstack.gt.nstack) pause 'nstack too small in sort'
         
         if(ir-i+1.ge.j-1) then
            istack(jstack)=ir
            istack(jstack-1)=i
            ir=j-1
         else
            istack(jstack)=j-1
            istack(jstack-1)=l
            l=i
         endif
      endif
      goto 1
      end

c	computes mean and variance (corrected two pass)
c	algorithm.
		
	subroutine moment(data,n,ave,var,sdev)
	real*8 ave,s, ep,p, sdev,var,data(n)
	s=0.d0
	do j=1,n
	s=s+data(j)
	enddo
	ave=s/n
	sdev=0.d0
	var=0.d0
	ep=0.d0
	do j=1,n
	s=data(j)-ave
	ep=ep+s
	p=s*s
	var=var+p
	enddo
	var=(var-ep**2/n)/(n-1)
	sdev=dsqrt(var)
	return
	end