LORENE: et_bin_bhns_extr_ylm.C Source File

00001 /*
00002  *  Method of class Et_bin_bhns_extr to construct spherical harmonics
00003  *
00004  *    (see file et_bin_bhns_extr.h for documentation).
00005  *
00006  */
00007 
00008 /*
00009  *   Copyright (c) 2004 Joshua A. Faber
00010  *
00011  *   This file is part of LORENE.
00012  *
00013  *   LORENE is free software; you can redistribute it and/or modify
00014  *   it under the terms of the GNU General Public License version 2
00015  *   as published by the Free Software Foundation.
00016  *
00017  *   LORENE is distributed in the hope that it will be useful,
00018  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00019  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00020  *   GNU General Public License for more details.
00021  *
00022  *   You should have received a copy of the GNU General Public License
00023  *   along with LORENE; if not, write to the Free Software
00024  *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00025  *
00026  */
00027 
00028 char et_bin_bhns_extr_ylm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_ylm.C,v 1.3 2005/02/28 23:18:07 k_taniguchi Exp $" ;
00029 
00030 /*
00031  * $Id: et_bin_bhns_extr_ylm.C,v 1.3 2005/02/28 23:18:07 k_taniguchi Exp $
00032  * $Log: et_bin_bhns_extr_ylm.C,v $
00033  * Revision 1.3  2005/02/28 23:18:07  k_taniguchi
00034  * Change the functions to constant ones
00035  *
00036  * Revision 1.2  2005/01/03 19:52:56  k_taniguchi
00037  * Change a factor multiplied/divided by sqrt(2).
00038  *
00039  * Revision 1.1  2004/12/29 16:30:46  k_taniguchi
00040  * *** empty log message ***
00041  *
00042  *
00043  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_ylm.C,v 1.3 2005/02/28 23:18:07 k_taniguchi Exp $
00044  *
00045  */
00046 
00047 // C headers
00048 #include <stdlib.h>
00049 #include <math.h>
00050 
00051 // Lorene headers
00052 #include "map.h"
00053 #include "tenseur.h"
00054 #include "et_bin_bhns_extr.h"
00055 
00056 void Et_bin_bhns_extr::get_ylm(int nylm, Cmp** ylmvec) const {
00057   
00058   //  IMPORTANT NOTE:
00059   // For Y_lm with m>=1, we have the real and imaginary parts, 
00060   // not Y_{l,m} and Y_{l,-m}.  This changes the normalization
00061   // properties.  In order to normalize properly, we multiply
00062   // all fields in get_integrals below by a factor of 2.0 when
00063   // m>=1.
00064 
00065   cout << "Constructing ylm" << endl;
00066 
00067   int nz = mp.get_mg()->get_nzone() ;
00068   int nr = mp.get_mg()->get_nr(0) ;  
00069   int np = mp.get_mg()->get_np(0) ;
00070   int nt = mp.get_mg()->get_nt(0) ;
00071 
00072   for (int l=0 ; l<nz ; l++) {    
00073 
00074     Mtbl Xabs (mp.x) ;
00075     Mtbl Yabs (mp.y) ;
00076     Mtbl Zabs (mp.z) ;
00077 
00078     for (int k=0 ; k<np ; k++) {
00079       for (int j=0 ; j<nt ; j++) {
00080         for (int i=0 ; i<nr ; i++) {
00081 
00082       double xval=Xabs(l,k,j,i);
00083       double yval=Yabs(l,k,j,i);
00084       double zval=Zabs(l,k,j,i);
00085       double rval=sqrt(xval*xval+yval*yval+zval*zval);
00086 
00087       //      cout <<l<<" "<<k<<" "<<j<<" "<<i<<endl;
00088 
00089       //l=0,m=0
00090       ylmvec[0]->set(l,k,j,i)=1.0*sqrt(1.0/4.0/M_PI);
00091       //      cout << " 0 " << endl;
00092       // l=1 included?
00093       if (nylm>1 ) {
00094         if (nylm <4) {abort();} else {
00095           //l=1,m=0
00096       ylmvec[1]->set(l,k,j,i)=zval*sqrt(3.0/4.0/M_PI);
00097       //l=1,m=1
00098       ylmvec[2]->set(l,k,j,i)=-1.0*xval*sqrt(3.0/8.0/M_PI);
00099       ylmvec[3]->set(l,k,j,i)=-1.0*yval*sqrt(3.0/8.0/M_PI);
00100         }
00101       }
00102       // l=2 included?
00103       if (nylm>4 ) {
00104         if (nylm <9) {abort();} else {
00105       //l=2,m=0
00106       ylmvec[4]->set(l,k,j,i)=(3.0*zval*zval-rval*rval)*sqrt(5.0/16.0/M_PI);
00107       //l=2,m=1
00108       ylmvec[5]->set(l,k,j,i)=-1.0*zval*xval*sqrt(15.0/8.0/M_PI);
00109       ylmvec[6]->set(l,k,j,i)=-1.0*zval*yval*sqrt(15.0/8.0/M_PI);
00110       //l=2,m=2
00111       ylmvec[7]->set(l,k,j,i)=(xval*xval-yval*yval)*sqrt(15.0/32.0/M_PI);
00112       ylmvec[8]->set(l,k,j,i)=2.0*xval*yval*sqrt(15.0/32.0/M_PI);
00113         }
00114       }
00115       // l=3 included?
00116       if (nylm>9 ) {
00117         if (nylm <16) {abort();} else {
00118       //l=3,m=0
00119       ylmvec[9]->set(l,k,j,i)=(5.0*pow(zval,3)-3.0*zval*rval*rval)*
00120         sqrt(7.0/16.0/M_PI);
00121       //l=3,m=1
00122       ylmvec[10]->set(l,k,j,i)=-1.0*(5.0*zval*zval-rval*rval)*xval*
00123         sqrt(21.0/64.0/M_PI);
00124       ylmvec[11]->set(l,k,j,i)=-1.0*(5.0*zval*zval-rval*rval)*yval*
00125         sqrt(21.0/64.0/M_PI);
00126       //l=3,m=2
00127       ylmvec[12]->set(l,k,j,i)=zval*(xval*xval-yval*yval)*
00128         sqrt(105./32.0/M_PI);
00129       ylmvec[13]->set(l,k,j,i)=zval*(2.0*xval*yval)*
00130         sqrt(105./32.0/M_PI);
00131       //l=3,m=3
00132       ylmvec[14]->set(l,k,j,i)=-1.0*(pow(xval,3)-3.0*xval*yval*yval)*
00133         sqrt(35.0/64.0/M_PI);
00134       ylmvec[15]->set(l,k,j,i)=-1.0*(3.0*xval*xval*yval-pow(yval,3))*
00135         sqrt(35.0/64.0/M_PI);
00136         }
00137       }
00138       // l=4 included?
00139       if (nylm>16 ) {
00140         if (nylm <25) {abort();} else {
00141       //l=4,m=0
00142       ylmvec[16]->set(l,k,j,i)=(35.0*pow(zval,4)-30.0*zval*zval*rval*rval+3*pow(rval,4))*
00143         sqrt(9.0/256.0/M_PI);
00144       //l=4,m=1
00145       ylmvec[17]->set(l,k,j,i)=-1.0*(7.0*pow(zval,3)-3*zval*rval*rval)*xval*
00146         sqrt(45.0/64.0/M_PI);
00147       ylmvec[18]->set(l,k,j,i)=-1.0*(7.0*pow(zval,3)-3*zval*rval*rval)*yval*
00148         sqrt(45.0/64.0/M_PI);
00149       //l=4,m=2
00150       ylmvec[19]->set(l,k,j,i)=(7.0*zval*zval-rval*rval)*(xval*xval-yval*yval)*
00151         sqrt(45./128.0/M_PI);
00152       ylmvec[20]->set(l,k,j,i)=(7.0*zval*zval-rval*rval)*(2.0*xval*yval)*
00153         sqrt(45./128.0/M_PI);
00154       //l=4,m=3
00155       ylmvec[21]->set(l,k,j,i)=-1.0*zval*(pow(xval,3)-3.0*xval*yval*yval)*
00156         sqrt(315.0/64.0/M_PI);
00157       ylmvec[22]->set(l,k,j,i)=-1.0*zval*(3.0*xval*xval*yval-pow(yval,3))*
00158         sqrt(315.0/64.0/M_PI);
00159       //l=4,m=4
00160       ylmvec[23]->set(l,k,j,i)=(pow(xval,4)-6*xval*xval*yval*yval+pow(yval,4))*
00161         sqrt(315.0/512.0/M_PI);
00162       ylmvec[24]->set(l,k,j,i)=4.0*xval*yval*(xval*xval-yval*yval)*
00163         sqrt(315.0/512.0/M_PI);
00164         }
00165       }
00166       // l=5 included?
00167       if (nylm>25 ) {
00168         if (nylm <36) {abort();} else {
00169       //l=5,m=0
00170       ylmvec[25]->set(l,k,j,i)=(63.0*pow(zval,5)-70.0*pow(zval,3)*rval*rval+15*zval*pow(rval,4))*
00171         sqrt(11.0/256.0/M_PI);
00172       //l=5,m=1
00173       ylmvec[26]->set(l,k,j,i)=-1.0*(21.0*pow(zval,4)-14*zval*zval*rval*rval+pow(rval,4))*xval*
00174         sqrt(165.0/512.0/M_PI);
00175       ylmvec[27]->set(l,k,j,i)=-1.0*(21.0*pow(zval,4)-14*zval*zval*rval*rval+pow(rval,4))*yval*
00176         sqrt(165.0/512.0/M_PI);
00177       //l=5,m=2
00178       ylmvec[28]->set(l,k,j,i)=(3.0*pow(zval,3)-zval*rval*rval)*(xval*xval-yval*yval)*
00179         sqrt(1155./128.0/M_PI);
00180       ylmvec[29]->set(l,k,j,i)=(3.0*pow(zval,3)-zval*rval*rval)*(2.0*xval*yval)*
00181         sqrt(1155./128.0/M_PI);
00182       //l=5,m=3
00183       ylmvec[30]->set(l,k,j,i)=-1.0*(9.0*zval*zval-rval*rval)*(pow(xval,3)-3.0*xval*yval*yval)*
00184         sqrt(385.0/1024.0/M_PI);
00185       ylmvec[31]->set(l,k,j,i)=-1.0*(9.0*zval*zval-rval*rval)*(3.0*xval*xval*yval-pow(yval,3))*
00186         sqrt(385.0/1024.0/M_PI);
00187       //l=5,m=4
00188       ylmvec[32]->set(l,k,j,i)=zval*(pow(xval,4)-6*xval*xval*yval*yval+pow(yval,4))*
00189         sqrt(3465.0/512.0/M_PI);
00190       ylmvec[33]->set(l,k,j,i)=zval*4.0*xval*yval*(xval*xval-yval*yval)*
00191         sqrt(3465.0/512.0/M_PI);
00192       //l=5,m=5
00193       ylmvec[34]->set(l,k,j,i)=-1.0*(pow(xval,5)-10.0*pow(xval,3)*yval*yval+5.0*xval*pow(yval,4))*
00194         sqrt(693.0/1024.0/M_PI);
00195       ylmvec[35]->set(l,k,j,i)=-1.0*(5.0*pow(xval,4)*yval-10.0*xval*xval*pow(yval,3)+pow(yval,5))*
00196         sqrt(693.0/1024.0/M_PI);
00197         }
00198       }
00199       // l=6 included?
00200       if (nylm>36 ) {
00201         if (nylm <49) {abort();} else {
00202       //l=6,m=0
00203       ylmvec[36]->set(l,k,j,i)=(231.0*pow(zval,6)-315.0*pow(zval,4)*rval*rval+105.0*zval*zval*pow(rval,4)-5.0*pow(rval,6))*
00204         sqrt(13.0/1024.0/M_PI);
00205       //l=6,m=1
00206       ylmvec[37]->set(l,k,j,i)=-1.0*(33.0*pow(zval,5)-30.0*pow(zval,3)*rval*rval+5.0*zval*pow(rval,4))*xval*
00207         sqrt(273.0/512.0/M_PI);
00208       ylmvec[38]->set(l,k,j,i)=-1.0*(33.0*pow(zval,5)-30.0*pow(zval,3)*rval*rval+5.0*zval*pow(rval,4))*yval*
00209         sqrt(273.0/512.0/M_PI);
00210       //l=6,m=2
00211       ylmvec[39]->set(l,k,j,i)=(33.0*pow(zval,4)-18.0*zval*zval*rval*rval+pow(rval,4))*(xval*xval-yval*yval)*
00212         sqrt(1365./4096.0/M_PI);
00213       ylmvec[40]->set(l,k,j,i)=(33.0*pow(zval,4)-18.0*zval*zval*rval*rval+pow(rval,4))*(2.0*xval*yval)*
00214         sqrt(1365./4096.0/M_PI);
00215       //l=6,m=3
00216       ylmvec[41]->set(l,k,j,i)=-1.0*(11.0*pow(zval,3)-3.0*zval*rval*rval)*(pow(xval,3)-3.0*xval*yval*yval)*
00217         sqrt(1365.0/1024.0/M_PI);
00218       ylmvec[42]->set(l,k,j,i)=-1.0*(11.0*pow(zval,3)-3.0*zval*rval*rval)*(3.0*xval*xval*yval-pow(yval,3))*
00219         sqrt(1365.0/1024.0/M_PI);
00220       //l=6,m=4
00221       ylmvec[43]->set(l,k,j,i)=(11.0*zval*zval-rval*rval)*(pow(xval,4)-6*xval*xval*yval*yval+pow(yval,4))*
00222         sqrt(819.0/2048.0/M_PI);
00223       ylmvec[44]->set(l,k,j,i)=(11.0*zval*zval-rval*rval)*4.0*xval*yval*(xval*xval-yval*yval)*
00224         sqrt(819.0/2048.0/M_PI);
00225       //l=6,m=5
00226       ylmvec[45]->set(l,k,j,i)=-1.0*zval*(pow(xval,5)-10.0*pow(xval,3)*yval*yval+5.0*xval*pow(yval,4))*
00227         sqrt(9009.0/1024.0/M_PI);
00228       ylmvec[46]->set(l,k,j,i)=-1.0*zval*(5.0*pow(xval,4)*yval-10.0*xval*xval*pow(yval,3)+pow(yval,5))*
00229         sqrt(9009.0/1024.0/M_PI);
00230       //l=6,m=6
00231       ylmvec[47]->set(l,k,j,i)=(pow(xval,6)-15.0*pow(xval,4)*yval*yval+15.0*xval*xval*pow(yval,4)-pow(yval,6))*
00232         sqrt(3003.0/4096.0/M_PI);
00233       ylmvec[48]->set(l,k,j,i)=(6.0*pow(xval,5)*yval-20.0*pow(xval,3)*pow(yval,3)+6.0*xval*pow(yval,5))*
00234         sqrt(3003.0/4096.0/M_PI);
00235         }
00236       }
00237       if(nylm >49) {
00238         cout << "l>6 not implemented!!!!!!!"<< endl;
00239         abort();
00240       }
00241     }
00242       }
00243     }
00244   }
00245 
00246 }
00247 
00248 void Et_bin_bhns_extr::get_integrals(int nylm, double* intvec, Cmp** ylmvec,
00249                      Cmp source) const {
00250 
00251   // As mentioned in the comment before get_ylm, our real/imaginary
00252   // representation of the Y_lm Cmp's does not agree with the normalization 
00253   // used to define
00254   // the spherical harmonic decomposition of Cmp arrays.  Thus, we multiply
00255   // all terms with m>=1 by a factor of 2.0 in order to
00256   // produce the correct decomposition.
00257 
00258   int nz=mp.get_mg()->get_nzone() ;
00259 
00260   Map_af mapping (mp);
00261   
00262   const double* a1 = mapping.get_alpha() ;
00263   const double* b1 = mapping.get_beta() ;
00264   
00265   double rlim=a1[nz-1]+b1[nz-1];
00266   
00267   int ll=0;
00268   int mm=0;
00269   int ncount=0;
00270   for (int n=0; n<nylm; n++) {
00271 
00272     Cmp ylmsource=(*ylmvec[n]*source);
00273     int symcheck=1;
00274     for (int l=0; l<nz; l++) {
00275       int symc=ylmsource.va.base.get_base_t(l);
00276       if(symc!=2304 && symc!=1280)symcheck=0;
00277     }
00278     if(symcheck==1) {
00279       intvec[n]=ylmsource.integrale()/(2.0*ll+1.0)/sqrt(2.0*M_PI)/pow(rlim,ll+1);
00280     } else {
00281       intvec[n]=0;
00282     }
00283     if(mm>=1)intvec[n]*=2.0;
00284 
00285     int lnew=0;
00286     int mnew=0;
00287     int nnew=0;
00288     if(mm<ll) {
00289       if(mm==0) {
00290     lnew=ll;
00291     mnew=mm+1;
00292     nnew=0;
00293       }
00294       if(mm>0&&ncount==0) {
00295     lnew=ll;
00296     mnew=mm;
00297     nnew=1;
00298       }
00299       if(mm>0&&ncount==1) {
00300     lnew=ll;
00301     mnew=mm+1;
00302     nnew=0;
00303       }
00304     }
00305     if(mm==ll) {
00306       if(mm==0) {
00307     lnew=ll+1;
00308     mnew=0;
00309     nnew=0;
00310       }
00311       if(mm>0&&ncount==0) {
00312     lnew=ll;
00313     mnew=mm;
00314     nnew=1;
00315       }
00316       if(mm>0&&ncount==1) {
00317     lnew=ll+1;
00318     mnew=0;
00319     nnew=0;
00320       }
00321     }
00322     ll=lnew;
00323     mm=mnew;
00324     ncount=nnew;
00325   }
00326 }