Index: trou_noir.c
===================================================================
--- trou_noir.c	(revision 213)
+++ trou_noir.c	(revision 214)
@@ -6,6 +6,7 @@
 		par Herve Aussel et Emmanuel Quemener
 
 	Conversion en C par Emmanuel Quemener en aout 1997
+        Modification par Emmanuel Quemener en aout 2019
 
 	Remerciements a :
 
@@ -15,31 +16,32 @@
 	- 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
+		gcc -O2 -o trou_noir trou_noir.c -lm
 */ 
 
 #include <stdio.h>
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
+#include <sys/time.h>
 
-#define nbr 200 /* Nombre de colonnes du spectre */
+#define nbr 256 /* Nombre de colonnes du spectre */
 
 #define PI 3.14159265359
 
-double atanp(double x,double y)
+#define TRACKPOINTS 2048
+
+#if TYPE == FP32
+#define MYFLOAT float
+#else
+#define MYFLOAT double
+#endif
+
+MYFLOAT atanp(MYFLOAT x,MYFLOAT y)
 {
-  double angle;
+  MYFLOAT angle;
 
   angle=atan2(y,x);
 
@@ -52,22 +54,21 @@
 }
 
 
-double f(double v)
+MYFLOAT f(MYFLOAT v)
 {
   return v;
 }
 
-double g(double u,double m,double b)
+MYFLOAT g(MYFLOAT u,MYFLOAT m,MYFLOAT b)
 {
-// return (3.*m/b*pow(u,2)-u);
   return (3.*m/b*pow(u,2)-u);
 }
 
 
-void calcul(double *us,double *vs,double up,double vp,
-	    double h,double m,double b)
+void calcul(MYFLOAT *us,MYFLOAT *vs,MYFLOAT up,MYFLOAT vp,
+	    MYFLOAT h,MYFLOAT m,MYFLOAT b)
 {
-  double c[4],d[4];
+  MYFLOAT c[4],d[4];
 
   c[0]=h*f(vp);
   c[1]=h*f(vp+c[0]/2.);
@@ -82,75 +83,80 @@
   *vs=vp+(d[0]+2.*d[1]+2.*d[2]+d[3])/6.;
 }
 
-void rungekutta(double *ps,double *us,double *vs,
-		double pp,double up,double vp,
-		double h,double m,double b)
+void rungekutta(MYFLOAT *ps,MYFLOAT *us,MYFLOAT *vs,
+		MYFLOAT pp,MYFLOAT up,MYFLOAT vp,
+		MYFLOAT h,MYFLOAT m,MYFLOAT b)
 {
   calcul(us,vs,up,vp,h,m,b);
   *ps=pp+h;
 }
 
 
-double decalage_spectral(double r,double b,double phi,
-			 double tho,double m)
+MYFLOAT decalage_spectral(MYFLOAT r,MYFLOAT b,MYFLOAT phi,
+			 MYFLOAT tho,MYFLOAT m)
 {
   return (sqrt(1-3*m/r)/(1+sqrt(m/pow(r,3))*b*sin(tho)*sin(phi)));
 }
 
-void spectre(int nt[nbr],double fx[nbr],double rf,int q,double b,double db,
-	     double h,double r,double m,double bss,double *flx)
+MYFLOAT spectre(MYFLOAT rf,MYFLOAT q,MYFLOAT b,MYFLOAT db,
+	     MYFLOAT h,MYFLOAT r,MYFLOAT m,MYFLOAT bss)
 {
-  int fi;
+  MYFLOAT flx;
 
-  fi=(int)(rf*nbr/bss);
-  nt[fi]+=1;
-  *flx=pow(r/m,q)*pow(rf,4)*b*db*h;
-  fx[fi]=fx[fi]+*flx;
+  flx=exp(q*log(r/m))*pow(rf,4)*b*db*h;
+  return(flx);
 }
 
-void spectre_cn(int nt[nbr],double fx[nbr],double rf,double b,double db,
-		double h,double r,double m,double bss,double *flx)
+MYFLOAT spectre_cn(MYFLOAT rf,MYFLOAT b,MYFLOAT db,
+		MYFLOAT h,MYFLOAT r,MYFLOAT m,MYFLOAT bss)
 {
-  double nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k;
+  
+  MYFLOAT flx;
+  MYFLOAT nu_rec,nu_em,qu,temp_em,flux_int;
   int fi,posfreq;
 
-  planck=6.62e-34;
-  k=1.38e-23;
-  temp=3.e7;
-  m_point=1.e14;
-  v=1.-3./r;
+#define planck 6.62e-34
+#define k 1.38e-23
+#define c2 9.e16
+#define temp 3.e7
+#define m_point 1.
 
-  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))));
+#define lplanck (log(6.62)-34.*log(10.))
+#define lk (log(1.38)-23.*log(10.))
+#define lc2 (log(9.)+16.*log(10.))
+  
+  MYFLOAT v=1.-3./r;
 
