Centre Blaise Pascal » Bench4GPU

      PROGRAM MODELISATION_NUMERIQUE

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c                                                                            c

c MODELISATION NUMERIQUE D'UN DISQUE D'ACCRETION AUTOUR D'UN TROU NOIR :     c

c                                                                            c

c APPARENCE & SPECTRE                                                        c

c                                                                            c

c Programme realise par : Herve AUSSEL                                       c

c                         Emmanuel QUEMENER                                  c

c                                                                            c

c Dans le cadre de l'unite de modelisation numerique du DEA                  c

c d'Astrophysique & Techniques Spatiales (Paris 7 & 11)                      c

c                                                                            c

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      IMPLICIT NONE

      REAL*8 pi

      PARAMETER (pi=3.141592654)

c

c PARAMETRES PHYSIQUES DU SYSTEME DISQUE, TROU-NOIR, OBSERVATEUR

c

      REAL*8 m                  ! masse reduite du trou noir

      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*4 bss                ! borne superieure du spectre

c

c PARAMETRES DE LA PLAQUE PHOTOGRAPHIQUE

c

      INTEGER*4 dim             ! dimension de la plaque photographique (en pixels)

      PARAMETER (dim=512)

      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

      REAL*8 bmx                ! parametre d'impact maximal

c

c PARAMETRES INITIAUX D' INTEGRATION

c

      REAL*8 db                 ! avancement du parametre d'impact

      REAL*8 b                  ! parametre d'impact du photon

      REAL*8 h                  ! avancement angulaire initial

c

c PARAMETRES D'INTEGRATION DU SYSTEME DIFFERENTIEL

c

      REAL*8 up,vp,pp           ! parametres du SED , precedents

      REAL*8 us,vs,ps           ! parametres du SED , suivants

c

c RESULTATS DE L'INTEGRATION DE L'EQUATION DIFFERENTIELLE

c

      REAL*8 rp(2000)           ! vecteur polaire

      REAL*8 nh                 ! nombre d'avancements angulaires

c

c RESULTATS DE LA SIMULATION

c

      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...

c

      INTEGER*4 pc              ! pourcentage du balayage realise

      REAL*8 phi                ! angle de la plaque photographique

      REAL*8 thi                ! angle trou noir-observateur & direction emission

      REAL*8 thx                ! angle thi maximal

      REAL*8 phd                ! angle donnant l'intersection des plans

      REAL*8 php                ! angle d'intersection intermediaire

      REAL*8 nr                 ! indice d'intersection reel du vecteur

      INTEGER*4 ni              ! indice d'intersection entier du vecteur

      INTEGER*4 ii              ! indice de 1' image (primaire, secondaire, ...)

      REAL*8 r                  ! rayon moyen d'impact sur le disque

      INTEGER*4 nbr             ! nombre de bandes de rapport de frequence

      INTEGER*4 i,imx,j,n

      LOGICAL tst

      REAL*8 atanp              ! fonction arctangente : resultat entre 0 et 2*pi

      LOGICAL raie            ! variable interactive de type de spectre

      LOGICAL dist            ! variable interactive de distance

c

c TYPE DE SPECTRE

c

      raie=.TRUE.               ! TRUE -> raie infiniment fine

c                               ! FALSE -> corps noir

c

c TYPE DE DISTANCE

c

      dist=.TRUE.               ! TRUE -> distance infinie

c                               ! 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.

            fp(i,j)=0.

         END DO

      END DO

c

c DEFINITION DES PARAMETRES PHYSIQUES

c

      m=1.

      rs=2.*m

      ri=3.*rs

      re=12.

      ro=100.

      tho=pi/180.*80.

      IF (raie .EQV. .TRUE.) THEN

         q=-2

      ELSE

         q=-3/4

      ENDIF

c

c DEFINITION DES PARAMETRES NUMERIQUES

c

      nmx=dim

      stp=dim/(2.*nmx)

      bmx=1.25*re

      b=0.

      thx=asin(bmx/ro)

      pc=0

      nbr=200

      IF (raie .EQV. .TRUE.) THEN

         bss=2.

      ELSE

         bss=1.e19

      ENDIF

c

c BOUCLE PRINCIPALE DE LA SIMULATION

c

      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.

         d=stp*n

         IF (dist .EQV. .TRUE.) THEN

            db=bmx/nmx

            b=db*n

            up=0.

            vp=1.

         ELSE

            thi=thx/nmx*n

            db=ro*DSIN(thi)-b

            b=ro*DSIN(thi)

            up=DSIN(thi)

            vp=DCOS(thi)

         ENDIF

         pp=0.

         nh=1

         CALL rungekutta(ps,us,vs,pp,up,vp,h,m,b)

         rp(nh)=ABS(b/us)

