Index: unk/SophyaLib/Samba/lambuilder.cc
===================================================================
--- /trunk/SophyaLib/Samba/lambuilder.cc	(revision 729)
+++ 	(revision )
@@ -1,320 +1,0 @@
-#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;
-  }
-}
Index: unk/SophyaLib/Samba/lambuilder.h
===================================================================
--- /trunk/SophyaLib/Samba/lambuilder.h	(revision 729)
+++ 	(revision )
@@ -1,220 +1,0 @@
-#ifndef LAMBUILDER_H_SEEN
-#define LAMBUILDER_H_SEEN
-
-#include <math.h>
-
-
- class LambdaBuilder{
- public:
-  LambdaBuilder();
-  LambdaBuilder(double theta,int lmax, int mmax);
-  LambdaBuilder(const LambdaBuilder& lb);
-  void reset();   
-  double lamlm() const {return lamlm_;}
-  double lamlm(int flip);
-  double lamlm(int l, int m);
-  double lamlm(int l, int m, int flip);
-  int l() const {return l_;}
-  int m() const {return m_;}
-  int lmax() const {return lmax_;}
-  int mmax() const {return mmax_;}
-  void step();
- protected:
-  double par_lm_;
-  int lmax_;
-  int mmax_;
-  int l_;
-  int m_;
-  double lamlm_;
-  double cth_;
-  double sth_;
-  double lammm_;
-  double lam_0_;
-  double lam_1_;
-  double a_rec_;
-  double b_rec_;
-  double fm2_;
-  //  static const double bignorm=1e-20*HUGE_VAL;
-//  static const double bignorm=1e-20*1.d288;
-  static const double bignorm;
-};
-
-class Lambda2Builder : public LambdaBuilder {
- public:
-  Lambda2Builder();
-  Lambda2Builder(double theta,int lmax, int mmax);
-  Lambda2Builder(const Lambda2Builder& lb);
-  void reset();
-  double lam2lmp() const {return -fact_*(w_-ix_);}
-  double lam2lmp(int flip);
-  double lam2lmp(int l, int m);
-  double lam2lmp(int l, int m, int flip);
-  double lam2lmm() const {return -fact_*(w_+ix_);}
-  double lam2lmm(int flip);
-  double lam2lmm(int l, int m);
-  double lam2lmm(int l, int m, int flip);
-  void step();
- private:
-  double fact_;
-  double w_;
-  double ix_;
-};
-#endif
-
-/*! \class LambdaBuilder
- \brief A lambda building class
-
-Objects in this class generate
-\f[
-            \lambda_l^m=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}
-	    P_l^m(\cos{\theta})
-\f]
-where \f$P_l^m\f$ are the associated Legendre polynomials.  The interest
-of course is that
-\f[
-            Y_{lm}(\cos{\theta})=\lambda_l^m(\cos{\theta}) e^{im\phi}.
-\f]
-
-Each object has a fixed theta, and maximum l and m to be calculated
-(lmax and mmax).
-
-The algorithm solves the \f$\lambda_l^m\f$'s iteratively so there
-is an optimal sequance of l and m's as follows...
-
-\code
-for (int m = 0; m <= nmmax; m++){
-  dum=lb.lamlm(m,m) 
-  for (int l = m+1; l<= nlmax; l++){
-    dum = lb.lamlm(l,m)
-  }
-}\endcode
-The function LambdaBuilder::step() automatically shifts the l and
-m's to the next optimal value.
-However, \f$\lambda_l^m\f$ can be calculated for any l and m less than
-or equal to lmax and lmax by going through the loop or restarting the
-loop.
-
-\f$\lambda_l^m(-\cos{\theta})\f$ can be quickly obtained using
-symmetry relations, which can be more efficient that creating a new
-LambdaBuilder object.
-
-The algorithm is based on the HEALPIX code.
-*/
-
-/*!  \fn LambdaBuilder::LambdaBuilder(double theta,int lmax, int mmax)
-  \param theta the desired \f$\theta\f$
-  \param lmax the maximum value of l
-  \param mmax the maximum value of m
-*/
-/*!
-  \fn LambdaBuilder::LambdaBuilder(const LambdaBuilder& lb)
-  \param lb the LambdaBuilder object to be copied
-*/
-/*! \fn void LambdaBuilder::reset()
-  \brief Resets the object to l=0, m=0
-*/
-/*! \fn double LambdaBuilder::lamlm() const
-  \returns \f$\lambda_l^m\f$ for current l and m
-*/
-/*!
-  \fn double LambdaBuilder::lamlm(int flip)
-  \param flip \f$\theta\f$ to be multiplied by sign(flip)
-  \returns \f$\lambda_l^m(\cos{(sign(flip)\theta}))\f$ for the current l and m.
-  Uses symmetry relations to get this value.
-*/
-/*! \fn double LambdaBuilder::lamlm(int l, int m)
-  \param l the desired l
-  \param m the desired m
-  \returns \f$\lambda_l^m\f$ for l and m
-*/
-/*!   \fn double LambdaBuilder::lamlm(int l, int m, int flip)
-  \param l the desired l
-  \param m the desired m
-   \param flip \f$\theta\f$ to be multiplied by sign(flip)
-  \returns \f$\lambda_l^m(\cos{(sign(flip)\theta}))\f$ for the given l and m.
-  Uses symmetry relations to get this value.
-*/
-/*!   \fn int LambdaBuilder::l() const
-  \returns current value of l
-*/
-/*!   \fn int LambdaBuilder::m() const
-  \returns current value of m
-*/
-/*!   \fn int LambdaBuilder::lmax() const
-  \returns the maximum value of l
-*/
-/*!   \fn int LambdaBuilder::mmax() const
-  \returns the maximum value of m
-*/
-/*!   \fn void LambdaBuilder::step()
-  \brief Steps to the next l and m
-*/
-
-/*! \class Lambda2Builder
- \brief A lambda building class
-
-Objects in this class generate
-\f[
-            _\pm\lambda_l^m=2\sqrt{\frac{(l-2)!}{(l+2)!}\frac{(2l+1)}{4\pi}\frac{(l-m)!}{(l+m)!}}\left( G^+_{lm} \mp G^-_{lm}\right)
-\f]
-where
-\f[G^+_{lm}(\cos{\theta})=-\left( \frac{l-m^2}{\sin^2{\theta}}+\frac{1}{2}l\left(l-1\right)\right)P_l^m(\cos{\theta})+\left(l+m\right)\frac{\cos{\theta}}{\sin^2{\theta}}P^m_{l-1}(\cos{\theta})
-\f]
-and
-\f[G^-_{lm}(\cos{\theta})=\frac{m}{\sin^2{\theta}}\left(\left(l-1\right)\cos{\theta}P^m_l(\cos{\theta})-\left(l+m\right)P^m_{l-1}(\cos{\theta})\right)
-\f]
-and \f$P_l^m\f$ are the associated Legendre polynomials.
-The interest in all this is that
-\f[_{\pm2}Y_l^m(\cos{\theta})=_\pm\lambda_l^m(\cos{\theta})e^{im\phi}
-\f]
-where \f$_{\pm2}Y_l^m\f$ are the spin-2 spherical harmonics.
-*/
-
-/*!   \fn  Lambda2Builder::Lambda2Builder(double theta,int lmax, int mmax)
-  \param theta the desired \f$\theta\f$
-  \param lmax the maximum value of l
-  \param mmax the maximum value of m
-*/
-/*!   \fn  Lambda2Builder::Lambda2Builder(const Lambda2Builder& lb)
-  \param lb the Lambda2Builder object to be copied
-*/
-/*!   \fn  double Lambda2Builder::lam2lmp() const
-  \returns \f$_+\lambda_l^m\f$ for current l and m
-*/
-/*!   \fn  double Lambda2Builder::lam2lmp(int flip);
-  \param flip \f$\theta\f$ to be multiplied by sign(flip)
-  \returns \f$_+\lambda_l^m(\cos{(sign(flip)\theta}))\f$ for the current l and m.
-  Uses symmetry relations to get this value.
-*/
-/*!   \fn  double Lambda2Builder::lam2lmp(int l, int m)
-  \param l the desired l
-  \param m the desired m
-  \returns \f$_+\lambda_l^m\f$ for l and m
-*/
-/*!   \fn  double Lambda2Builder::lam2lmp(int l, int m, int flip)
-  \param l the desired l
-  \param m the desired m
-   \param flip \f$\theta\f$ to be multiplied by sign(flip)
-  \returns \f$_+\lambda_l^m(\cos{(sign(flip)\theta}))\f$ for the given l and m.
-  Uses symmetry relations to get this value.
-*/
-/*!   \fn  double Lambda2Builder::lam2lmm() const
-  \returns \f$_-\lambda_l^m\f$ for current l and m
-*/
-/*!   \fn  double Lambda2Builder::lam2lmm(int flip)
-  \param flip \f$\theta\f$ to be multiplied by sign(flip)
-  \returns \f$_-\lambda_l^m(\cos{(sign(flip)\theta}))\f$ for the current l and m.
-  Uses symmetry relations to get this value.
-*/
-/*!   \fn  double Lambda2Builder::lam2lmm(int l, int m)
-  \param l the desired l
-  \param m the desired m
-  \returns \f$_-\lambda_l^m\f$ for l and m
-*/
-/*!   \fn  double Lambda2Builder::lam2lmm(int l, int m, int flip)
-  \param l the desired l
-  \param m the desired m
-   \param flip \f$\theta\f$ to be multiplied by sign(flip)
-  \returns \f$_-\lambda_l^m(\cos{(sign(flip)\theta}))\f$ for the given l and m.
-  Uses symmetry relations to get this value.
-)*/