-  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));
+  qu=1./sqrt((1.-3./r)*r)*(sqrt(r)-sqrt(6.)+sqrt(3.)/2.*log((sqrt(r)+sqrt(3.))/(sqrt(r)-sqrt(3.))* 0.17157287525380988 ));
 
-  flux_int=0;
-  *flx=0;
+  temp_em=temp*sqrt(m)*exp(0.25*log(m_point)-0.75*log(r)-0.125*log(v)+0.25*log(fabs(qu)));
 
-  for (fi=1;fi<nbr;fi++)
+  flux_int=0.;
+  flx=0.;
+
+  for (fi=0;fi<nbr;fi++)
     {
-      nu_em=bss*fi/nbr;
-      nu_rec=nu_em*rf; 
-      posfreq=1./bss*nu_rec*nbr;
+      nu_em=bss*(MYFLOAT)fi/(MYFLOAT)nbr;
+      nu_rec=nu_em*rf;
+      posfreq=(int)(nu_rec*(MYFLOAT)nbr/bss);
       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;
-	  fx[posfreq]+=flux_int;
-	  *flx+=flux_int;
-	  nt[posfreq]+=1;
-	}
+  	{
+	  flux_int=2.*planck/c2*pow(nu_em,3)/(exp(planck*nu_em/(k*temp_em))-1.)*exp(3.*log(rf))*b*db*h;
+  	  flx+=flux_int;
+  	}
     }
+
+  return((MYFLOAT)flx);
 }
 
-void impact(double d,double phi,int dim,double r,double b,double tho,double m,
-	    double **zp,double **fp,
-	    int nt[200],double fx[200],int q,double db,
-	    double h,double bss,int raie)
+void impact(MYFLOAT d,MYFLOAT phi,int dim,MYFLOAT r,MYFLOAT b,MYFLOAT tho,MYFLOAT m,
+	    MYFLOAT **zp,MYFLOAT **fp,
+	    MYFLOAT q,MYFLOAT db,
+	    MYFLOAT h,MYFLOAT bss,int raie)
 {
-  double xe,ye;
+  MYFLOAT xe,ye;
   int xi,yi;
-  double flx,rf;
+  MYFLOAT flx,rf;
   xe=d*sin(phi);
   ye=-d*cos(phi);
 
@@ -161,11 +167,13 @@
 
   if (raie==0)
     {
-      spectre(nt,fx,rf,q,b,db,h,r,m,bss,&flx);
+      bss=1.e19;
+      flx=spectre_cn(rf,b,db,h,r,m,bss);
     }
   else
     {
-      spectre_cn(nt,fx,rf,b,db,h,r,m,bss,&flx);
+      bss=2.;
+      flx=spectre(rf,q,b,db,h,r,m,bss);
     }
   
   if (zp[xi][yi]==0.)
@@ -177,6 +185,7 @@
     {
       fp[xi][yi]=flx;
     }
+
 }
 
 void sauvegarde_pgm(char nom[24],unsigned int **image,int dim)
@@ -198,123 +207,70 @@
   fclose(sortie);
 }
 
-void sauvegarde_dat(char nom[24],double tableau[3][nbr],int raie)
-{
-  FILE            *sortie;
-  unsigned long   i;
-  
-  sortie=fopen(nom,"w");
-
-  fprintf(sortie,"# Trou Noir entoure d'un Disque d'Accretion\n");
- 
-  if (raie==0)
-    {
-      fprintf(sortie,"# Colonne 1 : Frequence_Recue/Frequence_Emise\n");
-    }
-  else
-    {
-      fprintf(sortie,"# Colonne 1 : Fréquence d'Emission en Hertz\n");
-    }
-
-  fprintf(sortie,"# Colonne 2 : Intensite Normalisee\n");
-  fprintf(sortie,"# Colonne 3 : Nombre d'Impacts Normalise\n");
-  
-  for (i=1;i<nbr;i++)
-    {
-      fprintf(sortie,"%f\t%f\t%f\n",tableau[0][i],tableau[1][i],tableau[2][i]);
-    }
-
-  fclose(sortie);
-}
-
 int main(int argc,char *argv[])
 {
 
-  double m,rs,ri,re,tho,ro;
+  MYFLOAT m,rs,ri,re,tho;
   int q;
 
-  double bss,stp;
+  MYFLOAT bss,stp;
   int nmx,dim;
-  double d,bmx,db,b,h;
-  double up,vp,pp;
-  double us,vs,ps;
-  double rp[2000];
-  double **zp,**fp;
+  MYFLOAT d,bmx,db,b,h;
+  MYFLOAT up,vp,pp;
+  MYFLOAT us,vs,ps;
+  MYFLOAT rp[TRACKPOINTS];
+  MYFLOAT **zp,**fp;
   unsigned int **izp,**ifp;
-  double zmx,fmx,zen,fen;
-  double flux_tot,impc_tot;
-  double fx[nbr];
-  int nt[nbr];
-  double tableau[3][nbr];
-  double phi,thi,thx,phd,php,nr,r;
+  MYFLOAT zmx,fmx,zen,fen;
+  int zimx=0,zjmx=0,fimx=0,fjmx=0;
+  MYFLOAT flux_tot,impc_tot;
+  MYFLOAT phi,thi,thx,phd,php,nr,r;
   int ni,ii,i,imx,j,n,tst,dist,raie,pc,fcl,zcl;
-  double nh;
-
+  MYFLOAT nh;
+  struct timeval tv1,tv2;
+  struct timezone tz;  
+  double elapsed;
+  
   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) Nom du fichier contenant le spectre\n");
