#line 50 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
#include "gmtMirrors.h"
#line 2893 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__device__ rtd bending_modes_surface(vector *v, const int N_mode, rtd *b,
                                            double *BM,
                                            const int BM_N_SAMPLE, const double BM_radius)
{
  int ndx;
  rtd scale, S, s, t, fs, ft, onemt, onems;
  rtd *z;
  scale = 0.5/BM_radius;
  s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1);
  t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1);
  
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2904 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  z = BM;
  S = ( ( z[ndx]*onemt + z[ndx+1]*t )*onems +
	( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s );
  return -S;
}
#line 2912 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__device__ rtd partial_x_bending_modes_surface(vector *v, const int N_mode, rtd *b,
						      double *BM,
						      const int BM_N_SAMPLE, const double BM_radius)
{
  int ndx, ndx_a, ndx_b;
  rtd scale, S, S_, s, t, fs, ft, onemt, onems,
    s_a, s_b, onems_a, onems_b, h;
  rtd *z;
  scale = 0.5/BM_radius;
  s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1) + 1;
  t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1);
  
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2924 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  ndx_a = ndx;
  s_a = s;
  onems_a = onems;
  s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1) - 1;
  t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1);
  
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2930 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  ndx_b = ndx;
  s_b = s;
  onems_b = onems;
  h = BM_radius*2/(BM_N_SAMPLE-1);
  z = BM;
  S_ = ( ( z[ndx_a]*onemt + z[ndx_a+1]*t )*onems_a +
	 ( z[ndx_a+BM_N_SAMPLE]*onemt + z[ndx_a+BM_N_SAMPLE+1]*t )*s_a );
  S_ -= ( ( z[ndx_b]*onemt + z[ndx_b+1]*t )*onems_b +
	  ( z[ndx_b+BM_N_SAMPLE]*onemt + z[ndx_b+BM_N_SAMPLE+1]*t )*s_b );
  S_ /= 2*h;
  S = S_;
  return -S;
}
#line 2946 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__device__ rtd partial_y_bending_modes_surface(vector *v, const int N_mode, rtd *b,
						      double *BM,
						      const int BM_N_SAMPLE, const double BM_radius)
{
  int ndx, ndx_a, ndx_b;
  rtd scale, S, S_, s, t, fs, ft, onemt, onems,
    t_a, t_b, onemt_a, onemt_b, h;
  rtd *z;
  scale = 0.5/BM_radius;
  s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1);
  t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1) + 1;
  
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2958 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  ndx_a = ndx;
  t_a = t;
  onemt_a = onemt;
  s = (scale*v->x + 0.5)*(BM_N_SAMPLE-1);
  t = (scale*v->y + 0.5)*(BM_N_SAMPLE-1) - 1;
  
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2964 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  ndx_b = ndx;
  t_b = t;
  onemt_b = onemt;
  h = BM_radius*2/(BM_N_SAMPLE-1);
  z = BM;
  S_ = ( ( z[ndx_a]*onemt_a + z[ndx_a+1]*t_a )*onems +
	 ( z[ndx_a+BM_N_SAMPLE]*onemt_a + z[ndx_a+BM_N_SAMPLE+1]*t_a )*s );
  S_ -= ( ( z[ndx_b]*onemt_b + z[ndx_b+1]*t_b )*onems +
	  ( z[ndx_b+BM_N_SAMPLE]*onemt_b + z[ndx_b+BM_N_SAMPLE+1]*t_b )*s );
  S_ /= 2*h;
  S = S_;
  return -S;
}
#line 2980 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__device__ rtd partial_z_bending_modes_surface(void)
{
  return 0.0;
}

