Dear Carlos,
the scale_factor variable and the unit conversion functions were only implemented, to allow inexperienced users covert easily the energy unit of the source file.
Since you are writing your own code, and your energies are already in GeV, there is no need to add this feature.
I also made same changes to your routine: muonsatsurface.f (6.1 KB)
The biggest one is since you want to return two values, I converted the function to a subroutine
I removed the recursive call. I was there only to protect against returning non positive energy.
There is no problem using the probabilities from the file, but you need to keep in mind, that the sum of them has to be 1.0. You could implement a check, or add a normalization step.
Cheers,
David