-	  printf("  14) Valeur de normalisation des flux\n");
-	  printf("  15) Valeur de normalisation des decalages spectraux\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) Inclinaison par rapport au disque (en degres)\n");
+	  printf("  5) Rayonnement de disque MONOCHROMATIQUE ou CORPS_NOIR\n");
+	  printf("  6) Impression des images NEGATIVE ou POSITIVE\n"); 
+	  printf("  7) Nom de l'image des Flux\n");
+	  printf("  8) Nom de l'image 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) Nom du fichier contenant le spectre : spectre.dat\n");
-	  printf("  14) <non definie>\n");
-	  printf("  15) <non definie>\n");
+	  printf("  1) Dimension de l'image : 1024 pixels de cote\n");
+	  printf("  2) Masse relative du trou noir : 1\n");
+	  printf("  3) Dimension du disque exterieur : 12 \n");
+	  printf("  4) Inclinaison par rapport au disque (en degres) : 10\n");
+	  printf("  5) Rayonnement de disque CORPS_NOIR\n");
+	  printf("  6) Impression des images NEGATIVE ou POSITIVE\n"); 
+       	  printf("  7) Nom de l'image des flux : flux.pgm\n");
+	  printf("  8) Nom de l'image des z : z.pgm\n");
 	}
     }
   
-  if ((argc==14)||(argc==16))
+  if (argc==9)
     {
       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]));
+      tho=PI/180.*(90-atof(argv[4]));
       
       rs=2.*m;
       ri=3.*rs;
-      q=-2;
 
-      if (strcmp(argv[6],"FINIE")==0)
+      if (strcmp(argv[5],"CORPS_NOIR")==0)
 	{
-	  dist=0;
-	}
-      else
-	{
-	  dist=1;
-	}
-
-      if (strcmp(argv[7],"MONOCHROMATIQUE")==0)
-	{
 	  raie=0;
 	}
       else
@@ -322,52 +278,44 @@
 	  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;
+      dim=1024;
       m=1.;
       rs=2.*m;
       ri=3.*rs;
       re=12.;
-      ro=100.;
       tho=PI/180.*80;
-      q=-2;
-      dist=0;
+      // Corps noir
       raie=0;
     }
-      
+
+  if (raie==1)
+    {
+      bss=2.;
+      q=-2;
+    }
+  else
+    {
+      bss=1.e19;
+      q=-0.75;
+    }
+
       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);
   
-  for (i=0;i<nbr;i++)
-    {
-      fx[i]=0.;
-      nt[i]=0;
-    }  
-
-  zp=(double**)calloc(dim,sizeof(double*));
-  zp[0]=(double*)calloc(dim*dim,sizeof(double));
+  zp=(MYFLOAT**)calloc(dim,sizeof(MYFLOAT*));
+  zp[0]=(MYFLOAT*)calloc(dim*dim,sizeof(MYFLOAT));
   
-  fp=(double**)calloc(dim,sizeof(double*));
-  fp[0]=(double*)calloc(dim*dim,sizeof(double));
+  fp=(MYFLOAT**)calloc(dim,sizeof(MYFLOAT*));
+  fp[0]=(MYFLOAT*)calloc(dim*dim,sizeof(MYFLOAT));
 
   izp=(unsigned int**)calloc(dim,sizeof(unsigned int*));
   izp[0]=(unsigned int*)calloc(dim*dim,sizeof(unsigned int));
@@ -381,44 +329,26 @@
       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;