#line 56 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 922 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void m1_trace_chief_kernel(char *mask, ray *d__ray, int N_RAY,
				      rtd *d__conic_R, vector *d__conic_origin,
				      const rtd Fk, const rtd Fc,
				      vector *d__conic_self_origin,
				      rtd *d__rigid_body_R, vector *d__rigid_body_origin,
				      rtd *d__motion_R, vector *d__motion_origin)
{
  int j, ij, iCoordSys, iSource;
  rtd x1, y1, s0, k, l, m;
  rtd s1, S, K, L ,M, dSds;
  rtd G2, a;
  vector d__origin, xyz, klm, xyz_GS, klm_GS;

  iSource = blockIdx.z;
  ij = iSource*N_RAY;
  iCoordSys = 0;

  d__origin = d__conic_self_origin[iCoordSys];

  
#line 1002 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID BODY >>>
forward_transform(&xyz_GS, &(d__ray[ij].coordinates),
		  d__rigid_body_R+iCoordSys*9, d__rigid_body_origin+iCoordSys);
forward_transform_centered(&klm_GS, &(d__ray[ij].directions),
			   d__rigid_body_R+iCoordSys*9);

#line 943 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 1089 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC >>>
backward_transform_centered(&xyz, &xyz_GS,
			    d__conic_R+iCoordSys*9);
backward_transform_centered(&klm, &klm_GS,
			    d__conic_R+iCoordSys*9);
#line 944 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 1379 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k  = klm.x;
l  = klm.y;
m  = klm.z;
if (m==0) { return; }
s0 = -xyz.z/m;
x1 = xyz.x + k*s0;
y1 = xyz.y + l*s0;
d__ray[ij].optical_path_length = s0;
s0 = s1 = 0;
for (j=0; j<SNELL_N_ITERATION; j++)
{
  xyz.x = x1 + k*s0;
  xyz.y = y1 + l*s0;
  xyz.z = m*s0;
  
#line 1412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
S = conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
#line 1394 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 1415 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
M = partial_z_conic_surface();
#line 1395 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  dSds = K*k + L*l + M*m;
  if (dSds==0) { break; }
  s1 = s0 - S/dSds;
  if (abs(s1-s0)<TOL)
    {
      xyz.x = x1 + k*s1;
      xyz.y = y1 + l*s1;
      xyz.z = m*s1;
      
#line 1415 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
M = partial_z_conic_surface();
#line 1404 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      d__ray[ij].optical_path_length += s1;
      d__ray[ij].n_iteration = j;
      break;
    }
  s0 = s1;
}
#line 945 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 1096 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
G2 = K*K + L*L + M*M;
a  = k*K + l*L + m*M;

a *= -2.0/G2;
klm.x += a*K;
klm.y += a*L;
klm.z += a*M;
#line 946 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 1114 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 // CONIC <<<
 forward_transform_centered(&xyz_GS, &xyz,
			    d__conic_R+iCoordSys*9);
 forward_transform_centered(&klm_GS, &klm,
			    d__conic_R+iCoordSys*9);
#line 947 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 1128 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 // RIGID_BODY <<<
 backward_transform(&(d__ray[ij].coordinates), &xyz_GS,
		    d__rigid_body_R+iCoordSys*9,
		    d__rigid_body_origin+iCoordSys);
 backward_transform_centered(&(d__ray[ij].directions), &klm_GS,
			     d__rigid_body_R+iCoordSys*9);


#line 949 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 803 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void m1_trace_kernel(char *mask, ray *d__ray, int N_RAY,
                                rtd inner2, rtd outer2,
                                rtd *d__aperture_R, vector *d__aperture_origin,
                                rtd *d__conic_R, vector *d__conic_origin,
                                rtd *d__Fk, rtd *d__Fc,
                                vector *d__conic_self_origin,
                                rtd *d__rigid_body_R, vector *d__rigid_body_origin,
                                rtd *d__motion_R, vector *d__motion_origin,
                                int max_n, int n_mode,
                                rtd *d__a, rtd *d__cx, rtd *d__cy, rtd R,
                                ray *d__chief_ray, int *d__piston_mask,
                                char *d__valid_segments, float *reflectivity)
{
  int i, j, ij, iCoordSys, iSource, b_n_mode;
  rtd rho2;
  rtd x, y, z, x1, y1, s0, k, l, m;
  rtd s1, S, K, L ,M, dSds;
  rtd G2, a;
  rtd Fc, Fk;
  vector d__origin, xyz, klm, xyz_GS, klm_GS, Snormal;

  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = threadIdx.y;
  //  iCoordSys = blockIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  iSource = blockIdx.z;
  j = ij;
  ij += iSource*N_RAY;
  b_n_mode = max_n*(max_n+1)*0.5;
  if  ( ( j<N_RAY ) && (d__ray[ij].v) )// & (iCoordSys<1) )
  {

    for (iCoordSys=0; iCoordSys<7 ; iCoordSys++)
      {
	Fc = d__Fc[iCoordSys];
	Fk = d__Fk[iCoordSys];
        if (d__valid_segments[iCoordSys]>0) {
	d__origin = d__conic_self_origin[iCoordSys];

        
#line 1002 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID BODY >>>
forward_transform(&xyz_GS, &(d__ray[ij].coordinates),
		  d__rigid_body_R+iCoordSys*9, d__rigid_body_origin+iCoordSys);
forward_transform_centered(&klm_GS, &(d__ray[ij].directions),
			   d__rigid_body_R+iCoordSys*9);
#line 1009 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION >>>
forward_transform(&xyz_GS, &xyz_GS,
		  d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
forward_transform_centered(&klm_GS, &klm_GS,
			    d__motion_R+iCoordSys*9);

#line 844 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	
#line 1016 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// APERTURE >>>
forward_transform(&xyz, &xyz_GS,
		  d__aperture_R+iCoordSys*9,
		  d__aperture_origin+iCoordSys);
forward_transform_centered(&klm, &klm_GS,
			   d__aperture_R+iCoordSys*9);

	x = xyz.x;
	y = xyz.y;
	z = xyz.z;
	k = klm.x;
	l = klm.y;
	m = klm.z;

        
#line 1154 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (m==0) { return; }
s0 = -z/m;
x1 = x + k*s0;
y1 = y + l*s0;

#line 1032 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
        rho2 = x1*x1 + y1*y1;
	if ( (iCoordSys==6) && (rho2<inner2) )
	  {
	    return;
	  }

#line 846 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	if (rho2<=outer2)
	  {
	    mask[ij] = 1;
	    d__piston_mask[ij] = iCoordSys + 1;

            
#line 1050 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
klm.x = klm_GS.x;
klm.y =	klm_GS.y;
klm.z =	klm_GS.z;
xyz.x = xyz_GS.x;
xyz.y =	xyz_GS.y;
xyz.z =	xyz_GS.z;
#line 1163 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k  = klm.x;
l  = klm.y;
m  = klm.z;
if (m==0) { return; }
s0 = -xyz.z/m;
x1 = xyz.x + k*s0;
y1 = xyz.y + l*s0;
d__ray[ij].optical_path_length = s0;
s0 = s1 = 0;
for (j=0; j<SNELL_N_ITERATION; j++)
{
  xyz.x = x1 + k*s0;
  xyz.y = y1 + l*s0;
  xyz.z = m*s0;
  
#line 1232 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC >>>
backward_transform_centered(&xyz, &xyz,
			    d__conic_R+iCoordSys*9);
S = conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
#line 1334 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
M = partial_z_conic_surface();
#line 1241 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC <<<
forward_transform_centered(&xyz, &xyz,
                           d__conic_R+iCoordSys*9);
S += zernike_surface(&xyz, &d__origin, max_n,
                     d__a + iCoordSys*n_mode, R); 
#line 1338 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
Snormal.x = K;
Snormal.y = L;
Snormal.z = M;
forward_transform_centered(&Snormal, &Snormal,
			   d__conic_R+iCoordSys*9);
K = Snormal.x;
L = Snormal.y;
M = Snormal.z;
K += partial_x_zernike_surface(&xyz, &d__origin, max_n,
                               d__cx + iCoordSys*b_n_mode, R);
L += partial_y_zernike_surface(&xyz, &d__origin, max_n,
                               d__cy + iCoordSys*b_n_mode, R);
M += partial_z_zernike_surface();
#line 1178 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  dSds = K*k + L*l + M*m;
  if (dSds==0) { break; }
  s1 = s0 - S/dSds;
  if (abs(s1-s0)<TOL)
    {
      xyz.x = x1 + k*s1;
      xyz.y = y1 + l*s1;
      xyz.z = m*s1;
      
#line 1276 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC >>>
backward_transform_centered(&xyz, &xyz,
			    d__conic_R+iCoordSys*9);
K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
M = partial_z_conic_surface();
// CONIC <<<
forward_transform_centered(&xyz, &xyz,
                           d__conic_R+iCoordSys*9);
Snormal.x = K;
Snormal.y = L;
Snormal.z = M;
forward_transform_centered(&Snormal, &Snormal,
			   d__conic_R+iCoordSys*9);
K = Snormal.x;
L = Snormal.y;
M = Snormal.z;
K += partial_x_zernike_surface(&xyz, &d__origin, max_n,
                               d__cx + iCoordSys*b_n_mode, R);
L += partial_y_zernike_surface(&xyz, &d__origin, max_n,
                               d__cy + iCoordSys*b_n_mode, R);
M += partial_z_zernike_surface();
#line 1187 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      d__ray[ij].optical_path_length += s1;
      d__ray[ij].n_iteration = j;
      break;
    }
  s0 = s1;
}
#line 1096 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
G2 = K*K + L*L + M*M;
a  = k*K + l*L + m*M;

a *= -2.0/G2;
klm.x += a*K;
klm.y += a*L;
klm.z += a*M;
#line 1058 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
klm_GS.x = klm.x;
klm_GS.y = klm.y;
klm_GS.z = klm.z;
xyz_GS.x = xyz.x;
xyz_GS.y = xyz.y;
xyz_GS.z = xyz.z;
#line 1121 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 // MOTION <<<
 backward_transform(&xyz_GS, &xyz_GS,
		    d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
 backward_transform_centered(&klm_GS, &klm_GS,
			     d__motion_R+iCoordSys*9);
#line 1128 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 // RIGID_BODY <<<
 backward_transform(&(d__ray[ij].coordinates), &xyz_GS,
		    d__rigid_body_R+iCoordSys*9,
		    d__rigid_body_origin+iCoordSys);
 backward_transform_centered(&(d__ray[ij].directions), &klm_GS,
			     d__rigid_body_R+iCoordSys*9);

#line 852 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	    d__ray[ij].optical_path_difference +=
         d__ray[ij].optical_path_length - d__chief_ray[iSource].optical_path_length;
      d__ray[ij].throughput *= reflectivity[iCoordSys];
	    return;
	  }
        }
      }
  }
}
#line 863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void m1_bm_trace_kernel(char *mask, ray *d__ray, int N_RAY,
                                   rtd inner2, rtd outer2,
                                   rtd *d__aperture_R, vector *d__aperture_origin,
                                   rtd *d__conic_R, vector *d__conic_origin,
                                   rtd *d__Fk, rtd *d__Fc,
                                   vector *d__conic_self_origin,
                                   rtd *d__rigid_body_R, vector *d__rigid_body_origin,
                                   rtd *d__motion_R, vector *d__motion_origin,
                                   const int N_mode, rtd *b, double *BM,
                                   const int BM_N_SAMPLE, const double BM_radius,
                                   ray *d__chief_ray, int *d__piston_mask,
                                   char *d__valid_segments, float *reflectivity)
{
  int i, j, ij, iCoordSys, iSource;
  rtd rho2;
  rtd x, y, z, x1, y1, s0, k, l, m;
  rtd s1, S, K, L ,M, dSds;
  rtd G2, a;
  rtd Fc, Fk;
  vector d__origin, xyz, klm, xyz_GS, klm_GS, Snormal;

  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = threadIdx.y;
  //  iCoordSys = blockIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  iSource = blockIdx.z;
  j = ij;
  ij += iSource*N_RAY;
  if  ( ( j<N_RAY ) && (d__ray[ij].v) )// & (iCoordSys<1) )
  {

    for (iCoordSys=0; iCoordSys<7 ; iCoordSys++)
      {
	Fc = d__Fc[iCoordSys];
	Fk = d__Fk[iCoordSys];
        if (d__valid_segments[iCoordSys]>0) {
	d__origin = d__conic_self_origin[iCoordSys];

        
#line 1002 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID BODY >>>
forward_transform(&xyz_GS, &(d__ray[ij].coordinates),
		  d__rigid_body_R+iCoordSys*9, d__rigid_body_origin+iCoordSys);
forward_transform_centered(&klm_GS, &(d__ray[ij].directions),
			   d__rigid_body_R+iCoordSys*9);
#line 1009 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION >>>
forward_transform(&xyz_GS, &xyz_GS,
		  d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
forward_transform_centered(&klm_GS, &klm_GS,
			    d__motion_R+iCoordSys*9);

#line 903 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	
#line 1016 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// APERTURE >>>
forward_transform(&xyz, &xyz_GS,
		  d__aperture_R+iCoordSys*9,
		  d__aperture_origin+iCoordSys);
forward_transform_centered(&klm, &klm_GS,
			   d__aperture_R+iCoordSys*9);

	x = xyz.x;
	y = xyz.y;
	z = xyz.z;
	k = klm.x;
	l = klm.y;
	m = klm.z;

        
#line 1154 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (m==0) { return; }
s0 = -z/m;
x1 = x + k*s0;
y1 = y + l*s0;

#line 1032 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
        rho2 = x1*x1 + y1*y1;
	if ( (iCoordSys==6) && (rho2<inner2) )
	  {
	    return;
	  }

#line 905 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	if (rho2<=outer2)
	  {
	    mask[ij] = 1;
	    d__piston_mask[ij] = iCoordSys + 1;

            
#line 1069 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
klm.x = klm_GS.x;
klm.y =	klm_GS.y;
klm.z =	klm_GS.z;
xyz.x = xyz_GS.x;
xyz.y =	xyz_GS.y;
xyz.z =	xyz_GS.z;
#line 1196 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k  = klm.x;
l  = klm.y;
m  = klm.z;
if (m==0) { return; }
s0 = -xyz.z/m;
x1 = xyz.x + k*s0;
y1 = xyz.y + l*s0;
d__ray[ij].optical_path_length = s0;
s0 = s1 = 0;
for (j=0; j<SNELL_N_ITERATION; j++)
{
  xyz.x = x1 + k*s0;
  xyz.y = y1 + l*s0;
  xyz.z = m*s0;
  
#line 1232 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC >>>
backward_transform_centered(&xyz, &xyz,
			    d__conic_R+iCoordSys*9);
S = conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
#line 1355 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
M = partial_z_conic_surface();
#line 1249 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC <<<
forward_transform_centered(&xyz, &xyz,
                           d__conic_R+iCoordSys*9);

S += bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, 
			   BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE,
			   BM_N_SAMPLE, BM_radius); 

#line 1359 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
Snormal.x = K;
Snormal.y = L;
Snormal.z = M;
forward_transform_centered(&Snormal, &Snormal,
			   d__conic_R+iCoordSys*9);
K = Snormal.x;
L = Snormal.y;
M = Snormal.z;

K += partial_x_bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, 
				     BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE, 
				     BM_N_SAMPLE, BM_radius);
L += partial_y_bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, 
				     BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE,
				     BM_N_SAMPLE, BM_radius);
M += partial_z_bending_modes_surface();
#line 1211 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  dSds = K*k + L*l + M*m;
  if (dSds==0) { break; }
  s1 = s0 - S/dSds;
  if (abs(s1-s0)<TOL)
    {
      xyz.x = x1 + k*s1;
      xyz.y = y1 + l*s1;
      xyz.z = m*s1;
      
#line 1301 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC >>>
backward_transform_centered(&xyz, &xyz,
			    d__conic_R+iCoordSys*9);
K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
M = partial_z_conic_surface();
// CONIC <<<
forward_transform_centered(&xyz, &xyz,
                           d__conic_R+iCoordSys*9);
Snormal.x = K;
Snormal.y = L;
Snormal.z = M;
forward_transform_centered(&Snormal, &Snormal,
			   d__conic_R+iCoordSys*9);
K = Snormal.x;
L = Snormal.y;
M = Snormal.z;

K += partial_x_bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, 
				     BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE,
				     BM_N_SAMPLE, BM_radius);
L += partial_y_bending_modes_surface(&xyz, N_mode, b + iCoordSys*N_mode, 
				     BM + iCoordSys*BM_N_SAMPLE*BM_N_SAMPLE,
				     BM_N_SAMPLE, BM_radius);
M += partial_z_bending_modes_surface();

#line 1220 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      d__ray[ij].optical_path_length += s1;
      d__ray[ij].n_iteration = j;
      break;
    }
  s0 = s1;
}
#line 1096 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
G2 = K*K + L*L + M*M;
a  = k*K + l*L + m*M;

a *= -2.0/G2;
klm.x += a*K;
klm.y += a*L;
klm.z += a*M;
#line 1077 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
klm_GS.x = klm.x;
klm_GS.y = klm.y;
klm_GS.z = klm.z;
xyz_GS.x = xyz.x;
xyz_GS.y = xyz.y;
xyz_GS.z = xyz.z;
#line 1121 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 // MOTION <<<
 backward_transform(&xyz_GS, &xyz_GS,
		    d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
 backward_transform_centered(&klm_GS, &klm_GS,
			     d__motion_R+iCoordSys*9);
#line 1128 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 // RIGID_BODY <<<
 backward_transform(&(d__ray[ij].coordinates), &xyz_GS,
		    d__rigid_body_R+iCoordSys*9,
		    d__rigid_body_origin+iCoordSys);
 backward_transform_centered(&(d__ray[ij].directions), &klm_GS,
			     d__rigid_body_R+iCoordSys*9);

#line 911 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	    d__ray[ij].optical_path_difference +=
         d__ray[ij].optical_path_length - d__chief_ray[iSource].optical_path_length;
      d__ray[ij].throughput *= reflectivity[iCoordSys];
	    return;
	  }
        }
      }
  }
}
#line 1562 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void m1_trace_all_kernel(char *mask, ray *d__ray, int N_RAY,
                                rtd inner2, rtd outer2,
                                rtd *d__aperture_R, vector *d__aperture_origin,
                                rtd *d__conic_R, vector *d__conic_origin,
                                rtd *d__Fk, rtd *d__Fc,
                                vector *d__conic_self_origin,
                                rtd *d__rigid_body_R, vector *d__rigid_body_origin,
                                rtd *d__motion_R, vector *d__motion_origin,
                                int max_n, int n_mode,
                                rtd *d__a, rtd *d__cx, rtd *d__cy, rtd R,
                                ray *d__chief_ray, int *d__piston_mask,
                                char *d__valid_segments, float *reflectivity)
{
  int i, j, ij, iCoordSys, iSource, b_n_mode;
  rtd rho2;
  rtd x, y, z, x1, y1, s0, k, l, m;
  rtd s1, S, K, L ,M, dSds;
  rtd G2, a;
  rtd Fc, Fk;
  vector d__origin, xyz, klm, xyz_GS, klm_GS, Snormal;

  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = threadIdx.y;
  //  iCoordSys = blockIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  iSource = blockIdx.z;
  j = ij;
  ij += iSource*N_RAY;
  b_n_mode = max_n*(max_n+1)*0.5;
  if  ( ( j<N_RAY ) && (d__ray[ij].v) )// & (iCoordSys<1) )
  {
    iCoordSys = 6;
    // for (iCoordSys=0; iCoordSys<7 ; iCoordSys++)
      {
	Fc = d__Fc[iCoordSys];
	Fk = d__Fk[iCoordSys];
        if (d__valid_segments[iCoordSys]>0) {
	d__origin = d__conic_self_origin[iCoordSys];

        
#line 1002 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID BODY >>>
forward_transform(&xyz_GS, &(d__ray[ij].coordinates),
		  d__rigid_body_R+iCoordSys*9, d__rigid_body_origin+iCoordSys);
forward_transform_centered(&klm_GS, &(d__ray[ij].directions),
			   d__rigid_body_R+iCoordSys*9);
#line 1009 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION >>>
forward_transform(&xyz_GS, &xyz_GS,
		  d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
forward_transform_centered(&klm_GS, &klm_GS,
			    d__motion_R+iCoordSys*9);

#line 1603 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
		  {
	    mask[ij] = 1;
	    d__piston_mask[ij] = iCoordSys + 1;

            
#line 1050 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
klm.x = klm_GS.x;
klm.y =	klm_GS.y;
klm.z =	klm_GS.z;
xyz.x = xyz_GS.x;
xyz.y =	xyz_GS.y;
xyz.z =	xyz_GS.z;
#line 1163 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k  = klm.x;
l  = klm.y;
m  = klm.z;
if (m==0) { return; }
s0 = -xyz.z/m;
x1 = xyz.x + k*s0;
y1 = xyz.y + l*s0;
d__ray[ij].optical_path_length = s0;
s0 = s1 = 0;
for (j=0; j<SNELL_N_ITERATION; j++)
{
  xyz.x = x1 + k*s0;
  xyz.y = y1 + l*s0;
  xyz.z = m*s0;
  
#line 1232 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC >>>
backward_transform_centered(&xyz, &xyz,
			    d__conic_R+iCoordSys*9);
S = conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
#line 1334 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
M = partial_z_conic_surface();
#line 1241 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC <<<
forward_transform_centered(&xyz, &xyz,
                           d__conic_R+iCoordSys*9);
S += zernike_surface(&xyz, &d__origin, max_n,
                     d__a + iCoordSys*n_mode, R); 
#line 1338 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
Snormal.x = K;
Snormal.y = L;
Snormal.z = M;
forward_transform_centered(&Snormal, &Snormal,
			   d__conic_R+iCoordSys*9);
K = Snormal.x;
L = Snormal.y;
M = Snormal.z;
K += partial_x_zernike_surface(&xyz, &d__origin, max_n,
                               d__cx + iCoordSys*b_n_mode, R);
L += partial_y_zernike_surface(&xyz, &d__origin, max_n,
                               d__cy + iCoordSys*b_n_mode, R);
M += partial_z_zernike_surface();
#line 1178 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  dSds = K*k + L*l + M*m;
  if (dSds==0) { break; }
  s1 = s0 - S/dSds;
  if (abs(s1-s0)<TOL)
    {
      xyz.x = x1 + k*s1;
      xyz.y = y1 + l*s1;
      xyz.z = m*s1;
      
#line 1276 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// CONIC >>>
backward_transform_centered(&xyz, &xyz,
			    d__conic_R+iCoordSys*9);
K = partial_x_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
L = partial_y_conic_surface(&xyz, &d__origin, Fk, Fc, 0, &(xyz.x));
M = partial_z_conic_surface();
// CONIC <<<
forward_transform_centered(&xyz, &xyz,
                           d__conic_R+iCoordSys*9);
Snormal.x = K;
Snormal.y = L;
Snormal.z = M;
forward_transform_centered(&Snormal, &Snormal,
			   d__conic_R+iCoordSys*9);
K = Snormal.x;
L = Snormal.y;
M = Snormal.z;
K += partial_x_zernike_surface(&xyz, &d__origin, max_n,
                               d__cx + iCoordSys*b_n_mode, R);
L += partial_y_zernike_surface(&xyz, &d__origin, max_n,
                               d__cy + iCoordSys*b_n_mode, R);
M += partial_z_zernike_surface();
#line 1187 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      d__ray[ij].optical_path_length += s1;
      d__ray[ij].n_iteration = j;
      break;
    }
  s0 = s1;
}
#line 1096 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
G2 = K*K + L*L + M*M;
a  = k*K + l*L + m*M;

a *= -2.0/G2;
klm.x += a*K;
klm.y += a*L;
klm.z += a*M;
#line 1058 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
klm_GS.x = klm.x;
klm_GS.y = klm.y;
klm_GS.z = klm.z;
xyz_GS.x = xyz.x;
xyz_GS.y = xyz.y;
xyz_GS.z = xyz.z;
#line 1121 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 // MOTION <<<
 backward_transform(&xyz_GS, &xyz_GS,
		    d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
 backward_transform_centered(&klm_GS, &klm_GS,
			     d__motion_R+iCoordSys*9);
#line 1128 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 // RIGID_BODY <<<
 backward_transform(&(d__ray[ij].coordinates), &xyz_GS,
		    d__rigid_body_R+iCoordSys*9,
		    d__rigid_body_origin+iCoordSys);
 backward_transform_centered(&(d__ray[ij].directions), &klm_GS,
			     d__rigid_body_R+iCoordSys*9);

#line 1608 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	    d__ray[ij].optical_path_difference +=
         d__ray[ij].optical_path_length - d__chief_ray[iSource].optical_path_length;
      d__ray[ij].throughput *= reflectivity[iCoordSys];
	    return;
	  }
        }
      }
  }
}
#line 669 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void m1_blocking_kernel(char *mask, ray *d__ray, int N_RAY,
				rtd inner2, rtd outer2,
				rtd *d__aperture_R, vector *d__aperture_origin,
				rtd *d__rigid_body_R, vector *d__rigid_body_origin,
				rtd *d__motion_R, vector *d__motion_origin)
{
  int i, j, ij, iCoordSys, iSource;
  rtd rho2;
  rtd x, y, z, x1, y1, s0, k, l, m;
  vector xyz, klm, xyz_GS, klm_GS;

  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = threadIdx.y;
  //  iCoordSys = blockIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  iSource = blockIdx.z;
  j = ij;
  ij += iSource*N_RAY;
  if  ( ( j<N_RAY ) && (d__ray[ij].v) )// & (iCoordSys<1) )
  {

    for (iCoordSys=0; iCoordSys<7; iCoordSys++)
      {

        
#line 1002 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID BODY >>>
forward_transform(&xyz_GS, &(d__ray[ij].coordinates),
		  d__rigid_body_R+iCoordSys*9, d__rigid_body_origin+iCoordSys);
forward_transform_centered(&klm_GS, &(d__ray[ij].directions),
			   d__rigid_body_R+iCoordSys*9);
#line 1009 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION >>>
forward_transform(&xyz_GS, &xyz_GS,
		  d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
forward_transform_centered(&klm_GS, &klm_GS,
			    d__motion_R+iCoordSys*9);

#line 695 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	
#line 1016 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// APERTURE >>>
forward_transform(&xyz, &xyz_GS,
		  d__aperture_R+iCoordSys*9,
		  d__aperture_origin+iCoordSys);
forward_transform_centered(&klm, &klm_GS,
			   d__aperture_R+iCoordSys*9);

	x = xyz.x;
	y = xyz.y;
	z = xyz.z;
	k = klm.x;
	l = klm.y;
	m = klm.z;

        
#line 1154 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (m==0) { return; }
s0 = -z/m;
x1 = x + k*s0;
y1 = y + l*s0;

#line 1032 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
        rho2 = x1*x1 + y1*y1;
	if ( (iCoordSys==6) && (rho2<inner2) )
	  {
	    return;
	  }

#line 697 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
	if (rho2<=outer2)
	  {
	    mask[ij] = 0;
          }
       }
   }
}
#line 1682 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void track_kernel(float *x, float *y, float *z, const int N,
			     const int iCoordSys,
			     rtd *d__rigid_body_R, vector *d__rigid_body_origin,
			     rtd *d__motion_R, vector *d__motion_origin)
 {
   int i;
   vector v;
   i = blockIdx.x * blockDim.x + threadIdx.x;
   if (i<N) {
     v.x = x[i];
     v.y = y[i];
     v.z = z[i];
     // MOTION <<<
     backward_transform(&v, &v,
                        d__motion_R+iCoordSys*9,
                        d__motion_origin+iCoordSys);
     // RIGID_BODY <<<
     backward_transform(&v, &v,
                        d__rigid_body_R+iCoordSys*9,
                        d__rigid_body_origin+iCoordSys);
     x[i] = v.x;
     y[i] = v.y;
     z[i] = v.z;
   }
}
#line 1730 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void locate_kernel(float *x, float *y, float *z, const int N,
			     const int iCoordSys,
			     rtd *d__rigid_body_R, vector *d__rigid_body_origin,
			     rtd *d__motion_R, vector *d__motion_origin)
 {
   int i;
   vector v;
   i = blockIdx.x * blockDim.x + threadIdx.x;
   if (i<N) {
     v.x = x[i];
     v.y = y[i];
     v.z = z[i];
     // RIGID_BODY >>>
     forward_transform(&v, &v,
                       d__rigid_body_R+iCoordSys*9,
                       d__rigid_body_origin+iCoordSys);
     // MOTION >>>
     forward_transform(&v, &v,
                       d__motion_R+iCoordSys*9,
                       d__motion_origin+iCoordSys);
     x[i] = v.x;
     y[i] = v.y;
     z[i] = v.z;
   }
}
#line 2665 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void nearest_neighbor_kernel(rtd *BMi, rtd *partial_x_BMi, rtd *partial_y_BMi,
					const int NI, const rtd di,
					double *x_BM, double *y_BM, double *BM,
					const int BM_N_SAMPLE, const double BM_radius,
					const int k_mode)
{
  int i, j, k, k_BM, t, o,
    t_E, t_N, t_W, t_S;
  rtd xi, yi, rho_t, rho_k, rhoi, h,
    xi_E, yi_N, xi_W, yi_S, 
    rho_k_E, rho_k_N, rho_k_W, rho_k_S, 
    rho_t_E, rho_t_N, rho_t_W, rho_t_S;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (i<NI) && (j<NI)) {
    k = i + j*NI;

    
#line 2703 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
xi = i - 0.5*(NI-1);
yi = j - 0.5*(NI-1);
xi *= di;
yi *= di;
rhoi = hypot(xi,yi);

h = 4*5e-2;

xi_E = xi + h;
yi_N = yi + h;
xi_W = xi - h;
yi_S = yi - h;

#line 2684 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    BMi[k] = 0;
    partial_x_BMi[k] = 0;
    partial_y_BMi[k] = 0;

    if (rhoi<=BM_radius) {
      
#line 2718 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
t = 0;
rho_t = hypot(x_BM[t]-xi , y_BM[t]-yi);

t_E = 0;
rho_t_E = hypot(x_BM[t_E]-xi_E, y_BM[t_E]-yi);
t_N = 0;
rho_t_N = hypot(x_BM[t_N]-xi  , y_BM[t_N]-yi_N);
t_W = 0;
rho_t_W = hypot(x_BM[t_W]-xi_W, y_BM[t_W]-yi);
t_S = 0;
rho_t_S = hypot(x_BM[t_S]-xi  , y_BM[t_S]-yi_S);
for (k_BM=0; k_BM<BM_N_SAMPLE; k_BM++) {
  rho_k = hypot(xi-x_BM[k_BM] , yi-y_BM[k_BM]);
  if (rho_k<rho_t) {
    t = k_BM;
    rho_t = rho_k;
  }
  rho_k_E = hypot(x_BM[k_BM]-xi_E, y_BM[k_BM]-yi);
  if (rho_k_E<rho_t_E) {
    t_E = k_BM;
    rho_t_E = rho_k_E;
  }
  rho_k_N = hypot(x_BM[k_BM]-xi, y_BM[k_BM]-yi_N);
  if (rho_k_N<rho_t_N) {
    t_N = k_BM;
    rho_t_N = rho_k_N;
  }
  rho_k_W = hypot(x_BM[k_BM]-xi_W, y_BM[k_BM]-yi);
  if (rho_k_W<rho_t_W) {
    t_W = k_BM;
    rho_t_W = rho_k_W;
  }
  rho_k_S = hypot(x_BM[k_BM]-xi, y_BM[k_BM]-yi_S);
  if (rho_k_S<rho_t_S) {
    t_S = k_BM;
    rho_t_S = rho_k_S;
  }
}
#line 2690 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      o = (k_mode-1)*BM_N_SAMPLE;
      BMi[k] = BM[t+o];
      partial_x_BMi[k] = BM[t_W+o] - BM[t_E+o];
      partial_x_BMi[k] /= (x_BM[t_W] - x_BM[t_E]);
      partial_y_BMi[k] = BM[t_N+o] - BM[t_S+o];
      partial_y_BMi[k] /= (y_BM[t_N] - y_BM[t_S]);
    }

  }
}
#line 2788 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void bilinear_kernel(rtd *BMi, rtd *partial_x_BMi, rtd *partial_y_BMi,
				const int NI, const rtd di,
				double *BM,
				const int BM_N_SAMPLE, const double BM_radius,
				const int k_mode)
{
  int i, j, k, ndx, o;
  rtd xi, yi, rhoi, s, t, scale, fs, ft, onemt, onems, h;
  rtd *z;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (i<NI) && (j<NI)) {
    k = i*NI + j;

    xi = i - 0.5*(NI-1);
    yi = j - 0.5*(NI-1);
    xi *= di;
    yi *= di;
    rhoi = hypot(xi,yi);

    BMi[k] = 0;
    partial_x_BMi[k] = 0;
    partial_y_BMi[k] = 0;

    if (rhoi<=4.181) {

      o = (k_mode-1)*BM_N_SAMPLE*BM_N_SAMPLE;
      z = BM + o;
      scale = 0.5*di/BM_radius;
      h = BM_radius*2/(BM_N_SAMPLE-1);

      
#line 2829 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1);
t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1);
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2832 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
BMi[k] = ( z[ndx]*onemt + z[ndx+1]*t )*onems +
  ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s;
#line 2820 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      
#line 2836 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1) + 1;
t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1);
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2839 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
partial_x_BMi[k] = ( z[ndx]*onemt + z[ndx+1]*t )*onems +
    ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s;
s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1) - 1;
t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1);
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2844 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
partial_x_BMi[k] -= ( z[ndx]*onemt + z[ndx+1]*t )*onems +
    ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s;
partial_x_BMi[k] /= 2*h;
#line 2821 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      
#line 2849 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1);
t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1) + 1;
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2852 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
partial_y_BMi[k] = ( z[ndx]*onemt + z[ndx+1]*t )*onems +
    ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s;
s = (scale*(i + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1);
t = (scale*(j + (1-NI)*0.5) + 0.5)*(BM_N_SAMPLE-1) - 1;
#line 2863 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
if (s<0) { s=0.0; }
fs  = floorf(s);
if (t<0) { t=0.0; }
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*BM_N_SAMPLE );

if (s==(BM_N_SAMPLE-1)) { fs -= 1; ndx -= BM_N_SAMPLE; }
if (t==(BM_N_SAMPLE-1)) { ft -= 1; ndx -= 1; }
s -= fs;
t -= ft;
onemt = 1 - t;
onems = 1 - s;
#line 2857 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
partial_y_BMi[k] -= ( z[ndx]*onemt + z[ndx+1]*t )*onems +
    ( z[ndx+BM_N_SAMPLE]*onemt + z[ndx+BM_N_SAMPLE+1]*t )*s;
partial_y_BMi[k] /= 2*h;

#line 2823 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
   }
  }
}
#line 3403 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void edge_sensors_read(vector *v, const int N_OUT_SEGMENT,
				   const rtd D, const rtd alpha,
				   rtd *d__rigid_body_R, vector *d__rigid_body_origin,
				   rtd *d__motion_R, vector *d__motion_origin)
{
  int i, iCoordSys, k;
  rtd R, a, offset_angle;
  i = threadIdx.x;
  if (i<N_OUT_SEGMENT)
  {
    R = 0.5*D;
    a = alpha/180;

    
#line 3427 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
iCoordSys = 6;
k = i;
offset_angle = (3.0-2.0*i)/6.0 + a;
#line 3466 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x       *= R;
v[k].y       *= R;
v[k].z       = 0.0;
#line 3473 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
#line 3479 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys*9,
		   d__rigid_body_origin+iCoordSys);
#line 3434 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k             += N_OUT_SEGMENT;
offset_angle = (3.0-2.0*i)/6.0 - a;
#line 3466 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x       *= R;
v[k].y       *= R;
v[k].z       = 0.0;
#line 3473 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
#line 3479 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys*9,
		   d__rigid_body_origin+iCoordSys);

#line 3418 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    __syncthreads();

    
#line 3440 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
iCoordSys    = i;
k            = 2*N_OUT_SEGMENT + i;
offset_angle = -0.5 + a;
#line 3466 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x       *= R;
v[k].y       *= R;
v[k].z       = 0.0;
#line 3473 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
#line 3479 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys*9,
		   d__rigid_body_origin+iCoordSys);
