#line 181 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
#include "atmosphere.h"

#line 2051 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__device__ float sinc_atm(float x) {
  return (x==0) ? 1.0 : sin(x) / (x) ;
}

#line 2384 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void rayTracingKern(
const float* x_PUPIL,const float* y_PUPIL,
float* phase_screen_PUPIL,const int NXY_PUPIL,
const source* src, const int N_SOURCE,
const float r0,
const layer* layers,
float* d__phase_screen_LAYER,
const int N_LAYER, const float tau)
{

  
#line 2359 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int N_L, N_W, k_LAYER, k_SOURCE, ij, ndx, ix, iy, pitch;
float x, y, xi, yi, zi, s, c, t, fs, ft, g, onemt, xP, yP;
float *z;

#line 2396 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 2364 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ix = blockIdx.x * blockDim.x + threadIdx.x;
iy = blockIdx.y * blockDim.y + threadIdx.y;
pitch = gridDim.x * blockDim.x;
ij = iy * pitch + ix;
xP = x_PUPIL[ij];
yP = y_PUPIL[ij];
k_SOURCE = blockIdx.z;
pitch = 0;

#line 2398 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  if ((ij<NXY_PUPIL) && (k_SOURCE<N_SOURCE))
  {
    for (k_LAYER=0;k_LAYER<N_LAYER;k_LAYER++)
    {

      
#line 2291 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
g = 1 - layers[k_LAYER].altitude/src[k_SOURCE].height;
#line 2294 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x = g*xP;
y = g*yP;
#line 2298 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x += src[k_SOURCE].theta_x*layers[k_LAYER].altitude;
y += src[k_SOURCE].theta_y*layers[k_LAYER].altitude;
#line 2304 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sincosf(layers[k_LAYER].wind_direction, &s, &c);
xi = x*c + y*s;
yi = y*c - x*s;
#line 2309 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
xi -= tau*layers[k_LAYER].wind_speed;

#line 2405 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      
#line 2323 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_L = layers[k_LAYER].N_LENGTH_LAYER;
N_W = layers[k_LAYER].N_WIDTH_LAYER;
z = d__phase_screen_LAYER + pitch;
pitch += N_L*N_W;
#line 2331 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
s = (N_L - 1 )
  *(xi - layers[k_LAYER].WIDTH_LAYER*0.5 + layers[k_LAYER].LENGTH_LAYER)/layers[k_LAYER].LENGTH_LAYER;
t = (N_W - 1 )
  *(yi/layers[k_LAYER].WIDTH_LAYER + 0.5);
#line 2337 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*N_W );
#line 2344 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
if (fs==(N_L-1)) { s += 1 - fs; ndx -= N_W; } else { s -= fs; }
if (ft==(N_W-1)) { t += 1 - ft; ndx -= 1; }   else { t -= ft; }
#line 2348 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
onemt = 1 - t;
zi = ( z[ndx]*onemt + z[ndx+1]*t )*(1-s) +
        ( z[ndx+N_W]*onemt + z[ndx+N_W+1]*t )*s;



#line 2407 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	phase_screen_PUPIL[ij+k_SOURCE*NXY_PUPIL] += r0*zi;
    }
  }
}
#line 2413 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void rayTracingKernReg(
const float delta_x, const int N_X,
const float delta_y, const int N_Y,
float* phase_screen_PUPIL,
const source* src, const int N_SOURCE,
const float r0,
const layer* layers,
float* d__phase_screen_LAYER,
const int N_LAYER, const float tau)
{

  
#line 2359 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int N_L, N_W, k_LAYER, k_SOURCE, ij, ndx, ix, iy, pitch;
float x, y, xi, yi, zi, s, c, t, fs, ft, g, onemt, xP, yP;
float *z;

#line 2426 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 
#line 2374 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ix = blockIdx.x * blockDim.x + threadIdx.x;
iy = blockIdx.y * blockDim.y + threadIdx.y;
ij = iy + ix * N_Y;
xP = delta_x*( ix - (N_X-1)*0.5 );
yP = delta_y*( iy - (N_Y-1)*0.5 );
k_SOURCE = blockIdx.z;
pitch = 0;
int NXY_PUPIL = N_X*N_Y;

#line 2428 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
    if ( (ix<N_X) && (iy<N_Y) )
  {
    for (k_LAYER=0;k_LAYER<N_LAYER;k_LAYER++)
    {

      
#line 2291 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
g = 1 - layers[k_LAYER].altitude/src[k_SOURCE].height;
#line 2294 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x = g*xP;
y = g*yP;
#line 2298 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x += src[k_SOURCE].theta_x*layers[k_LAYER].altitude;
y += src[k_SOURCE].theta_y*layers[k_LAYER].altitude;
#line 2304 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sincosf(layers[k_LAYER].wind_direction, &s, &c);
xi = x*c + y*s;
yi = y*c - x*s;
#line 2309 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
xi -= tau*layers[k_LAYER].wind_speed;

#line 2435 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      
#line 2323 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_L = layers[k_LAYER].N_LENGTH_LAYER;
N_W = layers[k_LAYER].N_WIDTH_LAYER;
z = d__phase_screen_LAYER + pitch;
pitch += N_L*N_W;
#line 2331 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
s = (N_L - 1 )
  *(xi - layers[k_LAYER].WIDTH_LAYER*0.5 + layers[k_LAYER].LENGTH_LAYER)/layers[k_LAYER].LENGTH_LAYER;
t = (N_W - 1 )
  *(yi/layers[k_LAYER].WIDTH_LAYER + 0.5);
#line 2337 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*N_W );
#line 2344 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
if (fs==(N_L-1)) { s += 1 - fs; ndx -= N_W; } else { s -= fs; }
if (ft==(N_W-1)) { t += 1 - ft; ndx -= 1; }   else { t -= ft; }
#line 2348 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
onemt = 1 - t;
zi = ( z[ndx]*onemt + z[ndx+1]*t )*(1-s) +
        ( z[ndx+N_W]*onemt + z[ndx+N_W+1]*t )*s;



#line 2437 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	phase_screen_PUPIL[ij+k_SOURCE*NXY_PUPIL] += r0*zi;
    }
  }
}
__global__ void rayTracingKernReg_masked(
const float delta_x, const int N_X,
const float delta_y, const int N_Y,
float* phase_screen_PUPIL, char *mask,
const source* src, const int N_SOURCE,
const float r0,
const layer* layers,
float* d__phase_screen_LAYER,
const int N_LAYER, const float tau)
{

  
#line 2359 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int N_L, N_W, k_LAYER, k_SOURCE, ij, ndx, ix, iy, pitch;
float x, y, xi, yi, zi, s, c, t, fs, ft, g, onemt, xP, yP;
float *z;

#line 2454 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 
#line 2374 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ix = blockIdx.x * blockDim.x + threadIdx.x;
iy = blockIdx.y * blockDim.y + threadIdx.y;
ij = iy + ix * N_Y;
xP = delta_x*( ix - (N_X-1)*0.5 );
yP = delta_y*( iy - (N_Y-1)*0.5 );
k_SOURCE = blockIdx.z;
pitch = 0;
int NXY_PUPIL = N_X*N_Y;

#line 2456 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
    if ( (ix<N_X) && (iy<N_Y) && (mask[ij+k_SOURCE*NXY_PUPIL]) )
  {
    for (k_LAYER=0;k_LAYER<N_LAYER;k_LAYER++)
    {

      
#line 2291 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
g = 1 - layers[k_LAYER].altitude/src[k_SOURCE].height;
#line 2294 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x = g*xP;
y = g*yP;
#line 2298 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x += src[k_SOURCE].theta_x*layers[k_LAYER].altitude;
y += src[k_SOURCE].theta_y*layers[k_LAYER].altitude;
#line 2304 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sincosf(layers[k_LAYER].wind_direction, &s, &c);
xi = x*c + y*s;
yi = y*c - x*s;
#line 2309 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
xi -= tau*layers[k_LAYER].wind_speed;

#line 2463 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      
#line 2323 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_L = layers[k_LAYER].N_LENGTH_LAYER;
N_W = layers[k_LAYER].N_WIDTH_LAYER;
z = d__phase_screen_LAYER + pitch;
pitch += N_L*N_W;
#line 2331 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
s = (N_L - 1 )
  *(xi - layers[k_LAYER].WIDTH_LAYER*0.5 + layers[k_LAYER].LENGTH_LAYER)/layers[k_LAYER].LENGTH_LAYER;
t = (N_W - 1 )
  *(yi/layers[k_LAYER].WIDTH_LAYER + 0.5);
#line 2337 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*N_W );
#line 2344 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
if (fs==(N_L-1)) { s += 1 - fs; ndx -= N_W; } else { s -= fs; }
if (ft==(N_W-1)) { t += 1 - ft; ndx -= 1; }   else { t -= ft; }
#line 2348 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
onemt = 1 - t;
zi = ( z[ndx]*onemt + z[ndx+1]*t )*(1-s) +
        ( z[ndx+N_W]*onemt + z[ndx+N_W+1]*t )*s;



#line 2465 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	phase_screen_PUPIL[ij+k_SOURCE*NXY_PUPIL] += r0*zi;
    }
  }
}

#line 2082 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void circ_centroids_kernel(float *cx, float *cy, int N_xy,
				      const float R,
				      profile *turb, float r0, layer *layers, int N_LAYER,
				      float *zeta1, float *eta1, float *zeta2, float *eta2,
				      source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;

#line 2090 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  i_source = blockIdx.z;
  kl = j * gridDim.x * blockDim.x + i + i_source*N_xy;

  sum = 0;

  if (kl<N_xy)
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 2102 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      float o, so, co, c, _cx_, _cy_, phase_screen;
      c = -2.0/(R*N_xy);
      o = 2*PI*kl/N_xy;
      sincosf(o,&so,&co);
      x_kl0 = R*co;
      y_kl0 = R*so;
      cx[0] = cy[0] = 0.0;

      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 2113 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 2116 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1278 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl *= gl;
y_kl *= gl;
#line 1282 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 2119 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen = 1.4*r0*sum;
      _cx_ = c*co*phase_screen;
      _cy_ = c*so*phase_screen;
      atomicAdd(&cx[0], _cx_);
      atomicAdd(&cy[0], _cy_);
    }
}
#line 2165 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void circ_uplink_centroids_kernel(float *cx, float *cy, int N_xy,
				      const float R,
                                      profile *turb, float r0, layer *layers, int N_LAYER,
				      float *zeta1, float *eta1, float *zeta2, float *eta2,
				      source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;

#line 2173 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  i_source = blockIdx.z;
  kl = j * gridDim.x * blockDim.x + i + i_source*N_xy;

  sum = 0;

  if (kl<N_xy)
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 2185 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      float o, so, co, c, _cx_, _cy_, phase_screen;
      c = -2.0/(R*N_xy);
      o = 2*PI*kl/N_xy;
      sincosf(o,&so,&co);
      x_kl0 = R*co;
      y_kl0 = R*so;
      cx[0] = cy[0] = 0.0;

      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 2196 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 2199 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1282 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 2202 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
                  sum_l *= gl;
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen = 1.4*r0*sum;
      _cx_ = c*co*phase_screen;
      _cy_ = c*so*phase_screen;
      atomicAdd(&cx[0], _cx_);
      atomicAdd(&cy[0], _cy_);
    }
}
#line 2217 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void circ_focused_uplink_centroids_kernel(float *cx, float *cy, int N_xy,
				      const float R,
                                      profile *turb, float r0, layer *layers, int N_LAYER,
				      float *zeta1, float *eta1, float *zeta2, float *eta2,
				      source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;

#line 2225 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  i_source = blockIdx.z;
  kl = j * gridDim.x * blockDim.x + i + i_source*N_xy;

  sum = 0;

  if (kl<N_xy)
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 2237 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      float o, so, co, c, _cx_, _cy_, phase_screen;
      c = -2.0/(R*N_xy);
      o = 2*PI*kl/N_xy;
      sincosf(o,&so,&co);
      x_kl0 = R*co;
      y_kl0 = R*so;
      cx[0] = cy[0] = 0.0;

      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 2248 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 2251 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1278 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl *= gl;
y_kl *= gl;
#line 1282 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 2254 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
                  sum_l *= gl;
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen = 1.4*r0*sum;
      _cx_ = c*co*phase_screen;
      _cy_ = c*so*phase_screen;
      atomicAdd(&cx[0], _cx_);
      atomicAdd(&cy[0], _cy_);
    }
}

#line 1066 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void setupRandomSequence(curandState *state, int n_state, int seed)
{
  int id;
  id = blockIdx.x * blockDim.x + threadIdx.x;
  if (id<n_state)
    curand_init(seed, id, 0, &state[id]);
}

