#line 75 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
#include "segmentPistonSensor.h"

#line 560 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  __global__ void lenslet_trim(float *piecewise_phase, 
                               float *piecewise_amplitude, 
                               char *m,
                               const int n_px_lenslet, 
                               float *pupil_phase, 
                               float *pupil_amplitude, 
                               const int n_in, 
                               const float ri, const float ro,
                               const int N_LENSLET,
                               const float lenslet_size,
                               const float lenslet_height,
                               const float R,
                               source *src,
                               const float m2px)
{
  int i, j, k_out, k_in, iSource, k_LA,
    i_lenslet, j_lenslet, k_lenslet, 
    i_in, j_in;
  float theta, xc, yc, x0, y0, O, 
    s, c, h, x, y, xr, yr, pupil;

  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  iSource = blockIdx.z;
  if ( (i<n_px_lenslet) && (j<n_px_lenslet) ) {

    
#line 606 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
O = R - lenslet_size;
x0 = lenslet_height*src[iSource].zenith*cos(src[iSource].azimuth) + O;
y0 = lenslet_height*src[iSource].zenith*sin(src[iSource].azimuth) + O;

#line 588 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      k_lenslet = 12*iSource;

    for (k_LA=0;k_LA<6;k_LA++) {

      
#line 613 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
theta = k_LA*PI/3.0;
j_lenslet = k_lenslet/N_LENSLET;
i_lenslet = k_lenslet - N_LENSLET*j_lenslet;
sincosf(theta,&s,&c);
xc = ri*c + y0;
yc = ri*s + x0;
#line 593 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      
#line 638 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
i_in = (int) floor(yc*m2px);
j_in = (int) floor(xc*m2px);

k_out = lenslet2array( i,  j,  n_px_lenslet,
  		     i_lenslet,  j_lenslet,  N_LENSLET,  0);
k_in = j + i_in;
k_in += (i + j_in + iSource*n_in)*n_in;

h = lenslet_size*0.5;
x = 2*lenslet_size*(i - (n_px_lenslet-1)*0.5)/(n_px_lenslet-1);
y = 2*lenslet_size*(j - (n_px_lenslet-1)*0.5)/(n_px_lenslet-1);
xr = c*x + y*s;
yr = -s*x + y*c;
#line 659 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
pupil = ( ( (xr>=-h) && (xr<=h) ) && ( (yr>=-h) && (yr<=h) ) ) ? 1.0 : 0.0;
#line 653 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
piecewise_phase[k_out]     = pupil*pupil_phase[k_in];
piecewise_amplitude[k_out] = pupil*pupil_amplitude[k_in];
m[k_out]                   = (piecewise_amplitude[k_out]==1) ? 1 : 0;

#line 595 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      
#line 622 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
theta = (k_LA-0.5)*PI/3.0;
j_lenslet = (k_lenslet+6)/N_LENSLET;
i_lenslet = (k_lenslet+6) - N_LENSLET*j_lenslet;
sincosf(theta,&s,&c);
xc = ro*c + y0;
yc = ro*s + x0;

#line 596 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      
#line 638 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
i_in = (int) floor(yc*m2px);
j_in = (int) floor(xc*m2px);

k_out = lenslet2array( i,  j,  n_px_lenslet,
  		     i_lenslet,  j_lenslet,  N_LENSLET,  0);
k_in = j + i_in;
k_in += (i + j_in + iSource*n_in)*n_in;

h = lenslet_size*0.5;
x = 2*lenslet_size*(i - (n_px_lenslet-1)*0.5)/(n_px_lenslet-1);
y = 2*lenslet_size*(j - (n_px_lenslet-1)*0.5)/(n_px_lenslet-1);
xr = c*x + y*s;
yr = -s*x + y*c;
#line 659 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
pupil = ( ( (xr>=-h) && (xr<=h) ) && ( (yr>=-h) && (yr<=h) ) ) ? 1.0 : 0.0;
#line 653 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
piecewise_phase[k_out]     = pupil*pupil_phase[k_in];
piecewise_amplitude[k_out] = pupil*pupil_amplitude[k_in];
m[k_out]                   = (piecewise_amplitude[k_out]==1) ? 1 : 0;

#line 598 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
        ++k_lenslet;
    }
  }
}
#line 716 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
__global__ void lenslet_trim_with_middle_mask(float *piecewise_phase, 
                             float *piecewise_amplitude, 
                             char *m,
                             const int n_px_lenslet, 
                             float *pupil_phase, 
                             float *pupil_amplitude, 
                             const int n_in, 
                             const float ri, const float ro,
                             const int N_LENSLET,
                             const float lenslet_size,
                             const float lenslet_height,
                             const float R,
                             source *src,
                             const float m2px,
                             float w)
{
  int i, j, k_out, k_in, iSource, k_LA,
    i_lenslet, j_lenslet, k_lenslet, 
    i_in, j_in;
  float theta, xc, yc, x0, y0, O, 
    s, c, h, x, y, xr, yr, pupil;

  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  iSource = blockIdx.z;
  if ( (i<n_px_lenslet) && (j<n_px_lenslet) ) {

    
#line 606 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
O = R - lenslet_size;
x0 = lenslet_height*src[iSource].zenith*cos(src[iSource].azimuth) + O;
y0 = lenslet_height*src[iSource].zenith*sin(src[iSource].azimuth) + O;

#line 745 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      k_lenslet = 12*iSource;

    for (k_LA=0;k_LA<6;k_LA++) {

      
#line 613 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
theta = k_LA*PI/3.0;
j_lenslet = k_lenslet/N_LENSLET;
i_lenslet = k_lenslet - N_LENSLET*j_lenslet;
sincosf(theta,&s,&c);
xc = ri*c + y0;
yc = ri*s + x0;
#line 750 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      
#line 638 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
i_in = (int) floor(yc*m2px);
j_in = (int) floor(xc*m2px);

k_out = lenslet2array( i,  j,  n_px_lenslet,
  		     i_lenslet,  j_lenslet,  N_LENSLET,  0);
k_in = j + i_in;
k_in += (i + j_in + iSource*n_in)*n_in;

h = lenslet_size*0.5;
x = 2*lenslet_size*(i - (n_px_lenslet-1)*0.5)/(n_px_lenslet-1);
y = 2*lenslet_size*(j - (n_px_lenslet-1)*0.5)/(n_px_lenslet-1);
xr = c*x + y*s;
yr = -s*x + y*c;
#line 675 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  
#line 663 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  pupil = ( ( (yr>=-h) && (yr<=h) ) &&
            ( ( (xr>=-h) && (xr< -w) ) ||
              ( (xr> +w) && (xr<=+h) ) ) ) ? 1.0 : 0.0;
#line 653 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
piecewise_phase[k_out]     = pupil*pupil_phase[k_in];
piecewise_amplitude[k_out] = pupil*pupil_amplitude[k_in];
m[k_out]                   = (piecewise_amplitude[k_out]==1) ? 1 : 0;

#line 752 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      
#line 622 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
theta = (k_LA-0.5)*PI/3.0;
j_lenslet = (k_lenslet+6)/N_LENSLET;
i_lenslet = (k_lenslet+6) - N_LENSLET*j_lenslet;
sincosf(theta,&s,&c);
xc = ro*c + y0;
yc = ro*s + x0;

#line 753 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      
#line 638 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
i_in = (int) floor(yc*m2px);
j_in = (int) floor(xc*m2px);

k_out = lenslet2array( i,  j,  n_px_lenslet,
  		     i_lenslet,  j_lenslet,  N_LENSLET,  0);
k_in = j + i_in;
k_in += (i + j_in + iSource*n_in)*n_in;

h = lenslet_size*0.5;
x = 2*lenslet_size*(i - (n_px_lenslet-1)*0.5)/(n_px_lenslet-1);
y = 2*lenslet_size*(j - (n_px_lenslet-1)*0.5)/(n_px_lenslet-1);
xr = c*x + y*s;
yr = -s*x + y*c;
#line 668 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  pupil = ( ( (xr>=-h) && (xr<=h) ) &&
            ( ( (yr>=-h) && (yr< -w) ) ||
              ( (yr> +w) && (yr<=+h) ) ) ) ? 1.0 : 0.0;
#line 653 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
piecewise_phase[k_out]     = pupil*pupil_phase[k_in];
piecewise_amplitude[k_out] = pupil*pupil_amplitude[k_in];
m[k_out]                   = (piecewise_amplitude[k_out]==1) ? 1 : 0;

#line 755 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      ++k_lenslet;
    }
  }
}

