#line 112 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
#include "rayTracing.h"

#line 2049 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void transform_to_S_kernel(ray *d__ray, int N_RAY,
				      rtd *d__R, vector *d__origin)
{
  int i, j, ij;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  if ( (ij<N_RAY) && (d__ray[ij].v) )
  {
    forward_transform(&(d__ray[ij].coordinates), &(d__ray[ij].coordinates),
		      d__R, d__origin);
    forward_transform_centered(&(d__ray[ij].directions), &(d__ray[ij].directions),
			       d__R);
  }
}
#line 2112 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void transform_to_R_kernel(ray *d__ray, int N_RAY,
                                      rtd *d__R, vector *d__origin)
{
  int i, j, ij;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  if ( (ij<N_RAY) && (d__ray[ij].v) )
  {
    backward_transform(&(d__ray[ij].coordinates), &(d__ray[ij].coordinates),
		       d__R, d__origin);
    backward_transform_centered(&(d__ray[ij].directions), &(d__ray[ij].directions),
				d__R);
  }
}
#line 2180 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void intersect_chief_kernel(ray *d__ray, const int N_RAY,
				       const rtd Fk, const rtd Fc,
				       vector *d__origin,
				       const int N, rtd *a)
{
  int j, ij;
  rtd s0, s1, x1, y1, k, l, m, S, K, L ,M, dSds;
  vector vv;
  ij = 0;

  
#line 2244 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
k  = d__ray[ij].directions.x;
l  = d__ray[ij].directions.y;
m  = d__ray[ij].directions.z;
if (m==0) { return; }
s0 = -d__ray[ij].coordinates.z/m;
x1 = d__ray[ij].coordinates.x + k*s0;
y1 = d__ray[ij].coordinates.y + l*s0;
#line 2191 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  d__ray[ij].optical_path_length = s0;
  s0 = s1 = 0;
  for (j=0; j<SNELL_N_ITERATION; j++)
    {
      
#line 2273 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
vv.x = x1 + k*s0;
vv.y = y1 + l*s0;
vv.z = m*s0;
#line 2283 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
S = conic_surface(&vv, d__origin, Fk, Fc, N, a);
K = partial_x_conic_surface(&vv, d__origin, Fk, Fc, N, a);
L = partial_y_conic_surface(&vv, d__origin, Fk, Fc, N, a);
M = partial_z_conic_surface();
dSds = K*k + L*l + M*m;
if (dSds==0) { return; }
s1 = s0 - S/dSds;
if (abs(s1-s0)<TOL)
{
  d__ray[ij].coordinates.x = x1 + k*s1;
  d__ray[ij].coordinates.y = y1 + l*s1;
  d__ray[ij].coordinates.z = m*s1;
  d__ray[ij].surface_normal.x    =
    partial_x_conic_surface(&(d__ray[ij].coordinates), d__origin, Fk, Fc, N, a);
  d__ray[ij].surface_normal.y    =
    partial_y_conic_surface(&(d__ray[ij].coordinates), d__origin, Fk, Fc, N, a);
  d__ray[ij].surface_normal.z    = partial_z_conic_surface();
  d__ray[ij].optical_path_length += s1;
  d__ray[ij].n_iteration = j;
  break;
}
#line 2196 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      s0 = s1;
    }
}
__global__ void intersect_kernel(ray *d__ray, int N_RAY,
				 rtd Fk, const rtd Fc,
				 vector *d__origin,
                                 ray *d__chief_ray,
				 int N, rtd *a)
{
  int i, j, ij;
  rtd s0, s1, x1, y1, k, l, m, S, K, L ,M, dSds;
  vector vv;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  if ( (ij<N_RAY) && (d__ray[ij].v) )
  {
    
#line 2244 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
k  = d__ray[ij].directions.x;
l  = d__ray[ij].directions.y;
m  = d__ray[ij].directions.z;
if (m==0) { return; }
s0 = -d__ray[ij].coordinates.z/m;
x1 = d__ray[ij].coordinates.x + k*s0;
y1 = d__ray[ij].coordinates.y + l*s0;
#line 2214 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    d__ray[ij].optical_path_length = s0;
    s0 = s1 = 0;
    for (j=0; j<SNELL_N_ITERATION; j++)
    {
      
#line 2273 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
vv.x = x1 + k*s0;
vv.y = y1 + l*s0;
vv.z = m*s0;
#line 2283 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
S = conic_surface(&vv, d__origin, Fk, Fc, N, a);
K = partial_x_conic_surface(&vv, d__origin, Fk, Fc, N, a);
L = partial_y_conic_surface(&vv, d__origin, Fk, Fc, N, a);
M = partial_z_conic_surface();
dSds = K*k + L*l + M*m;
if (dSds==0) { return; }
s1 = s0 - S/dSds;
if (abs(s1-s0)<TOL)
{
  d__ray[ij].coordinates.x = x1 + k*s1;
  d__ray[ij].coordinates.y = y1 + l*s1;
  d__ray[ij].coordinates.z = m*s1;
  d__ray[ij].surface_normal.x    =
    partial_x_conic_surface(&(d__ray[ij].coordinates), d__origin, Fk, Fc, N, a);
  d__ray[ij].surface_normal.y    =
    partial_y_conic_surface(&(d__ray[ij].coordinates), d__origin, Fk, Fc, N, a);
  d__ray[ij].surface_normal.z    = partial_z_conic_surface();
  d__ray[ij].optical_path_length += s1;
  d__ray[ij].n_iteration = j;
  break;
}
#line 2219 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
	s0 = s1;
    }
  d__ray[ij].optical_path_difference +=
    d__ray[ij].optical_path_length - d__chief_ray->optical_path_length;
  }
}
#line 2487 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    __global__ void reflect_kernel(ray *d__ray, int N_RAY, float mu)
{
  int i, j, ij;
  rtd k, l, m, K, L ,M, G2, a;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  if ( (ij<N_RAY) && (d__ray[ij].v) )
  {
    
#line 2375 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
K = d__ray[ij].surface_normal.x;
L = d__ray[ij].surface_normal.y;
M = d__ray[ij].surface_normal.z;
G2 = K*K + L*L + M*M;
k = d__ray[ij].directions.x;
l = d__ray[ij].directions.y;
m = d__ray[ij].directions.z;
a = mu*(k*K + l*L + m*M)/G2;
#line 2497 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    a *= -2.0;
    d__ray[ij].directions.x = mu*k + a*K;
    d__ray[ij].directions.y = mu*l + a*L;
    d__ray[ij].directions.z = mu*m + a*M;
  }
}
#line 2427 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void refract_kernel(ray *d__ray, int N_RAY, const rtd mu)
{
  int i, j, ij;
  rtd k, l, m, K, L ,M, G2, a, b, gamma0, gamma1, gamma_relative_error;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  if ( (ij<N_RAY) && (d__ray[ij].v) )
  {
    
#line 2375 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
K = d__ray[ij].surface_normal.x;
L = d__ray[ij].surface_normal.y;
M = d__ray[ij].surface_normal.z;
G2 = K*K + L*L + M*M;
k = d__ray[ij].directions.x;
l = d__ray[ij].directions.y;
m = d__ray[ij].directions.z;
a = mu*(k*K + l*L + m*M)/G2;
#line 2437 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    
#line 2389 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
b = (mu*mu-1)/G2;
#line 2438 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    gamma0 = -0.5*b/a;
    for (j=0; j<SNELL_N_ITERATION; j++)
    {
      gamma1 = 0.5*(gamma0*gamma0-b)/(gamma0+a);
      gamma_relative_error = abs( gamma1/gamma0 - 1.0 );
        if (gamma_relative_error<1E-6)
      {
        d__ray[ij].directions.x = mu*k + gamma1*K;
        d__ray[ij].directions.y = mu*l + gamma1*L;
        d__ray[ij].directions.z = mu*m + gamma1*M;
        return;
      }
      gamma0 = gamma1;
    }
  }
}

#line 2066 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__host__ __device__  void forward_transform(vector *v_out, vector *v_in,
					 rtd *d__R, vector *d__origin)
{
  rtd u, v, w;
  u = v_in->x - d__origin->x;
  v = v_in->y - d__origin->y;
  w = v_in->z - d__origin->z;

  v_out->x = d__R[0]*u + d__R[3]*v + d__R[6]*w;
  v_out->y = d__R[1]*u + d__R[4]*v + d__R[7]*w;
  v_out->z = d__R[2]*u + d__R[5]*v + d__R[8]*w;
}
__host__ __device__  void forward_transform_centered(vector *v_out, vector *v_in,
						  rtd *d__R)
{
  rtd u, v, w;
  u = v_in->x;
  v = v_in->y;
  w = v_in->z;

  v_out->x = d__R[0]*u + d__R[3]*v + d__R[6]*w;
  v_out->y = d__R[1]*u + d__R[4]*v + d__R[7]*w;
  v_out->z = d__R[2]*u + d__R[5]*v + d__R[8]*w;
}
#line 2130 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__host__ __device__  void backward_transform(vector *v_out, vector *v_in,
                                    rtd *d__R, vector *d__origin)
{
  rtd u, v, w;

  u = v_in->x;
  v = v_in->y;
  w = v_in->z;

  v_out->x = d__R[0]*u + d__R[1]*v + d__R[2]*w;
  v_out->y = d__R[3]*u + d__R[4]*v + d__R[5]*w;
  v_out->z = d__R[6]*u + d__R[7]*v + d__R[8]*w;

  v_out->x += d__origin->x;
  v_out->y += d__origin->y;
  v_out->z += d__origin->z;
}
__host__ __device__  void backward_transform_centered(vector *v_out, vector *v_in,
						  rtd *d__R)
{
  rtd u, v, w;

  u = v_in->x;
  v = v_in->y;
  w = v_in->z;

  v_out->x = d__R[0]*u + d__R[1]*v + d__R[2]*w;
  v_out->y = d__R[3]*u + d__R[4]*v + d__R[5]*w;
  v_out->z = d__R[6]*u + d__R[7]*v + d__R[8]*w;
}

