I'm asking you for a help to solve a difficult (for myself) task. Now I will tell you what and why I'm trying to implement.

The task is to simulate the passage of elementary particles through a material consisting of nano-objects. To solve the task, you need to generate many random variables - the angles, on which an elementary particle will be scattered after the interaction with nano-objects. The random variables are distributed according to the law - the differential cross section of scattering, which has a rather complicated form.

In the beginning sets the scattering amplitude $fB[theta, v, npr]$, where $theta$ is the scattering angle, $v$ is the velocity of an elementary particle, $npr$ is the radius of the nanoparticle:

(* h - Planck constant *)
(* m - the mass of elementary particle *)
k[V_] := m V/h;
q[theta_, v_] := 2 k[v] Sin[theta/2];
fBintr[theta_, v_, npr_] := -((2 m)/h^2) Uvnyt (Sin[npr q[theta, v]] -
    npr Cos[npr q[theta, v]] q[theta, v])/q[theta, v]^3;
fBextr[theta_, v_, npr_] := -((2 m)/h^2) Uvnesh (Sin[npr q[theta, v]] -
    npr Cos[npr q[theta, v]] q[theta, v])/q[theta, v]^3;
fB[theta_, v_, npr_] := fBintr[theta, v, 0.75 npr] + fBextr[theta, v, npr];

As you can see, $fB[theta, v, npr]$ consists of two parts - internal and external, which describe a two-stage structure optical potential of the nano-object. Thus the nanoparticle consists of a solid core and a shell with different chemical composition and different optical interaction potential. The differential scattering cross-section is a square of the amplitude:

dsigma[theta_, v_, npr_] := Abs[fB[theta, v, npr]]^2;

Using the characteristic values for the task - $theta$ angle ranges [0, Pi], the velocity $v$ = 100 m/s, the radius of the nanoparticles $npr$ = 2.5 nm, we obtain the following graph: The differential scattering cross-section
And all looks good, but actually it is not. The differential scattering cross section can also be expressed using spherical Bessel functions of the first kind SphericalBesselJ that have instability at zero. It looks like this:
The instability at zero
Further, if we normalize this distribution to use as PDF $dsigma[t, 100, 2.5*10^{-9}]/a$:

mydist[v_, npr_] := With[{a = NIntegrate[dsigma[tt, v, npr], {tt, 0, Pi}]},
    ProbabilityDistribution[{"PDF", dsigma[t, v, npr]/a}, {t, 0, Pi}, 
    Assumptions -> {v > 0, npr > 0}]
RandomVariate[mydist[100, 2.5*10^-9], {10}]

Then RandomVariate it just refuses to work and to generate random angles. I think that this is due to the instability near zero, and because $dsigma$ has a large number of $theta$ cycles. So Mathematica is unable to verify that I specified a valid PDF because it loses too much precision in the integration and the integral $dsigma[theta, 100, 2.5*10^{-9}]/a$ will not be equal to 1.

I have tried several options for generating random angles.

  1. Using the RandomReal I generated uniformly distributed random value between the lower boundary $dsigma[theta, 100, 2.5 * 10 ^{-9}] / a$ (non-zero) to its upper boundary Limit[dsigma[theta, 100, 2.5*10^(-9)]/a, theta->0]. Then using FindRoot was looking for a matching $theta$. It works. But the problem is that such angles during the solution must generate millions, if not more. Accordingly, the method is very time-consuming in comparison with the generation of a large number of random variables using RandomVariate.
  2. Had a look at the nice tutorial how to Create Your Own Distribution with Oleksandr Pavlyk (and actually more confused). Learned another way to generate.

    Random`RandomRealRejection[ densityFunc, { hatFunc, hatRandomVector }, len, prec ] generates random variates from continuous univariate density. The density need not be normalized, as long as densityFunc <= hatFunc, and hatRandomVector draws from distribution with density ~ to hatFunc.

    a = NIntegrate[dsigma[t, 100, 2.5*10^(-9)], {t, 0, Pi},
        WorkingPrecision -> MachinePrecision];
    b = 0.01 + Limit[dsigma[x, 100, 2.5*10^(-9)], x -> 0]/a;
    c = dsigma[Pi, 100, 2.5*10^(-9)]/a;
    (* I am using piecewise *)
            Piecewise[{{Limit[dsigma[x, 100, 2.5*10^(-9)], x -> 0]/a, x == 0},
            {dsigma[x, 100, 2.5*10^(-9)]/a, x > 0}}], Listable],
        {b &, Function[{len, prec}, RandomReal[Pi, len, WorkingPrecision -> prec]]},
        10, MachinePrecision
    Show[Histogram[%, Automatic, "PDF"], 
        Plot[Piecewise[{{Limit[dsigma[x, 100, 2.5*10^(-9)], x -> 0]/a, x == 0},
        {dsigma[x, 100, 2.5*10^(-9)]/a, x > 0}}], {x, 0, Pi}, PlotStyle -> Red,
        PlotRange -> All], Plot[b, {x, 0, Pi}, PlotStyle -> Red, PlotRange -> All]

    And in the end successfully generated random angles: RandomRealRejection method
    The problem is that, as I said earlier, these angles need to have tens of millions, and the generation of only $10^5$ was approximately 8 minutes (+/- 1 minute). Thus, this option is not very suitable for a large number of calculations.

The question I would like to ask is how to specify the probability distribution in my case? I suspect that it is not enough to use the correct method of random variables generating, but it is necessary to find a suitable algorithm for the distribution assignment.

Thank you in advance for any advice and assistance!

closed as off-topic by MarcoB, corey979, m_goldberg, Kuba Feb 24 '17 at 7:04

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "The question is out of scope for this site. The answer to this question requires either advice from Wolfram support or the services of a professional consultant." – MarcoB, corey979, m_goldberg, Kuba
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 3
    Too much information. I would suggest honing your question down to the essence of the Mma issue ... in a couple of lines – wolfies Apr 25 '16 at 11:08
  • @wolfies, ok. So, how to generate a RandomVariate of a complex custom distribution? – AlexPlus Apr 25 '16 at 11:37
  • The starting point to your question is to define the distribution you are working with, in the simplest possible way. You just have a whole bunch of functions which depend on some unspecified unnamed parameters called Uvnyt and Uvnesh. You then NIntegrate the latter, which cannot possibly work since Uvnyt and Uvnesh are unstated symbolic constants, and not numerical. I cannot see how RandomVariate could possibly be expected to work on the latter: it requires that all parameters be given numerical values. – wolfies Apr 25 '16 at 14:42
  • @wolfies, thank you for your comment. You are right about Uvnyt and Uvnesh - they are not defined here, I forgot about it. But they just have the constant numerical values. So there is no problem with NIntegrate and RandomVariate in that case. – AlexPlus Apr 25 '16 at 17:15
  • I feel like your question is more about math.physics, than it's implementation. Bessel functions are not defined in zero. They are also not normalizable in regular sense. As far as I understand, which I am not sure, as I have not seen your derivations, you are trying to extract PDF's in a regular way. Well, I don't think it would work even analytically. Check the properties of Bessel functions. – MsTais Apr 25 '16 at 17:58