Index: TrouNoir/TrouNoir.py
===================================================================
--- TrouNoir/TrouNoir.py	(revision 223)
+++ TrouNoir/TrouNoir.py	(revision 224)
@@ -46,7 +46,11 @@
 #define EINSTEIN 0
 #define NEWTON 1
 
+#ifdef SETTRACKPOINTS
+#define TRACKPOINTS SETTRACKPOINTS
+#else
 #define TRACKPOINTS 2048
+#endif
 
 float atanp(float x,float y)
 {
@@ -297,7 +301,7 @@
 
 __kernel void Pixel(__global float *zImage,__global float *fImage,
                     __global float *Trajectories,__global int *IdLast,
-                    uint ImpactParameter,uint TrackPoints,
+                    uint ImpactParameter,
                     float Mass,float InternalRadius,
                     float ExternalRadius,float Angle,
                     int Line)
@@ -329,7 +333,7 @@
   bmx=1.25*re;
 
   // Angular step of integration
-  h=4.e0f*PI/(float)TrackPoints;
+  h=4.e0f*PI/(float)TRACKPOINTS;
 
   // Step of Impact Parameter
   db=bmx/(2.e0*(float)ImpactParameter);
@@ -360,11 +364,11 @@
 
       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.)+Trajectories[bi*TRACKPOINTS+ni];
       }
       else
       {
-        r=Trajectories[bi*TrackPoints+ni];
+        r=Trajectories[bi*TRACKPOINTS+ni];
       }
 	   
       if ((r<=re)&&(r>=ri))
@@ -386,7 +390,6 @@
 
 __kernel void Circle(__global float *Trajectories,__global int *IdLast,
                      __global float *zImage,__global float *fImage,
-                     int TrackPoints,
                      float Mass,float InternalRadius,
                      float ExternalRadius,float Angle,
                      int Line)
@@ -421,11 +424,11 @@
   bmx=1.25*re;
 
   // Angular step of integration
-  h=4.e0f*PI/(float)TrackPoints;
+  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;
   phd=atanp(cos(phi)*sin(tho),cos(tho));
@@ -443,11 +446,11 @@
 
      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.)+Trajectories[bi*TRACKPOINTS+ni];
      }
      else
      {
-	r=Trajectories[bi*TrackPoints+ni];
+	r=Trajectories[bi*TRACKPOINTS+ni];
      }
 	   
      if ((r<=re)&&(r>=ri))
@@ -466,8 +469,7 @@
 
 }
 
-__kernel void Trajectory(__global float *Trajectories,
-                         __global int *IdLast,int TrackPoints,
+__kernel void Trajectory(__global float *Trajectories,__global int *IdLast,
                          float Mass,float InternalRadius,
                          float ExternalRadius,float Angle,
                          int Line)
@@ -479,6 +481,71 @@
 
   // Perform trajectory for each pixel
 
+  float m,rs,re;
+
+  m=Mass;
+  rs=2.*m;
+  re=ExternalRadius;
+
+  float bmx,b,h;
+  int nh;
+
+  // 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;
+
+  barrier(CLK_GLOBAL_MEM_FENCE);
+ 
+}
+
+__kernel void Original(__global float *zImage,__global float *fImage,
+                       float Mass,float InternalRadius,
+                       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);
+
+  float Trajectories[2048];
+
+  // Perform trajectory for each pixel
+
   float m,rs,ri,re,tho;
   int raie,q;
 
@@ -490,19 +557,18 @@
   q=-2;
   raie=Line;
 
-  float d,bmx,db,b,h;
-  float phi,thi,phd,php,nr,r;
+  float bmx,db,b,h;
   int nh;
-  float zp,fp;
 
   // Autosize for image
-  bmx=1.25*re;
+  bmx=1.25e0f*re;
 
   // Angular step of integration
-  h=4.e0f*PI/(float)TrackPoints;
+  h=4.e0f*PI/(float)TRACKPOINTS;
 
   // impact parameter
   b=(float)bi/(float)bmaxi*bmx;
+  db=bmx/(2.e0f*(float)bmaxi);
 
   float up,vp,pp,us,vs,ps;
 
@@ -517,7 +583,7 @@
   // b versus us
   float bvus=fabs(b/us);
   float bvus0=bvus;
