Index: TrouNoir/TrouNoir.py
===================================================================
--- TrouNoir/TrouNoir.py	(revision 233)
+++ TrouNoir/TrouNoir.py	(revision 234)
@@ -1,12 +1,12 @@
 #!/usr/bin/env python
 #
-# TrouNoir model using PyOpenCL 
+# TrouNoir model using PyOpenCL
 #
-# CC BY-NC-SA 2019 : <emmanuel.quemener@ens-lyon.fr> 
+# CC BY-NC-SA 2019 : <emmanuel.quemener@ens-lyon.fr>
 #
 # Thanks to Andreas Klockner for PyOpenCL:
 # http://mathema.tician.de/software/pyopencl
-# 
+#
 # Original code programmed in Fortran 77 in mars 1994
 # for Practical Work of Numerical Simulation
 # DEA (old Master2) in astrophysics and spatial techniques in Meudon
@@ -17,7 +17,7 @@
 # GPUfication in CUDA under Python in august 2019
 #
 # Thanks to :
-# 
+#
 # - Herve Aussel for his part of code of black body spectrum
 # - Didier Pelat for his help to perform this work
 # - Jean-Pierre Luminet for his article published in 1979
@@ -102,7 +102,7 @@
 
 BlobOpenCL= """
 
-#define PI (float)3.14159265359
+#define PI (float)3.14159265359e0f
 #define nbr 256
 
 #define EINSTEIN 0
@@ -146,51 +146,51 @@
 #endif
 
 void calcul(float *us,float *vs,float up,float vp,
-	    float h,float m,float b)
+            float h,float m,float b)
 {
   float c0,c1,c2,c3,d0,d1,d2,d3;
 
   c0=h*f(vp);
-  c1=h*f(vp+c0/2.);
-  c2=h*f(vp+c1/2.);
+  c1=h*f(vp+c0/2.e0f);
+  c2=h*f(vp+c1/2.e0f);
   c3=h*f(vp+c2);
   d0=h*g(up,m,b);
-  d1=h*g(up+d0/2.,m,b);
-  d2=h*g(up+d1/2.,m,b);
+  d1=h*g(up+d0/2.e0f,m,b);
+  d2=h*g(up+d1/2.e0f,m,b);
   d3=h*g(up+d2,m,b);
 
-  *us=up+(c0+2.*c1+2.*c2+c3)/6.;
-  *vs=vp+(d0+2.*d1+2.*d2+d3)/6.;
+  *us=up+(c0+2.e0f*c1+2.e0f*c2+c3)/6.e0f;
+  *vs=vp+(d0+2.e0f*d1+2.e0f*d2+d3)/6.e0f;
 }
 
 void rungekutta(float *ps,float *us,float *vs,
-		float pp,float up,float vp,
-		float h,float m,float b)
+                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)
+                        float tho,float m)
 {
   return (sqrt(1-3*m/r)/(1+sqrt(m/pow(r,3))*b*sin(tho)*sin(phi)));
 }
 
 float spectre(float rf,int q,float b,float db,
-	     float h,float r,float m,float bss)
+              float h,float r,float m,float bss)
 {
   float flx;
 
 //  flx=exp(q*log(r/m))*pow(rf,4)*b*db*h;
-  flx=exp(q*log(r/m)+4.*log(rf))*b*db*h;
+  flx=exp(q*log(r/m)+4.e0f*log(rf))*b*db*h;
   return(flx);
 }
 
 float spectre_cn(float rf32,float b32,float db32,
-		 float h32,float r32,float m32,float bss32)
+                 float h32,float r32,float m32,float bss32)
 {
-  
+
 #define MYFLOAT float
 
   MYFLOAT rf=(MYFLOAT)(rf32);
@@ -205,24 +205,24 @@
   MYFLOAT nu_rec,nu_em,qu,temp_em,flux_int;
   int fi,posfreq;
 
-#define planck 6.62e-34
-#define k 1.38e-23
-#define c2 9.e16
-#define temp 3.e7
-#define m_point 1.
+#define planck 6.62e-34f
+#define k 1.38e-23f
+#define c2 9.e16f
+#define temp 3.e7f
+#define m_point 1.e0f
 
-#define lplanck (log(6.62)-34.*log(10.))
-#define lk (log(1.38)-23.*log(10.))
-#define lc2 (log(9.)+16.*log(10.))
+#define lplanck (log(6.62e0f)-34.e0f*log(10.e0f))
+#define lk (log(1.38e0f)-23.e0f*log(10.e0f))
+#define lc2 (log(9.e0f)+16.e0f*log(10.e0f))
 
-  MYFLOAT v=1.-3./r;
+  MYFLOAT v=1.e0f-3.e0f/r;
 
-  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 ));
+  qu=1.e0f/sqrt((1.e0f-3.e0f/r)*r)*(sqrt(r)-sqrt(6.e0f)+sqrt(3.e0f)/2.e0f*log((sqrt(r)+sqrt(3.e0f))/(sqrt(r)-sqrt(3.e0f))* 0.17157287525380988e0f ));
 
-  temp_em=temp*sqrt(m)*exp(0.25*log(m_point)-0.75*log(r)-0.125*log(v)+0.25*log(fabs(qu)));
+  temp_em=temp*sqrt(m)*exp(0.25e0f*log(m_point)-0.75e0f*log(r)-0.125e0f*log(v)+0.25e0f*log(fabs(qu)));
 
-  flux_int=0.;
-  flx=0.;
+  flux_int=0.e0f;
+  flx=0.e0f;
 
   for (fi=0;fi<nbr;fi++)
     {
@@ -230,26 +230,26 @@
       nu_rec=nu_em*rf;
       posfreq=(int)(nu_rec*(MYFLOAT)nbr/bss);
       if ((posfreq>0)&&(posfreq<nbr))
-	{
+        {
           // Initial version
-	  // flux_int=2.*planck/c2*pow(nu_em,3)/(exp(planck*nu_em/(k*temp_em))-1.);
+          // flux_int=2.*planck/c2*pow(nu_em,3)/(exp(planck*nu_em/(k*temp_em))-1.);
           // Version with log used
-	  //flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.);
-	  // flux_int*=pow(rf,3)*b*db*h;
-	  //flux_int*=exp(3.*log(rf))*b*db*h;
-	  flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.)*exp(3.*log(rf))*b*db*h;
+          //flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.);
+          // flux_int*=pow(rf,3)*b*db*h;
+          //flux_int*=exp(3.e0f*log(rf))*b*db*h;
+          flux_int=2.e0f*exp(lplanck-lc2+3.e0f*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.e0f)*exp(3.e0f*log(rf))*b*db*h;
 
-	  flx+=flux_int;
-	}
+          flx+=flux_int;
+        }
     }
 
   return((float)(flx));
 }
 
 void impact(float phi,float r,float b,float tho,float m,
-	    float *zp,float *fp,
-	    int q,float db,
-	    float h,int raie)
+            float *zp,float *fp,
+            int q,float db,
+            float h,int raie)
 {
   float flx,rf,bss;
 
@@ -257,16 +257,16 @@
 
   if (raie==0)
     {
-      bss=1.e19;
+      bss=1.e19f;
       flx=spectre_cn(rf,b,db,h,r,m,bss);
     }
   else
-    { 
-      bss=2.;
+    {
+      bss=2.e0f;
       flx=spectre(rf,q,b,db,h,r,m,bss);
     }
-  
-  *zp=1./rf;
+
+  *zp=1.e0f/rf;
   *fp=flx;
 
 }
@@ -295,7 +295,7 @@
   raie=Line;
 
   float bmx,db,b,h;
-  float rp0,rpp,rps;
+  float rp0,rps;
   float phi,phd;
   uint nh=0;
   float zp=0.e0f,fp=0.e0f;
@@ -356,7 +356,7 @@
   barrier(CLK_GLOBAL_MEM_FENCE);
 
   zImage[yi+sizex*xi]=(float)zp;
-  fImage[yi+sizex*xi]=(float)fp;  
+  fImage[yi+sizex*xi]=(float)fp;
 }
 
 __kernel void Pixel(__global float *zImage,__global float *fImage,
@@ -385,20 +385,20 @@
 
   float bmx,db,b,h;
   float phi,phd,php,nr,r;
-  float zp=0.,fp=0.;
+  float zp=0.e0f,fp=0.e0f;
 
   // Autosize for image, 25% greater than external radius
-  bmx=1.25*re;
+  bmx=1.25e0f*re;
 
   // Angular step of integration
   h=4.e0f*PI/(float)TRACKPOINTS;
 
   // Step of Impact Parameter
-  db=bmx/(2.e0*(float)ImpactParameter);
+  db=bmx/(2.e0f*(float)ImpactParameter);
 
   // set origin as center of image
-  float x=(float)xi-(float)(sizex/2)+(float)5e-1f;
-  float y=(float)yi-(float)(sizey/2)+(float)5e-1f;
+  float x=(float)xi-(float)(sizex/2)+(float)5.e-1f;
+  float y=(float)yi-(float)(sizey/2)+(float)5.e-1f;
 
   // angle extracted from cylindric symmetry
   phi=atanp(x,y);
@@ -422,19 +422,19 @@
 
       if (ni<IdLast[bi])
       {
-        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.)+Trajectories[bi*TRACKPOINTS+ni];
+        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.e0f)+Trajectories[bi*TRACKPOINTS+ni];
       }
       else
       {
         r=Trajectories[bi*TRACKPOINTS+ni];
       }
-	   
+
       if ((r<=re)&&(r>=ri))
       {
         ExitOnImpact=1;
         impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
       }
-	      
+
       HalfLap++;
     } while ((HalfLap<=2)&&(ExitOnImpact==0));
 
@@ -452,7 +452,7 @@
                      float ExternalRadius,float Angle,
                      int Line)
 {
-   // Integer Impact Parameter ID 
+   // Integer Impact Parameter ID
    int bi=get_global_id(0);
    // Integer points on circle
    int i=get_global_id(1);
@@ -474,10 +474,10 @@
 
   float bmx,db,b,h;
   float phi,phd;
-  float zp=0,fp=0;
+  float zp=0.e0f,fp=0.e0f;
 
   // Autosize for image
-  bmx=1.25*re;
+  bmx=1.25e0f*re;
 
   // Angular step of integration
   h=4.e0f*PI/(float)TRACKPOINTS;
@@ -486,7 +486,7 @@
   b=(float)bi/(float)bmaxi*bmx;
   db=bmx/(2.e0f*(float)bmaxi);
 
-  phi=2.*PI/(float)imx*(float)i;
+  phi=2.e0f*PI/(float)imx*(float)i;
   phd=atanp(cos(phi)*sin(tho),cos(tho));
   int yi=(int)((float)bi*sin(phi))+bmaxi;
   int xi=(int)((float)bi*cos(phi))+bmaxi;
@@ -502,24 +502,24 @@
 
      if (ni<IdLast[bi])
      {
-        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.)+Trajectories[bi*TRACKPOINTS+ni];
+        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.e0f)+Trajectories[bi*TRACKPOINTS+ni];
      }
      else
      {
-	r=Trajectories[bi*TRACKPOINTS+ni];
+        r=Trajectories[bi*TRACKPOINTS+ni];
      }
-	   
+
      if ((r<=re)&&(r>=ri))
      {
-	ExitOnImpact=1;
-	impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
+        ExitOnImpact=1;
+        impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
      }
-	      
+
      HalfLap++;
   } while ((HalfLap<=2)&&(ExitOnImpact==0));
-  
+
   zImage[yi+2*bmaxi*xi]=zp;
-  fImage[yi+2*bmaxi*xi]=fp;  
+  fImage[yi+2*bmaxi*xi]=fp;
 
   barrier(CLK_GLOBAL_MEM_FENCE);
 
@@ -530,7 +530,7 @@
                          float ExternalRadius,float Angle,
                          int Line)
 {
-  // Integer Impact Parameter ID 
+  // Integer Impact Parameter ID
   int bi=get_global_id(0);
   // Integer Impact Parameter Size (half of image)
   int bmaxi=get_global_size(0);
@@ -540,14 +540,14 @@
   float m,rs,re;
 
   m=Mass;
-  rs=2.*m;
+  rs=2.e0f*m;
   re=ExternalRadius;
 
   float bmx,b,h;
   int nh;
 
   // Autosize for image
-  bmx=1.25*re;
+  bmx=1.25e0f*re;
 
   // Angular step of integration
   h=4.e0f*PI/(float)TRACKPOINTS;
@@ -557,14 +557,14 @@
 
   float up,vp,pp,us,vs,ps;
 
-  up=0.;
-  vp=1.;
-      
-  pp=0.;
+  up=0.e0f;
+  vp=1.e0f;
+
+  pp=0.e0f;
   nh=0;
 
   rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-  
+
   // b versus us
   float bvus=fabs(b/us);
   float bvus0=bvus;
@@ -585,7 +585,7 @@
   IdLast[bi]=nh;
 
   barrier(CLK_GLOBAL_MEM_FENCE);
- 
+
 }
 
 __kernel void EachCircle(__global float *zImage,__global float *fImage,
@@ -593,18 +593,18 @@
                          float ExternalRadius,float Angle,
                          int Line)
 {
-  // Integer Impact Parameter ID 
-  int bi=get_global_id(0);
-  // Integer Impact Parameter Size (half of image)
-  int bmaxi=get_global_size(0);
+   // Integer Impact Parameter ID
+   uint bi=(uint)get_global_id(0);
+   // Integer Impact Parameter Size (half of image)
+   uint bmaxi=(uint)get_global_size(0);
 
-  float Trajectory[TRACKPOINTS];
+  private float Trajectory[TRACKPOINTS];
 
   float m,rs,ri,re,tho;
   int raie,q;
 
   m=Mass;
-  rs=2.*m;
+  rs=2.e0f*m;
   ri=InternalRadius;
   re=ExternalRadius;
   tho=Angle;
@@ -612,8 +612,9 @@
   raie=Line;
 
   float bmx,db,b,h;
-  int nh;
+  uint nh;
 
+
   // Autosize for image
   bmx=1.25e0f*re;
 
@@ -626,14 +627,14 @@
 
   float up,vp,pp,us,vs,ps;
 
-  up=0.;
-  vp=1.;
-      
-  pp=0.;
+  up=0.e0f;
+  vp=1.e0f;
+
+  pp=0.e0f;
   nh=0;
 
   rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-  
+
   // b versus us
   float bvus=fabs(b/us);
   float bvus0=bvus;
@@ -646,26 +647,28 @@
      up=us;
      vp=vs;
      rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-     bvus=fabs(b/us);
+     bvus=(float)fabs(b/us);
      Trajectory[nh]=bvus;
 
   } while ((bvus>=rs)&&(bvus<=bvus0));
 
+
   for (uint i=(uint)nh+1;i<TRACKPOINTS;i++) {
      Trajectory[i]=0.e0f;
   }
 
-  int imx=(int)(16*bi);
 
-  for (int i=0;i<imx;i++)
+  uint imx=(uint)(16*bi);
+
+  for (uint i=0;i<imx;i++)
   {
-     float zp=0.,fp=0.;
-     float phi=2.*PI/(float)imx*(float)i;
+     float zp=0.e0f,fp=0.e0f;
+     float phi=2.e0f*PI/(float)imx*(float)i;
      float phd=atanp(cos(phi)*sin(tho),cos(tho));
      uint yi=(uint)((float)bi*sin(phi)+bmaxi);
      uint xi=(uint)((float)bi*cos(phi)+bmaxi);
 
-     int HalfLap=0,ExitOnImpact=0,ni;
+     uint HalfLap=0,ExitOnImpact=0,ni;
      float php,nr,r;
 
      do
@@ -676,30 +679,30 @@
 
         if (ni<nh)
         {
-           r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni];
+           r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.e0f)+Trajectory[ni];
         }
         else
         {
            r=Trajectory[ni];
         }
-	   
+
         if ((r<=re)&&(r>=ri))
         {
            ExitOnImpact=1;
            impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
         }
-	      
+
         HalfLap++;
 
      } while ((HalfLap<=2)&&(ExitOnImpact==0));
- 
+
      zImage[yi+2*bmaxi*xi]=zp;
      fImage[yi+2*bmaxi*xi]=fp;
 
   }
 
   barrier(CLK_GLOBAL_MEM_FENCE);
- 
+
 }
 
 __kernel void Original(__global float *zImage,__global float *fImage,
@@ -718,7 +721,7 @@
    int raie,q;
 
    m=Mass;
-   rs=2.*m;
+   rs=2.e0f*m;
    ri=InternalRadius;
    re=ExternalRadius;
    tho=Angle;
@@ -726,7 +729,7 @@
    raie=Line;
 
    float bmx,db,b,h;
-   int nh;
+   uint nh;
 
    // Autosize for image
    bmx=1.25e0f*re;
@@ -734,7 +737,7 @@
    // Angular step of integration
    h=4.e0f*PI/(float)TRACKPOINTS;
 
-   // Integer Impact Parameter ID 
+   // Integer Impact Parameter ID
    for (int bi=0;bi<bmaxi;bi++)
    {
       // impact parameter
@@ -743,14 +746,14 @@
 
       float up,vp,pp,us,vs,ps;
 
-      up=0.;
-      vp=1.;
-      
-      pp=0.;
+      up=0.e0f;
+      vp=1.e0f;
+
+      pp=0.e0f;
       nh=0;
 
       rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-  
+
       // b versus us
       float bvus=fabs(b/us);
       float bvus0=bvus;
@@ -776,13 +779,13 @@
 
       for (int i=0;i<imx;i++)
       {
-         float zp=0,fp=0;
+         float zp=0.e0f,fp=0.e0f;
          float phi=2.e0f*PI/(float)imx*(float)i;
          float phd=atanp(cos(phi)*sin(tho),cos(tho));
          uint yi=(uint)((float)bi*sin(phi)+bmaxi);
          uint xi=(uint)((float)bi*cos(phi)+bmaxi);
 
-         int HalfLap=0,ExitOnImpact=0,ni;
+         uint HalfLap=0,ExitOnImpact=0,ni;
          float php,nr,r;
 
          do
@@ -793,32 +796,32 @@
 
             if (ni<nh)
             {
-               r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni];
+               r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.e0f)+Trajectory[ni];
             }
             else
             {
                r=Trajectory[ni];
             }
-	   
+
             if ((r<=re)&&(r>=ri))
             {
                ExitOnImpact=1;
                impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
             }
- 	      
+
             HalfLap++;
- 
+
          } while ((HalfLap<=2)&&(ExitOnImpact==0));
-  
+
          zImage[yi+2*bmaxi*xi]=zp;
          fImage[yi+2*bmaxi*xi]=fp;
-  
+
       }
 
-   } 
-   
+   }
+
    barrier(CLK_GLOBAL_MEM_FENCE);
- 
+
 }
 """
 