#line 481 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__host__ __device__ rtd conic_equation(vector *v, vector *v0, 
					      const rtd k, const rtd c,
					      const int N, rtd *a)
{
  rtd rho2, asphere;
  rho2    = v->rho2_shift(v0->x,v0->y);
  asphere = even_asphere(rho2, N, a);
  rho2   *= c;
  if (c==0)
    return 0.0 + asphere;
  else {
    if (k==0)
      return rho2*0.5 + asphere;
    else
      return rho2/( 1 + sqrt(1 - k*c*rho2) ) + asphere;
  }
}
#line 502 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd  conic_surface(vector *v, vector *v0, 
				     const rtd k, const rtd c,
				     const int N, rtd *a)
{
  return v->z - v0->z - conic_equation(v, v0, k, c, N, a);
}
#line 518 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
 __device__ rtd partial_x_conic_surface(vector *v, vector *v0,
					       const rtd k, const rtd c,
					       const int N, rtd *a)
{
  rtd rho2, asphere;
  rho2    = v->rho2_shift(v0->x,v0->y);
  asphere = partial_x_or_y_even_asphere(v->x - v0->x, rho2, N, a);
  rho2   *= c;
  if (c==0)
    return 0.0 - asphere;
  else
    if (k==0)
      return -(v->x - v0->x)*c - asphere;
    else {
      return -(v->x - v0->x)*c*rsqrt(1 - k*c*rho2) - asphere;
    }
}
#line 542 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd partial_y_conic_surface(vector *v, vector *v0,
					      const rtd k, const rtd c,
					      const int N, rtd *a)
{
  rtd rho2, asphere;
  rho2    = v->rho2_shift(v0->x,v0->y);
  asphere = partial_x_or_y_even_asphere(v->y - v0->y, rho2, N, a);
  rho2   *= c;
  if (c==0)
    return 0.0 - asphere;
  else
    if (k==0)
      return -(v->y - v0->y)*c - asphere;
    else {
      return -(v->y - v0->y)*c*rsqrt(1 - k*c*rho2) - asphere;
    }
}
#line 566 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd partial_z_conic_surface(void)
{
  return 1.0;
}
#line 583 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd aspheric_surface(vector *v, vector *v0,
                                       const rtd k, const rtd c,
				       const int max_n, rtd *a, const rtd R)
{
  return conic_surface(v, v0, k, c, 0, a) -
         zernike_surface(v->mag(R), v->angle(), max_n, a);
}
#line 593 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd zernike_surface(vector *v, vector *v0,
				       const int max_n, rtd *a, const rtd R)
{
  return -zernike_surface(v->mag(R), v->angle(), max_n, a);
}
#line 607 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd partial_x_aspheric_surface(vector *v, vector *v0,
					       const rtd k, const rtd c,
					       const int max_n, rtd *cx,
					       const rtd R)
{
  return partial_x_conic_surface(v, v0,	k, c, 0, cx) -
         zernike_surface(v->mag(R), v->angle(), max_n-1, cx);
}
#line 618 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd partial_x_zernike_surface(vector *v, vector *v0,
					       const int max_n, rtd *cx,
					       const rtd R)
{
  return -zernike_surface(v->mag(R), v->angle(), max_n-1, cx);
}
#line 631 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd partial_y_aspheric_surface(vector *v, vector *v0,
					      const rtd k, const rtd c,
					      const int max_n, rtd *cy,
					      const rtd R)
{
  return partial_y_conic_surface(v, v0,	k, c, 0, cy) -
         zernike_surface(v->mag(R), v->angle(), max_n-1, cy);
}
#line 642 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd partial_y_zernike_surface(vector *v, vector *v0,
					      const int max_n, rtd *cy,
					      const rtd R)
{
  return -zernike_surface(v->mag(R), v->angle(), max_n-1, cy);
}
#line 655 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd partial_z_aspheric_surface(void)
{
  return 1.0;
}
#line 662 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__device__ rtd partial_z_zernike_surface(void)
{
  return 0.0;
}
#line 713 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__host__ __device__ rtd even_asphere(rtd r2, int N, rtd *a)
{
  int k;
  rtd A, _r2_;
  if (N==0)
    return 0.0;
  else {
    _r2_ = r2;
    A = a[0]*_r2_;
    for (k=1;k<N;k++) {
      _r2_ *= r2;
      A += a[k]*_r2_;
    }
    return A;
  }
}
#line 745 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__host__ __device__ rtd partial_x_or_y_even_asphere(rtd x_or_y, rtd r2, int N, rtd *a)
{
  
#line 753 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
int k;
rtd A, _r2_;
if (N==0)
  A = 0.0;
else {
  _r2_ = 1;
  A = a[0];
  for (k=1;k<N;k++) {
    _r2_ *= r2;
    A += (k+1)*a[k]*_r2_;
  }
  A *= 2*x_or_y;
}
#line 748 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  return A;
}

#line 1314 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void surface_kern(rtd *S, rtd *r, rtd *o, const int N,
			     int max_n, rtd *a)
{
  int i, j, k;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  k = j * gridDim.x * blockDim.x + i;
  if (k<N) {
    S[k] = zernike_surface(r[k], o[k], max_n, a);
  }
}
#line 1743 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void surface_and_derivatives_kern(rtd *S, rtd *dSdx, rtd *dSdy,
                                             rtd *r, rtd *o, const int N,
                                             int max_n, rtd *a, rtd *ax, rtd *ay)
{
  int i, j, k;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  k = j * gridDim.x * blockDim.x + i;
  if (k<N) {
    zernike_surface_and_partials(S+k, dSdx+k, dSdy+k, r[k], o[k], max_n, a, ax, ay);
  }
}
#line 1785 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void projection_kern(rtd *a, float *phase, rtd *r, rtd *o, const int N,
				int max_n, rtd *ap, int mode)
{
  int i;
  rtd num, denom, z;
  num = denom = 0.0;
  ap[mode] = 1.0;
  for (i=0; i<N; i++) {
    if (phase[i]!=0)
      {
	z = zernike_surface(r[i], o[i], max_n, ap);
	num   += phase[i]*z;
	denom += z*z;
      }
  }
  a[mode]  = num/denom;
  ap[mode] = 0.0;
}
#line 1968 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void intersection(char *mask, ray *d__ray, int N_RAY)
{
  int i, j, ij, iSource;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = threadIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  iSource = blockIdx.z;
  if  ( ij<N_RAY )
  {
    ij += iSource*N_RAY;
    d__ray[ij].v &= mask[ij];
  }
}

#line 200 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::setup(vector _origin_, vector _euler_angles_)
{
  N = 1;
  HANDLE_ERROR( cudaMalloc((void**)&d__R, sizeof(rtd)*9*N ) );
  HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector)*N ) );
  int n_byte = sizeof(vector)*N;
  origin = (vector *)malloc( n_byte );
  euler_angles = (vector *)malloc( n_byte );
  memcpy( origin , &_origin_ , n_byte);
  memcpy(euler_angles,  &_euler_angles_, n_byte);
  
#line 307 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
R = (rtd *)malloc( sizeof(rtd)*9*N );
rtd ca, sa, cb, sb, cg, sg;
int k, idx;
for (k=0; k<N; k++)
  {
    idx = k;
    
#line 320 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
ca = cos(euler_angles[idx].x);
sa = sin(euler_angles[idx].x);
cb = cos(euler_angles[idx].y);
sb = sin(euler_angles[idx].y);
cg = cos(euler_angles[idx].z);
sg = sin(euler_angles[idx].z);
idx *= 9;
R[idx++] = cb*cg;
R[idx++] = sa*sb*cg - ca*sg;
R[idx++] = ca*sb*cg + sa*sg;
R[idx++] = cb*sg;
R[idx++] = sa*sb*sg + ca*cg;
R[idx++] = ca*sb*sg - sa*cg;
R[idx++] = -sb;
R[idx++] = sa*cb;
R[idx++] = ca*cb;
#line 314 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  }
#line 338 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
HANDLE_ERROR( cudaMemcpy( d__R, R,
			    sizeof(rtd)*9*N, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__origin, origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );
#line 316 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
/* for (k=0;k<9*N;k++) */
/*    *(float_R+k) = (float) R[k]; */
#line 211 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  strcpy(tag,"coordinate_system");
  info();
}
#line 218 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::setup(void)
{
  N = 1;
  HANDLE_ERROR( cudaMalloc((void**)&d__R, sizeof(rtd)*9*N ) );
  HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector)*N ) );
  int n_byte = sizeof(vector)*N;
  origin = (vector *)malloc( n_byte );
  euler_angles = (vector *)malloc( n_byte );
  origin[0].x = 0;
  origin[0].y = 0;
  origin[0].z = 0;
  euler_angles[0].x = 0;
  euler_angles[0].y = 0;
  euler_angles[0].z = 0;
  
