Index: TrouNoir.py
===================================================================
--- TrouNoir.py	(revision 200)
+++ TrouNoir.py	(revision 201)
@@ -86,45 +86,48 @@
 	     float h,float r,float m,float bss)
 {
   float flx;
-  int fi;
 
-  fi=(int)(rf*nbr/bss);
-  flx=pow(r/m,q)*pow(rf,4)*b*db*h;
+  flx=pow(r/m,q)*pow(rf,4)*b;
   return(flx);
 }
 
 float spectre_cn(float rf,float b,float db,
 		float h,float r,float m,float bss)
 {
-  float flx;
-  float nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k;
+  double flx;
+  double nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k,psc2,psk;
   int fi,posfreq;
 
   planck=6.62e-34;
   k=1.38e-23;
+  psk=5.38e-11;
+  psc2=7.35e-51;
   temp=3.e7;
   m_point=1.e14;
   v=1.-3./r;
 
-  qu=1/sqrt(1-3./r)/sqrt(r)*(sqrt(r)-sqrt(6.)+sqrt(3.)/2.*log((sqrt(r)+sqrt(3.))/(sqrt(r)-sqrt(3.))*(sqrt(6.)-sqrt(3.))/(sqrt(6.)+sqrt(3.))));
+  qu=1./sqrt(1.-3./r)/sqrt(r)*(sqrt(r)-sqrt(6.)+sqrt(3.)/2.*log((sqrt(r)+sqrt(3.))/(sqrt(r)-sqrt(3.))*(sqrt(6.)-sqrt(3.))/(sqrt(6.)+sqrt(3.))));
 
-  temp_em=temp*sqrt(m)*exp(0.25*log(m_point))*exp(-0.75*log(r))*
-    exp(-0.125*log(v))*exp(0.25*log(qu));
+  temp_em=temp*sqrt(m)*exp(0.25*log(m_point))*exp(-0.75*log(r))*exp(-0.125*log(v))*exp(0.25*log(qu));
 
   flux_int=0;
   flx=0;
 
-  for (fi=1;fi<nbr;fi++)
+  for (fi=0;fi<nbr;fi++)
     {
-      nu_em=bss*fi/nbr;
+      nu_em=bss*(float)fi/(float)nbr;
       nu_rec=nu_em*rf; 
-      posfreq=1./bss*nu_rec*nbr;
+      posfreq=(int)(1./bss*nu_rec*(float)nbr);
       if ((posfreq>0)&&(posfreq<nbr))
 	{
-	  flux_int=2*planck/9.e16f*pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b*db*h;
+	  //flux_int=2*planck/9.e16f*pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b*db*h;
+	  //flux_int=pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b;
+	  flux_int=2.*psc2*pow(nu_em,3)/(exp(psk*nu_em/temp_em)-1.)*pow(rf,3)*b;
 	  flx+=flux_int;
 	}
     }
+  //return(pow(rf,3)*b*temp_em);
+  return((float)flx);
   return(flx);
 }
 
@@ -151,7 +154,7 @@
 
 }
 
-__kernel void EachPixel(__global float *MyImage)
+__kernel void EachPixel(__global float *zImage,__global float *fImage)
 {
    uint xi=(uint)get_global_id(0);
    uint yi=(uint)get_global_id(1);
@@ -185,7 +188,7 @@
   //tho=PI/180.*80;
   tho=PI/180.*80.;
   q=-2;
-  dist=1;
+  dist=0;
   raie=0;
 
   nmx=(float)sizex/2.;
@@ -196,11 +199,20 @@
   thx=asin(bmx/ro);
   pc=0;
 
+  if (raie==0)
+    {
+      bss=2.;
+    }
+  else
+    {
+      bss=1.e6;
+    }
+
   h=PI/5.e2f;
 
   // set origin as center of image
   float x=(float)xi-(float)(sizex/2)+(float)5e-1f;
-  float y=(float)(sizey/2)-(float)yi-(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));
@@ -239,7 +251,6 @@
      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)))||((fmod(ps,PI)<fmod(phd,PI))&&(fmod(pp,PI)>fmod(phd,PI))))&&(rp[nh]>ri)&&(rp[nh]<re)?1:0;          
      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));
@@ -250,12 +261,14 @@
   else
   {
      zp=0.;
+     fp=0.;
   }
 
 //
-  MyImage[yi+sizex*xi]=(float)zp;
+  zImage[yi+sizex*xi]=(float)zp;
+  fImage[yi+sizex*xi]=(float)fp;
 
-  barrier(CLK_LOCAL_MEM_FENCE);
+  barrier(CLK_GLOBAL_MEM_FENCE);
  
 }
 """)
@@ -269,7 +282,7 @@
     image = Image.fromarray(SigmaInt)
     image.save("%s.jpg" % prefix)
     
-def BlackHole(MyImage,Device):
+def BlackHole(zImage,fImage,Device):
 
     kernel_params = {}
 
@@ -297,25 +310,30 @@
     # Je recupere les flag possibles pour les buffers
     mf = cl.mem_flags
 
-    print(MyImage)
+    print(zImage)
     
     # Program based on Kernel2
-    MyImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=MyImage)
+    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() 
 
-    CLLaunch=BlackHoleCL.EachPixel(queue,(MyImage.shape[0],MyImage.shape[1]),None,MyImageCL)
+    CLLaunch=BlackHoleCL.EachPixel(queue,(numpy.int32(zImage.shape[0]),numpy.int32(zImage.shape[1])),None,zImageCL,fImageCL)
                                      
     CLLaunch.wait()
     elapsed = time.time()-start_time
     print("Elapsed %f: " % elapsed)
     
-    cl.enqueue_copy(queue,MyImage,MyImageCL).wait()
-
-    print(MyImage)
-    MyImageCL.release()
+    cl.enqueue_copy(queue,zImage,zImageCL).wait()
+    cl.enqueue_copy(queue,fImage,fImageCL).wait()
     
-    ImageOutput(MyImage,"TrouNoir")
+    print(zImage.max())
+    print(fImage.max())
+    zImageCL.release()
+    fImageCL.release()
+    
+    ImageOutput(zImage,"TrouNoirZ")
+    ImageOutput(fImage,"TrouNoirF")
     return(elapsed)
 
 if __name__=='__main__':
@@ -460,10 +478,11 @@
 #    print(Devices,Alu)
 
 #    MyImage=numpy.where(numpy.random.zeros(Size,Size)>0,1,-1).astype(numpy.float32)
-    MyImage=numpy.zeros((Size,Size),dtype=numpy.float32)
+    zImage=numpy.zeros((Size,Size),dtype=numpy.float32)
+    fImage=numpy.zeros((Size,Size),dtype=numpy.float32)
 
-    print(MyImage)
+    print(zImage)
 	
-    duration=BlackHole(MyImage,Device)
+    duration=BlackHole(zImage,fImage,Device)
 
 
