The code of the IFS program:

      PROGRAM IFS

      IMPLICIT NONE

      integer ni,niter
      real pi
      parameter (pi=3.14159265359)
      integer i,j,i0,nf
      real xi(20,2),r(20,2),rho(20),theta(20),prob(0:20)
      real xn(2),xnp1(2),p,x,y(50000)
 

c     Lecture des parametres dans le fichier user.dat
      nf=1
      OPEN(nf,FILE='user.dat')

      read(nf,*) niter
      write(*,*) 'niter= ',niter

      read(nf,*) ni
      write(*,*) 'ni= ',ni

      do i=1,ni
         read(nf,*) (xi(i,j), j=1,2)
         write(*,*) 'x',i,'= ',xi(i,1),'   ','y',i,'= ',xi(i,2)
      enddo

      read(nf,*) (prob(i), i=1,ni)

c     prob represente en fait la somme cumulee des probabilites
      prob(0)=0
      do i=1,ni
         prob(i)=prob(i-1)+prob(i)
         write(*,*) 'p',i,'= ',prob(i)
      enddo

      do i=1,ni
         read(nf,*) rho(i), theta(i)
         write(*,*) 'rho',i,'= ',rho(i),'   ','theta',i,'= ',theta(i)
      enddo

      CLOSE(nf)

c     on passe en complexes cartesiens
      do i=1,ni
         r(i,1)=rho(i)*cos(theta(i)*pi/180)
         r(i,2)=rho(i)*sin(theta(i)*pi/180)
      enddo

      xn(1)=0.0
      xn(2)=0.0

      open(nf,FILE='attract.txt')

c     boucle principale
      do j=1,niter
c     sauvegarde des resultats
         write(nf,*) xn(1), xn(2)

c     on tire un nombre au hasard
         p=rand()

c    on cherche quelle fonction appliquer
         do i=1,ni
            if ((p.gt.prob(i-1)).and.(p.lt.prob(i))) then
               i0=i
            endif
         enddo

c     on fait le calcul du point suivant
         xnp1(1)=xi(i0,1)+r(i0,1)*(xn(1)-xi(i0,1))-
     &           r(i0,2)*(xn(2)-xi(i0,2))

         xnp1(2)=xi(i0,2)+r(i0,2)*(xn(1)-xi(i0,1))+
     &           r(i0,1)*(xn(2)-xi(i0,2))

c     on reactualise
         xn(1)=xnp1(1)
         xn(2)=xnp1(2)
      enddo

      close(nf)

      END
 

The code of the pdf calculator:

      PROGRAM proba
      integer niter,nz,i,j,nint
      real x,sup,inf,sup1,inf1,r
      real y(50000),z(50000)
      integer np(500000)
      character*20 filename
 

c     lecture des parametres
      open(1,FILE='proba.dat')
      read(1,*) filename
      read(1,*) niter
      read(1,*) inf,sup
      read(1,*) nint
      close(1)

c     lecture des donnes
      open(1,FILE=filename)
      do j=1,niter
         read(1,*) y(j)
      enddo
      close(1)

c     tri des points a traiter et calcul des bornes reels du domaine
      nz=0
      inf1=sup
      sup1=inf
      do j=1,niter
         if ((y(j).le.sup).and.(y(j).ge.inf)) then
            if (y(j).gt.sup1) then
               sup1=y(j)
            elseif (y(j).lt.inf1) then
               inf1=y(j)
            endif
            nz=nz+1
            z(nz)=y(j)
         endif
      enddo

      write(*,*) 'Nombre de points:', nz

      r=nint/(sup1-inf1)

      write(*,*) 'Nombre d''intervalles:', nint
      write(*,*) 'Borne Inf:', inf1
      write(*,*) 'Borne Sup:', sup1
      write(*,*) 'Taille d''echantillonage:',1.0/r

      do j=1,nz
         np(j)=0
      enddo

      do j=1,nz
         i=1+(z(j)-inf1)*r
         np(i)=np(i)+1
      enddo

      open(3,FILE='proba_'//filename)
      do i=1,nint
         x=inf1+(i-1)/r+0.5e0/r
         write(3,*) x, 1.0e0/nz*np(i)
      enddo
      close(3)
 

      END