//----------------------------------------------------------------------------------------------------------------------------------
//	p! = p.(p-1).(p-2) .. .3.2.1

#declare Fac = function(Val) { select(Val, 1, 1, prod(i, 1, Val, i)) } 

//----------------------------------------------------------------------------------------------------------------------------------
/*  ( p )      p!	
    (   )	=	--------
    ( q )   q!(p-q)!
*/

#declare Binom = function(P, Q) { Fac(P)/Fac(Q)/Fac(P - Q) }

//----------------------------------------------------------------------------------------------------------------------------------
/*
       1          l [     p  ( l )   l+m-1 {            }    2(l-p)-(l+m) ]
    ------- . Sigma [ (-1) . (   ) . Tau   { (2(l-p)-q) } . µ             ]		,	with µ = cos(Theta)  and  |m| !!
     l          p=0 [        ( p )   q=0   {            }                 ]
    2 * l!

	Derivative = 1 /2^l /Fac(l) * SUM(p=0 -> l) { (-1)^p*Binom(l, p)* PRODUCT(q=0 -> l+m-1) [2(l-p)-q] * µ^(2(l-p)-(l+m)) }
						 = Constant(l) * SUM(p=0 -> l) { Constant(l,m,p) * µ^(2(l-p)-(l+m)) }
						 = SUM(p=0 -> l) { Constant(l) * Constant(l,m,p) * µ^(2(l-p)-(l+m)) }
						 = SUM(p=0 -> l) { Constant_(l,m,p) * µ^(l-2p-m) }
*/

#macro Product(L, absM, P)

  #local Prod = 1;
  #local Q = 0;
  #while (Q <= L + absM - 1)
    #local Prod = Prod*(2*(L - P) - Q);
    #local Q = Q + 1;
  #end // while      

  Prod

#end // macro Product



#macro DerivateFunction(L, M)

  #if (L = 0)
    #local DerFn = function { 1 }
  #else
    #local absM = abs(M);
    #local ConstL = 1/pow(2, L)/Fac(L);
    #local DerFn = 
      function {
        #local P = 0;
        #while (P <= L)
          #local ConstLMP = ConstL*pow(-1, P)*Binom(L, P)*Product(L, absM, P);
          #local EE = L - 2*P - absM;
          #if (EE >= 0)
            + ConstLMP*pow(x/f_r(x, y, z), EE)
          #end // if
          #local P = P + 1;
        #end // while
      }
  #end // if

  function { DerFn(x, y, z) }
  
#end // macro DerivateFunction

//----------------------------------------------------------------------------------------------------------------------------------
//	Sigma_n_l = SUM(v=0 -> k) {Binom(n+l,k-v)/Fac(v)*(-2r/n)^v}
//						=   2^0/(0!.n^0).Binom(n+l,k-0).r^0 +
//							- 2^1/(1!.n^1).Binom(n+l,k-1).r^1 +
//							+ 2^2/(2!.n^2).Binom(n+l,k-2).r^2 +
//							- 2^3/(3!.n^3).Binom(n+l,k-3).r^3 + .....

#macro SigmaFunction(N, L)

  #local K = N - L - 1;
  #local SumFn =
    function(x, y, z, R) {
      #local V = 0;
      #while (V <= K)
        #local CS = Binom(N + L, K - V)/Fac(V)*pow(-2/N, V);
        + CS*pow(R, V)
        #local V = V + 1;
      #end // while
    }

  function { SumFn(x, y, z, f_r(x, y, z)) }

#end // macro SigmaFunction

//----------------------------------------------------------------------------------------------------------------------------------

#macro Exp_imPhiFunction(M)

  #local absM = abs(M);

  function {
    #if (absM <= 0)
      1
    #else
      cos(absM*atan2(y, z)) + sin(absM*atan2(y, z))
    #end // if
  }

#end // macro Exp_imPhiFunction

//----------------------------------------------------------------------------------------------------------------------------------

#declare Derivative = DerivateFunction(l, m)
    
#declare Sigma_n_l = SigmaFunction(n, l)

#declare Exp_imPhi = Exp_imPhiFunction(m)


#declare abs_m = abs(m);

#declare C_U = pow(2, l + 1) / pow(n, l + 2) / sqrt( Fac(2*l + 1)*Binom(n + l, n - l - 1) );

#declare C_Y = pow(-1, abs_m) * sqrt( (2*l + 1)/4/pi * Fac(l - abs_m)/Fac(l + abs_m) );


#declare Psi = 
	function	{
		select(f_r(x, y, z), 0, 0, // To prevent division by zero errors
			 C_U*pow(f_r(x, y, z), l) / exp(f_r(x, y, z)/n) * Sigma_n_l(x,y,z) 
			*C_Y*Exp_imPhi(x, y, z) * pow(f_r(0, y, z)/f_r(x, y, z), abs_m) * Derivative(x, y, z)
		)
	}

#if (Prob)
	#declare Func = function { 4*pi*(x*x + y*y + z*z)*pow(Psi(x, y, z), 2) }
#else
	#declare Func = function { Psi(x, y, z) + 0.5 }
#end	//

//	--------------------------------------------------------------------------------------------------------------------