PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoYlm.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoYlm - the class to calculate varous components of spherical        //
   4 //  harmonics                                                                 //
   5 //                                                                            //
   6 // Authors: Adam Kisiel kisiel@mps.ohio-state.edu                             //
   7 //                                                                            //
   8 ////////////////////////////////////////////////////////////////////////////////
   9
  10 #include "AliFemtoYlm.h"
  11 #include <TMath.h>
  12 #include <iostream>
  13
  14 double *AliFemtoYlm::fgPrefactors = 0x0;
  15 int    *AliFemtoYlm::fgPrefshift = 0x0;
  16 int    *AliFemtoYlm::fgPlmshift = 0x0;
  17
  18 AliFemtoYlm::AliFemtoYlm() {
  19   InitializeYlms();
  20 }
  21
  22 AliFemtoYlm::~AliFemtoYlm() {
  23   free(fgPrefactors);
  24   free(fgPrefshift);
  25   free(fgPlmshift);
  26 }
  27
  28
  29 AliFemtoYlm::AliFemtoYlm(const AliFemtoYlm& aYlm){
  30   fgPrefshift = aYlm.fgPrefshift;
  31   InitializeYlms();
  32 }
  33
  34 AliFemtoYlm& AliFemtoYlm::operator=(const AliFemtoYlm& aYlm){
  35   if (this == &aYlm)
  36     return *this;
  37
  38   InitializeYlms();
  39
  40   return *this;
  41 }
  42
  43 std::complex<double> AliFemtoYlm::Ceiphi(double phi){
  44   return std::complex<double>(cos(phi),sin(phi));
  45 }
  46
  47 double AliFemtoYlm::Legendre(int ell, int em, double ctheta){
  48   // Calculate a single Legendre value
  49   // *** Warning - NOT optimal - calculated all Plms up to L !!!
  50   double lbuf[36];
  51   AliFemtoYlm::LegendreUpToYlm(ell, ctheta, lbuf);
  52
  53   return lbuf[fgPlmshift[ell]-abs(em)];
  54 }
  55
  56 std::complex<double> AliFemtoYlm::Ylm(int ell,int m,double theta,double phi){
  57   // Calculate Ylm spherical input
  58   double ctheta;
  59   std::complex<double> answer;
  60   std::complex<double> ci(0.0,1.0);
  61   ctheta=cos(theta);
  62   answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi);
  63
  64   return answer;
  65 }
  66
  67 std::complex<double> AliFemtoYlm::Ylm(int ell, int m, double x, double y, double z){
  68   // Calculate Ylm cartesian input
  69   std::complex<double> answer;
  70   double ctheta,phi;
  71   double r = sqrt(x*x+y*y+z*z);
  72   if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0;
  73   else ctheta=z/r;
  74   phi=atan2(y,x);
  75   answer = (fgPrefactors[fgPrefshift[ell]+m]*Legendre(ell,m,ctheta))*Ceiphi(m*phi);
  76
  77   return answer;
  78 }
  79
  80 double AliFemtoYlm::ReYlm(int ell, int m, double theta, double phi){
  81         return real(AliFemtoYlm::Ylm(ell,m,theta,phi));
  82 }
  83
  84 double AliFemtoYlm::ImYlm(int ell, int m, double theta, double phi){
  85         return imag(AliFemtoYlm::Ylm(ell,m,theta,phi));
  86 }
  87
  88 double AliFemtoYlm::ReYlm(int ell, int m, double x,double y,double z){
  89         return real(AliFemtoYlm::Ylm(ell,m,x,y,z));
  90 }
  91
  92 double AliFemtoYlm::ImYlm(int ell, int m, double x,double y,double z){
  93         return imag(AliFemtoYlm::Ylm(ell,m,x,y,z));
  94 }
  95
  96 void AliFemtoYlm::InitializeYlms()
  97 {
  98   // Calculate prefactors for fast Ylm calculation
  99
 100   double oneoversqrtpi = 1.0/TMath::Sqrt(TMath::Pi());
 101
 102   if(fgPrefactors!=NULL)
 103     free(fgPrefactors);
 104   if(fgPrefshift!=NULL)
 105     free(fgPrefshift);
 106   if(fgPlmshift!=NULL)
 107     free(fgPlmshift);
 108
 109   fgPrefactors = (double *) malloc(sizeof(double) * 36);
 110   fgPrefshift  = (int *) malloc(sizeof(int) * 6);
 111   fgPlmshift   = (int *) malloc(sizeof(int) * 6);
 112
 113   // l=0 prefactors
 114   fgPrefactors[0] = 0.5*oneoversqrtpi;
 115
 116   // l=1 prefactors
 117   fgPrefactors[1] = 0.5*sqrt(3.0/2.0)*oneoversqrtpi;
 118   fgPrefactors[2] = 0.5*sqrt(3.0)*oneoversqrtpi;
 119   fgPrefactors[3] = -fgPrefactors[1];
 120
 121   // l=2 prefactors
 122   fgPrefactors[4] = 0.25*sqrt(15.0/2.0)*oneoversqrtpi;
 123   fgPrefactors[5] = 0.5*sqrt(15.0/2.0)*oneoversqrtpi;
 124   fgPrefactors[6] = 0.25*sqrt(5.0)*oneoversqrtpi;
 125   fgPrefactors[7] = -fgPrefactors[5];
 126   fgPrefactors[8] = fgPrefactors[4];
 127
 128   // l=3 prefactors
 129   fgPrefactors[9]  = 0.125*sqrt(35.0)*oneoversqrtpi;
 130   fgPrefactors[10] = 0.25*sqrt(105.0/2.0)*oneoversqrtpi;
 131   fgPrefactors[11] = 0.125*sqrt(21.0)*oneoversqrtpi;
 132   fgPrefactors[12] = 0.25*sqrt(7.0)*oneoversqrtpi;
 133   fgPrefactors[13] = -fgPrefactors[11];
 134   fgPrefactors[14] = fgPrefactors[10];
 135   fgPrefactors[15] = -fgPrefactors[9];
 136
 137   // l=4 prefactors
 138   fgPrefactors[16] = 3.0/16.0*sqrt(35.0/2.0)*oneoversqrtpi;
 139   fgPrefactors[17] = 3.0/8.0*sqrt(35.0)*oneoversqrtpi;
 140   fgPrefactors[18] = 3.0/8.0*sqrt(5.0/2.0)*oneoversqrtpi;
 141   fgPrefactors[19] = 3.0/8.0*sqrt(5.0)*oneoversqrtpi;
 142   fgPrefactors[20] = 3.0/16.0*oneoversqrtpi;
 143   fgPrefactors[21] = -fgPrefactors[19];
 144   fgPrefactors[22] = fgPrefactors[18];
 145   fgPrefactors[23] = -fgPrefactors[17];
 146   fgPrefactors[24] = fgPrefactors[16];
 147
 148   // l=5 prefactors
 149   fgPrefactors[25] = 3.0/32.0*sqrt(77.0)*oneoversqrtpi;
 150   fgPrefactors[26] = 3.0/16.0*sqrt(385.0/2.0)*oneoversqrtpi;
 151   fgPrefactors[27] = 1.0/32.0*sqrt(385.0)*oneoversqrtpi;
 152   fgPrefactors[28] = 1.0/8.0*sqrt(1155.0/2.0)*oneoversqrtpi;
 153   fgPrefactors[29] = 1.0/16.0*sqrt(165.0/2.0)*oneoversqrtpi;
 154   fgPrefactors[30] = 1.0/16.0*sqrt(11.0)*oneoversqrtpi;
 155   fgPrefactors[31] = -fgPrefactors[29];
 156   fgPrefactors[32] = fgPrefactors[28];
 157   fgPrefactors[33] = -fgPrefactors[27];
 158   fgPrefactors[34] = fgPrefactors[26];
 159   fgPrefactors[35] = -fgPrefactors[25];
 160
 161   fgPrefshift[0] = 0;
 162   fgPrefshift[1] = 2;
 163   fgPrefshift[2] = 6;
 164   fgPrefshift[3] = 12;
 165   fgPrefshift[4] = 20;
 166   fgPrefshift[5] = 30;
 167
 168   fgPlmshift[0] = 0;
 169   fgPlmshift[1] = 2;
 170   fgPlmshift[2] = 5;
 171   fgPlmshift[3] = 9;
 172   fgPlmshift[4] = 14;
 173   fgPlmshift[5] = 20;
 174 }
 175
 176 void AliFemtoYlm::LegendreUpToYlm(int lmax, double ctheta, double *lbuf)
 177 {
 178   // Calculate a set of legendre polynomials up to a given l
 179   // with spherical input
 180   double sins[6];
 181   double coss[6];
 182   sins[0] = 0.0;
 183   coss[0] = 1.0;
 184   sins[1] = sqrt(1-ctheta*ctheta);
 185   coss[1] = ctheta;
 186   for (int iter=2; iter<6; iter++) {
 187     sins[iter] = sins[iter-1]*sins[1];
 188     coss[iter] = coss[iter-1]*coss[1];
 189   }
 190
 191   // Legendre polynomials l=0
 192   lbuf[0] = 1.0;
 193
 194   // Legendre polynomials l=1
 195   if (lmax>0) {
 196     lbuf[1] = sins[1];
 197     lbuf[2] = coss[1];
 198   }
 199
 200   // Legendre polynomials l=2
 201   if (lmax>1) {
 202     lbuf[3] = sins[2];
 203     lbuf[4] = sins[1]*coss[1];
 204     lbuf[5] = 3*coss[2]-1;
 205   }
 206
 207   // Legendre polynomials l=3
 208   if (lmax>2) {
 209     lbuf[6] = sins[3];
 210     lbuf[7] = sins[2]*coss[1];
 211     lbuf[8] = (5*coss[2]-1)*sins[1];
 212     lbuf[9] = 5*coss[3]-3*coss[1];
 213   }
 214
 215   // Legendre polynomials l=4
 216   if (lmax>3) {
 217     lbuf[10] = sins[4];
 218     lbuf[11] = sins[3]*coss[1];
 219     lbuf[12] = (7*coss[2]-1)*sins[2];
 220     lbuf[13] = (7*coss[3]-3*coss[1])*sins[1];
 221     lbuf[14] = 35*coss[4]-30*coss[2]+3;
 222   }
 223
 224   // Legendre polynomials l=5
 225   if (lmax>4) {
 226     lbuf[15] = sins[5];
 227     lbuf[16] = sins[4]*coss[1];
 228     lbuf[17] = (9*coss[2]-1)*sins[3];
 229     lbuf[18] = (3*coss[3]-1*coss[1])*sins[2];
 230     lbuf[19] = (21*coss[4]-14*coss[2]+1)*sins[1];
 231     lbuf[20] = 63*coss[5]-70*coss[3]+15*coss[1];
 232   }
 233 }
 234
 235 void AliFemtoYlm::YlmUpToL(int lmax, double x, double y, double z, std::complex<double> *ylms)
 236 {
 237   // Calculate a set of Ylms up to a given l
 238   // with cartesian input
 239   double ctheta,phi;
 240
 241   double r = sqrt(x*x+y*y+z*z);
 242   if ( r < 1e-10 || fabs(z) < 1e-10 ) ctheta = 0.0;
 243   else ctheta=z/r;
 244   phi=atan2(y,x);
 245
 246   YlmUpToL(lmax, ctheta, phi, ylms);
 247
 248 }
 249
 250 void AliFemtoYlm::YlmUpToL(int lmax, double ctheta, double phi, std::complex<double> *ylms)
 251 {
 252   // Calculate a set of Ylms up to a given l
 253   // with spherical input
 254   int lcur = 0;
 255   double lpol;
 256
 257   double coss[6];
 258   double sins[6];
 259
 260   double lbuf[36];
 261   LegendreUpToYlm(lmax, ctheta, lbuf);
 262
 263   for (int iter=1; iter<=lmax; iter++) {
 264     coss[iter-1] = cos(iter*phi);
 265     sins[iter-1] = sin(iter*phi);
 266   }
 267   ylms[lcur++] = fgPrefactors[0]*lbuf[0] * std::complex<double>(1,0);
 268
 269   for (int il = 1; il<=lmax; il++) {
 270     // First im = 0
 271     ylms[lcur+il] = fgPrefactors[fgPrefshift[il]]*lbuf[fgPlmshift[il]]*std::complex<double>(1.0,0.0);
 272     // Im != 0
 273     for (int im=1; im<=il; im++) {
 274       lpol = lbuf[fgPlmshift[il]-im];
 275       ylms[lcur+il-im] = fgPrefactors[fgPrefshift[il]-im]*lpol*std::complex<double>(coss[im-1],-sins[im-1]);
 276       ylms[lcur+il+im] = fgPrefactors[fgPrefshift[il]+im]*lpol*std::complex<double>(coss[im-1],sins[im-1]);
 277     }
 278     lcur += 2*il + 1;
 279   }
 280 }