#line 79 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
#include "imaging.h"

#line 313 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void pointing_kern(float *theta_x, float *theta_y, 
			      float  *zenith, float *azimuth)
{
  int k;
  k = threadIdx.x;
  theta_x[k] = tanf(zenith[k])*cosf(azimuth[k]);
  theta_y[k] = tanf(zenith[k])*sinf(azimuth[k]);
}

#line 1437 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  __global__ void setupRandomSequence_imaging(curandState *state, const int n_data, 
					      const int LOCAL_RAND_SEED)
{
  int id = threadIdx.x + blockIdx.x*blockDim.x;
  if (id<n_data)
    curand_init(LOCAL_RAND_SEED, id, 0, &state[id]);
}
#line 1466 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__device__ float curand_poisson_knuth(curandState *state, float lambda)
{
  unsigned int k = 0;
  float p=exp(lambda);
  do{
      k++;
      p *= curand_uniform(state);
  }while (p > 1.0);
 return (float) k-1;
}
#line 1480 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__device__ float curand_poisson_atkinson(curandState *state, float lambda)
{
  /*
    The algorithm is ???method PA??? from ???The Computer Generation of Poisson Random Variables??? by A. C. Atkinson, Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979), pages 29-35. Method PA is on page 32.
   */
  float c, beta, alpha, k, u, x, n, v, y, lhs, rhs;
  c = 0.767 - 3.36/lambda;
  beta = PI/sqrt(3.0*lambda);
  alpha = beta*lambda;
  k = logf(c) - lambda - logf(beta);
  do{
    u = curand_uniform(state);
    x = (alpha - logf((1.0 - u)/u))/beta;
    n = floorf(x + 0.5);
    if (n >= 0)
      {
	v = curand_uniform(state);
	y = alpha - beta*x;
        lhs = 1.0 + expf(y);
	lhs *= lhs;
	lhs = y + logf(v/lhs);
	rhs = k + n*logf(lambda) - lgammaf(n+1);
      }
  }while(lhs > rhs);
  return n;
}
#line 1365 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void noisify(curandState *state, float *frame, const int n_frame,
			const float readOutNoiseRms, float exposureTime,
			const float nBackgroundPhoton,
			const float noiseFactor2,
			const float photoelectron_gain)
   {
     int i, j, id;
     float lambda;
     curandState localState;
     i = blockIdx.x * blockDim.x + threadIdx.x;
     j = blockIdx.y * blockDim.y + threadIdx.y;
     id = j * gridDim.x * blockDim.x + i;
     //int id = threadIdx.x + blockIdx.x*blockDim.x;
     if ( id < n_frame )
       {
         localState = state[id];
         frame[id] += nBackgroundPhoton;
         frame[id] *= photoelectron_gain;
         frame[id] *= noiseFactor2;
         frame[id] *= exposureTime;
         lambda = frame[id];
         //       frame[id] = (float) curand_poisson(&localState, (double) lambda);
           /*
         if (lambda>0) {
  	 frame[id] = (lambda<30) ?
  	   curand_poisson_knuth( &localState, lambda) :
  	   curand_poisson_atkinson( &localState, lambda);
         }
    */
	 frame[id] = (lambda<64) ?
	   curand_poisson_knuth( &localState, lambda) :
	   (lambda<4000) ?
	   (float) curand_poisson_gammainc( &localState, lambda) :
	   curand_normal(&localState)*sqrt(lambda) + lambda;
	 frame[id] += (1-noiseFactor2)*lambda/noiseFactor2;
	 frame[id] = (frame[id]<0) ? 0 : frame[id]; 
	 frame[id] += (readOutNoiseRms>0) ?
	   curand_normal(&localState)*readOutNoiseRms : 0;
	 //frame[id] -= nBackgroundPhoton*exposureTime*photoelectron_gain;
	 //frame[id] = (frame[id]>0) ? floorf(frame[id]) : ceilf(frame[id]);
	 state[id] = localState;
       }
     }
