Index: TrouNoir.py
===================================================================
--- TrouNoir.py	(revision 208)
+++ TrouNoir.py	(revision 209)
@@ -43,7 +43,7 @@
 KERNEL_CODE=string.Template("""
 
 #define PI (float)3.14159265359
-#define nbr 200
+#define nbr 256
 
 float atanp(float x,float y)
 {
@@ -107,47 +107,64 @@
 {
   float flx;
 
-  flx=pow(r/m,q)*pow(rf,4)*b;
+  flx=exp(q*log(r/m))*pow(rf,4)*b*db*h;
   return(flx);
 }
 
-float spectre_cn(float rf,float b,float db,float h,float r,float m,float bss)
+float spectre_cn(float rf32,float b32,float db32,
+		 float h32,float r32,float m32,float bss32)
 {
-  double flx;
-  double nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k,psc2,psk;
+  
+#define MYFLOAT double
+
+  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;
 
-  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;
+#define planck 6.62e-34
+#define k 1.38e-23
+#define c2 9.e16
+#define temp 3.e7
+#define m_point 1.
 
-  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.))));
+#define lplanck (log(6.62)-34.*log(10.))
+#define lk (log(1.38)-23.*log(10.))
+#define lc2 (log(9.)+16.*log(10.))
 
-  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));
+  MYFLOAT v=1.-3./r;
 
-  flux_int=0;
-  flx=0;
+  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*(float)fi/(float)nbr;
-      nu_rec=nu_em*rf; 
-      posfreq=(int)(1./bss*nu_rec*(float)nbr);
+      nu_em=bss*(MYFLOAT)fi/(MYFLOAT)nbr;
+      nu_rec=nu_em*rf;
+      posfreq=(int)(nu_rec*(MYFLOAT)nbr/bss);
       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=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;
+          // 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;
 	  flx+=flux_int;
 	}
     }
-  //return(pow(rf,3)*b*temp_em);
-  return((float)flx);
-  return(flx);
+
+  return((float)(flx));
 }
 
 void impact(float phi,float r,float b,float tho,float m,
@@ -161,14 +178,14 @@
 
   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);
     }
-  else
-    {
-      bss=1.e6;
-      flx=spectre_cn(rf,b,db,h,r,m,bss);
-    }
   
   *zp=1./rf;
   *fp=flx;
@@ -178,8 +195,7 @@
 __kernel void EachPixel(__global float *zImage,__global float *fImage,
                         float Mass,float InternalRadius,
                         float ExternalRadius,float Angle,
-                        float ObserverDistance,
-                        int BlackBody,int AngularCamera)
+                        int Line)
 {
    uint xi=(uint)get_global_id(0);
    uint yi=(uint)get_global_id(1);
@@ -188,18 +204,16 @@
 
    // Perform trajectory for each pixel, exit on hit
 
-  float m,rs,ri,re,tho,ro;
-  int q,dist,raie;
+  float m,rs,ri,re,tho;
+  int q,raie;
 
   m=Mass;
   rs=2.*m;
   ri=InternalRadius;
   re=ExternalRadius;
-  ro=ObserverDistance;
   tho=Angle;
   q=-2;
-  dist=AngularCamera;
-  raie=BlackBody;
+  raie=Line;
 
   float d,bmx,db,b,h;
   float rp[2048];
@@ -225,18 +239,11 @@
   // 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/2);
+  db=bmx/(float)(sizex);
 
-  if (dist==1) 
-  {
-     up=0.;
-     vp=1.;
-  }
-  else 
-  {
-     up=sin(thi);
-     vp=cos(thi);
-  }
+  up=0.;
+  vp=1.;
       
   pp=0.;
   nh=0;
@@ -279,8 +286,7 @@
                     uint ImpactParameter,uint TrackPoints,
                     float Mass,float InternalRadius,
                     float ExternalRadius,float Angle,
-                    float ObserverDistance,
-                    int BlackBody,int AngularCamera)
+                    int Line)
 {
    uint xi=(uint)get_global_id(0);
    uint yi=(uint)get_global_id(1);
@@ -289,32 +295,30 @@
 
    // Perform trajectory for each pixel
 
-  float m,rs,ri,re,tho,ro;
-  int q,dist,raie;
+  float m,rs,ri,re,tho;
+  int q,raie;
 
   m=Mass;
   rs=2.*m;
   ri=InternalRadius;
   re=ExternalRadius;
-  ro=ObserverDistance;
   tho=Angle;
   q=-2;
-  dist=AngularCamera;
-  raie=BlackBody;
+  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 externa radius
+  // 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/(float)ImpactParameter;
+  db=bmx/(2.e0*(float)ImpactParameter);
 
   // set origin as center of image
   float x=(float)xi-(float)(sizex/2)+(float)5e-1f;
@@ -371,8 +375,7 @@
                      int TrackPoints,
                      float Mass,float InternalRadius,
                      float ExternalRadius,float Angle,
-                     float ObserverDistance,
-                     int BlackBody,int AngularCamera)
+                     int Line)
 {
    // Integer Impact Parameter ID 
    int bi=get_global_id(0);
@@ -385,18 +388,15 @@
 
    // Perform trajectory for each pixel
 
-  float m,rs,ri,re,tho,ro;
-  int q,dist,raie;
+  float m,rs,ri,re,tho;
+  int q,raie;
 
   m=Mass;
   rs=2.*m;
   ri=InternalRadius;
   re=ExternalRadius;
-  ro=ObserverDistance;
   tho=Angle;
-  q=-2;
-  dist=AngularCamera;
-  raie=BlackBody;
+  raie=Line;
 
   float bmx,db,b,h;
   float phi,thi,phd;
@@ -411,7 +411,7 @@
 
   // impact parameter
   b=(float)bi/(float)bmaxi*bmx;
-  db=bmx/(float)bmaxi;
+  db=bmx/(2.e0*(float)bmaxi);
 
   phi=2.*PI/(float)imx*(float)i;
   phd=atanp(cos(phi)*sin(tho),cos(tho));
@@ -456,8 +456,7 @@
                          __global int *IdLast,int TrackPoints,
                          float Mass,float InternalRadius,
                          float ExternalRadius,float Angle,
-                         float ObserverDistance,
-                         int BlackBody,int AngularCamera)
+                         int Line)
 {
   // Integer Impact Parameter ID 
   int bi=get_global_id(0);
@@ -466,18 +465,16 @@
 
   // Perform trajectory for each pixel
 
-  float m,rs,ri,re,tho,ro;
-  int dist,raie,q;
+  float m,rs,ri,re,tho;
+  int raie,q;
 
   m=Mass;
   rs=2.*m;
   ri=InternalRadius;
   re=ExternalRadius;
-  ro=ObserverDistance;
   tho=Angle;
   q=-2;
-  dist=AngularCamera;
-  raie=BlackBody;
+  raie=Line;
 
   float d,bmx,db,b,h;
   float phi,thi,phd,php,nr,r;
@@ -495,17 +492,8 @@
 
   float up,vp,pp,us,vs,ps;
 
-  if (dist==1) 
-  {
-     up=0.;
-     vp=1.;
-  }
-  else 
-  {
-     thi=asin(b/ro);
-     up=sin(thi);
-     vp=cos(thi);
-  }
+  up=0.;
+  vp=1.;
       
   pp=0.;
   nh=0;
@@ -559,20 +547,16 @@
     InternalRadius=InputCL['InternalRadius']
     ExternalRadius=InputCL['ExternalRadius']
     Angle=InputCL['Angle']
-    ObserverDistance=InputCL['ObserverDistance']
     Method=InputCL['Method']
     TrackPoints=InputCL['TrackPoints']
 
     if InputCL['BlackBody']:
-        BlackBody=0
+        # Spectrum is Black Body one
+        Line=0
     else:
-        BlackBody=1
+        # Spectrum is Monochromatic Line one
+        Line=1
 
-    if InputCL['AngularCamera']:
-        AngularCamera=1
-    else:
-        AngularCamera=0
-        
     Trajectories=numpy.zeros((InputCL['Size']/2,InputCL['TrackPoints']),
                               dtype=numpy.float32)
     IdLast=numpy.zeros(InputCL['Size']/2,dtype=numpy.int32)
@@ -606,7 +590,6 @@
     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)
-    fImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=fImage)
     IdLastCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=IdLast)
     
     start_time=time.time() 
@@ -619,9 +602,7 @@
                                        numpy.float32(InternalRadius),
                                        numpy.float32(ExternalRadius),
                                        numpy.float32(Angle),
-                                       numpy.float32(ObserverDistance),
-                                       numpy.int32(BlackBody),
-                                       numpy.int32(AngularCamera))
+                                       numpy.int32(Line))
         CLLaunch.wait()
     elif Method=='TrajectoCircle':
         CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],),
@@ -631,9 +612,7 @@
                                         numpy.float32(InternalRadius),
                                         numpy.float32(ExternalRadius),
                                         numpy.float32(Angle),
-                                        numpy.float32(ObserverDistance),
-                                        numpy.int32(BlackBody),
-                                        numpy.int32(AngularCamera))
+                                        numpy.int32(Line))
         
         CLLaunch=BlackHoleCL.Circle(queue,(Trajectories.shape[0],
                                            zImage.shape[0]*4),None,
@@ -644,9 +623,7 @@
                                     numpy.float32(InternalRadius),
                                     numpy.float32(ExternalRadius),
                                     numpy.float32(Angle),
-                                    numpy.float32(ObserverDistance),
-                                    numpy.int32(BlackBody),
-                                    numpy.int32(AngularCamera))
+                                    numpy.int32(Line))
         CLLaunch.wait()
     else:
         CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],),
@@ -656,9 +633,7 @@
                                         numpy.float32(InternalRadius),
                                         numpy.float32(ExternalRadius),
                                         numpy.float32(Angle),
-                                        numpy.float32(ObserverDistance),
-                                        numpy.int32(BlackBody),
-                                        numpy.int32(AngularCamera))
+                                        numpy.int32(Line))
         
         CLLaunch=BlackHoleCL.Pixel(queue,(zImage.shape[0],
                                                   zImage.shape[1]),None,
@@ -669,9 +644,7 @@
                                         numpy.float32(InternalRadius),
                                         numpy.float32(ExternalRadius),
                                         numpy.float32(Angle),
-                                        numpy.float32(ObserverDistance),
-                                        numpy.int32(BlackBody),
-                                        numpy.int32(AngularCamera))
+                                        numpy.int32(Line))
         CLLaunch.wait()
     
     elapsed = time.time()-start_time
@@ -682,8 +655,8 @@
     cl.enqueue_copy(queue,Trajectories,TrajectoriesCL).wait()
     cl.enqueue_copy(queue,IdLast,IdLastCL).wait()
     
-    print(zImage.max())
-    print(fImage.max())
+    print(zImage.max(),numpy.where(zImage[:,:]==zImage.max()))
+    print(fImage.max(),numpy.where(fImage[:,:]==fImage.max()))
     zImageCL.release()
     fImageCL.release()
 
@@ -701,13 +674,10 @@
     #
     ExternalRadius=12.
     #
-    ObserverDistance=100.
     # Angle with normal to disc 10 degrees
     Angle = numpy.pi/180.*(90.-10.)
     # Radiation of disc : BlackBody or Monochromatic
-    BlackBody = True
-    # Camera : Angular Camera or plate with the size of camera
-    AngularCamera = True
+    BlackBody = False
     # Size of image
     Size=256
     # Variable Type
@@ -715,12 +685,12 @@
     # ?
     q=-2
     # Method of resolution
-    Method='Trajecto'
+    Method='TrajectoPixel'
           
-    HowToUse='%s -h [Help] -b [BlackBodyEmission] -c [AngularCamera] -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -o <ObservatorDistance> -a <AngleAboveDisc> -d <DeviceId> -g <CUDA/OpenCL> -t <EachPixel/TrajectoCircle/TrajectoPixel> -v <FP32/FP64>'
+    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>'
 
     try:
-        opts, args = getopt.getopt(sys.argv[1:],"hbcs:m:i:x:o:a:d:g:v:t:",["blackbody","camera","size=","mass=","internal=","external=","observer=","angle=","device=","gpustyle=","variabletype=","method="])
+        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="])
     except getopt.GetoptError:
         print(HowToUse % sys.argv[0])
         sys.exit(2)
@@ -777,14 +747,10 @@
             InternalRadius = float(arg)
         elif opt in ("-e", "--external"):
             ExternalRadius = float(arg)
-        elif opt in ("-o", "--observer"):
-            ObserverDistance = float(arg)
         elif opt in ("-a", "--angle"):
             Angle = numpy.pi/180.*(90.-float(arg))
         elif opt in ("-b", "--blackbody"):
             BlackBody = True
-        elif opt in ("-c", "--camera"):
-            AngularCamera = True
         elif opt in ("-t", "--method"):
             Method = arg
 
@@ -795,10 +761,8 @@
     print("Mass : %f" % Mass)
     print("Internal Radius : %f" % InternalRadius)
     print("External Radius : %f" % ExternalRadius)
-    print("Observer Distance : %f" % ObserverDistance)
     print("Angle with normal of (in radians) : %f" % Angle)
     print("Black Body Disc Emission (monochromatic instead) : %s" % BlackBody)
-    print("Angular Camera (dimension of object instead) : %s" % AngularCamera)
     print("Method of resolution : %s" % Method)
 
     if GpuStyle=='CUDA':
@@ -851,10 +815,8 @@
     InputCL['Mass']=Mass
     InputCL['InternalRadius']=InternalRadius
     InputCL['ExternalRadius']=ExternalRadius
-    InputCL['ObserverDistance']=ObserverDistance
     InputCL['Angle']=Angle
     InputCL['BlackBody']=BlackBody
-    InputCL['AngularCamera']=AngularCamera
     InputCL['Method']=Method
     InputCL['TrackPoints']=TrackPoints
     