@@ -874,7 +877,7 @@
 #endif
 
 __device__ void calcul(float *us,float *vs,float up,float vp,
-	               float h,float m,float b)
+                       float h,float m,float b)
 {
   float c0,c1,c2,c3,d0,d1,d2,d3;
 
@@ -892,21 +895,21 @@
 }
 
 __device__ void rungekutta(float *ps,float *us,float *vs,
-	                   float pp,float up,float vp,
-		           float h,float m,float b)
+                           float pp,float up,float vp,
+                           float h,float m,float b)
 {
   calcul(us,vs,up,vp,h,m,b);
   *ps=pp+h;
 }
 
 __device__ float decalage_spectral(float r,float b,float phi,
-			           float tho,float m)
+                                   float tho,float m)
 {
   return (sqrt(1-3*m/r)/(1+sqrt(m/pow(r,3))*b*sin(tho)*sin(phi)));
 }
 
 __device__ float spectre(float rf,int q,float b,float db,
-	                 float h,float r,float m,float bss)
+                         float h,float r,float m,float bss)
 {
   float flx;
 
@@ -916,9 +919,9 @@
 }
 
 __device__ float spectre_cn(float rf32,float b32,float db32,
-		            float h32,float r32,float m32,float bss32)
+                            float h32,float r32,float m32,float bss32)
 {
-  
+
 #define MYFLOAT float
 
   MYFLOAT rf=(MYFLOAT)(rf32);
@@ -958,26 +961,26 @@
       nu_rec=nu_em*rf;
       posfreq=(int)(nu_rec*(MYFLOAT)nbr/bss);
       if ((posfreq>0)&&(posfreq<nbr))
-	{
+        {
           // Initial version
-	  // flux_int=2.*planck/c2*pow(nu_em,3)/(exp(planck*nu_em/(k*temp_em))-1.);
+          // flux_int=2.*planck/c2*pow(nu_em,3)/(exp(planck*nu_em/(k*temp_em))-1.);
           // Version with log used
-	  //flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.);
-	  // flux_int*=pow(rf,3)*b*db*h;
-	  //flux_int*=exp(3.*log(rf))*b*db*h;
-	  flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.)*exp(3.*log(rf))*b*db*h;
+          //flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.);
+          // flux_int*=pow(rf,3)*b*db*h;
+          //flux_int*=exp(3.*log(rf))*b*db*h;
+          flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.)*exp(3.*log(rf))*b*db*h;
 
-	  flx+=flux_int;
-	}
+          flx+=flux_int;
+        }
     }
 
   return((float)(flx));
 }
 
 __device__ void impact(float phi,float r,float b,float tho,float m,
-	               float *zp,float *fp,
-	               int q,float db,
-	               float h,int raie)
+                       float *zp,float *fp,
+                       int q,float db,
+                       float h,int raie)
 {
   float flx,rf,bss;
 
@@ -989,11 +992,11 @@
       flx=spectre_cn(rf,b,db,h,r,m,bss);
     }
   else
