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}]