Index: TrouNoir.py
===================================================================
--- TrouNoir.py	(revision 216)
+++ TrouNoir.py	(revision 217)
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 #
 # TrouNoir model using PyOpenCL 
 #
@@ -45,6 +45,8 @@
 #define EINSTEIN 0
 #define NEWTON 1
 
+#define TRACKPOINTS 2048
+
 float atanp(float x,float y)
 {
   float angle;
@@ -213,8 +215,8 @@
 
    // Perform trajectory for each pixel, exit on hit
 
-  float m,rs,ri,re,tho;
-  int q,raie;
+  private float m,rs,ri,re,tho;
+  private int q,raie;
 
   m=Mass;
   rs=2.*m;
@@ -224,17 +226,17 @@
   q=-2;
   raie=Line;
 
-  float d,bmx,db,b,h;
-  float rp[2048];
-  float phi,thi,phd,php,nr,r;
-  int nh;
-  float zp,fp;
+  private float d,bmx,db,b,h;
+  private float rp[TRACKPOINTS];
+  private float phi,thi,phd,php,nr,r;
+  private int nh;
+  private float zp,fp;
 
   // Autosize for image
   bmx=1.25*re;
   b=0.;
 
-  h=4.e0f*PI/(float)2048;
+  h=4.e0f*PI/(float)TRACKPOINTS;
 
   // set origin as center of image
   float x=(float)xi-(float)(sizex/2)+(float)5e-1f;
@@ -533,6 +535,501 @@
 }
 """
 
+def KernelCodeCuda():
+    BlobCUDA= """
+
+#define PI (float)3.14159265359
+#define nbr 256
+
+#define EINSTEIN 0
+#define NEWTON 1
+
+#define TRACKPOINTS 2048
+
+__device__ float nothing(float x)
+{
+  return(x);
+}
+
+__device__ float atanp(float x,float y)
+{
+  float angle;
+
+  angle=atan2(y,x);
+
+  if (angle<0.e0f)
+    {
+      angle+=(float)2.e0f*PI;
+    }
+
+  return(angle);
+}
+
+__device__ float f(float v)
+{
+  return(v);
+}
+
+#if PHYSICS == NEWTON
+__device__ float g(float u,float m,float b)
+{
+  return (-u);
+}
+#else
+__device__ float g(float u,float m,float b)
+{
+  return (3.e0f*m/b*pow(u,2)-u);
+}
+#endif
+
+__device__ void calcul(float *us,float *vs,float up,float vp,
+	               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.);
+  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);
+  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.;
+}
+
+__device__ 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;
+}
+
+__device__ 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)));
+}
+
+__device__ float spectre(float rf,int q,float b,float db,
+	                 float h,float r,float m,float bss)
+{
+  float flx;
+
+  flx=exp(q*log(r/m))*pow(rf,4)*b*db*h;
+  return(flx);
+}
+
+__device__ float spectre_cn(float rf32,float b32,float db32,
+		            float h32,float r32,float m32,float bss32)
+{
+  
+#define MYFLOAT float
+
+  MYFLOAT rf=(MYFLOAT)(rf32);
+  MYFLOAT b=(MYFLOAT)(b32);
+  MYFLOAT db=(MYFLOAT)(db32);
+  MYFLOAT h=(MYFLOAT)(h32);
+  MYFLOAT r=(MYFLOAT)(r32);
+  MYFLOAT m=(MYFLOAT)(m32);
+  MYFLOAT bss=(MYFLOAT)(bss32);
+
+  MYFLOAT flx;
+  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 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;
+
+  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 ));
+
+  temp_em=temp*sqrt(m)*exp(0.25*log(m_point)-0.75*log(r)-0.125*log(v)+0.25*log(fabs(qu)));
+
+  flux_int=0.;
+  flx=0.;
+
+  for (fi=0;fi<nbr;fi++)
+    {
+      nu_em=bss*(MYFLOAT)fi/(MYFLOAT)nbr;
+      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.);
+          // 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;
+
+	  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 flx,rf,bss;
+
+  rf=decalage_spectral(r,b,phi,tho,m);
+
+  if (raie==0)
+    {
+      bss=1.e19;
+      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;
+
+}
+
+__global__ void EachPixel(float *zImage,float *fImage,
+                          float Mass,float InternalRadius,
+                          float ExternalRadius,float Angle,
+                          int Line)
+{
+   uint xi=(uint)blockIdx.x;
+   uint yi=(uint)blockIdx.y;
+   uint sizex=(uint)gridDim.x*blockDim.x;
+   uint sizey=(uint)gridDim.y*blockDim.y;
+
+   // Perform trajectory for each pixel, exit on hit
+
+  float m,rs,ri,re,tho;
+  int q,raie;
+
+  m=Mass;
+  rs=2.*m;
+  ri=InternalRadius;
+  re=ExternalRadius;
+  tho=Angle;
+  q=-2;
+  raie=Line;
+
+  float d,bmx,db,b,h;
+  float rp[TRACKPOINTS];
+  float phi,thi,phd,php,nr,r;
+  int nh;
+  float zp,fp;
+
+  // Autosize for image
+  bmx=1.25*re;
+  b=0.;
+
+  h=4.e0f*PI/(float)TRACKPOINTS;
+
+  // 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;
+  // angle extracted from cylindric symmetry
+  phi=atanp(x,y);
+  phd=atanp(cos(phi)*sin(tho),cos(tho));
+
+  float up,vp,pp,us,vs,ps;
+
+  // impact parameter
+  b=sqrt(x*x+y*y)*(float)2.e0f/(float)sizex*bmx;
+  // step of impact parameter;
+//  db=bmx/(float)(sizex/2);
+  db=bmx/(float)(sizex);
+
+  up=0.;
+  vp=1.;
+      
+  pp=0.;
+  nh=0;
+
+  rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
+    
+  rp[nh]=fabs(b/us);
+
+  int ExitOnImpact=0;
+
+  do
+  {
+     nh++;
+     pp=ps;
+     up=us;
+     vp=vs;
+     rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);	  
+     rp[nh]=fabs(b/us);
+     ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rp[nh]>ri)&&(rp[nh]<re)?1:0;          
+
+  } while ((rp[nh]>=rs)&&(rp[nh]<=rp[0])&&(ExitOnImpact==0));
+
+  if (ExitOnImpact==1) {
+     impact(phi,rp[nh-1],b,tho,m,&zp,&fp,q,db,h,raie);
+  }
+  else
+  {
+     zp=0.;
+     fp=0.;
+  }
+
+  zImage[yi+sizex*xi]=(float)zp;
+  fImage[yi+sizex*xi]=(float)fp;
+}
+
+__global__ void Pixel(float *zImage,float *fImage,
+                      float *Trajectories,int *IdLast,
+                      uint ImpactParameter,uint TrackPoints,
+                      float Mass,float InternalRadius,
+                      float ExternalRadius,float Angle,
+                      int Line)
+{
+  uint xi=(uint)blockIdx.x;
+  uint yi=(uint)blockIdx.y;
+  uint sizex=(uint)gridDim.x*blockDim.x;
+  uint sizey=(uint)gridDim.y*blockDim.y;
+
+  // Perform trajectory for each pixel
+
+  float m,rs,ri,re,tho;
+  int q,raie;
+
+  m=Mass;
+  rs=2.*m;
+  ri=InternalRadius;
+  re=ExternalRadius;
+  tho=Angle;
+  q=-2;
+  raie=Line;
+
+  float d,bmx,db,b,h;
+  float phi,thi,phd,php,nr,r;
+  int nh;
+  float zp=0,fp=0;
+  // Autosize for image, 25% greater than external radius
+  bmx=1.25*re;
+
+  // Angular step of integration
+  h=4.e0f*PI/(float)TrackPoints;
+
+  // Step of Impact Parameter
+  db=bmx/(2.e0*(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;
+  // angle extracted from cylindric symmetry
+  phi=atanp(x,y);
+  phd=atanp(cos(phi)*sin(tho),cos(tho));
+
+  // Real Impact Parameter
+  b=sqrt(x*x+y*y)*bmx/(float)ImpactParameter;
+
+  // Integer Impact Parameter
+  uint bi=(uint)sqrt(x*x+y*y);
+
+  int HalfLap=0,ExitOnImpact=0,ni;
+
+  if (bi<ImpactParameter)
+  {
+    do
+    {
+      php=phd+(float)HalfLap*PI;
+      nr=php/h;
+      ni=(int)nr;
+
+      if (ni<IdLast[bi])
+      {
+        r=(Trajectories[bi*TrackPoints+ni+1]-Trajectories[bi*TrackPoints+ni])*(nr-ni*1.)+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));
+
+  }
+
+  zImage[yi+sizex*xi]=zp;
+  fImage[yi+sizex*xi]=fp;
+}
+
+__global__ void Circle(float *Trajectories,int *IdLast,
+                       float *zImage,float *fImage,
+                       int TrackPoints,
+                       float Mass,float InternalRadius,
+                       float ExternalRadius,float Angle,
+                       int Line)
+{
+   // Integer Impact Parameter ID 
+   int bi=blockIdx.x;
+   // Integer points on circle
+   int i=blockIdx.y;
+   // Integer Impact Parameter Size (half of image)
+   int bmaxi=gridDim.x*blockDim.x;
+   // Integer Points on circle
+   int imx=gridDim.y*blockDim.y;
+
+   // Perform trajectory for each pixel
+
+  float m,rs,ri,re,tho;
+  int q,raie;
+
+  m=Mass;
+  rs=2.*m;
+  ri=InternalRadius;
+  re=ExternalRadius;
+  tho=Angle;
+  raie=Line;
+
+  float bmx,db,b,h;
+  float phi,thi,phd;
+  int nh;
+  float zp=0,fp=0;
+
+  // Autosize for image
+  bmx=1.25*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);
+
+  phi=2.*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;
+
+  int HalfLap=0,ExitOnImpact=0,ni;
+  float php,nr,r;
+
+  do
+  {
+     php=phd+(float)HalfLap*PI;
+     nr=php/h;
+     ni=(int)nr;
+
+     if (ni<IdLast[bi])
+     {
+        r=(Trajectories[bi*TrackPoints+ni+1]-Trajectories[bi*TrackPoints+ni])*(nr-ni*1.)+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));
+  
+  zImage[yi+2*bmaxi*xi]=zp;
+  fImage[yi+2*bmaxi*xi]=fp;  
+
+}
+
+__global__ void Trajectory(float *Trajectories,
+                           int *IdLast,int TrackPoints,
+                           float Mass,float InternalRadius,
+                           float ExternalRadius,float Angle,
+                           int Line)
+{
+  // Integer Impact Parameter ID 
+  int bi=blockIdx.x;
+  // Integer Impact Parameter Size (half of image)
+  int bmaxi=gridDim.x*blockDim.x;
+
+  // Perform trajectory for each pixel
+
+  float m,rs,ri,re,tho;
+  int raie,q;
+
+  m=Mass;
+  rs=2.*m;
+  ri=InternalRadius;
+  re=ExternalRadius;
+  tho=Angle;
+  q=-2;
+  raie=Line;
+
+  float d,bmx,db,b,h;
+  float phi,thi,phd,php,nr,r;
+  int nh;
+  float zp,fp;
+
+  // Autosize for image
+  bmx=1.25*re;
+
+  // Angular step of integration
+  h=4.e0f*PI/(float)TrackPoints;
+
+  // impact parameter
+  b=(float)bi/(float)bmaxi*bmx;
+
+  float up,vp,pp,us,vs,ps;
+
+  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;
+  Trajectories[bi*TrackPoints+nh]=bvus;
+
+  do
+  {
+     nh++;
+     pp=ps;
+     up=us;
+     vp=vs;
+     rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
+     bvus=fabs(b/us);
+     Trajectories[bi*TrackPoints+nh]=bvus;
+
+  } while ((bvus>=rs)&&(bvus<=bvus0));
+
+  IdLast[bi]=nh;
+
+}
+"""
+    return(BlobCUDA)
+    
 # def ImageOutput(sigma,prefix):
 #     from PIL import Image
 #     Max=sigma.max()
@@ -577,9 +1074,8 @@
         # Spectrum is Monochromatic Line one
         Line=1
 
-    Trajectories=numpy.zeros((InputCL['Size']/2,InputCL['TrackPoints']),
-                              dtype=numpy.float32)
-    IdLast=numpy.zeros(InputCL['Size']/2,dtype=numpy.int32)
+    Trajectories=numpy.zeros((int(InputCL['Size']/2),InputCL['TrackPoints']),dtype=numpy.float32)
+    IdLast=numpy.zeros(int(InputCL['Size']/2),dtype=numpy.int32)
 
     # Je detecte un peripherique GPU dans la liste des peripheriques
     Id=0
@@ -588,12 +1084,12 @@
         for device in platform.get_devices():
             if Id==Device:
                 XPU=device
-                print "CPU/GPU selected: ",device.name.lstrip()
+                print("CPU/GPU selected: ",device.name.lstrip())
                 HasXPU=True
             Id+=1
 
     if HasXPU==False:
-        print "No XPU #%i found in all of %i devices, sorry..." % (Device,Id-1)
+        print("No XPU #%i found in all of %i devices, sorry..." % (Device,Id-1))
         sys.exit()		
 	
     ctx = cl.Context([XPU])
@@ -658,7 +1154,7 @@
                                         numpy.int32(Line))
         
         CLLaunch=BlackHoleCL.Pixel(queue,(zImage.shape[0],
-                                                  zImage.shape[1]),None,
+                                          zImage.shape[1]),None,
                                    zImageCL,fImageCL,TrajectoriesCL,IdLastCL,
                                    numpy.uint32(Trajectories.shape[0]),
                                    numpy.uint32(Trajectories.shape[1]),
@@ -689,6 +1185,137 @@
 
     return(elapsed)
 
+def BlackHoleCUDA(zImage,fImage,InputCL):
+
+    kernel_params = {}
+
+    print(InputCL)
+
+    Device=InputCL['Device']
+    GpuStyle=InputCL['GpuStyle']
+    VariableType=InputCL['VariableType']
+    Size=InputCL['Size']
+    Mass=InputCL['Mass']
+    InternalRadius=InputCL['InternalRadius']
+    ExternalRadius=InputCL['ExternalRadius']
+    Angle=InputCL['Angle']
+    Method=InputCL['Method']
+    TrackPoints=InputCL['TrackPoints']
+    Physics=InputCL['Physics']
+
+    PhysicsList=DictionariesAPI()
+    
+    if InputCL['BlackBody']:
+        # Spectrum is Black Body one
+        Line=0
+    else:
+        # Spectrum is Monochromatic Line one
+        Line=1
+
+    Trajectories=numpy.zeros((int(InputCL['Size']/2),InputCL['TrackPoints']),dtype=numpy.float32)
+    IdLast=numpy.zeros(int(InputCL['Size']/2),dtype=numpy.int32)
+
+    try:
+        # For PyCUDA import
+        import pycuda.driver as cuda
+        from pycuda.compiler import SourceModule
+        
+        cuda.init()
+        for Id in range(cuda.Device.count()):
+            if Id==Device:
+                XPU=cuda.Device(Id)
+                print("GPU selected %s" % XPU.name())
+        print
+
+    except ImportError:
+        print("Platform does not seem to support CUDA")
+
+    Context=XPU.make_context()
+
+    try:
+        mod = SourceModule(KernelCodeCuda(),options=['--compiler-options','-DPHYSICS=%i' % (PhysicsList[Physics])])
+        print("Compilation seems to be OK")
+    except:
+        print("Compilation seems to break")
+
+    EachPixelCU=mod.get_function("EachPixel")
+    TrajectoryCU=mod.get_function("Trajectory")
+    PixelCU=mod.get_function("Pixel")
+    CircleCU=mod.get_function("Circle")
+    
+    TrajectoriesCU = cuda.mem_alloc(Trajectories.size*Trajectories.dtype.itemsize)
+    cuda.memcpy_htod(TrajectoriesCU, Trajectories)
+    zImageCU = cuda.mem_alloc(zImage.size*zImage.dtype.itemsize)
+    cuda.memcpy_htod(zImageCU, zImage)
+    fImageCU = cuda.mem_alloc(fImage.size*fImage.dtype.itemsize)
+    cuda.memcpy_htod(zImageCU, fImage)
+    IdLastCU = cuda.mem_alloc(IdLast.size*IdLast.dtype.itemsize)
+    cuda.memcpy_htod(IdLastCU, IdLast)
+
+    start_time=time.time() 
+
+    if Method=='EachPixel':
+        EachPixelCU(zImageCU,fImageCU,
+                    numpy.float32(Mass),
+                    numpy.float32(InternalRadius),
+                    numpy.float32(ExternalRadius),
+                    numpy.float32(Angle),
+                    numpy.int32(Line),
+                    grid=(zImage.shape[0],zImage.shape[1]),block=(1,1,1))
+    elif Method=='TrajectoCircle':
+        TrajectoryCU(TrajectoriesCU,IdLastCU,
+                     numpy.uint32(Trajectories.shape[1]),
+                     numpy.float32(Mass),
+                     numpy.float32(InternalRadius),
+                     numpy.float32(ExternalRadius),
+                     numpy.float32(Angle),
+                     numpy.int32(Line),
+                     grid=(Trajectories.shape[0],1),block=(1,1,1))
+        
+        CircleCU(TrajectoriesCU,IdLastCU,zImageCU,fImageCU,
+                 numpy.uint32(Trajectories.shape[1]),
+                 numpy.float32(Mass),
+                 numpy.float32(InternalRadius),
+                 numpy.float32(ExternalRadius),
+                 numpy.float32(Angle),
+                 numpy.int32(Line),
+                 grid=(Trajectories.shape[0],zImage.shape[0]*4),block=(1,1,1))
+    else:
+        TrajectoryCU(TrajectoriesCU,IdLastCU,
+                     numpy.uint32(Trajectories.shape[1]),
+                     numpy.float32(Mass),
+                     numpy.float32(InternalRadius),
+                     numpy.float32(ExternalRadius),
+                     numpy.float32(Angle),
+                     numpy.int32(Line),
+                     grid=(Trajectories.shape[0],1),block=(1,1,1))
+        
+        PixelCU(zImageCU,fImageCU,TrajectoriesCU,IdLastCU,
+                numpy.uint32(Trajectories.shape[0]),
+                numpy.uint32(Trajectories.shape[1]),
+                numpy.float32(Mass),
+                numpy.float32(InternalRadius),
+                numpy.float32(ExternalRadius),
+                numpy.float32(Angle),
+                numpy.int32(Line),
+                grid=(zImage.shape[0],zImage.shape[1],1),block=(1,1,1))
+    
+
+    elapsed = time.time()-start_time
+    print("\nElapsed Time : %f" % elapsed)
+
+    cuda.memcpy_dtoh(zImage,zImageCU)
+    cuda.memcpy_dtoh(fImage,fImageCU)
+
+    zMaxPosition=numpy.where(zImage[:,:]==zImage.max())
+    fMaxPosition=numpy.where(fImage[:,:]==fImage.max())
+    print("Z max @(%i,%i) : %f" % (zMaxPosition[0][0],zMaxPosition[1][0],zImage.max()))
+    print("Flux max @(%i,%i) : %f\n" % (fMaxPosition[0][0],fMaxPosition[1][0],fImage.max()))
+
+    Context.detach()
+
+    return(elapsed)
+
 if __name__=='__main__':
 	
     GpuStyle = 'OpenCL'
@@ -858,9 +1485,12 @@
     InputCL['Method']=Method
     InputCL['TrackPoints']=TrackPoints
     InputCL['Physics']=Physics
-    
-    duration=BlackHoleCL(zImage,fImage,InputCL)
 
+    if GpuStyle=='OpenCL':
+        duration=BlackHoleCL(zImage,fImage,InputCL)
+    else:
+        duration=BlackHoleCUDA(zImage,fImage,InputCL)
+        
     Hostname=gethostname()
     Date=time.strftime("%Y%m%d_%H%M%S")
     ImageInfo="%s_Device%i_%s_%s" % (Method,Device,Hostname,Date)