-    { 
+    {
       bss=2.;
       flx=spectre(rf,q,b,db,h,r,m,bss);
     }
-  
+
   *zp=1./rf;
   *fp=flx;
 
@@ -1051,12 +1054,11 @@
 
   up=0.;
   vp=1.;
-      
   pp=0.;
   nh=0;
 
   rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-    
+
   rps=fabs(b/us);
   rp0=rps;
 
@@ -1069,7 +1071,7 @@
      up=us;
      vp=vs;
      rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-     rpp=rps;	  
+     rpp=rps;
      rps=fabs(b/us);
      ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rps>ri)&&(rps<re)?1:0;          
 
@@ -1080,8 +1082,8 @@
   }
   else
   {
-     zp=0.;
-     fp=0.;
+     zp=0.e0f;
+     fp=0.e0f;
   }
 
   __syncthreads();
@@ -1108,7 +1110,7 @@
   int q,raie;
 
   m=Mass;
-  rs=2.*m;
+  rs=2.e0f*m;
   ri=InternalRadius;
   re=ExternalRadius;
   tho=Angle;
@@ -1120,13 +1122,13 @@
   int nh;
   float zp=0,fp=0;
   // Autosize for image, 25% greater than external radius
-  bmx=1.25*re;
+  bmx=1.25e0f*re;
 
   // Angular step of integration
   h=4.e0f*PI/(float)TRACKPOINTS;
 
   // Step of Impact Parameter