__global__ void flux_integral(float *frame, const int n_frame,
                        float exposureTime,
                        const float photoelectron_gain)
{
  int i, j, id;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  id = j * gridDim.x * blockDim.x + i;
  //int id = threadIdx.x + blockIdx.x*blockDim.x;
  if ( id < n_frame )
    {
      frame[id] *= photoelectron_gain;
      frame[id] *= exposureTime;
    }
}
#line 1004 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void complexAmplitude(float2 *wave_PUPIL, int N_DFT,
                                 float *wavefront_amplitude,
				 float *wavefront_phase,
				 float wavenumber,
				 int N_PX_PUPIL, int N_PX_IMAGE, int M,
				 int N_SIDE_LENSLET, source *src,
				 float pixel_scale,
				 float *theta_x, float *theta_y,
				 char absolute_pointing, char DFT_SHIFT)
{
  int i, j, iLenslet, jLenslet, ij_PUPIL, ij_DFT, iSource;
  float c ,s, phasor, dx, dy, dx_int, dy_int;
  dx = dy = 0.0;
#line 1062 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  // source index
  iSource = blockIdx.z;
#line 1066 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  //  wave pixel coordinate (N_DFTLxN_SIDE_LENSLET)^2
  if ( threads2lenslet(threadIdx, blockIdx, blockDim,
		       &i, &j, N_DFT,
		       &iLenslet, &jLenslet, N_SIDE_LENSLET) ) {
#line 1072 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    ij_DFT = lenslet2row(i,j,N_DFT,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
    if ((i<N_PX_PUPIL) && (j<N_PX_PUPIL))
      {
#line 1077 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
ij_PUPIL = lenslet2array_overlap(i,j,N_PX_PUPIL,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
#line 1090 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        phasor = (DFT_SHIFT==0) ? 0.0 : PI*(i+j)*(1-N_DFT + (N_PX_IMAGE & 1))/N_DFT;
        
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 1092 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
	phasor -= 2*PI*(i*dx+j*dy)/N_DFT;
#line 1094 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        phasor += wavefront_phase[ij_PUPIL]*wavenumber;
        sincosf(phasor,&c,&s);
#line 1105 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = wavefront_amplitude[ij_PUPIL];
        wave_PUPIL[ij_DFT].x *= c;
        wave_PUPIL[ij_DFT].y *= s;
        return;
      }
    if ((i<N_DFT) && (j<N_DFT))
      wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = 0;
  }
#line 1020 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
__global__ void complexAmplitudeTT7(float2 *wave_PUPIL, int N_DFT,
                                 float *wavefront_amplitude,
				 float *wavefront_phase,
				 float wavenumber,
                                 int *piston_mask,
				 int N_PX_PUPIL, int N_PX_IMAGE, int M,
				 int N_SIDE_LENSLET, source *src,
				 float pixel_scale,
				 float *theta_x, float *theta_y,
				 char absolute_pointing, char DFT_SHIFT)
{
  int i, j, iLenslet, jLenslet, ij_PUPIL, ij_DFT, iSource;
  float c ,s, phasor, dx, dy, dx_int, dy_int;
  dx = dy = 0.0;
#line 1062 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  // source index
  iSource = blockIdx.z;
#line 1066 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  //  wave pixel coordinate (N_DFTLxN_SIDE_LENSLET)^2
  if ( threads2lenslet(threadIdx, blockIdx, blockDim,
		       &i, &j, N_DFT,
		       &iLenslet, &jLenslet, N_SIDE_LENSLET) ) {
#line 1072 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    ij_DFT = lenslet2row(i,j,N_DFT,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
    if ((i<N_PX_PUPIL) && (j<N_PX_PUPIL))
      {
#line 1036 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
ij_PUPIL = lenslet2array_overlap(i,j,N_PX_PUPIL,iLenslet,jLenslet,N_SIDE_LENSLET,0);
  
#line 1090 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        phasor = (DFT_SHIFT==0) ? 0.0 : PI*(i+j)*(1-N_DFT + (N_PX_IMAGE & 1))/N_DFT;
        
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 1092 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
	phasor -= 2*PI*(i*dx+j*dy)/N_DFT;
#line 1094 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        phasor += wavefront_phase[ij_PUPIL]*wavenumber;
        sincosf(phasor,&c,&s);
#line 1114 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = wavefront_amplitude[ij_PUPIL]*
                                               (piston_mask[ij_PUPIL]==(iSource+1));
        wave_PUPIL[ij_DFT].x *= c;
        wave_PUPIL[ij_DFT].y *= s;
        return;
      }
    if ((i<N_DFT) && (j<N_DFT))
      wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = 0;
  }
#line 1038 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 1041 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void complexAmplitudeWithModulation(float2 *wave_PUPIL, int N_DFT,
                                 float *wavefront_amplitude,
				 float *wavefront_phase,
				 float wavenumber,
				 int N_PX_PUPIL, int N_PX_IMAGE, int M,
				 int N_SIDE_LENSLET, source *src,
				 float pixel_scale,
				 float *theta_x, float *theta_y,
				 char absolute_pointing, char DFT_SHIFT,
                                 float modulation_r, float modulation_o)
{
  int i, j, iLenslet, jLenslet, ij_PUPIL, ij_DFT, iSource;
  float c ,s, phasor, dx, dy, dx_int, dy_int;
  dx = dy = 0.0;
#line 1062 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  // source index
  iSource = blockIdx.z;
#line 1066 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  //  wave pixel coordinate (N_DFTLxN_SIDE_LENSLET)^2
  if ( threads2lenslet(threadIdx, blockIdx, blockDim,
		       &i, &j, N_DFT,
		       &iLenslet, &jLenslet, N_SIDE_LENSLET) ) {
#line 1072 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    ij_DFT = lenslet2row(i,j,N_DFT,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
    if ((i<N_PX_PUPIL) && (j<N_PX_PUPIL))
      {
#line 1077 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
ij_PUPIL = lenslet2array_overlap(i,j,N_PX_PUPIL,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
#line 1090 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        phasor = (DFT_SHIFT==0) ? 0.0 : PI*(i+j)*(1-N_DFT + (N_PX_IMAGE & 1))/N_DFT;
        
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 1092 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
	phasor -= 2*PI*(i*dx+j*dy)/N_DFT;
#line 1099 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        sincosf(modulation_o,&s,&c);
        phasor -= 8*PI*modulation_r*(i*c+j*s)/N_DFT;
#line 1094 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        phasor += wavefront_phase[ij_PUPIL]*wavenumber;
        sincosf(phasor,&c,&s);
#line 1105 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = wavefront_amplitude[ij_PUPIL];
        wave_PUPIL[ij_DFT].x *= c;
        wave_PUPIL[ij_DFT].y *= s;
        return;
      }
    if ((i<N_DFT) && (j<N_DFT))
      wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = 0;
  }
#line 1058 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 1126 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void complexAmplitudeNoOverlap(float2 *wave_PUPIL, int N_DFT,
                                 float *wavefront_amplitude,
				 float *wavefront_phase,
				 float wavenumber,
				 int N_PX_PUPIL, int N_PX_IMAGE, int M,
				 int N_SIDE_LENSLET, source *src,
				 float pixel_scale,
				 float *theta_x, float *theta_y,
				 char absolute_pointing, char DFT_SHIFT)
{
  int i, j, iLenslet, jLenslet, ij_PUPIL, ij_DFT, iSource;
  float c ,s, phasor, dx, dy, dx_int, dy_int;
  dx = dy = 0.0;
#line 1062 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  // source index
  iSource = blockIdx.z;
#line 1066 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  //  wave pixel coordinate (N_DFTLxN_SIDE_LENSLET)^2
  if ( threads2lenslet(threadIdx, blockIdx, blockDim,
		       &i, &j, N_DFT,
		       &iLenslet, &jLenslet, N_SIDE_LENSLET) ) {
#line 1072 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    ij_DFT = lenslet2row(i,j,N_DFT,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
    if ((i<N_PX_PUPIL) && (j<N_PX_PUPIL))
      {
#line 1146 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
ij_PUPIL = lenslet2array(i,j,N_PX_PUPIL,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
#line 1090 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        phasor = (DFT_SHIFT==0) ? 0.0 : PI*(i+j)*(1-N_DFT + (N_PX_IMAGE & 1))/N_DFT;
        
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 1092 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
	phasor -= 2*PI*(i*dx+j*dy)/N_DFT;
#line 1094 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
        phasor += wavefront_phase[ij_PUPIL]*wavenumber;
        sincosf(phasor,&c,&s);
#line 1105 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = wavefront_amplitude[ij_PUPIL];
        wave_PUPIL[ij_DFT].x *= c;
        wave_PUPIL[ij_DFT].y *= s;
        return;
      }
    if ((i<N_DFT) && (j<N_DFT))
      wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = 0;
  }
#line 1142 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
} 
#line 1149 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void complexAmplitudeNoOverlapSPS(float2 *wave_PUPIL, int N_DFT,
                                 float *wavefront_amplitude,
				 float *wavefront_phase,
				 float wavenumber,
				 int N_PX_PUPIL, int N_PX_IMAGE, int M,
				 int N_SIDE_LENSLET, source *src,
				 float d, char DFT_SHIFT)
{
  int i, j, iLenslet, jLenslet, kLenslet,
    ij_PUPIL, ij_DFT, iSource;
  float c ,s, phasor, theta, ctheta, stheta;
  
#line 1062 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  // source index
  iSource = blockIdx.z;
#line 1066 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  //  wave pixel coordinate (N_DFTLxN_SIDE_LENSLET)^2
  if ( threads2lenslet(threadIdx, blockIdx, blockDim,
		       &i, &j, N_DFT,
		       &iLenslet, &jLenslet, N_SIDE_LENSLET) ) {
#line 1072 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    ij_DFT = lenslet2row(i,j,N_DFT,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
    if ((i<N_PX_PUPIL) && (j<N_PX_PUPIL))
      {
#line 1161 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  
#line 1146 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
ij_PUPIL = lenslet2array(i,j,N_PX_PUPIL,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
#line 1162 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  kLenslet = jLenslet*N_SIDE_LENSLET + iLenslet;
  kLenslet = kLenslet%12;
  if (kLenslet<6)
    theta = (3-kLenslet)*PI/3.0;
  else
    theta = (2-(kLenslet-6))*PI/3.0;
  sincosf(theta,&ctheta,&stheta);
  phasor = (DFT_SHIFT==0) ? 0.0 : PI*(i+j)*(1-N_DFT + (N_PX_IMAGE & 1))/N_DFT;
  phasor -= 2*PI*d*(i*ctheta+j*stheta)/N_DFT;
  phasor += wavefront_phase[ij_PUPIL]*wavenumber;
  sincosf(phasor,&c,&s);
  wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = wavefront_amplitude[ij_PUPIL];
  wave_PUPIL[ij_DFT].x *= c;
  wave_PUPIL[ij_DFT].y *= s;
  return;
}
if ((i<N_DFT) && (j<N_DFT))
  wave_PUPIL[ij_DFT].x = wave_PUPIL[ij_DFT].y = 0;
}
}
#line 701 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void wave2intensity(float2 *wave_PUPIL, const int N, float intensity_norm)
{
  int i, j, ij;
  float x ,y;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  if ( ij < N )
    {
      x =  wave_PUPIL[ij].x;
      y =  wave_PUPIL[ij].y;
      wave_PUPIL[ij].x = x*x + y*y;
      wave_PUPIL[ij].y = 0.0;
      wave_PUPIL[ij].x *= intensity_norm;
    }
}
#line 782 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void apply_field_stop(float2 *d__wave_PUPIL, int N_DFT,
				 int N_PX_PUPIL, float field_stop_diam)
{
  int i, j, k, w, half;
  float xi, yj, rij, l, c, s, phasor, a, b;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (i<N_DFT) && (j<N_DFT) ) {
    w = N_DFT & 1;
    half = (N_DFT-2+w)/2;
    if (i>half)
      xi = (float) (i - N_DFT);
    else
      xi = (float) i;
    if (j>half)
      yj = (float) (j - N_DFT);
    else
      yj = (float) j;
    rij = hypot(xi,yj);
    l = 0.5*field_stop_diam;
    k = i*N_DFT + j;
    if (rij>l) {
      d__wave_PUPIL[k].x = d__wave_PUPIL[k].y = 0.0;
    } else {
      phasor = 2.*(i+j)*(N_PX_PUPIL-1)/N_DFT;
      sincospif(phasor,&s,&c);
      a = d__wave_PUPIL[k].x;
      b = d__wave_PUPIL[k].y;
      d__wave_PUPIL[k].x = a*c - b*s;
      d__wave_PUPIL[k].y = b*c + a*s;
    }
  }
}
#line 858 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void apply_pyramid(float2 *d__wave_PUPIL, int N_DFT, int N_PX_PUPIL, float alpha)
{
  int i, j, k, w, half, iSource;
  float xi, yj, c, s, phasor, a, b;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  iSource = blockIdx.z;
  if ( (i<N_DFT) && (j<N_DFT) ) {
    k = i*N_DFT + j + iSource*N_DFT*N_DFT;
    w = N_PX_PUPIL & 1;
    half = (N_DFT-2+w)/2;
    if (i>half)
      i -= N_DFT;
    if (j>half)
      j -= N_DFT;
    xi = (float) i;
    yj = (float) j;
    phasor = -2*(i+j)*(N_PX_PUPIL*1.5 + 1 - w*0.5)/N_DFT;
    phasor -= alpha*(fabs(xi)+fabs(yj));
    sincospif(phasor,&s,&c);
    a = d__wave_PUPIL[k].x;
    b = d__wave_PUPIL[k].y;
    d__wave_PUPIL[k].x = (a*c - b*s)/N_DFT;
    d__wave_PUPIL[k].y = (b*c + a*s)/N_DFT;
  }
}
#line 634 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void apply_irradiance(float2 *wave_PUPIL, const int N_DFT,
				   const int N_SIDE_LENSLET, const float fwhm)
{
  int i, j, iLenslet, jLenslet, ij_DFT, iSource, w, half, N2;
  float xi, yj, sigma2, fft_I;
  iSource = blockIdx.z;
  if ( threads2lenslet(threadIdx, blockIdx, blockDim,
                       &i, &j, N_DFT,
                       &iLenslet, &jLenslet, N_SIDE_LENSLET) ) {
    ij_DFT = lenslet2row(i,j,N_DFT,iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
    w = N_DFT & 1;
    half = (N_DFT-2+w)/2;
    if (i>half)
      xi = (float) (i - N_DFT);
    else
      xi = (float) i;
    if (j>half)
      yj = (float) (j - N_DFT);
    else
      yj = (float) j;
    xi /= N_DFT;
    yj /= N_DFT;
    sigma2 = 8*logf(2);
    sigma2 = fwhm*fwhm/sigma2;
    fft_I = 2*PI*PI*sigma2;
    fft_I *= xi*xi + yj*yj;
    fft_I = expf( -fft_I);
    N2 = N_DFT*N_DFT;
    wave_PUPIL[ij_DFT].x *= fft_I/N2;
    wave_PUPIL[ij_DFT].y *= fft_I/N2;
  }
}

#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void binning(float *frame, const float2 *wave, int N_DFT,
			source *src, float pixel_scale,
			float *theta_x, float *theta_y, char absolute_pointing,
			int N_PX_IMAGE, int N_PX_CAMERA, int N_SIDE_LENSLET)
{
  int u, v, iLenslet, jLenslet, ij_DFT,
    M, offset, iPxL, jPxL, i, j, ij_CAM, ii, jj,iSource;
  float dx, dy, dx_int, dy_int;
  dx_int = dy_int = 0.0;

  // source index
  iSource = blockIdx.z;

  // BINNING
  if ( threads2lenslet(threadIdx, blockIdx, blockDim,
		       &u, &v, N_PX_CAMERA,
		       &iLenslet, &jLenslet, N_SIDE_LENSLET) )
    {

      M = N_PX_IMAGE/N_PX_CAMERA;
      offset = (N_PX_IMAGE-N_DFT)/2;

      __shared__ float bin[16][16];
      ii = threadIdx.x;
      jj = threadIdx.y;

      bin[ii][jj] = 0.0;
      
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 1242 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
      for (i=0;i<M;i++)
      	{
      	  iPxL = u*M + i - offset - dx_int;
      	  if ( iPxL>=0 && iPxL<N_DFT)
      	    for (j=0;j<M;j++)
      	      {
      		jPxL = v*M + j - offset - dy_int;
      		if ( jPxL>=0 && jPxL<N_DFT)
      		  {
      		    ij_DFT = lenslet2row(iPxL,jPxL,N_DFT,
					 iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
      		    bin[ii][jj] += src[iSource].N_PHOTON*wave[ij_DFT].x;
      		  }
      	      }
      	}
      /* u += dx_int; */
      /* v += dy_int; */
      ij_CAM = lenslet2array(u,v,N_PX_CAMERA,
			     iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
      frame[ij_CAM] +=   bin[ii][jj]/(N_DFT*N_DFT);
    }

}
#line 1267 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
__global__ void binningTT7(float *frame, const float2 *wave, int N_DFT,
			source *src, float pixel_scale,
			float *theta_x, float *theta_y, char absolute_pointing,
			int N_PX_IMAGE, int N_PX_CAMERA, int N_SIDE_LENSLET)
{
  int u, v, iLenslet, jLenslet, ij_DFT,
    M, offset, iPxL, jPxL, i, j, ij_CAM, ii, jj,iSource;
  float dx, dy, dx_int, dy_int;
  dx_int = dy_int = 0.0;

  // source index
  iSource = blockIdx.z;

  // BINNING
  if ( threads2lenslet(threadIdx, blockIdx, blockDim,
		       &u, &v, N_PX_CAMERA,
		       &iLenslet, &jLenslet, N_SIDE_LENSLET) )
    {

      M = N_PX_IMAGE/N_PX_CAMERA;
      offset = (N_PX_IMAGE-N_DFT)/2;

      __shared__ float bin[16][16];
      ii = threadIdx.x;
      jj = threadIdx.y;

      bin[ii][jj] = 0.0;
      
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 1295 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
      for (i=0;i<M;i++)
      	{
      	  iPxL = u*M + i - offset - dx_int;
      	  if ( iPxL>=0 && iPxL<N_DFT)
      	    for (j=0;j<M;j++)
      	      {
      		jPxL = v*M + j - offset - dy_int;
      		if ( jPxL>=0 && jPxL<N_DFT)
      		  {
      		    ij_DFT = lenslet2row(iPxL,jPxL,N_DFT,
					 iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
      		    bin[ii][jj] += src[0].N_PHOTON*wave[ij_DFT].x;
      		  }
      	      }
      	}
      /* u += dx_int; */
      /* v += dy_int; */
      ij_CAM = lenslet2array(u,v,N_PX_CAMERA,
			     iLenslet,jLenslet,N_SIDE_LENSLET,iSource);
      frame[ij_CAM] +=   bin[ii][jj]/(N_DFT*N_DFT);
    }

}
#line 1326 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::readout(float exposureTime, float readOutNoiseRms) {
#line 1334 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
float unit_exposureTime;
unit_exposureTime = (N_FRAME>0) ? exposureTime/N_FRAME : exposureTime;
N_PHOTON_PER_FRAME *= unit_exposureTime;
N = N_PX_CAMERA * N_PX_CAMERA * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
dim3 blockDim(16,16);
dim3 gridDim(Nb,Nb);
#line 1328 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
noisify <<< gridDim, blockDim >>> (devStates, d__frame, N,
				   readOutNoiseRms, unit_exposureTime, 
				   0.0, 1.0, photoelectron_gain);
}
#line 1425 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  void imaging::noiseless_readout(float exposureTime) {
       
#line 1334 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
float unit_exposureTime;
unit_exposureTime = (N_FRAME>0) ? exposureTime/N_FRAME : exposureTime;
N_PHOTON_PER_FRAME *= unit_exposureTime;
N = N_PX_CAMERA * N_PX_CAMERA * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
dim3 blockDim(16,16);
dim3 gridDim(Nb,Nb);
#line 1427 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
       flux_integral <<< gridDim, blockDim >>> (d__frame, N,
					  unit_exposureTime,
					  photoelectron_gain);
   }
#line 1345 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::readout(float exposureTime, float readOutNoiseRms, 
		      float nBackgroundPhoton,  float noiseFactor) {
#line 1334 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
float unit_exposureTime;
unit_exposureTime = (N_FRAME>0) ? exposureTime/N_FRAME : exposureTime;
N_PHOTON_PER_FRAME *= unit_exposureTime;
N = N_PX_CAMERA * N_PX_CAMERA * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
dim3 blockDim(16,16);
dim3 gridDim(Nb,Nb);
#line 1348 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
float nbgp, nf2;
int n;
nbgp = (N_FRAME>0) ? nBackgroundPhoton*N_FRAME : nBackgroundPhoton;
/*
n = N_SIDE_LENSLET*N_PX_CAMERA;
n *= n;
nbgp /= n;
*/
nf2 = noiseFactor*noiseFactor;
noisify <<< gridDim, blockDim >>> (devStates, d__frame, N,
				   readOutNoiseRms, unit_exposureTime,
				   nbgp, nf2,
				   photoelectron_gain);
}

#line 166 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::setup(int __N_PX_PUPIL, int __N_SIDE_LENSLET,
		    int DFT_osf, float IMAGE_osf, float CAMERA_osf) {
  N_SOURCE    = 1;
  
#line 181 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
BIN_IMAGE   = 1;
N_PX_PUPIL  = __N_PX_PUPIL;
N_DFT       = DFT_osf*N_PX_PUPIL;//*round_up_to_nhp2(N_PX_PUPIL);
N_PX_IMAGE  = IMAGE_osf*N_PX_PUPIL;
N_PX_CAMERA = CAMERA_osf*N_PX_IMAGE;
N_SIDE_LENSLET = __N_SIDE_LENSLET;
N_LENSLET   = N_SIDE_LENSLET*N_SIDE_LENSLET;
N_FRAME     = 0;
#line 170 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  
#line 225 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 248 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int n_byte = sizeof(float) * N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__zenith  , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__azimuth , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_x , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_y , n_byte) );      
HANDLE_ERROR( cudaMemset(d__zenith, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__azimuth, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_x, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_y, 0, n_byte) );
#line 227 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
pixel_scale = 0.0;
absolute_pointing = 0;
photoelectron_gain = 1.0;
INFO("@(CEO)>imaging: Device memory allocation\n");
int N_DFT_bytes = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__wave_PUPIL , sizeof(float2) * N_DFT_bytes) );
HANDLE_ERROR( cudaMalloc( (void**)&d__frame ,
			  sizeof(float) *N_PX_CAMERA*N_PX_CAMERA*N_LENSLET*N_SOURCE ) );
