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
![]() |
![]() |
![]() |