Index: TrouNoir/trou_noir.f
===================================================================
--- TrouNoir/trou_noir.f	(revision 204)
+++ TrouNoir/trou_noir.f	(revision 205)
@@ -26,7 +26,7 @@
       REAL*8 rs,ri,re           ! rayons de Schwarzschild,interieur et exterieur
       REAL*8 tho                ! angle normale-disque & direction-observateur
       REAL*8 ro                 ! distance entre l'observateur et le trou noir
-      INTEGER*4 q               ! parametre de spectre
+      REAL*8 q               ! parametre de spectre
       REAL*4 bss                ! borne superieure du spectre
 
 c
@@ -34,7 +34,7 @@
 c
 
       INTEGER*4 dim             ! dimension de la plaque photographique (en pixels)
-      PARAMETER (dim=512)
+      PARAMETER (dim=1024)
       REAL*4 stp                ! pas d'iteration en pixels sur la plaque
       INTEGER*4 nmx             ! nombre total d'iterations
       REAL*8 d                  ! distance exploree du centre de la plaque
@@ -59,8 +59,8 @@
 c RESULTATS DE L'INTEGRATION DE L'EQUATION DIFFERENTIELLE
 c
 
-      REAL*8 rp(2000)           ! vecteur polaire
-      REAL*8 nh                 ! nombre d'avancements angulaires
+      REAL*8 rp(2048)    ! vecteur polaire
+      INTEGER*4 nh              ! nombre d'avancements angulaires
 
 c
 c RESULTATS DE LA SIMULATION
@@ -68,8 +68,6 @@
 
       REAL*4 zp(dim,dim) 	! image en decalage spectral
       REAL*4 fp(dim,dim) 	! image en flux
-      REAL*8 fx(200)            ! flux par tranche de decalage de frequence
-      INTEGER*4 nt(200) 	! nombre de photon par tranche de decalage
 
 c
 c AUTRES VARIABLES...
@@ -92,29 +90,25 @@
       LOGICAL raie            ! variable interactive de type de spectre
       LOGICAL dist            ! variable interactive de distance
 
+      REAL*4 t_start,t_stop
 c
 c TYPE DE SPECTRE
 c
 
-      raie=.TRUE.               ! TRUE -> raie infiniment fine
-c                               ! FALSE -> corps noir
+c      raie=.TRUE.               ! TRUE -> raie infiniment fine
+      raie=.FALSE.               ! FALSE -> corps noir
 
 c
 c TYPE DE DISTANCE
 c
 
-      dist=.TRUE.               ! TRUE -> distance infinie
-c                               ! FALSE -> distance finie
+c      dist=.TRUE.               ! TRUE -> distance infinie
+      dist=.FALSE.               ! FALSE -> distance finie
 
 c
 c INITIALISATION DES TABLEAUX
 c
 
-      DO i=1,200
-         fx(i)=0.
-         nt(i)=0
-      END DO
-      
       DO i=1,dim
          DO j=1,dim
             zp(i,j)=0.
@@ -148,7 +142,7 @@
       b=0.
       thx=asin(bmx/ro)
       pc=0
-      nbr=200
+      nbr=256
       IF (raie .EQV. .TRUE.) THEN
          bss=2.
       ELSE
@@ -159,15 +153,15 @@
 c BOUCLE PRINCIPALE DE LA SIMULATION
 c
 
+      CALL TIMER(t_start)
+      
       DO n=1,nmx