-  db=bmx/(2.e0*(float)ImpactParameter);
+  db=bmx/(2.e0f*(float)ImpactParameter);
 
   // set origin as center of image
   float x=(float)xi-(float)(sizex/2)+(float)5e-1f;
@@ -1153,19 +1155,19 @@
 
       if (ni<IdLast[bi])
       {
-        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.)+Trajectories[bi*TRACKPOINTS+ni];
+        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.e0f)+Trajectories[bi*TRACKPOINTS+ni];
       }
       else
       {
         r=Trajectories[bi*TRACKPOINTS+ni];
       }
-	   
+
       if ((r<=re)&&(r>=ri))
       {
         ExitOnImpact=1;
         impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
       }
-	      
+
       HalfLap++;
     } while ((HalfLap<=2)&&(ExitOnImpact==0));
 
@@ -1181,7 +1183,7 @@
                        float ExternalRadius,float Angle,
                        int Line)
 {
-   // Integer Impact Parameter ID 
+   // Integer Impact Parameter ID
    int bi=blockIdx.x*blockDim.x+threadIdx.x;
    // Integer points on circle
    int i=blockIdx.y*blockDim.y+threadIdx.y;
@@ -1196,7 +1198,7 @@
   int q,raie;
 
   m=Mass;
-  rs=2.*m;
+  rs=2.e0f*m;
   ri=InternalRadius;
   re=ExternalRadius;
   tho=Angle;
@@ -1208,16 +1210,16 @@
   float zp=0,fp=0;
 
   // Autosize for image
-  bmx=1.25*re;
+  bmx=1.25e0f*re;
 
   // Angular step of integration
   h=4.e0f*PI/(float)TRACKPOINTS;
 
   // impact parameter
   b=(float)bi/(float)bmaxi*bmx;
-  db=bmx/(2.e0*(float)bmaxi);
+  db=bmx/(2.e0f*(float)bmaxi);
 
-  phi=2.*PI/(float)imx*(float)i;
+  phi=2.e0f*PI/(float)imx*(float)i;
   phd=atanp(cos(phi)*sin(tho),cos(tho));
   int yi=(int)((float)bi*sin(phi))+bmaxi;
   int xi=(int)((float)bi*cos(phi))+bmaxi;
@@ -1233,24 +1235,24 @@
 
      if (ni<IdLast[bi])
      {
-        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.)+Trajectories[bi*TRACKPOINTS+ni];
+        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.e0f)+Trajectories[bi*TRACKPOINTS+ni];
      }
      else
      {
-	r=Trajectories[bi*TRACKPOINTS+ni];
+        r=Trajectories[bi*TRACKPOINTS+ni];
      }
