   const int num_y = (int)(4 * (simulation->esp) * (float)simulation->grid_size);
   const int num_z = (int)((simulation->esp) * (float)simulation->grid_size);
   const int fhd_size = num_x * num_y * num_z;
+
+  // Output temporal detectors
+  const int num_x_tdet = simulations[0].det.x_temp_numdets;
+  const int num_y_tdet = simulations[0].det.y_temp_numdets;
+  const long num_tbins = simulations[0].det.temp_bins;
+  const long timegrid_size = num_x_tdet * num_y_tdet * num_tbins;
 
   cudaError_t cudastat;
   clock_t time1, time2;

    tempfhd[xyz] = ((double)HostMem.fhd[xyz]/(0xFFFFFFFF*photons_finished));
  }

  printf ("Photons simulated: %llu\n\n", photons_finished);
  FreeMemStructs(&HostMem, &DeviceMem);
  return photons_finished;


  const int fhd_size = num_x * num_y * num_z; //x + HEIGHT* (y + WIDTH* z)


  // Store in local variables the number of time detectors

  const int num_x_tdet = simulations[0].det.x_temp_numdets;
  const int num_y_tdet = simulations[0].det.y_temp_numdets;
  const long num_tbins = simulations[0].det.temp_bins;

  const long timegrid_size = num_x_tdet * num_y_tdet * num_tbins;

  // FHD simulation
  // Run a simulation

