#include <stdlib.h>
#include <stdio.h>
#include "machdefs.h"
#include <iostream>
#include <float.h>
#include <math.h>
#include "harmspher.h"

using namespace std;

/*
----------------------------------------------------------------
 ----- Calcul des fonctions de Legendre Plm(cos(Theta)) ------ 
----------------------------------------------------------------
       (partie theta des harmoniques spheriques)

 HarmSph_array(int lmax,int m,double cost,double *Plm)
 HarmSph_array_teta(int lmax,int m,double teta,double *Plm)

 - Input:
     lmax, m and teta or cos(teta)    (lmax>=m)
     array Plm with size at least "lmax-m+1"
     (calling with HarmSph_array_teta for small teta is more accurate)

 - Action:
     compute Psph for fixed "m" and all "l" from "m" to "lmax"

 - Return array: Psph(l,m)
     index Psph(l,m) = Plm[l-m]

 - Return value:
     Psph(lmax,m) = Plm[lmax-m]

 - Rappel:    x = cos(teta)
   Polynomes de Legendre associes: P(l,m)
   Polynomes de Legendre pour les harmoniques spheriques: Psph(l,m)
                (souvent appelles lambda_lm)
      Psph(l,m) = sqrt((2l+1)/4Pi) * sqrt((l-m)!/(l+m)!) * P(l,m)
----------------------------------------------------------------
*/

double HarmSph_array_teta(int lmax,int m,double teta,double *xlm)
{
 // Debug level
 const int lp=0;
 // Si une valeur finale est <10^EXPO_MINI on retourne zero
 const double EXPO_MINI = DBL_MIN_10_EXP +5;
 // Si une valeur dans la reccurence est <FACTMINI on renormalise
 const double FACTMINI = 1.e-50;

 if(m<0 || m>lmax) {
   cout<<"HarmSph_array_Error: bad arguments lmax="<<lmax<<" m="<<m<<endl;
   return(0.);
 }

 double x = cos(teta);
 if(lp) cout<<"HarmSph_array_teta: x="<<x<<" lmax="<<lmax<<" m="<<m<<endl;

 /* ---- calcul Psph(0,0) = 1/sqrt(4Pi) */
 xlm[0] = M_2_SQRTPI/4.;
 if(m==0 && lmax==0) return xlm[0];

 /* ---- calcul Psph(m,m):
 Recurrence:
   P(l,m) = (-1)^m * (2m-1)!! * (1-x^2)^(m/2)
   Psph(m+1,m+1) = -Psph(m,m) * sin(teta)*sqrt((2m+3)/(2m+2))
   ou Psph(m,m) = - Psph(m-1,m-1) * sin(teta)*sqrt((2m+1)/2m)
   et Psph(0,0) = 1/sqrt(4Pi)
 Pb: pour m grand et sin(teta) petit sin(teta)^m devient
     plus petit que la limite machine (1e-300) et on obtient zero
 */
 double expos = 0.;
 if(m>0) {          // Psph(0,0) == xlm[0]
   //double st = sqrt((1.-x)*(1.+x));
   double st = fabs(sin(teta));
   for(int i=1;i<=m;i++) {
     xlm[0] *= -sqrt(((double)i+0.5)/(double)i) *st;
     double fxlm = fabs(xlm[0]);
     if(fxlm<FACTMINI && fxlm>0.) {
       if(lp>1) cout<<i<<" renorm: xlm[0]="<<xlm[0]<<" e="<<log10(fxlm);
       xlm[0] = (xlm[0]>0.)? 1.: -1.;
       expos += log10(fxlm);
       if(lp>1) cout<<" --> "<<xlm[0]<<", expos="<<expos<<endl;
     }
   }
 }
 if(lp) cout<<"End computing Pmm="<<xlm[0]<<", expos="<<expos<<endl;

 /* ---- calcul Psph(m+1,m):
 Formule:
   P(m+1,m) = cos(teta)*sqrt(2m+1)*P(m,m)
   Psph(m+1,m) = cos(teta)*sqrt(2m+3)*Psph(m,m)
 En fait cela s'integre sans probleme dans la recurrence generale ci-dessous
 */
 /* ---- calcul de Psph(l,m)  pour l>m+1 :
 Recurrence:
   P(l,m)  = [ (2l-1)*x * P(l-1,m) - (l-1+m) * P(l-2,m) ] / (l-m)
   Psph(l,m) = sqrt((4*l^2-1)/(l^2-m^2))
             * [ cos(teta) * Psph(l-1,m)
               - sqrt(((l-1)^2-m^2)/(4*(l-1)^2-1)) * Psph(l-2,m) ]
   Psph(l,m) = sqrt((2l-1)*(2l+1)/((l-m)*(l+m)))
             * [ cos(teta) * Psph(l-1,m)
               - sqrt((l-m-1)*(l-1+m)/((2l-3)*(2l-1))) * Psph(l-2,m) ]

 */
 if(lmax>m) {
   const double expo_factmini = log10(FACTMINI);
   double factmaxi = 1./FACTMINI;
   double pm2=0., pm1=xlm[0], b_rec=0.;
   for(int ll=m+1;ll<=lmax;ll++) {
     int il = ll-m;
     double a_rec = sqrt( (2.*ll-1.)*(2.*ll+1.)/((ll-m)*(ll+m)) );
     xlm[il] = a_rec*(x*pm1-b_rec*pm2);
     double fxlm = fabs(xlm[il]);
     // Par construction "expos" est negatif ou nul (car FACTMINI<1)
     if(expos<0. && fxlm>factmaxi) {
       // On renormalise par factmaxi si on a encore de la reserve d'exposant
       double fact,efact;
       if(expos<expo_factmini) {
         fact  = factmaxi;
         efact = -expo_factmini;
         expos += efact;
       } else {
         fact = pow(10.,expos);
         efact = expos;
         expos = 0.;
         if(lp>1) cout<<"!! expos==0 !!"<<endl;
       }
       if(lp>1) cout<<ll<<" renorm: xlm["<<il<<"]="<<xlm[il]<<", expos="<<expos<<endl;
       for(int i=0;i<=ll-m;i++) xlm[i] /= fact;
     }
     pm2 = xlm[il-1];
     pm1 = xlm[il];
     b_rec = 1./a_rec;
   }
 }

 if(lp) cout<<"Return: Final expos="<<expos<<endl;
 if(expos!=0.) {
   for(int i=0;i<=lmax-m;i++) {
     double fxml=fabs(xlm[i]);
     if(lp>1) cout<<"     xlm["<<i<<"]="<<xlm[i];
     if(fxml>0.) {
       double e = log10(fxml);
       if(lp>1) cout<<" (e="<<e<<")";
       if(expos+e<EXPO_MINI) xlm[i] = 0.;
       else if(expos>EXPO_MINI) xlm[i] *= pow(10.,expos);
       else xlm[i] = pow(10.,expos+e);
     }
     if(lp>1) cout<<" ---> "<<xlm[i]<<endl;
   }
 }
 return xlm[lmax-m];

}