-  Trajectories[bi*TrackPoints+nh]=bvus;
+  Trajectories[nh]=bvus;
 
   do
   {
@@ -527,12 +593,53 @@
      vp=vs;
      rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
      bvus=fabs(b/us);
-     Trajectories[bi*TrackPoints+nh]=bvus;
+     Trajectories[nh]=bvus;
 
   } while ((bvus>=rs)&&(bvus<=bvus0));
 
-  IdLast[bi]=nh;
+  int imx=(int)(16*bi);
 
+  for (int i=0;i<imx;i++)
+  {
+     float zp=0,fp=0;
+     float phi=2.*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;
+     float php,nr,r;
+
+     do
+     {
+        php=phd+(float)HalfLap*PI;
+        nr=php/h;
+        ni=(int)nr;
+
+        if (ni<nh)
+        {
+           r=(Trajectories[ni+1]-Trajectories[ni])*(nr-ni*1.)+Trajectories[ni];
+        }
+        else
+        {
+           r=Trajectories[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);
  
 }
@@ -547,7 +654,11 @@
 #define EINSTEIN 0
 #define NEWTON 1
 
-#define TRACKPOINTS 2048
+#ifdef SETTRACKPOINTS
+#define TRACKPOINTS SETTRACKPOINTS
+#else
+#define TRACKPOINTS
+#endif
 
 __device__ float nothing(float x)
 {
@@ -804,7 +915,7 @@
 
 __global__ void Pixel(float *zImage,float *fImage,
                       float *Trajectories,int *IdLast,
-                      uint ImpactParameter,uint TrackPoints,
+                      uint ImpactParameter,
                       float Mass,float InternalRadius,
                       float ExternalRadius,float Angle,
                       int Line)
@@ -835,7 +946,7 @@
   bmx=1.25*re;
 
   // Angular step of integration
-  h=4.e0f*PI/(float)TrackPoints;
+  h=4.e0f*PI/(float)TRACKPOINTS;
 
   // Step of Impact Parameter
   db=bmx/(2.e0*(float)ImpactParameter);
@@ -865,11 +976,11 @@
 
       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.)+Trajectories[bi*TRACKPOINTS+ni];
       }
       else
       {
-        r=Trajectories[bi*TrackPoints+ni];
+        r=Trajectories[bi*TRACKPOINTS+ni];
       }
 	   
       if ((r<=re)&&(r>=ri))
@@ -889,7 +1000,6 @@
 
 __global__ void Circle(float *Trajectories,int *IdLast,
                        float *zImage,float *fImage,
-                       int TrackPoints,
                        float Mass,float InternalRadius,
                        float ExternalRadius,float Angle,
                        int Line)
@@ -924,7 +1034,7 @@
   bmx=1.25*re;
 
   // Angular step of integration
-  h=4.e0f*PI/(float)TrackPoints;
+  h=4.e0f*PI/(float)TRACKPOINTS;
 
   // impact parameter
   b=(float)bi/(float)bmaxi*bmx;
@@ -946,11 +1056,11 @@
 
      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.)+Trajectories[bi*TRACKPOINTS+ni];
      }
      else
      {
-	r=Trajectories[bi*TrackPoints+ni];
+	r=Trajectories[bi*TRACKPOINTS+ni];
      }
 	   
      if ((r<=re)&&(r>=ri))
@@ -967,8 +1077,7 @@
 
 }
 
-__global__ void Trajectory(float *Trajectories,
-                           int *IdLast,int TrackPoints,
+__global__ void Trajectory(float *Trajectories,int *IdLast,
                            float Mass,float InternalRadius,
                            float ExternalRadius,float Angle,
                            int Line)
@@ -1000,7 +1109,7 @@
   bmx=1.25*re;
 
   // Angular step of integration
-  h=4.e0f*PI/(float)TrackPoints;
+  h=4.e0f*PI/(float)TRACKPOINTS;
 
   // impact parameter
   b=(float)bi/(float)bmaxi*bmx;
@@ -1018,7 +1127,7 @@
   // b versus us
   float bvus=fabs(b/us);
   float bvus0=bvus;
-  Trajectories[bi*TrackPoints+nh]=bvus;
+  Trajectories[bi*TRACKPOINTS+nh]=bvus;
 
   do
   {
@@ -1028,13 +1137,126 @@
      vp=vs;
      rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
      bvus=fabs(b/us);
-     Trajectories[bi*TrackPoints+nh]=bvus;
+     Trajectories[bi*TRACKPOINTS+nh]=bvus;
 
   } while ((bvus>=rs)&&(bvus<=bvus0));
 
   IdLast[bi]=nh;
 
 }