-         pc=INT(100./nmx*n)
-         IF (pc.NE.INT(100./nmx*(n-1))) WRITE (*,*) pc, ' %'
 
 c
 c Affectation des conditions initiales de la resolution Runge-Kutta
 c ( Premier etage d'imbrication -> n )
 c
-         h=pi/500.
+         h=4*pi/2048.
          d=stp*n
          IF (dist .EQV. .TRUE.) THEN
             db=bmx/nmx
@@ -206,7 +200,7 @@
 c ( Deuxieme etage d'imbrication -> n & 1 )
 c
 
-         DO i=nh+1,2000
+         DO i=nh+1,2048
             rp(i)=0. 
          END DO
 
@@ -261,8 +255,8 @@
 c S'il y a impact calcul - du rapport de frequence
 c - du spectre de photons et de flux
 
-                  CALL impact(raie,d,phi,dim,r,b,tho,m,zp,fp,nt,
-     & fx,q,db,h,nbr,bss)
+                  CALL impact(raie,d,phi,dim,r,b,tho,m,zp,fp,
+     & q,db,h,nbr,bss)
 
                END IF
 
@@ -271,16 +265,21 @@
          END DO
       END DO
 
+
+      CALL TIMER(t_stop)
+
+      WRITE (*,*) 'Elapsed time = ', t_stop - t_start
+      
 c Appel de la routine d'affichage
 
-      CALL affichage(zp,fp,nt,fx,dim,nbr,bss,bmx,m,ri,re,tho,pi)
+      CALL affichage(zp,fp,dim,nbr,bss,bmx,m,ri,re,tho,pi)
 
       WRITE(*,*)
       WRITE(*,*) 'C''est TERMINE ! '
       WRITE(*,*) 
 
       STOP
-      END
+      END PROGRAM
 
 c
 c Cette fonction permet de trouver un arctangente entre 0 et 2*pi
@@ -311,6 +310,19 @@
       END
 
 c
+c
+c
+      SUBROUTINE timer(ytime)
+      REAL*4 ytime
+      INTEGER*4 clk_max,clk_rate,clk_read
+
+      CALL SYSTEM_CLOCK(clk_read,clk_rate,clk_max)
+
+      ytime=FLOAT(clk_read)/FLOAT(clk_rate)
+
+      RETURN
+      END
+c
 c Premiere fonction du systeme -> ( d(u)/d(phi)=f(phi,u,v) )
 c
 
@@ -366,9 +378,9 @@
 c
 
       SUBROUTINE impact(raie,d,phi,dim,r,b,tho,m,zp,fp, 
-     & nt,fx,q,db,h,nbr,bss)
-      REAL*8 d,phi,r,b,tho,m,fx(200),db,h
-      INTEGER*4 dim,nt(200),q,nbr
+     & q,db,h,nbr,bss)
+      REAL*8 d,phi,r,b,tho,m,db,h,q
+      INTEGER*4 dim
       REAL*4 zp(dim,dim),fp(dim,dim),bss
       REAL*4 xe,ye
       INTEGER*4 xi,yi
@@ -390,9 +402,9 @@
 c pour le cas d'une raie infiniment fine ou un corps noir
 
       IF (raie .EQV. .TRUE.) THEN
-         CALL spectre(nt,fx,rf,q,b,db,h,r,m,nbr,bss,flx)
+         CALL spectre(rf,q,b,db,h,r,m,bss,flx)
       ELSE
-         CALL spectre_bb(nt,fx,rf,q,b,db,h,r,m,nbr,bss,flx)
+         CALL spectre_bb(rf,b,db,h,r,m,nbr,bss,flx)
       END IF
 
 c Affectation sur chacune des images du decalage spectral et du flux
@@ -413,32 +425,25 @@
 
 c Calcul du flux puis du spectre pour une raie infiniment fine
 
-      SUBROUTINE spectre(nt,fx,rf,q,b,db,h,r,m,nbr,bss,flx)
-      REAL*8 fx(200)
-      REAL*8 rf,b,db,h,r,m,flx
+      SUBROUTINE spectre(rf,q,b,db,h,r,m,bss,flx)
+      REAL*8 rf,b,db,h,r,m,flx,q
       REAL*4 bss
-      INTEGER*4 nt(200),q,nbr
-      INTEGER*4 fi
-      fi=INT(rf*nbr/bss)
-      nt(fi)=nt(fi)+1
       flx=(r/m)**q*rf**4*b*db*h
-      fx(fi)=fx(fi)+flx
       RETURN
       END
 
 c Calcul du flux puis du spectre pour un corps noir
 
-      SUBROUTINE spectre_bb(nr,fx,rf,q,b,db,h,r,m,nbr,bss,flx)
-      REAL*8 fx(200),nu_rec,nu_em,qu,v
+      SUBROUTINE spectre_bb(rf,b,db,h,r,m,nbr,bss,flx)
+      REAL*8 nu_rec,nu_em,qu,v
       REAL*8 rf,b,db,h,r,m,flx,temp,temp_em,flux_int
       REAL*8 m_point
       REAL*4 bss
-      INTEGER*4 nr(200),q,nbr
+      INTEGER*4 nbr
       INTEGER*4 fi,posfreq
-      REAL*8 planck,c,k
-      PARAMETER (planck=6.62e-34) ! definition des constantes
-      PARAMETER (c=3.e8)
-      PARAMETER (k=1.38e-23)
+      REAL*8, PARAMETER :: planck=6.62e-34 ! definition des constantes
+      REAL*8, PARAMETER :: c=3.e8
+      REAL*8, PARAMETER :: k=1.38e-23
 
 c
 c Definition des parametres pour le calcul de la temperature d'emission
@@ -454,14 +459,14 @@
 
 c qu -> Q du rapport
 
-      qu=((1-3./r)**-.5)*(r**-.5)*((r**.5-6.**.5)+
-     & ((3.**.5)/2.)*log ((r**.5)-(3.**.5))*0.171517)
+      qu=((1-3./r)**-0.5)*(r**-0.5)*(r**0.5-6.**0.5+
+     & (3.**0.5/2.)*log((r**0.5+3.**0.5)/(r**0.5-3**0.5)*0.1715728752))
       qu=abs(qu)
-      temp_em=temp*m**.5
-      temp_em=temp_em*m_point**.25
-      temp_em=temp_em*r**(-.75)
-      temp_em=temp_em*v**(-1./8.)
-      temp_em=temp_em*qu**.25
+      temp_em=temp*m**0.5
+      temp_em=temp_em*m_point**0.25
+      temp_em=temp_em*r**(-0.75)
+      temp_em=temp_em*v**(-0.125)
+      temp_em=temp_em*qu**0.25
 
 c 
 c initialisation des compteurs de flux
@@ -479,7 +484,7 @@
       DO fi=1,nbr
          nu_em=bss*fi/nbr 	!frequence d' emission
          nu_rec=nu_em*rf 	!frequence de reception
-         posfreq=INT(nu_rec*200./bss) !case ou se trouve nu rec
+         posfreq=INT(nu_rec*FLOAT(nbr)/bss) !case ou se trouve nu rec
 
 c
 c test pour savoir si la frequence de reception est bien dans le
@@ -487,26 +492,22 @@
 c
 
          if ((posfreq.GE.1).AND.(posfreq.LE.nbr)) THEN
-            nr(posfreq)=nr(posfreq)+1
 
 c Loi du corps noir
 
-            flux_int=(2.*h*(nu_em**3/c**2))*
-     & (1./(exp(planck*nu_em/(k*temp_em))-1.))
+            flux_int=2.*nu_em**3/c**2*planck/
+     & (exp(planck*nu_em/(k*temp_em))-1.)
 
 c Integration sur le pixel
 
             flux_int=flux_int*rf**3*b*db*h
 
-c Remplissage du spectre
+c Integration bolometrique
 
-            fx(posfreq)=flux_int+fx(posfreq)
-
-c Intergration bolometrique
-
             flx=flx+flux_int
          ENDIF
       END DO
+
       RETURN
       END
 
@@ -521,20 +522,20 @@
 c - Dessin du spectre en nombre de photons
 c
 
-      SUBROUTINE affichage(zp,fp,nt,fx,dim,nbr,bss,bmx,m,ri,re,tho,pi)
-      INTEGER*4 nt(200),dim,nbr
+      SUBROUTINE affichage(zp,fp,dim,nbr,bss,bmx,m,ri,re,tho,pi)
+      INTEGER*4 dim
       REAL*4 zp(dim,dim),fp(dim,dim),bss
-      REAL*8 fx(200),bmx,m,ri,re,tho,pi
+      REAL*8 bmx,m,ri,re,tho,pi
       INTEGER*4 np,i,j
       REAL*8 fm,nm
-      REAL*4 fmx,zmx,x(200),y(200),tr(6),c(20),bmp
+      REAL*4 fmx,zmx,x(1024),y(1024),tr(6),c(20),bmp
       LOGICAL tg
       LOGICAL ps
 
       tg=.TRUE.                 ! TRUE : affichage en tons de gris
 c                               ! FALSE : affichage en contours
 
-      ps=.FALSE.                ! TRUE : affichage postscript 
+      ps=.TRUE.                ! TRUE : affichage postscript 
 c                               ! FALSE : affichage fenetre Xwindow
 
 
@@ -542,16 +543,6 @@
 c TRANSFORMATION DE L'IMAGE DES FLUX : FLUX -> FLUX RELATIFS
 c
 
-      np=0
-      fm=0.
-      DO i=1,dim
-         DO j=1,dim
-            IF (fp(i,j).NE.0.) np=np+1
-            fm=fm+fp(i,j)
-         END DO
-      END DO
-      fm=fm/np
-
 c
 c DETERMINATION DU DECALAGE SPECTRAL MAXIMUM ET DU FLUX MAXIMUM
 c
@@ -560,12 +551,15 @@
       zmx=0.
       DO i=1,dim
          DO j=1,dim
-            fp(i,j)=fp(i,j)/fm
             IF (fmx.LT.fp(i,j)) fmx=fp(i,j)
             IF (zmx.LT.zp(i,j)) zmx=zp(i,j)
          END DO
       END DO
 
+      write(*,*) 'Flux max :',fmx
+      write(*,*) 'Z max :',zmx
+      
+      
       DATA tr/0.,1.,0.,0.,0.,1./
 
 c
@@ -573,7 +567,7 @@
 c
 
       if (ps .EQV. .TRUE.) THEN
-         CALL PGBEGIN (0,'image1.ps/ps',1,1)
+         CALL PGBEGIN (0,'z.ps/ps',1,1)
       ELSE
          CALL PGBEGIN (0,'/xwindow',1,1)
       END IF
@@ -591,7 +585,7 @@
          DO i=1,20
             c(i)=0.2*(i-1)
          END DO
-         CALL PGCONS(zp,dim,dim,1,dim,1,dim,c,12,tr)
+         CALL PGCONS(zp,dim,dim,1,dim,1,dim,c,20,tr)
       END IF
       CALL PGEND
 
@@ -600,7 +594,7 @@
 c
 
       if (ps .EQV. .TRUE.) THEN
-         CALL PGBEGIN (0,'image2.ps/ps',1,1)
+         CALL PGBEGIN (0,'flux.ps/ps',1,1)
       ELSE
          CALL PGBEGIN (0,'/xwindow',1,1)
       END IF
@@ -622,114 +616,7 @@
       END IF
       CALL PGEND
 
-c
-c AFFICHAGE DES DONNEES PHYSIQUES
-c
-
-      bmp=bmx
-
-      if (ps .EQV. .TRUE.) THEN
-         CALL PGBEGIN (0,'image3.ps/ps',1,1)
-      ELSE
-         CALL PGBEGIN (0,'/xwindow',1,1)
-      END IF
-      CALL PGENV (-bmp,bmp,-bmp,bmp,1,0)
-      CALL PGLABEL('Rayon de Schwarzschild &
-     &Parametre d''impact critique','',
-     &'DISQUE D''ACCRETION AUTOUR D''UN TROU NOIR')
-      DO i=1,200
-         x(i)=2.*m*cos(2.*pi*i/200.)
-         y(i)=2.*m*sin(2.*pi*i/200.)
-      END DO
-      CALL PGLINE(200,x,y)
-      DO i=1,200
-         x(i)=m*SQRT(27.)*cos(2.*pi*i/200.)
-         y(i)=m*SQRT(27.)*sin(2.*pi*i/200.)
-      END DO 
-      CALL PGLINE(200,x,y)
-      CALL PGEND
-
-c
-c AFFICHAGE DES CONTOURS DU DISQUE
-c
-
-      if (ps .EQV. .TRUE.) THEN
-         CALL PGBEGIN (0,'image4.ps/ps',1,1)
-      ELSE
-         CALL PGBEGIN (0,'/xwindow',1,1)
-      END IF
-      CALL PGENV(-bmp,bmp,-bmp,bmp,1,0)
-      CALL PGLABEL('Limites interieure & 
-     &exterieure du disque a l''infini','',
-     &'DISQUE D''ACCRETION AUTOUR D''UN TROU NOIR')
-      DO i=1,200
-         x(i)=ri*cos(2.*pi*i/200.)
-         y(i)=ri*cos(tho)*sin (2.*pi*i/200.)
-      END DO
-      CALL PGLINE(200,x,y)
-      DO i=1,200
-         x(i)=re*cos(2.*pi*i/200.)
-         y(i)=re*cos(tho)*sin(2.*pi*i/200.)
-      END DO
-      CALL PGLINE(200,x,y)
-      CALL PGEND
-
-c
-c DETERMINATION DES FLUX ET NOMBRE D'IMPACTS MOYENS SUR LES SPECTRES
-c
-
-      fm=0.
-      nm=0.
-      stp=bss/nbr
-      j=0
-      do i=1,nbr
-         x(i)=i*stp
-         fm=fm+fx(i)
-         nm=nm+1.*nt(i)
-         IF (nt(i).GT.0) j=j+1
-      END DO
-      fm=fm/j
-      nm=nm/j
-
-c
-c AFFICHAGE DU SPECTRE EN FLUX
-c
-
-      DO i=1,nbr
-         y(i)=fx(i)/fm
-      END DO
-      if (ps .EQV. .TRUE.) THEN
-         CALL PGBEGIN (0,'image5.ps/ps',1,1)
-      ELSE
-         CALL PGBEGIN (0,'/xwindow',1,1)
-      END IF
-      CALL PGENV (0.,2.,0.,10.,0,1)
-      CALL PGLABEL('Rapport des frequences recue/emise',
-     &'Flux relatif','SPECTRE DU DISQUE')
-      CALL PGLINE(200,x,y)
-      CALL PGEND
-
-c
-c AFFICHAGE DU SPECTRE EN NOMBRE D'IMPACTS
-c
-
-      DO i=1,nbr
-         y(i)=1.*nt(i)/nm
-      END DO
-      if (ps .EQV. .TRUE.) THEN
-         CALL PGBEGIN (0,'image6.ps/ps',1,1)
-      ELSE
-         CALL PGBEGIN (0,'/xwindow',1,1)
-      END IF
-      CALL PGENV(0.,2.,0.,10.,0,1)
-      CALL PGLABEL('Rapport des frequences recue/emise',
-     &' Nombre relatif d''impacts','SPECTRE DU DISQUE')
-      CALL PGLINE(200,x,y)
-      CALL PGEND
-
-      RETURN
       END
-
 c compilation : 
 c gfortran --std=gnu -o trou_noir_Fortran trou_noir.f -lpgplot
 c SORTIE GRAPHIQUE : -> Envoi sur l'ecran : /xwindow
