black-hole-accretion.py

#!/usr/bin/env python3
#
# TrouNoir model using PyOpenCL or PyCUDA
#
# CC BY-NC-SA 2019 : <emmanuel.quemener@ens-lyon.fr>
#
# Part of matrix programs from: https://forge.cbp.ens-lyon.fr/svn/bench4gpu/
#
# Thanks to Andreas Klockner for PyOpenCL and PyCUDA:
# http://mathema.tician.de/software/pyopencl
#
# Original code programmed in Fortran 77 in mars 1994
# for Practical Work of Numerical Simulation
# DEA (old Master2) in astrophysics and spatial techniques in Meudon
# by Herve Aussel & Emmanuel Quemener
#
# Conversion in C done by Emmanuel Quemener in august 1997
# GPUfication in OpenCL under Python in july 2019
# GPUfication in CUDA under Python in august 2019
#
# Thanks to :
#
# - Herve Aussel for his part of code of black body spectrum
# - Didier Pelat for his help to perform this work
# - Jean-Pierre Luminet for his article published in 1979
# - Numerical Recipes for Runge Kutta recipes
# - Luc Blanchet for his disponibility about my questions in General Relativity
# - Pierre Lena for his passion about science and vulgarisation

# If crash on OpenCL Intel implementation, add following options and force
# export PYOPENCL_COMPILER_OUTPUT=1
# export CL_CONFIG_USE_VECTORIZER=True
# export CL_CONFIG_CPU_VECTORIZER_MODE=16

import pyopencl as cl
import numpy
import time
import sys
import getopt
from socket import gethostname


def DictionariesAPI():
    PhysicsList = {"Einstein": 0, "Newton": 1}
    return PhysicsList


#
# Blank space below to simplify debugging on OpenCL code
#


