ylmgen.h

00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*
00028  *  Code for efficient calculation of Y_lm(theta,phi=0)
00029  *
00030  *  Copyright (C) 2005, 2006 Max-Planck-Society
00031  *  Author: Martin Reinecke
00032  */
00033 
00034 #ifndef PLANCK_YLMGEN_H
00035 #define PLANCK_YLMGEN_H
00036 
00037 #include <cmath>
00038 #include "arr.h"
00039 #include "lsconstants.h"
00040 
00041 /*! Class for efficient calculation of Y_lm(theta,phi=0) */
00042 class Ylmgen
00043   {
00044   private:
00045     double fsmall, fbig, eps, cth_crit;
00046     int lmax, mmax, m_last, m_crit;
00047     arr<double> cf;
00048     arr<double[2]> recfac;
00049     arr<double> mfac;
00050     arr<double> t1fac, t2fac;
00051 
00052     enum { large_exponent2 = 90, minscale=-4 };
00053 
00054     void recalc_recfac (int m)
00055       {
00056       using namespace std;
00057 
00058       if (m_last==m) return;
00059 
00060       double f_old=1;
00061       for (int l=m; l<recfac.size(); ++l)
00062         {
00063         recfac[l][0] = t1fac[l]*t2fac[l+m]*t2fac[l-m];
00064         recfac[l][1] = recfac[l][0]/f_old;
00065         f_old = recfac[l][0];
00066         }
00067 
00068       m_last=m;
00069       }
00070 
00071   public:
00072     /*! Creates a generator which will calculate Y_lm(theta,phi=0)
00073         up to \a l=l_max and \a m=m_max. It may regard Y_lm whose absolute
00074         magnitude is smaller than \a epsilon as zero. */
00075     Ylmgen (int l_max, int m_max, double epsilon=1e-30)
00076       : eps(epsilon), cth_crit(2.), lmax(l_max), mmax(m_max), m_last(-1),
00077         m_crit(mmax+1), cf(-minscale+11), recfac(lmax+1), mfac(mmax+1),
00078         t1fac(lmax+1), t2fac(lmax+mmax+1)
00079       {
00080       using namespace std;
00081 
00082       fsmall = ldexp(1.,-large_exponent2);
00083       fbig   = ldexp(1., large_exponent2);
00084       for (int m=0; m<cf.size(); ++m)
00085         cf[m] = ldexp(1.,(m+minscale)*large_exponent2);
00086 
00087       mfac[0] = 1;
00088       for (int m=1; m<mfac.size(); ++m)
00089         mfac[m] = mfac[m-1]*sqrt((2*m+1.)/(2*m));
00090       for (int m=0; m<mfac.size(); ++m)
00091         mfac[m] = inv_ln2*log(inv_sqrt4pi*mfac[m]);
00092       for (int l=0; l<t1fac.size(); ++l)
00093         t1fac[l] = sqrt(4.*(l+1)*(l+1)-1.);
00094       for (int i=0; i<t2fac.size(); ++i)
00095         t2fac[i] = 1./sqrt(i+1.);
00096       }
00097 
00098     /*! For a colatitude given by \a cth and \a sth (representing cos(theta)
00099         and sin(theta)) and a multipole moment \a m, calculate the
00100         Y_lm(theta,phi=0) for \a m<=l<=lmax and return in it \a result[l].
00101         On exit, \a firstl is the \a l index of the first Y_lm with an
00102         absolute magnitude larger than \a epsilon. If \a firstl>lmax, all
00103         absolute values are smaller than \a epsilon.
00104         \a result[l] is undefined for all \a l<firstl. */
00105     void get_Ylm (double cth, double sth, int m, arr<double> &result,
00106       int &firstl)
00107       {
00108       using namespace std;
00109 
00110       planck_assert (m<=mmax, "get_Ylm: m larger than mmax");
00111 
00112       if (((m>=m_crit)&&(abs(cth)>=cth_crit)) || ((m>0)&&(sth==0)))
00113         { firstl=lmax+1; return; }
00114 
00115       recalc_recfac(m);
00116       result.alloc(lmax+1);
00117 
00118       double logval = mfac[m];
00119       if (m>0) logval += m*inv_ln2*log(sth);
00120       int scale = int (logval/large_exponent2)-minscale;
00121       double corfac = (scale<0) ? 0. : cf[scale];
00122       double lam_1 = 0;
00123       double lam_2 = exp(ln2*(logval-(scale+minscale)*large_exponent2));
00124       if (m&1) lam_2 = -lam_2;
00125 
00126       int l=m;
00127       while (true)
00128         {
00129         if (abs(lam_2*corfac)>eps) break;
00130         if (++l>lmax) break;
00131         double lam_0 = cth*lam_2*recfac[l-1][0] - lam_1*recfac[l-1][1];
00132         if (abs(lam_0*corfac)>eps) { lam_1=lam_2; lam_2=lam_0; break; }
00133         if (++l>lmax) break;
00134         lam_1 = cth*lam_0*recfac[l-1][0] - lam_2*recfac[l-1][1];
00135         if (abs(lam_1*corfac)>eps) { lam_2=lam_1; lam_1=lam_0; break; }
00136         if (++l>lmax) break;
00137         lam_2 = cth*lam_1*recfac[l-1][0] - lam_0*recfac[l-1][1];
00138 
00139         while (abs(lam_2)>fbig)
00140           {
00141           lam_1 *= fsmall;
00142           lam_2 *= fsmall;
00143           ++scale;
00144           corfac = (scale<0) ? 0. : cf[scale];
00145           }
00146         }
00147 
00148       firstl=l;
00149       if (l>lmax)
00150         { m_crit=m; cth_crit=abs(cth); return; }
00151 
00152       lam_1*=corfac;
00153       lam_2*=corfac;
00154 
00155       for (;l<lmax-2;l+=3)
00156         {
00157         result[l]=lam_2;
00158         double lam_0 = cth*lam_2*recfac[l][0] - lam_1*recfac[l][1];
00159         result[l+1] = lam_0;
00160         lam_1 = cth*lam_0*recfac[l+1][0] - lam_2*recfac[l+1][1];
00161         result[l+2] = lam_1;
00162         lam_2 = cth*lam_1*recfac[l+2][0] - lam_0*recfac[l+2][1];
00163         }
00164       while (true)
00165         {
00166         result[l]=lam_2;
00167         if (++l>lmax) break;
00168         double lam_0 = cth*lam_2*recfac[l-1][0] - lam_1*recfac[l-1][1];
00169         result[l] = lam_0;
00170         if (++l>lmax) break;
00171         lam_1 = cth*lam_0*recfac[l-1][0] - lam_2*recfac[l-1][1];
00172         result[l] = lam_1;
00173         if (++l>lmax) break;
00174         lam_2 = cth*lam_1*recfac[l-1][0] - lam_0*recfac[l-1][1];
00175         }
00176       }
00177 
00178   };
00179 
00180 #endif