#line 307 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
R = (rtd *)malloc( sizeof(rtd)*9*N );
rtd ca, sa, cb, sb, cg, sg;
int k, idx;
for (k=0; k<N; k++)
  {
    idx = k;
    
#line 320 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
ca = cos(euler_angles[idx].x);
sa = sin(euler_angles[idx].x);
cb = cos(euler_angles[idx].y);
sb = sin(euler_angles[idx].y);
cg = cos(euler_angles[idx].z);
sg = sin(euler_angles[idx].z);
idx *= 9;
R[idx++] = cb*cg;
R[idx++] = sa*sb*cg - ca*sg;
R[idx++] = ca*sb*cg + sa*sg;
R[idx++] = cb*sg;
R[idx++] = sa*sb*sg + ca*cg;
R[idx++] = ca*sb*sg - sa*cg;
R[idx++] = -sb;
R[idx++] = sa*cb;
R[idx++] = ca*cb;
#line 314 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  }
#line 338 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
HANDLE_ERROR( cudaMemcpy( d__R, R,
			    sizeof(rtd)*9*N, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__origin, origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );
#line 316 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
/* for (k=0;k<9*N;k++) */
/*    *(float_R+k) = (float) R[k]; */
#line 233 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  strcpy(tag,"coordinate_system");
  info();
}
#line 242 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::setup(vector *_origin_, vector *_euler_angles_, int _N_)
{
  strcpy(tag,"coordinate_system");
  
#line 261 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
N = _N_;
HANDLE_ERROR( cudaMalloc((void**)&d__R, sizeof(rtd)*9*N ) );
HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector)*N ) );
int n_byte = sizeof(vector)*N;
origin = (vector *)malloc( n_byte );
euler_angles = (vector *)malloc( n_byte );
memcpy( origin , _origin_ , n_byte);
memcpy(euler_angles,  _euler_angles_, n_byte);
#line 307 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
R = (rtd *)malloc( sizeof(rtd)*9*N );
rtd ca, sa, cb, sb, cg, sg;
int k, idx;
for (k=0; k<N; k++)
  {
    idx = k;
    
#line 320 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
ca = cos(euler_angles[idx].x);
sa = sin(euler_angles[idx].x);
cb = cos(euler_angles[idx].y);
sb = sin(euler_angles[idx].y);
cg = cos(euler_angles[idx].z);
sg = sin(euler_angles[idx].z);
idx *= 9;
R[idx++] = cb*cg;
R[idx++] = sa*sb*cg - ca*sg;
R[idx++] = ca*sb*cg + sa*sg;
R[idx++] = cb*sg;
R[idx++] = sa*sb*sg + ca*cg;
R[idx++] = ca*sb*sg - sa*cg;
R[idx++] = -sb;
R[idx++] = sa*cb;
R[idx++] = ca*cb;
#line 314 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  }
#line 338 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
HANDLE_ERROR( cudaMemcpy( d__R, R,
			    sizeof(rtd)*9*N, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__origin, origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );
#line 316 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
/* for (k=0;k<9*N;k++) */
/*    *(float_R+k) = (float) R[k]; */
#line 270 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
info();
#line 246 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 251 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::setup(vector *_origin_, vector *_euler_angles_,
                              int _N_, char *_tag_)
{
  strcpy(tag,_tag_);
  
#line 261 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
N = _N_;
HANDLE_ERROR( cudaMalloc((void**)&d__R, sizeof(rtd)*9*N ) );
HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector)*N ) );
int n_byte = sizeof(vector)*N;
origin = (vector *)malloc( n_byte );
euler_angles = (vector *)malloc( n_byte );
memcpy( origin , _origin_ , n_byte);
memcpy(euler_angles,  _euler_angles_, n_byte);
#line 307 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
R = (rtd *)malloc( sizeof(rtd)*9*N );
rtd ca, sa, cb, sb, cg, sg;
int k, idx;
for (k=0; k<N; k++)
  {
    idx = k;
    
#line 320 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
ca = cos(euler_angles[idx].x);
sa = sin(euler_angles[idx].x);
cb = cos(euler_angles[idx].y);
sb = sin(euler_angles[idx].y);
cg = cos(euler_angles[idx].z);
sg = sin(euler_angles[idx].z);
idx *= 9;
R[idx++] = cb*cg;
R[idx++] = sa*sb*cg - ca*sg;
R[idx++] = ca*sb*cg + sa*sg;
R[idx++] = cb*sg;
R[idx++] = sa*sb*sg + ca*cg;
R[idx++] = ca*sb*sg - sa*cg;
R[idx++] = -sb;
R[idx++] = sa*cb;
R[idx++] = ca*cb;
#line 314 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  }
#line 338 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
HANDLE_ERROR( cudaMemcpy( d__R, R,
			    sizeof(rtd)*9*N, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__origin, origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );
#line 316 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
/* for (k=0;k<9*N;k++) */
/*    *(float_R+k) = (float) R[k]; */
#line 270 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
info();
#line 256 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 345 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::cleanup(void)
{
  INFO("@(CEO)>coordinate_system: freeing memory!\n");
  free( origin );
  free( euler_angles );
  free(R);
  HANDLE_ERROR( cudaFree( d__R ) );
  HANDLE_ERROR( cudaFree( d__origin ) );
}
#line 358 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::update(void)
{
  int k, idx;
  rtd ca, sa, cb, sb, cg, sg;
  for (k=0; k<N; k++)
    {
      idx = k;
      
#line 320 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
ca = cos(euler_angles[idx].x);
sa = sin(euler_angles[idx].x);
cb = cos(euler_angles[idx].y);
sb = sin(euler_angles[idx].y);
cg = cos(euler_angles[idx].z);
sg = sin(euler_angles[idx].z);
idx *= 9;
R[idx++] = cb*cg;
R[idx++] = sa*sb*cg - ca*sg;
R[idx++] = ca*sb*cg + sa*sg;
R[idx++] = cb*sg;
R[idx++] = sa*sb*sg + ca*cg;
R[idx++] = ca*sb*sg - sa*cg;
R[idx++] = -sb;
R[idx++] = sa*cb;
R[idx++] = ca*cb;
#line 366 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    }
  
#line 338 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
HANDLE_ERROR( cudaMemcpy( d__R, R,
			    sizeof(rtd)*9*N, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__origin, origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );
#line 368 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 371 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::update(vector _origin_, vector _euler_angles_,int idx)
{
  int n_byte = sizeof(vector);
  memcpy( origin + idx , &_origin_ , n_byte);
  memcpy(euler_angles + idx,  &_euler_angles_, n_byte);
  rtd ca, sa, cb, sb, cg, sg;
  
#line 320 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
ca = cos(euler_angles[idx].x);
sa = sin(euler_angles[idx].x);
cb = cos(euler_angles[idx].y);
sb = sin(euler_angles[idx].y);
cg = cos(euler_angles[idx].z);
sg = sin(euler_angles[idx].z);
idx *= 9;
R[idx++] = cb*cg;
R[idx++] = sa*sb*cg - ca*sg;
R[idx++] = ca*sb*cg + sa*sg;
R[idx++] = cb*sg;
R[idx++] = sa*sb*sg + ca*cg;
R[idx++] = ca*sb*sg - sa*cg;
R[idx++] = -sb;
R[idx++] = sa*cb;
R[idx++] = ca*cb;
#line 378 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  
#line 338 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
HANDLE_ERROR( cudaMemcpy( d__R, R,
			    sizeof(rtd)*9*N, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__origin, origin,
			    sizeof(vector)*N, cudaMemcpyHostToDevice ) );
#line 379 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 387 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::info_details(void)
{
  INFO("\n\x1B[1;42m@(CEO)>%s:\x1B[;42m\n",tag);
  INFO(" ID (Tx          Ty        Tz)[m]  (Ox         Oy         Oz)[deg]\n");
  float c = 180.0/PI;
  int k;//, i, j, idx;
  for (k=0; k<N; k++)
    {
      INFO(" %2d  %+.3e   %+.3e   %+.3e      %+.3e   %+.3e   %+.3e\n", k,
	      origin[k].x, origin[k].y, origin[k].z,
	      euler_angles[k].x*c, euler_angles[k].y*c, euler_angles[k].z*c);
      /*
      INFO("Direction cosine matrix:\n");
      for (i=0;i<3;i++)
	{
	  INFO("||  ");
	  for (j=0;j<3;j++)
	    {
	      idx = j + i*3 + k*9;
	      INFO("%+.3e  ",R[idx]);
	    }
	  INFO("||\n");
	}
      */
    }
  INFO("----------------------------------------------------\x1B[0m\n");
}
#line 417 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void coordinate_system::info(void)
{
  //info_details();
/*
  INFO("\n\x1B[1;42m@(CEO)>%s:\x1B[;42m\n",tag);
  INFO(" %d coordinate systems\n",N);
  INFO("----------------------------------------------------\x1B[0m\n");
*/
}

#line 773 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void conic::setup(rtd _c_, rtd _k_, vector _origin_, vector _euler_angles_)
{
  c = _c_;
  k = _k_;
  origin.x = 0;
  origin.y = 0;
  origin.z = 0;
  ref_frame.setup(_origin_, _euler_angles_);
  
#line 847 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
refractive_index = 1.0;
HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector) ) );
HANDLE_ERROR( cudaMemcpy( d__origin, &origin,
			    sizeof(vector), cudaMemcpyHostToDevice ) );
even_asphere_N = 0;
d__even_asphere_a = NULL;
#line 782 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 786 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void conic::setup(rtd _c_, rtd _k_, vector _origin_, vector _euler_angles_, vector conic_origin)
{
  c = _c_;
  k = _k_;
  origin.x = conic_origin.x;
  origin.y = conic_origin.y;
  origin.z = conic_origin.z;
  ref_frame.setup(_origin_, _euler_angles_);
  
#line 847 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
refractive_index = 1.0;
HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector) ) );
HANDLE_ERROR( cudaMemcpy( d__origin, &origin,
			    sizeof(vector), cudaMemcpyHostToDevice ) );
even_asphere_N = 0;
d__even_asphere_a = NULL;
#line 795 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  info();
}
#line 799 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void conic::setup(rtd _c_, rtd _k_, vector _origin_, vector _euler_angles_, 
                  vector conic_origin, rtd _refractive_index_)
{
  c = _c_;
  k = _k_;
  origin.x = conic_origin.x;
  origin.y = conic_origin.y;
  origin.z = conic_origin.z;
  ref_frame.setup(_origin_, _euler_angles_);
  
#line 847 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
refractive_index = 1.0;
HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector) ) );
HANDLE_ERROR( cudaMemcpy( d__origin, &origin,
			    sizeof(vector), cudaMemcpyHostToDevice ) );
