A better way to create a response function for CSD

This implementation includes the recursive calibration of the response function for CSD, Using an FA threshold is not a very robust method.It is dependent on the dataset (non-informed used subjectivity), and still depends on the diffusion tensor (FA and first eigenvector), which has low accuracy at high b-value. This function recursively calibrates the response function, for more information see Tax et al. NeuroImage 2014.

This implementation is slow only on one line (788) which calls real_sph_harm. I would appreciate any help to optimize this function. Thx!
Thanks @pv for your feedback. I added these lines according to the function real_sph_harm in dipy, which uses different convention from scipy:

```
def real_sph_harm(m, n, theta, phi):
    # dipy uses a convention for theta and phi that is reversed with respect to
    # function signature of scipy.special.sph_harm
    sh = sph_harm(np.abs(m), n, phi, theta)
    real_sh = np.where(m > 0, sh.imag, sh.real)
    real_sh *= np.where(m == 0, 1., np.sqrt(2))
    return real_sh
```

When I use only the sph_harm function from scipy and yours I get two different results: I tried with ```m = np.array([ 0 -2 -1  0  1  2 -4 -3 -2 -1  0  1  2  3  4 -6 -5 -4 -3 -2 -1  0  1  2  3
  4  5  6 -8 -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8])```, ```n = np.array([0 2 2 2 2 2 4 4 4 4 4 4 4 4 4 6 6 6 6 6 6 6 6 6 6 6 6 6 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8])```, and lists of theta and phi
```
theta = np.array([ 3.08471252  1.5187582   1.57217178  2.14791083  1.29522236  0.42860416
  2.49054837  1.29093715  0.35181337  2.05196433  1.88205117  2.0172383
  1.02859042  2.19811741  2.17987986  0.54873935  0.45649432  1.84598342
  1.63990503  2.45273092  0.98544454  2.45363653  2.77217512  1.30317169
  1.58706864  1.73514323  2.39954954  1.72457508  1.67126966  0.43811831
  2.14015788  1.1752897   1.87887165  1.21314693  0.22387675  2.8249272
  1.98664626  2.40358562  0.76385825  0.79364915  0.77623336  1.38302801
  2.38465215  1.63795237  2.14682278  1.45391706  0.93587114  2.52734611
  1.89946007  1.72251569  0.71732319  1.07701534  1.01218582  1.38454091
  1.07637993  1.91337867  0.74734146  1.25701982  1.55382064  1.55546254
  2.17049165  1.8352862   0.72752261  1.25029512])
```
```
phi = np.array([ 3.12102262 -2.84760175 -1.98409072  3.13050127  2.43523068  2.42322314
  1.73045006 -2.54577098 -1.55268133 -2.65769859 -2.4172876  -2.08889762
  2.31478875 -2.33433495  1.78212456  3.05074976  1.68468091 -1.89277974
 -2.60209419  2.65297338 -2.63136679  2.18608633  2.78594995 -2.20327949
  1.68485126 -3.05338411 -3.1373576  -1.6529307   2.71042227 -2.35366419
  2.19039837  1.69626492  1.6792922   1.99498978 -3.00995926 -2.24411209
 -2.92893596 -2.7347392  -1.84608784  1.79238318  2.66000342  2.69323866
 -1.94707106 -2.25352281  2.71805959  2.21382628  3.12755056 -2.33513684
  2.53355843  2.29259891 -2.85015235  2.72232179 -2.17820482 -1.8102771
 -1.74767739  2.94928089  2.19218307  3.03103661  1.95027264  3.00054962
 -1.6910932   2.03966614 -2.312912   -2.93454625])
```
Does this give different results for you too?