+
+__global__ void Original(float *zImage,float *fImage,
+                         float Mass,float InternalRadius,
+                         float ExternalRadius,float Angle,
+                         int Line)
+{
+  // 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;
+
+  float Trajectories[2048];
+
+  // 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 bmx,db,b,h;
+  int nh;
+
+  // Autosize for image
+  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.e0f*(float)bmaxi);
+
+  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[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[nh]=bvus;
+
+  } while ((bvus>=rs)&&(bvus<=bvus0));
+
+  int imx=(int)(16*bi);
+
+  for (int i=0;i<imx;i++)
+  {
+     float zp=0,fp=0;
+     float phi=2.*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;
+     float php,nr,r;
+
+     do
+     {
+        php=phd+(float)HalfLap*PI;
+        nr=php/h;
+        ni=(int)nr;
+
+        if (ni<nh)
+        {
+           r=(Trajectories[ni+1]-Trajectories[ni])*(nr-ni*1.)+Trajectories[ni];
+        }
+        else
+        {
+           r=Trajectories[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));
+
+   __syncthreads();
+
+   zImage[yi+2*bmaxi*xi]=zp;
+   fImage[yi+2*bmaxi*xi]=fp;
+
+  }
+ 
+}
 """
     return(BlobCUDA)
     
@@ -1110,7 +1332,7 @@
     
     #    BlackHoleCL = cl.Program(ctx,KERNEL_CODE.substitute(kernel_params)).build()
 
-    BuildOptions="-cl-mad-enable -DPHYSICS=%i " % (PhysicsList[Physics])
+    BuildOptions="-cl-mad-enable -DPHYSICS=%i -DSETTRACKPOINTS=%i " % (PhysicsList[Physics],InputCL['TrackPoints'])
 
     BlackHoleCL = cl.Program(ctx,BlobOpenCL).build(options = BuildOptions)
     
@@ -1133,10 +1355,18 @@
                                        numpy.float32(Angle),
                                        numpy.int32(Line))
         CLLaunch.wait()
+    elif Method=='Original':
+        CLLaunch=BlackHoleCL.Original(queue,(zImage.shape[0]/2,),
+                                      None,zImageCL,fImageCL,
+                                      numpy.float32(Mass),
+                                      numpy.float32(InternalRadius),
+                                      numpy.float32(ExternalRadius),
+                                      numpy.float32(Angle),
+                                      numpy.int32(Line))
+        CLLaunch.wait()
     elif Method=='TrajectoCircle':
         CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],),
                                         None,TrajectoriesCL,IdLastCL,
-                                        numpy.uint32(Trajectories.shape[1]),
                                         numpy.float32(Mass),
                                         numpy.float32(InternalRadius),
                                         numpy.float32(ExternalRadius),
@@ -1147,7 +1377,6 @@
                                            zImage.shape[0]*4),None,
                                     TrajectoriesCL,IdLastCL,
                                     zImageCL,fImageCL,
-                                    numpy.uint32(Trajectories.shape[1]),
                                     numpy.float32(Mass),
                                     numpy.float32(InternalRadius),
                                     numpy.float32(ExternalRadius),
@@ -1157,7 +1386,6 @@
     else:
         CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],),
                                         None,TrajectoriesCL,IdLastCL,
-                                        numpy.uint32(Trajectories.shape[1]),
                                         numpy.float32(Mass),
                                         numpy.float32(InternalRadius),
                                         numpy.float32(ExternalRadius),
@@ -1168,7 +1396,6 @@
                                           zImage.shape[1]),None,
                                    zImageCL,fImageCL,TrajectoriesCL,IdLastCL,
                                    numpy.uint32(Trajectories.shape[0]),
-                                   numpy.uint32(Trajectories.shape[1]),
                                         numpy.float32(Mass),
                                         numpy.float32(InternalRadius),
                                         numpy.float32(ExternalRadius),
@@ -1246,12 +1473,13 @@
     Context=XPU.make_context()
 
     try:
-        mod = SourceModule(KernelCodeCuda(),options=['--compiler-options','-DPHYSICS=%i' % (PhysicsList[Physics])])
+        mod = SourceModule(KernelCodeCuda(),options=['--compiler-options','-DPHYSICS=%i -DSETTRACKPOINTS=%i' % (PhysicsList[Physics],TrackPoints)])
         print("Compilation seems to be OK")
     except:
         print("Compilation seems to break")
 
     EachPixelCU=mod.get_function("EachPixel")
+    OriginalCU=mod.get_function("Original")
     TrajectoryCU=mod.get_function("Trajectory")
     PixelCU=mod.get_function("Pixel")
     CircleCU=mod.get_function("Circle")
@@ -1276,9 +1504,17 @@
                     numpy.int32(Line),
                     grid=(zImage.shape[0]/32,zImage.shape[1]/32),
                     block=(32,32,1))
+    elif Method=='Original':
+        OriginalCU(zImageCU,fImageCU,
+                   numpy.float32(Mass),
+                   numpy.float32(InternalRadius),
+                   numpy.float32(ExternalRadius),
+                   numpy.float32(Angle),
+                   numpy.int32(Line),
+                   grid=(zImage.shape[0]/32/2,1),
+                   block=(32,1,1))
     elif Method=='TrajectoCircle':
         TrajectoryCU(TrajectoriesCU,IdLastCU,
-                     numpy.uint32(Trajectories.shape[1]),
                      numpy.float32(Mass),
                      numpy.float32(InternalRadius),
                      numpy.float32(ExternalRadius),
@@ -1288,7 +1524,6 @@
                      block=(32,1,1))
         
         CircleCU(TrajectoriesCU,IdLastCU,zImageCU,fImageCU,
-                 numpy.uint32(Trajectories.shape[1]),
                  numpy.float32(Mass),
                  numpy.float32(InternalRadius),
                  numpy.float32(ExternalRadius),
@@ -1298,7 +1533,6 @@
                  block=(32,32,1))
     else:
         TrajectoryCU(TrajectoriesCU,IdLastCU,
-                     numpy.uint32(Trajectories.shape[1]),
                      numpy.float32(Mass),
                      numpy.float32(InternalRadius),
                      numpy.float32(ExternalRadius),
@@ -1309,7 +1543,6 @@
         
         PixelCU(zImageCU,fImageCU,TrajectoriesCU,IdLastCU,
                 numpy.uint32(Trajectories.shape[0]),
-                numpy.uint32(Trajectories.shape[1]),
                 numpy.float32(Mass),
                 numpy.float32(InternalRadius),
                 numpy.float32(ExternalRadius),
@@ -1368,7 +1601,7 @@
     # No output as image
     NoImage = False
     
-    HowToUse='%s -h [Help] -b [BlackBodyEmission] -n [NoImage] -p <Einstein/Newton> -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -a <AngleAboveDisc> -d <DeviceId> -c <Greyscale/Red2Yellow> -g <CUDA/OpenCL> -t <EachPixel/TrajectoCircle/TrajectoPixel> -v <FP32/FP64>'
+    HowToUse='%s -h [Help] -b [BlackBodyEmission] -n [NoImage] -p <Einstein/Newton> -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -a <AngleAboveDisc> -d <DeviceId> -c <Greyscale/Red2Yellow> -g <CUDA/OpenCL> -t <EachPixel/TrajectoCircle/TrajectoPixel/Original> -v <FP32/FP64>'
 
     try:
         opts, args = getopt.getopt(sys.argv[1:],"hbns:m:i:x:a:d:g:v:t:c:p:",["blackbody","noimage","camera","size=","mass=","internal=","external=","angle=","device=","gpustyle=","variabletype=","method=","colors=","physics="])
Index: TrouNoir/BenchTrouNoirC.sh
===================================================================
--- TrouNoir/BenchTrouNoirC.sh	(revision 223)
+++ TrouNoir/BenchTrouNoirC.sh	(revision 224)
@@ -1,5 +1,6 @@
 #!/bin/bash
-SEQ="6 1 14"
+#SEQ="6 1 14"
+SEQ="6 1 13"
 HOST=$(hostname)
 DATE=$(date "+%Y%m%d-%H%M")
 LOGFILE=trou_noir_FP32-${HOST}-${DATE}.log