#line 254 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::setup(gmt_m1 *M1, source *src, 
                                float _dispersion_,
                                float _field_of_view_,
                                float _nyquist_factor_)
{
  int k, N_DFT, DFT_osf;
  BIN_IMAGE=2;
  lenslet_size   = 1.5;
  
#line 332 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_height = 82.5;
ri             = (M1->D_full+0.357)*0.5 ;
ro             = M1->L*sqrt(3)*0.5;
m2px           = src->rays.N_L/src->rays.L;
R              = src->rays.L*0.5;
D_px2          = D_px = src->rays.N_L;
D_px2         *= D_px;
dispersion     = _dispersion_;
field_of_view  = _field_of_view_;
nyquist_factor = _nyquist_factor_;
N_GS           = src->N_SRC;
#line 346 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_LENSLET2 =  N_PX_LENSLET = 2*ceil(lenslet_size*m2px);;
N_PX_LENSLET2 *= N_PX_LENSLET;
#line 351 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_LENSLET2 =  N_LENSLET = ceil(sqrt(12*N_GS));
N_LENSLET2 *= N_LENSLET;
#line 356 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX2 =  N_PX = N_PX_LENSLET*N_LENSLET;
N_PX2 *= N_PX;
lenslet.setup(N_PX2, 1);
lenslet_mask.setup(N_PX2);
HANDLE_ERROR( cudaMemset(lenslet_mask.m,  0, sizeof(char)*N_PX2 ) );
lenslet.M = &lenslet_mask;
HANDLE_ERROR( cudaMemset( lenslet.amplitude, 0 , sizeof(float)*lenslet.N_PX) );
#line 366 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
DFT_osf            = (int) rint(2*nyquist_factor);
if (dispersion>0) {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength() - spectral_bandwidth*0.5;
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = (int) 2*dispersion*spectral_bandwidth/pixel_scale;
 } else {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength();
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = 1;
 }
#line 382 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_IMAGE  = (int) ceil(field_of_view/pixel_scale);
N_PX_IMAGE += N_PX_IMAGE%2;
camera.setup(N_PX_LENSLET-1,N_LENSLET,DFT_osf,N_PX_IMAGE,BIN_IMAGE,1);
FFT.setup(camera.N_PX_CAMERA-1,N_LENSLET,2,camera.N_PX_CAMERA*2,1,1);
#line 407 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_src.setup(src[0].photometric_band,0.0,0.0,INFINITY);
lenslet_src.wavefront.N         = lenslet.N;
lenslet_src.wavefront.N_PX      = lenslet.N_PX;
lenslet_src.wavefront.amplitude = lenslet.amplitude;
lenslet_src.wavefront.phase     = lenslet.phase;
lenslet_src.wavefront.M         = lenslet.M;
#line 417 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
fft_src.setup("V",0.0,0.0,INFINITY);
fft_src.wavefront.N         = 1;
fft_src.wavefront.N_PX      = camera.N_PX_CAMERA*camera.N_SIDE_LENSLET;
fft_src.wavefront.N_PX     *= fft_src.wavefront.N_PX;
fft_src.wavefront.amplitude = camera.d__frame;

HANDLE_ERROR( cudaMalloc( (void**)&fft_phase, sizeof(float)*fft_src.wavefront.N_PX ) );
HANDLE_ERROR( cudaMemset( fft_phase, 0 , sizeof(float)*fft_src.wavefront.N_PX ) );
fft_src.wavefront.phase     = fft_phase;

fft_mask.setup(fft_src.wavefront.N_PX);
fft_mask.area               = 1.0/N_LAMBDA;
fft_src.wavefront.M         = &fft_mask;
#line 263 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  
#line 390 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
camera_array = (imaging *) malloc( sizeof(imaging)*N_LAMBDA) ;
fprintf(stdout,"@(CEO)>segmentPistonSensor: camera array setup!");
if (N_LAMBDA>1) {
for (k=0;k<N_LAMBDA;k++) {
  N_DFT = (int) DFT_osf*N_PX_LENSLET*(1 + k*(spectral_bandwidth/lambda0)/(N_LAMBDA-1));
  camera_array[k].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
}
 } else {
  N_DFT = (int) DFT_osf*N_PX_LENSLET;
  camera_array[k].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
 }
#line 264 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  info();
}
void segmentPistonSensor::setup(gmt_m1 *M1, source *src,
                                float _lenslet_size_,
                                float _dispersion_,
                                float _field_of_view_,
                                float _nyquist_factor_)
{
  int k, N_DFT, DFT_osf;
  BIN_IMAGE=2;
  lenslet_size   = _lenslet_size_;
  
#line 332 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_height = 82.5;
ri             = (M1->D_full+0.357)*0.5 ;
ro             = M1->L*sqrt(3)*0.5;
m2px           = src->rays.N_L/src->rays.L;
R              = src->rays.L*0.5;
D_px2          = D_px = src->rays.N_L;
D_px2         *= D_px;
dispersion     = _dispersion_;
field_of_view  = _field_of_view_;
nyquist_factor = _nyquist_factor_;
N_GS           = src->N_SRC;
#line 346 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_LENSLET2 =  N_PX_LENSLET = 2*ceil(lenslet_size*m2px);;
N_PX_LENSLET2 *= N_PX_LENSLET;
#line 351 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_LENSLET2 =  N_LENSLET = ceil(sqrt(12*N_GS));
N_LENSLET2 *= N_LENSLET;
#line 356 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX2 =  N_PX = N_PX_LENSLET*N_LENSLET;
N_PX2 *= N_PX;
lenslet.setup(N_PX2, 1);
lenslet_mask.setup(N_PX2);
HANDLE_ERROR( cudaMemset(lenslet_mask.m,  0, sizeof(char)*N_PX2 ) );
lenslet.M = &lenslet_mask;
HANDLE_ERROR( cudaMemset( lenslet.amplitude, 0 , sizeof(float)*lenslet.N_PX) );
#line 366 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
DFT_osf            = (int) rint(2*nyquist_factor);
if (dispersion>0) {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength() - spectral_bandwidth*0.5;
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = (int) 2*dispersion*spectral_bandwidth/pixel_scale;
 } else {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength();
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = 1;
 }
#line 382 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_IMAGE  = (int) ceil(field_of_view/pixel_scale);
N_PX_IMAGE += N_PX_IMAGE%2;
camera.setup(N_PX_LENSLET-1,N_LENSLET,DFT_osf,N_PX_IMAGE,BIN_IMAGE,1);
FFT.setup(camera.N_PX_CAMERA-1,N_LENSLET,2,camera.N_PX_CAMERA*2,1,1);
#line 407 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_src.setup(src[0].photometric_band,0.0,0.0,INFINITY);
lenslet_src.wavefront.N         = lenslet.N;
lenslet_src.wavefront.N_PX      = lenslet.N_PX;
lenslet_src.wavefront.amplitude = lenslet.amplitude;
lenslet_src.wavefront.phase     = lenslet.phase;
lenslet_src.wavefront.M         = lenslet.M;
#line 417 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
fft_src.setup("V",0.0,0.0,INFINITY);
fft_src.wavefront.N         = 1;
fft_src.wavefront.N_PX      = camera.N_PX_CAMERA*camera.N_SIDE_LENSLET;
fft_src.wavefront.N_PX     *= fft_src.wavefront.N_PX;
fft_src.wavefront.amplitude = camera.d__frame;

HANDLE_ERROR( cudaMalloc( (void**)&fft_phase, sizeof(float)*fft_src.wavefront.N_PX ) );
HANDLE_ERROR( cudaMemset( fft_phase, 0 , sizeof(float)*fft_src.wavefront.N_PX ) );
fft_src.wavefront.phase     = fft_phase;

fft_mask.setup(fft_src.wavefront.N_PX);
fft_mask.area               = 1.0/N_LAMBDA;
fft_src.wavefront.M         = &fft_mask;
#line 276 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  
#line 390 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
camera_array = (imaging *) malloc( sizeof(imaging)*N_LAMBDA) ;
fprintf(stdout,"@(CEO)>segmentPistonSensor: camera array setup!");
if (N_LAMBDA>1) {
for (k=0;k<N_LAMBDA;k++) {
  N_DFT = (int) DFT_osf*N_PX_LENSLET*(1 + k*(spectral_bandwidth/lambda0)/(N_LAMBDA-1));
  camera_array[k].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
}
 } else {
  N_DFT = (int) DFT_osf*N_PX_LENSLET;
  camera_array[k].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
 }
#line 277 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      info();
}
#line 281 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::setup(gmt_m1 *M1, source *src, 
                                float _dispersion_,
                                float _field_of_view_,
                                float _nyquist_factor_,
                                int _BIN_IMAGE_)
{
  int k, N_DFT, DFT_osf;
  BIN_IMAGE=_BIN_IMAGE_;
  lenslet_size   = 1.5;
  
#line 332 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_height = 82.5;
ri             = (M1->D_full+0.357)*0.5 ;
ro             = M1->L*sqrt(3)*0.5;
m2px           = src->rays.N_L/src->rays.L;
R              = src->rays.L*0.5;
D_px2          = D_px = src->rays.N_L;
D_px2         *= D_px;
dispersion     = _dispersion_;
field_of_view  = _field_of_view_;
nyquist_factor = _nyquist_factor_;
N_GS           = src->N_SRC;
#line 346 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_LENSLET2 =  N_PX_LENSLET = 2*ceil(lenslet_size*m2px);;
N_PX_LENSLET2 *= N_PX_LENSLET;
#line 351 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_LENSLET2 =  N_LENSLET = ceil(sqrt(12*N_GS));
N_LENSLET2 *= N_LENSLET;
#line 356 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX2 =  N_PX = N_PX_LENSLET*N_LENSLET;
N_PX2 *= N_PX;
lenslet.setup(N_PX2, 1);
lenslet_mask.setup(N_PX2);
HANDLE_ERROR( cudaMemset(lenslet_mask.m,  0, sizeof(char)*N_PX2 ) );
lenslet.M = &lenslet_mask;
HANDLE_ERROR( cudaMemset( lenslet.amplitude, 0 , sizeof(float)*lenslet.N_PX) );
#line 366 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
DFT_osf            = (int) rint(2*nyquist_factor);
if (dispersion>0) {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength() - spectral_bandwidth*0.5;
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = (int) 2*dispersion*spectral_bandwidth/pixel_scale;
 } else {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength();
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = 1;
 }
#line 382 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_IMAGE  = (int) ceil(field_of_view/pixel_scale);
N_PX_IMAGE += N_PX_IMAGE%2;
camera.setup(N_PX_LENSLET-1,N_LENSLET,DFT_osf,N_PX_IMAGE,BIN_IMAGE,1);
FFT.setup(camera.N_PX_CAMERA-1,N_LENSLET,2,camera.N_PX_CAMERA*2,1,1);
#line 407 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_src.setup(src[0].photometric_band,0.0,0.0,INFINITY);
lenslet_src.wavefront.N         = lenslet.N;
lenslet_src.wavefront.N_PX      = lenslet.N_PX;
lenslet_src.wavefront.amplitude = lenslet.amplitude;
lenslet_src.wavefront.phase     = lenslet.phase;
lenslet_src.wavefront.M         = lenslet.M;
#line 417 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
fft_src.setup("V",0.0,0.0,INFINITY);
fft_src.wavefront.N         = 1;
fft_src.wavefront.N_PX      = camera.N_PX_CAMERA*camera.N_SIDE_LENSLET;
fft_src.wavefront.N_PX     *= fft_src.wavefront.N_PX;
fft_src.wavefront.amplitude = camera.d__frame;

HANDLE_ERROR( cudaMalloc( (void**)&fft_phase, sizeof(float)*fft_src.wavefront.N_PX ) );
HANDLE_ERROR( cudaMemset( fft_phase, 0 , sizeof(float)*fft_src.wavefront.N_PX ) );
fft_src.wavefront.phase     = fft_phase;

fft_mask.setup(fft_src.wavefront.N_PX);
fft_mask.area               = 1.0/N_LAMBDA;
fft_src.wavefront.M         = &fft_mask;
#line 291 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  
#line 390 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
camera_array = (imaging *) malloc( sizeof(imaging)*N_LAMBDA) ;
fprintf(stdout,"@(CEO)>segmentPistonSensor: camera array setup!");
if (N_LAMBDA>1) {
for (k=0;k<N_LAMBDA;k++) {
  N_DFT = (int) DFT_osf*N_PX_LENSLET*(1 + k*(spectral_bandwidth/lambda0)/(N_LAMBDA-1));
  camera_array[k].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
}
 } else {
  N_DFT = (int) DFT_osf*N_PX_LENSLET;
  camera_array[k].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
 }
#line 292 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  info();
}
void segmentPistonSensor::setup(gmt_m1 *M1, source *src, 
                                float _lenslet_size_,
                                float _dispersion_,
                                float _field_of_view_,
                                float _nyquist_factor_,
                                int _BIN_IMAGE_)
{
  int k, N_DFT, DFT_osf;
  BIN_IMAGE=_BIN_IMAGE_;
  lenslet_size   = _lenslet_size_;
  
#line 332 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_height = 82.5;
ri             = (M1->D_full+0.357)*0.5 ;
ro             = M1->L*sqrt(3)*0.5;
m2px           = src->rays.N_L/src->rays.L;
R              = src->rays.L*0.5;
D_px2          = D_px = src->rays.N_L;
D_px2         *= D_px;
dispersion     = _dispersion_;
field_of_view  = _field_of_view_;
nyquist_factor = _nyquist_factor_;
N_GS           = src->N_SRC;
#line 346 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_LENSLET2 =  N_PX_LENSLET = 2*ceil(lenslet_size*m2px);;
N_PX_LENSLET2 *= N_PX_LENSLET;
#line 351 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_LENSLET2 =  N_LENSLET = ceil(sqrt(12*N_GS));
N_LENSLET2 *= N_LENSLET;
#line 356 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX2 =  N_PX = N_PX_LENSLET*N_LENSLET;
N_PX2 *= N_PX;
lenslet.setup(N_PX2, 1);
lenslet_mask.setup(N_PX2);
HANDLE_ERROR( cudaMemset(lenslet_mask.m,  0, sizeof(char)*N_PX2 ) );
lenslet.M = &lenslet_mask;
HANDLE_ERROR( cudaMemset( lenslet.amplitude, 0 , sizeof(float)*lenslet.N_PX) );
#line 366 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
DFT_osf            = (int) rint(2*nyquist_factor);
if (dispersion>0) {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength() - spectral_bandwidth*0.5;
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = (int) 2*dispersion*spectral_bandwidth/pixel_scale;
 } else {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength();
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = 1;
 }
#line 382 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_IMAGE  = (int) ceil(field_of_view/pixel_scale);
N_PX_IMAGE += N_PX_IMAGE%2;
camera.setup(N_PX_LENSLET-1,N_LENSLET,DFT_osf,N_PX_IMAGE,BIN_IMAGE,1);
FFT.setup(camera.N_PX_CAMERA-1,N_LENSLET,2,camera.N_PX_CAMERA*2,1,1);
#line 407 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_src.setup(src[0].photometric_band,0.0,0.0,INFINITY);
lenslet_src.wavefront.N         = lenslet.N;
lenslet_src.wavefront.N_PX      = lenslet.N_PX;
lenslet_src.wavefront.amplitude = lenslet.amplitude;
lenslet_src.wavefront.phase     = lenslet.phase;
lenslet_src.wavefront.M         = lenslet.M;
#line 417 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
fft_src.setup("V",0.0,0.0,INFINITY);
fft_src.wavefront.N         = 1;
fft_src.wavefront.N_PX      = camera.N_PX_CAMERA*camera.N_SIDE_LENSLET;
fft_src.wavefront.N_PX     *= fft_src.wavefront.N_PX;
fft_src.wavefront.amplitude = camera.d__frame;

HANDLE_ERROR( cudaMalloc( (void**)&fft_phase, sizeof(float)*fft_src.wavefront.N_PX ) );
HANDLE_ERROR( cudaMemset( fft_phase, 0 , sizeof(float)*fft_src.wavefront.N_PX ) );
fft_src.wavefront.phase     = fft_phase;

fft_mask.setup(fft_src.wavefront.N_PX);
fft_mask.area               = 1.0/N_LAMBDA;
fft_src.wavefront.M         = &fft_mask;
#line 305 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  
#line 390 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
camera_array = (imaging *) malloc( sizeof(imaging)*N_LAMBDA) ;
fprintf(stdout,"@(CEO)>segmentPistonSensor: camera array setup!");
if (N_LAMBDA>1) {
for (k=0;k<N_LAMBDA;k++) {
  N_DFT = (int) DFT_osf*N_PX_LENSLET*(1 + k*(spectral_bandwidth/lambda0)/(N_LAMBDA-1));
  camera_array[k].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
}
 } else {
  N_DFT = (int) DFT_osf*N_PX_LENSLET;
  camera_array[k].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
 }
#line 306 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
      info();
}
#line 316 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::setup_alt(gmt_m1 *M1, source *src, 
                                    float _dispersion_,
                                    float _field_of_view_,
                                    float _nyquist_factor_,
                                    int _BIN_IMAGE_)
{
  int DFT_osf;
  BIN_IMAGE=_BIN_IMAGE_;
  lenslet_size   = 1.5;
  
#line 332 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_height = 82.5;
ri             = (M1->D_full+0.357)*0.5 ;
ro             = M1->L*sqrt(3)*0.5;
m2px           = src->rays.N_L/src->rays.L;
R              = src->rays.L*0.5;
D_px2          = D_px = src->rays.N_L;
D_px2         *= D_px;
dispersion     = _dispersion_;
field_of_view  = _field_of_view_;
nyquist_factor = _nyquist_factor_;
N_GS           = src->N_SRC;
#line 346 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_LENSLET2 =  N_PX_LENSLET = 2*ceil(lenslet_size*m2px);;
N_PX_LENSLET2 *= N_PX_LENSLET;
#line 351 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_LENSLET2 =  N_LENSLET = ceil(sqrt(12*N_GS));
N_LENSLET2 *= N_LENSLET;
#line 356 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX2 =  N_PX = N_PX_LENSLET*N_LENSLET;
N_PX2 *= N_PX;
lenslet.setup(N_PX2, 1);
lenslet_mask.setup(N_PX2);
HANDLE_ERROR( cudaMemset(lenslet_mask.m,  0, sizeof(char)*N_PX2 ) );
lenslet.M = &lenslet_mask;
HANDLE_ERROR( cudaMemset( lenslet.amplitude, 0 , sizeof(float)*lenslet.N_PX) );
#line 366 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
DFT_osf            = (int) rint(2*nyquist_factor);
if (dispersion>0) {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength() - spectral_bandwidth*0.5;
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = (int) 2*dispersion*spectral_bandwidth/pixel_scale;
 } else {
  spectral_bandwidth = src[0].spectral_bandwidth();
  lambda0            = src[0].wavelength();
  pixel_scale = lambda0/(DFT_osf*2*lenslet_size);
  N_LAMBDA    = 1;
 }
#line 382 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
N_PX_IMAGE  = (int) ceil(field_of_view/pixel_scale);
N_PX_IMAGE += N_PX_IMAGE%2;
camera.setup(N_PX_LENSLET-1,N_LENSLET,DFT_osf,N_PX_IMAGE,BIN_IMAGE,1);
FFT.setup(camera.N_PX_CAMERA-1,N_LENSLET,2,camera.N_PX_CAMERA*2,1,1);
#line 407 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
lenslet_src.setup(src[0].photometric_band,0.0,0.0,INFINITY);
lenslet_src.wavefront.N         = lenslet.N;
lenslet_src.wavefront.N_PX      = lenslet.N_PX;
lenslet_src.wavefront.amplitude = lenslet.amplitude;
lenslet_src.wavefront.phase     = lenslet.phase;
lenslet_src.wavefront.M         = lenslet.M;
#line 417 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
fft_src.setup("V",0.0,0.0,INFINITY);
fft_src.wavefront.N         = 1;
fft_src.wavefront.N_PX      = camera.N_PX_CAMERA*camera.N_SIDE_LENSLET;
fft_src.wavefront.N_PX     *= fft_src.wavefront.N_PX;
fft_src.wavefront.amplitude = camera.d__frame;

HANDLE_ERROR( cudaMalloc( (void**)&fft_phase, sizeof(float)*fft_src.wavefront.N_PX ) );
HANDLE_ERROR( cudaMemset( fft_phase, 0 , sizeof(float)*fft_src.wavefront.N_PX ) );
fft_src.wavefront.phase     = fft_phase;

fft_mask.setup(fft_src.wavefront.N_PX);
fft_mask.area               = 1.0/N_LAMBDA;
fft_src.wavefront.M         = &fft_mask;
#line 326 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  
#line 402 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
camera_array = (imaging *) malloc( sizeof(imaging)) ;
#line 327 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  info();
}
#line 434 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::cleanup(void)
{
  fprintf(stdout,"@(CEO)>segmentPistonSensor: freeing memory!\n");
  lenslet.cleanup();
  lenslet_mask.cleanup();
  camera.cleanup();
  for (int k=0;k<N_LAMBDA;k++) {
    camera_array[k].cleanupSegmentPistonSensor();
  }
  free( camera_array );
  FFT.cleanup();
  fft_mask.cleanup();
  HANDLE_ERROR( cudaFree( fft_phase ) );
}
#line 451 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::cleanup_alt(void)
{
  fprintf(stdout,"@(CEO)>segmentPistonSensor: freeing memory!\n");
  lenslet.cleanup();
  lenslet_mask.cleanup();
  camera.cleanup();
  free( camera_array );
  FFT.cleanup();
  fft_mask.cleanup();
  HANDLE_ERROR( cudaFree( fft_phase ) );
}
#line 469 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::propagate(source *gs)
{
  int k;
  float d, wavenumber;
  
#line 480 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
dim3 blockDim(16,16);
dim3 gridDim(N_PX_LENSLET/16+1,N_PX_LENSLET/16+1,N_GS);
lenslet_trim <<< gridDim, blockDim >>> (lenslet.phase, 
					lenslet.amplitude, 
					lenslet_mask.m,
					N_PX_LENSLET,
					gs->wavefront.phase, 
					gs->wavefront.amplitude, 
					D_px, ri, ro, N_LENSLET,
					lenslet_size, 
					lenslet_height, R,
					gs->dev_ptr, m2px);
lenslet_mask.set_filter_quiet();
lenslet_mask.area = N_GS*lenslet_mask.nnz*
  gs->wavefront.M->area/gs->wavefront.M->nnz/N_LAMBDA;
#line 474 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  
#line 498 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  //lenslet_src.magnitude = gs[0].magnitude;
lenslet_src.copy_magnitude(gs);
lenslet_src.fwhm = gs->fwhm;
if (N_LAMBDA>1) {
for (k=0;k<N_LAMBDA;k++) {
  camera_array[k].d__frame = camera.d__frame;
  d = spectral_bandwidth*( k - 0.5*(N_LAMBDA-1))/(N_LAMBDA-1);
  wavenumber = 2*PI/(gs[0].wavelength() + d);
  d *= (dispersion/pixel_scale);
  camera_array[k].propagateNoOverlapSPS(&lenslet_src, d, wavenumber);
}
 } else {
  camera_array[k].d__frame = camera.d__frame;
  d = 0.0;
  wavenumber = 2*PI/(gs[0].wavelength() + d);
  camera_array[k].propagateNoOverlapSPS(&lenslet_src, d, wavenumber);
 }
  camera.N_FRAME++;
#line 475 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
}
#line 519 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::propagate_alt(source *gs)
{
  int k, N_DFT, DFT_osf;
  float d, wavenumber;
  
#line 480 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
dim3 blockDim(16,16);
dim3 gridDim(N_PX_LENSLET/16+1,N_PX_LENSLET/16+1,N_GS);
lenslet_trim <<< gridDim, blockDim >>> (lenslet.phase, 
					lenslet.amplitude, 
					lenslet_mask.m,
					N_PX_LENSLET,
					gs->wavefront.phase, 
					gs->wavefront.amplitude, 
					D_px, ri, ro, N_LENSLET,
					lenslet_size, 
					lenslet_height, R,
					gs->dev_ptr, m2px);
lenslet_mask.set_filter_quiet();
lenslet_mask.area = N_GS*lenslet_mask.nnz*
  gs->wavefront.M->area/gs->wavefront.M->nnz/N_LAMBDA;
#line 524 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  DFT_osf            = (int) rint(2*nyquist_factor);
  
#line 528 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  //lenslet_src.magnitude = gs[0].magnitude;
lenslet_src.copy_magnitude(gs);
lenslet_src.fwhm = gs->fwhm;
if (N_LAMBDA>1) {
for (k=0;k<N_LAMBDA;k++) {
  N_DFT = (int) DFT_osf*N_PX_LENSLET*(1 + k*(spectral_bandwidth/lambda0)/(N_LAMBDA-1));
  camera_array[0].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
  camera_array[0].d__frame = camera.d__frame;
  d = spectral_bandwidth*( k - 0.5*(N_LAMBDA-1))/(N_LAMBDA-1);
  wavenumber = 2*PI/(gs[0].wavelength() + d);
  d *= (dispersion/pixel_scale);
  camera_array[0].propagateNoOverlapSPS(&lenslet_src, d, wavenumber);
  camera_array[0].cleanupSegmentPistonSensor();

}
 } else {
  N_DFT = (int) DFT_osf*N_PX_LENSLET;
  camera_array[0].setupSegmentPistonSensor(N_PX_LENSLET-1,N_LENSLET,N_DFT,N_PX_IMAGE,BIN_IMAGE,1);
  camera_array[0].d__frame = camera.d__frame;
  d = 0.0;
  wavenumber = 2*PI/(gs[0].wavelength() + d);
  camera_array[0].propagateNoOverlapSPS(&lenslet_src, d, wavenumber);
  camera_array[0].cleanupSegmentPistonSensor();
 }
  camera.N_FRAME++;
#line 526 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
}
#line 692 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::propagate(source *gs, float middle_mask_width)
{
  int k;
  float d, wavenumber;
  dim3 blockDim(16,16);
  dim3 gridDim(N_PX_LENSLET/16+1,N_PX_LENSLET/16+1,N_GS);
  lenslet_trim_with_middle_mask <<< gridDim, blockDim >>> (lenslet.phase, 
                                                      lenslet.amplitude, 
                                                      lenslet_mask.m,
                                                      N_PX_LENSLET,
                                                      gs->wavefront.phase, 
                                                      gs->wavefront.amplitude, 
                                                      D_px, ri, ro, N_LENSLET,
                                                      lenslet_size, 
                                                      lenslet_height, R,
                                                      gs->dev_ptr, m2px,
                                                      0.5*middle_mask_width);
  lenslet_mask.set_filter_quiet();
  lenslet_mask.area = N_GS*lenslet_mask.nnz*
  gs->wavefront.M->area/gs->wavefront.M->nnz/N_LAMBDA;
  
#line 498 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
  //lenslet_src.magnitude = gs[0].magnitude;
lenslet_src.copy_magnitude(gs);
lenslet_src.fwhm = gs->fwhm;
if (N_LAMBDA>1) {
for (k=0;k<N_LAMBDA;k++) {
  camera_array[k].d__frame = camera.d__frame;
  d = spectral_bandwidth*( k - 0.5*(N_LAMBDA-1))/(N_LAMBDA-1);
  wavenumber = 2*PI/(gs[0].wavelength() + d);
  d *= (dispersion/pixel_scale);
  camera_array[k].propagateNoOverlapSPS(&lenslet_src, d, wavenumber);
}
 } else {
  camera_array[k].d__frame = camera.d__frame;
  d = 0.0;
  wavenumber = 2*PI/(gs[0].wavelength() + d);
  camera_array[k].propagateNoOverlapSPS(&lenslet_src, d, wavenumber);
 }
  camera.N_FRAME++;
#line 713 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
}
#line 766 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::readout(float exposureTime, float readOutNoiseRms, 
				  float backgroundMagnitude)
{
  float p2, nBackgroundPhoton;
  p2 = pixel_scale*camera.BIN_IMAGE/ARCSEC(1);
  p2 *= p2;
  fprintf(stdout,"p2=%g\n",p2);
  /*
  lenslet_src.magnitude = backgroundMagnitude;
  nBackgroundPhoton = p2*lenslet_src.n_photon();
  nBackgroundPhoton *= lenslet_mask.area*N_LAMBDA/N_GS/12;
  */
  nBackgroundPhoton = lenslet_src.n_background_photon(backgroundMagnitude*p2)/N_LAMBDA/12;
  fprintf(stdout,"nBackgroundPhoton=%g\n",exposureTime*nBackgroundPhoton);
  camera.readout(exposureTime, readOutNoiseRms, nBackgroundPhoton, 1.0);
}
#line 789 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::fft(void)
{
  fft_src.wavefront.masked();
  FFT.propagateNoOverlapBare( &fft_src );
}
#line 801 "/home/ubuntu/projects/crseo/sys/CEO/segmentPistonSensor/segmentPistonSensor.nw"
void segmentPistonSensor::info(void)
{
  fprintf(stdout,"\n\x1B[1;42m@(CEO)>segmentPistonSensor:\x1B[;42m\n");
  fprintf(stdout," . lenslet size and conjugation height: %3.1fm and %4.1fm\n",
	  lenslet_size, lenslet_height);
  fprintf(stdout," . center wavelength, spectral bandwidth and resolution: %6.3fmicron, %6.3fmicron, %d\n",
	  (lambda0+spectral_bandwidth*0.5)*1e6,
	  spectral_bandwidth*1e6, N_LAMBDA);
  fprintf(stdout," . pixel scale       : %6.3farcsec\n",pixel_scale*camera.BIN_IMAGE/ARCSEC(1));
  fprintf(stdout," . field-of-view     : %6.3farcsec\n",N_PX_IMAGE*pixel_scale/ARCSEC(1));
  fprintf(stdout,"----------------------------------------------------\x1B[0m\n"); 
}