BlobOpenCL = """

#define PI (float)3.14159265359e0f
#define nbr 256

#define EINSTEIN 0
#define NEWTON 1

#ifdef SETTRACKPOINTS
#define TRACKPOINTS SETTRACKPOINTS
#else
#define TRACKPOINTS 2048
#endif

float atanp(float x,float y)
{
  float angle;

  angle=atan2(y,x);

  if (angle<0.e0f)
    {
      angle+=(float)2.e0f*PI;
    }

  return angle;
}

float f(float v)
{
  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)
{
  float c0,c1,c2,c3,d0,d1,d2,d3;

  c0=h*f(vp);
  c1=h*f(vp+c0/2.e0f);
  c2=h*f(vp+c1/2.e0f);
  c3=h*f(vp+c2);
  d0=h*g(up,m,b);
  d1=h*g(up+d0/2.e0f,m,b);
  d2=h*g(up+d1/2.e0f,m,b);
  d3=h*g(up+d2,m,b);

  *us=up+(c0+2.e0f*c1+2.e0f*c2+c3)/6.e0f;
  *vs=vp+(d0+2.e0f*d1+2.e0f*d2+d3)/6.e0f;
}

void rungekutta(float *ps,float *us,float *vs,
                float pp,float up,float vp,
                float h,float m,float b)
{
  calcul(us,vs,up,vp,h,m,b);
  *ps=pp+h;
}

float decalage_spectral(float r,float b,float phi,
                        float tho,float m)
{
  return (sqrt(1-3*m/r)/(1+sqrt(m/pow(r,3))*b*sin(tho)*sin(phi)));
}

float spectre(float rf,int q,float b,float db,
              float h,float r,float m,float bss)
{
  float flx;

//  flx=exp(q*log(r/m))*pow(rf,4)*b*db*h;
  flx=exp(q*log(r/m)+4.e0f*log(rf))*b*db*h;
  return(flx);
}

float spectre_cn(float rf32,float b32,float db32,
                 float h32,float r32,float m32,float bss32)
{

#define MYFLOAT float

  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;

#define planck 6.62e-34f
#define k 1.38e-23f
#define c2 9.e16f
#define temp 3.e7f
#define m_point 1.e0f

#define lplanck (log(6.62e0f)-34.e0f*log(10.e0f))
#define lk (log(1.38e0f)-23.e0f*log(10.e0f))
#define lc2 (log(9.e0f)+16.e0f*log(10.e0f))

  MYFLOAT v=1.e0f-3.e0f/r;

  qu=1.e0f/sqrt((1.e0f-3.e0f/r)*r)*(sqrt(r)-sqrt(6.e0f)+sqrt(3.e0f)/2.e0f*log((sqrt(r)+sqrt(3.e0f))/(sqrt(r)-sqrt(3.e0f))* 0.17157287525380988e0f ));  // # noqa: E501

  temp_em=temp*sqrt(m)*exp(0.25e0f*log(m_point)-0.75e0f*log(r)-0.125e0f*log(v)+0.25e0f*log(fabs(qu)));

  flux_int=0.e0f;
  flx=0.e0f;

  for (fi=0;fi<nbr;fi++)
    {
      nu_em=bss*(MYFLOAT)fi/(MYFLOAT)nbr;
      nu_rec=nu_em*rf;
      posfreq=(int)(nu_rec*(MYFLOAT)nbr/bss);
      if ((posfreq>0)&&(posfreq<nbr))
        {
          // 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;
          //flux_int*=exp(3.e0f*log(rf))*b*db*h;
          flux_int=2.e0f*exp(lplanck-lc2+3.e0f*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.e0f)*exp(3.e0f*log(rf))*b*db*h;

          flx+=flux_int;
        }
    }

  return((float)(flx));
}

void impact(float phi,float r,float b,float tho,float m,
            float *zp,float *fp,
            int q,float db,
            float h,int raie)
{
  float flx,rf,bss;

  rf=decalage_spectral(r,b,phi,tho,m);

  if (raie==0)
    {
      bss=1.e19f;
      flx=spectre_cn(rf,b,db,h,r,m,bss);
    }
  else
    {
      bss=2.e0f;
      flx=spectre(rf,q,b,db,h,r,m,bss);
    }

  *zp=1.e0f/rf;
  *fp=flx;

}

__kernel void EachPixel(__global float *zImage,__global float *fImage,
                        float Mass,float InternalRadius,
                        float ExternalRadius,float Angle,
                        int Line)
{
   uint xi=(uint)get_global_id(0);
   uint yi=(uint)get_global_id(1);
   uint sizex=(uint)get_global_size(0);
   uint sizey=(uint)get_global_size(1);

   // Perform trajectory for each pixel, exit on hit

  float m,rs,ri,re,tho;
  int q,raie;

  m=Mass;
  rs=2.e0f*m;
  ri=InternalRadius;
  re=ExternalRadius;
  tho=Angle;
  q=-2;
  raie=Line;

  float bmx,db,b,h;
  float rp0,rps;
  float phi,phd;
  uint nh=0;
  float zp=0.e0f,fp=0.e0f;

  // Autosize for image
  bmx=1.25e0f*re;

  h=4.e0f*PI/(float)TRACKPOINTS;

  // set origin as center of image
  float x=(float)xi-(float)(sizex/2)+(float)5.e-1f;
  float y=(float)yi-(float)(sizey/2)+(float)5.e-1f;
  // angle extracted from cylindric symmetry
  phi=atanp(x,y);
  phd=atanp(cos(phi)*sin(tho),cos(tho));


  float up,vp,pp,us,vs,ps;

  // impact parameter
  b=sqrt(x*x+y*y)*(float)2.e0f/(float)sizex*bmx;
  // step of impact parameter;
  db=bmx/(float)(sizex);

  up=0.e0f;
  vp=1.e0f;
  pp=0.e0f;

  rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);

  rps=fabs(b/us);
  rp0=rps;

  int ExitOnImpact=0;

  do
  {
     nh++;
     pp=ps;
     up=us;
     vp=vs;
     rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
     rps=fabs(b/us);
     ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rps>=ri)&&(rps<=re)?1:0;

  } while ((rps>=rs)&&(rps<=rp0)&&(ExitOnImpact==0)&&(nh<TRACKPOINTS));


  if (ExitOnImpact==1) {
     impact(phi,rps,b,tho,m,&zp,&fp,q,db,h,raie);
  }
  else
  {
     zp=0.e0f;
     fp=0.e0f;
  }

  barrier(CLK_GLOBAL_MEM_FENCE);

  zImage[yi+sizex*xi]=(float)zp;
  fImage[yi+sizex*xi]=(float)fp;
}

__kernel void Pixel(__global float *zImage,__global float *fImage,
                    __global float *Trajectories,__global int *IdLast,
                    uint ImpactParameter,
                    float Mass,float InternalRadius,
                    float ExternalRadius,float Angle,
                    int Line)
{
   uint xi=(uint)get_global_id(0);
   uint yi=(uint)get_global_id(1);
   uint sizex=(uint)get_global_size(0);
   uint sizey=(uint)get_global_size(1);

   // Perform trajectory for each pixel

  float m,ri,re,tho;
  int q,raie;

  m=Mass;
  ri=InternalRadius;
  re=ExternalRadius;
  tho=Angle;
  q=-2;
  raie=Line;

  float bmx,db,b,h;
  float phi,phd,php,nr,r;
  float zp=0.e0f,fp=0.e0f;

  // Autosize for image, 25% greater than external radius
  bmx=1.25e0f*re;

  // Angular step of integration
  h=4.e0f*PI/(float)TRACKPOINTS;

  // Step of Impact Parameter
  db=bmx/(2.e0f*(float)ImpactParameter);

  // set origin as center of image
  float x=(float)xi-(float)(sizex/2)+(float)5.e-1f;
  float y=(float)yi-(float)(sizey/2)+(float)5.e-1f;

  // angle extracted from cylindric symmetry
  phi=atanp(x,y);
  phd=atanp(cos(phi)*sin(tho),cos(tho));

  // Real Impact Parameter
  b=sqrt(x*x+y*y)*bmx/(float)ImpactParameter;

  // Integer Impact Parameter
  uint bi=(uint)sqrt(x*x+y*y);

  int HalfLap=0,ExitOnImpact=0,ni;

  if (bi<ImpactParameter)
  {
    do
    {
      php=phd+(float)HalfLap*PI;
      nr=php/h;
      ni=(int)nr;

      if (ni<IdLast[bi])
      {
        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.e0f)+Trajectories[bi*TRACKPOINTS+ni];
      }
      else
      {
        r=Trajectories[bi*TRACKPOINTS+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));

  }

  barrier(CLK_GLOBAL_MEM_FENCE);

  zImage[yi+sizex*xi]=zp;
  fImage[yi+sizex*xi]=fp;
}

__kernel void Circle(__global float *Trajectories,__global int *IdLast,
                     __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 points on circle
   int i=get_global_id(1);
   // Integer Impact Parameter Size (half of image)
   int bmaxi=get_global_size(0);
   // Integer Points on circle
   int imx=get_global_size(1);

   // Perform trajectory for each pixel

  float m,ri,re,tho;
  int q,raie;

  m=Mass;
  ri=InternalRadius;
  re=ExternalRadius;
  tho=Angle;
  raie=Line;

  float bmx,db,b,h;
  float phi,phd;
  float zp=0.e0f,fp=0.e0f;

  // 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);

  phi=2.e0f*PI/(float)imx*(float)i;
  phd=atanp(cos(phi)*sin(tho),cos(tho));
  int yi=(int)((float)bi*sin(phi))+bmaxi;
  int xi=(int)((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<IdLast[bi])
     {
        r=(Trajectories[bi*TRACKPOINTS+ni+1]-Trajectories[bi*TRACKPOINTS+ni])*(nr-ni*1.e0f)+Trajectories[bi*TRACKPOINTS+ni];
     }
     else
     {
        r=Trajectories[bi*TRACKPOINTS+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);

}

__kernel void Trajectory(__global float *Trajectories,__global int *IdLast,
                         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);

  // Perform trajectory for each pixel

  float m,rs,re;

  m=Mass;
  rs=2.e0f*m;
  re=ExternalRadius;

  float bmx,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;

  float up,vp,pp,us,vs,ps;

  up=0.e0f;
  vp=1.e0f;

  pp=0.e0f;
  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 EachCircle(__global float *zImage,__global float *fImage,
                         float Mass,float InternalRadius,
                         float ExternalRadius,float Angle,
                         int Line)
{
   // Integer Impact Parameter ID
   uint bi=(uint)get_global_id(0);
   // Integer Impact Parameter Size (half of image)
   uint bmaxi=(uint)get_global_size(0);

  private float Trajectory[TRACKPOINTS];

  float m,rs,ri,re,tho;
  int raie,q;

  m=Mass;
  rs=2.e0f*m;
  ri=InternalRadius;
  re=ExternalRadius;
  tho=Angle;
  q=-2;
  raie=Line;

  float bmx,db,b,h;
  uint 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.e0f;
  vp=1.e0f;

  pp=0.e0f;
  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=(float)fabs(b/us);
     Trajectory[nh]=bvus;

  } while ((bvus>=rs)&&(bvus<=bvus0));


  for (uint i=(uint)nh+1;i<TRACKPOINTS;i++) {
     Trajectory[i]=0.e0f;
  }


  uint imx=(uint)(16*bi);

  for (uint i=0;i<imx;i++)
  {
     float zp=0.e0f,fp=0.e0f;
     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);

     uint 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.e0f)+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);

}

__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.e0f*m;
   ri=InternalRadius;
   re=ExternalRadius;
   tho=Angle;
   q=-2;
   raie=Line;

   float bmx,db,b,h;
   uint 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.e0f;
      vp=1.e0f;

      pp=0.e0f;
      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.e0f,fp=0.e0f;
         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);

         uint 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.e0f)+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():
    BlobCUDA = """

