#include <math.h>
#include "lambuilder.h"
#include <iostream.h>
#ifdef __MWERKS__
#include "unixmac.h"
#endif

//++ 
// Class	Lambda2Builder
//
//    (see separate documentation)
//--

const double LambdaBuilder::bignorm = 1.e268; // = 1e-20*1.d288

Lambda2Builder::Lambda2Builder(){}

Lambda2Builder::Lambda2Builder(double theta, int lmax, int mmax):
  LambdaBuilder(theta, lmax, mmax){
    fact_=0;
    w_=0;
    ix_=0;
}

Lambda2Builder::Lambda2Builder(const Lambda2Builder& lb):
  LambdaBuilder(lb){
    fact_=lb.fact_;
    w_=lb.w_;
    ix_=lb.ix_;
}

void Lambda2Builder::reset(){
  LambdaBuilder::reset();
  fact_=0;
  w_=0;
  ix_=0;
}

double Lambda2Builder::lam2lmp(int flip){
  if (flip < 0){
    return -fact_*par_lm_*(w_+ix_);
  } else{
    return lam2lmp();
  }
}
 
double Lambda2Builder::lam2lmp(int l, int m){
  if ((l <= lmax_ && l==l_+1 && m==m_) ||
      (m <= mmax_ && m==m_+1 && l==m && l_==lmax_)){
    step();
    return lam2lmp();
  } else if (l==l_ && m==m_){
    return lam2lmp();
  } else if ((l < l_ && m <= mmax_) || (m == m_ && l < l_) || m<m_){
    reset();
    while (l != l_ || m != m_){
      step();
    }
    return lam2lmp();
  } else if (l <= lmax_ && m <= mmax_ && (m > m_ || (m==m_ && l > l_))){ 
    while (l != l_ || m != m_){
      step();
    }
    return lam2lmp();
  } else{
    return -1;
  }
}

double Lambda2Builder::lam2lmp(int l, int m, int flip){

  if ((l <= lmax_ && l==l_+1 && m==m_) ||
      (m <= mmax_ && m==m_+1 && l==m && l_==lmax_)){
    step();
    return lam2lmp(flip);
  } else if (l==l_ && m==m_){
    return lam2lmp(flip);
  } else if ((l < l_ && m <= mmax_) || (m == m_ && l < l_) || m<m_){
    reset();
    while (l != l_ || m != m_){
      step();
    }
    return lam2lmp(flip);
  } else if (l <= lmax_ && m <= mmax_ && (m > m_ || (m==m_ && l > l_))){ 
    while (l != l_ || m != m_){
      step();
    }
    return lam2lmp(flip);
  } else{
    return -1;
  }
}

double Lambda2Builder::lam2lmm(int flip){
  if (flip < 0){
    return -fact_*par_lm_*(w_-ix_);
  } else{
    return lam2lmp();
  }
}

double Lambda2Builder::lam2lmm(int l, int m){
  if ((l <= lmax_ && l==l_+1 && m==m_) ||
      (m <= mmax_ && m==m_+1 && l==m && l_==lmax_)){
    step();
    return lam2lmm();
  } else if (l==l_ && m==m_){
    return lam2lmm();
  } else if ((l < l_ && m <= mmax_) || (m == m_ && l < l_) || m<m_){
    reset();
    while (l != l_ || m != m_){
      step();
    }
    return lam2lmm();
  } else if (l <= lmax_ && m <= mmax_ && (m > m_ || (m==m_ && l > l_))){ 
    while (l != l_ || m != m_){
      step();
    }
    return lam2lmm();
  } else{
    return -1;
  }
}

double Lambda2Builder::lam2lmm(int l, int m, int flip){
  if ((l <= lmax_ && l==l_+1 && m==m_) ||
      (m <= mmax_ && m==m_+1 && l==m && l_==lmax_)){
    step();
    return lam2lmm(flip);
  } else if (l==l_ && m==m_){
    return lam2lmm(flip);
  } else if ((l < l_ && m <= mmax_) || (m == m_ && l < l_) || m<m_){
    reset();
    while (l != l_ || m != m_){
      step();
    }
    return lam2lmm(flip);
  } else if (l <= lmax_ && m <= mmax_ && (m > m_ || (m==m_ && l > l_))){ 
    while (l != l_ || m != m_){
      step();
    }
    return lam2lmm(flip);
  } else{
    return -1;
  }
}

void Lambda2Builder::step(){
  if (l_>=lmax_ && m_>=mmax_){return;}
  if (l_==lmax_){//move on to the next m and calculate l=m
    LambdaBuilder::step();
    if (l_>=2){
      w_=-2.*(m_*(2.*m_+1.)-m_*(m_-1.)/sth_/sth_)*lamlm_;
      ix_=-2.*m_*(m_-1.)*cth_/sth_/sth_*lamlm_;
      fact_=sqrt(1./(m_-1.)/m_/(m_+1.)/(m_+2.));
    } else{
      w_=0;
      ix_=0;
      fact_=0;
    }
  } else{
    double lastlamlm=lamlm_;
    LambdaBuilder::step();
    if (l_>=2){    
      w_=-2.*(-((l_-m_*m_)/sth_/sth_+.5*l_*(l_-1.))*lamlm_+
	      (l_+m_)*cth_/sth_/sth_*sqrt((2.*l_+1.)/(2.*l_-1.)*(l_-m_)/(l_+m_))*
	      lastlamlm);
      ix_=-2.*m_/sth_/sth_*((l_-1.)*cth_*lamlm_-(l_+m_)*
			    sqrt((2.*l_+1.)/(2.*l_-1.)*(l_-m_)/(l_+m_))*
			    lastlamlm);
      fact_=sqrt(1./(l_-1.)/l_/(l_+1.)/(l_+2.));
    } else{
      w_=0;
      ix_=0;
      fact_=0;
    }
  }
}