-	   
+
      if ((r<=re)&&(r>=ri))
      {
-	ExitOnImpact=1;
-	impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
+        ExitOnImpact=1;
+        impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
      }
-	      
+
      HalfLap++;
   } while ((HalfLap<=2)&&(ExitOnImpact==0));
-  
+
   zImage[yi+2*bmaxi*xi]=zp;
-  fImage[yi+2*bmaxi*xi]=fp;  
+  fImage[yi+2*bmaxi*xi]=fp;
 
 }
 
@@ -1259,7 +1261,7 @@
                            float ExternalRadius,float Angle,
                            int Line)
 {
-  // Integer Impact Parameter ID 
+  // Integer Impact Parameter ID
   int bi=blockIdx.x*blockDim.x+threadIdx.x;
   // Integer Impact Parameter Size (half of image)
   int bmaxi=gridDim.x*blockDim.x;
@@ -1270,7 +1272,7 @@
   int raie,q;
 
   m=Mass;
-  rs=2.*m;
+  rs=2.e0f*m;
   ri=InternalRadius;
   re=ExternalRadius;
   tho=Angle;
@@ -1283,7 +1285,7 @@
   float zp,fp;
 
   // Autosize for image
-  bmx=1.25*re;
+  bmx=1.25e0f*re;
 
   // Angular step of integration
   h=4.e0f*PI/(float)TRACKPOINTS;
@@ -1293,14 +1295,13 @@
 
   float up,vp,pp,us,vs,ps;
 
-  up=0.;
-  vp=1.;
-      
-  pp=0.;
+  up=0.e0f;
+  vp=1.e0f;
+  pp=0.e0f;
   nh=0;
 
   rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-  
+
   // b versus us
   float bvus=fabs(b/us);
   float bvus0=bvus;
@@ -1327,7 +1328,7 @@
                            float ExternalRadius,float Angle,
                            int Line)
 {
-  // Integer Impact Parameter ID 
+  // Integer Impact Parameter ID
   int bi=blockIdx.x*blockDim.x+threadIdx.x;
 
   // Integer Impact Parameter Size (half of image)
@@ -1365,12 +1366,11 @@
 
   up=0.;
   vp=1.;
-      
   pp=0.;
   nh=0;
 
   rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-  
+
   // b versus us
   float bvus=fabs(b/us);
   float bvus0=bvus;
@@ -1415,13 +1415,13 @@
         {
            r=Trajectory[ni];
         }
-	   
+
         if ((r<=re)&&(r>=ri))
         {
            ExitOnImpact=1;
            impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
         }
-	      
+
         HalfLap++;
 
      } while ((HalfLap<=2)&&(ExitOnImpact==0));