#line 3444 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k += 6;
offset_angle = -0.5 - a;
#line 3466 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x       *= R;
v[k].y       *= R;
v[k].z       = 0.0;
#line 3473 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
#line 3479 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys*9,
		   d__rigid_body_origin+iCoordSys);
#line 3449 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k += 6;
offset_angle = -1.0/6.0 + a;
#line 3466 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x       *= R;
v[k].y       *= R;
v[k].z       = 0.0;
#line 3473 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
#line 3479 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys*9,
		   d__rigid_body_origin+iCoordSys);
#line 3452 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k += 6;
offset_angle = -1.0/6.0 - a;
#line 3466 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x       *= R;
v[k].y       *= R;
v[k].z       = 0.0;
#line 3473 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
#line 3479 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys*9,
		   d__rigid_body_origin+iCoordSys);
#line 3457 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k += 6;
offset_angle = -5.0/6.0 + a;
#line 3466 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x       *= R;
v[k].y       *= R;
v[k].z       = 0.0;
#line 3473 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
#line 3479 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys*9,
		   d__rigid_body_origin+iCoordSys);
#line 3460 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k += 6;
offset_angle = -5.0/6.0 - a;
#line 3466 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x       *= R;
v[k].y       *= R;
v[k].z       = 0.0;
#line 3473 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys*9, d__motion_origin+iCoordSys);
#line 3479 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys*9,
		   d__rigid_body_origin+iCoordSys);

#line 3422 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  }
}
#line 3522 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void edge_sensors_data_kernel(vector *d__dv,
					vector *d__v,
					vector *d__dv0,
					const int N_OUT_SEGMENT)
{
  int k;
  k = threadIdx.x;
  if (k<N_OUT_SEGMENT)
    {
      d__dv[k]                 =
      	d__v[k+3*N_OUT_SEGMENT] -
        d__v[k] -
        d__dv0[k];
      d__dv[k+  N_OUT_SEGMENT] =
      	d__v[k+2*N_OUT_SEGMENT] -
        d__v[k+  N_OUT_SEGMENT] -
        d__dv0[k+N_OUT_SEGMENT];
      d__dv[k+2*N_OUT_SEGMENT] =
      	d__v[(k+1)%6+7*N_OUT_SEGMENT] -
        d__v[k+4*N_OUT_SEGMENT] -
        d__dv0[k+2*N_OUT_SEGMENT];
      d__dv[k+3*N_OUT_SEGMENT] =
      	d__v[(k+1)%6+6*N_OUT_SEGMENT] -
        d__v[k+5*N_OUT_SEGMENT] -
        d__dv0[k+3*N_OUT_SEGMENT];
    }
}
#line 3801 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void laser_coordinates(vector *v, const int N_OUT_SEGMENT,
				  const rtd D,
				  const rtd alpha,
				  const rtd height,
				  rtd *d__rigid_body_R, vector *d__rigid_body_origin,
				  rtd *d__motion_R, vector *d__motion_origin)
{
  int i, j, iCoordSys_laser, iCoordSys_cam, k;
  rtd R, a, offset_angle;
  i = threadIdx.x;
  j = threadIdx.y;
  if (i<N_OUT_SEGMENT)
  {
    R = 0.5*D;
    a = alpha/180;
    
#line 3823 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
iCoordSys_laser = 6;
iCoordSys_cam   = i;
k = i + 24*j;
offset_angle = (3.0-2.0*i)/6.0 + a;
#line 3843 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x      *= R;
v[k].y      *= R;
v[k].z       = height;
#line 3849 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v[k].z      *= (j==0) ? 1 : -1;
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys_laser*9,
                   d__motion_origin+iCoordSys_laser);
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys_laser*9,
		   d__rigid_body_origin+iCoordSys_laser);
// RIGID BODY >>>
forward_transform(v+k, v+k,
		  d__rigid_body_R+iCoordSys_cam*9,
                  d__rigid_body_origin+iCoordSys_cam);
// MOTION >>>
forward_transform(v+k, v+k,
		  d__motion_R+iCoordSys_cam*9,
                  d__motion_origin+iCoordSys_cam);
#line 3828 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k             += N_OUT_SEGMENT;
offset_angle = (3.0-2.0*i)/6.0 - a;
#line 3843 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x      *= R;
v[k].y      *= R;
v[k].z       = height;
#line 3849 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v[k].z      *= (j==0) ? 1 : -1;
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys_laser*9,
                   d__motion_origin+iCoordSys_laser);
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys_laser*9,
		   d__rigid_body_origin+iCoordSys_laser);
// RIGID BODY >>>
forward_transform(v+k, v+k,
		  d__rigid_body_R+iCoordSys_cam*9,
                  d__rigid_body_origin+iCoordSys_cam);
// MOTION >>>
forward_transform(v+k, v+k,
		  d__motion_R+iCoordSys_cam*9,
                  d__motion_origin+iCoordSys_cam);
#line 3817 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    __syncthreads();
    
#line 3833 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
iCoordSys_laser = i;
iCoordSys_cam   = (i+1)%N_OUT_SEGMENT;
k = 2*N_OUT_SEGMENT + i + 24*j;
offset_angle = -1.0/6.0 + a;
#line 3843 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x      *= R;
v[k].y      *= R;
v[k].z       = height;
#line 3849 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v[k].z      *= (j==0) ? 1 : -1;
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys_laser*9,
                   d__motion_origin+iCoordSys_laser);
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys_laser*9,
		   d__rigid_body_origin+iCoordSys_laser);
// RIGID BODY >>>
forward_transform(v+k, v+k,
		  d__rigid_body_R+iCoordSys_cam*9,
                  d__rigid_body_origin+iCoordSys_cam);
// MOTION >>>
forward_transform(v+k, v+k,
		  d__motion_R+iCoordSys_cam*9,
                  d__motion_origin+iCoordSys_cam);
#line 3838 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k += 6;
offset_angle = -1.0/6.0 - a;
#line 3843 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x      *= R;
v[k].y      *= R;
v[k].z       = height;
#line 3849 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v[k].z      *= (j==0) ? 1 : -1;
// MOTION <<<
backward_transform(v+k, v+k,
		   d__motion_R+iCoordSys_laser*9,
                   d__motion_origin+iCoordSys_laser);
// RIGID_BODY <<<
backward_transform(v+k, v+k,
		   d__rigid_body_R+iCoordSys_laser*9,
		   d__rigid_body_origin+iCoordSys_laser);
// RIGID BODY >>>
forward_transform(v+k, v+k,
		  d__rigid_body_R+iCoordSys_cam*9,
                  d__rigid_body_origin+iCoordSys_cam);
// MOTION >>>
forward_transform(v+k, v+k,
		  d__motion_R+iCoordSys_cam*9,
                  d__motion_origin+iCoordSys_cam);
#line 3819 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  }
}
#line 3868 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void laser_coordinates_to_GCS(vector *v0, vector *v, const int N_OUT_SEGMENT,
				  rtd *d__rigid_body_R, vector *d__rigid_body_origin)
{
  int i, j, k;
  i = threadIdx.x;
  j = threadIdx.y;
  if (i<N_OUT_SEGMENT)
  {
    k = i + 24*j;
    // RIGID_BODY <<<
    backward_transform(v0+k, v+k,
                       d__rigid_body_R+i*9,
                       d__rigid_body_origin+i);
    k += N_OUT_SEGMENT;
    // RIGID_BODY <<<
    backward_transform(v0+k, v+k,
                       d__rigid_body_R+i*9,
                       d__rigid_body_origin+i);
    __syncthreads();
    k = 2*N_OUT_SEGMENT + i + 24*j;
    i = (i+1)%N_OUT_SEGMENT;
    // RIGID_BODY <<<
    backward_transform(v0+k, v+k,
                       d__rigid_body_R+i*9,
                       d__rigid_body_origin+i);
    k += N_OUT_SEGMENT;
    // RIGID_BODY <<<
    backward_transform(v0+k, v+k,
                       d__rigid_body_R+i*9,
                       d__rigid_body_origin+i);
  }
}
#line 3913 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void camera_coordinates(vector *v, vector *v0,
                                   const int N_OUT_SEGMENT,
				   const rtd D,
				   const rtd alpha,
				   const rtd height,
				   rtd *d__rigid_body_R,
                                   vector *d__rigid_body_origin)
{
  int i, j, k;
  rtd R, a, offset_angle;
  i = threadIdx.x;
  j = threadIdx.y;
  if (i<N_OUT_SEGMENT)
  {
    R = 0.5*D;
    a = alpha/180;
    
#line 3936 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k = i + 24*j;
offset_angle = -0.5 - a;
#line 3843 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x      *= R;
v[k].y      *= R;
v[k].z       = height;
#line 3953 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v[k].z      *= (j==0) ? -1 : 1;
// RIGID_BODY <<<
backward_transform(v0+k, v+k,
		   d__rigid_body_R+i*9,
		   d__rigid_body_origin+i);
#line 3939 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k += N_OUT_SEGMENT;
offset_angle = -0.5 + a;
#line 3843 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x      *= R;
v[k].y      *= R;
v[k].z       = height;
#line 3953 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v[k].z      *= (j==0) ? -1 : 1;
// RIGID_BODY <<<
backward_transform(v0+k, v+k,
		   d__rigid_body_R+i*9,
		   d__rigid_body_origin+i);
#line 3930 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    __syncthreads();
    
#line 3944 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k = 2*N_OUT_SEGMENT + i + 24*j;
offset_angle = -5.0/6.0 - a;
#line 3843 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x      *= R;
v[k].y      *= R;
v[k].z       = height;
#line 3953 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v[k].z      *= (j==0) ? -1 : 1;
// RIGID_BODY <<<
backward_transform(v0+k, v+k,
		   d__rigid_body_R+i*9,
		   d__rigid_body_origin+i);
#line 3947 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
k += N_OUT_SEGMENT;
offset_angle = -5.0/6.0 + a;
#line 3843 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sincospi( offset_angle, &(v[k].y), &(v[k].x) );
v[k].x      *= R;
v[k].y      *= R;
v[k].z       = height;
#line 3953 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v[k].z      *= (j==0) ? -1 : 1;
// RIGID_BODY <<<
backward_transform(v0+k, v+k,
		   d__rigid_body_R+i*9,
		   d__rigid_body_origin+i);
#line 3932 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  }
}
#line 3992 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void edge_sensors_camera_vector_kernel(vector *d__k,
                                                vector *d__A, vector*d__B,
                                                vector *d__uP, vector *d__vP,
                                                const int N_OUT_SEGMENT)
{
  int i, j, k;
  i = threadIdx.x;
  j = 24*threadIdx.y;
  vector n;
  n.x = n.y = 0.0;
  n.z = 1.0;
  if (i<N_OUT_SEGMENT)
    {
      k = i + j;
      d__k[k] = d__B[k] - d__A[k];
      
#line 4026 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
d__k[k].left_cross(&d__uP[k],&n);
d__uP[k].unit();
d__uP[k].left_cross(&d__vP[k], &d__k[k]);
d__vP[k].unit();

#line 4009 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      k += N_OUT_SEGMENT;
      d__k[k] =	d__B[k] - d__A[k];
      
#line 4026 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
d__k[k].left_cross(&d__uP[k],&n);
d__uP[k].unit();
d__uP[k].left_cross(&d__vP[k], &d__k[k]);
d__vP[k].unit();

#line 4013 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      __syncthreads();

      k = 2*N_OUT_SEGMENT + i + j;
      d__k[k] =	d__B[k-i+(i+1)%6] - d__A[k];
      
#line 4026 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
d__k[k].left_cross(&d__uP[k],&n);
d__uP[k].unit();
d__uP[k].left_cross(&d__vP[k], &d__k[k]);
d__vP[k].unit();

#line 4019 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      k += N_OUT_SEGMENT;
      d__k[k] =	d__B[k-i+(i+1)%6] - d__A[k];
      
#line 4026 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
d__k[k].left_cross(&d__uP[k],&n);
d__uP[k].unit();
d__uP[k].left_cross(&d__vP[k], &d__k[k]);
d__vP[k].unit();
#line 4022 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    }
}
#line 4052 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void edge_sensors_laser_vector_kernel(vector *d__k_laser,
                                                vector *d__k_cam,
                                                const int N_OUT_SEGMENT,
						 rtd *d__rigid_body_R,
						 rtd *d__motion_R)
{
  int i, j, iCoordSys_laser, iCoordSys_cam, k;
  i = threadIdx.x;
  j = 24*threadIdx.y;
  if (i<N_OUT_SEGMENT)
  {
    iCoordSys_laser = 6;
    iCoordSys_cam   = i;
    k = i + j;
    
#line 4080 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform_centered(d__k_laser+k, d__k_cam+k,
			    d__rigid_body_R+iCoordSys_cam*9);
// RIGID BODY >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__rigid_body_R+iCoordSys_laser*9);
// MOTION <<<
backward_transform_centered(d__k_laser+k, d__k_laser+k,
			    d__motion_R+iCoordSys_laser*9);
// RIGID_BODY <<<
backward_transform_centered(d__k_laser+k, d__k_laser+k,
			    d__rigid_body_R+iCoordSys_laser*9);
// RIGID BODY >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__rigid_body_R+iCoordSys_cam*9);
// MOTION >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__motion_R+iCoordSys_cam*9);
#line 4067 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    k             += N_OUT_SEGMENT;
    
#line 4080 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform_centered(d__k_laser+k, d__k_cam+k,
			    d__rigid_body_R+iCoordSys_cam*9);
// RIGID BODY >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__rigid_body_R+iCoordSys_laser*9);
// MOTION <<<
backward_transform_centered(d__k_laser+k, d__k_laser+k,
			    d__motion_R+iCoordSys_laser*9);
// RIGID_BODY <<<
backward_transform_centered(d__k_laser+k, d__k_laser+k,
			    d__rigid_body_R+iCoordSys_laser*9);
// RIGID BODY >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__rigid_body_R+iCoordSys_cam*9);
// MOTION >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__motion_R+iCoordSys_cam*9);
#line 4069 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    __syncthreads();
    iCoordSys_laser = i;
    iCoordSys_cam   = (i+1)%N_OUT_SEGMENT;
    k = 2*N_OUT_SEGMENT + i + j;
    
#line 4080 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform_centered(d__k_laser+k, d__k_cam+k,
			    d__rigid_body_R+iCoordSys_cam*9);
// RIGID BODY >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__rigid_body_R+iCoordSys_laser*9);
// MOTION <<<
backward_transform_centered(d__k_laser+k, d__k_laser+k,
			    d__motion_R+iCoordSys_laser*9);
// RIGID_BODY <<<
backward_transform_centered(d__k_laser+k, d__k_laser+k,
			    d__rigid_body_R+iCoordSys_laser*9);
// RIGID BODY >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__rigid_body_R+iCoordSys_cam*9);
// MOTION >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__motion_R+iCoordSys_cam*9);
#line 4074 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    k             += N_OUT_SEGMENT;
    
#line 4080 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
// RIGID_BODY <<<
backward_transform_centered(d__k_laser+k, d__k_cam+k,
			    d__rigid_body_R+iCoordSys_cam*9);
// RIGID BODY >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__rigid_body_R+iCoordSys_laser*9);
// MOTION <<<
backward_transform_centered(d__k_laser+k, d__k_laser+k,
			    d__motion_R+iCoordSys_laser*9);
// RIGID_BODY <<<
backward_transform_centered(d__k_laser+k, d__k_laser+k,
			    d__rigid_body_R+iCoordSys_laser*9);
// RIGID BODY >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__rigid_body_R+iCoordSys_cam*9);
// MOTION >>>
forward_transform_centered(d__k_laser+k, d__k_laser+k,
			   d__motion_R+iCoordSys_cam*9);
#line 4076 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  }
}
#line 4152 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  __global__ void laser_camera_intersection(rtd *x, rtd *y, rtd *d,
					  vector *B, vector *A,
					  vector *k_cam, vector *k_laser,
					  vector *uP, vector *vP,
					  const int N_OUT_SEGMENT)
{
  int i, j, i_cam, i_laser;
  rtd t;
  vector w, Ilp;
  i = threadIdx.x;
  j = 24*threadIdx.y;
  if (i<N_OUT_SEGMENT)
    {
      i_cam = i_laser = i + j;
      
#line 4185 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
w          = B[i_cam] - A[i_laser];
d[i_laser] = w.norm();
t          = (k_cam[i_cam]*w)/(k_cam[i_cam]*k_laser[i_laser]);
Ilp        = k_laser[i_laser]*t + A[i_laser];
w          = Ilp - B[i_cam];
x[i_cam]   = w*uP[i_cam];
y[i_cam]   = w*vP[i_cam];

#line 4168 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      i_cam   += N_OUT_SEGMENT;
      i_laser += N_OUT_SEGMENT;
      
#line 4185 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
w          = B[i_cam] - A[i_laser];
d[i_laser] = w.norm();
t          = (k_cam[i_cam]*w)/(k_cam[i_cam]*k_laser[i_laser]);
Ilp        = k_laser[i_laser]*t + A[i_laser];
w          = Ilp - B[i_cam];
x[i_cam]   = w*uP[i_cam];
y[i_cam]   = w*vP[i_cam];

#line 4172 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      __syncthreads();

      i_cam   = 2*N_OUT_SEGMENT + (i+1)%6 + j;
      i_laser = 2*N_OUT_SEGMENT + i + j;
      
#line 4185 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
w          = B[i_cam] - A[i_laser];
d[i_laser] = w.norm();
t          = (k_cam[i_cam]*w)/(k_cam[i_cam]*k_laser[i_laser]);
Ilp        = k_laser[i_laser]*t + A[i_laser];
w          = Ilp - B[i_cam];
x[i_cam]   = w*uP[i_cam];
y[i_cam]   = w*vP[i_cam];

#line 4178 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      i_cam   += N_OUT_SEGMENT;
      i_laser += N_OUT_SEGMENT;
      
#line 4185 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
w          = B[i_cam] - A[i_laser];
d[i_laser] = w.norm();
t          = (k_cam[i_cam]*w)/(k_cam[i_cam]*k_laser[i_laser]);
Ilp        = k_laser[i_laser]*t + A[i_laser];
w          = Ilp - B[i_cam];
x[i_cam]   = w*uP[i_cam];
y[i_cam]   = w*vP[i_cam];
#line 4181 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    }
}
#line 616 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void m1_preset_kernel(ray *d__ray, int N_RAY,
				 rtd L, int N_L, rtd R)
{
  int i, j, k, iSource;
  rtd x, y, r;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  iSource = blockIdx.z;
  if ( (i<N_L) && (j<N_L) ) {
    k = i + j*N_L + iSource*N_RAY;
    x   = L*(i - (N_L-1)*0.5)/(N_L-1);
    y   = L*(j - (N_L-1)*0.5)/(N_L-1);
    r = hypot(x,y);
    if (r>R)
      d__ray[k].v = 0;
  }
}

#line 2398 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void surface_update(double *BMS, const int N_mode, rtd *b,
					     double *BM, int N_SAMPLE)
{
  int i, j, k, kz, l, k_mode, o;
  rtd *z;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  l = blockIdx.z;
  if ( (i<N_SAMPLE) && (j<N_SAMPLE)) {
    k = kz = i*N_SAMPLE + j;
    k += l*N_SAMPLE*N_SAMPLE;
    BMS[k] = 0.0;
    for (k_mode=0; k_mode<N_mode; k_mode++) {
      o = k_mode*N_SAMPLE*N_SAMPLE;
      z = BM + o;
      BMS[k] += b[k_mode + l*N_mode]*z[kz];
    }
  }
}
#line 2300 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::setup(int _n_mode_)
{
  N = 1;
  n_mode = _n_mode_;
  
#line 2355 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte = sizeof(rtd)*n_mode*N;
b = (rtd *)malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__b, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__b, b,
              n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemset(d__b, 0, n_byte ) );
#line 2305 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  load_reg();
  