c

c Resolution de l'equation differentielle pour un parametre b donne

c Les resultats sont stockes dans le vecteur polaire rp(n)

c ( Deuxieme etage d'imbrication -> n & nh )

c

         DO WHILE ((rp(nh).GE.rs).AND.(rp(nh).LE.rp(1)))

            nh=nh+1

            pp=ps

            up=us

            vp=vs

            CALL rungekutta(ps,us,vs,pp,up,vp,h,m,b)

            rp(nh)=b/us

         END DO

c

c Mise a zero de tous les termes non redefinis du vecteur polaire

c ( Deuxieme etage d'imbrication -> n & 1 )

c

         DO i=nh+1,2000

            rp(i)=0. 

         END DO

c

c Determination des points de la trajectoire coupant le plan du disque

c par echantillonnage de 1' angle phi sur un cercle, a pas radial constant

c ( Deuxieme etage d'imbrication -> n & 1 )

c

         imx=INT(8*d)

         DO i=0,imx

            phi=2*pi/imx*i

c Calcul de l'angle "brut" d'intersection dans la base polaire

            phd=atanp(DCOS(phi)*DSIN(tho),DCOS(tho))

            phd=phd-INT(phd/pi)*pi

            ii=0

            tst=.FALSE.

c

c Boucle tant que - que l'image tertiaire n'a pas ete teste

c - le rayon n'est pas dans le domaine du disque

c ( Troisieme etage d'imbrication -> n , 1 , (ii, tst) , r )

c

            DO WHILE((ii.LE.2).AND.(tst.EQV..FALSE.))

c Calcul de l'angle d'intersection pour chacune des images

               php=phd+ii*pi

               nr=php/h

               ni=INT(nr)

c Interpolation lineaire de la distance photon-trou noir

               IF (ni.LT.nh) THEN

                  r=(rp(ni+1)-rp(ni))*(nr-ni*1.)+rp(ni)

               ELSE

                  r=rp(ni)

               ENDIF

c

c Test d'impact sur le disque

c ( Quatrieme etage d'imbrication -> n , i , (ii,tst) , r )

c

               IF ((r.LE.re).AND.(r.GE.ri)) then

                  tst=.TRUE.

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)

               END IF

               ii=ii+1

            END DO

         END DO

      END DO

c Appel de la routine d'affichage

      CALL affichage(zp,fp,nt,fx,dim,nbr,bss,bmx,m,ri,re,tho,pi)

      WRITE(*,*)

      WRITE(*,*) 'C''est TERMINE ! '

      WRITE(*,*) 

      STOP

      END

c

c Cette fonction permet de trouver un arctangente entre 0 et 2*pi

c

      REAL*8 function atanp(x,y)

      REAL*8 x,y,pi,eps

      pi=3.141592654

      eps=1.e-20

      IF (ABS(x).LE.eps) then

         IF (ABS(y).EQ.Y) then

            atanp=pi/2.

         ELSE

            atanp=3.*pi/2.

         END IF

      ELSE

         IF (ABS(x).EQ.x) then

            IF (ABS(y).EQ.Y) then

               atanp=DATAN(y/x)

            ELSE

               atanp=DATAN(y/x)+2.*pi

            END IF

         ELSE

            atanp=DATAN(y/x)+pi

         END IF

      END IF

      RETURN

      END

c

c Premiere fonction du systeme -> ( d(u)/d(phi)=f(phi,u,v) )

c

      REAL*8 function f(v)

      REAL*8 v

      f=v

      RETURN

      END

c

c Deuxieme fonction du systeme -> ( d(v)/d(phi)=g(phi,u,v) )

c

      REAL*8 function g(u,m,b)

      REAL*8 u, m, b

      g=3.*m/b*u**2-u

      RETURN

      END

c

c Routine de d'intergration par la methode de Runge-Kutta

c

      SUBROUTINE rungekutta(ps,us,vs,pp,up,vp,h,m,b)

      REAL*8 ps,us,vs,pp,up,vp,h,m,b

      CALL calcul(us,vs,up,vp,h,m,b)

      ps=pp+h

      RETURN

      END

      SUBROUTINE calcul(us,vs,up,vp,h,m,b)

      REAL*8 us,vs,up,vp,h,m,b

      REAL*8 c(4),d(4),f,g

      c(1)=h*f(vp)

      c(2)=h* f(vp+c(1)/2.)

      c(3)=h*f(vp+c(2)/2.)

      c(4)=h*f(vp+c(3))

      d(1)=h*g(up,m,b)

      d(2)=h*g(up+d(1)/2.,m,b)

      d(3)=h*g(up+d(2)/2.,m,b)

      d(4)=h*g(up+d(3),m,b)

      us=up+(c(1)+2.*c(2)+2.*c(3)+c(4))/6.

      vs=vp+(d(1)+2.*d(2)+2.*d(3)+d(4))/6.

      RETURN

      END

c

c En cas d'impact, cette procedure :

c - etablit la position du pixel consisdere sur les images

c - appelle la routine de calcul de decalage spectral

c - appelle la routine de calcul du flux et du spectre

c - affecte a chacun des pixels sa valeur calculee

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

      REAL*4 zp(dim,dim),fp(dim,dim),bss

      REAL*4 xe,ye

      INTEGER*4 xi,yi

      REAL*8 flx,rf

      LOGICAL raie

      xe=d*DSIN(phi)

      ye=-d*DCOS(phi)

c Calcul des coordonnees (x,y) du pixel sur chacune des images

      xi=INT(xe)+INT(dim/2.)

      yi=INT(ye)+INT(dim/2.)

c Appel de la routine de decalage spectral

      CALL decalage_spectral(rf,r,b,phi,tho,m)

c Appel de la routine de calcul de flux et de modification de spectre

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)

      ELSE

         CALL spectre_bb(nt,fx,rf,q,b,db,h,r,m,nbr,bss,flx)

      END IF

c Affectation sur chacune des images du decalage spectral et du flux

      IF(zp(xi,yi).EQ.0.) zp(xi,yi)=1./rf

      IF(fp(xi,yi).EQ.0.) fp(xi,yi)=flx

      RETURN

      END

c Calcul du rapport entre la frequence recue et la frequence emise

      SUBROUTINE decalage_spectral(rf,r,b,phi,tho,m)

      REAL*8 rf,r,b,phi,tho,m

      rf=sqrt(1-3*m/r)/(1+sqrt(m/r**3)*b*sin(tho)*sin(phi))

      RETURN

      END

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

      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

      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 fi,posfreq

      REAL*8 planck,c,k

      PARAMETER (planck=6.62e-34) ! definition des constantes

      PARAMETER (c=3.e8)

      PARAMETER (k=1.38e-23)

c

c Definition des parametres pour le calcul de la temperature d'emission

c puis calcul

c

      temp=3.e7

      m_point=1.

c v -> C du rapport

      v=1.-3./r

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=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

c 

c initialisation des compteurs de flux

c flux int : flux recu pour le pixel dans la tranque posfreq

c flx : flux "bolometrique", integre de 0. a bss

c

      flux_int=0.

      flx=0.

c

c Balayage en frequence du spectre

c

      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

c

c test pour savoir si la frequence de reception est bien dans le

c domaine du spectre calcule

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.))

c Integration sur le pixel

            flux_int=flux_int*rf**3*b*db*h

c Remplissage du spectre

            fx(posfreq)=flux_int+fx(posfreq)

c Intergration bolometrique

            flx=flx+flux_int

         ENDIF

      END DO

      RETURN

      END

c

c AFFICHAGE DES RESULTATS : 6 fenetres de sortie

c

c - Image en decalage spectral

c - Image en flux

c - Dessin de parametres physiques interessants

c - Dessin des limites du disque a l'infini

c - Dessin du spectre en flux

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

      REAL*4 zp(dim,dim),fp(dim,dim),bss

      REAL*8 fx(200),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

      LOGICAL tg

      LOGICAL ps

      tg=.TRUE.                 ! TRUE : affichage en tons de gris

c                               ! FALSE : affichage en contours

      ps=.FALSE.                ! TRUE : affichage postscript 

c                               ! FALSE : affichage fenetre Xwindow

c

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

      fmx=0.

      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

      DATA tr/0.,1.,0.,0.,0.,1./

c

c AFFICHAGE DE L'IMAGE DES DECALAGES SPECTRAUX

c

      if (ps .EQV. .TRUE.) THEN

         CALL PGBEGIN (0,'image1.ps/ps',1,1)

      ELSE

         CALL PGBEGIN (0,'/xwindow',1,1)

      END IF

      CALL PGENV (1.,dim/1.,1.,dim/1.,1,-1)

      CALL PGLABEL('Image des decalages spectraux','',

     &'DISQUE D''ACCRETION AUTOUR D''UN TROU NOIR')

      IF (tg .EQV. .TRUE.) THEN

         DO i=0,25

            DO j=0,25

               zp(i+25,j+25)=1.

            END DO

         END DO

         CALL PGGRAY(zp,dim,dim,2,dim-1,2,dim-1,zmx,0.,tr)

      ELSE

         DO i=1,20

            c(i)=0.2*(i-1)

         END DO

         CALL PGCONS(zp,dim,dim,1,dim,1,dim,c,12,tr)

      END IF

      CALL PGEND

c

c AFFICHAGE DE L'IMAGE DES FLUX

c

      if (ps .EQV. .TRUE.) THEN

         CALL PGBEGIN (0,'image2.ps/ps',1,1)

      ELSE

         CALL PGBEGIN (0,'/xwindow',1,1)

      END IF

      CALL PGENV(1.,dim/1.,1.,dim/1.,1,-1)

      CALL PGLABEL(' Image des flux','',

     &'DISQUE D''ACCRETION AUTOUR D''UN TROU NOIR')

      if (tg) THEN

         DO i=0,25

            DO j=0,25

               fp(i+25,j+25)=1.

            END DO

         END DO

         CALL PGGRAY(fp,dim,dim,2,dim-1,2,dim-1,fmx,0.,tr)

      ELSE

         DO i=1,8

            c(i)=0.0625*2.**(i-1)

         END DO

         CALL PGCONS(fp,dim,dim,1,dim,1,dim,c,8,tr)

      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

c ( dans PGBEGIN ) -> Envoi en format postscript : nom image.ps/ps

/*

	Programme original realise en Fortran 77 en mars 1994

	pour les Travaux Pratiques de Modelisation Numerique

	DEA d'astrophysique et techniques spatiales de Meudon

		par Herve Aussel et Emmanuel Quemener

	Conversion en C par Emmanuel Quemener en aout 1997

	Remerciements a :

	- Herve Aussel pour sa procedure sur le spectre de corps noir

	- Didier Pelat pour l'aide lors de ce travail

	- Jean-Pierre Luminet pour son article de 1979

	- Le Numerical Recipies pour ses recettes de calcul

	- Luc Blanchet pour sa disponibilite lors de mes interrogations en RG

	Mes Coordonnees :	Emmanuel Quemener

				Departement Optique

				ENST de Bretagne

				BP 832

				29285 BREST Cedex

				Emmanuel.Quemener@enst-bretagne.fr

	Compilation sous gcc ( Compilateur GNU sous Linux ) :

		gcc -O6 -m486 -o trou_noir trou_noir.c -lm

*/ 

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#include <string.h>

#define nbr 200 /* Nombre de colonnes du spectre */

#define PI 3.14159265359

float atanp(float x,float y)

{

  float angle;

  angle=atan2(y,x);

  if (angle<0)

    {

      angle+=2*PI;

    }

  return angle;

}

float f(float v)

{

  return v;

}

float g(float u,float m,float b)

{

  return (3.*m/b*pow(u,2)-u);

}

void calcul(float *us,float *vs,float up,float vp,

	    float h,float m,float b)

{

  float c[4],d[4];

  c[0]=h*f(vp);

  c[1]=h*f(vp+c[0]/2.);

  c[2]=h*f(vp+c[1]/2.);

  c[3]=h*f(vp+c[2]);

  d[0]=h*g(up,m,b);

  d[1]=h*g(up+d[0]/2.,m,b);

  d[2]=h*g(up+d[1]/2.,m,b);

  d[3]=h*g(up+d[2],m,b);

  *us=up+(c[0]+2.*c[1]+2.*c[2]+c[3])/6.;

  *vs=vp+(d[0]+2.*d[1]+2.*d[2]+d[3])/6.;

}

void rungekutta(float *ps,float *us,float *vs,

		float pp,float up,float vp,

		float h,float m,float b)

{

  calcul(us,vs,up,vp,h,m,b);

  *ps=pp+h;

}

float decalage_spectral(float r,float b,float phi,

			 float tho,float m)

{

  return (sqrt(1-3*m/r)/(1+sqrt(m/pow(r,3))*b*sin(tho)*sin(phi)));

}

/* void spectre(float rf,int q,float b,float db, */

/* 	     float h,float r,float m,float bss,float *flx) */

float spectre(float rf,int q,float b,float db,

	     float h,float r,float m,float bss)

{

  float flx;

  int fi;

  fi=(int)(rf*nbr/bss);

  flx=pow(r/m,q)*pow(rf,4)*b*db*h;

  return(flx);

}

/* void spectre_cn(float rf,float b,float db, */

/* 		float h,float r,float m,float bss,float *flx) */

float spectre_cn(float rf,float b,float db,

		float h,float r,float m,float bss)

{

  float flx;

  float nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k;

  int fi,posfreq;

  planck=6.62e-34;

  k=1.38e-23;

  temp=3.e7;

  m_point=1.e14;

  v=1.-3./r;

  qu=1/sqrt(1-3./r)/sqrt(r)*(sqrt(r)-sqrt(6)+sqrt(3)/2*log((sqrt(r)+sqrt(3))/(sqrt(r)-sqrt(3))*(sqrt(6)-sqrt(3))/(sqrt(6)+sqrt(3))));

  temp_em=temp*sqrt(m)*exp(0.25*log(m_point))*exp(-0.75*log(r))*

    exp(-0.125*log(v))*exp(0.25*log(qu));

  flux_int=0;

  flx=0;

  for (fi=1;fi<nbr;fi++)

    {

      nu_em=bss*fi/nbr;

      nu_rec=nu_em*rf; 

      posfreq=1./bss*nu_rec*nbr;

      if ((posfreq>0)&&(posfreq<nbr))

	{

	  flux_int=2*planck/9e16*pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b*db*h;

	  flx+=flux_int;

	}

    }

  return(flx);

}

void impact(float d,float phi,int dim,float r,float b,float tho,float m,

	    float **zp,float **fp,

	    int q,float db,

	    float h,float bss,int raie)

{

  float xe,ye;

  int xi,yi;

  float flx,rf;

  xe=d*sin(phi);

  ye=-d*cos(phi);

  xi=(int)xe+dim/2;

  yi=(int)ye+dim/2;

  rf=decalage_spectral(r,b,phi,tho,m);

  if (raie==0)

    {

      flx=spectre(rf,q,b,db,h,r,m,bss);

    }

  else

    {

      flx=spectre_cn(rf,b,db,h,r,m,bss);

    }

  if (zp[xi][yi]==0.)

    {

      zp[xi][yi]=1./rf;

    }

  if (fp[xi][yi]==0.)

    {

      fp[xi][yi]=flx;

    }

}

void sauvegarde_pgm(char nom[24],unsigned int **image,int dim)

{

  FILE            *sortie;

  unsigned long   i,j;

  sortie=fopen(nom,"w");

  fprintf(sortie,"P5\n");

  fprintf(sortie,"%i %i\n",dim,dim);

  fprintf(sortie,"255\n");

  for (j=0;j<dim;j++) for (i=0;i<dim;i++)

    {

      fputc(image[i][j],sortie);

    }

  fclose(sortie);

}

int main(int argc,char *argv[])

{

  float m,rs,ri,re,tho,ro;

  int q;

  float bss,stp;

  int nmx,dim;

  float d,bmx,db,b,h;

  float up,vp,pp;

  float us,vs,ps;

  float rp[2000];

  float **zp,**fp;

  unsigned int **izp,**ifp;

  float zmx,fmx,zen,fen;

  float flux_tot,impc_tot;

  float phi,thi,thx,phd,php,nr,r;

  int ni,ii,i,imx,j,n,tst,dist,raie,pc,fcl,zcl;

  float nh;

  if (argc==2)

    {

      if (strcmp(argv[1],"-aide")==0)

	{

	  printf("\nSimulation d'un disque d'accretion autour d'un trou noir\n");

	  printf("\nParametres a definir :\n\n");

	  printf("   1) Dimension de l'Image\n");

	  printf("   2) Masse relative du trou noir\n");

	  printf("   3) Dimension du disque exterieur\n");

	  printf("   4) Distance de l'observateur\n");

	  printf("   5) Inclinaison par rapport au disque (en degres)\n");

	  printf("   6) Observation a distance FINIE ou INFINIE\n");

	  printf("   7) Rayonnement de disque MONOCHROMATIQUE ou CORPS_NOIR\n");

	  printf("   8) Normalisation des flux INTERNE ou EXTERNE\n");

	  printf("   9) Normalisation de z INTERNE ou EXTERNE\n"); 

	  printf("  10) Impression des images NEGATIVE ou POSITIVE\n"); 

	  printf("  11) Nom de l'image des Flux\n");

	  printf("  12) Nom de l'image des decalages spectraux\n");

	  printf("  13) Valeur de normalisation des flux\n");

	  printf("  14) Valeur de normalisation des decalages spectraux\n");

	  printf("\nSi aucun parametre defini, parametres par defaut :\n\n");

	  printf("   1) Dimension de l'image : 256 pixels de cote\n");

	  printf("   2) Masse relative du trou noir : 1\n");

	  printf("   3) Dimension du disque exterieur : 12 \n");

	  printf("   4) Distance de l'observateur : 100 \n");

	  printf("   5) Inclinaison par rapport au disque (en degres) : 10\n");

	  printf("   6) Observation a distance FINIE\n");

	  printf("   7) Rayonnement de disque MONOCHROMATIQUE\n");

	  printf("   8) Normalisation des flux INTERNE\n");

	  printf("   9) Normalisation des z INTERNE\n");

	  printf("  10) Impression des images NEGATIVE ou POSITIVE\n"); 

       	  printf("  11) Nom de l'image des flux : flux.pgm\n");

	  printf("  12) Nom de l'image des z : z.pgm\n");

	  printf("  13) <non definie>\n");

	  printf("  14) <non definie>\n");

	}

    }

  if ((argc==13)||(argc==15))

    {

      printf("# Utilisation les valeurs definies par l'utilisateur\n");

      dim=atoi(argv[1]);

      m=atof(argv[2]);

      re=atof(argv[3]);

      ro=atof(argv[4]);

      tho=PI/180.*(90-atof(argv[5]));

      rs=2.*m;

      ri=3.*rs;

      q=-2;

      if (strcmp(argv[6],"FINIE")==0)

	{

	  dist=0;

	}

      else

	{

	  dist=1;

	}

      if (strcmp(argv[7],"MONOCHROMATIQUE")==0)

	{

	  raie=0;

	}

      else

	{

	  raie=1;

	}

      if (strcmp(argv[8],"EXTERNE")==0)

	{

	  fen=atof(argv[14]);

	}

      if (strcmp(argv[9],"EXTERNE")==0)

	{

	  zen=atof(argv[15]);

	}

    }

  else

    {

      printf("# Utilisation les valeurs par defaut\n");

      dim=256;

      m=1.;

      rs=2.*m;

      ri=3.*rs;

      re=12.;

      ro=100.;

      tho=PI/180.*80;

      q=-2;

      dist=0;

      raie=0;

    }

      printf("# Dimension de l'image : %i\n",dim);

      printf("# Masse : %f\n",m);

      printf("# Rayon singularite : %f\n",rs);

      printf("# Rayon interne : %f\n",ri);

      printf("# Rayon externe : %f\n",re);

      printf("# Distance de l'observateur : %f\n",ro);

      printf("# Inclinaison a la normale en radian : %f\n",tho);

  zp=(float**)calloc(dim,sizeof(float*));

  zp[0]=(float*)calloc(dim*dim,sizeof(float));

  fp=(float**)calloc(dim,sizeof(float*));

  fp[0]=(float*)calloc(dim*dim,sizeof(float));

  izp=(unsigned int**)calloc(dim,sizeof(unsigned int*));

  izp[0]=(unsigned int*)calloc(dim*dim,sizeof(unsigned int));

  ifp=(unsigned int**)calloc(dim,sizeof(unsigned int*));

  ifp[0]=(unsigned int*)calloc(dim*dim,sizeof(unsigned int));

  for (i=1;i<dim;i++)

    {

      zp[i]=zp[i-1]+dim;

      fp[i]=fp[i-1]+dim;

      izp[i]=izp[i-1]+dim;

      ifp[i]=ifp[i-1]+dim;

    }

  nmx=dim;

  stp=dim/(2.*nmx);

  bmx=1.25*re;

  b=0.;

  thx=asin(bmx/ro);

  pc=0;

  if (raie==0)

    {

      bss=2;

    }

  else

    {

      bss=3e21;

    }

  for (n=1;n<=nmx;n++)

    {     

      h=PI/500.;

      d=stp*n;

      if (dist==1)

	{

	  db=bmx/(float)nmx;

	  b=db*(float)n;

	  up=0.;

	  vp=1.;

	}

      else

	{

	  thi=thx/(float)nmx*(float)n;

	  db=ro*sin(thi)-b;

	  b=ro*sin(thi);

	  up=sin(thi);

	  vp=cos(thi);

	}

      pp=0.;

      nh=1;

      rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);

      rp[(int)nh]=fabs(b/us);

      do

	{

	  nh++;

	  pp=ps;

	  up=us;

	  vp=vs;

	  rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);

	  rp[(int)nh]=b/us;

	} while ((rp[(int)nh]>=rs)&&(rp[(int)nh]<=rp[1]));

      for (i=nh+1;i<2000;i++)

	{

	  rp[i]=0.; 

	}

      imx=(int)(8*d);

      for (i=0;i<=imx;i++)

	{

	  phi=2.*PI/(float)imx*(float)i;

	  phd=atanp(cos(phi)*sin(tho),cos(tho));

	  phd=fmod(phd,PI);

	  ii=0;

	  tst=0;

	  do

	    {

	      php=phd+(float)ii*PI;

	      nr=php/h;

	      ni=(int)nr;

	      if ((float)ni<nh)

		{

		  r=(rp[ni+1]-rp[ni])*(nr-ni*1.)+rp[ni];

		}

	      else

		{

		  r=rp[ni];

		}

	      if ((r<=re)&&(r>=ri))

		{

		  tst=1;

		  impact(d,phi,dim,r,b,tho,m,zp,fp,q,db,h,bss,raie);

		}

	      ii++;

	    } while ((ii<=2)&&(tst==0));

	}

    }

  fmx=fp[0][0];

  zmx=zp[0][0];

  for (i=0;i<dim;i++) for (j=0;j<dim;j++)

    {

      if (fmx<fp[i][j])

	{

	  fmx=fp[i][j];

	}

      if (zmx<zp[i][j])

	{

	  zmx=zp[i][j];

	}

    }

  printf("\nLe flux maximal detecte est de %f",fmx);

  printf("\nLe decalage spectral maximal detecte est de %f\n\n",zmx);

  if (strcmp(argv[8],"EXTERNE")==0)

    {

      fmx=fen;

    }

  if (strcmp(argv[9],"EXTERNE")==0)

    {  

      zmx=zen;

    }

  for (i=0;i<dim;i++) for (j=0;j<dim;j++)

    {

      zcl=(int)(255/zmx*zp[i][dim-1-j]);

      fcl=(int)(255/fmx*fp[i][dim-1-j]);

      if (strcmp(argv[8],"NEGATIVE")==0)

	{

	  if (zcl>255)

	    {

	      izp[i][j]=0;

	    }

	  else

	    {

	      izp[i][j]=255-zcl;

	    }

	  if (fcl>255)

	    {

	      ifp[i][j]=0;

	    }

	  else

	    {

	      ifp[i][j]=255-fcl;

	    } 

	}

      else

	{

	  if (zcl>255)

	    {

	      izp[i][j]=255;

	    }

	  else

	    {

	      izp[i][j]=zcl;

	    }

	  if (fcl>255)

	    {

	      ifp[i][j]=255;

	    }

	  else

	    {

	      ifp[i][j]=fcl;

	    } 

	}

    }

  if ((argc==14)||(argc==16))

   {

     sauvegarde_pgm(argv[11],ifp,dim);

     sauvegarde_pgm(argv[12],izp,dim);

   }

  else

    {

      sauvegarde_pgm("z.pgm",izp,dim);

      sauvegarde_pgm("flux.pgm",ifp,dim);

    }

  free(zp[0]);

  free(zp);

  free(fp[0]);

  free(fp);

  free(izp[0]);

  free(izp);

  free(ifp[0]);

  free(ifp);

}

