Parameter niter=1000 gives the maximum number of iterations made before the final check for periodicity of solution. Diverging solutions are caught earlier than that. I've tried niter=10000 with only several values of k in the output different. Possible optimization would include catching the common k=3 outcomes earlier than at iteration niter. As is, this program runs usually for a few hours on a workstation.
c c The Complex Universe Survey c implicit double precision (a-h,o-z) parameter (niter=1000) c------------current parameter (NX=620, NY=310) parameter (x0=-4.2d0, x1=-4.0d0, y0=-.08d0, y1=0d0) c------------previous * parameter (NX=1000, NY=500) * parameter (x0=-.85d0, x1=.15d0, y0=-.5d0, y1=0d0) c------------previous * parameter (NX=800, NY=400) * parameter (x0=-5d0, x1=15d0, y0=-10d0, y1=0d0) c------------------ complex*16 z,z0 open(1,file='surv620.dat',status='new') write(1,*)NX,NY write(1,*)x0,x1 write(1,*)y0,y1 do 100 ix=1,NX do 100 iy=1,NY x=x0+(x1-x0)*ix/dble(NX) y=y0+(y1-y0)*iy/dble(NY) z=dcmplx(x,y) z0=z call object(z,z0,niter,k) write(1,2)k 2 format(i4) 100 continue close(1) stop end c c c subroutine object(z,z0,niter,k) parameter(nc=41,eps=1d-4,big=1d33) complex*16 z,z0,zz,zn zz=z0 do 10 j=1,niter if(j.eq.50.or.j.gt.niter-nc)then c---------check for NaN in the 50th step & the last nc iterations zre=real(zz) zim=imag(zz) c-----------k=0 means NaN result, k=-10 finite result k=0 if(zre.gt.-big.and.zre.lt.big.and.zim.gt.-big.and.zim.lt.big)k=-10 if(k.eq.0)return endif c---------iteration 10 zz=z**zz c---------the last 3*nc points are analysed for cyclic pattern do 20 jj=1,3 zn=zz do 20 j=1,nc zz=z**zz if(abs(zz-zn).lt.eps)then k=j return endif 20 continue c---------k=-10 no cycle found, k>0 period of the cycle return endBack to the home page of Pawel Artymowicz