Index: TrouNoir/TrouNoir.py
===================================================================
--- TrouNoir/TrouNoir.py	(revision 210)
+++ TrouNoir/TrouNoir.py	(revision 211)
@@ -26,25 +26,25 @@
 
 import pyopencl as cl
 import numpy
-from PIL import Image
 import time,string
 from numpy.random import randint as nprnd
 import sys
 import getopt
 import matplotlib.pyplot as plt
+from socket import gethostname
 
+def DictionariesAPI():
+    PhysicsList={'Einstein':0,'Newton':1}
+    return(PhysicsList)
 
+BlobOpenCL= """
 
-
-
-
-
-
-KERNEL_CODE=string.Template("""
-
 #define PI (float)3.14159265359
 #define nbr 256
 
+#define EINSTEIN 0
+#define NEWTON 1
+
 float atanp(float x,float y)
 {
   float angle;
@@ -64,12 +64,18 @@
   return v;
 }
 
+#if PHYSICS == NEWTON
 float g(float u,float m,float b)
 {
+  return (-u);
+}
+#else
+float g(float u,float m,float b)
+{
   return (3.e0f*m/b*pow(u,2)-u);
 }
+#endif
 
-
 void calcul(float *us,float *vs,float up,float vp,
 	    float h,float m,float b)
 {
@@ -156,10 +162,13 @@
       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=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;
 	}
     }
@@ -522,17 +531,25 @@
   barrier(CLK_GLOBAL_MEM_FENCE);
  
 }
-""")
+"""
 
-def ImageOutput(sigma,prefix):
-    Max=sigma.max()
-    Min=sigma.min()
+# 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)
-    image.save("%s.jpg" % prefix)
-    
+#     # Normalize value as 8bits Integer
+#     SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
+#     image = Image.fromarray(SigmaInt)
+#     image.save("%s.jpg" % prefix)
+
+def ImageOutput(sigma,prefix,Colors):
+    import matplotlib.pyplot as plt
+    if Colors == 'Red2Yellow':
+        plt.imsave("%s.png" % prefix, sigma, cmap='afmhot')
+    else:
+        plt.imsave("%s.png" % prefix, sigma, cmap='Greys_r')
+
 def BlackHoleCL(zImage,fImage,InputCL):
 
     kernel_params = {}
@@ -549,7 +566,10 @@
     Angle=InputCL['Angle']
     Method=InputCL['Method']
     TrackPoints=InputCL['TrackPoints']
+    Physics=InputCL['Physics']
 
+    PhysicsList=DictionariesAPI()
+    
     if InputCL['BlackBody']:
         # Spectrum is Black Body one
         Line=0
@@ -580,13 +600,16 @@
     queue = cl.CommandQueue(ctx,
 			    properties=cl.command_queue_properties.PROFILING_ENABLE)
 
-    BlackHoleCL = cl.Program(ctx,KERNEL_CODE.substitute(kernel_params)).build()
     
+    #    BlackHoleCL = cl.Program(ctx,KERNEL_CODE.substitute(kernel_params)).build()
+
+    BuildOptions="-cl-mad-enable -DPHYSICS=%i " % (PhysicsList[Physics])
+
+    BlackHoleCL = cl.Program(ctx,BlobOpenCL).build(options = BuildOptions)
+    
     # Je recupere les flag possibles pour les buffers
     mf = cl.mem_flags
 
-    print(zImage)
-    
     TrajectoriesCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=Trajectories)
     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)
@@ -594,7 +617,6 @@
     
     start_time=time.time() 
 
-    print(Trajectories.shape)
     if Method=='EachPixel':
         CLLaunch=BlackHoleCL.EachPixel(queue,(zImage.shape[0],zImage.shape[1]),
                                        None,zImageCL,fImageCL,
@@ -648,15 +670,17 @@
         CLLaunch.wait()
     
     elapsed = time.time()-start_time
-    print("Elapsed %f: " % elapsed)
+    print("\nElapsed Time : %f" % elapsed)
     
     cl.enqueue_copy(queue,zImage,zImageCL).wait()
     cl.enqueue_copy(queue,fImage,fImageCL).wait()
     cl.enqueue_copy(queue,Trajectories,TrajectoriesCL).wait()
     cl.enqueue_copy(queue,IdLast,IdLastCL).wait()
-    
-    print(zImage.max(),numpy.where(zImage[:,:]==zImage.max()))
-    print(fImage.max(),numpy.where(fImage[:,:]==fImage.max()))
+
+    zMaxPosition=numpy.where(zImage[:,:]==zImage.max())
+    fMaxPosition=numpy.where(fImage[:,:]==fImage.max())
+    print("Z max @(%i,%i) : %f" % (zMaxPosition[0],zMaxPosition[1],zImage.max()))
+    print("Flux max @(%i,%i) : %f\n" % (fMaxPosition[0],fMaxPosition[1],fImage.max()))
     zImageCL.release()
     fImageCL.release()
 
@@ -686,11 +710,17 @@
     q=-2
     # Method of resolution
     Method='TrajectoPixel'
-          
-    HowToUse='%s -h [Help] -b [BlackBodyEmission] -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -a <AngleAboveDisc> -d <DeviceId> -g <CUDA/OpenCL> -t <EachPixel/TrajectoCircle/TrajectoPixel> -v <FP32/FP64>'
+    # Colors for output image
+    Colors='Greyscale'
+    # Physics
+    Physics='Einstein'
+    # 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>'
 
     try:
-        opts, args = getopt.getopt(sys.argv[1:],"hbs:m:i:x:a:d:g:v:t:",["blackbody","camera","size=","mass=","internal=","external=","angle=","device=","gpustyle=","variabletype=","method="])
+        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="])
     except getopt.GetoptError:
         print(HowToUse % sys.argv[0])
         sys.exit(2)
@@ -751,8 +781,14 @@
             Angle = numpy.pi/180.*(90.-float(arg))
         elif opt in ("-b", "--blackbody"):
             BlackBody = True
+        elif opt in ("-n", "--noimage"):
+            NoImage = True
         elif opt in ("-t", "--method"):
             Method = arg
+        elif opt in ("-c", "--colors"):
+            Colors = arg
+        elif opt in ("-p", "--physics"):
+            Physics = arg
 
     print("Device Identification selected : %s" % Device)
     print("GpuStyle used : %s" % GpuStyle)
@@ -764,6 +800,8 @@
     print("Angle with normal of (in radians) : %f" % Angle)
     print("Black Body Disc Emission (monochromatic instead) : %s" % BlackBody)
     print("Method of resolution : %s" % Method)
+    print("Colors for output images : %s" % Colors)
+    print("Physics used for Trajectories : %s" % Physics)
 
     if GpuStyle=='CUDA':
         print("\nSelection of CUDA device") 
@@ -819,8 +857,15 @@
     InputCL['BlackBody']=BlackBody
     InputCL['Method']=Method
     InputCL['TrackPoints']=TrackPoints
+    InputCL['Physics']=Physics
     
     duration=BlackHoleCL(zImage,fImage,InputCL)
 
-    ImageOutput(zImage,"TrouNoirZ_%s" % Method)
-    ImageOutput(fImage,"TrouNoirF_%s" % Method)
+    Hostname=gethostname()
+    Date=time.strftime("%Y%m%d_%H%M%S")
+    ImageInfo="%s_Device%i_%s_%s" % (Method,Device,Hostname,Date)
+    
+    
+    if not NoImage:
+        ImageOutput(zImage,"TrouNoirZ_%s" % ImageInfo,Colors)
+        ImageOutput(fImage,"TrouNoirF_%s" % ImageInfo,Colors)