-    }
+  // Set start timer
+  gettimeofday(&tv1, &tz);
 
   for (n=1;n<=nmx;n++)
     {     
-      h=PI/500.;
+      h=4.*PI/(MYFLOAT)TRACKPOINTS;
       d=stp*n;
 
-      if (dist==1)
-	{
-	  db=bmx/(double)nmx;
-	  b=db*(double)n;
-	  up=0.;
-	  vp=1.;
-	}
-      else
-	{
-	  thi=thx/(double)nmx*(double)n;
-	  db=ro*sin(thi)-b;
-	  b=ro*sin(thi);
-	  up=sin(thi);
-	  vp=cos(thi);
-	}
+      db=bmx/(MYFLOAT)nmx;
+      b=db*(MYFLOAT)n;
+      up=0.;
+      vp=1.;
       
       pp=0.;
       nh=1;
@@ -439,7 +369,7 @@
 	  
 	} while ((rp[(int)nh]>=rs)&&(rp[(int)nh]<=rp[1]));
       
-      for (i=nh+1;i<2000;i++)
+      for (i=nh+1;i<TRACKPOINTS;i++)
 	{
 	  rp[i]=0.; 
 	}
@@ -448,7 +378,7 @@
       
       for (i=0;i<=imx;i++)
 	{
-	  phi=2.*PI/(double)imx*(double)i;
+	  phi=2.*PI/(MYFLOAT)imx*(MYFLOAT)i;
 	  phd=atanp(cos(phi)*sin(tho),cos(tho));
 	  phd=fmod(phd,PI);
 	  ii=0;
@@ -456,11 +386,11 @@
 	  
 	  do
 	    {
-	      php=phd+(double)ii*PI;
+	      php=phd+(MYFLOAT)ii*PI;
 	      nr=php/h;
 	      ni=(int)nr;
 
-	      if ((double)ni<nh)
+	      if ((MYFLOAT)ni<nh)
 		{
 		  r=(rp[ni+1]-rp[ni])*(nr-ni*1.)+rp[ni];
 		}
@@ -472,7 +402,7 @@
 	      if ((r<=re)&&(r>=ri))
 		{
 		  tst=1;
-		  impact(d,phi,dim,r,b,tho,m,zp,fp,nt,fx,q,db,h,bss,raie);
+		  impact(d,phi,dim,r,b,tho,m,zp,fp,q,db,h,bss,raie);
 		}
 	      
 	      ii++;
@@ -480,6 +410,12 @@
 	}
     }
 
+  // Set stop timer
+  gettimeofday(&tv2, &tz);
+  
+  elapsed=(double)((tv2.tv_sec-tv1.tv_sec) * 1000000L +
+			  (tv2.tv_usec-tv1.tv_usec))/1000000;  
+  
   fmx=fp[0][0];
   zmx=zp[0][0];
   
@@ -487,34 +423,29 @@
     {
       if (fmx<fp[i][j])
 	{
+	  fimx=i;
+	  fjmx=j;
 	  fmx=fp[i][j];
 	}
       
       if (zmx<zp[i][j])
 	{
+	  zimx=i;
+	  zjmx=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);
+  printf("\nElapsed Time : %lf",(double)elapsed);
+  printf("\nZ max @(%i,%i) : %f",zimx,zjmx,zmx);
+  printf("\nFlux max @(%i,%i) : %f\n\n",fimx,fjmx,fmx);
 
-  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 (strcmp(argv[6],"NEGATIVE")==0)
 	{
 	  if (zcl>255)
 	    {
@@ -559,33 +490,15 @@
 	
     }
 
-  flux_tot=0;
-  impc_tot=0;
-
-  for (i=1;i<nbr;i++)
-    {
-      flux_tot+=fx[i];
-      impc_tot+=nt[i];
-    }
-
-  for (i=1;i<nbr;i++)
-    {
-      tableau[0][i]=bss*i/nbr;
-      tableau[1][i]=fx[i]/flux_tot;
-      tableau[2][i]=(double)nt[i]/(double)impc_tot;
-    }
-
-  if ((argc==14)||(argc==16))
+  if (argc==9)
    {
-     sauvegarde_pgm(argv[11],ifp,dim);
-     sauvegarde_pgm(argv[12],izp,dim);
-     sauvegarde_dat(argv[13],tableau,raie);
+     sauvegarde_pgm(argv[7],ifp,dim);
+     sauvegarde_pgm(argv[8],izp,dim);
    }
   else
     {
       sauvegarde_pgm("z.pgm",izp,dim);
       sauvegarde_pgm("flux.pgm",ifp,dim);
-      sauvegarde_dat("spectre.dat",tableau,raie);
     }
 
   free(zp[0]);