even_asphere_N = 0;
d__even_asphere_a = NULL;
#line 809 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  refractive_index = _refractive_index_;
}
#line 813 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void conic::setup(rtd _c_, rtd _k_, vector _origin_, vector _euler_angles_, 
                  vector conic_origin, rtd _refractive_index_,
                  int asphere_N, rtd *asphere_a)
{
  c = _c_;
  k = _k_;
  origin.x = conic_origin.x;
  origin.y = conic_origin.y;
  origin.z = conic_origin.z;
  ref_frame.setup(_origin_, _euler_angles_);
  
#line 847 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
refractive_index = 1.0;
HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector) ) );
HANDLE_ERROR( cudaMemcpy( d__origin, &origin,
			    sizeof(vector), cudaMemcpyHostToDevice ) );
even_asphere_N = 0;
d__even_asphere_a = NULL;
#line 824 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  refractive_index = _refractive_index_;
  even_asphere_N = asphere_N;
  HANDLE_ERROR( cudaMalloc((void**)&d__even_asphere_a, sizeof(rtd)*even_asphere_N ) );
  HANDLE_ERROR( cudaMemcpy( d__even_asphere_a, &asphere_a,
			    sizeof(rtd)*even_asphere_N, cudaMemcpyHostToDevice ) );
}
#line 834 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void conic::setup(rtd _c_, rtd _k_)
{
  c = _c_;
  k = _k_;
  origin.x = 0;
  origin.y = 0;
  origin.z = 0;
  ref_frame.setup();
  
#line 847 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
refractive_index = 1.0;
HANDLE_ERROR( cudaMalloc((void**)&d__origin, sizeof(vector) ) );
HANDLE_ERROR( cudaMemcpy( d__origin, &origin,
			    sizeof(vector), cudaMemcpyHostToDevice ) );
even_asphere_N = 0;
d__even_asphere_a = NULL;
#line 843 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 888 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void conic::cleanup(void)
{
  INFO("@(CEO)>conic: freeing memory!\n");
  ref_frame.cleanup();
  HANDLE_ERROR( cudaFree( d__origin ) );
  if (d__even_asphere_a!=NULL)
    HANDLE_ERROR( cudaFree( d__even_asphere_a ) );
}
#line 900 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void conic::trace(bundle *rays)
{
  transform_to_S_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1,
                                       ref_frame.d__R, ref_frame.d__origin);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  transform_to_S_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
					            ref_frame.d__R, ref_frame.d__origin);
  intersect_chief_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1,
                                        k, c, d__origin,
					even_asphere_N, d__even_asphere_a);
  intersect_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
					       k, c, d__origin,
					       rays->d__chief_ray,
					       even_asphere_N, d__even_asphere_a);
  if (rays->refractive_index==refractive_index) 
     reflect(rays);
   else { 
     refract(rays, rays->refractive_index/refractive_index);
     rays->refractive_index = refractive_index;
   }
   transform_to_R_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1,
				       ref_frame.d__R, ref_frame.d__origin);
  transform_to_R_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
        	 				    ref_frame.d__R, ref_frame.d__origin);
}
#line 928 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void conic::info(void)
{
  INFO("\n\x1B[1;42m@(CEO)>conic:\x1B[;42m\n");
  INFO(" . curvature       : %g\n",c);
  INFO(" . conic constant  : %g\n",k);
  INFO(" . refractive index: %g\n",refractive_index);
  INFO("----------------------------------------------------\x1B[0m\n");  
  ref_frame.info();
}

#line 1044 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::setup(int _max_n_, rtd *_a_)
{
  N = 1;
  max_n = _max_n_;
  n_mode = (max_n+1.0)*(max_n+2.0)*0.5;
  int n_byte = sizeof(rtd)*n_mode;
  a = (rtd *)malloc( n_byte );
  memcpy( a, _a_,  n_byte);
  HANDLE_ERROR( cudaMalloc((void**)&d__a, n_byte ) );
  HANDLE_ERROR( cudaMemcpy( d__a, a,
			    n_byte, cudaMemcpyHostToDevice ) );

  ref_frame.setup();

  
#line 1122 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
if (max_n==0)
  return;

INFO("Computing Zernike derivative coefficients: ....");

int nel = n_mode*n_mode/2;
int n_col_ptr = (max_n+1.0)*max_n*0.5 + 1;

bx = (rtd *) malloc( sizeof(rtd)*nel);
bx_row_idx = (unsigned int *) malloc( sizeof(unsigned int)*nel);
bx_col_ptr = (unsigned int *) malloc( sizeof(unsigned int)*n_col_ptr);
memset( bx_col_ptr, 0 , sizeof(unsigned int)*n_col_ptr );

by = (rtd *) malloc( sizeof(rtd)*nel);
by_row_idx = (unsigned int *) malloc( sizeof(unsigned int)*nel);
by_col_ptr = (unsigned int *) malloc( sizeof(unsigned int)*n_col_ptr);
memset( by_col_ptr, 0 , sizeof(unsigned int)*n_col_ptr );

#line 1179 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
int j, n, m, jp, np, mp, kx=0, ky=0;
for (np=0;np<max_n;np++) {
  jp = np*(np+1)/2 + 1;
  mp = np%2;
  while (mp<=np) {
    
#line 1199 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
for (n=np+1;n<=max_n;n++) {
  j = n*(n+1)/2 + 1;
  m = n%2;
  while (m<=n) {
    
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1204 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    ++j;
    if (m>0) {
      
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++j;
    }
    m += 2;
  }
 }
bx_col_ptr[jp] = kx;
by_col_ptr[jp] = ky;
#line 1185 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++jp;
    if (mp>0) {
       
#line 1199 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
for (n=np+1;n<=max_n;n++) {
  j = n*(n+1)/2 + 1;
  m = n%2;
  while (m<=n) {
    
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1204 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    ++j;
    if (m>0) {
      
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++j;
    }
    m += 2;
  }
 }
bx_col_ptr[jp] = kx;
by_col_ptr[jp] = ky;
#line 1188 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
	++jp;
    }
    mp += 2;
  }
 }
bx_col_ptr[jp-1] = kx;
bx_nnz = kx;
by_col_ptr[jp-1] = ky;
by_nnz = ky;
#line 1141 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
INFO("\b\b");
/*
INFO("bx_nnz=%d\n",bx_nnz);
INFO("b ; row_idx ; col_ptr\n");
for (int k=0;k<bx_nnz;k++) {
  if (k<n_col_ptr)
  INFO("b=%+4.2f ; %2d ; %2d\n",bx[k],bx_row_idx[k],bx_col_ptr[k]);
  else
  INFO("b=%+4.2f ; %2d\n",bx[k],bx_row_idx[k]);
}
INFO("by_nnz=%d\n",by_nnz);
INFO("b ; row_idy ; col_ptr\n");
for (int k=0;k<by_nnz;k++) {
  if (k<n_col_ptr)
  INFO("b=%+4.2f ; %2d ; %2d\n",by[k],by_row_idx[k],by_col_ptr[k]);
  else
  INFO("b=%+4.2f ; %2d\n",by[k],by_row_idx[k]);
}
*/

n_byte = N*sizeof(rtd)*(n_col_ptr-1);
cx = (rtd *) malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__cx, n_byte ) );

cy = (rtd *) malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__cy, n_byte ) );

surface_derivative_coefs();
INFO("\b\b");

HANDLE_ERROR( cudaMemcpy( d__cx, cx,
			  n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__cy, cy,
			  n_byte, cudaMemcpyHostToDevice ) );
