 __device__ void Spin(PhotonStruct*, float,unsigned long long*,unsigned int*);
 __device__ unsigned int Reflect(PhotonStruct*, int, unsigned long long*, unsigned int*, int);
 __device__ unsigned int PhotonSurvive(PhotonStruct*, unsigned long long*, unsigned int*);
-
+__device__ unsigned int MoveToFirstBoundary(PhotonStruct*, unsigned short, short*);
 
 __global__ void MCd(MemStruct DeviceMem)
 {

	for(;ii<NUMSTEPS_GPU;ii++) {
    // Main while loop

  	if(bulks_dc[p.bulkpos].mutr!=FLT_MAX)
			//s = -__logf(rand_MWC_oc(&x,&a))*bulks_dc[p.bulkpos].mutr;//sample step length [cm] //HERE AN OPEN_OPEN FUNCTION WOULD BE APPRECIATED
      s = bulks_dc[p.bulkpos].mutr;
    else
			//s = 100.0f; //temporary, say the step in glass is 100 cm.
      s = 2.0f;

		//Check for layer transitions and in case, calculate s
		new_bulk = p.bulkpos;

    //int side_scape = 0;

    //Check for upwards reflection/transmission and move to surface
		if(p.z+s*p.dz<0.) {
      new_bulk = 0;
      //side_scape = 0;
      s = __fdividef(-p.z,p.dz);
    }

    //Check for downward reflection/transmission and move to surface
		if(p.z+s*p.dz>(*esp_dc)){
      new_bulk = last_bulk;
      //side_scape = 0;
      s = __fdividef((*esp_dc)-p.z,p.dz);
    }

    // Move photon TODO:here?
    p.x += p.dx*s;
    p.y += p.dy*s;
    p.z += p.dz*s;

    //if(p.z>(*esp_dc)) p.z=(*esp_dc);//needed? TODO
		//if(p.z<0.) p.z=0.;//needed? TODO

    // Retrieve bulk position
    if(new_bulk!=0 && new_bulk!=last_bulk) {
      if (fabsf(p.x)<2*(*esp_dc) && fabsf(p.y)<2*(*esp_dc) && p.z<(*esp_dc)){
        // Inside space of 3D matrix. Calculate index of the current voxel
        // Use round to zero so there are no over sampled voxels (for ex: (max_x,0,0) and (0,1,0) should not map to the same voxel)
        index = __float2uint_rz((p.x+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
              + num_x * (__float2uint_rz((p.y+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
              + num_y * __float2uint_rz((p.z)*__int2float_rn(*grid_size_dc)));
        new_bulk = DeviceMem.bulk_info[index];
      }
      else {
        // Photon scaped to the sides
        // Outside space of 3D matrix, assume inside homogeneous medium (should we assume it is outside the bulk (bulkpos=0)? TODO)
        new_bulk = 1;
        }
    }

    unsigned int reflected;


	float s;	// step length
	unsigned int index; // temporal variable to store indexes to arrays
	unsigned int w; // photon weight
  //bool in_glass = FALSE;

	PhotonStruct p = DeviceMem.p[begin+tx];

