I got an algorithm which works well in regard of what I need. Now, I should want to get an analytic form too.
$$(x,y) \in [ {-\pi \over 2} , {\pi \over 2} ]^2 , h \in [ {0} , {\pi \over 8} ] , c \in [ {0 } , {\pi } ] ,
n \in N_+^* , m \in N_+^* , m \simeq 3 n$$ . n and m are static parameters.
The Monte-Carlo algorithm uses 3 classical states f, g and z which have probability functions depending on x and h: $$k(x,h) = { { ( 1 - ( 1 - Cos^{2n}( x -h )-Sin^{2n}( x -h ) )^m) } \over { (Cos^{2n}( x -h )+Sin^{2n}( x -h ))} } \\
f(x,h) = {k(x,h) \quad Cos^{2n}( x -h ) } \\
g(x,h) = {k(x,h) \quad Sin^{2n}( x -h ) } \\
z(x,h) = { 1 - f(x,h) - g(x,h) }$$
I built the distribution with a markov chain describing a photon in some perfect crystal with a little dissipation. m and n acting as dissipation parameters.
It takes 3 random values x,h and y, assumes that c=y-x and computes the function P(x,h,c) which is $$ = { f(x,h) f(x+c,h) + g(x,h) g(x+c,h) }$$ and constructs p(c) as a probability of c which is independent of x and h. For symmetry reasons, c>0.
$$p(c) ={ {8 \over {\pi ( \pi -c ) }} { \int_{-\pi \over 2}^{{\pi \over 2}-c} { \int_{0}^{\pi \over 8}
{P(x,h,c) } \quad dh } \quad dx } } \\ \\$$
Question #1 : Is it the correct construction in regard of the above distribution ? I assumed that the order doesn't matter when computing p(c), this allowing to sum first on the common h and then on x. But I walk on eggs ! y the other random value is absent. Then I bind the 1st integral to get only allowed values. For example, if $$c=\pi/2$$ and $$x = \pi/4$$, then $$y=x+c > \pi/2$$ which is is not allowed.
Question #2 :How to compute the integral ? I always failed.
Facing difficulties to get pretty and easy expressions, I abandoned the previous f,g,z distribution to build handly another distribution less complicated. There is no more the parameter m , other variables are the same. The results fit less my hopes but it is enough.
the new probability distribution is:
$$f(x,h) = {{Cos^{2n+2}( x -h )} \over { Cos^{2n}( x -h )+Sin^{2n}( x -h ) } }\\
g(x,h) = {{Sin^{2n+2}( x -h )} \over { Cos^{2n}( x -h )+Sin^{2n}( x -h ) } } \\
z(x,h) = { 1 - f(x,h) - g(x,h) }$$
This time $$P(x,h,c) = { f(x,h) f(x+c,h) + g(x,h) g(x+c,h) } \\
={ { Cos^{2n+2}( x -h ) Cos^{2n+2}( x+c -h ) + Sin^{2n+2}( x -h ) Sin^{2n+2}( x+c -h ) }
\over { (Cos^{2n}( x+c -h )+Sin^{2n}( x+c -h ) ) (Cos^{2n}( x -h )+Sin^{2n}( x -h ) )} }$$
And p:
$$p(c) ={ {8 \over {\pi ( \pi -c ) }}
{ \int_{-\pi \over 2}^{{\pi \over 2}-c}
{ \int_{0}^{\pi \over 8}
{ { Cos^{2n+2}( x -h ) Cos^{2n+2}( x+c -h ) + Sin^{2n+2}( x -h ) Sin^{2n+2}( x+c -h ) }
\over { (Cos^{2n}( x+c -h )+Sin^{2n}( x+c -h ) ) (Cos^{2n}( x -h )+Sin^{2n}( x -h ) )} } dh
} \quad dx
} }$$
I found painful methods until n=3 giving very long results but I failed to do better.
The simulation seems to show that the best curve is got with n=7 or n=9 for the 2 distributions ( n=7, m=22 for the original one).
Question #3 :How to compute the integral ?
Question #4 : to check them in simulations: do you know another probability functions f and g with 1 or 2 arguments, possibly well known, which can replace the above ugly attempts?