code: 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
    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);
  
    up=0.;
    vp=1.;
    pp=0.;
    nh=0;
  
    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);
       rpp=rps;
       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));
  
    if (ExitOnImpact==1) {
       impact(phi,rpp,b,tho,m,&zp,&fp,q,db,h,raie);
    }
    else
    {
       zp=0.e0f;
       fp=0.e0f;
    }
  
    __syncthreads();
  
    zImage[yi+sizex*xi]=(float)zp;
    fImage[yi+sizex*xi]=(float)fp;
  }
  
  __global__ void Pixel(float *zImage,float *fImage,
                        float *Trajectories,int *IdLast,
                        uint ImpactParameter,
                        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
  
    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,fp=0;
    // 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)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));
  
    // 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));
  
    }
  
    zImage[yi+sizex*xi]=zp;
    fImage[yi+sizex*xi]=fp;
  }
  
  __global__ void Circle(float *Trajectories,int *IdLast,
                         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;
     // Integer points on circle
     int i=blockIdx.y*blockDim.y+threadIdx.y;
     // Integer Impact Parameter Size (half of image)
     int bmaxi=gridDim.x*blockDim.x;
     // Integer Points on circle
     int imx=gridDim.y*blockDim.y;
  
     // 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,fp=0;
  
    // 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;
  
  }
  
  __global__ void Trajectory(float *Trajectories,int *IdLast,
                             float Mass,float InternalRadius,
                             float ExternalRadius,float Angle,
                             int Line)
  {
    // Integer Impact Parameter ID
    int bi=blockIdx.x*blockDim.x+threadIdx.x;
    // Integer Impact Parameter Size (half of image)
    int bmaxi=gridDim.x*blockDim.x;
  
    // 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;
  
  }
  
  __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;
  
    // Integer Impact Parameter Size (half of image)
    int bmaxi=gridDim.x*blockDim.x;
  
    float Trajectory[2048];
  
    // 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;
  
    // 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));
  
    int imx=(int)(16*bi);
  
    for (int i=0;i<imx;i++)
    {
       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);
       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));
  
     __syncthreads();
  
     zImage[yi+2*bmaxi*xi]=zp;
     fImage[yi+2*bmaxi*xi]=fp;
  
    }
  
  }
  
  __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.e0f*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
  
file_path: pyopencl/examples/black-hole-accretion.py
imports: import pyopencl
  import numpy
  import time
  import sys
  import getopt
  from socket import gethostname
label_desc: Use image node_id 2 for calling a global function with example usage: KernelCodeCuda() and returns: BlobCUDA