reset();
#line 686 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
INFO("@(CEO)>imaging: Creating a 2D FFT plan\n");
int n_DFT[NRANK] = {N_DFT, N_DFT};
int iodist = N_DFT*N_DFT;
int BATCH = N_LENSLET*N_SOURCE;
/* Create a 2D FFT plan. */
HANDLE_ERROR_CUFFT( cufftPlanMany(&plan, NRANK, n_DFT,
				  NULL, 1, iodist,
				  NULL, 1, iodist,
				  CUFFT_C2C,BATCH),
		    "Unable to create plan");
/* HANDLE_ERROR_CUFFT( cufftSetCompatibilityMode(plan, CUFFT_COMPATIBILITY_NATIVE), */
/* 		    "Unable to set compatibility mode to native"); */
#line 237 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
INFO("@(CEO)> Setting up random generator\n");
int n_data;
n_data = N_SOURCE*N_LENSLET*N_PX_CAMERA*N_PX_CAMERA;
HANDLE_ERROR( cudaMalloc( (void**)&devStates, n_data*sizeof(curandState)) );
dim3 blockDim(256);
dim3 gridDim(n_data/256+1);
setupRandomSequence_imaging <<< gridDim,blockDim >>> (devStates, n_data, LOCAL_RAND_SEED);
info();
#line 171 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 173 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::setup(int __N_PX_PUPIL, int __N_SIDE_LENSLET,
		    int DFT_osf, float IMAGE_osf, float CAMERA_osf,
		    int __N_SOURCE) {
  N_SOURCE    = __N_SOURCE;
  
#line 181 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
BIN_IMAGE   = 1;
N_PX_PUPIL  = __N_PX_PUPIL;
N_DFT       = DFT_osf*N_PX_PUPIL;//*round_up_to_nhp2(N_PX_PUPIL);
N_PX_IMAGE  = IMAGE_osf*N_PX_PUPIL;
N_PX_CAMERA = CAMERA_osf*N_PX_IMAGE;
N_SIDE_LENSLET = __N_SIDE_LENSLET;
N_LENSLET   = N_SIDE_LENSLET*N_SIDE_LENSLET;
N_FRAME     = 0;
#line 178 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  
#line 225 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 248 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int n_byte = sizeof(float) * N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__zenith  , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__azimuth , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_x , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_y , n_byte) );      
HANDLE_ERROR( cudaMemset(d__zenith, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__azimuth, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_x, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_y, 0, n_byte) );
#line 227 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
pixel_scale = 0.0;
absolute_pointing = 0;
photoelectron_gain = 1.0;
INFO("@(CEO)>imaging: Device memory allocation\n");
int N_DFT_bytes = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__wave_PUPIL , sizeof(float2) * N_DFT_bytes) );
HANDLE_ERROR( cudaMalloc( (void**)&d__frame ,
			  sizeof(float) *N_PX_CAMERA*N_PX_CAMERA*N_LENSLET*N_SOURCE ) );