INFO("\b\b!\n");
#line 1059 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 1064 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::setup(int _max_n_, rtd *_a_, vector _origin_, vector _euler_angles_)
{
  N = 1;
  max_n = _max_n_;
  n_mode = (max_n+1.0)*(max_n+2.0)*0.5;
  int n_byte = sizeof(rtd)*n_mode;
  a = (rtd *)malloc( n_byte );
  memcpy( a, _a_,  n_byte);
  HANDLE_ERROR( cudaMalloc((void**)&d__a, n_byte ) );
  HANDLE_ERROR( cudaMemcpy( d__a, a,
			    n_byte, cudaMemcpyHostToDevice ) );

  ref_frame.setup(_origin_, _euler_angles_);

  
#line 1122 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
if (max_n==0)
  return;

INFO("Computing Zernike derivative coefficients: ....");

int nel = n_mode*n_mode/2;
int n_col_ptr = (max_n+1.0)*max_n*0.5 + 1;

bx = (rtd *) malloc( sizeof(rtd)*nel);
bx_row_idx = (unsigned int *) malloc( sizeof(unsigned int)*nel);
bx_col_ptr = (unsigned int *) malloc( sizeof(unsigned int)*n_col_ptr);
memset( bx_col_ptr, 0 , sizeof(unsigned int)*n_col_ptr );

by = (rtd *) malloc( sizeof(rtd)*nel);
by_row_idx = (unsigned int *) malloc( sizeof(unsigned int)*nel);
by_col_ptr = (unsigned int *) malloc( sizeof(unsigned int)*n_col_ptr);
memset( by_col_ptr, 0 , sizeof(unsigned int)*n_col_ptr );

#line 1179 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
int j, n, m, jp, np, mp, kx=0, ky=0;
for (np=0;np<max_n;np++) {
  jp = np*(np+1)/2 + 1;
  mp = np%2;
  while (mp<=np) {
    
#line 1199 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
for (n=np+1;n<=max_n;n++) {
  j = n*(n+1)/2 + 1;
  m = n%2;
  while (m<=n) {
    
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1204 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    ++j;
    if (m>0) {
      
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++j;
    }
    m += 2;
  }
 }
bx_col_ptr[jp] = kx;
by_col_ptr[jp] = ky;
#line 1185 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++jp;
    if (mp>0) {
       
#line 1199 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
for (n=np+1;n<=max_n;n++) {
  j = n*(n+1)/2 + 1;
  m = n%2;
  while (m<=n) {
    
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1204 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    ++j;
    if (m>0) {
      
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++j;
    }
    m += 2;
  }
 }
bx_col_ptr[jp] = kx;
by_col_ptr[jp] = ky;
#line 1188 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
	++jp;
    }
    mp += 2;
  }
 }
bx_col_ptr[jp-1] = kx;
bx_nnz = kx;
by_col_ptr[jp-1] = ky;
by_nnz = ky;
#line 1141 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
INFO("\b\b");
/*
INFO("bx_nnz=%d\n",bx_nnz);
INFO("b ; row_idx ; col_ptr\n");
for (int k=0;k<bx_nnz;k++) {
  if (k<n_col_ptr)
  INFO("b=%+4.2f ; %2d ; %2d\n",bx[k],bx_row_idx[k],bx_col_ptr[k]);
  else
  INFO("b=%+4.2f ; %2d\n",bx[k],bx_row_idx[k]);
}
INFO("by_nnz=%d\n",by_nnz);
INFO("b ; row_idy ; col_ptr\n");
for (int k=0;k<by_nnz;k++) {
  if (k<n_col_ptr)
  INFO("b=%+4.2f ; %2d ; %2d\n",by[k],by_row_idx[k],by_col_ptr[k]);
  else
  INFO("b=%+4.2f ; %2d\n",by[k],by_row_idx[k]);
}
*/

n_byte = N*sizeof(rtd)*(n_col_ptr-1);
cx = (rtd *) malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__cx, n_byte ) );

cy = (rtd *) malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__cy, n_byte ) );

surface_derivative_coefs();
INFO("\b\b");

HANDLE_ERROR( cudaMemcpy( d__cx, cx,
			  n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__cy, cy,
			  n_byte, cudaMemcpyHostToDevice ) );
INFO("\b\b!\n");
#line 1079 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 1084 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::setup(int _max_n_, rtd *_a_, int _N_)
{
  N = _N_;
  max_n = _max_n_;
  n_mode = (max_n+1.0)*(max_n+2.0)*0.5;
  int n_byte = sizeof(rtd)*n_mode*N;
  a = (rtd *)malloc( n_byte );
  memcpy( a, _a_,  n_byte);
  HANDLE_ERROR( cudaMalloc((void**)&d__a, n_byte ) );
  HANDLE_ERROR( cudaMemcpy( d__a, a,
			    n_byte, cudaMemcpyHostToDevice ) );

  ref_frame.setup();

  
#line 1122 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
if (max_n==0)
  return;

INFO("Computing Zernike derivative coefficients: ....");

int nel = n_mode*n_mode/2;
int n_col_ptr = (max_n+1.0)*max_n*0.5 + 1;

bx = (rtd *) malloc( sizeof(rtd)*nel);
bx_row_idx = (unsigned int *) malloc( sizeof(unsigned int)*nel);
bx_col_ptr = (unsigned int *) malloc( sizeof(unsigned int)*n_col_ptr);
memset( bx_col_ptr, 0 , sizeof(unsigned int)*n_col_ptr );

by = (rtd *) malloc( sizeof(rtd)*nel);
by_row_idx = (unsigned int *) malloc( sizeof(unsigned int)*nel);
by_col_ptr = (unsigned int *) malloc( sizeof(unsigned int)*n_col_ptr);
memset( by_col_ptr, 0 , sizeof(unsigned int)*n_col_ptr );

#line 1179 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
int j, n, m, jp, np, mp, kx=0, ky=0;
for (np=0;np<max_n;np++) {
  jp = np*(np+1)/2 + 1;
  mp = np%2;
  while (mp<=np) {
    
#line 1199 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
for (n=np+1;n<=max_n;n++) {
  j = n*(n+1)/2 + 1;
  m = n%2;
  while (m<=n) {
    
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1204 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    ++j;
    if (m>0) {
      
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++j;
    }
    m += 2;
  }
 }
bx_col_ptr[jp] = kx;
by_col_ptr[jp] = ky;
#line 1185 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++jp;
    if (mp>0) {
       
#line 1199 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
for (n=np+1;n<=max_n;n++) {
  j = n*(n+1)/2 + 1;
  m = n%2;
  while (m<=n) {
    
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1204 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    ++j;
    if (m>0) {
      
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++j;
    }
    m += 2;
  }
 }
bx_col_ptr[jp] = kx;
by_col_ptr[jp] = ky;
#line 1188 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
	++jp;
    }
    mp += 2;
  }
 }
bx_col_ptr[jp-1] = kx;
bx_nnz = kx;
by_col_ptr[jp-1] = ky;
by_nnz = ky;
#line 1141 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
INFO("\b\b");
/*
INFO("bx_nnz=%d\n",bx_nnz);
INFO("b ; row_idx ; col_ptr\n");
for (int k=0;k<bx_nnz;k++) {
  if (k<n_col_ptr)
  INFO("b=%+4.2f ; %2d ; %2d\n",bx[k],bx_row_idx[k],bx_col_ptr[k]);
  else
  INFO("b=%+4.2f ; %2d\n",bx[k],bx_row_idx[k]);
}
INFO("by_nnz=%d\n",by_nnz);
INFO("b ; row_idy ; col_ptr\n");
for (int k=0;k<by_nnz;k++) {
  if (k<n_col_ptr)
  INFO("b=%+4.2f ; %2d ; %2d\n",by[k],by_row_idx[k],by_col_ptr[k]);
  else
  INFO("b=%+4.2f ; %2d\n",by[k],by_row_idx[k]);
}
*/

n_byte = N*sizeof(rtd)*(n_col_ptr-1);
cx = (rtd *) malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__cx, n_byte ) );

cy = (rtd *) malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__cy, n_byte ) );

surface_derivative_coefs();
INFO("\b\b");

HANDLE_ERROR( cudaMemcpy( d__cx, cx,
			  n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__cy, cy,
			  n_byte, cudaMemcpyHostToDevice ) );
