2 weeks over! That was fast..
And I have sent my pull request finishing the first of 4 sub-parts of my project.
The link to the pull request is https://github.com/scipy/scipy/pull/3692.
Here the function to implement spherical harmonics has been rewritten in Cython. It uses lpmv fortran sub-routine (with the help of the existing wrapper function pmv_wrap) to calculate associated Legendre Function.
Initially the function to calculate spherical harmonics was written in python which called lpmn function to calculate the associated Legendre function. Lpmn function in turn uses lpmn fortran routine to calculate the function. The routine iteratively evaluates associated Legendre function for all orders and degrees m and n. While the routine lpmv avoids the iteration by using recursion.
The usage of a different routine and rewriting the code in cython has given a huge leap in the speed of the function which is shown in a simple benchmark performed by Pauli,
Before:
 In [4]: %timeit sph_harm(50, 50, 0.5, 0.9)
10000 loops, best of 3:
81.3 µs per loop
 In [5]: x = np.linspace(0, 1, 2000)
 In [6]: %timeit sph_harm(50, 50, 0.5, x)
 10 loops, best of 3:
127 ms per loop

After:
 In [2]: %timeit sph_harm(50, 50, 0.5, 0.9)
100000 loops, best of 3:
 4.52 µs per loop
 In [3]: x = np.linspace(0, 1, 2000)
In [4]: %timeit sph_harm(50, 50, 0.5, x)
1000 loops, best of 3: 1.87 ms per loop

Now the next part of the project is to implement the ellipsoidal harmonic functions. It is kind of challenging as they are complex and they havent been implemented much. My next 2-3 weeks are dedicated to the implementation of the both kinds of ellipsoidal functions, thus meeting my mid sem deadline.