reset();
#line 686 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
INFO("@(CEO)>imaging: Creating a 2D FFT plan\n");
int n_DFT[NRANK] = {N_DFT, N_DFT};
int iodist = N_DFT*N_DFT;
int BATCH = N_LENSLET*N_SOURCE;
/* Create a 2D FFT plan. */
HANDLE_ERROR_CUFFT( cufftPlanMany(&plan, NRANK, n_DFT,
				  NULL, 1, iodist,
				  NULL, 1, iodist,
				  CUFFT_C2C,BATCH),
		    "Unable to create plan");
/* HANDLE_ERROR_CUFFT( cufftSetCompatibilityMode(plan, CUFFT_COMPATIBILITY_NATIVE), */
/* 		    "Unable to set compatibility mode to native"); */
#line 237 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
INFO("@(CEO)> Setting up random generator\n");
int n_data;
n_data = N_SOURCE*N_LENSLET*N_PX_CAMERA*N_PX_CAMERA;
HANDLE_ERROR( cudaMalloc( (void**)&devStates, n_data*sizeof(curandState)) );
dim3 blockDim(256);
dim3 gridDim(n_data/256+1);
setupRandomSequence_imaging <<< gridDim,blockDim >>> (devStates, n_data, LOCAL_RAND_SEED);
info();
#line 179 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 191 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::setup(int __N_PX_PUPIL, int __N_SIDE_LENSLET,
		    int DFT_osf, int N_PX_IMAGE_, float CAMERA_osf,
		    int __N_SOURCE) {
  N_SOURCE    = __N_SOURCE;
  
#line 199 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
BIN_IMAGE   = 1;
N_PX_PUPIL  = __N_PX_PUPIL;
N_DFT       = DFT_osf*N_PX_PUPIL;//round_up_to_nhp2(N_PX_PUPIL);
N_PX_IMAGE  = N_PX_IMAGE_;
N_PX_CAMERA = CAMERA_osf*N_PX_IMAGE;
N_SIDE_LENSLET = __N_SIDE_LENSLET;
N_LENSLET   = N_SIDE_LENSLET*N_SIDE_LENSLET;
N_FRAME     = 0;
#line 196 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  
#line 225 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 248 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int n_byte = sizeof(float) * N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__zenith  , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__azimuth , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_x , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_y , n_byte) );      
HANDLE_ERROR( cudaMemset(d__zenith, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__azimuth, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_x, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_y, 0, n_byte) );
#line 227 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
pixel_scale = 0.0;
absolute_pointing = 0;
photoelectron_gain = 1.0;
INFO("@(CEO)>imaging: Device memory allocation\n");
int N_DFT_bytes = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__wave_PUPIL , sizeof(float2) * N_DFT_bytes) );
HANDLE_ERROR( cudaMalloc( (void**)&d__frame ,
			  sizeof(float) *N_PX_CAMERA*N_PX_CAMERA*N_LENSLET*N_SOURCE ) );