#line 1075 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void generateVariates(float *zeta1, float *eta1, float *zeta2, float *eta2,
			    curandState *state, int n_state)
{
  int id;
  curandState localState;
  id = blockIdx.x * blockDim.x + threadIdx.x;
  if (id<n_state)
    {
      localState = state[id*4];
      zeta1[id] = curand_uniform(&localState);
      zeta1[id] = sqrtf( -logf( zeta1[id] ) );
      state[id*4] = localState;

      localState = state[id*4+1];
      eta1[id] = curand_uniform(&localState);
      eta1[id] = 2*PI*eta1[id];
      state[id*4+1] = localState;

      localState = state[id*4+2];
      zeta2[id] = curand_uniform(&localState);
      zeta2[id] = sqrtf( -logf( zeta2[id] ) );
      state[id*4+2] = localState;

      localState = state[id*4+3];
      eta2[id] = curand_uniform(&localState);
      eta2[id] = 2*PI*eta2[id];
      state[id*4+3] = localState;

    }
}

#line 1108 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::reset(void)
{
  int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
  dim3 blockDim(64,1);
  dim3 gridDim( 1+nel/64 , 1 );
  generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
}

#line 1171 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void plps(float *phase_screen,
		     float const *x, float const *y, int N_xy,
		     profile *turb, float r0, layer *layers, int N_LAYER,
		     float *zeta1, float *eta1, float *zeta2, float *eta2,
		     source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;

#line 1179 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 1218 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = j * gridDim.x * blockDim.x + i + i_source*N_xy;
x_kl0 = x[kl];
y_kl0 = y[kl];

#line 1181 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  sum = 0;

  if (kl<N_xy)
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 1188 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1191 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 1194 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1278 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl *= gl;
y_kl *= gl;
#line 1282 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 1197 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen[kl] = 1.4*r0*sum;
    }
}

#line 1402 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  __global__ void square_plps(float *phase_screen,
		       float const delta_x, int N_x, float const delta_y, int N_y,
       		       profile *turb, float r0, layer *layers, int N_LAYER,
		       float *zeta1, float *eta1, float *zeta2, float *eta2,
		       source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;

#line 1410 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 1612 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * N_y + j + i_source*N_x*N_y;
x_kl0 = i*delta_x - (N_x-1)*delta_x*0.5;
y_kl0 = j*delta_y - (N_y-1)*delta_y*0.5;

#line 1412 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  sum = 0;
  if ( (i<N_x) && (j<N_y) )
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 1418 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1421 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 1424 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1278 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl *= gl;
y_kl *= gl;
#line 1282 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 1427 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen[kl] += 1.4*r0*sum;
    }
}
#line 1436 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  __global__ void square_plps_layers(float *phase_screen,
		       float const delta_x, int N_x, float const delta_y, int N_y,
       		       profile *turb, float r0, layer *layers, int n_LAYER,
		       float *zeta1, float *eta1, float *zeta2, float *eta2,
		       source *src, float time)
{
  float v;
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;
#line 1444 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 1620 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * N_y + j + i_source*N_x*N_y;
x_kl0 = delta_x*(i - (N_y-1)*0.5 - N_x + N_y);
y_kl0 = delta_y*(j - (N_y-1)*0.5);

#line 1446 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  sum = 0;
  if ( (i<N_x) && (j<N_y) )
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 1452 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1455 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 1458 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		// for (l=0;l<N_LAYER;l++)
		l = n_LAYER;
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1278 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl *= gl;
y_kl *= gl;
#line 1267 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
v = hypotf(layers[l].vx,layers[l].vy);
x_kl += layers[l].altitude*src[i_source].theta_x - v*time;
y_kl += layers[l].altitude*src[i_source].theta_y;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 1462 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen[kl] = 1.4*r0*sum;
    }
}
#line 1471 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void square_plps_masked(float *phase_screen, char *mask,
		       float const delta_x, int N_x, float const delta_y, int N_y,
       		       profile *turb, float r0, layer *layers, int N_LAYER,
		       float *zeta1, float *eta1, float *zeta2, float *eta2,
		       source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;

#line 1479 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 1612 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * N_y + j + i_source*N_x*N_y;
x_kl0 = i*delta_x - (N_x-1)*delta_x*0.5;
y_kl0 = j*delta_y - (N_y-1)*delta_y*0.5;

#line 1481 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  sum = 0;
  if ( ( (i<N_x) && (j<N_y) ) && (mask[kl]) )
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 1487 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1490 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 1493 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1278 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl *= gl;
y_kl *= gl;
#line 1282 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 1496 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen[kl] += 1.4*r0*sum;
    }
}
#line 1535 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void square_gplps(float *phase_screen,
			       float const delta_x, int N_x, float const delta_y, int N_y,
			       profile *turb, float r0, layer *layers, int N_LAYER,
			       float *zeta1, float *eta1, float *zeta2, float *eta2,
			       source *src, float time, float exponent)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;

#line 1543 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 1612 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * N_y + j + i_source*N_x*N_y;
x_kl0 = i*delta_x - (N_x-1)*delta_x*0.5;
y_kl0 = j*delta_y - (N_y-1)*delta_y*0.5;

#line 1545 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  sum = 0;
  if ( (i<N_x) && (j<N_y) )
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 1551 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1603 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -0.5*exponent)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1554 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 1557 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1278 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl *= gl;
y_kl *= gl;
#line 1282 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 1560 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen[kl] = 1.4*r0*sum;
    }
}
#line 1569 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void square_gplps_masked(float *phase_screen, char *mask,
				    float const delta_x, int N_x, float const delta_y, int N_y,
				    profile *turb, float r0, layer *layers, int N_LAYER,
				    float *zeta1, float *eta1, float *zeta2, float *eta2,
				    source *src, float time, float exponent)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1214 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
float sum, sum_l;

#line 1577 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 1612 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * N_y + j + i_source*N_x*N_y;
x_kl0 = i*delta_x - (N_x-1)*delta_x*0.5;
y_kl0 = j*delta_y - (N_y-1)*delta_y*0.5;

#line 1579 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  sum = 0;
  if ( ( (i<N_x) && (j<N_y) ) && (mask[i * N_x + j]) )
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 1585 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1603 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -0.5*exponent)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1588 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1250 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 sum_l = 0;
#line 1591 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 1272 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = x_kl0;
y_kl = y_kl0;
#line 1278 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl *= gl;
y_kl *= gl;
#line 1282 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
#line 1286 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
sum_l += sqrt(layers[l].xi0)*
   ( zeta1[ijl]*cosf( eta1[ijl] +
		     freq_mag*( x_kl*cos_freq_ang + y_kl*sin_freq_ang ) ) +
     zeta2[ijl]*cosf( eta2[ijl] -
		     freq_mag*( x_kl*sin_freq_ang - y_kl*cos_freq_ang ) ));
#line 1594 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum += sqrt_spectrum_kernel*sum_l;
	    }
	}
      phase_screen[kl] = 1.4*r0*sum;
    }
}

#line 1692 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void plps_xy_gradient(float *sx, float *sy,
                                 float *x, float *y, int Nxy,
                                 profile *turb, float wavelength, float r0,
                                 layer *layers, int N_LAYER,
                                 float *zeta1, float *eta1, float *zeta2, float *eta2,
                                 source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1700 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
    float sum_x, sum_l_x, sum_y, sum_l_y, red0, red1, red2;

  i = blockIdx.x * blockDim.x + threadIdx.x;
  kl=i;

  sum_x = sum_y = 0;
  if (i<Nxy)
    {

      x_kl0 = x[i];
      y_kl0 = y[i];

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 1714 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1717 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1720 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		cos_freq_ang *= freq_mag;
		sin_freq_ang *= freq_mag;
		sum_l_x = sum_l_y = 0;
	      for (l=0;l<N_LAYER;l++)
		{
        ijl = ij + l*turb->N_a*turb->N_k;
      gl = 1 - layers[l].altitude/src[i_source].height;
      x_kl = gl*x_kl0;
      y_kl = gl*y_kl0;
      x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
      y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
      red0 = sqrt(layers[l].xi0);
      red1 = zeta1[ijl]*sin( eta1[ijl] + cos_freq_ang*x_kl + sin_freq_ang*y_kl );
      red2 = zeta2[ijl]*sin( eta2[ijl] - sin_freq_ang*x_kl + cos_freq_ang*y_kl );
      sum_l_x += red0*( red1*cos_freq_ang - red2*sin_freq_ang);
      sum_l_y += red0*( red1*sin_freq_ang + red2*cos_freq_ang);
		}
	      sum_x += sqrt_spectrum_kernel*sum_l_x;
	      sum_y += sqrt_spectrum_kernel*sum_l_y;
	    }
	}
      red0 = 1.4*powf(r0,-5.0/6.0)*wavelength/(2*PI);
      sx[kl] = red0*sum_x;
      sy[kl] = red0*sum_y;
    }
}
#line 1835 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void plps_gradient(float *sx, float *sy, int NL, float const d,
			      profile *turb, float wavelength, float r0,
			      layer *layers, int N_LAYER,
			      float *zeta1, float *eta1, float *zeta2, float *eta2,
			      source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1842 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
    float sum_x, sum_l_x, sum_y, sum_l_y, red0, red1, red2;

  
#line 2029 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * NL + j + 2*i_source*NL*NL;
x_kl0 = (i+0.5)*d - NL*d*0.5;
y_kl0 = (j+0.5)*d - NL*d*0.5;

#line 1846 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  sum_x = sum_y = 0;
  if ( (i<NL) && (j<NL) )
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;

#line 1852 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1855 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1858 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		cos_freq_ang *= freq_mag;
		sin_freq_ang *= freq_mag;
		sum_l_x = sum_l_y = 0;
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 2038 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = gl*x_kl0;
y_kl = gl*y_kl0;
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
red0 = sqrt(layers[l].xi0)*
  sinc_atm(0.5*gl*cos_freq_ang*d)*sinc_atm(0.5*gl*sin_freq_ang*d);
red1 = zeta1[ijl]*sin( eta1[ijl] + cos_freq_ang*x_kl + sin_freq_ang*y_kl );
red2 = zeta2[ijl]*sin( eta2[ijl] - sin_freq_ang*x_kl + cos_freq_ang*y_kl );
sum_l_x += red0*( red1*cos_freq_ang - red2*sin_freq_ang);
sum_l_y += red0*( red1*sin_freq_ang + red2*cos_freq_ang);
#line 1864 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum_x += sqrt_spectrum_kernel*sum_l_x;
	      sum_y += sqrt_spectrum_kernel*sum_l_y;
	    }
	}
      red0 = 1.4*powf(r0,-5.0/6.0)*wavelength/(2*PI);
      sx[kl] = red0*sum_x;
      sy[kl] = red0*sum_y;
    }
}
__global__ void plps_gradient_rolling_shutter(float *sx, float *sy, int NL, float const d,
			      profile *turb, float wavelength, float r0,
			      layer *layers, int N_LAYER,
			      float *zeta1, float *eta1, float *zeta2, float *eta2,
			      source *src, float _time_, float delay)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1881 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  float sum_x, sum_l_x, sum_y, sum_l_y, red0, red1, red2;
  float time, tau;
  tau = delay/NL/NL;

  
#line 2029 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * NL + j + 2*i_source*NL*NL;
x_kl0 = (i+0.5)*d - NL*d*0.5;
y_kl0 = (j+0.5)*d - NL*d*0.5;

#line 1887 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  sum_x = sum_y = 0;
  if ( (i<NL) && (j<NL) )
    {

      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;
#line 1892 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      time = _time_ + (i*NL+j)*tau;

      for (i=0;i<turb->N_k;i++)
	{
	  
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1897 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
          for (j=0;j<turb->N_a;j++)
	    {
	      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1900 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		cos_freq_ang *= freq_mag;
		sin_freq_ang *= freq_mag;
		sum_l_x = sum_l_y = 0;
	      for (l=0;l<N_LAYER;l++)
		{
		  
#line 2038 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = gl*x_kl0;
y_kl = gl*y_kl0;
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
red0 = sqrt(layers[l].xi0)*
  sinc_atm(0.5*gl*cos_freq_ang*d)*sinc_atm(0.5*gl*sin_freq_ang*d);
red1 = zeta1[ijl]*sin( eta1[ijl] + cos_freq_ang*x_kl + sin_freq_ang*y_kl );
red2 = zeta2[ijl]*sin( eta2[ijl] - sin_freq_ang*x_kl + cos_freq_ang*y_kl );
sum_l_x += red0*( red1*cos_freq_ang - red2*sin_freq_ang);
sum_l_y += red0*( red1*sin_freq_ang + red2*cos_freq_ang);
#line 1906 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		}
	      sum_x += sqrt_spectrum_kernel*sum_l_x;
	      sum_y += sqrt_spectrum_kernel*sum_l_y;
	    }
	}
      red0 = 1.4*powf(r0,-5.0/6.0)*wavelength/(2*PI);
      sx[kl] = red0*sum_x;
      sy[kl] = red0*sum_y;
    }
}
#line 1933 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
__global__ void plps_gradient_mask(float *sx, float *sy, int NL,
				   char *valid_lenslet, float const d,
				   profile *turb, float wavelength, float r0,
				   layer *layers, int N_LAYER,
				   float *zeta1, float *eta1, float *zeta2, float *eta2,
				   source *src, float time)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1941 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
    float sum_x, sum_l_x, sum_y, sum_l_y, red0, red1, red2;

  