INFO("\b\b!\n");
#line 1099 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 1103 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::setup(int _max_n_, rtd *_a_,
		     vector _origin_, vector _euler_angles_, int _N_)
{
  N = _N_;
  max_n = _max_n_;
  n_mode = (max_n+1.0)*(max_n+2.0)*0.5;
  int n_byte = sizeof(rtd)*n_mode*N;
  a = (rtd *)malloc( n_byte );
  memcpy( a, _a_,  n_byte);
  HANDLE_ERROR( cudaMalloc((void**)&d__a, n_byte ) );
  HANDLE_ERROR( cudaMemcpy( d__a, a,
			    n_byte, cudaMemcpyHostToDevice ) );

  ref_frame.setup(_origin_, _euler_angles_);

  
#line 1122 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
if (max_n==0)
  return;

INFO("Computing Zernike derivative coefficients: ....");

int nel = n_mode*n_mode/2;
int n_col_ptr = (max_n+1.0)*max_n*0.5 + 1;

bx = (rtd *) malloc( sizeof(rtd)*nel);
bx_row_idx = (unsigned int *) malloc( sizeof(unsigned int)*nel);
bx_col_ptr = (unsigned int *) malloc( sizeof(unsigned int)*n_col_ptr);
memset( bx_col_ptr, 0 , sizeof(unsigned int)*n_col_ptr );

by = (rtd *) malloc( sizeof(rtd)*nel);
by_row_idx = (unsigned int *) malloc( sizeof(unsigned int)*nel);
by_col_ptr = (unsigned int *) malloc( sizeof(unsigned int)*n_col_ptr);
memset( by_col_ptr, 0 , sizeof(unsigned int)*n_col_ptr );

#line 1179 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
int j, n, m, jp, np, mp, kx=0, ky=0;
for (np=0;np<max_n;np++) {
  jp = np*(np+1)/2 + 1;
  mp = np%2;
  while (mp<=np) {
    
#line 1199 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
for (n=np+1;n<=max_n;n++) {
  j = n*(n+1)/2 + 1;
  m = n%2;
  while (m<=n) {
    
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1204 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    ++j;
    if (m>0) {
      
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++j;
    }
    m += 2;
  }
 }
bx_col_ptr[jp] = kx;
by_col_ptr[jp] = ky;
#line 1185 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++jp;
    if (mp>0) {
       
#line 1199 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
for (n=np+1;n<=max_n;n++) {
  j = n*(n+1)/2 + 1;
  m = n%2;
  while (m<=n) {
    
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1204 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    ++j;
    if (m>0) {
      
#line 1216 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (zern_dx_coef(bx+kx,j,n,m,jp,np,mp)) {
    bx_row_idx[kx] = j-1;
    ++kx;
  }
  if (zern_dy_coef(by+ky,j,n,m,jp,np,mp)) {
    by_row_idx[ky] = j-1;
    ++ky;
  }
#line 1207 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      ++j;
    }
    m += 2;
  }
 }
bx_col_ptr[jp] = kx;
by_col_ptr[jp] = ky;
#line 1188 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
	++jp;
    }
    mp += 2;
  }
 }
bx_col_ptr[jp-1] = kx;
bx_nnz = kx;
by_col_ptr[jp-1] = ky;
by_nnz = ky;
#line 1141 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
INFO("\b\b");
/*
INFO("bx_nnz=%d\n",bx_nnz);
INFO("b ; row_idx ; col_ptr\n");
for (int k=0;k<bx_nnz;k++) {
  if (k<n_col_ptr)
  INFO("b=%+4.2f ; %2d ; %2d\n",bx[k],bx_row_idx[k],bx_col_ptr[k]);
  else
  INFO("b=%+4.2f ; %2d\n",bx[k],bx_row_idx[k]);
}
INFO("by_nnz=%d\n",by_nnz);
INFO("b ; row_idy ; col_ptr\n");
for (int k=0;k<by_nnz;k++) {
  if (k<n_col_ptr)
  INFO("b=%+4.2f ; %2d ; %2d\n",by[k],by_row_idx[k],by_col_ptr[k]);
  else
  INFO("b=%+4.2f ; %2d\n",by[k],by_row_idx[k]);
}
*/

n_byte = N*sizeof(rtd)*(n_col_ptr-1);
cx = (rtd *) malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__cx, n_byte ) );

cy = (rtd *) malloc( n_byte );
HANDLE_ERROR( cudaMalloc((void**)&d__cy, n_byte ) );

surface_derivative_coefs();
INFO("\b\b");

HANDLE_ERROR( cudaMemcpy( d__cx, cx,
			  n_byte, cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( d__cy, cy,
			  n_byte, cudaMemcpyHostToDevice ) );
INFO("\b\b!\n");
#line 1119 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
}
#line 1228 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::surface_derivative_coefs(void)
{
  int k,l,n,p, idx_a, idx_c;
  n = (max_n+1.0)*max_n*0.5;
  for (p=0;p<N;p++) {
    idx_a = p*n_mode;
    idx_c = p*n;
    for (k=0;k<n;k++) {
      cx[idx_c] = 0.0;
      for (l=bx_col_ptr[k];l<bx_col_ptr[k+1];l++)
	cx[idx_c] += a[ bx_row_idx[l] + idx_a ] * bx[l];
      /* INFO("cx = %4.2f\n",cx[k]); */
      cy[idx_c] = 0.0;
      for (l=by_col_ptr[k];l<by_col_ptr[k+1];l++)
	cy[idx_c] += a[ by_row_idx[l] + idx_a ] * by[l];
      /* INFO("cy = %4.2f\n",cy[k]); */
      idx_c++;
    }
  }
}
#line 1252 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::update(rtd *_a_)
{
  int n_byte = sizeof(rtd)*n_mode*N;
  memcpy( a, _a_,  n_byte);
  HANDLE_ERROR( cudaMemcpy( d__a, a,
			    n_byte, cudaMemcpyHostToDevice ) );
  surface_derivative_coefs();

  int n_col_ptr = (max_n+1.0)*max_n*0.5 + 1;
  n_byte = N*sizeof(rtd)*(n_col_ptr-1);
  HANDLE_ERROR( cudaMemcpy( d__cx, cx,
			    n_byte, cudaMemcpyHostToDevice ) );
  HANDLE_ERROR( cudaMemcpy( d__cy, cy,
			    n_byte, cudaMemcpyHostToDevice ) );
}
#line 1271 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::cleanup(void)
{
  INFO("@(CEO)>zernikeS: freeing memory!\n");
  free( a );
  HANDLE_ERROR( cudaFree( d__a) );
  if (max_n>0) {
    free( bx );
    free( bx_row_idx );
    free( bx_col_ptr );
    free( by );
    free( by_row_idx );
    free( by_col_ptr );
    free( cx );
    free( cy );
    HANDLE_ERROR( cudaFree( d__cx) );
    HANDLE_ERROR( cudaFree( d__cy) );
  }
}
#line 1295 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::surface(rtd *d__S, rtd *d__r, rtd *d__o, const int N)
{
    dim3 blockDim(N_THREAD,N_THREAD);
    dim3 gridDim(N/N_THREAD2+1,1);
    surface_kern <<< gridDim, blockDim >>> (d__S, d__r, d__o, N, max_n, d__a);
}
#line 1304 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::surface(rtd *d__S, rtd *d__r, rtd *d__o, const int N, int surf_id)
{
  //INFO("a offset = %d\n",surf_id*n_mode);
    dim3 blockDim(N_THREAD,N_THREAD);
    dim3 gridDim(N/N_THREAD2+1,1);
    surface_kern <<< gridDim, blockDim >>> (d__S, d__r, d__o, N, max_n, d__a + surf_id*n_mode);

}
#line 1661 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::surface_derivatives(rtd *d__dSdx, rtd *d__dSdy,
				   rtd *d__r, rtd *d__o, const int N)
{
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(N/N_THREAD2+1,1);
  surface_kern <<< gridDim, blockDim >>> (d__dSdx, d__r, d__o, N, max_n-1, d__cx);
  surface_kern <<< gridDim, blockDim >>> (d__dSdy, d__r, d__o, N, max_n-1, d__cy);
}
#line 1716 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::surface_and_derivatives(rtd *d__S, rtd *d__dSdx, rtd *d__dSdy,
                                       rtd *d__r, rtd *d__o, const int N)
{
    dim3 blockDim(N_THREAD,N_THREAD);
    dim3 gridDim(N/N_THREAD2+1,1);
    surface_and_derivatives_kern <<< gridDim, blockDim >>> (d__S, d__dSdx, d__dSdy,
                                                            d__r, d__o, N, max_n,
                                                            d__a, d__cx, d__cy);
}
#line 1728 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::surface_and_derivatives(rtd *d__S, rtd *d__dSdx, rtd *d__dSdy,
                                       rtd *d__r, rtd *d__o, const int N, int surf_id)
{
  //INFO("a offset = %d\n",surf_id*n_mode);
    dim3 blockDim(N_THREAD,N_THREAD);
    dim3 gridDim(N/N_THREAD2+1,1);
    surface_and_derivatives_kern <<< gridDim, blockDim >>> (d__S, d__dSdx, d__dSdy,
                                                            d__r, d__o, N, max_n,
                                                            d__a + surf_id*n_mode,
                                                            d__cx + surf_id*n_mode,
                                                            d__cy + surf_id*n_mode);

}
#line 1768 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void zernikeS::projection(float *d__phase, rtd *d__r, rtd *d__o, const int N)
{
  int k;
  int n_byte = sizeof(rtd)*n_mode;
  rtd *d__ap;
  HANDLE_ERROR( cudaMalloc((void**)&d__ap, n_byte ) );
  HANDLE_ERROR( cudaMemset(d__ap, 0, n_byte ) );
  dim3 blockDim(1,1);
  dim3 gridDim(1,1);
  for (k=0;k<n_mode;k++)
    projection_kern <<< gridDim, blockDim >>> (d__a, d__phase, d__r, d__o, N, max_n, d__ap, k);
  HANDLE_ERROR( cudaFree( d__ap) );
  HANDLE_ERROR( cudaMemcpy( a, d__a,
			    sizeof(rtd)*n_mode, cudaMemcpyDeviceToHost ) );
}