reset();
#line 686 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
INFO("@(CEO)>imaging: Creating a 2D FFT plan\n");
int n_DFT[NRANK] = {N_DFT, N_DFT};
int iodist = N_DFT*N_DFT;
int BATCH = N_LENSLET*N_SOURCE;
/* Create a 2D FFT plan. */
HANDLE_ERROR_CUFFT( cufftPlanMany(&plan, NRANK, n_DFT,
				  NULL, 1, iodist,
				  NULL, 1, iodist,
				  CUFFT_C2C,BATCH),
		    "Unable to create plan");
/* HANDLE_ERROR_CUFFT( cufftSetCompatibilityMode(plan, CUFFT_COMPATIBILITY_NATIVE), */
/* 		    "Unable to set compatibility mode to native"); */
#line 237 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
INFO("@(CEO)> Setting up random generator\n");
int n_data;
n_data = N_SOURCE*N_LENSLET*N_PX_CAMERA*N_PX_CAMERA;
HANDLE_ERROR( cudaMalloc( (void**)&devStates, n_data*sizeof(curandState)) );
dim3 blockDim(256);
dim3 gridDim(n_data/256+1);
setupRandomSequence_imaging <<< gridDim,blockDim >>> (devStates, n_data, LOCAL_RAND_SEED);
info();
#line 197 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 209 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::setup(int __N_PX_PUPIL, int __N_SIDE_LENSLET,
		    int DFT_osf, int N_PX_IMAGE_, int _BIN_IMAGE_,
		    int __N_SOURCE) {
  N_SOURCE    = __N_SOURCE;
  BIN_IMAGE   = _BIN_IMAGE_;
  N_PX_PUPIL  = __N_PX_PUPIL+1;
  N_DFT       = DFT_osf*N_PX_PUPIL;//round_up_to_nhp2(N_PX_PUPIL);
  N_PX_IMAGE  = N_PX_IMAGE_;
  N_PX_CAMERA = N_PX_IMAGE/BIN_IMAGE;
  N_SIDE_LENSLET = __N_SIDE_LENSLET;
  N_LENSLET   = N_SIDE_LENSLET*N_SIDE_LENSLET;
  N_FRAME     = 0;
  
#line 225 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 248 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int n_byte = sizeof(float) * N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__zenith  , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__azimuth , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_x , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_y , n_byte) );      
HANDLE_ERROR( cudaMemset(d__zenith, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__azimuth, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_x, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_y, 0, n_byte) );
#line 227 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
pixel_scale = 0.0;
absolute_pointing = 0;
photoelectron_gain = 1.0;
INFO("@(CEO)>imaging: Device memory allocation\n");
int N_DFT_bytes = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__wave_PUPIL , sizeof(float2) * N_DFT_bytes) );
HANDLE_ERROR( cudaMalloc( (void**)&d__frame ,
			  sizeof(float) *N_PX_CAMERA*N_PX_CAMERA*N_LENSLET*N_SOURCE ) );
reset();
#line 686 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
INFO("@(CEO)>imaging: Creating a 2D FFT plan\n");
int n_DFT[NRANK] = {N_DFT, N_DFT};
int iodist = N_DFT*N_DFT;
int BATCH = N_LENSLET*N_SOURCE;
/* Create a 2D FFT plan. */
HANDLE_ERROR_CUFFT( cufftPlanMany(&plan, NRANK, n_DFT,
				  NULL, 1, iodist,
				  NULL, 1, iodist,
				  CUFFT_C2C,BATCH),
		    "Unable to create plan");
/* HANDLE_ERROR_CUFFT( cufftSetCompatibilityMode(plan, CUFFT_COMPATIBILITY_NATIVE), */
/* 		    "Unable to set compatibility mode to native"); */
#line 237 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
INFO("@(CEO)> Setting up random generator\n");
int n_data;
n_data = N_SOURCE*N_LENSLET*N_PX_CAMERA*N_PX_CAMERA;
HANDLE_ERROR( cudaMalloc( (void**)&devStates, n_data*sizeof(curandState)) );
dim3 blockDim(256);
dim3 gridDim(n_data/256+1);
setupRandomSequence_imaging <<< gridDim,blockDim >>> (devStates, n_data, LOCAL_RAND_SEED);
info();
#line 222 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 259 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::setupSegmentPistonSensor(int __N_PX_PUPIL, int __N_SIDE_LENSLET,
                                       int _N_DFT_, int N_PX_IMAGE_, int _BIN_IMAGE_,
                                       int __N_SOURCE) {
  N_SOURCE    = __N_SOURCE;
  BIN_IMAGE   = _BIN_IMAGE_;
  N_PX_PUPIL  = __N_PX_PUPIL+1;
  N_DFT       = _N_DFT_;
  N_PX_IMAGE  = N_PX_IMAGE_;
  N_PX_CAMERA = N_PX_IMAGE/BIN_IMAGE;
  N_SIDE_LENSLET = __N_SIDE_LENSLET;
  N_LENSLET   = N_SIDE_LENSLET*N_SIDE_LENSLET;
  N_FRAME     = 0;

  
#line 248 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int n_byte = sizeof(float) * N_SOURCE;
HANDLE_ERROR( cudaMalloc( (void**)&d__zenith  , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__azimuth , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_x , n_byte) );      
HANDLE_ERROR( cudaMalloc( (void**)&d__theta_y , n_byte) );      
HANDLE_ERROR( cudaMemset(d__zenith, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__azimuth, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_x, 0, n_byte) );
HANDLE_ERROR( cudaMemset(d__theta_y, 0, n_byte) );

#line 274 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  pixel_scale = 0.0;
  absolute_pointing = 0;
  photoelectron_gain = 1.0;

  int N_DFT_bytes = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
  cudaMalloc( (void**)&d__wave_PUPIL , sizeof(float2) * N_DFT_bytes);

  int n_DFT[NRANK] = {N_DFT, N_DFT};
  int iodist = N_DFT*N_DFT;
  int BATCH = N_LENSLET*N_SOURCE;
  /* Create a 2D FFT plan. */
  HANDLE_ERROR_CUFFT( cufftPlanMany(&plan, NRANK, n_DFT,
				    NULL, 1, iodist,
				    NULL, 1, iodist,
				    CUFFT_C2C,BATCH),
		      "Unable to create plan");
  /* HANDLE_ERROR_CUFFT( cufftSetCompatibilityMode(plan, CUFFT_COMPATIBILITY_NATIVE), */
  /* 		      "Unable to set compatibility mode to native"); */
}

#line 325 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::cleanup(void) {
  INFO("@(CEO)>imaging: Device memory free\n");
  cufftDestroy(plan);
  cudaFree( d__wave_PUPIL );
  cudaFree( d__frame );
  cudaFree( devStates );
  cudaFree(d__zenith);
  cudaFree(d__azimuth);
  cudaFree(d__theta_x);
  cudaFree(d__theta_y);
}
void imaging::cleanupSegmentPistonSensor(void) {
  cufftDestroy(plan);
  cudaFree( d__wave_PUPIL );
}


