No hope for Mathematica

There is no hope as such for Mathematica problem. I would have to reply solely on IDL and probable integration of which with Mathematica. I still have one more option open, to write iteration method inside Mathematica. One such fine example can be found here:
pi[n_, theta_] := 
Module[{p1,p2}, 
   p1=0; 
   p2=1; 
   Do[ 
     p1=p2*(2i-1)/(i-1) Cos[theta] -  i/(i-1) p1; 
     tmp=p2; 
     p2=p1; 
     p1=tmp,{i,2,n}] 
   p2 
]

But it doesn't make a sense to try too much with Mathematica, when IDL is working pretty good.

Problem with Mathematica

The problem with Mathematica intensified. It seems that recursion is not working properly with Mathematica. The question has been submitted in the forum. It should be seen, what others understand for it.
But IDL is working OK. I have found the error with the understanding that if “nang” has been taken for granted for “q” wasn't a good idea, since algorithm relies heavily on value of “cos(theta)” which has reflective values on either side of 90 degrees. So now the conversion from q to angle has been introduced. Product of which can be seen next time. No intense involvment has been on my side due to time constraints.

Major IDL Changes

IDL programme has be reedited. All the mistakes relating to misunderstanding of role of “j” and “jj” has been cleared. Now the next step would be to try to integrate the relation q=2k sin(theta/2) must be introduced. All the values related to “nang” must be taken care after. Optimism is there to be able to implement it with “q” vector. The attempt with Mathematica had not been so successful yet. Trouble remains with exact execution of recursive function.

It is written now as,
pi[i_,theta_] := (2i-1)/(i-1) Cos[theta] pi [i-1, theta] - n/(n-1) pi[i-2, theta]
tau[i_,theta_] := I Cos[theta] pi[i,theta] - (n+1) pi[i-1,theta]

The problem still lies with initialization pi[0,theta] = 0, pi[1,theta] = 1 for all 'theta'.

Then with writing S1 & S2 function in loop, which are, in fact, Sum over 1-nstop, probably can be written like,
S1Temp[l_,theta_] := (2l+1)/(l(l+1)) a[l] pi[l,theta] + b[l] tau[l,theta]
then
S1[theta_] = Sum[S1Temp[l,theta], {l,1,nstop}]

New plan

Preparation:
Read given 2 papers & Jens-Nielson for units-problem
Jens-Nielson & Klaus's Arbeit for Chapter I & II
Write:
Debye Series
Programming:
Implementation of loops in Mathematica

Regarding Debye Series:
some papers and van den Hulst and original Debye paper

-----
done for today:
Modification in IDL programme. S1(theta) & S2(theta) ranges over n(1->nstop) and “nang” belongs to “theta”. The array for S1 & S2 are limited for allowed number of angles, so should be array of “Mie”.
In Mathematica, all the initialization should be made, then S1 & S2 must be calculated, where they have Table form.

Understanding Recursion

nang: number of angles between 0-90 degrees
dang: defines angles for calculation (kleinste Einheit)

For j: 1 -- nang
theta = (j-1) * nang
mu_j = cos(theta)
where dang = 0.5*Pi * 1/(nang-1)

mxnang: maximum allowable “nang”

Possibility I:
mu = Table[0,{mxnang}]
For [i=o, i<= nang, theta = (i-1)*dang; mu[[i+1]] = cos(theta); i++]

Possibility II:
mu = Table[theta=(i-1)*dang; cos(theta), {i,0,nang}]
mu = (theta=(i-1)*dang; cos(theta))&/@ Range[0,nang]
mu = Array[theta=(i-1)*dang; cos(theta)&, nang, 0]

Mathematica Reworking

Writing
◆ Explanation w.r.t. Debye Series

Programming
◆ Implementation of recursive algorithm in Mathematica
◆ Resurrecting the idea of saving data for Legendre Polynome
◆ Improving IDL code

Saving the data of Legendre Polynome and later reusing it for the Mathematica code doesn't save necessarily lotta time. Just 36% of time has been saved, i.e. in case of 100 mulitiplicity where it takes 20613 seconds, I can complete this with 2168.3865 min. i.e. 36 Hrs. instead of 57 Hrs.
The next step would be to try recursive algorithm in Mathematica.