#line 1328 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__host__ __device__ rtd zernike_surface(rtd r, rtd o, int max_n, rtd *a)
{
  int n, m, mH, a_idx,j;
  rtd R, R2, R4, A, S_n, S,
    H1, H2, H3, sm, cm, s2, c2, s ,c, sqrt2;
  if (max_n<0)//( (r>1.0) || (max_n<0) )
    S = 0.0;
  else {
    R = 1.0;
    S = a[0]*R;
    
#line 1370 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
sqrt2 = sqrt(2.0);
sincos(2*o,&s2,&c2);
for (n=1;n<=max_n;n++) {
  j = (n+1)*(n+2)/2;
  a_idx = j - 1;

  m = n;
  R = pow(r,n);
  sincos(m*o, &sm, &cm);
  if (j%2) {
    A = a[a_idx--]*sm;
    A += a[a_idx--]*cm;
  } else {
    A = a[a_idx--]*cm;
    A += a[a_idx--]*sm;
  }
  j -= 2;
  A *= sqrt2;
  S_n = R*A;

  m -= 2;
  if (m>=0) {
    R2 = n*R - (n-1.0)*pow(r,m);
    
#line 1426 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (m>0) {
    c = cm*c2 + sm*s2;
    s = sm*c2 - cm*s2;
    if (j%2) {
      A = a[a_idx--]*s;
      A += a[a_idx--]*c;
    } else {
      A = a[a_idx--]*c;
      A += a[a_idx--]*s;
    }
    j -= 2;
    A *= sqrt2;
    cm = c;
    sm = s;
  } else {
    A = a[a_idx--];
    --j;
  }
#line 1394 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    S_n += R2*A;
  }

  m -= 2;
  while (m>=0) {
    mH = m + 4;
    if (r>0) {
      
#line 1485 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H3 = -4.0*(mH-2.0)*(mH-3.0)/((n+mH-2.0)*(n-mH+4.0));
#line 1402 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      
#line 1477 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H2 = H3*0.25*(n+mH)*(n-mH+2.0)/(mH-1.0) + mH - 2.0;
#line 1403 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      
#line 1470 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H1 = mH*(mH-1.0)*0.5 - mH*H2 + H3*(n+mH+2.0)*(n-mH)*0.125;
#line 1404 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      
#line 1451 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
R4 = H1*R + ( H2 + H3/(r*r) )*R2;
#line 1405 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      
#line 1426 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (m>0) {
    c = cm*c2 + sm*s2;
    s = sm*c2 - cm*s2;
    if (j%2) {
      A = a[a_idx--]*s;
      A += a[a_idx--]*c;
    } else {
      A = a[a_idx--]*c;
      A += a[a_idx--]*s;
    }
    j -= 2;
    A *= sqrt2;
    cm = c;
    sm = s;
  } else {
    A = a[a_idx--];
    --j;
  }
#line 1406 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    } else {
      R4 = (m==0) ? pow(-1.0,0.5*n) : 0.0;
      A = 1.0;
    }
    S_n += R4*A;
    R  = R2;
    R2 = R4;
    m -= 2;
  }
  S += sqrt(n+1.0)*S_n;
}
#line 1339 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    }
  return S;
}
#line 1672 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
int zern_dx_coef(rtd *b, int j, int n, int m, int jp, int np, int mp)
{
  int delta_m, rule_0, rule_1, rule_2, rule_3;
  delta_m = mp - m;
  rule_0 = (delta_m==1) || (delta_m==-1);
  rule_1 = (j-jp)%2 == 0 && !( (m==0) || (mp==0) );
  rule_2 = ( (m==0) && (mp!=0) ) && (jp%2 == 0);
  rule_3 = ( (m!=0) && (mp==0) ) && (j%2  == 0);
  if ( rule_0 && (rule_1||rule_2||rule_3) )
    *b = sqrt( (n+1.0)*(np+1.0) );
  else
    return 0;
  if ( (m==0) || (mp==0) )
    *b *= sqrt(2.0);
  return 1;
}
#line 1690 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
int zern_dy_coef(rtd *b, int j, int n, int m, int jp, int np, int mp)
{
  int delta_m, rule_0, rule_1, rule_2, rule_3, sign_rule_1, sign_rule_2;
  delta_m = mp - m;
  rule_0 = (delta_m==1) || (delta_m==-1);
  rule_1 = (j-jp)%2 == 1 && !( (m==0) || (mp==0) );
  rule_2 = ( (m==0) && (mp!=0) ) && (jp%2 == 1);
  rule_3 = ( (m!=0) && (mp==0) ) && (j%2  == 1);
  if ( rule_0 && (rule_1||rule_2||rule_3) )
    *b = sqrt( (n+1.0)*(np+1.0) );
  else
    return 0;
  sign_rule_1 = (delta_m==+1) && (j%2 == 1);
  sign_rule_2 = (delta_m==-1) && (j%2 == 0);
  if ( ( (m!=0) && (mp!=0) ) && ( sign_rule_1 || sign_rule_2 ) )
    *b *= -1;
  if ( (m==0) || (mp==0) )
    *b *= sqrt(2.0);
  return 1;
}
#line 1345 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__host__ __device__ void zernike_surface_and_partials(rtd *ZS, rtd *dZSdx, rtd *dZSdy,
						     rtd r, rtd o, int max_n,
						     rtd *a, rtd *ax, rtd *ay)
{
  int n, m, mH, a_idx,j;
  rtd R, R2, R4, A, Ax, Ay, S_n, Sx_n, Sy_n, S, Sx, Sy,
    H1, H2, H3, sm, cm, s2, c2, s ,c, sqrt2;
  if (max_n==0) {//( (r>1.0) || (max_n<0) )
    S = a[0];
    Sx = Sy = 0.0;
  } else {
    R = 1.0;
    S =  a[0]*R;
    Sx = ax[0]*R;
    Sy = ay[0]*R;
    
#line 1489 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
sqrt2 = sqrt(2.0);
sincos(2*o,&s2,&c2);
for (n=1;n<=max_n-1;n++) {
  j = (n+1)*(n+2)/2;
  a_idx = j - 1;

  m = n;
  R = pow(r,n);
  sincos(m*o, &sm, &cm);
  if (j%2) {
    Ax = ax[a_idx]*sm;
    Ay = ay[a_idx]*sm;
    A  = a[a_idx--]*sm;
    Ax += ax[a_idx]*cm;
    Ay += ay[a_idx]*cm;
    A  += a[a_idx--]*cm;
  } else {
    Ax = ax[a_idx]*cm;
    Ay = ay[a_idx]*cm;
    A  = a[a_idx--]*cm;
    Ax += ax[a_idx]*sm;
    Ay += ay[a_idx]*sm;
    A  += a[a_idx--]*sm;
  }
  j -= 2;
  Ax *= sqrt2;
  Ay *= sqrt2;
  A *= sqrt2;
  Sx_n = R*Ax;
  Sy_n = R*Ay;
  S_n  = R*A;

  m -= 2;
  if (m>=0) {
    R2 = n*R - (n-1.0)*pow(r,m);
    
#line 1603 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (m>0) {
    c = cm*c2 + sm*s2;
    s = sm*c2 - cm*s2;
    if (j%2) {
      Ax = ax[a_idx]*s;
      Ay = ay[a_idx]*s;
      A  = a[a_idx--]*s;
      Ax += ax[a_idx]*c;
      Ay += ay[a_idx]*c;
      A  += a[a_idx--]*c;
    } else {
      Ax = ax[a_idx]*c;
      Ay = ay[a_idx]*c;
      A  = a[a_idx--]*c;
      Ax += ax[a_idx]*s;
      Ay += ay[a_idx]*s;
      A  += a[a_idx--]*s;
    }
    j -= 2;
    Ax *= sqrt2;
    Ay *= sqrt2;
    A  *= sqrt2;
    cm = c;
    sm = s;
  } else {
    Ax = ax[a_idx];
    Ay = ay[a_idx];
    A  = a[a_idx--];
    --j;
  }
#line 1525 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    Sx_n += R2*Ax;
    Sy_n += R2*Ay;
    S_n  += R2*A;
  }

  m -= 2;
  while (m>=0) {
    mH = m + 4;
    if (r>0) {
      
#line 1485 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H3 = -4.0*(mH-2.0)*(mH-3.0)/((n+mH-2.0)*(n-mH+4.0));
#line 1535 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      
#line 1477 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H2 = H3*0.25*(n+mH)*(n-mH+2.0)/(mH-1.0) + mH - 2.0;
#line 1536 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      
#line 1470 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H1 = mH*(mH-1.0)*0.5 - mH*H2 + H3*(n+mH+2.0)*(n-mH)*0.125;
#line 1537 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      
#line 1451 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
R4 = H1*R + ( H2 + H3/(r*r) )*R2;
#line 1538 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      
#line 1603 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (m>0) {
    c = cm*c2 + sm*s2;
    s = sm*c2 - cm*s2;
    if (j%2) {
      Ax = ax[a_idx]*s;
      Ay = ay[a_idx]*s;
      A  = a[a_idx--]*s;
      Ax += ax[a_idx]*c;
      Ay += ay[a_idx]*c;
      A  += a[a_idx--]*c;
    } else {
      Ax = ax[a_idx]*c;
      Ay = ay[a_idx]*c;
      A  = a[a_idx--]*c;
      Ax += ax[a_idx]*s;
      Ay += ay[a_idx]*s;
      A  += a[a_idx--]*s;
    }
    j -= 2;
    Ax *= sqrt2;
    Ay *= sqrt2;
    A  *= sqrt2;
    cm = c;
    sm = s;
  } else {
    Ax = ax[a_idx];
    Ay = ay[a_idx];
    A  = a[a_idx--];
    --j;
  }
#line 1539 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    } else {
      R4 = (m==0) ? pow(-1.0,0.5*n) : 0.0;
      Ax = 1.0;
      Ay = 1.0;
      A  = 1.0;
    }
    Sx_n += R4*Ax;
    Sy_n += R4*Ay;
    S_n  += R4*A;
    R  = R2;
    R2 = R4;
    m -= 2;
  }
  Sx += sqrt(n+1.0)*Sx_n;
  Sy += sqrt(n+1.0)*Sy_n;
  S  += sqrt(n+1.0)*S_n;
}