#line 297 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::set_pointing_direction(float *zen, float *azim) {
  int n_byte = sizeof(float) * N_SOURCE;
  HANDLE_ERROR( cudaMemcpy( d__zenith, zen, n_byte,
			    cudaMemcpyHostToDevice ) );
  HANDLE_ERROR( cudaMemcpy( d__azimuth, azim, n_byte,
			    cudaMemcpyHostToDevice ) );
  pointing_kern <<< 1, N_SOURCE >>> (d__theta_x, d__theta_y, d__zenith, d__azimuth);
  /*
  zenith = zen;
  azimuth = azim;
  theta_x = tanf(zenith)*cosf(azimuth);
  theta_y = tanf(zenith)*sinf(azimuth);
  */
}

#line 732 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::reset(void) {
  HANDLE_ERROR( cudaMemset(d__frame, 0,
			   sizeof(float)*N_PX_CAMERA*N_PX_CAMERA*N_LENSLET*N_SOURCE) );
  N_FRAME = 0;
}

#line 1448 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::reset_rng(int SEED) {
  int n_data;
  n_data = N_SOURCE*N_LENSLET*N_PX_CAMERA*N_PX_CAMERA;
  LOCAL_RAND_SEED = SEED;
  dim3 blockDim(256);
  dim3 gridDim(n_data/256+1);
  setupRandomSequence_imaging <<< gridDim,blockDim >>> (devStates, n_data, LOCAL_RAND_SEED);  
}


#line 1540 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::info(void)
{
#ifndef SILENT
  fprintf(stdout,"\n\x1B[1;42m@(CEO)>imaging:\x1B[;42m\n");
  fprintf(stdout," . imaging device    : %4d\n",N_LENSLET);
  fprintf(stdout," . pupil sampling    : %4d\n",N_PX_PUPIL);
  fprintf(stdout," . DFT sampling      : %4d\n",N_DFT);
  fprintf(stdout," . imagelet sampling : %4d\n",N_PX_IMAGE);
  fprintf(stdout," . frame sampling    : %4d\n",N_PX_CAMERA);
  fprintf(stdout," . source #          : %4d\n",N_SOURCE);
  fprintf(stdout,"----------------------------------------------------\x1B[0m\n");
#endif
}

#line 346 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::propagate_cpx(source *src) {
  int M = N_PX_IMAGE/N_PX_CAMERA;
  dim3 gridDim;
  dim3 blockDim(16,16);
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  complexAmplitude <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
                                              src->wavefront.amplitude,
                                              src->wavefront.phase,
                                              src->wavenumber(),
                                              N_PX_PUPIL, N_PX_IMAGE, M,
                                              N_SIDE_LENSLET, src->dev_ptr,
                                              pixel_scale, d__theta_x, d__theta_y,
                                              absolute_pointing,1);

  HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
                      "Unable to execute plan forward FT");
}
#line 368 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::propagate(source *src) {
  int M = N_PX_IMAGE/N_PX_CAMERA;
  /* if ( (pixel_scale>0) && (zenith<0) ) */
  /*   set_pointing_direction(0.0,0.0); */

  /*
  if (pixel_scale>0 ) {
      pixel_scale /= BIN_IMAGE;
      printf(" pixel scale: %.2E\n",pixel_scale);
      float dx, dy, dx_int, dy_int;
      int iSource=0;
      
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 380 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
      printf("ANGULAR OFFSET:\n");
      printf(" . x: %3.0f + %4.3f\n",dx_int,dx);
      printf(" . y: %3.0f + %4.3f\n",dy_int,dy);
    }
  */

  //absolute_pointing = 1;
  //  fprintf(stdout," pixel scale: %.2E\n",pixel_scale);
  //  fprintf(stdout,"Absolute pointing: %d\n",absolute_pointing);
  dim3 gridDim;
  dim3 blockDim(16,16);
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  complexAmplitude <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
                                              src->wavefront.amplitude,
                                              src->wavefront.phase,
                                              src->wavenumber(),
                                              N_PX_PUPIL, N_PX_IMAGE, M,
                                              N_SIDE_LENSLET, src->dev_ptr,
                                              pixel_scale, d__theta_x, d__theta_y,
                                              absolute_pointing,1);
  
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 721 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
N_PHOTON_PER_FRAME = N_PHOTON_PER_SECOND_PER_FRAME = src->n_photon()*src->wavefront.M->area;
float intensity_norm;
intensity_norm = src->wavefront.M->area/src->wavefront.M->nnz;
/*
fprintf(stdout,"Source # of photon: %g\n",src->n_photon());
fprintf(stdout,"Wavefront mask area and nnz: %f & %f\n",src->wavefront.M->area,src->wavefront.M->nnz);
fprintf(stdout,"# of photon per second: %g\n",N_PHOTON_PER_SECOND_PER_FRAME);
*/
#line 598 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
gridDim  = dim3(Nb,Nb);
wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, intensity_norm);
if (src->fwhm>0) {
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  
#line 622 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_INVERSE),
		    "Unable to execute plan forward FT");
#line 626 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
apply_irradiance <<< gridDim , blockDim >>> (d__wave_PUPIL, N_DFT, N_SIDE_LENSLET, src->fwhm);
#line 629 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 606 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 669 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                N_SOURCE);
binning <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                   src->dev_ptr,
				   pixel_scale, d__theta_x, d__theta_y,
				   absolute_pointing,
                                   N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);

++N_FRAME;
#line 401 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 405 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::propagateNoOverlap(source *src) {
  int M = N_PX_IMAGE/N_PX_CAMERA;
  dim3 gridDim;
  dim3 blockDim(16,16);
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  complexAmplitudeNoOverlap <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
                                                       src->wavefront.amplitude,
                                                       src->wavefront.phase,
                                                       src->wavenumber(),
                                                       N_PX_PUPIL, N_PX_IMAGE, M,
                                                       N_SIDE_LENSLET, src->dev_ptr,
                                                       pixel_scale, d__theta_x, d__theta_y,
                                                       absolute_pointing,1);
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 721 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
N_PHOTON_PER_FRAME = N_PHOTON_PER_SECOND_PER_FRAME = src->n_photon()*src->wavefront.M->area;
float intensity_norm;
intensity_norm = src->wavefront.M->area/src->wavefront.M->nnz;
/*
fprintf(stdout,"Source # of photon: %g\n",src->n_photon());
fprintf(stdout,"Wavefront mask area and nnz: %f & %f\n",src->wavefront.M->area,src->wavefront.M->nnz);
fprintf(stdout,"# of photon per second: %g\n",N_PHOTON_PER_SECOND_PER_FRAME);
*/
#line 420 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
gridDim  = dim3(Nb,Nb);
wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, intensity_norm);
gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                N_SOURCE);
binning <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                   src->dev_ptr,
                                   pixel_scale, d__theta_x, d__theta_y,
                                   absolute_pointing,
                                   N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);
++N_FRAME;
}
#line 437 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::propagateNoOverlapBare(source *src) {
  int M = N_PX_IMAGE/N_PX_CAMERA;
  dim3 gridDim;
  dim3 blockDim(16,16);
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  complexAmplitudeNoOverlap <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
					      src->wavefront.amplitude,
					      src->wavefront.phase,
                                              src->wavenumber(),
					      N_PX_PUPIL, N_PX_IMAGE, M,
					      N_SIDE_LENSLET, src->dev_ptr,
                                              pixel_scale, d__theta_x, d__theta_y,
					      absolute_pointing,1);
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 451 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
gridDim  = dim3(Nb,Nb);
wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, 1.0);
gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                N_SOURCE);
binning <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                   src->dev_ptr,
                                   pixel_scale, d__theta_x, d__theta_y,
                                   absolute_pointing,
                                   N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);