#line 2362 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte = BM_N_SAMPLE*BM_N_SAMPLE*N*sizeof(double);
HANDLE_ERROR( cudaMalloc((void**)&d__BMS, n_byte ) );
HANDLE_ERROR( cudaMemset(d__BMS, 0, n_byte ) );
#line 2307 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2311 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::setup(int _n_mode_, int N_SET)
{
  N = N_SET;
  n_mode = _n_mode_;
  
#line 2355 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte = sizeof(rtd)*n_mode*N;
b = (rtd *)malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__b, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__b, b,
              n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemset(d__b, 0, n_byte ) );
#line 2316 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  load_reg();
  
#line 2362 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte = BM_N_SAMPLE*BM_N_SAMPLE*N*sizeof(double);
HANDLE_ERROR( cudaMalloc((void**)&d__BMS, n_byte ) );
HANDLE_ERROR( cudaMemset(d__BMS, 0, n_byte ) );
#line 2318 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
} 
#line 2322 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::setupKL(int _n_mode_)
{
  N = 1;
  n_mode = _n_mode_;
  
#line 2355 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte = sizeof(rtd)*n_mode*N;
b = (rtd *)malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__b, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__b, b,
              n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemset(d__b, 0, n_byte ) );
#line 2327 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  load_KL();
  
#line 2362 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte = BM_N_SAMPLE*BM_N_SAMPLE*N*sizeof(double);
HANDLE_ERROR( cudaMalloc((void**)&d__BMS, n_byte ) );
HANDLE_ERROR( cudaMemset(d__BMS, 0, n_byte ) );
#line 2329 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2333 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::setupKL(int _n_mode_, int N_SET)
{
  N = N_SET;
  n_mode = _n_mode_;
  
#line 2355 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte = sizeof(rtd)*n_mode*N;
b = (rtd *)malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__b, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__b, b,
              n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemset(d__b, 0, n_byte ) );
#line 2338 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  load_KL();
  
#line 2362 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte = BM_N_SAMPLE*BM_N_SAMPLE*N*sizeof(double);
HANDLE_ERROR( cudaMalloc((void**)&d__BMS, n_byte ) );
HANDLE_ERROR( cudaMemset(d__BMS, 0, n_byte ) );
#line 2340 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
} 
#line 2344 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::setupPolish(int _n_mode_, int N_SET)
{
  N = N_SET;
  n_mode = _n_mode_;
  
#line 2355 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte = sizeof(rtd)*n_mode*N;
b = (rtd *)malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__b, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__b, b,
              n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemset(d__b, 0, n_byte ) );
#line 2349 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  load_polish();
  
#line 2362 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte = BM_N_SAMPLE*BM_N_SAMPLE*N*sizeof(double);
HANDLE_ERROR( cudaMalloc((void**)&d__BMS, n_byte ) );
HANDLE_ERROR( cudaMemset(d__BMS, 0, n_byte ) );
#line 2351 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
} 
#line 2369 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::
#line 2371 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
cleanup(void)
{
  INFO("@(CEO)>bending_modes: freeing memory!\n");
  free( b );
  HANDLE_ERROR( cudaFree( d__b ) );
  HANDLE_ERROR( cudaFree( d__BM_buffer ) );
  HANDLE_ERROR( cudaFree( d__BMS ) );
}
#line 2385 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::update(rtd *b)
{
  int n_byte = sizeof(rtd)*n_mode*N;
  HANDLE_ERROR( cudaMemcpy( d__b, b,
			    n_byte, cudaMemcpyHostToDevice ) );
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(BM_N_SAMPLE/N_THREAD+1,BM_N_SAMPLE/N_THREAD+1,N);
  
  surface_update <<< gridDim, blockDim >>> 
    (d__BMS, n_mode, d__b, d__BM, BM_N_SAMPLE);
}
#line 2429 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::load(void)
{
  FILE * pFile;
  long lSize;
  double * buffer;
  size_t result;

  char path[256]; 
  sprintf(path, "%sgmtMirrors/bendingModes.bin",getenv("CEOPATH"));
  INFO("@(CEO)>bending_modes: loading bending modes from %s\n",path);

  pFile = fopen ( path , "rb" );
  if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

  // obtain file size:
  fseek (pFile , 0 , SEEK_END);
  lSize = ftell (pFile)/sizeof(double);
  rewind (pFile);

  //  fprintf(stdout,"lSize=%ld\n",lSize);

  // allocate memory to contain the whole file:
  buffer = (double*) malloc (lSize*sizeof(double));
  if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);}

  // copy the file into the buffer:
  result = fread (buffer,sizeof(double),lSize,pFile);
  //  fprintf(stdout,"result=%ld\n",result);
  if (result != lSize) {fputs ("Reading error",stderr); exit (3);}

  /* the whole file is now loaded in the memory buffer. */
  BM_radius = 4.181;
  BM_N_SAMPLE = 27685;
  HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) );
  HANDLE_ERROR( cudaMemcpy( d__BM_buffer  , buffer, lSize*sizeof(double),
			    cudaMemcpyHostToDevice ) );
  d__x_BM = d__BM_buffer;
  d__y_BM = d__BM_buffer + BM_N_SAMPLE;
  d__BM   = d__BM_buffer + BM_N_SAMPLE*2;

  // terminate
  fclose (pFile);
  free (buffer);
}
#line 2483 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::load_reg(void)
{
  FILE * pFile;
  long lSize;
  double * buffer;
  size_t result;

  char path[256]; 
  sprintf(path, "%sgmtMirrors/bendingModesReg.bin",getenv("CEOPATH"));
  INFO("@(CEO)>bending_modes: loading bending modes from %s\n",path);

  pFile = fopen ( path , "rb" );
  if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

  // obtain file size:
  fseek (pFile , 0 , SEEK_END);
  lSize = ftell (pFile)/sizeof(double);
  rewind (pFile);

  //  fprintf(stdout,"lSize=%ld\n",lSize);

  // allocate memory to contain the whole file:
  buffer = (double*) malloc (lSize*sizeof(double));
  if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);}

  // copy the file into the buffer:
  result = fread (buffer,sizeof(double),lSize,pFile);
  //  fprintf(stdout,"result=%ld\n",result);
  if (result != lSize) {fputs ("Reading error",stderr); exit (3);}

  /* the whole file is now loaded in the memory buffer. */
  BM_radius = 0.5*8.462;//4.181;
  BM_N_SAMPLE = 169;
  HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) );
  HANDLE_ERROR( cudaMemcpy( d__BM_buffer  , buffer, lSize*sizeof(double),
			    cudaMemcpyHostToDevice ) );
  d__BM   = d__BM_buffer;

  // terminate
  fclose (pFile);
  free (buffer);
}
#line 2535 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::load_polish(void)
{
  FILE * pFile;
  long lSize;
  double * buffer;
  size_t result;

  char path[256]; 
  sprintf(path, "%sgmtMirrors/polishingMap.bin",getenv("CEOPATH"));
  INFO("@(CEO)>bending_modes: loading bending modes from %s\n",path);

  pFile = fopen ( path , "rb" );
  if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

  // obtain file size:
  fseek (pFile , 0 , SEEK_END);
  lSize = ftell (pFile)/sizeof(double);
  rewind (pFile);

  //  fprintf(stdout,"lSize=%ld\n",lSize);

  // allocate memory to contain the whole file:
  buffer = (double*) malloc (lSize*sizeof(double));
  if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);}

  // copy the file into the buffer:
  result = fread (buffer,sizeof(double),lSize,pFile);
  //  fprintf(stdout,"result=%ld\n",result);
  if (result != lSize) {fputs ("Reading error",stderr); exit (3);}

  /* the whole file is now loaded in the memory buffer. */
  BM_radius = 0.5*8.52;
  BM_N_SAMPLE = 483;
  HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) );
  HANDLE_ERROR( cudaMemcpy( d__BM_buffer  , buffer, lSize*sizeof(double),
			    cudaMemcpyHostToDevice ) );
  d__BM   = d__BM_buffer;

  // terminate
  fclose (pFile);
  free (buffer);
}
#line 2587 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::load_KL(void)
{
  FILE * pFile;
  long lSize;
  double * buffer;
  size_t result;

  char path[256]; 
  sprintf(path, "%sgmtMirrors/KarhunenLoeveModes.bin",getenv("CEOPATH"));
  INFO("@(CEO)>bending_modes: loading bending modes from %s\n",path);

  pFile = fopen ( path , "rb" );
  if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

  // obtain file size:
  fseek (pFile , 0 , SEEK_END);
  lSize = ftell (pFile)/sizeof(double);
  rewind (pFile);

  //  fprintf(stdout,"lSize=%ld\n",lSize);

  // allocate memory to contain the whole file:
  buffer = (double*) malloc (lSize*sizeof(double));
  if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);}

  // copy the file into the buffer:
  result = fread (buffer,sizeof(double),lSize,pFile);
  //  fprintf(stdout,"result=%ld\n",result);
  if (result != lSize) {fputs ("Reading error",stderr); exit (3);}

  /* the whole file is now loaded in the memory buffer. */
  BM_radius = 0.5*1.0425*1.1;//4.25;
  BM_N_SAMPLE = 256;
  HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) );
  HANDLE_ERROR( cudaMemcpy( d__BM_buffer  , buffer, lSize*sizeof(double),
			    cudaMemcpyHostToDevice ) );
  d__BM   = d__BM_buffer;

  // terminate
  fclose (pFile);
  free (buffer);
}
#line 2647 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::nearest_neighbor(rtd *d__BMi,
			      rtd *d__partial_x_BMi,
			      rtd *d__partial_y_BMi,
			      int NI, rtd di, int k_mode)
{
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(NI/N_THREAD+1,NI/N_THREAD+1);
  nearest_neighbor_kernel <<< gridDim, blockDim >>> (d__BMi,
						     d__partial_x_BMi,
						     d__partial_y_BMi,
						     NI, di,
						     d__x_BM, d__y_BM,
						     d__BM, BM_N_SAMPLE,
						     BM_radius, k_mode);
}
#line 2771 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void bending_modes::bilinear(rtd *d__BMi,
			     rtd *d__partial_x_BMi,
			     rtd *d__partial_y_BMi,
			     int NI, rtd di, int k_mode)
{
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(NI/N_THREAD+1,NI/N_THREAD+1);
  bilinear_kernel <<< gridDim, blockDim >>> (d__BMi,
					     d__partial_x_BMi,
					     d__partial_y_BMi,
					     NI, di,
					     d__BM, BM_N_SAMPLE,
					     BM_radius, k_mode);
}
#line 3210 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
__global__ void modes2surface(double *BMS, const int n_mode,
                              rtd *b, int *s2b,
                              double *BM, const int N_MODE,
                              const int N_SAMPLE)
{
  int i, j, k, kz, l, k_mode, o;
  rtd *z;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  l = blockIdx.z;
  if ( (i<N_SAMPLE) && (j<N_SAMPLE)) {
    k = kz = i*N_SAMPLE + j;
    k += l*N_SAMPLE*N_SAMPLE;
    BMS[k] = 0.0;
    for (k_mode=0; k_mode<n_mode; k_mode++) {
      o = (s2b[l]*N_MODE+k_mode)*N_SAMPLE*N_SAMPLE;
      z = BM + o;
      BMS[k] += b[k_mode + l*n_mode]*z[kz];
    }
  }
}
#line 3041 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void modes::setup(char *_filename_, int _N_, int _n_mode_)
{
  N = _N_;
  n_mode = _n_mode_;
  strcpy(filename,_filename_);
  load();
  
#line 2355 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte = sizeof(rtd)*n_mode*N;
b = (rtd *)malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__b, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__b, b,
              n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemset(d__b, 0, n_byte ) );
#line 3048 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 2362 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte = BM_N_SAMPLE*BM_N_SAMPLE*N*sizeof(double);
HANDLE_ERROR( cudaMalloc((void**)&d__BMS, n_byte ) );
HANDLE_ERROR( cudaMemset(d__BMS, 0, n_byte ) );
#line 3049 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
} 
void modes::setup(int _BM_N_SAMPLE_, double _BM_radius_,
                  int _N_SET_, int _N_MODE_, int *s2b,
                  double *buffer, int _N_, int _n_mode_)
{
  N = _N_;
  n_mode = _n_mode_;
  load(_BM_N_SAMPLE_, _BM_radius_,
       _N_SET_, _N_MODE_, s2b,
       buffer);
  
#line 2355 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte = sizeof(rtd)*n_mode*N;
b = (rtd *)malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__b, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__b, b,
              n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemset(d__b, 0, n_byte ) );
#line 3060 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 2362 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte = BM_N_SAMPLE*BM_N_SAMPLE*N*sizeof(double);
HANDLE_ERROR( cudaMalloc((void**)&d__BMS, n_byte ) );
HANDLE_ERROR( cudaMemset(d__BMS, 0, n_byte ) );
#line 3061 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
} 
#line 3066 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void modes::
#line 2371 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
cleanup(void)
{
  INFO("@(CEO)>bending_modes: freeing memory!\n");
  free( b );
  HANDLE_ERROR( cudaFree( d__b ) );
  HANDLE_ERROR( cudaFree( d__BM_buffer ) );
  HANDLE_ERROR( cudaFree( d__BMS ) );
}
#line 3085 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void modes::load(void)
{
  FILE * pFile;
  long lSize;
  double * buffer;
  int *s2b;
  size_t result;

  //char path[256]; 
  //  sprintf(path, "%sgmtMirrors/%s.ceo",getenv("CEOPATH"),filename);
  INFO("@(CEO)>modes: loading modes from %s\n",filename);

  pFile = fopen ( filename, "rb" );
  if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

  // read modal basis info
  result = fread (&BM_N_SAMPLE,sizeof(int)   ,1,pFile);
  result = fread (&BM_radius  ,sizeof(double),1,pFile);
  result = fread (&N_SET,sizeof(int)   ,1,pFile);
  result = fread (&N_MODE,sizeof(int)   ,1,pFile);
  INFO("@(CEO)>modes: grid sampling and size: %d pixel, %f meter\n",BM_N_SAMPLE,BM_radius);
  INFO("@(CEO)>modes: number of set: %d and number of modes per set: %d\n",N_SET,N_MODE);

  // read surface 2 mode array
  s2b = (int*) malloc (N*sizeof(int));
  result = fread (s2b  ,sizeof(int),N,pFile);
  HANDLE_ERROR( cudaMalloc((void**)&d__s2b, N*sizeof(int)) );
  HANDLE_ERROR( cudaMemcpy( d__s2b  , s2b, N*sizeof(int),
			    cudaMemcpyHostToDevice ) );
  

  // allocate memory to contain the whole file:
  lSize = N_SET*N_MODE*BM_N_SAMPLE*BM_N_SAMPLE;
  buffer = (double*) malloc (lSize*sizeof(double));
  if (buffer == NULL) {fputs ("Memory error",stderr); exit (2);}

  // copy the file into the buffer:
  result = fread (buffer  ,sizeof(double),lSize,pFile);
  if (result != lSize) {fputs ("Reading error",stderr); exit (3);}

  /* the whole file is now loaded in the memory buffer. */
  BM_radius *= 0.5;
  HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) );
  HANDLE_ERROR( cudaMemcpy( d__BM_buffer  , buffer, lSize*sizeof(double),
			    cudaMemcpyHostToDevice ) );
  d__BM   = d__BM_buffer;

  // terminate
  fclose (pFile);
  free(s2b);
  free (buffer);
}
#line 3145 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void modes::load(int _BM_N_SAMPLE_, double _BM_radius_,
                 int _N_SET_, int _N_MODE_, int *s2b,
                 double *buffer)
{
  long lSize;

  BM_N_SAMPLE = _BM_N_SAMPLE_;
  BM_radius   = _BM_radius_;
  N_SET       = _N_SET_;
  N_MODE      = _N_MODE_;

  HANDLE_ERROR( cudaMalloc((void**)&d__s2b, N*sizeof(int)) );
  HANDLE_ERROR( cudaMemcpy( d__s2b  , s2b, N*sizeof(int),
			    cudaMemcpyHostToDevice ) );

  lSize = N_SET*N_MODE*BM_N_SAMPLE*BM_N_SAMPLE;

  BM_radius *= 0.5;
  HANDLE_ERROR( cudaMalloc((void**)&d__BM_buffer, lSize*sizeof(double)) );
  HANDLE_ERROR( cudaMemcpy( d__BM_buffer  , buffer, lSize*sizeof(double),
			    cudaMemcpyHostToDevice ) );
  d__BM   = d__BM_buffer;

}
#line 3177 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void modes::reset_modes(double *buffer)
{
  long lSize = N_SET*N_MODE*BM_N_SAMPLE*BM_N_SAMPLE;
  HANDLE_ERROR( cudaMemcpy( d__BM_buffer  , buffer, lSize*sizeof(double),
			    cudaMemcpyHostToDevice ) );
  d__BM   = d__BM_buffer;
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(BM_N_SAMPLE/N_THREAD+1,BM_N_SAMPLE/N_THREAD+1,N);
  
  modes2surface <<< gridDim, blockDim >>> 
    (d__BMS, n_mode, d__b, d__s2b, d__BM, N_MODE, BM_N_SAMPLE);

}
#line 3197 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void modes::update(rtd *b)
{
  int n_byte = sizeof(rtd)*n_mode*N;
  HANDLE_ERROR( cudaMemcpy( d__b, b,
			    n_byte, cudaMemcpyHostToDevice ) );
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(BM_N_SAMPLE/N_THREAD+1,BM_N_SAMPLE/N_THREAD+1,N);
  
  modes2surface <<< gridDim, blockDim >>> 
    (d__BMS, n_mode, d__b, d__s2b, d__BM, N_MODE, BM_N_SAMPLE);
}
#line 256 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::setup(void)
{
  //  BS = NULL;
  ZS = NULL;
  
#line 293 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
M_ID    = 1;
D_assembly = 25.498;
D_full  = 8.417;
D_clear = 8.365;
area    = 357;
area0   = 357;
ri      = 2.875/8.417;
beta    = 13.601685*PI/180.0;
L       = 8.710;
N       = 7;
conic_c = 1.0/36.0;
conic_k = 1-0.9982857;
height  = 0.0;
idx_offset = 0;
V = NULL;
#line 261 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 310 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte;
float segment_reflectivity[7] = {1.,1.,1.,1.,1.,1.,1.};
n_byte = 7*sizeof(float);
HANDLE_ERROR( cudaMalloc((void**)&d__segment_reflectivity, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, segment_reflectivity,
                          n_byte, cudaMemcpyHostToDevice ) );

rtd *h__conic_c;
n_byte = 7*sizeof(rtd);
h__conic_c = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_c[k] = conic_c;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_c, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_c, h__conic_c,
              n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_c);

rtd *h__conic_k;
n_byte = 7*sizeof(rtd);
h__conic_k = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_k[k] = conic_k;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_k, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_k, h__conic_k,
                          n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_k);

area_fraction = 1.0;
vector __v0(0.0,0.0,0.0),
       __v(D_full*0.5,0.0,0.0),
       origin[12],
       euler_angles[12];
rtd D_c, o, zo;
 int k, idx;
char tag[16];

depth = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);
D_c = L;
__v.x = L;
zo = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);

HANDLE_ERROR( cudaMalloc((void**)&d__valid_segments, sizeof(char)*7) );
HANDLE_ERROR( cudaMemset( d__valid_segments, 1, sizeof(char)*7) );

#line 385 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = height;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
for (k=0; k<N-1; k++)
{
  o = PI*(3-2*k)/6.0;
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 394 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = D_c*cos(o);
  origin[idx].y = D_c*sin(o);
  origin[idx].z = height + zo;
  if (M_ID==1) {
    euler_angles[idx].x = beta;
    euler_angles[idx].y = 0.0;
    euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }
}
sprintf(tag,"M%d RIGID BODY",M_ID);
rigid_body_CS.setup(origin, euler_angles, N, tag);

#line 421 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = (M_ID==1) ? depth: -depth;
euler_angles[N-1].x = euler_angles[N-1].y = euler_angles[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 427 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = (M_ID==1) ? depth: -depth;
  euler_angles[idx].x = 0.0;
  euler_angles[idx].y = 0.0;
  euler_angles[idx].z = 0.0;
}
sprintf(tag,"M%d APERTURE",M_ID);
aperture_CS.setup(origin, euler_angles, N, tag);