/*

	Programme original realise en Fortran 77 en mars 1994

	pour les Travaux Pratiques de Modelisation Numerique

	DEA d'astrophysique et techniques spatiales de Meudon

		par Herve Aussel et Emmanuel Quemener

	Conversion en C par Emmanuel Quemener en aout 1997

	Remerciements a :

	- Herve Aussel pour sa procedure sur le spectre de corps noir

	- Didier Pelat pour l'aide lors de ce travail

	- Jean-Pierre Luminet pour son article de 1979

	- Le Numerical Recipies pour ses recettes de calcul

	- Luc Blanchet pour sa disponibilite lors de mes interrogations en RG

	Mes Coordonnees :	Emmanuel Quemener

				Departement Optique

				ENST de Bretagne

				BP 832

				29285 BREST Cedex

				Emmanuel.Quemener@enst-bretagne.fr

	Compilation sous gcc ( Compilateur GNU sous Linux ) :

		gcc -O6 -m486 -o trou_noir trou_noir.c -lm

*/ 

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#include <string.h>

#define nbr 200 /* Nombre de colonnes du spectre */

#define PI 3.14159265359

float atanp(float x,float y)

{

  float angle;

  angle=atan2(y,x);

  if (angle<0)

    {

      angle+=2*PI;

    }

  return angle;

}