@@ -1432,7 +1432,7 @@
    fImage[yi+2*bmaxi*xi]=fp;
 
   }
- 
+
 }
 
 __global__ void Original(float *zImage,float *fImage,
@@ -1451,7 +1451,7 @@
    int raie,q;
 
    m=Mass;
-   rs=2.*m;
+   rs=2.e0f*m;
    ri=InternalRadius;
    re=ExternalRadius;
    tho=Angle;
@@ -1467,7 +1467,7 @@
    // Angular step of integration
    h=4.e0f*PI/(float)TRACKPOINTS;
 
-   // Integer Impact Parameter ID 
+   // Integer Impact Parameter ID
    for (int bi=0;bi<bmaxi;bi++)
    {
       // impact parameter
@@ -1478,12 +1478,11 @@
 
       up=0.;
       vp=1.;
-      
       pp=0.;
       nh=0;
 
       rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
-  
+
       // b versus us
       float bvus=fabs(b/us);
       float bvus0=bvus;
@@ -1532,33 +1531,32 @@
             {
                r=Trajectory[ni];
             }
-	   
+
             if ((r<=re)&&(r>=ri))
             {
                ExitOnImpact=1;
                impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie);
             }
- 	      
+
             HalfLap++;
- 
+
          } while ((HalfLap<=2)&&(ExitOnImpact==0));