#line 463 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = origin[N-1].z = 0.0;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
conic_origin[N-1].x = conic_origin[N-1].y = conic_origin[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 471 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  o = PI*(3-2*k)/6.0;
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = 0.0;
   if (M_ID==1) {
     euler_angles[idx].x = beta;
     euler_angles[idx].y = 0.0;
     euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }     
  conic_origin[idx].x = -D_c*cos(o);
  conic_origin[idx].y = -D_c*sin(o);
  conic_origin[idx].z = -zo;
}
//euler_angles[0].x = 30.0*PI/180.0/3600.0;
sprintf(tag,"M%d CONIC",M_ID);
conic_CS.setup(origin, euler_angles, N, tag);
HANDLE_ERROR( cudaMalloc((void**)&d__conic_origin, sizeof(vector)*N ) );
HANDLE_ERROR( cudaMemcpy( d__conic_origin, conic_origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );

#line 503 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
for (k=0; k<N; k++)
{
  origin[k].x = 0.0;
  origin[k].y = 0.0;
  origin[k].z = 0.0;
  euler_angles[k].x = 0.0;
  euler_angles[k].y = 0.0;
  euler_angles[k].z = 0.0;
}
sprintf(tag,"M%d MOTION",M_ID);
motion_CS.setup(origin, euler_angles, N, tag);

#line 1894 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sprintf(tag,"GLOBAL TIP-TILT");
vector TT_CS_origin, TT_CS_euler_angles;
TT_CS_origin.x = TT_CS_origin.y = TT_CS_origin.z = 0.0;
TT_CS_euler_angles.x = TT_CS_euler_angles.y = TT_CS_euler_angles.z = 0.0;
TT_CS.setup(&TT_CS_origin, &TT_CS_euler_angles, 1, tag);
cublasCreate(&handle);
HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) );
HANDLE_ERROR( cudaMalloc((void**)&d__D, sizeof(double)*9 ) );
#line 262 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 274 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::setup(char *_filename_, int _N_, int _n_mode_)
{
  BS.setup(_filename_, _N_, _n_mode_);
  ZS = NULL;
  
#line 293 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
M_ID    = 1;
D_assembly = 25.498;
D_full  = 8.417;
D_clear = 8.365;
area    = 357;
area0   = 357;
ri      = 2.875/8.417;
beta    = 13.601685*PI/180.0;
L       = 8.710;
N       = 7;
conic_c = 1.0/36.0;
conic_k = 1-0.9982857;
height  = 0.0;
idx_offset = 0;
V = NULL;
#line 279 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 310 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte;
float segment_reflectivity[7] = {1.,1.,1.,1.,1.,1.,1.};
n_byte = 7*sizeof(float);
HANDLE_ERROR( cudaMalloc((void**)&d__segment_reflectivity, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, segment_reflectivity,
                          n_byte, cudaMemcpyHostToDevice ) );

rtd *h__conic_c;
n_byte = 7*sizeof(rtd);
h__conic_c = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_c[k] = conic_c;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_c, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_c, h__conic_c,
              n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_c);

rtd *h__conic_k;
n_byte = 7*sizeof(rtd);
h__conic_k = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_k[k] = conic_k;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_k, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_k, h__conic_k,
                          n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_k);

area_fraction = 1.0;
vector __v0(0.0,0.0,0.0),
       __v(D_full*0.5,0.0,0.0),
       origin[12],
       euler_angles[12];
rtd D_c, o, zo;
 int k, idx;
char tag[16];

depth = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);
D_c = L;
__v.x = L;
zo = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);

HANDLE_ERROR( cudaMalloc((void**)&d__valid_segments, sizeof(char)*7) );
HANDLE_ERROR( cudaMemset( d__valid_segments, 1, sizeof(char)*7) );

#line 385 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = height;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
for (k=0; k<N-1; k++)
{
  o = PI*(3-2*k)/6.0;
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 394 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = D_c*cos(o);
  origin[idx].y = D_c*sin(o);
  origin[idx].z = height + zo;
  if (M_ID==1) {
    euler_angles[idx].x = beta;
    euler_angles[idx].y = 0.0;
    euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }
}
sprintf(tag,"M%d RIGID BODY",M_ID);
rigid_body_CS.setup(origin, euler_angles, N, tag);

#line 421 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = (M_ID==1) ? depth: -depth;
euler_angles[N-1].x = euler_angles[N-1].y = euler_angles[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 427 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = (M_ID==1) ? depth: -depth;
  euler_angles[idx].x = 0.0;
  euler_angles[idx].y = 0.0;
  euler_angles[idx].z = 0.0;
}
sprintf(tag,"M%d APERTURE",M_ID);
aperture_CS.setup(origin, euler_angles, N, tag);

#line 463 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = origin[N-1].z = 0.0;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
conic_origin[N-1].x = conic_origin[N-1].y = conic_origin[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 471 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  o = PI*(3-2*k)/6.0;
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = 0.0;
   if (M_ID==1) {
     euler_angles[idx].x = beta;
     euler_angles[idx].y = 0.0;
     euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }     
  conic_origin[idx].x = -D_c*cos(o);
  conic_origin[idx].y = -D_c*sin(o);
  conic_origin[idx].z = -zo;
}
//euler_angles[0].x = 30.0*PI/180.0/3600.0;
sprintf(tag,"M%d CONIC",M_ID);
conic_CS.setup(origin, euler_angles, N, tag);
HANDLE_ERROR( cudaMalloc((void**)&d__conic_origin, sizeof(vector)*N ) );
HANDLE_ERROR( cudaMemcpy( d__conic_origin, conic_origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );

#line 503 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
for (k=0; k<N; k++)
{
  origin[k].x = 0.0;
  origin[k].y = 0.0;
  origin[k].z = 0.0;
  euler_angles[k].x = 0.0;
  euler_angles[k].y = 0.0;
  euler_angles[k].z = 0.0;
}
sprintf(tag,"M%d MOTION",M_ID);
motion_CS.setup(origin, euler_angles, N, tag);

#line 1894 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sprintf(tag,"GLOBAL TIP-TILT");
vector TT_CS_origin, TT_CS_euler_angles;
TT_CS_origin.x = TT_CS_origin.y = TT_CS_origin.z = 0.0;
TT_CS_euler_angles.x = TT_CS_euler_angles.y = TT_CS_euler_angles.z = 0.0;
TT_CS.setup(&TT_CS_origin, &TT_CS_euler_angles, 1, tag);
cublasCreate(&handle);
HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) );
HANDLE_ERROR( cudaMalloc((void**)&d__D, sizeof(double)*9 ) );
#line 280 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 284 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::setup( zernikeS *_ZS_)
{
  ZS = _ZS_;
  //BS = NULL;
  
#line 293 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
M_ID    = 1;
D_assembly = 25.498;
D_full  = 8.417;
D_clear = 8.365;
area    = 357;
area0   = 357;
ri      = 2.875/8.417;
beta    = 13.601685*PI/180.0;
L       = 8.710;
N       = 7;
conic_c = 1.0/36.0;
conic_k = 1-0.9982857;
height  = 0.0;
idx_offset = 0;
V = NULL;
#line 289 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 310 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte;
float segment_reflectivity[7] = {1.,1.,1.,1.,1.,1.,1.};
n_byte = 7*sizeof(float);
HANDLE_ERROR( cudaMalloc((void**)&d__segment_reflectivity, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, segment_reflectivity,
                          n_byte, cudaMemcpyHostToDevice ) );

rtd *h__conic_c;
n_byte = 7*sizeof(rtd);
h__conic_c = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_c[k] = conic_c;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_c, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_c, h__conic_c,
              n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_c);

rtd *h__conic_k;
n_byte = 7*sizeof(rtd);
h__conic_k = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_k[k] = conic_k;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_k, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_k, h__conic_k,
                          n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_k);

area_fraction = 1.0;
vector __v0(0.0,0.0,0.0),
       __v(D_full*0.5,0.0,0.0),
       origin[12],
       euler_angles[12];
rtd D_c, o, zo;
 int k, idx;
char tag[16];

depth = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);
D_c = L;
__v.x = L;
zo = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);

HANDLE_ERROR( cudaMalloc((void**)&d__valid_segments, sizeof(char)*7) );
HANDLE_ERROR( cudaMemset( d__valid_segments, 1, sizeof(char)*7) );

#line 385 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = height;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
for (k=0; k<N-1; k++)
{
  o = PI*(3-2*k)/6.0;
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 394 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = D_c*cos(o);
  origin[idx].y = D_c*sin(o);
  origin[idx].z = height + zo;
  if (M_ID==1) {
    euler_angles[idx].x = beta;
    euler_angles[idx].y = 0.0;
    euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }
}
sprintf(tag,"M%d RIGID BODY",M_ID);
rigid_body_CS.setup(origin, euler_angles, N, tag);

#line 421 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = (M_ID==1) ? depth: -depth;
euler_angles[N-1].x = euler_angles[N-1].y = euler_angles[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 427 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = (M_ID==1) ? depth: -depth;
  euler_angles[idx].x = 0.0;
  euler_angles[idx].y = 0.0;
  euler_angles[idx].z = 0.0;
}
sprintf(tag,"M%d APERTURE",M_ID);
aperture_CS.setup(origin, euler_angles, N, tag);

#line 463 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = origin[N-1].z = 0.0;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
conic_origin[N-1].x = conic_origin[N-1].y = conic_origin[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 471 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  o = PI*(3-2*k)/6.0;
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = 0.0;
   if (M_ID==1) {
     euler_angles[idx].x = beta;
     euler_angles[idx].y = 0.0;
     euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }     
  conic_origin[idx].x = -D_c*cos(o);
  conic_origin[idx].y = -D_c*sin(o);
  conic_origin[idx].z = -zo;
}
//euler_angles[0].x = 30.0*PI/180.0/3600.0;
sprintf(tag,"M%d CONIC",M_ID);
conic_CS.setup(origin, euler_angles, N, tag);
HANDLE_ERROR( cudaMalloc((void**)&d__conic_origin, sizeof(vector)*N ) );
HANDLE_ERROR( cudaMemcpy( d__conic_origin, conic_origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );

#line 503 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
for (k=0; k<N; k++)
{
  origin[k].x = 0.0;
  origin[k].y = 0.0;
  origin[k].z = 0.0;
  euler_angles[k].x = 0.0;
  euler_angles[k].y = 0.0;
  euler_angles[k].z = 0.0;
}
sprintf(tag,"M%d MOTION",M_ID);
motion_CS.setup(origin, euler_angles, N, tag);

#line 1894 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sprintf(tag,"GLOBAL TIP-TILT");
vector TT_CS_origin, TT_CS_euler_angles;
TT_CS_origin.x = TT_CS_origin.y = TT_CS_origin.z = 0.0;
TT_CS_euler_angles.x = TT_CS_euler_angles.y = TT_CS_euler_angles.z = 0.0;
TT_CS.setup(&TT_CS_origin, &TT_CS_euler_angles, 1, tag);
cublasCreate(&handle);
HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) );
HANDLE_ERROR( cudaMalloc((void**)&d__D, sizeof(double)*9 ) );
#line 290 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 523 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::cleanup(void)
{
  INFO("@(CEO)>gmt_m1: freeing memory!\n");
  
#line 530 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
INFO(" |-");
BS.cleanup();
INFO(" |-");
#line 1904 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
TT_CS.cleanup();
cublasDestroy(handle);
HANDLE_ERROR( cudaFree( d__C ) );
HANDLE_ERROR( cudaFree( d__D ) );
#line 534 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
INFO(" |-");
aperture_CS.cleanup();
INFO(" |-");
conic_CS.cleanup();
INFO(" |-");
rigid_body_CS.cleanup();
INFO(" |-");
motion_CS.cleanup();
HANDLE_ERROR( cudaFree( d__conic_origin ) );
HANDLE_ERROR( cudaFree( d__valid_segments ) );
HANDLE_ERROR( cudaFree( d__conic_c ) );
HANDLE_ERROR( cudaFree( d__conic_k ) );
HANDLE_ERROR( cudaFree( d__segment_reflectivity ) );
#line 527 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 551 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::update_conic_c(rtd *_conic_c_)
{
  HANDLE_ERROR( cudaMemcpy( d__conic_c, _conic_c_,
              sizeof(rtd)*7, cudaMemcpyHostToDevice ) );
}
#line 560 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  void gmt_m1::update_conic_k(rtd *_conic_k_)
{
  HANDLE_ERROR( cudaMemcpy( d__conic_k, _conic_k_,
                            sizeof(rtd)*7, cudaMemcpyHostToDevice ) );
}
#line 569 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::update(vector _origin_, vector _euler_angles_,int idx)
{
    
#line 575 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  //fprintf(stdout,"\n\x1B[31m");
motion_CS.update( _origin_, _euler_angles_, --idx);
//motion_CS.info();
#line 572 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 582 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::reset(void)
{
    
#line 588 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
vector zero;
zero.x = 0.0;
zero.y = 0.0;
zero.z = 0.0;
//fprintf(stdout,"\n\x1B[31m");
for (int k=0; k<N; k++)
   motion_CS.update( zero, zero, k);
//motion_CS.info();
#line 585 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 606 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::preset(bundle *rays, rtd margin)
{
 dim3 blockDim(N_THREAD,N_THREAD);
 dim3 gridDim(rays->N_L/N_THREAD+1,rays->N_L/N_THREAD+1, rays->N_BUNDLE);
 m1_preset_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
					      rays->L, rays->N_L,
					      0.5*(margin+D_assembly));
}
#line 640 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::blocking(bundle *rays)
{
  
#line 646 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 rtd R2, Rri2;
 R2 = D_clear*D_clear*0.25;
 Rri2 = R2*ri*ri;
 V = &(rays->V);
 // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2);
 // printf("N_BUNDLE=%d\n",rays->N_BUNDLE);
 int nel = rays->N_RAY*rays->N_BUNDLE;
 dim3 blockDim(N_THREAD2);
 dim3 gridDim(nel/N_THREAD2+1);
 fill_ones_char <<< gridDim,blockDim >>> (V->m,nel);
 blockDim = dim3(N_THREAD,N_THREAD);
 gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE);
 m1_blocking_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
					     Rri2, R2,
					     aperture_CS.d__R,
					     aperture_CS.d__origin,
					     rigid_body_CS.d__R,
					     rigid_body_CS.d__origin,
					     motion_CS.d__R,
					     motion_CS.d__origin);
 intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY);
#line 643 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 718 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::trace(bundle *rays)
{
  rays->V.area = area0*rays->N_BUNDLE*area_fraction;
  dim3 blockDim, gridDim;
  rtd R2, Rri2;
  R2 = D_clear*D_clear*0.25;
  Rri2 = R2*ri*ri;
  
#line 729 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 /* if (V!=NULL) { */
 /*    blockDim = dim3(N_THREAD,N_THREAD); */
 /*    gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); */
 /*    intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY); */
 /* } */
 V = &(rays->V);
 // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2);
 blockDim = dim3(1,1);
 gridDim  = dim3(1,1,rays->N_BUNDLE);
 m1_trace_chief_kernel <<< gridDim , blockDim >>> (V->m, rays->d__chief_ray, 1,
						   conic_CS.d__R,
						   conic_CS.d__origin,
						   conic_k, conic_c,
						   d__conic_origin,
						   rigid_body_CS.d__R,
						   rigid_body_CS.d__origin,
						   motion_CS.d__R,
						   motion_CS.d__origin);
 blockDim = dim3(N_THREAD,N_THREAD);
 gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE);
 HANDLE_ERROR( cudaMemset(V->m, 0, sizeof(char)*rays->N_RAY*rays->N_BUNDLE ) );
 HANDLE_ERROR( cudaMemset(rays->d__piston_mask, 0,
			  sizeof(int)*rays->N_RAY*rays->N_BUNDLE ) );
/*if (ZS!=NULL)
  m1_trace_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
                                              Rri2, R2,
                                              aperture_CS.d__R,
                                              aperture_CS.d__origin,
                                              conic_CS.d__R,
                                              conic_CS.d__origin,
                                              d__conic_k, d__conic_c,
                                              d__conic_origin,
                                              rigid_body_CS.d__R,
                                              rigid_body_CS.d__origin,
                                              motion_CS.d__R,
                                              motion_CS.d__origin,
                                              ZS->max_n, ZS->n_mode,
                                              ZS->d__a, ZS->d__cx, ZS->d__cy,
                                              0.5*D_full,
                                              rays->d__chief_ray,
                                              rays->d__piston_mask,
                                              d__valid_segments,
                                              d__segment_reflectivity);

					      if (BS!=NULL)*/
  m1_bm_trace_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
                                                 Rri2, R2,
                                                 aperture_CS.d__R,
                                                 aperture_CS.d__origin,
                                                 conic_CS.d__R,
                                                 conic_CS.d__origin,
                                                 d__conic_k, d__conic_c,
                                                 d__conic_origin,
                                                 rigid_body_CS.d__R,
                                                 rigid_body_CS.d__origin,
                                                 motion_CS.d__R,
                                                 motion_CS.d__origin,
                                                 BS.n_mode, BS.d__b,
                                                 BS.d__BMS,
                                                 BS.BM_N_SAMPLE, BS.BM_radius,
                                                 rays->d__chief_ray,
                                                 rays->d__piston_mask,
                                                 d__valid_segments,
                                                 d__segment_reflectivity);
intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY);
float previous_nnz = rays->V.nnz;
//fprintf(stdout,"M%d previous nnz rays: %f\n",M_ID,previous_nnz);
rays->V.set_filter_quiet();
//fprintf(stdout,"M%d current nnz rays: %f\n",M_ID,rays->V.nnz);
if (previous_nnz<rays->N_RAY_TOTAL)
  rays->V.area *= rays->V.nnz/previous_nnz;
//fprintf(stdout,"M%d current area: %f\n",M_ID,rays->V.area);
#line 726 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 952 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::test_ray_tracing(void)
{
  int iCoordSys;
  vector in, out;

  in.x = 8.71*cos(PI/6);
  in.y = 8.71*sin(PI/6);
  in.z = 0.0;

  printf("Output vector:\n");
  for (iCoordSys = 0; iCoordSys<N; iCoordSys++)
    {
      printf(" #%d\n",iCoordSys);
        forward_transform(&out, &in,
                          rigid_body_CS.R+iCoordSys*9,
                          rigid_body_CS.origin+iCoordSys);
        printf(". x = %.2e ; y = %.2e ; z = %.2e\n",out.x,out.y,out.z);

        forward_transform(&out, &out,
                          motion_CS.R+iCoordSys*9,
                          motion_CS.origin+iCoordSys);
        printf(". x = %.2e ; y = %.2e ; z = %.2e\n",out.x,out.y,out.z);

	backward_transform_centered(&out, &out,
				    conic_CS.R+iCoordSys*9);
        printf(". x = %.2e ; y = %.2e ; z = %.2e\n",out.x,out.y,out.z);

	forward_transform_centered(&out, &out,
				    conic_CS.R+iCoordSys*9);
        printf(". x = %.2e ; y = %.2e ; z = %.2e\n",out.x,out.y,out.z);

        backward_transform(&out, &out,
                          motion_CS.R+iCoordSys*9,
                          motion_CS.origin+iCoordSys);
        printf(". x = %.2e ; y = %.2e ; z = %.2e\n",out.x,out.y,out.z);

        backward_transform(&out, &out,
                          rigid_body_CS.R+iCoordSys*9,
                          rigid_body_CS.origin+iCoordSys);
        printf(". x = %.2e ; y = %.2e ; z = %.2e\n",out.x-in.x,out.y-in.y,out.z-in.z);
    }

}
#line 1474 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::traceall(bundle *rays)
{
  rays->V.area = area0*rays->N_BUNDLE*area_fraction;
  dim3 blockDim, gridDim;
  rtd R2, Rri2;
  R2 = D_clear*D_clear*0.25;
  Rri2 = R2*ri*ri*0;
  
#line 1485 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 /* if (V!=NULL) { */
 /*    blockDim = dim3(N_THREAD,N_THREAD); */
 /*    gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); */
 /*    intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY); */
 /* } */
 V = &(rays->V);
 // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2);
 blockDim = dim3(1,1);
 gridDim  = dim3(1,1,rays->N_BUNDLE);
 m1_trace_chief_kernel <<< gridDim , blockDim >>> (V->m, rays->d__chief_ray, 1,
						   conic_CS.d__R,
						   conic_CS.d__origin,
						   conic_k, conic_c,
						   d__conic_origin,
						   rigid_body_CS.d__R,
						   rigid_body_CS.d__origin,
						   motion_CS.d__R,
						   motion_CS.d__origin);
 blockDim = dim3(N_THREAD,N_THREAD);
 gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE);