float f(float v)

{

  return v;

}

float g(float u,float m,float b)

{

  return (3.*m/b*pow(u,2)-u);

}

void calcul(float *us,float *vs,float up,float vp,

	    float h,float m,float b)

{

  float c[4],d[4];

  c[0]=h*f(vp);

  c[1]=h*f(vp+c[0]/2.);

  c[2]=h*f(vp+c[1]/2.);

  c[3]=h*f(vp+c[2]);

  d[0]=h*g(up,m,b);

  d[1]=h*g(up+d[0]/2.,m,b);

  d[2]=h*g(up+d[1]/2.,m,b);

  d[3]=h*g(up+d[2],m,b);

  *us=up+(c[0]+2.*c[1]+2.*c[2]+c[3])/6.;

  *vs=vp+(d[0]+2.*d[1]+2.*d[2]+d[3])/6.;

}

void rungekutta(float *ps,float *us,float *vs,

		float pp,float up,float vp,

		float h,float m,float b)

{

  calcul(us,vs,up,vp,h,m,b);

  *ps=pp+h;

}

float decalage_spectral(float r,float b,float phi,

			 float tho,float m)

{

  return (sqrt(1-3*m/r)/(1+sqrt(m/pow(r,3))*b*sin(tho)*sin(phi)));

}