#line 2029 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * NL + j + 2*i_source*NL*NL;
x_kl0 = (i+0.5)*d - NL*d*0.5;
y_kl0 = (j+0.5)*d - NL*d*0.5;

#line 1945 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      sum_x = sum_y = 0;
  if ( (i<NL) && (j<NL) )
    {
      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;
#line 1949 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
	if (valid_lenslet[i*NL+j]>0)
	  {

	    for (i=0;i<turb->N_k;i++)
	      {
		
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 1955 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		  for (j=0;j<turb->N_a;j++)
		    {
		      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 1958 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
			cos_freq_ang *= freq_mag;
		      sin_freq_ang *= freq_mag;
		      sum_l_x = sum_l_y = 0;
		      for (l=0;l<N_LAYER;l++)
			{
			  
#line 2038 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = gl*x_kl0;
y_kl = gl*y_kl0;
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
red0 = sqrt(layers[l].xi0)*
  sinc_atm(0.5*gl*cos_freq_ang*d)*sinc_atm(0.5*gl*sin_freq_ang*d);
red1 = zeta1[ijl]*sin( eta1[ijl] + cos_freq_ang*x_kl + sin_freq_ang*y_kl );
red2 = zeta2[ijl]*sin( eta2[ijl] - sin_freq_ang*x_kl + cos_freq_ang*y_kl );
sum_l_x += red0*( red1*cos_freq_ang - red2*sin_freq_ang);
sum_l_y += red0*( red1*sin_freq_ang + red2*cos_freq_ang);
#line 1964 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
			    }
		      sum_x += sqrt_spectrum_kernel*sum_l_x;
		      sum_y += sqrt_spectrum_kernel*sum_l_y;
		    }
	      }
	    red0 = 1.4*powf(r0,-5.0/6.0)*wavelength/(2*PI);
	    sx[kl] = red0*sum_x;
	    sy[kl] = red0*sum_y;
	  } else {
	  sx[kl] = 0.0;
	  sy[kl] = 0.0;
	}
    }
}
__global__ void plps_gradient_mask_rolling_shutter(float *sx, float *sy, int NL,
				   char *valid_lenslet, float const d,
				   profile *turb, float wavelength, float r0,
				   layer *layers, int N_LAYER,
				   float *zeta1, float *eta1, float *zeta2, float *eta2,
				   source *src, float _time_, float delay)
{
  
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  int i, j, l, ij, ijl, kl, i_source;
float freq_mag0, delta_freq_mag0, f_red0, f_red,
  freq_L0_square, x_kl, y_kl, x_kl0, y_kl0, gl,
  freq_mag, delta_freq_mag, sqrt_spectrum_kernel,
  freq_ang, cos_freq_ang, sin_freq_ang;
#line 1986 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  float sum_x, sum_l_x, sum_y, sum_l_y, red0, red1, red2;
  float time, tau;
  tau = delay/NL/NL;

  
#line 2029 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
i_source = blockIdx.z;
kl = i * NL + j + 2*i_source*NL*NL;
x_kl0 = (i+0.5)*d - NL*d*0.5;
y_kl0 = (j+0.5)*d - NL*d*0.5;

#line 1992 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      sum_x = sum_y = 0;
  if ( (i<NL) && (j<NL) )
    {
      
#line 1227 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red0 = powf(turb->f,1/turb->N_k);
freq_mag0 = 0.5*turb->kmin*( f_red0 + 1 )/f_red0;
delta_freq_mag0 = turb->kmin*( f_red0 - 1 )/f_red0;
f_red  = 1;
freq_L0_square = 2*PI/turb->L0;
freq_L0_square *= freq_L0_square;
#line 1996 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
      time = _time_ + (i*NL+j)*tau;
	if (valid_lenslet[i*NL+j]>0)
	  {

	    for (i=0;i<turb->N_k;i++)
	      {
		
#line 1236 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
f_red *= f_red0;
freq_mag = freq_mag0*f_red;
delta_freq_mag = delta_freq_mag0*f_red;
sqrt_spectrum_kernel =
  powf( freq_mag*freq_mag + freq_L0_square, -11.0/12.0)*
  sqrt(freq_mag*delta_freq_mag*turb->delta);
#line 2003 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
		  for (j=0;j<turb->N_a;j++)
		    {
		      
#line 1245 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
 freq_ang = (j+0.5)*turb->delta;
 sincosf(freq_ang, &sin_freq_ang, &cos_freq_ang);
 ij = i*turb->N_a + j;
#line 2006 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
			cos_freq_ang *= freq_mag;
		      sin_freq_ang *= freq_mag;
		      sum_l_x = sum_l_y = 0;
		      for (l=0;l<N_LAYER;l++)
			{
			  
#line 2038 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
ijl = ij + l*turb->N_a*turb->N_k;
gl = 1 - layers[l].altitude/src[i_source].height;
x_kl = gl*x_kl0;
y_kl = gl*y_kl0;
x_kl += layers[l].altitude*src[i_source].theta_x - layers[l].vx*time;
y_kl += layers[l].altitude*src[i_source].theta_y - layers[l].vy*time;
red0 = sqrt(layers[l].xi0)*
  sinc_atm(0.5*gl*cos_freq_ang*d)*sinc_atm(0.5*gl*sin_freq_ang*d);
red1 = zeta1[ijl]*sin( eta1[ijl] + cos_freq_ang*x_kl + sin_freq_ang*y_kl );
red2 = zeta2[ijl]*sin( eta2[ijl] - sin_freq_ang*x_kl + cos_freq_ang*y_kl );
sum_l_x += red0*( red1*cos_freq_ang - red2*sin_freq_ang);
sum_l_y += red0*( red1*sin_freq_ang + red2*cos_freq_ang);
#line 2012 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
			    }
		      sum_x += sqrt_spectrum_kernel*sum_l_x;
		      sum_y += sqrt_spectrum_kernel*sum_l_y;
		    }
	      }
	    red0 = 1.4*powf(r0,-5.0/6.0)*wavelength/(2*PI);
	    sx[kl] = red0*sum_x;
	    sy[kl] = red0*sum_y;
	  } else {
	  sx[kl] = 0.0;
	  sy[kl] = 0.0;
	}
    }
}

#line 984 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int layer::setup(float _altitude, float _xi0,
		 float _wind_speed, float _wind_direction,
		 float W, int N_W,
		 float field_size,
		 int OSF, float duration)
{
  int N = 0;
  float c, d, tan_wind;
  altitude       = _altitude;
  xi0            = _xi0;
  wind_speed     = _wind_speed;
  wind_direction = _wind_direction;
  vx = wind_speed*cosf(wind_direction);
  vy = wind_speed*sinf(wind_direction);
  if (N_W>0) {
    d = W/(N_W-1);
    d /= OSF;
    tan_wind = tanf(wind_direction);
    c = ( fabsf(tan_wind) + 1.0)/sqrtf(tan_wind*tan_wind + 1.0);
    WIDTH_LAYER    = (W + 2*altitude*tanf(0.5*field_size))*c;
    LENGTH_LAYER   = WIDTH_LAYER + wind_speed*duration;
    N_WIDTH_LAYER  = (int) ceilf( WIDTH_LAYER/d ) + 1;
    N_LENGTH_LAYER = (int) ceilf( LENGTH_LAYER/d ) + 1;
    WIDTH_LAYER    = d*(N_WIDTH_LAYER-1);
    LENGTH_LAYER   = d*(N_LENGTH_LAYER-1);
    N = N_WIDTH_LAYER*N_LENGTH_LAYER;
    INFO("(%5.2fX%5.2f) meters & (%dX%d) pixels ==> %.2fMB \n",
	    LENGTH_LAYER,WIDTH_LAYER,
	    N_LENGTH_LAYER,N_WIDTH_LAYER,N*sizeof(float)/powf(2,20));
  }
  return N;
}
#line 816 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void profile::setup(float *altitude_, float *xi0_,
		    float *wind_speed_, float *wind_direction_, int N_LAYER) {
  altitude       = (float*)malloc(sizeof(float)*N_LAYER);
  xi0            = (float*)malloc(sizeof(float)*N_LAYER);
  wind_speed     = (float*)malloc(sizeof(float)*N_LAYER);
  wind_direction = (float*)malloc(sizeof(float)*N_LAYER);
  for (int k=0;k<N_LAYER;k++)
    {
      altitude[k]       = altitude_[k];
      xi0[k]            = xi0_[k];
      wind_speed[k]     = wind_speed_[k];
      wind_direction[k] = wind_direction_[k];
    }
  l0    = 1e-3;
  L     = 50e3;
  f     = MAX(L,3*L0);
  kmin  = 2*PI/f;
  f     = f/l0;
  delta = 0.05; // TODO: like f, should depends on the input parameters but 5% seems optimal from a statistical stand point
  N_k = logf(f)/log( (2+delta)/(2-delta) );
  N_a = powf( f , 1/N_k ) ;
  N_a = 0.25*PI*( (N_a + 1)/(N_a - 1) );
  N_k = ceilf( N_k );
  N_a = ceilf( N_a );
  f = powf( (4*N_a + PI)/(4*N_a - PI) , N_k );
  delta = 0.5*PI/N_a;
  INFO("\n@(CEO)>profile:\n");
  INFO(" . N_k = %f\n",N_k);
  INFO(" . N_a = %f\n",N_a);
  INFO(" . frequency range      = %8.2f\n",f);
  INFO(" . frequency resolution = %6.4f\n",delta);
  INFO(" . minimum frequency    = %6.4f\n",kmin);
}
#line 852 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void profile::cleanup(void) {
  free(altitude);
  free(xi0);
  free(wind_speed);
  free(wind_direction);
}

