Index: TrouNoir.py
===================================================================
--- TrouNoir.py	(revision 224)
+++ TrouNoir.py	(revision 225)
@@ -532,20 +532,18 @@
  
 }
 
-__kernel void Original(__global float *zImage,__global float *fImage,
-                       float Mass,float InternalRadius,
-                       float ExternalRadius,float Angle,
-                       int Line)
+__kernel void EachCircle(__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];
+  float Trajectory[TRACKPOINTS];
 
-  // Perform trajectory for each pixel
-
   float m,rs,ri,re,tho;
   int raie,q;
 
@@ -583,7 +581,7 @@
   // b versus us
   float bvus=fabs(b/us);
   float bvus0=bvus;
-  Trajectories[nh]=bvus;
+  Trajectory[nh]=bvus;
 
   do
   {
@@ -593,15 +591,19 @@
      vp=vs;
      rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
      bvus=fabs(b/us);
-     Trajectories[nh]=bvus;
+     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++)
   {
-     float zp=0,fp=0;
+     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);
@@ -618,11 +620,11 @@
 
         if (ni<nh)
         {
-           r=(Trajectories[ni+1]-Trajectories[ni])*(nr-ni*1.)+Trajectories[ni];
+           r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni];
         }
         else
         {
-           r=Trajectories[ni];
+           r=Trajectory[ni];
         }
 	   
         if ((r<=re)&&(r>=ri))
@@ -643,6 +645,125 @@
   barrier(CLK_GLOBAL_MEM_FENCE);
  
 }
+
+__kernel void Original(__global float *zImage,__global float *fImage,
+                       uint Size,float Mass,float InternalRadius,
+                       float ExternalRadius,float Angle,
+                       int Line)
+{
+   // Integer Impact Parameter Size (half of image)
+   uint bmaxi=(uint)Size;
+
+   float Trajectory[TRACKPOINTS];
+
+   // 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;
+
+   // Integer Impact Parameter ID 
+   for (int bi=0;bi<bmaxi;bi++)
+   {
+      // 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;
+      Trajectory[nh]=bvus;
+
+      do
+      {
+         nh++;
+         pp=ps;
+         up=us;
+         vp=vs;
+         rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
+         bvus=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++)
+      {
+         float zp=0,fp=0;
+         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;
+         float php,nr,r;
+
+         do
+         {
+            php=phd+(float)HalfLap*PI;
+            nr=php/h;
+            ni=(int)nr;
+
+            if (ni<nh)
+            {
+               r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+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);
+ 
+}
 """
 
 def KernelCodeCuda():
@@ -1145,10 +1266,10 @@
 
 }
 
-__global__ void Original(float *zImage,float *fImage,
-                         float Mass,float InternalRadius,
-                         float ExternalRadius,float Angle,
-                         int Line)
+__global__ void EachCircle(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;
@@ -1156,7 +1277,7 @@
   // Integer Impact Parameter Size (half of image)
   int bmaxi=gridDim.x*blockDim.x;
 
-  float Trajectories[2048];
+  float Trajectory[2048];
 
   // Perform trajectory for each pixel
 
@@ -1197,7 +1318,7 @@
   // b versus us
   float bvus=fabs(b/us);
   float bvus0=bvus;
-  Trajectories[nh]=bvus;
+  Trajectory[nh]=bvus;
 
   do
   {
@@ -1207,7 +1328,7 @@
      vp=vs;
      rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
      bvus=fabs(b/us);
-     Trajectories[nh]=bvus;
+     Trajectory[nh]=bvus;
 
   } while ((bvus>=rs)&&(bvus<=bvus0));
 
@@ -1232,11 +1353,11 @@
 
         if (ni<nh)
         {
-           r=(Trajectories[ni+1]-Trajectories[ni])*(nr-ni*1.)+Trajectories[ni];
+           r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni];
         }
         else
         {
-           r=Trajectories[ni];
+           r=Trajectory[ni];
         }
 	   
         if ((r<=re)&&(r>=ri))
@@ -1257,6 +1378,123 @@
   }
  
 }
+
+__global__ void Original(float *zImage,float *fImage,
+                         uint Size,float Mass,float InternalRadius,
+                         float ExternalRadius,float Angle,
+                         int Line)
+{
+   // Integer Impact Parameter Size (half of image)
+   uint bmaxi=(uint)Size;
+
+   float Trajectory[TRACKPOINTS];
+
+   // 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;
+
+   // Integer Impact Parameter ID 
+   for (int bi=0;bi<bmaxi;bi++)
+   {
+      // 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;
+      Trajectory[nh]=bvus;
+
+      do
+      {
+         nh++;
+         pp=ps;
+         up=us;
+         vp=vs;
+         rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
+         bvus=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++)
+      {
+         float zp=0,fp=0;
+         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;
+         float php,nr,r;
+
+         do
+         {
+            php=phd+(float)HalfLap*PI;
+            nr=php/h;
+            ni=(int)nr;
+
+            if (ni<nh)
+            {
+               r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+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;
+  
+      }
+
+   } 
+   
+}
 """
     return(BlobCUDA)
     
@@ -1356,14 +1594,24 @@
                                        numpy.int32(Line))
         CLLaunch.wait()
     elif Method=='Original':
-        CLLaunch=BlackHoleCL.Original(queue,(zImage.shape[0]/2,),
+        CLLaunch=BlackHoleCL.Original(queue,(1,),
                                       None,zImageCL,fImageCL,
+                                      numpy.uint32(zImage.shape[0]/2),
                                       numpy.float32(Mass),
                                       numpy.float32(InternalRadius),
                                       numpy.float32(ExternalRadius),
                                       numpy.float32(Angle),
                                       numpy.int32(Line))
         CLLaunch.wait()
+    elif Method=='EachCircle':
+        CLLaunch=BlackHoleCL.EachCircle(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,
@@ -1480,6 +1728,7 @@
 
     EachPixelCU=mod.get_function("EachPixel")
     OriginalCU=mod.get_function("Original")
+    EachCircleCU=mod.get_function("EachCircle")
     TrajectoryCU=mod.get_function("Trajectory")
     PixelCU=mod.get_function("Pixel")
     CircleCU=mod.get_function("Circle")
@@ -1504,15 +1753,25 @@
                     numpy.int32(Line),
                     grid=(zImage.shape[0]/32,zImage.shape[1]/32),
                     block=(32,32,1))
+    elif Method=='EachCircle':
+        EachCircleCU(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=='Original':
         OriginalCU(zImageCU,fImageCU,
+                   numpy.uint32(zImage.shape[0]/2),
                    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))
+                   grid=(1,1),
+                   block=(1,1,1))
     elif Method=='TrajectoCircle':
         TrajectoryCU(TrajectoriesCU,IdLastCU,
                      numpy.float32(Mass),
@@ -1601,7 +1860,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/Original> -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/EachCircle/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="])