LambdaBuilder::LambdaBuilder(){}

LambdaBuilder::LambdaBuilder(double theta, int lmax, int mmax){
  lmax_=lmax;
  mmax_=mmax;
  l_=0;
  m_=0;
  lamlm_=1. / sqrt(4.*M_PI);
  cth_=cos(theta);
  sth_=sin(theta);
  lammm_=1. / sqrt(4.*M_PI) *bignorm;
  par_lm_=1.;
  lam_0_=0.;
  lam_1_=1./bignorm;
  a_rec_=sqrt(3.);
  b_rec_=0;
  fm2_=0;
}

LambdaBuilder::LambdaBuilder(const LambdaBuilder& lb){ 
  lmax_=lb.lmax_;
  mmax_=lb.mmax_;
  l_=lb.l_;
  m_=lb.m_;
  lamlm_=lb.lamlm_;
  cth_=lb.cth_;
  sth_=lb.sth_;
  lammm_=lb.lammm_;
  par_lm_=lb.par_lm_;
  lam_0_=lb.lam_0_;
  lam_1_=lb.lam_1_;
  a_rec_=lb.a_rec_;
  b_rec_=lb.b_rec_;
  fm2_=lb.fm2_;
}

void LambdaBuilder::reset(){
  l_=0;
  m_=0;
  lamlm_=1. / sqrt(4.*M_PI);
  lammm_=1. / sqrt(4.*M_PI) *bignorm;
  par_lm_=1.;
  lam_0_=0.;
  lam_1_=1./bignorm;
  a_rec_=sqrt(3.);
  b_rec_=0;
  fm2_=0;
}

double LambdaBuilder::lamlm(int flip){
  if (flip < 0){
    return par_lm_*lamlm();
  } else{
    return lamlm();
  }
}

double LambdaBuilder::lamlm(int l, int m){
  if ((l <= lmax_ && l==l_+1 && m==m_) ||
      (m <= mmax_ && m==m_+1 && l==m && l_==lmax_)){
    step();
    return lamlm();
  } else if (l==l_ && m==m_){
    return lamlm_;
  } else if ((l < l_ && m <= mmax_) || (m == m_ && l < l_) || m<m_){
    reset();
    while (l != l_ || m != m_){
      step();
    }
    return lamlm();
  } else if (l <= lmax_ && m <= mmax_ && (m > m_ || (m==m_ && l > l_))){ 
    while (l != l_ || m != m_){
      step();
    }
    return lamlm();
  } else{
    return -1;
  }
}

double LambdaBuilder::lamlm(int l, int m, int flip){
  if ((l <= lmax_ && l==l_+1 && m==m_) ||
      (m <= mmax_ && m==m_+1 && l==m && l_==lmax_)){
    step();
    return lamlm(flip);
  } else if (l==l_ && m==m_){
    return lamlm(flip);
  } else if ((l < l_ && m <= mmax_) || (m == m_ && l < l_) || m<m_){
    reset();
    while (l != l_ || m != m_){
      step();
    }
    return lamlm(flip);
  } else if (l <= lmax_ && m <= mmax_ && (m > m_ || (m==m_ && l > l_))){ 
    while (l != l_ || m != m_){
      step();
    }
    return lamlm(flip);
  } else{
    return -1;
  }
}

void LambdaBuilder::step(){
  /*-----------------------------------------------------------------------
    computes lambda_lm(theta) in an optimal fashion
    looping over m
    For example, if nlmax =3 and nmmax=3 the order is
    (l,m)=(0,0),(1,0),(2,0),(3,0),(1,1),(2,1),(3,1),(2,2),(3,2),(3,3)
    -----------------------------------------------------------------------*/
  if (l_>=lmax_ && m_>=mmax_){return;}
  if (l_==lmax_){//move on to the next m and calculate l=m
    m_=m_+1;
    l_=m_;
    double f2m = 2. * m_;
    fm2_ = ((double) m_)*m_;
    par_lm_ = 1.;  // = (-1)^(l+m)
    if (m_ >= 1) { /* lambda_0_0 for m>0*/
      lammm_ = - lammm_ * sth_ * sqrt( (f2m+1.)/ f2m );
    }
    lamlm_ = lammm_ / bignorm;
    lam_0_ = 0.;
    lam_1_ = 1. / bignorm;   /*small number --> to avoid overflow*/
    a_rec_ = sqrt( f2m + 3.);
    b_rec_ = 0.;
  } else{
    l_=l_+1;
    par_lm_ = - par_lm_; /*  = (-1)^(l+m)*/
    double lam_2 = (cth_ * lam_1_ - b_rec_ * lam_0_) * a_rec_;
    //    cout<<lam_2<<" "<<cth_<<" "<<lam_1_<<" "<<b_rec_<<" "<<lam_0_<<" "<<a_rec_<<endl;
    //actual lambda_lm (small and big numbers cancel out)
    lamlm_ = lam_2 * lammm_;
    
    b_rec_ = 1. / a_rec_;
    double fl2 = (l_+1.)*(l_+1.);
    a_rec_ = sqrt( (4. * fl2 - 1.) / (fl2-fm2_) );
    //    cout <<fl2<<" "<<fm2<<endl;
    lam_0_ = lam_1_;
    lam_1_ = lam_2;
  }
}