#line 595 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_set_id(int _ID_) {
  ID = _ID_;
}
#line 589 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_set_eph(float _EPH_) {
  EPH = _EPH_;
}
#line 460 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::setup(float r0_, float L0, int _NLAYER,
		       float *altitude, float *xi0,
		       float *wind_speed, float *wind_direction) {
  wavelength = 500e-9;
  wavenumber = wavelength/2/PI;
  r0 = r0_;
  N_LAYER = _NLAYER;
  turbulence.L0 = L0;
  turbulence.setup(altitude, xi0,
		   wind_speed, wind_direction, N_LAYER);
  field_size = 0.0;
  N_PHASE_LAYER = 0;
  W = 0.0;
  N_W = 0;
  layers_OSF            = 2;
  layers_duration       = 0.0;
  layers_tau0           = 0.0;
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 784 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 478 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  phase_screen_LAYER = NULL;
  d__phase_screen_LAYER = NULL;
  info();
}
#line 496 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::setup(float r0_, float L0, int _NLAYER,
		       float *altitude, float *xi0,
		       float *wind_speed, float *wind_direction,
                       float _W_, int _N_W_,
		       float _field_size_,
		       float duration) {
  wavelength = 500e-9;
  wavenumber = wavelength/2/PI;
  r0 = r0_;
  N_LAYER = _NLAYER;
  turbulence.L0 = L0;
  turbulence.setup(altitude, xi0,
		   wind_speed, wind_direction, N_LAYER);
  field_size = _field_size_;
  W = _W_;
  N_W = _N_W_;
  layers_OSF            = 2;
  layers_duration       = duration;
  layers_tau0           = 0.0;
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 784 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 516 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  HANDLE_ERROR( cudaMalloc( (void**)&d__phase_screen_LAYER, N_PHASE_LAYER*sizeof(float) ) );
  
#line 1036 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
dim3 blockDimLayers(16,16);
dim3 gridDimLayers;
source layers_src;
layers_src.setup("V",0.0,0.0,INFINITY);
int N = 0, i_LAYER, N_L, __N_W__;
float delta_W = W/(N_W-1);
delta_W /= layers_OSF;
tid.tic();
INFO(" . Computing layer phase screen: ");
for (i_LAYER=0;i_LAYER<N_LAYER;i_LAYER++) {
  INFO("%d,",i_LAYER);
  N_L = layers[i_LAYER].N_LENGTH_LAYER;
  __N_W__ = layers[i_LAYER].N_WIDTH_LAYER;

  gridDimLayers = dim3( 1+N_L/16 , 1+__N_W__/16);
  square_plps_layers <<< gridDimLayers , blockDimLayers >>> (d__phase_screen_LAYER + N,
							     delta_W, N_L,
							     delta_W, __N_W__,
							     d__turbulence,
                   wavenumber,//*powf(r0,-5.0/6.0),
							     d__layers, i_LAYER,
							     zeta1, eta1, zeta2, eta2,
							     layers_src.dev_ptr, layers_tau0);
  N += N_L*__N_W__;
 }
INFO("\b\n");
tid.toc();
#line 518 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
}
#line 530 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::setup(float r0_, float L0, int _NLAYER,
		       float *altitude, float *xi0,
		       float *wind_speed, float *wind_direction,
                       float _W_, int _N_W_,
		       float _field_size_,
		       float duration,
		       const char *fullpath_to_phasescreens,
		       int _N_DURATION) {
  wavelength = 500e-9;
  wavenumber = wavelength/2/PI;
  r0 = r0_;
  N_LAYER = _NLAYER;
  turbulence.L0 = L0;
  turbulence.setup(altitude, xi0,
		   wind_speed, wind_direction, N_LAYER);
  field_size = _field_size_;
  W = _W_;
  N_W = _N_W_;
  layers_OSF            = 2;
  layers_duration       = duration;
  layers_tau0           = 0.0;
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 784 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 552 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 556 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_screen_LAYER, N_PHASE_LAYER*sizeof(float) ) );
N_DURATION = _N_DURATION;
struct stat sb;
if (stat(fullpath_to_phasescreens,&sb)==0) {
  
#line 935 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
FILE *layers_fid;
size_t N_READ;
layers_fid = fopen(fullpath_to_phasescreens,"rb");

N_READ = fread(&N_DURATION,       sizeof(int),  1,layers_fid);
if (N_READ==0)
  ERROR("Cannot read file!");

fclose(layers_fid);

struct stat sb;
int fd;
off_t offset, pa_offset;

INFO("\n@(CEO)>atmosphere: Memory-mapping the phase screens in the layers from file %s...\n",fullpath_to_phasescreens);
tid.tic();
fd = open(fullpath_to_phasescreens, O_RDONLY);
if (fd == -1)
  handle_error("open");

if (fstat(fd, &sb) == -1)           /* To obtain file size */
  handle_error("fstat");
offset = sizeof(int);
pa_offset = offset & ~(sysconf(_SC_PAGE_SIZE) - 1); 
mmap_size = sb.st_size - pa_offset;
//phase_screen_LAYER = (float*)malloc(sizeof(float)*N_PHASE_LAYER*N_DURATION);
//N_READ     = fread(phase_screen_LAYER,sizeof(float),N_PHASE_LAYER*N_DURATION,layers_fid);
//if (N_READ==0)
//  ERROR("Cannot read file!");
phase_screen_LAYER = (float *) mmap(NULL, mmap_size, PROT_READ, MAP_PRIVATE, fd, pa_offset);
if (phase_screen_LAYER == MAP_FAILED)
  handle_error("mmap");

HANDLE_ERROR( cudaMemcpy( d__phase_screen_LAYER,
                          phase_screen_LAYER,
                          sizeof(float)*N_PHASE_LAYER,
                          cudaMemcpyHostToDevice ) );
tid.toc();
#line 561 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
} else {
  save_layer_phasescreens(fullpath_to_phasescreens, N_DURATION);
}
#line 553 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
}
#line 568 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_setup() {
  wavelength = 500e-9;
  wavenumber = wavelength/2/PI;
  r0 = 15e-2;
  turbulence.L0 = 60; 
  ID = 0;
  