n = max_n;
j = (n+1)*(n+2)/2;
a_idx = j - 1;

m = n;
R = pow(r,n);
sincos(m*o, &sm, &cm);
if (j%2) {
  A = a[a_idx--]*sm;
  A += a[a_idx--]*cm;
 } else {
  A = a[a_idx--]*cm;
  A += a[a_idx--]*sm;
 }
j -= 2;
A *= sqrt2;
S_n = R*A;

m -= 2;
if (m>=0) {
  R2 = n*R - (n-1.0)*pow(r,m);
  
#line 1426 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (m>0) {
    c = cm*c2 + sm*s2;
    s = sm*c2 - cm*s2;
    if (j%2) {
      A = a[a_idx--]*s;
      A += a[a_idx--]*c;
    } else {
      A = a[a_idx--]*c;
      A += a[a_idx--]*s;
    }
    j -= 2;
    A *= sqrt2;
    cm = c;
    sm = s;
  } else {
    A = a[a_idx--];
    --j;
  }
#line 1579 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
      S_n += R2*A;
 }

m -= 2;
while (m>=0) {
  mH = m + 4;
  if (r>0) {
    
#line 1485 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H3 = -4.0*(mH-2.0)*(mH-3.0)/((n+mH-2.0)*(n-mH+4.0));
#line 1587 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    
#line 1477 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H2 = H3*0.25*(n+mH)*(n-mH+2.0)/(mH-1.0) + mH - 2.0;
#line 1588 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    
#line 1470 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
H1 = mH*(mH-1.0)*0.5 - mH*H2 + H3*(n+mH+2.0)*(n-mH)*0.125;
#line 1589 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    
#line 1451 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
R4 = H1*R + ( H2 + H3/(r*r) )*R2;
#line 1590 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    
#line 1426 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  if (m>0) {
    c = cm*c2 + sm*s2;
    s = sm*c2 - cm*s2;
    if (j%2) {
      A = a[a_idx--]*s;
      A += a[a_idx--]*c;
    } else {
      A = a[a_idx--]*c;
      A += a[a_idx--]*s;
    }
    j -= 2;
    A *= sqrt2;
    cm = c;
    sm = s;
  } else {
    A = a[a_idx--];
    --j;
  }
#line 1591 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
  } else {
    R4 = (m==0) ? pow(-1.0,0.5*n) : 0.0;
    A = 1.0;
  }
  S_n += R4*A;
  R  = R2;
  R2 = R4;
  m -= 2;
 }
S += sqrt(n+1.0)*S_n;
#line 1361 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
    }
  ZS[0]    = S;
  dZSdx[0] = Sx;
  dZSdy[0] = Sy;
}

#line 1916 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
__global__ void vignetting_kernel(char *mask, ray *d__ray, int N_RAY,
				  rtd inner2, rtd outer2,
				  rtd *d__R, vector *d__origin)
{
  int i, j, ij, iCoordSys, idx;
  rtd rho2;
  rtd x, y, z, u, v, w, x1, y1, s0, k, l, m;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = threadIdx.y;
  iCoordSys = blockIdx.y;
  ij = j * gridDim.x * blockDim.x + i;
  if ( ij<N_RAY )
  {

    idx = iCoordSys;

    u = d__ray[ij].coordinates.x - d__origin[idx].x;
    v = d__ray[ij].coordinates.y - d__origin[idx].y;
    w = d__ray[ij].coordinates.z - d__origin[idx].z;

    idx = iCoordSys*9;

    x = d__R[0+idx]*u + d__R[3+idx]*v + d__R[6+idx]*w;
    y = d__R[1+idx]*u + d__R[4+idx]*v + d__R[7+idx]*w;
    z = d__R[2+idx]*u + d__R[5+idx]*v + d__R[8+idx]*w;

    u = d__ray[ij].directions.x;
    v = d__ray[ij].directions.y;
    w = d__ray[ij].directions.z;

    k = d__R[0+idx]*u + d__R[3+idx]*v + d__R[6+idx]*w;
    l = d__R[1+idx]*u + d__R[4+idx]*v + d__R[7+idx]*w;
    m = d__R[2+idx]*u + d__R[5+idx]*v + d__R[8+idx]*w;

    if (m==0) { return; }
    s0 = -z/m;
    x1 = x + k*s0;
    y1 = y + l*s0;
    rho2 = x1*x1 + y1*y1;
    if ( (iCoordSys==0) && (rho2<inner2) )
      {
	return;
      }
    if (rho2<=outer2)
      {
	mask[ij] = 1;
      }
  }
}
#line 1845 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void aperture::setup(rtd _D_, int _D_px_,
		     vector _origin_, vector _euler_angles_)
{
  D_px = _D_px_;
  D_seg = D = _D_;
  ri   = 0.0;
  N = 1;
  V.setup(D_px*D_px);
  ref_frame.setup(_origin_, _euler_angles_);
}
#line 1858 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void aperture::setup(rtd _D_, rtd _ri_, int _D_px_,
		     vector _origin_, vector _euler_angles_)
{
  D_px = _D_px_;
  D_seg = D = _D_;
  ri   = _ri_;
  N = 1;
  V.setup(D_px*D_px);
  ref_frame.setup(_origin_, _euler_angles_);
}
#line 1872 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void aperture::setup(rtd _D_, rtd _ri_, rtd _D_seg_, int _D_px_,
		     vector *_origin_, vector *_euler_angles_, int N_segment)
{
  D_px  = _D_px_;
  D     = _D_;
  D_seg = _D_seg_;
  ri    = _ri_;
  N = N_segment;
  V.setup(D_px*D_px*N);
  ref_frame.setup(_origin_, _euler_angles_, N);
}
#line 1889 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void aperture::cleanup(void)
{
  INFO("@(CEO)>aperture: freeing memory!\n");
  V.cleanup();
  ref_frame.cleanup();
}
#line 1899 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void aperture::vignetting(bundle *rays)
{
  rtd R2, Rri2;
  R2 = D_seg*D_seg*0.25;
  Rri2 = R2*ri*ri;
  //INFO("R2=%g - Rri2=%5.2f\n",R2,Rri2);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  HANDLE_ERROR( cudaMemset(V.m, 0, sizeof(char)*rays->N_RAY ) );
  vignetting_kernel <<< gridDim , blockDim >>> (V.m, rays->d__ray, rays->N_RAY,
						Rri2, R2,
						ref_frame.d__R,
						ref_frame.d__origin);
  intersection <<< gridDim , blockDim >>> (V.m, rays->d__ray, rays->N_RAY);
}

#line 2029 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void transform_to_S(bundle *rays, conic *F)
{
  transform_to_S_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1,
				   F->ref_frame.d__R, F->ref_frame.d__origin);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  transform_to_S_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
					     F->ref_frame.d__R, F->ref_frame.d__origin);
}
void transform_to_S(bundle *rays, aperture *A)
{
  transform_to_S_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1,
				   A->ref_frame.d__R, A->ref_frame.d__origin);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  transform_to_S_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
					     A->ref_frame.d__R, A->ref_frame.d__origin);
}
#line 2092 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void transform_to_R(bundle *rays, conic *F)
{
  transform_to_R_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1,
				   F->ref_frame.d__R, F->ref_frame.d__origin);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  transform_to_R_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
        	 				     F->ref_frame.d__R, F->ref_frame.d__origin);
}
void transform_to_R(bundle *rays, aperture *A)
{
  transform_to_R_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1,
				   A->ref_frame.d__R, A->ref_frame.d__origin);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  transform_to_R_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
        	 				     A->ref_frame.d__R, A->ref_frame.d__origin);
}
#line 2166 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void intersect(bundle *rays, conic *F)
{
  intersect_chief_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1,
					F->k, F->c, F->d__origin,
					F->even_asphere_N, F->d__even_asphere_a);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  intersect_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY,
					       F->k, F->c, F->d__origin,
					       rays->d__chief_ray,
					       F->even_asphere_N, F->d__even_asphere_a);
}

#line 2469 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void reflect(bundle *rays)
{
  reflect_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1, 1.0);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  reflect_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY, 1.0);
}
#line 2418 "/home/ubuntu/projects/crseo/sys/CEO/rayTracing/rayTracing.nw"
void refract(bundle *rays, const rtd mu)
{
  refract_kernel <<< 1 , 1 >>> (rays->d__chief_ray, 1, mu);
  dim3 blockDim(N_THREAD,N_THREAD);
  dim3 gridDim(rays->N_RAY/N_THREAD2+1,1);
  refract_kernel <<< gridDim , blockDim >>> (rays->d__ray, rays->N_RAY, mu);
}