++N_FRAME;
}
#line 468 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  void imaging::propagateNoOverlapSPS(source *src, float d, float wavenumber) {
  int M = N_PX_IMAGE/N_PX_CAMERA;
  dim3 gridDim;
  dim3 blockDim(16,16);
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  complexAmplitudeNoOverlapSPS <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
					      src->wavefront.amplitude,
					      src->wavefront.phase,
                                              wavenumber,
					      N_PX_PUPIL, N_PX_IMAGE, M,
					      N_SIDE_LENSLET, src->dev_ptr,
                                              d,1);
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");

#line 721 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
N_PHOTON_PER_FRAME = N_PHOTON_PER_SECOND_PER_FRAME = src->n_photon()*src->wavefront.M->area;
float intensity_norm;
intensity_norm = src->wavefront.M->area/src->wavefront.M->nnz;
/*
fprintf(stdout,"Source # of photon: %g\n",src->n_photon());
fprintf(stdout,"Wavefront mask area and nnz: %f & %f\n",src->wavefront.M->area,src->wavefront.M->nnz);
fprintf(stdout,"# of photon per second: %g\n",N_PHOTON_PER_SECOND_PER_FRAME);
*/

#line 484 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
gridDim  = dim3(Nb,Nb);
wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, intensity_norm);
gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                N_SOURCE);
// fprintf(stdout,"FWHM=%f\n",src->fwhm);
if (src->fwhm>0) {
  INFO("CONVOLUTION!\n");
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  
#line 622 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_INVERSE),
		    "Unable to execute plan forward FT");
#line 626 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
apply_irradiance <<< gridDim , blockDim >>> (d__wave_PUPIL, N_DFT, N_SIDE_LENSLET, src->fwhm);
#line 629 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 496 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
binning <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                   src->dev_ptr,
                                   pixel_scale, d__theta_x, d__theta_y,
                                   absolute_pointing,
                                   N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);
++N_FRAME;
}
#line 507 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::propagateTT7(source *src) {
  
#line 536 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int M = N_PX_IMAGE/N_PX_CAMERA;
/* if ( (pixel_scale>0) && (zenith<0) ) */
/*   set_pointing_direction(0.0,0.0); */
/*
if (pixel_scale>0 ) {
    printf(" pixel scale: %.2E\n",pixel_scale);
float dx, dy, dx_int, dy_int;
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 544 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
printf("ANGULAR OFFSET:\n");
printf(" . x: %3.0f + %4.3f\n",dx_int,dx);
printf(" . y: %3.0f + %4.3f\n",dy_int,dy);
  }
*/
//absolute_pointing = 1;
//fprintf(stdout,"Absolute pointing: %d\n",absolute_pointing);
dim3 gridDim;
dim3 blockDim(16,16);
gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
#line 509 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  complexAmplitudeTT7 <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
                                                 src->wavefront.amplitude,
                                                 src->wavefront.phase,
                                                 src->wavenumber(),
                                                 src->rays.d__piston_mask,
                                                 N_PX_PUPIL, N_PX_IMAGE, M,
                                                 N_SIDE_LENSLET, src->dev_ptr,
                                                 pixel_scale, d__theta_x, d__theta_y,
                                                 absolute_pointing,1);
  
#line 557 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
                    "Unable to execute plan forward FT");
N_PHOTON_PER_FRAME = N_PHOTON_PER_SECOND_PER_FRAME = src->n_photon()*src->wavefront.M->area;
float intensity_norm;
intensity_norm = src->wavefront.M->area/src->wavefront.M->nnz;

fprintf(stdout,"Source # of photon: %g\n",src->n_photon());
fprintf(stdout,"Wavefront mask area and nnz: (%f,%f) & %f\n",src->rays.V.area,src->wavefront.M->area,src->wavefront.M->nnz);
fprintf(stdout,"# of photon per second: %g\n",N_PHOTON_PER_SECOND_PER_FRAME);

int N, Nb;
N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
gridDim  = dim3(Nb,Nb);
wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, intensity_norm);
if (src->fwhm>0) {
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_INVERSE),
                      "Unable to execute plan forward FT");
  apply_irradiance <<< gridDim , blockDim >>> (d__wave_PUPIL, N_DFT, N_SIDE_LENSLET, src->fwhm);
  HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
                      "Unable to execute plan forward FT");
}
gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                N_SOURCE);
binningTT7 <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                   src->dev_ptr,
                                   pixel_scale, d__theta_x, d__theta_y,
                                   absolute_pointing,
                                   N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);

++N_FRAME;
#line 519 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
void imaging::propagateTT7(source *src, int *d__piston_mask) {
  
#line 536 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int M = N_PX_IMAGE/N_PX_CAMERA;
/* if ( (pixel_scale>0) && (zenith<0) ) */
/*   set_pointing_direction(0.0,0.0); */
/*
if (pixel_scale>0 ) {
    printf(" pixel scale: %.2E\n",pixel_scale);
float dx, dy, dx_int, dy_int;
#line 990 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
if (pixel_scale>0) {
  dx  = absolute_pointing*src[iSource].theta_x - theta_x[iSource];
  dx /= pixel_scale;
  dx_int = rintf( dx );
  dx -= dx_int;
  //dx *= M;
  dy  = absolute_pointing*src[iSource].theta_y - theta_y[iSource];
  dy /= pixel_scale;
  dy_int = rintf( dy );
  dy -= dy_int;
  //dy *= M;
}
#line 544 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
printf("ANGULAR OFFSET:\n");
printf(" . x: %3.0f + %4.3f\n",dx_int,dx);
printf(" . y: %3.0f + %4.3f\n",dy_int,dy);
  }
*/
//absolute_pointing = 1;
//fprintf(stdout,"Absolute pointing: %d\n",absolute_pointing);
dim3 gridDim;
dim3 blockDim(16,16);
gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
#line 522 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    complexAmplitudeTT7 <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
                                                   src->wavefront.amplitude,
                                                   src->wavefront.phase,
                                                   src->wavenumber(),
                                                   d__piston_mask,
                                                   N_PX_PUPIL, N_PX_IMAGE, M,
                                                   N_SIDE_LENSLET, src->dev_ptr,
                                                   pixel_scale, d__theta_x, d__theta_y,
                                                   absolute_pointing,1);
  
#line 557 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
                    "Unable to execute plan forward FT");
N_PHOTON_PER_FRAME = N_PHOTON_PER_SECOND_PER_FRAME = src->n_photon()*src->wavefront.M->area;
float intensity_norm;
intensity_norm = src->wavefront.M->area/src->wavefront.M->nnz;

fprintf(stdout,"Source # of photon: %g\n",src->n_photon());
fprintf(stdout,"Wavefront mask area and nnz: (%f,%f) & %f\n",src->rays.V.area,src->wavefront.M->area,src->wavefront.M->nnz);
fprintf(stdout,"# of photon per second: %g\n",N_PHOTON_PER_SECOND_PER_FRAME);

int N, Nb;
N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
gridDim  = dim3(Nb,Nb);
wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, intensity_norm);
if (src->fwhm>0) {
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_INVERSE),
                      "Unable to execute plan forward FT");
  apply_irradiance <<< gridDim , blockDim >>> (d__wave_PUPIL, N_DFT, N_SIDE_LENSLET, src->fwhm);
  HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
                      "Unable to execute plan forward FT");
}
gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                N_SOURCE);
binningTT7 <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                   src->dev_ptr,
                                   pixel_scale, d__theta_x, d__theta_y,
                                   absolute_pointing,
                                   N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);

++N_FRAME;
#line 532 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
      }
#line 758 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::propagateThroughFieldStop(source *src, float field_stop_diam) {
  int M = N_PX_IMAGE/N_PX_CAMERA;
  dim3 gridDim;
  dim3 blockDim(16,16);
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  complexAmplitude <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
					      src->wavefront.amplitude,
					      src->wavefront.phase,
                                              src->wavenumber(),
					      N_PX_PUPIL, N_PX_IMAGE, M,
					      N_SIDE_LENSLET, src->dev_ptr,
                                              pixel_scale, d__theta_x, d__theta_y,
					      absolute_pointing,0);
  
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");


