static double IMP_PROTEIN_DENSITY=0.826446; // 0.826446=1/1.21   don't know exactly where that constant comes from
static double protein_density_ = IMP_PROTEIN_DENSITY;
double get_used_protein_density(){
  return protein_density_;
 };

double set_used_protein_density(double densityValue){
  IMP_USAGE_CHECK(densityValue > 0,
        "A protein density value is most probably a positive value");
  protein_density_=densityValue;
  return protein_density_;
 }

double set_used_protein_density(ProteinDensityReference densityReference){
  double density=0.0;
  switch (densityReference) {
    case HARPAZ:
      density = IMP_PROTEIN_DENSITY;
      break;
    case ANDERSSON :
      density = 0.7347;  // Andersson and Hovmöller (1998) Theoretical 1.22 g/cm3
      break;
    case TSAI:
      density = 0.84309; // Tsai et al. (1999) Theoretical 1.40 g/cm3
      break;
    case QUILLIN:
      density = 0.86116; // Quillin and Matthews (2000) Theoretical 1.43 g/cm3
      break;
    case SQUIRE:
      density = 0.82503; // Squire and Himmel (1979) and Gekko and Noguchi (1979) Experimental 1.37
      break;
    default :
      // should warn IMP_WARNING("unknown density reference... density set to default");
      density = IMP_PROTEIN_DENSITY;
   }
  density=set_used_protein_density(density);
  return density;
 }










Float compute_threshold_for_approximate_volume(DensityMap* d, Float volumeDesired) {
  Float voxelVolume     = d->get_spacing()*d->get_spacing()*d->get_spacing();
  long numVoxelsNeeded  = volumeDesired / voxelVolume;
  long mapSizeInVoxels  = d->get_number_of_voxels();
  std::vector<emreal> data(mapSizeInVoxels);
  for (long l=0;l<mapSizeInVoxels;++l) {
    // This ugly loop should be replaced by a memcopy if only I had access to d->data_
    data[l]=(d->get_value(l));
   }
  std::sort(data.begin(),data.end());
  emreal threshold = data[mapSizeInVoxels-numVoxelsNeeded];
  return static_cast<Float>(threshold);
 }

Float compute_threshold_for_approximate_mass(DensityMap* d, Float desiredMass){
  Float desiredVolume = desiredMass / get_used_protein_density();
  return compute_threshold_for_approximate_volume(d,desiredVolume);
}

Float approximate_molecular_mass(DensityMap* d, Float threshold) {
  long counter=0;//number of voxels above the threshold
  for(long l=0;l<d->get_number_of_voxels();l++) {
    if (d->get_value(l) > threshold) {
      ++counter;
     }
   }
  Float s=d->get_spacing();
  return  s*s*s*counter*get_used_protein_density();
 }