#define PI (float)3.14159265359
#define nbr 256

#define EINSTEIN 0
#define NEWTON 1

#ifdef SETTRACKPOINTS
#define TRACKPOINTS SETTRACKPOINTS
#else
#define TRACKPOINTS
#endif
__device__ float nothing(float x)
{
  return(x);
}

__device__ float atanp(float x,float y)
{
  float angle;

  angle=atan2(y,x);

  if (angle<0.e0f)
    {
      angle+=(float)2.e0f*PI;
    }

  return(angle);
}

__device__ float f(float v)
{
  return(v);
}

#if PHYSICS == NEWTON
__device__ float g(float u,float m,float b)
{
  return (-u);
}
#else
__device__ float g(float u,float m,float b)
{
  return (3.e0f*m/b*pow(u,2)-u);
}
#endif

__device__ void calcul(float *us,float *vs,float up,float vp,
                       float h,float m,float b)
{
  float c0,c1,c2,c3,d0,d1,d2,d3;

  c0=h*f(vp);
  c1=h*f(vp+c0/2.);
  c2=h*f(vp+c1/2.);
  c3=h*f(vp+c2);
  d0=h*g(up,m,b);
  d1=h*g(up+d0/2.,m,b);
  d2=h*g(up+d1/2.,m,b);
  d3=h*g(up+d2,m,b);

  *us=up+(c0+2.*c1+2.*c2+c3)/6.;
  *vs=vp+(d0+2.*d1+2.*d2+d3)/6.;
}