/*HANDLE_ERROR( cudaMemset(V->m, 0, sizeof(char)*rays->N_RAY*rays->N_BUNDLE ) );
 HANDLE_ERROR( cudaMemset(rays->d__piston_mask, 0,
			  sizeof(int)*rays->N_RAY*rays->N_BUNDLE ) );
*/
if (ZS!=NULL)
  m1_trace_all_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
                                              Rri2, R2,
                                              aperture_CS.d__R,
                                              aperture_CS.d__origin,
                                              conic_CS.d__R,
                                              conic_CS.d__origin,
                                              d__conic_k, d__conic_c,
                                              d__conic_origin,
                                              rigid_body_CS.d__R,
                                              rigid_body_CS.d__origin,
                                              motion_CS.d__R,
                                              motion_CS.d__origin,
                                              ZS->max_n, ZS->n_mode,
                                              ZS->d__a, ZS->d__cx, ZS->d__cy,
                                              0.5*D_full,
                                              rays->d__chief_ray,
                                              rays->d__piston_mask,
                                              d__valid_segments,
                                              d__segment_reflectivity);
/*
if (BS!=NULL)
  m1_bm_trace_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
                                                 Rri2, R2,
                                                 aperture_CS.d__R,
                                                 aperture_CS.d__origin,
                                                 conic_CS.d__R,
                                                 conic_CS.d__origin,
                                                 d__conic_k, d__conic_c,
                                                 d__conic_origin,
                                                 rigid_body_CS.d__R,
                                                 rigid_body_CS.d__origin,
                                                 motion_CS.d__R,
                                                 motion_CS.d__origin,
                                                 BS.n_mode, BS.d__b,
                                                 BS.d__BMS,
                                                 BS.BM_N_SAMPLE, BS.BM_radius,
                                                 rays->d__chief_ray,
                                                 rays->d__piston_mask,
                                                 d__valid_segments,
                                                 d__segment_reflectivity);

intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY);
float previous_nnz = rays->V.nnz;
//fprintf(stdout,"M%d previous nnz rays: %f\n",M_ID,previous_nnz);
rays->V.set_filter_quiet();
//fprintf(stdout,"M%d current nnz rays: %f\n",M_ID,rays->V.nnz);
if (previous_nnz<rays->N_RAY_TOTAL)
  rays->V.area *= rays->V.nnz/previous_nnz;
//fprintf(stdout,"M%d current area: %f\n",M_ID,rays->V.area);
*/
#line 1482 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 1624 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::remove(int *seg_ID, int N_ID)
{
  
#line 1631 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int k;
char valid[7];
area_fraction = (7-N_ID)/7.0;
memset(valid, 1, sizeof(char)*7);
for (k=0;k<N_ID;k++)
  valid[seg_ID[k]-1] = 0;
HANDLE_ERROR( cudaMemcpy( d__valid_segments, valid,
                          sizeof(char)*7, cudaMemcpyHostToDevice ) );
#line 1627 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 1643 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::keep(int *seg_ID, int N_ID)
{
  
#line 1650 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int k;
char valid[7];
area_fraction = N_ID/7.0;
memset(valid, 0, sizeof(char)*7);
for (k=0;k<N_ID;k++)
  valid[seg_ID[k]-1] = 1;
HANDLE_ERROR( cudaMemcpy( d__valid_segments, valid,
                          sizeof(char)*7, cudaMemcpyHostToDevice ) );
#line 1646 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 1666 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::track(float *d__x, float *d__y, float *d__z, int N, int idx)
{
  
#line 1672 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
dim3 blockDim(N_THREAD,1);
dim3 gridDim(N/N_THREAD+1,1);
track_kernel <<< gridDim, blockDim >>> (d__x, d__y, d__z, N,
					idx,
					rigid_body_CS.d__R,
					rigid_body_CS.d__origin,
					motion_CS.d__R,
					motion_CS.d__origin);
#line 1669 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 1714 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::locate(float *d__x, float *d__y, float *d__z, int N, int idx)
{
  
#line 1720 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
dim3 blockDim(N_THREAD,1);
dim3 gridDim(N/N_THREAD+1,1);
locate_kernel <<< gridDim, blockDim >>> (d__x, d__y, d__z, N,
					idx,
					rigid_body_CS.d__R,
					rigid_body_CS.d__origin,
					motion_CS.d__R,
					motion_CS.d__origin);
#line 1717 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 1888 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::global_tiptilt(float tip, float tilt)
{
  
#line 1910 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
double alpha, beta;
alpha = 1.0;
beta  = 0.0;
int k, k9;
vector origin, euler_angles, o;
origin.x = origin.y = origin.z = 0.0;
euler_angles.x = tip;
euler_angles.y = tilt;
euler_angles.z = 0.0;
TT_CS.update(origin,euler_angles,0);
for (k=0;k<N;k++)
  {
    k9 = k*9;
    
#line 1794 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_T,
	     3,3,3,
	     &alpha,
	     TT_CS.d__R, 3,
	     rigid_body_CS.d__R+k9, 3,
	     &beta,
	     d__C, 3);
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_N,
	     3,3,3,
	     &alpha,
	     rigid_body_CS.d__R+k9, 3,
	     d__C, 3,
	     &beta,
	     d__D, 3);
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_N,
	     3,3,3,
	     &alpha,
	     motion_CS.d__R+k9, 3,
	     d__D, 3,
	     &beta,
	     d__C, 3);

 HANDLE_ERROR( cudaMemcpy( motion_CS.d__R+k9, d__C,
			   sizeof(double)*9, cudaMemcpyDeviceToDevice ) );
 HANDLE_ERROR( cudaMemcpy( motion_CS.R+k9, motion_CS.d__R+k9,
			   sizeof(double)*9, cudaMemcpyDeviceToHost ) );
 /*
 int i, j, idx = 0;
for (i=0;i<3;i++)
   {
     fprintf(stdout,"||  ");
     for (j=0;j<3;j++)
       {
	 idx = j + i*3 + k9;
	 fprintf(stdout,"%+.4e  ",motion_CS.R[idx]);
       }
     fprintf(stdout,"||\n");
   }
 */
#line 1845 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
motion_CS.euler_angles[k].x = atan2(motion_CS.R[7+k9],motion_CS.R[8+k9]);
motion_CS.euler_angles[k].y = asin(-1.0*motion_CS.R[6+k9]);
motion_CS.euler_angles[k].z = atan2(motion_CS.R[3+k9],motion_CS.R[k9]);
#line 1856 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
vector o;
forward_transform(&o, rigid_body_CS.origin + k,
		  TT_CS.R, TT_CS.origin);
o.x = rigid_body_CS.origin[k].x - TT_CS.origin->x - o.x;
o.y = rigid_body_CS.origin[k].y - TT_CS.origin->y - o.y;
o.z = rigid_body_CS.origin[k].z - TT_CS.origin->z - o.z;
#line 1865 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
forward_transform_centered(&o, &o, rigid_body_CS.R + k9);
motion_CS.origin[k].x += o.x; 
motion_CS.origin[k].y += o.y;
motion_CS.origin[k].z += o.z;
HANDLE_ERROR( cudaMemcpy( motion_CS.d__origin+k, motion_CS.origin+k,
			  sizeof(vector), cudaMemcpyHostToDevice ) );

/*
 fprintf(stdout,"#%d >> MOTION CS origins [micron]  : %+.3e, %+.3e, %+.3e\n",k,
	 motion_CS.origin[k].x*1e6,
	 motion_CS.origin[k].y*1e6,
	 motion_CS.origin[k].z*1e6);

 float r2d;
 r2d = 1000*3600*180.0/PI;
 fprintf(stdout,"#%d >> MOTION CS Euler angles [mas]: %+.3e, %+.3e, %+.3e\n",k,
	 r2d*motion_CS.euler_angles[k].x,
	 r2d*motion_CS.euler_angles[k].y,
	 r2d*motion_CS.euler_angles[k].z);
*/
#line 1924 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  }
#line 1891 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 1932 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m1::set_reflectivity(float *reflectivity)
{
  
#line 1938 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte;
n_byte = 7*sizeof(float);
HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, reflectivity,
                          n_byte, cudaMemcpyHostToDevice ) );
#line 1935 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 1961 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::setup(void)
{
  //  BS = NULL;
  ZS = NULL;
  
#line 1988 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 M_ID    = 2;
 D_assembly = 3.168;
 D_full  = 1.0425*1.005; // slightly larger to always encompass M1; needed for the Zernike
 D_clear = 1.0415;
 ri      = 0.0;
 beta    = -14.777498*PI/180.0;
 L       = 1.08774;
 N       = 7;
 conic_c = -1.0/4.1639009;
 conic_k = 1-0.71692784;
 height  = 20.26247614;
 idx_offset = 3;
#line 1966 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 310 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte;
float segment_reflectivity[7] = {1.,1.,1.,1.,1.,1.,1.};
n_byte = 7*sizeof(float);
HANDLE_ERROR( cudaMalloc((void**)&d__segment_reflectivity, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, segment_reflectivity,
                          n_byte, cudaMemcpyHostToDevice ) );

rtd *h__conic_c;
n_byte = 7*sizeof(rtd);
h__conic_c = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_c[k] = conic_c;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_c, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_c, h__conic_c,
              n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_c);

rtd *h__conic_k;
n_byte = 7*sizeof(rtd);
h__conic_k = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_k[k] = conic_k;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_k, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_k, h__conic_k,
                          n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_k);

area_fraction = 1.0;
vector __v0(0.0,0.0,0.0),
       __v(D_full*0.5,0.0,0.0),
       origin[12],
       euler_angles[12];
rtd D_c, o, zo;
 int k, idx;
char tag[16];

depth = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);
D_c = L;
__v.x = L;
zo = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);

HANDLE_ERROR( cudaMalloc((void**)&d__valid_segments, sizeof(char)*7) );
HANDLE_ERROR( cudaMemset( d__valid_segments, 1, sizeof(char)*7) );

#line 385 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = height;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
for (k=0; k<N-1; k++)
{
  o = PI*(3-2*k)/6.0;
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 394 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = D_c*cos(o);
  origin[idx].y = D_c*sin(o);
  origin[idx].z = height + zo;
  if (M_ID==1) {
    euler_angles[idx].x = beta;
    euler_angles[idx].y = 0.0;
    euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }
}
sprintf(tag,"M%d RIGID BODY",M_ID);
rigid_body_CS.setup(origin, euler_angles, N, tag);

#line 421 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = (M_ID==1) ? depth: -depth;
euler_angles[N-1].x = euler_angles[N-1].y = euler_angles[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 427 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = (M_ID==1) ? depth: -depth;
  euler_angles[idx].x = 0.0;
  euler_angles[idx].y = 0.0;
  euler_angles[idx].z = 0.0;
}
sprintf(tag,"M%d APERTURE",M_ID);
aperture_CS.setup(origin, euler_angles, N, tag);

#line 463 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = origin[N-1].z = 0.0;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
conic_origin[N-1].x = conic_origin[N-1].y = conic_origin[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 471 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  o = PI*(3-2*k)/6.0;
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = 0.0;
   if (M_ID==1) {
     euler_angles[idx].x = beta;
     euler_angles[idx].y = 0.0;
     euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }     
  conic_origin[idx].x = -D_c*cos(o);
  conic_origin[idx].y = -D_c*sin(o);
  conic_origin[idx].z = -zo;
}
//euler_angles[0].x = 30.0*PI/180.0/3600.0;
sprintf(tag,"M%d CONIC",M_ID);
conic_CS.setup(origin, euler_angles, N, tag);
HANDLE_ERROR( cudaMalloc((void**)&d__conic_origin, sizeof(vector)*N ) );
HANDLE_ERROR( cudaMemcpy( d__conic_origin, conic_origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );

#line 503 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
for (k=0; k<N; k++)
{
  origin[k].x = 0.0;
  origin[k].y = 0.0;
  origin[k].z = 0.0;
  euler_angles[k].x = 0.0;
  euler_angles[k].y = 0.0;
  euler_angles[k].z = 0.0;
}
sprintf(tag,"M%d MOTION",M_ID);
motion_CS.setup(origin, euler_angles, N, tag);

#line 1894 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sprintf(tag,"GLOBAL TIP-TILT");
vector TT_CS_origin, TT_CS_euler_angles;
TT_CS_origin.x = TT_CS_origin.y = TT_CS_origin.z = 0.0;
TT_CS_euler_angles.x = TT_CS_euler_angles.y = TT_CS_euler_angles.z = 0.0;
TT_CS.setup(&TT_CS_origin, &TT_CS_euler_angles, 1, tag);
cublasCreate(&handle);
HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) );
HANDLE_ERROR( cudaMalloc((void**)&d__D, sizeof(double)*9 ) );
#line 1967 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 1979 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::setup(char *_filename_, int _N_, int _n_mode_)
{
  BS.setup(_filename_, _N_, _n_mode_);
  ZS = NULL;
  
#line 1988 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 M_ID    = 2;
 D_assembly = 3.168;
 D_full  = 1.0425*1.005; // slightly larger to always encompass M1; needed for the Zernike
 D_clear = 1.0415;
 ri      = 0.0;
 beta    = -14.777498*PI/180.0;
 L       = 1.08774;
 N       = 7;
 conic_c = -1.0/4.1639009;
 conic_k = 1-0.71692784;
 height  = 20.26247614;
 idx_offset = 3;
#line 1984 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 310 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte;
float segment_reflectivity[7] = {1.,1.,1.,1.,1.,1.,1.};
n_byte = 7*sizeof(float);
HANDLE_ERROR( cudaMalloc((void**)&d__segment_reflectivity, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, segment_reflectivity,
                          n_byte, cudaMemcpyHostToDevice ) );

rtd *h__conic_c;
n_byte = 7*sizeof(rtd);
h__conic_c = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_c[k] = conic_c;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_c, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_c, h__conic_c,
              n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_c);

rtd *h__conic_k;
n_byte = 7*sizeof(rtd);
h__conic_k = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_k[k] = conic_k;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_k, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_k, h__conic_k,
                          n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_k);

area_fraction = 1.0;
vector __v0(0.0,0.0,0.0),
       __v(D_full*0.5,0.0,0.0),
       origin[12],
       euler_angles[12];
rtd D_c, o, zo;
 int k, idx;
char tag[16];

depth = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);
D_c = L;
__v.x = L;
zo = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);

HANDLE_ERROR( cudaMalloc((void**)&d__valid_segments, sizeof(char)*7) );
HANDLE_ERROR( cudaMemset( d__valid_segments, 1, sizeof(char)*7) );

#line 385 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = height;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
for (k=0; k<N-1; k++)
{
  o = PI*(3-2*k)/6.0;
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 394 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = D_c*cos(o);
  origin[idx].y = D_c*sin(o);
  origin[idx].z = height + zo;
  if (M_ID==1) {
    euler_angles[idx].x = beta;
    euler_angles[idx].y = 0.0;
    euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }
}
sprintf(tag,"M%d RIGID BODY",M_ID);
rigid_body_CS.setup(origin, euler_angles, N, tag);

#line 421 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = (M_ID==1) ? depth: -depth;
euler_angles[N-1].x = euler_angles[N-1].y = euler_angles[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 427 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = (M_ID==1) ? depth: -depth;
  euler_angles[idx].x = 0.0;
  euler_angles[idx].y = 0.0;
  euler_angles[idx].z = 0.0;
}
sprintf(tag,"M%d APERTURE",M_ID);
aperture_CS.setup(origin, euler_angles, N, tag);

#line 463 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = origin[N-1].z = 0.0;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
conic_origin[N-1].x = conic_origin[N-1].y = conic_origin[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 471 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  o = PI*(3-2*k)/6.0;
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = 0.0;
   if (M_ID==1) {
     euler_angles[idx].x = beta;
     euler_angles[idx].y = 0.0;
     euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }     
  conic_origin[idx].x = -D_c*cos(o);
  conic_origin[idx].y = -D_c*sin(o);
  conic_origin[idx].z = -zo;
}
//euler_angles[0].x = 30.0*PI/180.0/3600.0;
sprintf(tag,"M%d CONIC",M_ID);
conic_CS.setup(origin, euler_angles, N, tag);
HANDLE_ERROR( cudaMalloc((void**)&d__conic_origin, sizeof(vector)*N ) );
HANDLE_ERROR( cudaMemcpy( d__conic_origin, conic_origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );

#line 503 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
for (k=0; k<N; k++)
{
  origin[k].x = 0.0;
  origin[k].y = 0.0;
  origin[k].z = 0.0;
  euler_angles[k].x = 0.0;
  euler_angles[k].y = 0.0;
  euler_angles[k].z = 0.0;
}
sprintf(tag,"M%d MOTION",M_ID);
motion_CS.setup(origin, euler_angles, N, tag);

#line 1894 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sprintf(tag,"GLOBAL TIP-TILT");
vector TT_CS_origin, TT_CS_euler_angles;
TT_CS_origin.x = TT_CS_origin.y = TT_CS_origin.z = 0.0;
TT_CS_euler_angles.x = TT_CS_euler_angles.y = TT_CS_euler_angles.z = 0.0;
TT_CS.setup(&TT_CS_origin, &TT_CS_euler_angles, 1, tag);
cublasCreate(&handle);
HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) );
HANDLE_ERROR( cudaMalloc((void**)&d__D, sizeof(double)*9 ) );
#line 1985 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2003 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::setup(zernikeS *_ZS_)
{
  ZS = _ZS_;
  //BS = NULL;
  M_ID    = 2;
  
#line 1988 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 M_ID    = 2;
 D_assembly = 3.168;
 D_full  = 1.0425*1.005; // slightly larger to always encompass M1; needed for the Zernike
 D_clear = 1.0415;
 ri      = 0.0;
 beta    = -14.777498*PI/180.0;
 L       = 1.08774;
 N       = 7;
 conic_c = -1.0/4.1639009;
 conic_k = 1-0.71692784;
 height  = 20.26247614;
 idx_offset = 3;
#line 2009 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  
#line 310 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte;
float segment_reflectivity[7] = {1.,1.,1.,1.,1.,1.,1.};
n_byte = 7*sizeof(float);
HANDLE_ERROR( cudaMalloc((void**)&d__segment_reflectivity, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, segment_reflectivity,
                          n_byte, cudaMemcpyHostToDevice ) );

rtd *h__conic_c;
n_byte = 7*sizeof(rtd);
h__conic_c = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_c[k] = conic_c;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_c, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_c, h__conic_c,
              n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_c);

rtd *h__conic_k;
n_byte = 7*sizeof(rtd);
h__conic_k = (rtd *)malloc( n_byte );
for (int k=0; k<7; k++)
  h__conic_k[k] = conic_k;
HANDLE_ERROR( cudaMalloc((void**)&d__conic_k, n_byte ) );
HANDLE_ERROR( cudaMemcpy( d__conic_k, h__conic_k,
                          n_byte, cudaMemcpyHostToDevice ) );
free(h__conic_k);

area_fraction = 1.0;
vector __v0(0.0,0.0,0.0),
       __v(D_full*0.5,0.0,0.0),
       origin[12],
       euler_angles[12];
rtd D_c, o, zo;
 int k, idx;
char tag[16];

depth = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);
D_c = L;
__v.x = L;
zo = conic_equation(&__v,&__v0,conic_k,conic_c, 0, &D_c);

HANDLE_ERROR( cudaMalloc((void**)&d__valid_segments, sizeof(char)*7) );
HANDLE_ERROR( cudaMemset( d__valid_segments, 1, sizeof(char)*7) );