-  
+
          zImage[yi+2*bmaxi*xi]=zp;
          fImage[yi+2*bmaxi*xi]=fp;
-  
+
       }
 
-   } 
-   
+   }
+
 }
 """
     return(BlobCUDA)
-    
+
 # def ImageOutput(sigma,prefix):
 #     from PIL import Image
 #     Max=sigma.max()
 #     Min=sigma.min()
-    
 #     # Normalize value as 8bits Integer
 #     SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
 #     image = Image.fromarray(SigmaInt)
@@ -1566,12 +1564,13 @@
 
 def ImageOutput(sigma,prefix,Colors):
     import matplotlib.pyplot as plt
-    start_time=time.time() 
+    start_time=time.time()
     if Colors == 'Red2Yellow':
         plt.imsave("%s.png" % prefix, sigma, cmap='afmhot')
     else:
         plt.imsave("%s.png" % prefix, sigma, cmap='Greys_r')
     save_time = time.time()-start_time
+    print("Save image as %s.png file" % prefix)
     print("Save Time : %f" % save_time)
 
 def BlackHoleCL(zImage,fImage,InputCL):
@@ -1610,6 +1609,7 @@
     for platform in cl.get_platforms():
         for device in platform.get_devices():
             if Id==Device:
+                PF4XPU=platform.name
                 XPU=device
                 print("CPU/GPU selected: ",device.name.lstrip())
                 HasXPU=True
@@ -1617,8 +1617,8 @@
 
     if HasXPU==False:
         print("No XPU #%i found in all of %i devices, sorry..." % (Device,Id-1))
-        sys.exit()		
-	
+        sys.exit()
+
     ctx = cl.Context([XPU])
     queue = cl.CommandQueue(ctx,
 			    properties=cl.command_queue_properties.PROFILING_ENABLE)
@@ -1628,9 +1628,9 @@
 
     BuildOptions="-DPHYSICS=%i -DSETTRACKPOINTS=%i " % (PhysicsList[Physics],InputCL['TrackPoints'])
 
-    print('My Platform is ',platform.name)
+    print('My Platform is ',PF4XPU)
     
-    if 'Intel' in platform.name or 'Experimental' in platform.name or 'Clover' in platform.name or 'Portable' in platform.name :
+    if 'Intel' in PF4XPU or 'Experimental' in PF4XPU or 'Clover' in PF4XPU or 'Portable' in PF4XPU :
         print('No extra options for Intel and Clover!')
     else:
         BuildOptions = BuildOptions+" -cl-mad-enable"
@@ -1649,7 +1649,7 @@
     zImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=zImage)
     fImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=fImage)
     
-    start_time=time.time() 
+    start_time=time.time()
 
     if Method=='EachPixel':
         CLLaunch=BlackHoleCL.EachPixel(queue,(zImage.shape[0],zImage.shape[1]),
@@ -1810,7 +1810,7 @@
     IdLastCU = cuda.mem_alloc(IdLast.size*IdLast.dtype.itemsize)
     cuda.memcpy_htod(IdLastCU, IdLast)
 
-    start_time=time.time() 
+    start_time=time.time()
 
     if Method=='EachPixel':
         EachPixelCU(zImageCU,fImageCU,