__device__ void rungekutta(float *ps,float *us,float *vs,
                           float pp,float up,float vp,
                           float h,float m,float b)
{
  calcul(us,vs,up,vp,h,m,b);
  *ps=pp+h;
}

__device__ float decalage_spectral(float r,float b,float phi,
                                   float tho,float m)
{
  return (sqrt(1-3*m/r)/(1+sqrt(m/pow(r,3))*b*sin(tho)*sin(phi)));
}

__device__ float spectre(float rf,int q,float b,float db,
                         float h,float r,float m,float bss)
{
  float flx;

//  flx=exp(q*log(r/m))*pow(rf,4)*b*db*h;
  flx=exp(q*log(r/m)+4.*log(rf))*b*db*h;
  return(flx);
}

__device__ float spectre_cn(float rf32,float b32,float db32,
                            float h32,float r32,float m32,float bss32)
{

#define MYFLOAT float

  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;

#define planck 6.62e-34
#define k 1.38e-23
#define c2 9.e16
#define temp 3.e7
#define m_point 1.

#define lplanck (log(6.62)-34.*log(10.))
#define lk (log(1.38)-23.*log(10.))
#define lc2 (log(9.)+16.*log(10.))

  MYFLOAT v=1.-3./r;

  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 ));  // # noqa: #051

  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*(MYFLOAT)fi/(MYFLOAT)nbr;
      nu_rec=nu_em*rf;
      posfreq=(int)(nu_rec*(MYFLOAT)nbr/bss);
      if ((posfreq>0)&&(posfreq<nbr))
        {
          // 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;
          //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;
        }
    }

  return((float)(flx));
}

__device__ void impact(float phi,float r,float b,float tho,float m,
                       float *zp,float *fp,
                       int q,float db,
                       float h,int raie)
{
  float flx,rf,bss;

  rf=decalage_spectral(r,b,phi,tho,m);

  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);
    }

  *zp=1./rf;
  *fp=flx;

}

__global__ void EachPixel(float *zImage,float *fImage,
                          float Mass,float InternalRadius,
                          float ExternalRadius,float Angle,
                          int Line)
{
   uint xi=(uint)(blockIdx.x*blockDim.x+threadIdx.x);
   uint yi=(uint)(blockIdx.y*blockDim.y+threadIdx.y);
   uint sizex=(uint)gridDim.x*blockDim.x;
   uint sizey=(uint)gridDim.y*blockDim.y;


   // Perform trajectory for each pixel, exit on hit

  float m,rs,ri,re,tho;
  int q,raie;

  m=Mass;
  rs=2.*m;
  ri=InternalRadius;
  re=ExternalRadius;
  tho=Angle;
  q=-2;
  raie=Line;

  float bmx,db,b,h;
  float rp0,rpp,rps;
  float phi,phd;
  int nh;
  float zp,fp;

  // Autosize for image
  bmx=1.25*re;
  b=0.;

  h=4.e0f*PI/(float)TRACKPOINTS;

  // set origin as center of image
  float x=(float)xi-(float)(sizex/2)+(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));

  float up,vp,pp,us,vs,ps;

  // impact parameter