#line 385 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = height;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
for (k=0; k<N-1; k++)
{
  o = PI*(3-2*k)/6.0;
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 394 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = D_c*cos(o);
  origin[idx].y = D_c*sin(o);
  origin[idx].z = height + zo;
  if (M_ID==1) {
    euler_angles[idx].x = beta;
    euler_angles[idx].y = 0.0;
    euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }
}
sprintf(tag,"M%d RIGID BODY",M_ID);
rigid_body_CS.setup(origin, euler_angles, N, tag);

#line 421 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = 0.0;
origin[N-1].z = (M_ID==1) ? depth: -depth;
euler_angles[N-1].x = euler_angles[N-1].y = euler_angles[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 427 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = (M_ID==1) ? depth: -depth;
  euler_angles[idx].x = 0.0;
  euler_angles[idx].y = 0.0;
  euler_angles[idx].z = 0.0;
}
sprintf(tag,"M%d APERTURE",M_ID);
aperture_CS.setup(origin, euler_angles, N, tag);

#line 463 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
origin[N-1].x = origin[N-1].y = origin[N-1].z = 0.0;
euler_angles[N-1].x = 0.0;
euler_angles[N-1].y = (M_ID==1) ? 0.0 : PI;
euler_angles[N-1].z = (idx_offset/2)*PI;
conic_origin[N-1].x = conic_origin[N-1].y = conic_origin[N-1].z = 0.0;
for (k=0; k<N-1; k++)
{
  
#line 412 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
idx = (k + idx_offset)%6;
#line 471 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  o = PI*(3-2*k)/6.0;
  origin[idx].x = 0.0;
  origin[idx].y = 0.0;
  origin[idx].z = 0.0;
   if (M_ID==1) {
     euler_angles[idx].x = beta;
     euler_angles[idx].y = 0.0;
     euler_angles[idx].z = -PI*k/3.0;
  }
  if (M_ID==2) {
    euler_angles[idx].x = -beta;
    euler_angles[idx].y = PI;
    euler_angles[idx].z = -PI*k/3.0;
  }     
  conic_origin[idx].x = -D_c*cos(o);
  conic_origin[idx].y = -D_c*sin(o);
  conic_origin[idx].z = -zo;
}
//euler_angles[0].x = 30.0*PI/180.0/3600.0;
sprintf(tag,"M%d CONIC",M_ID);
conic_CS.setup(origin, euler_angles, N, tag);
HANDLE_ERROR( cudaMalloc((void**)&d__conic_origin, sizeof(vector)*N ) );
HANDLE_ERROR( cudaMemcpy( d__conic_origin, conic_origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );

#line 503 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
for (k=0; k<N; k++)
{
  origin[k].x = 0.0;
  origin[k].y = 0.0;
  origin[k].z = 0.0;
  euler_angles[k].x = 0.0;
  euler_angles[k].y = 0.0;
  euler_angles[k].z = 0.0;
}
sprintf(tag,"M%d MOTION",M_ID);
motion_CS.setup(origin, euler_angles, N, tag);

#line 1894 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
sprintf(tag,"GLOBAL TIP-TILT");
vector TT_CS_origin, TT_CS_euler_angles;
TT_CS_origin.x = TT_CS_origin.y = TT_CS_origin.z = 0.0;
TT_CS_euler_angles.x = TT_CS_euler_angles.y = TT_CS_euler_angles.z = 0.0;
TT_CS.setup(&TT_CS_origin, &TT_CS_euler_angles, 1, tag);
cublasCreate(&handle);
HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) );
HANDLE_ERROR( cudaMalloc((void**)&d__D, sizeof(double)*9 ) );
#line 2010 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2015 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::cleanup(void)
{
  INFO("@(CEO)>gmt_m2: freeing memory!\n");
  
#line 530 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
INFO(" |-");
BS.cleanup();
INFO(" |-");
#line 1904 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
TT_CS.cleanup();
cublasDestroy(handle);
HANDLE_ERROR( cudaFree( d__C ) );
HANDLE_ERROR( cudaFree( d__D ) );
#line 534 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
INFO(" |-");
aperture_CS.cleanup();
INFO(" |-");
conic_CS.cleanup();
INFO(" |-");
rigid_body_CS.cleanup();
INFO(" |-");
motion_CS.cleanup();
HANDLE_ERROR( cudaFree( d__conic_origin ) );
HANDLE_ERROR( cudaFree( d__valid_segments ) );
HANDLE_ERROR( cudaFree( d__conic_c ) );
HANDLE_ERROR( cudaFree( d__conic_k ) );
HANDLE_ERROR( cudaFree( d__segment_reflectivity ) );
#line 2019 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2022 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::update_conic_c(rtd *_conic_c_)
{
  HANDLE_ERROR( cudaMemcpy( d__conic_c, _conic_c_,
              sizeof(rtd)*7, cudaMemcpyHostToDevice ) );
}
#line 2029 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  void gmt_m2::update_conic_k(rtd *_conic_k_)
{
  HANDLE_ERROR( cudaMemcpy( d__conic_k, _conic_k_,
                            sizeof(rtd)*7, cudaMemcpyHostToDevice ) );
}
#line 2041 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::blocking(bundle *rays)
{
  
#line 646 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 rtd R2, Rri2;
 R2 = D_clear*D_clear*0.25;
 Rri2 = R2*ri*ri;
 V = &(rays->V);
 // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2);
 // printf("N_BUNDLE=%d\n",rays->N_BUNDLE);
 int nel = rays->N_RAY*rays->N_BUNDLE;
 dim3 blockDim(N_THREAD2);
 dim3 gridDim(nel/N_THREAD2+1);
 fill_ones_char <<< gridDim,blockDim >>> (V->m,nel);
 blockDim = dim3(N_THREAD,N_THREAD);
 gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE);
 m1_blocking_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
					     Rri2, R2,
					     aperture_CS.d__R,
					     aperture_CS.d__origin,
					     rigid_body_CS.d__R,
					     rigid_body_CS.d__origin,
					     motion_CS.d__R,
					     motion_CS.d__origin);
 intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY);
#line 2044 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2052 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::trace(bundle *rays)
{
  dim3 blockDim, gridDim;
  rtd R2, Rri2;
  R2 = D_clear*D_clear*0.25;
  Rri2 = R2*ri*ri;
  
#line 729 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 /* if (V!=NULL) { */
 /*    blockDim = dim3(N_THREAD,N_THREAD); */
 /*    gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); */
 /*    intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY); */
 /* } */
 V = &(rays->V);
 // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2);
 blockDim = dim3(1,1);
 gridDim  = dim3(1,1,rays->N_BUNDLE);
 m1_trace_chief_kernel <<< gridDim , blockDim >>> (V->m, rays->d__chief_ray, 1,
						   conic_CS.d__R,
						   conic_CS.d__origin,
						   conic_k, conic_c,
						   d__conic_origin,
						   rigid_body_CS.d__R,
						   rigid_body_CS.d__origin,
						   motion_CS.d__R,
						   motion_CS.d__origin);
 blockDim = dim3(N_THREAD,N_THREAD);
 gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE);
 HANDLE_ERROR( cudaMemset(V->m, 0, sizeof(char)*rays->N_RAY*rays->N_BUNDLE ) );
 HANDLE_ERROR( cudaMemset(rays->d__piston_mask, 0,
			  sizeof(int)*rays->N_RAY*rays->N_BUNDLE ) );
/*if (ZS!=NULL)
  m1_trace_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
                                              Rri2, R2,
                                              aperture_CS.d__R,
                                              aperture_CS.d__origin,
                                              conic_CS.d__R,
                                              conic_CS.d__origin,
                                              d__conic_k, d__conic_c,
                                              d__conic_origin,
                                              rigid_body_CS.d__R,
                                              rigid_body_CS.d__origin,
                                              motion_CS.d__R,
                                              motion_CS.d__origin,
                                              ZS->max_n, ZS->n_mode,
                                              ZS->d__a, ZS->d__cx, ZS->d__cy,
                                              0.5*D_full,
                                              rays->d__chief_ray,
                                              rays->d__piston_mask,
                                              d__valid_segments,
                                              d__segment_reflectivity);

					      if (BS!=NULL)*/
  m1_bm_trace_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
                                                 Rri2, R2,
                                                 aperture_CS.d__R,
                                                 aperture_CS.d__origin,
                                                 conic_CS.d__R,
                                                 conic_CS.d__origin,
                                                 d__conic_k, d__conic_c,
                                                 d__conic_origin,
                                                 rigid_body_CS.d__R,
                                                 rigid_body_CS.d__origin,
                                                 motion_CS.d__R,
                                                 motion_CS.d__origin,
                                                 BS.n_mode, BS.d__b,
                                                 BS.d__BMS,
                                                 BS.BM_N_SAMPLE, BS.BM_radius,
                                                 rays->d__chief_ray,
                                                 rays->d__piston_mask,
                                                 d__valid_segments,
                                                 d__segment_reflectivity);
intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY);
float previous_nnz = rays->V.nnz;
//fprintf(stdout,"M%d previous nnz rays: %f\n",M_ID,previous_nnz);
rays->V.set_filter_quiet();
//fprintf(stdout,"M%d current nnz rays: %f\n",M_ID,rays->V.nnz);
if (previous_nnz<rays->N_RAY_TOTAL)
  rays->V.area *= rays->V.nnz/previous_nnz;
//fprintf(stdout,"M%d current area: %f\n",M_ID,rays->V.area);
#line 2059 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
void gmt_m2::traceall(bundle *rays)
{
  rays->V.area = area0*rays->N_BUNDLE*area_fraction;
  dim3 blockDim, gridDim;
  rtd R2, Rri2;
  R2 = D_clear*D_clear*0.25;
  Rri2 = R2*ri*ri*0;
  
#line 1485 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 /* if (V!=NULL) { */
 /*    blockDim = dim3(N_THREAD,N_THREAD); */
 /*    gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE); */
 /*    intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY); */
 /* } */
 V = &(rays->V);
 // fprintf(stdout,"R2=%5.2f - Rri2=%5.2f\n",R2,Rri2);
 blockDim = dim3(1,1);
 gridDim  = dim3(1,1,rays->N_BUNDLE);
 m1_trace_chief_kernel <<< gridDim , blockDim >>> (V->m, rays->d__chief_ray, 1,
						   conic_CS.d__R,
						   conic_CS.d__origin,
						   conic_k, conic_c,
						   d__conic_origin,
						   rigid_body_CS.d__R,
						   rigid_body_CS.d__origin,
						   motion_CS.d__R,
						   motion_CS.d__origin);
 blockDim = dim3(N_THREAD,N_THREAD);
 gridDim  = dim3(rays->N_RAY/N_THREAD2+1,1,rays->N_BUNDLE);
/*HANDLE_ERROR( cudaMemset(V->m, 0, sizeof(char)*rays->N_RAY*rays->N_BUNDLE ) );
 HANDLE_ERROR( cudaMemset(rays->d__piston_mask, 0,
			  sizeof(int)*rays->N_RAY*rays->N_BUNDLE ) );
*/
if (ZS!=NULL)
  m1_trace_all_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
                                              Rri2, R2,
                                              aperture_CS.d__R,
                                              aperture_CS.d__origin,
                                              conic_CS.d__R,
                                              conic_CS.d__origin,
                                              d__conic_k, d__conic_c,
                                              d__conic_origin,
                                              rigid_body_CS.d__R,
                                              rigid_body_CS.d__origin,
                                              motion_CS.d__R,
                                              motion_CS.d__origin,
                                              ZS->max_n, ZS->n_mode,
                                              ZS->d__a, ZS->d__cx, ZS->d__cy,
                                              0.5*D_full,
                                              rays->d__chief_ray,
                                              rays->d__piston_mask,
                                              d__valid_segments,
                                              d__segment_reflectivity);
/*
if (BS!=NULL)
  m1_bm_trace_kernel <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY,
                                                 Rri2, R2,
                                                 aperture_CS.d__R,
                                                 aperture_CS.d__origin,
                                                 conic_CS.d__R,
                                                 conic_CS.d__origin,
                                                 d__conic_k, d__conic_c,
                                                 d__conic_origin,
                                                 rigid_body_CS.d__R,
                                                 rigid_body_CS.d__origin,
                                                 motion_CS.d__R,
                                                 motion_CS.d__origin,
                                                 BS.n_mode, BS.d__b,
                                                 BS.d__BMS,
                                                 BS.BM_N_SAMPLE, BS.BM_radius,
                                                 rays->d__chief_ray,
                                                 rays->d__piston_mask,
                                                 d__valid_segments,
                                                 d__segment_reflectivity);

intersection <<< gridDim , blockDim >>> (V->m, rays->d__ray, rays->N_RAY);
float previous_nnz = rays->V.nnz;
//fprintf(stdout,"M%d previous nnz rays: %f\n",M_ID,previous_nnz);
rays->V.set_filter_quiet();
//fprintf(stdout,"M%d current nnz rays: %f\n",M_ID,rays->V.nnz);
if (previous_nnz<rays->N_RAY_TOTAL)
  rays->V.area *= rays->V.nnz/previous_nnz;
//fprintf(stdout,"M%d current area: %f\n",M_ID,rays->V.area);
*/
#line 2068 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
      }
#line 2073 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::update(vector _origin_, vector _euler_angles_,int idx)
{
    
#line 575 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  //fprintf(stdout,"\n\x1B[31m");
motion_CS.update( _origin_, _euler_angles_, --idx);
//motion_CS.info();
#line 2076 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2081 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::reset(void)
{
    
#line 588 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
vector zero;
zero.x = 0.0;
zero.y = 0.0;
zero.z = 0.0;
//fprintf(stdout,"\n\x1B[31m");
for (int k=0; k<N; k++)
   motion_CS.update( zero, zero, k);
//motion_CS.info();
#line 2084 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2092 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::remove(int *seg_ID, int N_ID)
{
  
#line 1631 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int k;
char valid[7];
area_fraction = (7-N_ID)/7.0;
memset(valid, 1, sizeof(char)*7);
for (k=0;k<N_ID;k++)
  valid[seg_ID[k]-1] = 0;
HANDLE_ERROR( cudaMemcpy( d__valid_segments, valid,
                          sizeof(char)*7, cudaMemcpyHostToDevice ) );
#line 2095 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2100 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::keep(int *seg_ID, int N_ID)
{
  
#line 1650 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int k;
char valid[7];
area_fraction = N_ID/7.0;
memset(valid, 0, sizeof(char)*7);
for (k=0;k<N_ID;k++)
  valid[seg_ID[k]-1] = 1;
HANDLE_ERROR( cudaMemcpy( d__valid_segments, valid,
                          sizeof(char)*7, cudaMemcpyHostToDevice ) );
#line 2103 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2112 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::track(float *d__x, float *d__y, float *d__z, int N, int idx)
{
  
#line 1672 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
dim3 blockDim(N_THREAD,1);
dim3 gridDim(N/N_THREAD+1,1);
track_kernel <<< gridDim, blockDim >>> (d__x, d__y, d__z, N,
					idx,
					rigid_body_CS.d__R,
					rigid_body_CS.d__origin,
					motion_CS.d__R,
					motion_CS.d__origin);
#line 2115 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2123 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::locate(float *d__x, float *d__y, float *d__z, int N, int idx)
{
  
#line 1720 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
dim3 blockDim(N_THREAD,1);
dim3 gridDim(N/N_THREAD+1,1);
locate_kernel <<< gridDim, blockDim >>> (d__x, d__y, d__z, N,
					idx,
					rigid_body_CS.d__R,
					rigid_body_CS.d__origin,
					motion_CS.d__R,
					motion_CS.d__origin);
#line 2126 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2134 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::global_tiptilt(float tip, float tilt)
{
  
#line 1910 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
double alpha, beta;
alpha = 1.0;
beta  = 0.0;
int k, k9;
vector origin, euler_angles, o;
origin.x = origin.y = origin.z = 0.0;
euler_angles.x = tip;
euler_angles.y = tilt;
euler_angles.z = 0.0;
TT_CS.update(origin,euler_angles,0);
for (k=0;k<N;k++)
  {
    k9 = k*9;
    
#line 1794 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_T,
	     3,3,3,
	     &alpha,
	     TT_CS.d__R, 3,
	     rigid_body_CS.d__R+k9, 3,
	     &beta,
	     d__C, 3);
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_N,
	     3,3,3,
	     &alpha,
	     rigid_body_CS.d__R+k9, 3,
	     d__C, 3,
	     &beta,
	     d__D, 3);
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_N,
	     3,3,3,
	     &alpha,
	     motion_CS.d__R+k9, 3,
	     d__D, 3,
	     &beta,
	     d__C, 3);

 HANDLE_ERROR( cudaMemcpy( motion_CS.d__R+k9, d__C,
			   sizeof(double)*9, cudaMemcpyDeviceToDevice ) );
 HANDLE_ERROR( cudaMemcpy( motion_CS.R+k9, motion_CS.d__R+k9,
			   sizeof(double)*9, cudaMemcpyDeviceToHost ) );
 /*
 int i, j, idx = 0;
for (i=0;i<3;i++)
   {
     fprintf(stdout,"||  ");
     for (j=0;j<3;j++)
       {
	 idx = j + i*3 + k9;
	 fprintf(stdout,"%+.4e  ",motion_CS.R[idx]);
       }
     fprintf(stdout,"||\n");
   }
 */
#line 1845 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
motion_CS.euler_angles[k].x = atan2(motion_CS.R[7+k9],motion_CS.R[8+k9]);
motion_CS.euler_angles[k].y = asin(-1.0*motion_CS.R[6+k9]);
motion_CS.euler_angles[k].z = atan2(motion_CS.R[3+k9],motion_CS.R[k9]);
#line 1856 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
vector o;
forward_transform(&o, rigid_body_CS.origin + k,
		  TT_CS.R, TT_CS.origin);
o.x = rigid_body_CS.origin[k].x - TT_CS.origin->x - o.x;
o.y = rigid_body_CS.origin[k].y - TT_CS.origin->y - o.y;
o.z = rigid_body_CS.origin[k].z - TT_CS.origin->z - o.z;
#line 1865 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
forward_transform_centered(&o, &o, rigid_body_CS.R + k9);
motion_CS.origin[k].x += o.x; 
motion_CS.origin[k].y += o.y;
motion_CS.origin[k].z += o.z;
HANDLE_ERROR( cudaMemcpy( motion_CS.d__origin+k, motion_CS.origin+k,
			  sizeof(vector), cudaMemcpyHostToDevice ) );

/*
 fprintf(stdout,"#%d >> MOTION CS origins [micron]  : %+.3e, %+.3e, %+.3e\n",k,
	 motion_CS.origin[k].x*1e6,
	 motion_CS.origin[k].y*1e6,
	 motion_CS.origin[k].z*1e6);

 float r2d;
 r2d = 1000*3600*180.0/PI;
 fprintf(stdout,"#%d >> MOTION CS Euler angles [mas]: %+.3e, %+.3e, %+.3e\n",k,
	 r2d*motion_CS.euler_angles[k].x,
	 r2d*motion_CS.euler_angles[k].y,
	 r2d*motion_CS.euler_angles[k].z);
*/
#line 1924 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
  }
#line 2137 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 2145 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::pointing_neutral(float tip, float tilt)
{
  double alpha, beta, *d__C;
  vector origin, euler_angles;
  cublasHandle_t handle;
  coordinate_system TT_CS;

  origin.x = origin.y = 0.0;
  origin.z = height - 4390.312E-03;
  euler_angles.x = tip;
  euler_angles.y = tilt;
  euler_angles.z = 0.0;

  TT_CS.setup(origin, euler_angles);
  cublasCreate(&handle);

  alpha = 1;
  beta  = 0;
  HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) );
  int k, k9;
  for (k=0;k<N;k++)
    {
      k9 = k*9;
      
#line 1794 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_T,
	     3,3,3,
	     &alpha,
	     TT_CS.d__R, 3,
	     rigid_body_CS.d__R+k9, 3,
	     &beta,
	     d__C, 3);
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_N,
	     3,3,3,
	     &alpha,
	     rigid_body_CS.d__R+k9, 3,
	     d__C, 3,
	     &beta,
	     d__D, 3);
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_N,
	     3,3,3,
	     &alpha,
	     motion_CS.d__R+k9, 3,
	     d__D, 3,
	     &beta,
	     d__C, 3);

 HANDLE_ERROR( cudaMemcpy( motion_CS.d__R+k9, d__C,
			   sizeof(double)*9, cudaMemcpyDeviceToDevice ) );
 HANDLE_ERROR( cudaMemcpy( motion_CS.R+k9, motion_CS.d__R+k9,
			   sizeof(double)*9, cudaMemcpyDeviceToHost ) );
 /*
 int i, j, idx = 0;
for (i=0;i<3;i++)
   {
     fprintf(stdout,"||  ");
     for (j=0;j<3;j++)
       {
	 idx = j + i*3 + k9;
	 fprintf(stdout,"%+.4e  ",motion_CS.R[idx]);
       }
     fprintf(stdout,"||\n");
   }
 */
#line 1845 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
motion_CS.euler_angles[k].x = atan2(motion_CS.R[7+k9],motion_CS.R[8+k9]);
motion_CS.euler_angles[k].y = asin(-1.0*motion_CS.R[6+k9]);
motion_CS.euler_angles[k].z = atan2(motion_CS.R[3+k9],motion_CS.R[k9]);
#line 1856 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
vector o;
forward_transform(&o, rigid_body_CS.origin + k,
		  TT_CS.R, TT_CS.origin);
o.x = rigid_body_CS.origin[k].x - TT_CS.origin->x - o.x;
o.y = rigid_body_CS.origin[k].y - TT_CS.origin->y - o.y;
o.z = rigid_body_CS.origin[k].z - TT_CS.origin->z - o.z;
#line 1865 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
forward_transform_centered(&o, &o, rigid_body_CS.R + k9);
motion_CS.origin[k].x += o.x; 
motion_CS.origin[k].y += o.y;
motion_CS.origin[k].z += o.z;
HANDLE_ERROR( cudaMemcpy( motion_CS.d__origin+k, motion_CS.origin+k,
			  sizeof(vector), cudaMemcpyHostToDevice ) );

/*
 fprintf(stdout,"#%d >> MOTION CS origins [micron]  : %+.3e, %+.3e, %+.3e\n",k,
	 motion_CS.origin[k].x*1e6,
	 motion_CS.origin[k].y*1e6,
	 motion_CS.origin[k].z*1e6);

 float r2d;
 r2d = 1000*3600*180.0/PI;
 fprintf(stdout,"#%d >> MOTION CS Euler angles [mas]: %+.3e, %+.3e, %+.3e\n",k,
	 r2d*motion_CS.euler_angles[k].x,
	 r2d*motion_CS.euler_angles[k].y,
	 r2d*motion_CS.euler_angles[k].z);
*/
#line 2169 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    }

  cublasDestroy(handle);
  TT_CS.cleanup();
  HANDLE_ERROR( cudaFree( d__C ) );
}
#line 2181 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::coma_neutral(float tip, float tilt)
{
  double alpha, beta, *d__C;
  vector origin, euler_angles;
  cublasHandle_t handle;
  coordinate_system TT_CS;

  origin.x = origin.y = 0.0;
  origin.z = height - 2246.410E-03;
  euler_angles.x = tip;
  euler_angles.y = tilt;
  euler_angles.z = 0.0;

  TT_CS.setup(origin, euler_angles);
  cublasCreate(&handle);

  alpha = 1;
  beta  = 0;
  HANDLE_ERROR( cudaMalloc((void**)&d__C, sizeof(double)*9 ) );
  int k, k9;
  for (k=0;k<N;k++)
    {
      k9 = k*9;
      
#line 1794 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_T,
	     3,3,3,
	     &alpha,
	     TT_CS.d__R, 3,
	     rigid_body_CS.d__R+k9, 3,
	     &beta,
	     d__C, 3);
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_N,
	     3,3,3,
	     &alpha,
	     rigid_body_CS.d__R+k9, 3,
	     d__C, 3,
	     &beta,
	     d__D, 3);
 cublasDgemm(handle,
	     CUBLAS_OP_N, CUBLAS_OP_N,
	     3,3,3,
	     &alpha,
	     motion_CS.d__R+k9, 3,
	     d__D, 3,
	     &beta,
	     d__C, 3);

 HANDLE_ERROR( cudaMemcpy( motion_CS.d__R+k9, d__C,
			   sizeof(double)*9, cudaMemcpyDeviceToDevice ) );
 HANDLE_ERROR( cudaMemcpy( motion_CS.R+k9, motion_CS.d__R+k9,
			   sizeof(double)*9, cudaMemcpyDeviceToHost ) );
 /*
 int i, j, idx = 0;
for (i=0;i<3;i++)
   {
     fprintf(stdout,"||  ");
     for (j=0;j<3;j++)
       {
	 idx = j + i*3 + k9;
	 fprintf(stdout,"%+.4e  ",motion_CS.R[idx]);
       }
     fprintf(stdout,"||\n");
   }
 */