/* void spectre(float rf,int q,float b,float db, */

/* 	     float h,float r,float m,float bss,float *flx) */

float spectre(float rf,int q,float b,float db,

	     float h,float r,float m,float bss)

{

  float flx;

  int fi;

  fi=(int)(rf*nbr/bss);

  flx=pow(r/m,q)*pow(rf,4)*b*db*h;

  return(flx);

}

/* void spectre_cn(float rf,float b,float db, */

/* 		float h,float r,float m,float bss,float *flx) */

float spectre_cn(float rf,float b,float db,

		float h,float r,float m,float bss)

{

  float flx;

  float nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k;

  int fi,posfreq;

  planck=6.62e-34;

  k=1.38e-23;

  temp=3.e7;

  m_point=1.e14;

  v=1.-3./r;

  qu=1/sqrt(1-3./r)/sqrt(r)*(sqrt(r)-sqrt(6)+sqrt(3)/2*log((sqrt(r)+sqrt(3))/(sqrt(r)-sqrt(3))*(sqrt(6)-sqrt(3))/(sqrt(6)+sqrt(3))));

  temp_em=temp*sqrt(m)*exp(0.25*log(m_point))*exp(-0.75*log(r))*

    exp(-0.125*log(v))*exp(0.25*log(qu));

  flux_int=0;

  flx=0;

  for (fi=1;fi<nbr;fi++)

    {

      nu_em=bss*fi/nbr;

      nu_rec=nu_em*rf; 

      posfreq=1./bss*nu_rec*nbr;

      if ((posfreq>0)&&(posfreq<nbr))

	{

	  flux_int=2*planck/9e16*pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b*db*h;

	  flx+=flux_int;

	}

    }

  return(flx);

}

void impact(float d,float phi,int dim,float r,float b,float tho,float m,

	    float **zp,float **fp,

	    int q,float db,

	    float h,float bss,int raie)

{

  float xe,ye;

  int xi,yi;

  float flx,rf;

  xe=d*sin(phi);

  ye=-d*cos(phi);

  xi=(int)xe+dim/2;

  yi=(int)ye+dim/2;

  rf=decalage_spectral(r,b,phi,tho,m);

  if (raie==0)

    {

      flx=spectre(rf,q,b,db,h,r,m,bss);

    }

  else

    {

      flx=spectre_cn(rf,b,db,h,r,m,bss);

    }

  if (zp[xi][yi]==0.)

    {

      zp[xi][yi]=1./rf;

    }

  if (fp[xi][yi]==0.)

    {

      fp[xi][yi]=flx;

    }

}

void sauvegarde_pgm(char nom[24],unsigned int **image,int dim)

{

  FILE            *sortie;

  unsigned long   i,j;

  sortie=fopen(nom,"w");

  fprintf(sortie,"P5\n");

  fprintf(sortie,"%i %i\n",dim,dim);

  fprintf(sortie,"255\n");

  for (j=0;j<dim;j++) for (i=0;i<dim;i++)

    {