The original wave function of Hydrogen:

///////////////////////////////////////////////////////////////////////////////////////////////////////////

           1              m                         -m    m *
Psi      = - * U   (r) * Y (Theta,Phi)        and  Y   = Y      * = imaginary
   n,l,m   r    n,l       l                         l     l
                                          2*r   
          a = 2*l+1  ,  k = n-l-1  ,  t = --- :
                                           n
                                ___________
                    _____      /   1                     |          (    )v |
                   / 2*r      /-----------     -t/2   k  |( k+a )   ( -t )  |
          U   =   /  --- *   /     ( k+a )  * e    * SUM |(     ) * ( -- )  |
           n,l   V    n     / a! * (     )           v=0 |( k-v )   ( v! )  |
                           V       (  k  )


          u = cos(Theta):
                          ________________                           m  |         l  { (       )l } |
           m       m     / 2*l+1   (l-m)!    i*m*Phi     m          d   |  1     d   { (  2    )  } |
          Y  = (-1) *   /  ----- * ------ * e       * sin (Theta) * --- |----- * --- { ( u - 1 )  } |
           l           V   4*pi    (l+m)!                             m | l        l {            } |
                                                                    du  |2 * l!  du  {            } |



//////////////////////////////////////////////////////////////////////////////////////////////////////////

Imaginary part:

    i*m*Phi
   e        = cos(m*Phi) + i*sin(m*Phi)
 
   This was my bottle neck at the start of my coding, for I coundn't find anything about the interpretation of 'i'.
   At last I decided that it must be something like the following and  ... yes it worked!!

   If m > 0 : cos(|m|*Phi) + sin(|m|*Phi)
   If m < 0 : cos(|m|*Phi) - sin(|m|*Phi)

   or:

    i*m*Phi
   e        = cos(|m|*Phi) + m /|m|*sin(|m|*Phi)
   
When POV-Ray got the sum() and prod() functions, I could do the following:

Differentation:
    p
   d   (  n )                                  n-p   p-1  |       |    n-p
   --- ( x  ) = n * (n-1) * ... * {n-(p-1)} * x    = PROD | n - v | * x
     p (    )                                        v=0  |       |
   dx

Series expansion of (u^2 - 1)^l :

        n     n  | ( n )    n-i    i |            (  2    )l    l  | ( l )    2*(l-p)      p |
   (a+b)   = SUM | (   ) * a    * b  |     gives  ( u - 1 )  = SUM | (   ) * u       * (-1)  |
             i=0 | ( i )             |                         p=0 | ( p )                   |


////////////////////////////////////////////////////////////////////////////////////////////////////////////