#line 1845 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
motion_CS.euler_angles[k].x = atan2(motion_CS.R[7+k9],motion_CS.R[8+k9]);
motion_CS.euler_angles[k].y = asin(-1.0*motion_CS.R[6+k9]);
motion_CS.euler_angles[k].z = atan2(motion_CS.R[3+k9],motion_CS.R[k9]);
#line 1856 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
vector o;
forward_transform(&o, rigid_body_CS.origin + k,
		  TT_CS.R, TT_CS.origin);
o.x = rigid_body_CS.origin[k].x - TT_CS.origin->x - o.x;
o.y = rigid_body_CS.origin[k].y - TT_CS.origin->y - o.y;
o.z = rigid_body_CS.origin[k].z - TT_CS.origin->z - o.z;
#line 1865 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
forward_transform_centered(&o, &o, rigid_body_CS.R + k9);
motion_CS.origin[k].x += o.x; 
motion_CS.origin[k].y += o.y;
motion_CS.origin[k].z += o.z;
HANDLE_ERROR( cudaMemcpy( motion_CS.d__origin+k, motion_CS.origin+k,
			  sizeof(vector), cudaMemcpyHostToDevice ) );

/*
 fprintf(stdout,"#%d >> MOTION CS origins [micron]  : %+.3e, %+.3e, %+.3e\n",k,
	 motion_CS.origin[k].x*1e6,
	 motion_CS.origin[k].y*1e6,
	 motion_CS.origin[k].z*1e6);

 float r2d;
 r2d = 1000*3600*180.0/PI;
 fprintf(stdout,"#%d >> MOTION CS Euler angles [mas]: %+.3e, %+.3e, %+.3e\n",k,
	 r2d*motion_CS.euler_angles[k].x,
	 r2d*motion_CS.euler_angles[k].y,
	 r2d*motion_CS.euler_angles[k].z);
*/
#line 2205 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
    }

  cublasDestroy(handle);
  TT_CS.cleanup();
  HANDLE_ERROR( cudaFree( d__C ) );
}
#line 2218 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void gmt_m2::set_reflectivity(float *reflectivity)
{
  
#line 1938 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
int n_byte;
n_byte = 7*sizeof(float);
HANDLE_ERROR( cudaMemcpy( d__segment_reflectivity, reflectivity,
                          n_byte, cudaMemcpyHostToDevice ) );
#line 2221 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}

#line 3308 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void stereoscopic_edge_sensors::setup(gmt_m1 *_mirror_)
{
  
#line 3327 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
N = 48;
N_DATA = 24;
n_byte      = sizeof(vector)*N;
n_byte_DATA = sizeof(vector)*N_DATA;
#line 3335 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
alpha = 19.5;
#line 3340 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
mirror = _mirror_;
#line 3346 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
v0 = (vector *)malloc(n_byte);
v  = (vector *)malloc(n_byte);
HANDLE_ERROR( cudaMalloc((void**)&d__v0, n_byte ) );
HANDLE_ERROR( cudaMalloc((void**)&d__v,  n_byte ) );
#line 3361 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
dv0  = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__dv0, n_byte_DATA ) );
dv  = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__dv, n_byte_DATA ) );
#line 3373 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
INFO("@(CEO)>gmt_m1: Initializing edge sensors coordinates!\n");
dim3 blockDim(6,1);
dim3 gridDim(1,1);
edge_sensors_read <<< gridDim, blockDim >>> (d__v0, mirror->N-1, mirror->D_full, alpha,
					     mirror->rigid_body_CS.d__R,
                                             mirror->rigid_body_CS.d__origin,
                                             mirror->motion_CS.d__R,
                                             mirror->motion_CS.d__origin);
HANDLE_ERROR( cudaMemcpy( v0, d__v0,
                          n_byte, cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaMemcpy( d__v, d__v0,
                          n_byte, cudaMemcpyDeviceToDevice ) );
HANDLE_ERROR( cudaMemcpy( v, d__v,
                          n_byte, cudaMemcpyDeviceToHost ) );
#line 3389 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
HANDLE_ERROR( cudaMemset(d__dv0, 0, n_byte_DATA ) );
edge_sensors_data_kernel <<< gridDim, blockDim >>> (d__dv,
                                                   d__v,
                                                   d__dv0,
                                                   mirror->N-1);
HANDLE_ERROR( cudaMemcpy( d__dv0, d__dv,
                          n_byte_DATA, cudaMemcpyDeviceToDevice ) );
HANDLE_ERROR( cudaMemcpy( dv0, d__dv0,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaMemset(d__dv, 0, n_byte_DATA ) );
HANDLE_ERROR( cudaMemcpy( dv, d__dv,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
#line 3311 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
void stereoscopic_edge_sensors::cleanup(void)
{
  INFO("@(CEO)>stereoscopic_edge_sensors: freeing memory!\n");
  
#line 3351 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
free( v0 );
free( v );
HANDLE_ERROR( cudaFree( d__v0 ) );
HANDLE_ERROR( cudaFree( d__v ) );
#line 3366 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
free( dv0 );
HANDLE_ERROR( cudaFree( d__dv0 ) );
free( dv );
HANDLE_ERROR( cudaFree( d__dv ) );
#line 3316 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
void stereoscopic_edge_sensors::data(void)
{
  
#line 3506 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
dim3 blockDim(6,1);
dim3 gridDim(1,1);
edge_sensors_read <<< gridDim, blockDim >>> (d__v, mirror->N-1, mirror->D_full, alpha,
       					     mirror->rigid_body_CS.d__R,
                                             mirror->rigid_body_CS.d__origin,
                                             mirror->motion_CS.d__R,
                                             mirror->motion_CS.d__origin);
edge_sensors_data_kernel <<< gridDim, blockDim >>> (d__dv,
						     d__v,
						     d__dv0,
						     mirror->N-1);
HANDLE_ERROR( cudaMemcpy( dv, d__dv,
			    sizeof(vector)*N_DATA, cudaMemcpyDeviceToHost ) );
#line 3320 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
#line 3569 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
void lateral_edge_sensors::setup(gmt_m1 *_mirror_)
{
  N_HEIGHT = 1;
  height = 0.0;
  
#line 3602 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
//N = 48;
N_DATA = 24*N_HEIGHT;
#line 3608 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
alpha = 19.5;
#line 3617 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
mirror = _mirror_;
#line 3633 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte_DATA = sizeof(vector)*N_DATA;
k_cam  = (vector *)malloc(n_byte_DATA);
k_laser  = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__k_cam, n_byte_DATA ) );
HANDLE_ERROR( cudaMalloc((void**)&d__k_laser, n_byte_DATA ) );
uP  = (vector *)malloc(n_byte_DATA);
vP  = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__uP, n_byte_DATA ) );
HANDLE_ERROR( cudaMalloc((void**)&d__vP, n_byte_DATA ) );
x  = (rtd*)malloc(sizeof(rtd)*N_DATA);
y  = (rtd*)malloc(sizeof(rtd)*N_DATA);
d  = (rtd*)malloc(sizeof(rtd)*N_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__x, sizeof(rtd)*N_DATA ) );
HANDLE_ERROR( cudaMalloc((void**)&d__y, sizeof(rtd)*N_DATA ) );
HANDLE_ERROR( cudaMalloc((void**)&d__d, sizeof(rtd)*N_DATA ) );
#line 3757 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
B = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__B, n_byte_DATA ) );
B0 = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__B0, n_byte_DATA ) );
A = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__A, n_byte_DATA ) );
A0 = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__A0, n_byte_DATA ) );
#line 3782 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
INFO("@(CEO)>gmt_m1: Initializing edge sensors coordinates!\n");
dim3 blockDim(6,N_HEIGHT);
dim3 gridDim(1,1);
laser_coordinates <<< gridDim, blockDim >>> (d__A, mirror->N-1, mirror->D_full,
					     alpha, height,
					     mirror->rigid_body_CS.d__R,
                                             mirror->rigid_body_CS.d__origin,
                                             mirror->motion_CS.d__R,
                                             mirror->motion_CS.d__origin);
HANDLE_ERROR( cudaMemcpy( A, d__A,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
laser_coordinates_to_GCS <<< gridDim, blockDim >>> (d__A0, d__A, mirror->N-1,
                                                    mirror->rigid_body_CS.d__R,
                                                    mirror->rigid_body_CS.d__origin);
HANDLE_ERROR( cudaMemcpy( A0, d__A0,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
#line 3903 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
camera_coordinates <<< gridDim, blockDim >>> (d__B, d__B0, mirror->N-1,
                                              mirror->D_full, alpha, height,
					      mirror->rigid_body_CS.d__R,
					      mirror->rigid_body_CS.d__origin);
HANDLE_ERROR( cudaMemcpy(B , d__B,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaMemcpy(B0 , d__B0,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
#line 3985 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
edge_sensors_camera_vector_kernel <<< gridDim, blockDim >>> (d__k_cam, d__A, d__B,
                                                           d__uP, d__vP,
                                                           mirror->N-1);
HANDLE_ERROR( cudaMemcpy(k_cam , d__k_cam,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
#line 3574 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
void lateral_edge_sensors::setup(gmt_m1 *_mirror_, rtd _height_)
{
  N_HEIGHT = 2;
  height = _height_;
  
#line 3602 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
//N = 48;
N_DATA = 24*N_HEIGHT;
#line 3608 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
alpha = 19.5;
#line 3617 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
mirror = _mirror_;
#line 3633 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
n_byte_DATA = sizeof(vector)*N_DATA;
k_cam  = (vector *)malloc(n_byte_DATA);
k_laser  = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__k_cam, n_byte_DATA ) );
HANDLE_ERROR( cudaMalloc((void**)&d__k_laser, n_byte_DATA ) );
uP  = (vector *)malloc(n_byte_DATA);
vP  = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__uP, n_byte_DATA ) );
HANDLE_ERROR( cudaMalloc((void**)&d__vP, n_byte_DATA ) );
x  = (rtd*)malloc(sizeof(rtd)*N_DATA);
y  = (rtd*)malloc(sizeof(rtd)*N_DATA);
d  = (rtd*)malloc(sizeof(rtd)*N_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__x, sizeof(rtd)*N_DATA ) );
HANDLE_ERROR( cudaMalloc((void**)&d__y, sizeof(rtd)*N_DATA ) );
HANDLE_ERROR( cudaMalloc((void**)&d__d, sizeof(rtd)*N_DATA ) );
#line 3757 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
B = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__B, n_byte_DATA ) );
B0 = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__B0, n_byte_DATA ) );
A = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__A, n_byte_DATA ) );
A0 = (vector *)malloc(n_byte_DATA);
HANDLE_ERROR( cudaMalloc((void**)&d__A0, n_byte_DATA ) );
#line 3782 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
INFO("@(CEO)>gmt_m1: Initializing edge sensors coordinates!\n");
dim3 blockDim(6,N_HEIGHT);
dim3 gridDim(1,1);
laser_coordinates <<< gridDim, blockDim >>> (d__A, mirror->N-1, mirror->D_full,
					     alpha, height,
					     mirror->rigid_body_CS.d__R,
                                             mirror->rigid_body_CS.d__origin,
                                             mirror->motion_CS.d__R,
                                             mirror->motion_CS.d__origin);
HANDLE_ERROR( cudaMemcpy( A, d__A,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
laser_coordinates_to_GCS <<< gridDim, blockDim >>> (d__A0, d__A, mirror->N-1,
                                                    mirror->rigid_body_CS.d__R,
                                                    mirror->rigid_body_CS.d__origin);
HANDLE_ERROR( cudaMemcpy( A0, d__A0,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
#line 3903 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
camera_coordinates <<< gridDim, blockDim >>> (d__B, d__B0, mirror->N-1,
                                              mirror->D_full, alpha, height,
					      mirror->rigid_body_CS.d__R,
					      mirror->rigid_body_CS.d__origin);
HANDLE_ERROR( cudaMemcpy(B , d__B,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaMemcpy(B0 , d__B0,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
#line 3985 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
edge_sensors_camera_vector_kernel <<< gridDim, blockDim >>> (d__k_cam, d__A, d__B,
                                                           d__uP, d__vP,
                                                           mirror->N-1);
HANDLE_ERROR( cudaMemcpy(k_cam , d__k_cam,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
#line 3580 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
void lateral_edge_sensors::cleanup(void)
{
  INFO("@(CEO)>lateral_edge_sensors: freeing memory!\n");
  
#line 3650 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
free( k_cam );
free( k_laser );
HANDLE_ERROR( cudaFree( d__k_cam ) );
HANDLE_ERROR( cudaFree( d__k_laser ) );
free( uP );
free( vP );
HANDLE_ERROR( cudaFree( d__uP ) );
HANDLE_ERROR( cudaFree( d__vP ) );
free( x );
free( y );
free( d );
HANDLE_ERROR( cudaFree( d__x ) );
HANDLE_ERROR( cudaFree( d__y ) );
HANDLE_ERROR( cudaFree( d__d ) );
#line 3766 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
free( A );
free( A0 );
free( B );
free( B0 );
HANDLE_ERROR( cudaFree( d__A ) );
HANDLE_ERROR( cudaFree( d__A0 ) );
HANDLE_ERROR( cudaFree( d__B ) );
HANDLE_ERROR( cudaFree( d__B0 ) );
#line 3585 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}
void lateral_edge_sensors::data(void)
{
  
#line 4042 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
dim3 blockDim(6,N_HEIGHT);
dim3 gridDim(1,1);
edge_sensors_laser_vector_kernel <<< gridDim, blockDim >>> (d__k_laser, d__k_cam,
                                                            mirror->N-1,
							    mirror->rigid_body_CS.d__R,
							    mirror->motion_CS.d__R);
HANDLE_ERROR( cudaMemcpy(k_laser , d__k_laser,
                          N_DATA*sizeof(vector), cudaMemcpyDeviceToHost ) );
#line 4105 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
laser_coordinates <<< gridDim, blockDim >>> (d__A, mirror->N-1, mirror->D_full,
					     alpha, height,
					     mirror->rigid_body_CS.d__R,
                                             mirror->rigid_body_CS.d__origin,
                                             mirror->motion_CS.d__R,
                                             mirror->motion_CS.d__origin);
int n_byte_DATA;
n_byte_DATA = sizeof(vector)*N_DATA;
HANDLE_ERROR( cudaMemcpy( A, d__A,
                          n_byte_DATA, cudaMemcpyDeviceToHost ) );
#line 4138 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
laser_camera_intersection <<< gridDim, blockDim >>> (d__x, d__y, d__d,
						     d__B, d__A,
						     d__k_cam, d__k_laser,
						     d__uP, d__vP,
						     mirror->N-1);
HANDLE_ERROR( cudaMemcpy( x, d__x,
                          sizeof(rtd)*N_DATA, cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaMemcpy( y, d__y,
                          sizeof(rtd)*N_DATA, cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaMemcpy( d, d__d,
                          sizeof(rtd)*N_DATA, cudaMemcpyDeviceToHost ) );
#line 3589 "/home/ubuntu/projects/crseo/sys/CEO/gmtMirrors/gmtMirrors.nw"
}