#line 600 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_LAYER = 7;
float altitude[7], xi0[7], wind_speed[7], wind_direction[7];
int n_bytes = sizeof(float)*7;
switch(ID) {
  case 2 :
  {
    INFO("@(CEO)>atmosphere: GMT case 2: STRONG GROUND LAYER (Goodwin LCO Jan08, case 7)\n");
    float altitude_2[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_2[]            = {0.2039, 0.1116, 0.1706, 0.2281, 0.1694, 0.0615, 0.054},
    wind_speed_2[]     = {9.4092, 9.6384, 9.7905, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_2[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_2, n_bytes);
    memcpy(xi0, xi0_2, n_bytes);
    memcpy(wind_speed, wind_speed_2, n_bytes);
    memcpy(wind_direction, wind_direction_2, n_bytes);
    break;
  }
  case 3 :
  {
    INFO("@(CEO)>atmosphere: GMT case 3: GOOD CONDITIONS (Goodwin LCO Jan08, case 1)\n");
    float altitude_3[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_3[]            = {0.1410, 0.0381, 0.0542, 0.3403, 0.2528, 0.0917, 0.0819},
    wind_speed_3[]     = {2.1953, 2.2575, 2.3029, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_3[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_3, n_bytes);
    memcpy(xi0, xi0_3, n_bytes);
    memcpy(wind_speed, wind_speed_3, n_bytes);
    memcpy(wind_direction, wind_direction_3, n_bytes);
    break;
  }
  case 4 :
  {
    INFO("@(CEO)>atmosphere: GMT case 4: POOR CONDITIONS (Goodwin LCO Jan08, case 9)\n");
    float altitude_4[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_4[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_4[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_4[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_4, n_bytes);
    memcpy(xi0, xi0_4, n_bytes);
    memcpy(wind_speed, wind_speed_4, n_bytes);
    memcpy(wind_direction, wind_direction_4, n_bytes);
    break;
  }
  case 5 :
  {
    INFO("@(CEO)>atmosphere: GMT case 5: THIN CLOUDS (Goodwin LCO Jan08, case 5)\n");
    float altitude_5[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_5[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_5[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_5[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_5, n_bytes);
    memcpy(xi0, xi0_5, n_bytes);
    memcpy(wind_speed, wind_speed_5, n_bytes);
    memcpy(wind_direction, wind_direction_5, n_bytes);
    break;
  }
  default :
  {
    INFO("@(CEO)>atmosphere: GMT case 1: MEDIAN (Goodwin LCO Jan08, case 5)\n");
    float altitude_1[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_1[]            = {0.1257, 0.0874, 0.0666, 0.3498, 0.2273, 0.0681, 0.0751},
    wind_speed_1[]     = {5.6540, 5.7964, 5.8942, 6.6370, 13.2925, 34.8250, 29.4187},
    wind_direction_1[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_1, n_bytes);
    memcpy(xi0, xi0_1, n_bytes);
    memcpy(wind_speed, wind_speed_1, n_bytes);
    memcpy(wind_direction, wind_direction_1, n_bytes);
  }
}
for (int k=0;k<N_LAYER;k++) {
  altitude[k] -= EPH;
 }
turbulence.setup(altitude, xi0,
                 wind_speed, wind_direction, N_LAYER);
#line 575 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  field_size = 0.0;
  N_PHASE_LAYER = 0;
  W = 0.0;
  N_W = 0;
  layers_OSF            = 2;
  layers_duration       = 0.0;
  layers_tau0           = 0.0;
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 784 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 583 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  phase_screen_LAYER = NULL;
  d__phase_screen_LAYER = NULL;
}
#line 676 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_setup(float r0_, float L0) {
  wavelength = 500e-9;
  wavenumber = wavelength/2/PI;
  r0 = r0_;
  turbulence.L0 = L0;
  
#line 600 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_LAYER = 7;
float altitude[7], xi0[7], wind_speed[7], wind_direction[7];
int n_bytes = sizeof(float)*7;
switch(ID) {
  case 2 :
  {
    INFO("@(CEO)>atmosphere: GMT case 2: STRONG GROUND LAYER (Goodwin LCO Jan08, case 7)\n");
    float altitude_2[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_2[]            = {0.2039, 0.1116, 0.1706, 0.2281, 0.1694, 0.0615, 0.054},
    wind_speed_2[]     = {9.4092, 9.6384, 9.7905, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_2[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_2, n_bytes);
    memcpy(xi0, xi0_2, n_bytes);
    memcpy(wind_speed, wind_speed_2, n_bytes);
    memcpy(wind_direction, wind_direction_2, n_bytes);
    break;
  }
  case 3 :
  {
    INFO("@(CEO)>atmosphere: GMT case 3: GOOD CONDITIONS (Goodwin LCO Jan08, case 1)\n");
    float altitude_3[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_3[]            = {0.1410, 0.0381, 0.0542, 0.3403, 0.2528, 0.0917, 0.0819},
    wind_speed_3[]     = {2.1953, 2.2575, 2.3029, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_3[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_3, n_bytes);
    memcpy(xi0, xi0_3, n_bytes);
    memcpy(wind_speed, wind_speed_3, n_bytes);
    memcpy(wind_direction, wind_direction_3, n_bytes);
    break;
  }
  case 4 :
  {
    INFO("@(CEO)>atmosphere: GMT case 4: POOR CONDITIONS (Goodwin LCO Jan08, case 9)\n");
    float altitude_4[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_4[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_4[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_4[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_4, n_bytes);
    memcpy(xi0, xi0_4, n_bytes);
    memcpy(wind_speed, wind_speed_4, n_bytes);
    memcpy(wind_direction, wind_direction_4, n_bytes);
    break;
  }
  case 5 :
  {
    INFO("@(CEO)>atmosphere: GMT case 5: THIN CLOUDS (Goodwin LCO Jan08, case 5)\n");
    float altitude_5[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_5[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_5[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_5[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_5, n_bytes);
    memcpy(xi0, xi0_5, n_bytes);
    memcpy(wind_speed, wind_speed_5, n_bytes);
    memcpy(wind_direction, wind_direction_5, n_bytes);
    break;
  }
  default :
  {
    INFO("@(CEO)>atmosphere: GMT case 1: MEDIAN (Goodwin LCO Jan08, case 5)\n");
    float altitude_1[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_1[]            = {0.1257, 0.0874, 0.0666, 0.3498, 0.2273, 0.0681, 0.0751},
    wind_speed_1[]     = {5.6540, 5.7964, 5.8942, 6.6370, 13.2925, 34.8250, 29.4187},
    wind_direction_1[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_1, n_bytes);
    memcpy(xi0, xi0_1, n_bytes);
    memcpy(wind_speed, wind_speed_1, n_bytes);
    memcpy(wind_direction, wind_direction_1, n_bytes);
  }
}
for (int k=0;k<N_LAYER;k++) {
  altitude[k] -= EPH;
 }
turbulence.setup(altitude, xi0,
                 wind_speed, wind_direction, N_LAYER);
#line 682 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  field_size = 0.0;
  N_PHASE_LAYER = 0;
  W = 0.0;
  N_W = 0;
  layers_OSF            = 2;
  layers_duration       = 0.0;
  layers_tau0           = 0.0;
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 784 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 690 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  phase_screen_LAYER = NULL;
  d__phase_screen_LAYER = NULL;
}
#line 716 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_setup(float r0_, float L0,
                           float _W_, int _N_W_,
                           float _field_size_,
			   float duration) {
  
#line 727 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
wavelength = 500e-9;
wavenumber = wavelength/2/PI;
r0 = r0_;
turbulence.L0 = L0;
#line 600 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_LAYER = 7;
float altitude[7], xi0[7], wind_speed[7], wind_direction[7];
int n_bytes = sizeof(float)*7;
switch(ID) {
  case 2 :
  {
    INFO("@(CEO)>atmosphere: GMT case 2: STRONG GROUND LAYER (Goodwin LCO Jan08, case 7)\n");
    float altitude_2[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_2[]            = {0.2039, 0.1116, 0.1706, 0.2281, 0.1694, 0.0615, 0.054},
    wind_speed_2[]     = {9.4092, 9.6384, 9.7905, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_2[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_2, n_bytes);
    memcpy(xi0, xi0_2, n_bytes);
    memcpy(wind_speed, wind_speed_2, n_bytes);
    memcpy(wind_direction, wind_direction_2, n_bytes);
    break;
  }
  case 3 :
  {
    INFO("@(CEO)>atmosphere: GMT case 3: GOOD CONDITIONS (Goodwin LCO Jan08, case 1)\n");
    float altitude_3[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_3[]            = {0.1410, 0.0381, 0.0542, 0.3403, 0.2528, 0.0917, 0.0819},
    wind_speed_3[]     = {2.1953, 2.2575, 2.3029, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_3[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_3, n_bytes);
    memcpy(xi0, xi0_3, n_bytes);
    memcpy(wind_speed, wind_speed_3, n_bytes);
    memcpy(wind_direction, wind_direction_3, n_bytes);
    break;
  }
  case 4 :
  {
    INFO("@(CEO)>atmosphere: GMT case 4: POOR CONDITIONS (Goodwin LCO Jan08, case 9)\n");
    float altitude_4[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_4[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_4[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_4[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_4, n_bytes);
    memcpy(xi0, xi0_4, n_bytes);
    memcpy(wind_speed, wind_speed_4, n_bytes);
    memcpy(wind_direction, wind_direction_4, n_bytes);
    break;
  }
  case 5 :
  {
    INFO("@(CEO)>atmosphere: GMT case 5: THIN CLOUDS (Goodwin LCO Jan08, case 5)\n");
    float altitude_5[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_5[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_5[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_5[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_5, n_bytes);
    memcpy(xi0, xi0_5, n_bytes);
    memcpy(wind_speed, wind_speed_5, n_bytes);
    memcpy(wind_direction, wind_direction_5, n_bytes);
    break;
  }
  default :
  {
    INFO("@(CEO)>atmosphere: GMT case 1: MEDIAN (Goodwin LCO Jan08, case 5)\n");
    float altitude_1[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_1[]            = {0.1257, 0.0874, 0.0666, 0.3498, 0.2273, 0.0681, 0.0751},
    wind_speed_1[]     = {5.6540, 5.7964, 5.8942, 6.6370, 13.2925, 34.8250, 29.4187},
    wind_direction_1[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_1, n_bytes);
    memcpy(xi0, xi0_1, n_bytes);
    memcpy(wind_speed, wind_speed_1, n_bytes);
    memcpy(wind_direction, wind_direction_1, n_bytes);
  }
}
for (int k=0;k<N_LAYER;k++) {
  altitude[k] -= EPH;
 }
turbulence.setup(altitude, xi0,
                 wind_speed, wind_direction, N_LAYER);
#line 732 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
field_size = _field_size_;
W = _W_;
N_W = _N_W_;
N_PHASE_LAYER = 0;
layers_OSF            = 2;
layers_duration       = duration;
layers_tau0           = 0.0;
#line 721 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 784 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 722 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  HANDLE_ERROR( cudaMalloc( (void**)&d__phase_screen_LAYER, N_PHASE_LAYER*sizeof(float) ) );
  
#line 1036 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
dim3 blockDimLayers(16,16);
dim3 gridDimLayers;
source layers_src;
layers_src.setup("V",0.0,0.0,INFINITY);
int N = 0, i_LAYER, N_L, __N_W__;
float delta_W = W/(N_W-1);
delta_W /= layers_OSF;
tid.tic();
INFO(" . Computing layer phase screen: ");
for (i_LAYER=0;i_LAYER<N_LAYER;i_LAYER++) {
  INFO("%d,",i_LAYER);
  N_L = layers[i_LAYER].N_LENGTH_LAYER;
  __N_W__ = layers[i_LAYER].N_WIDTH_LAYER;

  gridDimLayers = dim3( 1+N_L/16 , 1+__N_W__/16);
  square_plps_layers <<< gridDimLayers , blockDimLayers >>> (d__phase_screen_LAYER + N,
							     delta_W, N_L,
							     delta_W, __N_W__,
							     d__turbulence,
                   wavenumber,//*powf(r0,-5.0/6.0),
							     d__layers, i_LAYER,
							     zeta1, eta1, zeta2, eta2,
							     layers_src.dev_ptr, layers_tau0);
  N += N_L*__N_W__;
 }
INFO("\b\n");
tid.toc();
#line 724 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
}
#line 756 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_setup(float r0_, float L0,
			     float _W_, int _N_W_,
			     float _field_size_,
			     float duration,
			     const char *fullpath_to_phasescreens,
			     int _N_DURATION) {
  
#line 727 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
wavelength = 500e-9;
wavenumber = wavelength/2/PI;
r0 = r0_;
turbulence.L0 = L0;
#line 600 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_LAYER = 7;
float altitude[7], xi0[7], wind_speed[7], wind_direction[7];
int n_bytes = sizeof(float)*7;
switch(ID) {
  case 2 :
  {
    INFO("@(CEO)>atmosphere: GMT case 2: STRONG GROUND LAYER (Goodwin LCO Jan08, case 7)\n");
    float altitude_2[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_2[]            = {0.2039, 0.1116, 0.1706, 0.2281, 0.1694, 0.0615, 0.054},
    wind_speed_2[]     = {9.4092, 9.6384, 9.7905, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_2[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_2, n_bytes);
    memcpy(xi0, xi0_2, n_bytes);
    memcpy(wind_speed, wind_speed_2, n_bytes);
    memcpy(wind_direction, wind_direction_2, n_bytes);
    break;
  }
  case 3 :
  {
    INFO("@(CEO)>atmosphere: GMT case 3: GOOD CONDITIONS (Goodwin LCO Jan08, case 1)\n");
    float altitude_3[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_3[]            = {0.1410, 0.0381, 0.0542, 0.3403, 0.2528, 0.0917, 0.0819},
    wind_speed_3[]     = {2.1953, 2.2575, 2.3029, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_3[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_3, n_bytes);
    memcpy(xi0, xi0_3, n_bytes);
    memcpy(wind_speed, wind_speed_3, n_bytes);
    memcpy(wind_direction, wind_direction_3, n_bytes);
    break;
  }
  case 4 :
  {
    INFO("@(CEO)>atmosphere: GMT case 4: POOR CONDITIONS (Goodwin LCO Jan08, case 9)\n");
    float altitude_4[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_4[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_4[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_4[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_4, n_bytes);
    memcpy(xi0, xi0_4, n_bytes);
    memcpy(wind_speed, wind_speed_4, n_bytes);
    memcpy(wind_direction, wind_direction_4, n_bytes);
    break;
  }
  case 5 :
  {
    INFO("@(CEO)>atmosphere: GMT case 5: THIN CLOUDS (Goodwin LCO Jan08, case 5)\n");
    float altitude_5[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_5[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_5[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_5[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_5, n_bytes);
    memcpy(xi0, xi0_5, n_bytes);
    memcpy(wind_speed, wind_speed_5, n_bytes);
    memcpy(wind_direction, wind_direction_5, n_bytes);
    break;
  }
  default :
  {
    INFO("@(CEO)>atmosphere: GMT case 1: MEDIAN (Goodwin LCO Jan08, case 5)\n");
    float altitude_1[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_1[]            = {0.1257, 0.0874, 0.0666, 0.3498, 0.2273, 0.0681, 0.0751},
    wind_speed_1[]     = {5.6540, 5.7964, 5.8942, 6.6370, 13.2925, 34.8250, 29.4187},
    wind_direction_1[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_1, n_bytes);
    memcpy(xi0, xi0_1, n_bytes);
    memcpy(wind_speed, wind_speed_1, n_bytes);
    memcpy(wind_direction, wind_direction_1, n_bytes);
  }
}
for (int k=0;k<N_LAYER;k++) {
  altitude[k] -= EPH;
 }
turbulence.setup(altitude, xi0,
                 wind_speed, wind_direction, N_LAYER);
#line 732 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
field_size = _field_size_;
W = _W_;
N_W = _N_W_;
N_PHASE_LAYER = 0;
layers_OSF            = 2;
layers_duration       = duration;
layers_tau0           = 0.0;
#line 763 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 784 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
LOCAL_RAND_SEED = RAND_SEED;
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 764 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 556 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_screen_LAYER, N_PHASE_LAYER*sizeof(float) ) );
N_DURATION = _N_DURATION;
struct stat sb;
if (stat(fullpath_to_phasescreens,&sb)==0) {
  
#line 935 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
FILE *layers_fid;
size_t N_READ;
layers_fid = fopen(fullpath_to_phasescreens,"rb");

N_READ = fread(&N_DURATION,       sizeof(int),  1,layers_fid);
if (N_READ==0)
  ERROR("Cannot read file!");

fclose(layers_fid);

struct stat sb;
int fd;
off_t offset, pa_offset;

INFO("\n@(CEO)>atmosphere: Memory-mapping the phase screens in the layers from file %s...\n",fullpath_to_phasescreens);
tid.tic();
fd = open(fullpath_to_phasescreens, O_RDONLY);
if (fd == -1)
  handle_error("open");

if (fstat(fd, &sb) == -1)           /* To obtain file size */
  handle_error("fstat");
offset = sizeof(int);
pa_offset = offset & ~(sysconf(_SC_PAGE_SIZE) - 1); 
mmap_size = sb.st_size - pa_offset;
//phase_screen_LAYER = (float*)malloc(sizeof(float)*N_PHASE_LAYER*N_DURATION);
//N_READ     = fread(phase_screen_LAYER,sizeof(float),N_PHASE_LAYER*N_DURATION,layers_fid);
//if (N_READ==0)
//  ERROR("Cannot read file!");
phase_screen_LAYER = (float *) mmap(NULL, mmap_size, PROT_READ, MAP_PRIVATE, fd, pa_offset);
if (phase_screen_LAYER == MAP_FAILED)
  handle_error("mmap");

HANDLE_ERROR( cudaMemcpy( d__phase_screen_LAYER,
                          phase_screen_LAYER,
                          sizeof(float)*N_PHASE_LAYER,
                          cudaMemcpyHostToDevice ) );
tid.toc();
#line 561 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
} else {
  save_layer_phasescreens(fullpath_to_phasescreens, N_DURATION);
}
#line 765 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
}
#line 695 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_setup(float r0_, float L0, int _RAND_SEED_) {
  wavelength = 500e-9;
  wavenumber = wavelength/2/PI;
  r0 = r0_;
  turbulence.L0 = L0;
  
#line 600 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_LAYER = 7;
float altitude[7], xi0[7], wind_speed[7], wind_direction[7];
int n_bytes = sizeof(float)*7;
switch(ID) {
  case 2 :
  {
    INFO("@(CEO)>atmosphere: GMT case 2: STRONG GROUND LAYER (Goodwin LCO Jan08, case 7)\n");
    float altitude_2[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_2[]            = {0.2039, 0.1116, 0.1706, 0.2281, 0.1694, 0.0615, 0.054},
    wind_speed_2[]     = {9.4092, 9.6384, 9.7905, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_2[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_2, n_bytes);
    memcpy(xi0, xi0_2, n_bytes);
    memcpy(wind_speed, wind_speed_2, n_bytes);
    memcpy(wind_direction, wind_direction_2, n_bytes);
    break;
  }
  case 3 :
  {
    INFO("@(CEO)>atmosphere: GMT case 3: GOOD CONDITIONS (Goodwin LCO Jan08, case 1)\n");
    float altitude_3[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_3[]            = {0.1410, 0.0381, 0.0542, 0.3403, 0.2528, 0.0917, 0.0819},
    wind_speed_3[]     = {2.1953, 2.2575, 2.3029, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_3[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_3, n_bytes);
    memcpy(xi0, xi0_3, n_bytes);
    memcpy(wind_speed, wind_speed_3, n_bytes);
    memcpy(wind_direction, wind_direction_3, n_bytes);
    break;
  }
  case 4 :
  {
    INFO("@(CEO)>atmosphere: GMT case 4: POOR CONDITIONS (Goodwin LCO Jan08, case 9)\n");
    float altitude_4[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_4[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_4[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_4[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_4, n_bytes);
    memcpy(xi0, xi0_4, n_bytes);
    memcpy(wind_speed, wind_speed_4, n_bytes);
    memcpy(wind_direction, wind_direction_4, n_bytes);
    break;
  }
  case 5 :
  {
    INFO("@(CEO)>atmosphere: GMT case 5: THIN CLOUDS (Goodwin LCO Jan08, case 5)\n");
    float altitude_5[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_5[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_5[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_5[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_5, n_bytes);
    memcpy(xi0, xi0_5, n_bytes);
    memcpy(wind_speed, wind_speed_5, n_bytes);
    memcpy(wind_direction, wind_direction_5, n_bytes);
    break;
  }
  default :
  {
    INFO("@(CEO)>atmosphere: GMT case 1: MEDIAN (Goodwin LCO Jan08, case 5)\n");
    float altitude_1[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_1[]            = {0.1257, 0.0874, 0.0666, 0.3498, 0.2273, 0.0681, 0.0751},
    wind_speed_1[]     = {5.6540, 5.7964, 5.8942, 6.6370, 13.2925, 34.8250, 29.4187},
    wind_direction_1[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_1, n_bytes);
    memcpy(xi0, xi0_1, n_bytes);
    memcpy(wind_speed, wind_speed_1, n_bytes);
    memcpy(wind_direction, wind_direction_1, n_bytes);
  }
}
for (int k=0;k<N_LAYER;k++) {
  altitude[k] -= EPH;
 }
turbulence.setup(altitude, xi0,
                 wind_speed, wind_direction, N_LAYER);
#line 701 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  field_size = 0.0;
  N_PHASE_LAYER = 0;
  W = 0.0;
  N_W = 0;
  layers_OSF            = 2;
  layers_duration       = 0.0;
  layers_tau0           = 0.0;
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 709 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  LOCAL_RAND_SEED = _RAND_SEED_;
  
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 711 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  phase_screen_LAYER = NULL;
  d__phase_screen_LAYER = NULL;
}
#line 741 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_setup(float r0_, float L0,
                           float _W_, int _N_W_,
                           float _field_size_,
			   float duration,
                           int _RAND_SEED_) {
  
#line 727 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
wavelength = 500e-9;
wavenumber = wavelength/2/PI;
r0 = r0_;
turbulence.L0 = L0;
#line 600 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_LAYER = 7;
float altitude[7], xi0[7], wind_speed[7], wind_direction[7];
int n_bytes = sizeof(float)*7;
switch(ID) {
  case 2 :
  {
    INFO("@(CEO)>atmosphere: GMT case 2: STRONG GROUND LAYER (Goodwin LCO Jan08, case 7)\n");
    float altitude_2[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_2[]            = {0.2039, 0.1116, 0.1706, 0.2281, 0.1694, 0.0615, 0.054},
    wind_speed_2[]     = {9.4092, 9.6384, 9.7905, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_2[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_2, n_bytes);
    memcpy(xi0, xi0_2, n_bytes);
    memcpy(wind_speed, wind_speed_2, n_bytes);
    memcpy(wind_direction, wind_direction_2, n_bytes);
    break;
  }
  case 3 :
  {
    INFO("@(CEO)>atmosphere: GMT case 3: GOOD CONDITIONS (Goodwin LCO Jan08, case 1)\n");
    float altitude_3[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_3[]            = {0.1410, 0.0381, 0.0542, 0.3403, 0.2528, 0.0917, 0.0819},
    wind_speed_3[]     = {2.1953, 2.2575, 2.3029, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_3[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_3, n_bytes);
    memcpy(xi0, xi0_3, n_bytes);
    memcpy(wind_speed, wind_speed_3, n_bytes);
    memcpy(wind_direction, wind_direction_3, n_bytes);
    break;
  }
  case 4 :
  {
    INFO("@(CEO)>atmosphere: GMT case 4: POOR CONDITIONS (Goodwin LCO Jan08, case 9)\n");
    float altitude_4[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_4[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_4[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_4[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_4, n_bytes);
    memcpy(xi0, xi0_4, n_bytes);
    memcpy(wind_speed, wind_speed_4, n_bytes);
    memcpy(wind_direction, wind_direction_4, n_bytes);
    break;
  }
  case 5 :
  {
    INFO("@(CEO)>atmosphere: GMT case 5: THIN CLOUDS (Goodwin LCO Jan08, case 5)\n");
    float altitude_5[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_5[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_5[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_5[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_5, n_bytes);
    memcpy(xi0, xi0_5, n_bytes);
    memcpy(wind_speed, wind_speed_5, n_bytes);
    memcpy(wind_direction, wind_direction_5, n_bytes);
    break;
  }
  default :
  {
    INFO("@(CEO)>atmosphere: GMT case 1: MEDIAN (Goodwin LCO Jan08, case 5)\n");
    float altitude_1[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_1[]            = {0.1257, 0.0874, 0.0666, 0.3498, 0.2273, 0.0681, 0.0751},
    wind_speed_1[]     = {5.6540, 5.7964, 5.8942, 6.6370, 13.2925, 34.8250, 29.4187},
    wind_direction_1[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_1, n_bytes);
    memcpy(xi0, xi0_1, n_bytes);
    memcpy(wind_speed, wind_speed_1, n_bytes);
    memcpy(wind_direction, wind_direction_1, n_bytes);
  }
}
for (int k=0;k<N_LAYER;k++) {
  altitude[k] -= EPH;
 }
turbulence.setup(altitude, xi0,
                 wind_speed, wind_direction, N_LAYER);
#line 732 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
field_size = _field_size_;
W = _W_;
N_W = _N_W_;
N_PHASE_LAYER = 0;
layers_OSF            = 2;
layers_duration       = duration;
layers_tau0           = 0.0;
#line 747 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 748 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  LOCAL_RAND_SEED = _RAND_SEED_;
  
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 750 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  HANDLE_ERROR( cudaMalloc( (void**)&d__phase_screen_LAYER, N_PHASE_LAYER*sizeof(float) ) );
  
#line 1036 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
dim3 blockDimLayers(16,16);
dim3 gridDimLayers;
source layers_src;
layers_src.setup("V",0.0,0.0,INFINITY);
int N = 0, i_LAYER, N_L, __N_W__;
float delta_W = W/(N_W-1);
delta_W /= layers_OSF;
tid.tic();
INFO(" . Computing layer phase screen: ");
for (i_LAYER=0;i_LAYER<N_LAYER;i_LAYER++) {
  INFO("%d,",i_LAYER);
  N_L = layers[i_LAYER].N_LENGTH_LAYER;
  __N_W__ = layers[i_LAYER].N_WIDTH_LAYER;

  gridDimLayers = dim3( 1+N_L/16 , 1+__N_W__/16);
  square_plps_layers <<< gridDimLayers , blockDimLayers >>> (d__phase_screen_LAYER + N,
							     delta_W, N_L,
							     delta_W, __N_W__,
							     d__turbulence,
                   wavenumber,//*powf(r0,-5.0/6.0),
							     d__layers, i_LAYER,
							     zeta1, eta1, zeta2, eta2,
							     layers_src.dev_ptr, layers_tau0);
  N += N_L*__N_W__;
 }
INFO("\b\n");
tid.toc();
#line 752 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
}
#line 768 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::gmt_setup(float r0_, float L0,
			     float _W_, int _N_W_,
			     float _field_size_,
			     float duration,
			     const char *fullpath_to_phasescreens,
			     int _N_DURATION,
                             int _RAND_SEED_) {
  
#line 727 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
wavelength = 500e-9;
wavenumber = wavelength/2/PI;
r0 = r0_;
turbulence.L0 = L0;
#line 600 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
N_LAYER = 7;
float altitude[7], xi0[7], wind_speed[7], wind_direction[7];
int n_bytes = sizeof(float)*7;
switch(ID) {
  case 2 :
  {
    INFO("@(CEO)>atmosphere: GMT case 2: STRONG GROUND LAYER (Goodwin LCO Jan08, case 7)\n");
    float altitude_2[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_2[]            = {0.2039, 0.1116, 0.1706, 0.2281, 0.1694, 0.0615, 0.054},
    wind_speed_2[]     = {9.4092, 9.6384, 9.7905, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_2[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_2, n_bytes);
    memcpy(xi0, xi0_2, n_bytes);
    memcpy(wind_speed, wind_speed_2, n_bytes);
    memcpy(wind_direction, wind_direction_2, n_bytes);
    break;
  }
  case 3 :
  {
    INFO("@(CEO)>atmosphere: GMT case 3: GOOD CONDITIONS (Goodwin LCO Jan08, case 1)\n");
    float altitude_3[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_3[]            = {0.1410, 0.0381, 0.0542, 0.3403, 0.2528, 0.0917, 0.0819},
    wind_speed_3[]     = {2.1953, 2.2575, 2.3029, 2.7028, 8.2883, 30.1824, 13.0364},
    wind_direction_3[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_3, n_bytes);
    memcpy(xi0, xi0_3, n_bytes);
    memcpy(wind_speed, wind_speed_3, n_bytes);
    memcpy(wind_direction, wind_direction_3, n_bytes);
    break;
  }
  case 4 :
  {
    INFO("@(CEO)>atmosphere: GMT case 4: POOR CONDITIONS (Goodwin LCO Jan08, case 9)\n");
    float altitude_4[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_4[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_4[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_4[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_4, n_bytes);
    memcpy(xi0, xi0_4, n_bytes);
    memcpy(wind_speed, wind_speed_4, n_bytes);
    memcpy(wind_direction, wind_direction_4, n_bytes);
    break;
  }
  case 5 :
  {
    INFO("@(CEO)>atmosphere: GMT case 5: THIN CLOUDS (Goodwin LCO Jan08, case 5)\n");
    float altitude_5[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_5[]            = {0.1335, 0.0730, 0.1117, 0.3311, 0.2170, 0.0613, 0.0724},
    wind_speed_5[]     = {9.4092, 9.6384, 9.7905, 10.8527, 18.2550, 39.1520, 43.7936},
    wind_direction_5[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_5, n_bytes);
    memcpy(xi0, xi0_5, n_bytes);
    memcpy(wind_speed, wind_speed_5, n_bytes);
    memcpy(wind_direction, wind_direction_5, n_bytes);
    break;
  }
  default :
  {
    INFO("@(CEO)>atmosphere: GMT case 1: MEDIAN (Goodwin LCO Jan08, case 5)\n");
    float altitude_1[]       = {25, 275, 425, 1250, 4000, 8000, 13000},
    xi0_1[]            = {0.1257, 0.0874, 0.0666, 0.3498, 0.2273, 0.0681, 0.0751},
    wind_speed_1[]     = {5.6540, 5.7964, 5.8942, 6.6370, 13.2925, 34.8250, 29.4187},
    wind_direction_1[] = {0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462};
    memcpy(altitude, altitude_1, n_bytes);
    memcpy(xi0, xi0_1, n_bytes);
    memcpy(wind_speed, wind_speed_1, n_bytes);
    memcpy(wind_direction, wind_direction_1, n_bytes);
  }
}
for (int k=0;k<N_LAYER;k++) {
  altitude[k] -= EPH;
 }
turbulence.setup(altitude, xi0,
                 wind_speed, wind_direction, N_LAYER);
#line 732 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
field_size = _field_size_;
W = _W_;
N_W = _N_W_;
N_PHASE_LAYER = 0;
layers_OSF            = 2;
layers_duration       = duration;
layers_tau0           = 0.0;
#line 776 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 788 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
info();
N_DURATION = 0;
#line 1018 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
layers = (layer*) malloc( sizeof(layer)*N_LAYER );
HANDLE_ERROR( cudaMalloc( (void**)&d__layers, sizeof(layer)*N_LAYER ) );
N_PHASE_LAYER = 0;
INFO("\n@(CEO)>atmosphere: __ Layer Phase Screen __\n");
for (int i_LAYER=0; i_LAYER<N_LAYER; i_LAYER++) {
  INFO(" #%d: ",i_LAYER);
  N_PHASE_LAYER += layers[i_LAYER].setup(altitude[i_LAYER], xi0[i_LAYER],
					wind_speed[i_LAYER], wind_direction[i_LAYER],
					 W, N_W, field_size,
					 layers_OSF, layers_duration);
  HANDLE_ERROR( cudaMemcpy( d__layers + i_LAYER, layers + i_LAYER,
			    sizeof(layer),
			    cudaMemcpyHostToDevice ) );
}
if (N_PHASE_LAYER>0)
  INFO(" Total: %.2fMB\n",N_PHASE_LAYER*sizeof(float)/powf(2,20));
#line 791 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__turbulence, sizeof(profile) ) );
HANDLE_ERROR( cudaMemcpy( d__turbulence, &turbulence,
			  sizeof(profile), cudaMemcpyHostToDevice ) );
#line 777 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  LOCAL_RAND_SEED = _RAND_SEED_;
  
#line 796 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
int nel = N_LAYER*turbulence.N_k*turbulence.N_a;
stopwatch tid;
tid.tic();
INFO("\n@(CEO)>atmosphere: initializing %d variates with SEED=%d ...",nel,LOCAL_RAND_SEED);
HANDLE_ERROR( cudaMalloc( (void**)&devStates, 4*nel*sizeof(curandState)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta1, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta1,  nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&zeta2, nel*sizeof(float)) );
HANDLE_ERROR( cudaMalloc( (void**)&eta2,  nel*sizeof(float)) );
dim3 blockDim(64,1);
dim3 gridDim( 1+4*nel/64 , 1 );
setupRandomSequence <<< gridDim,blockDim >>> (devStates, 4*nel, LOCAL_RAND_SEED);
HANDLE_ERROR( cudaDeviceSynchronize() );
generateVariates <<< gridDim,blockDim >>> (zeta1, eta1, zeta2, eta2, devStates, nel);
HANDLE_ERROR( cudaDeviceSynchronize() );
INFO("done!\n");
tid.toc();
#line 779 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 556 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_screen_LAYER, N_PHASE_LAYER*sizeof(float) ) );
N_DURATION = _N_DURATION;
struct stat sb;
if (stat(fullpath_to_phasescreens,&sb)==0) {
  
#line 935 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
FILE *layers_fid;
size_t N_READ;
layers_fid = fopen(fullpath_to_phasescreens,"rb");

N_READ = fread(&N_DURATION,       sizeof(int),  1,layers_fid);
if (N_READ==0)
  ERROR("Cannot read file!");

fclose(layers_fid);

struct stat sb;
int fd;
off_t offset, pa_offset;

INFO("\n@(CEO)>atmosphere: Memory-mapping the phase screens in the layers from file %s...\n",fullpath_to_phasescreens);
tid.tic();
fd = open(fullpath_to_phasescreens, O_RDONLY);
if (fd == -1)
  handle_error("open");

if (fstat(fd, &sb) == -1)           /* To obtain file size */
  handle_error("fstat");
offset = sizeof(int);
pa_offset = offset & ~(sysconf(_SC_PAGE_SIZE) - 1); 
mmap_size = sb.st_size - pa_offset;
//phase_screen_LAYER = (float*)malloc(sizeof(float)*N_PHASE_LAYER*N_DURATION);
//N_READ     = fread(phase_screen_LAYER,sizeof(float),N_PHASE_LAYER*N_DURATION,layers_fid);
//if (N_READ==0)
//  ERROR("Cannot read file!");
phase_screen_LAYER = (float *) mmap(NULL, mmap_size, PROT_READ, MAP_PRIVATE, fd, pa_offset);
if (phase_screen_LAYER == MAP_FAILED)
  handle_error("mmap");

HANDLE_ERROR( cudaMemcpy( d__phase_screen_LAYER,
                          phase_screen_LAYER,
                          sizeof(float)*N_PHASE_LAYER,
                          cudaMemcpyHostToDevice ) );
tid.toc();
#line 561 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
} else {
  save_layer_phasescreens(fullpath_to_phasescreens, N_DURATION);
}
#line 780 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
}

#line 1118 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::cleanup(void) {
  INFO("@(CEO)>atmosphere: freeing memory!\n");
  turbulence.cleanup();
  HANDLE_ERROR( cudaFree( devStates ) );
  HANDLE_ERROR( cudaFree( zeta1  ));
  HANDLE_ERROR( cudaFree( eta1 ) );
  HANDLE_ERROR( cudaFree( zeta2 ) );
  HANDLE_ERROR( cudaFree( eta2 ) );
  if (d__phase_screen_LAYER!=NULL) {
    HANDLE_ERROR( cudaFree( d__phase_screen_LAYER ) );
    free( layers );
    HANDLE_ERROR( cudaFree( d__layers ) );
  }
  if (phase_screen_LAYER!=NULL)
    munmap(phase_screen_LAYER, mmap_size);
  //    free( phase_screen_LAYER );
}

#line 339 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::info(void)
{
INFO("\n\x1B[1;42m@(CEO)>atmosphere: Von Karman atmospheric turbulence model\x1B[;42m\n");
INFO(" . wavelength = %5.2fnm\n . r0 = %5.2fcm\n . L0 = %5.2fm\n",
       wavelength*1e9,r0*1e2,turbulence.L0);
INFO("----------------------------------------------------\n");
INFO("  Layer   Altitude[m]   fr0    wind([m/s] [deg])\n");
for (int kLayer=0;kLayer<N_LAYER;kLayer++) {
  INFO("  %2d      %8.2f      %4.2f    (%5.2f %6.2f)\n",
  kLayer,
  turbulence.altitude[kLayer],
  turbulence.xi0[kLayer],
  turbulence.wind_speed[kLayer],
  turbulence.wind_direction[kLayer]*180/PI);
}
INFO("----------------------------------------------------\x1B[0m\n");
}

#line 861 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::save_layer_phasescreens(const char* _fullpath_to_phasescreens, int _N_DURATION)
{
  int k_DURATION;
  size_t N_WRITE, N_DATA;
  stopwatch tid;
  time_t now;
  struct tm beg;
  char buff[128];

  FILE *fid;
  fid = fopen(_fullpath_to_phasescreens,"wb");
  fwrite(&_N_DURATION,sizeof(int),1,fid);

  N_DATA = N_PHASE_LAYER*_N_DURATION;
  phase_screen_LAYER = (float*)malloc(sizeof(float)*N_DATA);

  for (k_DURATION=0;k_DURATION<_N_DURATION;k_DURATION++) {

    time(&now);
    beg = *localtime(&now);

    INFO("\n * LAYER PHASE SCREENS COMPUTING: %d/%d *\n",k_DURATION+1,_N_DURATION);
    layers_tau0 = k_DURATION*layers_duration;
    
#line 1036 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
dim3 blockDimLayers(16,16);
dim3 gridDimLayers;
source layers_src;
layers_src.setup("V",0.0,0.0,INFINITY);
int N = 0, i_LAYER, N_L, __N_W__;
float delta_W = W/(N_W-1);
delta_W /= layers_OSF;
tid.tic();
INFO(" . Computing layer phase screen: ");
for (i_LAYER=0;i_LAYER<N_LAYER;i_LAYER++) {
  INFO("%d,",i_LAYER);
  N_L = layers[i_LAYER].N_LENGTH_LAYER;
  __N_W__ = layers[i_LAYER].N_WIDTH_LAYER;

  gridDimLayers = dim3( 1+N_L/16 , 1+__N_W__/16);
  square_plps_layers <<< gridDimLayers , blockDimLayers >>> (d__phase_screen_LAYER + N,
							     delta_W, N_L,
							     delta_W, __N_W__,
							     d__turbulence,
                   wavenumber,//*powf(r0,-5.0/6.0),
							     d__layers, i_LAYER,
							     zeta1, eta1, zeta2, eta2,
							     layers_src.dev_ptr, layers_tau0);
  N += N_L*__N_W__;
 }
INFO("\b\n");
tid.toc();

#line 886 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
    HANDLE_ERROR( cudaMemcpy( phase_screen_LAYER + k_DURATION*N_PHASE_LAYER,
			      d__phase_screen_LAYER,
			      sizeof(float)*N_PHASE_LAYER,
			      cudaMemcpyDeviceToHost ) );
    N_WRITE = fwrite(phase_screen_LAYER + k_DURATION*N_PHASE_LAYER,
		     sizeof(float),N_PHASE_LAYER,fid);
    INFO(" ==>> Written %lu elements out of %lu elements to file %s\n",
	    N_WRITE,N_PHASE_LAYER,_fullpath_to_phasescreens);
    if (N_WRITE<N_PHASE_LAYER) { break; }

    beg.tm_sec += tid.elapsedTime*1e-3*(_N_DURATION-k_DURATION);
    mktime(&beg);
    strftime(buff, sizeof buff, "%c", &beg);

    INFO("\n * ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ * ETA: %s\n",buff);
  }
  fclose(fid);

  layers_tau0 = 0.0;
  HANDLE_ERROR( cudaMemcpy( d__phase_screen_LAYER,
			    phase_screen_LAYER,
			    sizeof(float)*N_PHASE_LAYER,
			    cudaMemcpyHostToDevice ) );
}

#line 1139 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen(float *phase_screen,
		       float const *x, float const *y, int N_xy,
		       source *src, float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( ceilf(sqrt(N_xy)/16) , ceilf(sqrt(N_xy)/16));
  plps <<< gridDim , blockDim >>> (phase_screen,
				   x, y, N_xy,
				   d__turbulence,
				   wavenumber*powf(r0,-5.0/6.0),
				   d__layers, N_LAYER,
				   zeta1, eta1, zeta2, eta2,
				   src->dev_ptr, time);
}
#line 1156 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen(float const *x, float const *y, int N_xy,
		       source *src, float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( ceilf(sqrt(N_xy)/16) , ceilf(sqrt(N_xy)/16));
  plps <<< gridDim , blockDim >>> (src->wavefront.phase,
				   x, y, N_xy,
				   d__turbulence,
				   wavenumber*powf(r0,-5.0/6.0),
				   d__layers, N_LAYER,
				   zeta1, eta1, zeta2, eta2,
				   src->dev_ptr, time);
}
#line 1301 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen(float const delta_x, int N_x, float const delta_y, int N_y,
		       source *src, float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+N_x/16 , 1+N_y/16);
  square_plps <<< gridDim , blockDim >>> (src->wavefront.phase,
				   delta_x, N_x, delta_y, N_y,
				   d__turbulence,
				   wavenumber*powf(r0,-5.0/6.0),
				   d__layers, N_LAYER,
				   zeta1, eta1, zeta2, eta2,
				   src->dev_ptr, time);
}
#line 1315 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen(float *phase_screen,
                       float const delta_x, int N_x, float const delta_y, int N_y,
		       source *src, float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+N_x/16 , 1+N_y/16);
  square_plps <<< gridDim , blockDim >>> (phase_screen,
				   delta_x, N_x, delta_y, N_y,
				   d__turbulence,
				   wavenumber*powf(r0,-5.0/6.0),
				   d__layers, N_LAYER,
				   zeta1, eta1, zeta2, eta2,
				   src->dev_ptr, time);
}
#line 1330 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen(source *src,
                       float const delta_x, int N_x, float const delta_y, int N_y,
		       float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+N_x/16 , 1+N_y/16, src->N_SRC);
  if (src->wavefront.M==NULL) {
    square_plps <<< gridDim , blockDim >>>
      (src->wavefront.phase,
       delta_x, N_x, delta_y, N_y,
       d__turbulence,
       wavenumber*powf(r0,-5.0/6.0),
       d__layers, N_LAYER,
       zeta1, eta1, zeta2, eta2,
       src->dev_ptr, time);
  } else {
    square_plps_masked <<< gridDim , blockDim >>>
      (src->wavefront.phase, src->wavefront.M->m,
       delta_x, N_x, delta_y, N_y,
       d__turbulence,
       wavenumber*powf(r0,-5.0/6.0),
       d__layers, N_LAYER,
       zeta1, eta1, zeta2, eta2,
       src->dev_ptr, time);
  }
}
#line 1357 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen(source *src, float time)
{
  complex_amplitude *wavefront;
  wavefront = &(src->wavefront);
  dim3 blockDim(16,16);
  dim3 gridDim( 1+wavefront->M->size_px[0]/16 ,
		1+wavefront->M->size_px[1]/16, src->N_SRC);
  if (wavefront->M==NULL)
    square_plps <<< gridDim , blockDim >>>
      (wavefront->phase,
       wavefront->M->delta[0], wavefront->M->size_px[0],
       wavefront->M->delta[1], wavefront->M->size_px[1],
       d__turbulence,
       wavenumber*powf(r0,-5.0/6.0),
       d__layers, N_LAYER,
       zeta1, eta1, zeta2, eta2,
       src->dev_ptr, time);
  else
    square_plps_masked <<< gridDim , blockDim >>>
      (wavefront->phase, wavefront->M->m,
       wavefront->M->delta[0], wavefront->M->size_px[0],
       wavefront->M->delta[1], wavefront->M->size_px[1],
       d__turbulence,
       wavenumber*powf(r0,-5.0/6.0),
       d__layers, N_LAYER,
       zeta1, eta1, zeta2, eta2,
       src->dev_ptr, time);
}
#line 1386 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen(source *src, int N_SRC,
                       float const delta_x, int N_x, float const delta_y, int N_y,
		       float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+N_x/16 , 1+N_y/16, N_SRC);
  square_plps <<< gridDim , blockDim >>> (src->wavefront.phase,
				   delta_x, N_x, delta_y, N_y,
				   d__turbulence,
				   wavenumber*powf(r0,-5.0/6.0),
                                   d__layers, N_LAYER,
				   zeta1, eta1, zeta2, eta2,
				   src->dev_ptr, time);
}
#line 1506 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen(source *src,
				  float const delta_x, int N_x,
				  float const delta_y, int N_y,
				  float time, float exponent)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+N_x/16 , 1+N_y/16, src->N_SRC);
  if (src->wavefront.M==NULL) {
    square_gplps <<< gridDim , blockDim >>>
      (src->wavefront.phase,
       delta_x, N_x, delta_y, N_y,
       d__turbulence,
       wavenumber*powf(r0,-5.0/6.0),
       d__layers, N_LAYER,
       zeta1, eta1, zeta2, eta2,
       src->dev_ptr, time, exponent);
  } else {
    square_gplps_masked <<< gridDim , blockDim >>>
      (src->wavefront.phase, src->wavefront.M->m,
       delta_x, N_x, delta_y, N_y,
       d__turbulence,
       wavenumber*powf(r0,-5.0/6.0),
       d__layers, N_LAYER,
       zeta1, eta1, zeta2, eta2,
       src->dev_ptr, time, exponent);
  }
}
#line 240 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
  
#line 1677 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen_gradient(float *sx, float *sy,
                                           float *x, float *y, int Nxy,
                                           source *src, float time)
{
  dim3 blockDim(256);
  dim3 gridDim( 1+Nxy/256);
  plps_xy_gradient <<< gridDim , blockDim >>> (sx, sy,
                                               x,y,Nxy,
                                               d__turbulence, wavelength, r0,
                                               d__layers, N_LAYER,
                                               zeta1, eta1, zeta2, eta2,
                                               src->dev_ptr, time);
}
#line 1750 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen_gradient(float *sx, float *sy, int NL,
                                           float const d, source *src, float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+NL/16 , 1+NL/16);
  plps_gradient <<< gridDim , blockDim >>> (sx, sy, NL, d,
					    d__turbulence, wavelength, r0,
					    d__layers, N_LAYER,
					    zeta1, eta1, zeta2, eta2,
					    src->dev_ptr, time);
}
#line 1762 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen_gradient(centroiding *cog, int NL, float const d,
					   source *src, float time)
{
  dim3 blockDim(16,16);
    dim3 gridDim( 1+NL/16 , 1+NL/16, src->N_SRC);
  if (cog->MASK_SET)
    plps_gradient_mask <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy, NL,
						   cog->lenslet_mask, d,
						   d__turbulence, wavelength, r0,
						   d__layers, N_LAYER,
						   zeta1, eta1, zeta2, eta2,
						   src->dev_ptr, time);
  else
    plps_gradient <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy, NL, d,
					      d__turbulence, wavelength, r0,
					      d__layers, N_LAYER,
					      zeta1, eta1, zeta2, eta2,
					      src->dev_ptr, time);
}
void atmosphere::get_phase_screen_gradient_rolling_shutter(centroiding *cog, int NL, float const d,
							   source *src, float time, float delay)
{
  dim3 blockDim(16,16);
    dim3 gridDim( 1+NL/16 , 1+NL/16, src->N_SRC);
  if (cog->MASK_SET)
    plps_gradient_mask_rolling_shutter <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy, NL,
								   cog->lenslet_mask, d,
								   d__turbulence, wavelength, r0,
								   d__layers, N_LAYER,
								   zeta1, eta1, zeta2, eta2,
								   src->dev_ptr, time, delay);
  else
    plps_gradient_rolling_shutter <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy, NL, d,
							      d__turbulence, wavelength, r0,
							      d__layers, N_LAYER,
							      zeta1, eta1, zeta2, eta2,
							      src->dev_ptr, time,delay);
}
#line 1801 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen_gradient(centroiding *cog, int NL, float const d,
					   source *src, int N_SRC, float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+NL/16 , 1+NL/16, N_SRC);
  if (cog->MASK_SET)
    plps_gradient_mask <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy, NL,
						   cog->lenslet_mask, d,
						   d__turbulence, wavelength, r0,
						   d__layers, N_LAYER,
						   zeta1, eta1, zeta2, eta2,
						   src->dev_ptr, time);
  else
    plps_gradient <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy, NL, d,
					      d__turbulence, wavelength, r0,
					      d__layers, N_LAYER,
					      zeta1, eta1, zeta2, eta2,
					      src->dev_ptr, time);
}
#line 1821 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen_gradient(centroiding *cog, int NL, char *valid_lenslet,
                                           float const d, source *src, int N_SRC, float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+NL/16 , 1+NL/16, N_SRC);
  plps_gradient_mask <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy, NL,
						 valid_lenslet, d,
						 d__turbulence, wavelength, r0,
						 d__layers, N_LAYER,
						 zeta1, eta1, zeta2, eta2,
						 src->dev_ptr, time);
}
#line 1920 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen_gradient(float *sx, float *sy, int NL, char *valid_lenslet,
					   float const d, source *src, float time)
{
  dim3 blockDim(16,16);
  dim3 gridDim( 1+NL/16 , 1+NL/16);
  plps_gradient_mask <<< gridDim , blockDim >>> (sx, sy, NL, valid_lenslet, d,
						 d__turbulence, wavelength, r0,
						 d__layers, N_LAYER,
						 zeta1, eta1, zeta2, eta2,
						 src->dev_ptr, time);
}
#line 2068 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen_circ_centroids(centroiding *cog, const float R,
                                                 source *src, int N_SRC, float time)
{
  int N_o = 360;
  dim3 blockDim(16,16);
  dim3 gridDim( ceilf(sqrt(N_o)/16) , ceilf(sqrt(N_o)/16));
  circ_centroids_kernel <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy,  N_o, R,
						    d__turbulence,
						    wavenumber*powf(r0,-5.0/6.0),
						    d__layers, N_LAYER,
						    zeta1, eta1, zeta2, eta2,
						    src->dev_ptr, time);
}
#line 2141 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::get_phase_screen_circ_uplink_centroids(centroiding *cog, const float R,
                                                 source *src, int N_SRC, float time, char focused)
{
  int N_o = 360;
  dim3 blockDim(16,16);
  dim3 gridDim( ceilf(sqrt(N_o)/16) , ceilf(sqrt(N_o)/16));
  if (focused==1)
    circ_focused_uplink_centroids_kernel <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy,  N_o, R,
								     d__turbulence,
								     wavenumber*powf(r0,-5.0/6.0),
								     d__layers, N_LAYER,
								     zeta1, eta1, zeta2, eta2,
								     src->dev_ptr, time);
  else
    circ_uplink_centroids_kernel <<< gridDim , blockDim >>> (cog->d__cx, cog->d__cy,  N_o, R,
							     d__turbulence,
							     wavenumber*powf(r0,-5.0/6.0),
							     d__layers, N_LAYER,
							     zeta1, eta1, zeta2, eta2,
							     src->dev_ptr, time);
}

#line 2473 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::rayTracing(
const float* x_PUPIL,const float* y_PUPIL,
float* phase_screen_PUPIL, const int NXY_PUPIL,
const source *src, const float tau)
{
  int k_DURATION;
  float new_tau = tau - layers_tau0;
  stopwatch tid;

  //INFO(" . tau=%f\n . layers_tau=%f\n . layers_duration=%f\n",tau,layers_tau0,layers_duration);

  if ( (fabs(tau-layers_tau0)>layers_duration) || (new_tau<0) ) {

    k_DURATION = (int) (tau/layers_duration);
    layers_tau0 = k_DURATION*layers_duration;
    new_tau     = tau - k_DURATION*layers_duration;

    if (N_DURATION>0) {
      INFO("\n@(CEO)>atmosphere: Loading phase screens (%.0fs) to device...\n",layers_tau0);
      tid.tic();
      HANDLE_ERROR( cudaMemcpy( d__phase_screen_LAYER,
				phase_screen_LAYER + k_DURATION*N_PHASE_LAYER,
				sizeof(float)*N_PHASE_LAYER,
				cudaMemcpyHostToDevice ) );
      tid.toc();
    } else {
      
#line 1036 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
dim3 blockDimLayers(16,16);
dim3 gridDimLayers;
source layers_src;
layers_src.setup("V",0.0,0.0,INFINITY);
int N = 0, i_LAYER, N_L, __N_W__;
float delta_W = W/(N_W-1);
delta_W /= layers_OSF;
tid.tic();
INFO(" . Computing layer phase screen: ");
for (i_LAYER=0;i_LAYER<N_LAYER;i_LAYER++) {
  INFO("%d,",i_LAYER);
  N_L = layers[i_LAYER].N_LENGTH_LAYER;
  __N_W__ = layers[i_LAYER].N_WIDTH_LAYER;

  gridDimLayers = dim3( 1+N_L/16 , 1+__N_W__/16);
  square_plps_layers <<< gridDimLayers , blockDimLayers >>> (d__phase_screen_LAYER + N,
							     delta_W, N_L,
							     delta_W, __N_W__,
							     d__turbulence,
                   wavenumber,//*powf(r0,-5.0/6.0),
							     d__layers, i_LAYER,
							     zeta1, eta1, zeta2, eta2,
							     layers_src.dev_ptr, layers_tau0);
  N += N_L*__N_W__;
 }
INFO("\b\n");
tid.toc();
#line 2500 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
    }

  }
  dim3 blockDim(16,16);
  dim3 gridDim( ceilf(sqrtf(NXY_PUPIL)/16) , ceilf(sqrtf(NXY_PUPIL)/16), src->N_SRC);
  rayTracingKern <<< gridDim, blockDim >>>(
  x_PUPIL, y_PUPIL,
  phase_screen_PUPIL, NXY_PUPIL,
  src->dev_ptr, src->N_SRC,
  powf(r0,-5.0/6.0),
  d__layers,
  d__phase_screen_LAYER,
  N_LAYER, tau-layers_tau0);
}
#line 2524 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
void atmosphere::rayTracing(source *src,
			    const float delta_x, const int N_X,
			    const float delta_y, const int N_Y,
			    const float tau)
{
  int k_DURATION;
  float new_tau = tau - layers_tau0;
  stopwatch tid;

  //INFO(" . tau=%f\n . layers_tau=%f\n . layers_duration=%f\n",tau,layers_tau0,layers_duration);

  if ( (fabs(tau-layers_tau0)>layers_duration) || (new_tau<0) ) {

    k_DURATION = (int) (tau/layers_duration);
    layers_tau0 = k_DURATION*layers_duration;
    new_tau     = tau - k_DURATION*layers_duration;

    if (N_DURATION>0) {
      INFO("\n@(CEO)>atmosphere: Loading phase screens (%.0fs) to device...\n",layers_tau0);
      tid.tic();
      HANDLE_ERROR( cudaMemcpy( d__phase_screen_LAYER,
				phase_screen_LAYER + k_DURATION*N_PHASE_LAYER,
				sizeof(float)*N_PHASE_LAYER,
				cudaMemcpyHostToDevice ) );
      tid.toc();
    } else {
      
#line 1036 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
dim3 blockDimLayers(16,16);
dim3 gridDimLayers;
source layers_src;
layers_src.setup("V",0.0,0.0,INFINITY);
int N = 0, i_LAYER, N_L, __N_W__;
float delta_W = W/(N_W-1);
delta_W /= layers_OSF;
tid.tic();
INFO(" . Computing layer phase screen: ");
for (i_LAYER=0;i_LAYER<N_LAYER;i_LAYER++) {
  INFO("%d,",i_LAYER);
  N_L = layers[i_LAYER].N_LENGTH_LAYER;
  __N_W__ = layers[i_LAYER].N_WIDTH_LAYER;

  gridDimLayers = dim3( 1+N_L/16 , 1+__N_W__/16);
  square_plps_layers <<< gridDimLayers , blockDimLayers >>> (d__phase_screen_LAYER + N,
							     delta_W, N_L,
							     delta_W, __N_W__,
							     d__turbulence,
                   wavenumber,//*powf(r0,-5.0/6.0),
							     d__layers, i_LAYER,
							     zeta1, eta1, zeta2, eta2,
							     layers_src.dev_ptr, layers_tau0);
  N += N_L*__N_W__;
 }
INFO("\b\n");
tid.toc();
#line 2551 "/home/ubuntu/projects/crseo/sys/CEO/atmosphere/atmosphere.nw"
    }

  }

  dim3 blockDim(16,16);
  dim3 gridDim( 1+N_X/16 , 1+N_Y/16, src->N_SRC);
  if (src->wavefront.M==NULL) {
    rayTracingKernReg <<< gridDim, blockDim >>>(delta_x, N_X,
						delta_y, N_Y,
						src->wavefront.phase,
						src->dev_ptr, src->N_SRC,
						powf(r0,-5.0/6.0),
						d__layers,
						d__phase_screen_LAYER,
						N_LAYER, new_tau);
  } else {
    rayTracingKernReg_masked <<< gridDim, blockDim >>>(delta_x, N_X,
						       delta_y, N_Y,
						       src->wavefront.phase,
						       src->wavefront.M->m,
						       src->dev_ptr, src->N_SRC,
						       powf(r0,-5.0/6.0),
						       d__layers,
						       d__phase_screen_LAYER,
						       N_LAYER, new_tau);
  }
}