#line 774 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  field_stop_diam *= N_DFT/N_PX_PUPIL;
  gridDim = dim3(1+N_DFT/16,1+N_DFT/16);
  apply_field_stop <<< gridDim , blockDim >>>( d__wave_PUPIL, N_DFT, N_PX_PUPIL,
					       field_stop_diam);
  
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 721 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
N_PHOTON_PER_FRAME = N_PHOTON_PER_SECOND_PER_FRAME = src->n_photon()*src->wavefront.M->area;
float intensity_norm;
intensity_norm = src->wavefront.M->area/src->wavefront.M->nnz;
/*
fprintf(stdout,"Source # of photon: %g\n",src->n_photon());
fprintf(stdout,"Wavefront mask area and nnz: %f & %f\n",src->wavefront.M->area,src->wavefront.M->nnz);
fprintf(stdout,"# of photon per second: %g\n",N_PHOTON_PER_SECOND_PER_FRAME);
*/
#line 598 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
gridDim  = dim3(Nb,Nb);
wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, intensity_norm);
if (src->fwhm>0) {
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  
#line 622 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_INVERSE),
		    "Unable to execute plan forward FT");
#line 626 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
apply_irradiance <<< gridDim , blockDim >>> (d__wave_PUPIL, N_DFT, N_SIDE_LENSLET, src->fwhm);
#line 629 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 606 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 669 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                N_SOURCE);
binning <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                   src->dev_ptr,
				   pixel_scale, d__theta_x, d__theta_y,
				   absolute_pointing,
                                   N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);

++N_FRAME;
#line 779 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 836 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::propagateThroughPyramid(source *src, float alpha) {
  int M = N_PX_IMAGE/N_PX_CAMERA;
  dim3 gridDim;
  dim3 blockDim(16,16);
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  complexAmplitude <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
					      src->wavefront.amplitude,
					      src->wavefront.phase,
                                              src->wavenumber(),
					      N_PX_PUPIL, N_PX_IMAGE, M,
					      N_SIDE_LENSLET, src->dev_ptr,
                                              pixel_scale, d__theta_x, d__theta_y,
					      absolute_pointing,0);
  
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");

#line 851 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
  gridDim = dim3(1+N_DFT/16,1+N_DFT/16, N_SOURCE);
  apply_pyramid <<< gridDim , blockDim >>>( d__wave_PUPIL, N_DFT, N_PX_PUPIL, alpha);

  
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 721 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
N_PHOTON_PER_FRAME = N_PHOTON_PER_SECOND_PER_FRAME = src->n_photon()*src->wavefront.M->area;
float intensity_norm;
intensity_norm = src->wavefront.M->area/src->wavefront.M->nnz;
/*
fprintf(stdout,"Source # of photon: %g\n",src->n_photon());
fprintf(stdout,"Wavefront mask area and nnz: %f & %f\n",src->wavefront.M->area,src->wavefront.M->nnz);
fprintf(stdout,"# of photon per second: %g\n",N_PHOTON_PER_SECOND_PER_FRAME);
*/
#line 598 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
int N, Nb;
N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
Nb = (int) ceilf( sqrtf(N)/16.0 );
gridDim  = dim3(Nb,Nb);
wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, intensity_norm);
if (src->fwhm>0) {
  gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
  
#line 622 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_INVERSE),
		    "Unable to execute plan forward FT");
#line 626 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
apply_irradiance <<< gridDim , blockDim >>> (d__wave_PUPIL, N_DFT, N_SIDE_LENSLET, src->fwhm);
#line 629 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 606 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 669 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                N_SOURCE);
binning <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                   src->dev_ptr,
				   pixel_scale, d__theta_x, d__theta_y,
				   absolute_pointing,
                                   N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);

++N_FRAME;
#line 855 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
}
#line 888 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::propagateThroughModulatedPyramid(source *src,
                                               float modulation,
                                               int modulation_sampling,
                                               float alpha) {
  int M, nTheta, kTheta, N, Nb;
  float theta, intensity_norm;
  nTheta = modulation_sampling;
  //printf("nTheta=%d\n",nTheta);
  M = N_PX_IMAGE/N_PX_CAMERA;
  dim3 gridDim;
  dim3 blockDim(16,16);
  for (kTheta = 0;kTheta<nTheta;kTheta++) {
    gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
    theta = 2*PI*kTheta/nTheta;
    complexAmplitudeWithModulation <<< gridDim , blockDim >>>(d__wave_PUPIL, N_DFT,
                                                              src->wavefront.amplitude,
                                                              src->wavefront.phase,
                                                              src->wavenumber(),
                                                              N_PX_PUPIL, N_PX_IMAGE, M,
                                                              N_SIDE_LENSLET, src->dev_ptr,
                                                              pixel_scale, d__theta_x, d__theta_y,
                                                              absolute_pointing,0,
                                                              modulation,theta);
    
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");

#line 913 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    gridDim = dim3(1+N_DFT/16,1+N_DFT/16,N_SOURCE);
    apply_pyramid <<< gridDim , blockDim >>>( d__wave_PUPIL, N_DFT, N_PX_PUPIL, alpha);

    
#line 681 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 917 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    N_PHOTON_PER_FRAME = N_PHOTON_PER_SECOND_PER_FRAME = src->n_photon()*src->wavefront.M->area;
    intensity_norm = src->wavefront.M->area/src->wavefront.M->nnz/nTheta;
    N = N_DFT * N_DFT * N_LENSLET*N_SOURCE;
    Nb = (int) ceilf( sqrtf(N)/16.0 );
    gridDim  = dim3(Nb,Nb);
    wave2intensity <<< gridDim , blockDim >>>( d__wave_PUPIL, N, intensity_norm);
    if (src->fwhm>0) {
      gridDim  = dim3(1+N_DFT*N_SIDE_LENSLET/16,1+N_DFT*N_SIDE_LENSLET/16,N_SOURCE);
      
#line 622 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_INVERSE),
		    "Unable to execute plan forward FT");
#line 626 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
apply_irradiance <<< gridDim , blockDim >>> (d__wave_PUPIL, N_DFT, N_SIDE_LENSLET, src->fwhm);
#line 629 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
HANDLE_ERROR_CUFFT( cufftExecC2C(plan, d__wave_PUPIL, d__wave_PUPIL, CUFFT_FORWARD),
		    "Unable to execute plan forward FT");
#line 926 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
    }
    gridDim  = dim3(1+N_PX_CAMERA*N_SIDE_LENSLET/16,1+N_PX_CAMERA*N_SIDE_LENSLET/16,
                    N_SOURCE);
    binning <<< gridDim , blockDim >>>(d__frame, d__wave_PUPIL, N_DFT,
                                       src->dev_ptr,
                                       pixel_scale, d__theta_x, d__theta_y,
                                       absolute_pointing,
                                       N_PX_IMAGE, N_PX_CAMERA, N_SIDE_LENSLET);

  }
  ++N_FRAME;
}

#line 1519 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
float imaging::strehl_ratio(imaging *ref)
{
  int idx;
  idx = (N_PX_CAMERA - 1)/2;
  float psf0, ref_psf0, sr;
  HANDLE_ERROR( cudaMemcpy( &psf0, d__frame+idx*(1+N_PX_CAMERA),sizeof(float),
			    cudaMemcpyDeviceToHost ) );
  HANDLE_ERROR( cudaMemcpy( &ref_psf0, ref->d__frame+idx*(1+N_PX_CAMERA),sizeof(float),
			    cudaMemcpyDeviceToHost ) );

  sr = (psf0*ref->N_FRAME)/(ref_psf0*N_FRAME);
  return sr;
}


#line 1556 "/home/ubuntu/projects/crseo/sys/CEO/imaging/imaging.nw"
void imaging::frame2file(const char *filename){
  int n;
  n = N_PX_CAMERA*N_SIDE_LENSLET;
  n *= n;
  printf("n=%d\n",n);
  dev2file(filename,d__frame,n*N_SOURCE,n,N_SOURCE);
}