| 1 | #include <math.h>
|
|---|
| 2 | #include <vector>
|
|---|
| 3 | #include <fftserver.h>
|
|---|
| 4 | #include <complex>
|
|---|
| 5 | #include "ana2fast.h"
|
|---|
| 6 | #include "lambuilder.h"
|
|---|
| 7 |
|
|---|
| 8 | /*extern "C" {
|
|---|
| 9 | void fft_gpd_(long double* ,int& ,int& ,int& ,int& ,long double*);
|
|---|
| 10 | }*/
|
|---|
| 11 |
|
|---|
| 12 | void map2a2lm(int nsmax,int nlmax,int nmmax,const vector<float>& mapq,
|
|---|
| 13 | const vector<float>& mapu,
|
|---|
| 14 | vector< vector< complex<double> > >& a2lme,
|
|---|
| 15 | vector< vector< complex<double> > >& a2lmb,
|
|---|
| 16 | double cos_theta_cut){
|
|---|
| 17 |
|
|---|
| 18 | // REAL*4 powspec(0:nlmax)
|
|---|
| 19 |
|
|---|
| 20 | // integer npmiss,npmt,id_miss(10000)
|
|---|
| 21 |
|
|---|
| 22 | //create the maps for which there are nice basis functions
|
|---|
| 23 |
|
|---|
| 24 | vector< complex<float> > mapp(mapq.size());
|
|---|
| 25 | vector< complex<float> > mapm(mapq.size());
|
|---|
| 26 | for (int i=0;i< (signed) mapq.size();i++){
|
|---|
| 27 | mapp[i]=complex<float>(mapq[i],mapu[i]);
|
|---|
| 28 | mapm[i]=complex<float>(mapq[i],-mapu[i]);
|
|---|
| 29 | //cout <<"the maps"<< mapp[i]<<" "<<mapm[i]<<endl;
|
|---|
| 30 | }
|
|---|
| 31 |
|
|---|
| 32 | vector< vector< complex<double> > > a2lmp;
|
|---|
| 33 | vector< vector< complex<double> > > a2lmm;
|
|---|
| 34 | a2lmp.resize(nlmax+1);
|
|---|
| 35 | for (int i=0; i< (signed) a2lmp.size();i++){
|
|---|
| 36 | a2lmp[i].resize(nmmax+1);
|
|---|
| 37 | for (int j=0; j< (signed) a2lmp[i].size();j++)a2lmp[i][j]=0;
|
|---|
| 38 | }
|
|---|
| 39 | a2lmm.resize(nlmax+1);
|
|---|
| 40 | for (int i=0; i< (signed) a2lmm.size();i++){
|
|---|
| 41 | a2lmm[i].resize(nmmax+1);
|
|---|
| 42 | for (int j=0; j< (signed) a2lmm[i].size();j++)a2lmm[i][j]=0;
|
|---|
| 43 | }
|
|---|
| 44 |
|
|---|
| 45 | /*-----------------------------------------------------------------------
|
|---|
| 46 | computes the integral in phi : phas_m(theta)
|
|---|
| 47 | for each parallele from north to south pole
|
|---|
| 48 | -----------------------------------------------------------------------*/
|
|---|
| 49 |
|
|---|
| 50 | int istart_north = 0;
|
|---|
| 51 | int istart_south = 12*nsmax*nsmax;
|
|---|
| 52 |
|
|---|
| 53 | double dth1 = 1. / (3.*nsmax*nsmax);
|
|---|
| 54 | double dth2 = 2. / (3.*nsmax);
|
|---|
| 55 | double dst1 = 1. / (sqrt(6.) * nsmax);
|
|---|
| 56 |
|
|---|
| 57 | vector< complex<double> > phas_np(nmmax+1), phas_sp(nmmax+1),
|
|---|
| 58 | phas_nm(nmmax+1),phas_sm(nmmax+1);
|
|---|
| 59 |
|
|---|
| 60 | for (int ith = 1; ith <= 2*nsmax;ith++){
|
|---|
| 61 | int nph, kphi0;
|
|---|
| 62 | double cth, sth, sth2;
|
|---|
| 63 | //assign doesn't seem to exist in our compiler
|
|---|
| 64 | //phas_n.assign(nmmax+1,(complex<float>) 0);
|
|---|
| 65 | //phas_s.assign(nmmax+1,(complex<float>) 0);
|
|---|
| 66 | for (int i=0;i< nmmax+1;i++){
|
|---|
| 67 | phas_np[i]=0; phas_sp[i]=0;phas_nm[i]=0;phas_sm[i]=0;
|
|---|
| 68 | }
|
|---|
| 69 |
|
|---|
| 70 | if (ith <= nsmax-1){ /* north polar cap */
|
|---|
| 71 | nph = 4*ith;
|
|---|
| 72 | kphi0 = 1;
|
|---|
| 73 | cth = 1. - dth1*ith*ith; /* cos(theta) */
|
|---|
| 74 | sth = sin( 2. * asin( ith * dst1 ) ) ; /* sin(theta) */
|
|---|
| 75 | sth2 = sth*sth;
|
|---|
| 76 | } else { /* tropical band + equat. */
|
|---|
| 77 | nph = 4*nsmax;
|
|---|
| 78 | kphi0 = (ith+1-nsmax) % 2;
|
|---|
| 79 | cth = (2.*nsmax-ith) * dth2;
|
|---|
| 80 | sth = sqrt((1.-cth)*(1.+cth)); /* ! sin(theta)*/
|
|---|
| 81 | sth2=(1.-cth)*(1.+cth);
|
|---|
| 82 | }
|
|---|
| 83 |
|
|---|
| 84 | //part of the sky out of the symetric cut
|
|---|
| 85 | bool keep_it = (abs(cth) >= cos_theta_cut);
|
|---|
| 86 |
|
|---|
| 87 | //make sure that map is well defined
|
|---|
| 88 | if (keep_it){
|
|---|
| 89 | comp_phas2_2(nsmax,nlmax,nmmax,mapp,mapm,istart_north,nph,phas_np,
|
|---|
| 90 | phas_nm,kphi0);
|
|---|
| 91 | }
|
|---|
| 92 | istart_north = istart_north + nph;
|
|---|
| 93 |
|
|---|
| 94 | istart_south = istart_south - nph;
|
|---|
| 95 | if (ith < 2*nsmax && keep_it){
|
|---|
| 96 | comp_phas2_2(nsmax,nlmax,nmmax,mapp,mapm,istart_south,nph,phas_sp,
|
|---|
| 97 | phas_sm,kphi0);
|
|---|
| 98 | }
|
|---|
| 99 | /*-----------------------------------------------------------------------
|
|---|
| 100 | computes the a_lm by integrating over theta
|
|---|
| 101 | lambda_lm(theta) * phas_m(theta)
|
|---|
| 102 | for each m and l
|
|---|
| 103 | -----------------------------------------------------------------------*/
|
|---|
| 104 | Lambda2Builder l2b(acos(cth),nlmax,nmmax);
|
|---|
| 105 | //cout << "fft:"<<phas_np[0]<<" "<<phas_sp[0]<<" "<<phas_nm[0]<<" "<<phas_sm[0]<<endl;
|
|---|
| 106 | for (int m = 0; m <= nmmax; m++){
|
|---|
| 107 | cout << phas_np[m]<<" "<<phas_sp[m]<<" "<<phas_nm[m]<<" "<<phas_sm[m]<<endl;
|
|---|
| 108 | a2lmp[m][m]+=l2b.lam2lmp(m,m)*phas_np[m]+l2b.lam2lmp(m,m,-1)*phas_sp[m];
|
|---|
| 109 | a2lmm[m][m]+=l2b.lam2lmm(m,m)*phas_nm[m]+l2b.lam2lmm(m,m,-1)*phas_sm[m];
|
|---|
| 110 | for (int l = m+1; l<= nlmax; l++){
|
|---|
| 111 | a2lmp[l][m]+=
|
|---|
| 112 | l2b.lam2lmp(l,m)*phas_np[m]+l2b.lam2lmp(l,m,-1)*phas_sp[m];
|
|---|
| 113 | a2lmm[l][m]+=
|
|---|
| 114 | l2b.lam2lmm(l,m)*phas_nm[m]+l2b.lam2lmm(l,m,-1)*phas_sm[m];
|
|---|
| 115 | }
|
|---|
| 116 | }
|
|---|
| 117 | }
|
|---|
| 118 | complex<double> im(0,1);
|
|---|
| 119 | a2lme.resize(nlmax+1);
|
|---|
| 120 | for (int i=0; i< (signed) a2lme.size();i++){
|
|---|
| 121 | a2lme[i].resize(nmmax+1);
|
|---|
| 122 | }
|
|---|
| 123 | a2lmb.resize(nlmax+1);
|
|---|
| 124 | for (int i=0; i< (signed) a2lmb.size();i++){
|
|---|
| 125 | a2lmb[i].resize(nmmax+1);
|
|---|
| 126 | }
|
|---|
| 127 | float domega=(4.*M_PI)/(12.*nsmax*nsmax);
|
|---|
| 128 | for (int m = 0; m <= nmmax; m++){
|
|---|
| 129 | a2lme[m][m]=-(a2lmp[m][m]+a2lmm[m][m])/2.*static_cast<double>(domega);
|
|---|
| 130 | a2lmb[m][m]=im*(a2lmp[m][m]-a2lmm[m][m])/2.*static_cast<double>(domega);
|
|---|
| 131 | for (int l = m+1; l<= nlmax; l++){
|
|---|
| 132 | a2lme[l][m]=-(a2lmp[l][m]+a2lmm[l][m])/2.*static_cast<double>(domega);
|
|---|
| 133 | a2lmb[l][m]=im*(a2lmp[l][m]-a2lmm[l][m])/2.*static_cast<double>(domega);
|
|---|
| 134 | }
|
|---|
| 135 | }
|
|---|
| 136 | //for (int l = 2; l<= nlmax; l++){
|
|---|
| 137 | //cout << "calc almp,m"<<a2lmp[l][0]<<" "<<a2lmm[l][0]<<endl;}
|
|---|
| 138 | }
|
|---|
| 139 |
|
|---|
| 140 | void comp_phas2_2(int nsmax,int nlmax,int nmmax,
|
|---|
| 141 | const vector< complex<float> >& datain,
|
|---|
| 142 | const vector< complex<float> >& datain2,
|
|---|
| 143 | int start,int nph,vector< complex<double> >& dataout,
|
|---|
| 144 | vector< complex<double> >& dataout2, int kphi0){
|
|---|
| 145 | /*=======================================================================
|
|---|
| 146 | integrates (data * phi-dependence-of-Ylm) over phi
|
|---|
| 147 | --> function of m can be computed by FFT
|
|---|
| 148 | with 0<= m <= npoints/2 (: Nyquist)
|
|---|
| 149 | because the data is real the negative m are the conjugate of the
|
|---|
| 150 | positive ones
|
|---|
| 151 |
|
|---|
| 152 | arguments d'appels : GLM
|
|---|
| 153 | =======================================================================*/
|
|---|
| 154 |
|
|---|
| 155 | int ksign = -1;
|
|---|
| 156 | double phi0 = kphi0*M_PI/nph;
|
|---|
| 157 |
|
|---|
| 158 | complex<double>* data= new complex<double>[4*nsmax];
|
|---|
| 159 | complex<double>* data2= new complex<double>[4*nsmax];
|
|---|
| 160 | for (int i = 0; i< nph;i++){
|
|---|
| 161 | data[i] = datain[i+start];
|
|---|
| 162 | data2[i] = datain2[i+start];
|
|---|
| 163 | }
|
|---|
| 164 | for (int i = nph; i< 4*nsmax;i++){
|
|---|
| 165 | data[i] = 0;
|
|---|
| 166 | data2[i] = 0;
|
|---|
| 167 | }
|
|---|
| 168 |
|
|---|
| 169 | FFTServer fft;
|
|---|
| 170 | fft.fftb(nph,data);
|
|---|
| 171 | fft.fftb(nph,data2);
|
|---|
| 172 |
|
|---|
| 173 | //in the output the frequencies are respectively 0,1,2,..,nph/2,-nph/2+1,..,-2,-1
|
|---|
| 174 | // only the first nph/2+1 (positive freq.) are interesting
|
|---|
| 175 | int im_max = min(nph/2,nmmax);
|
|---|
| 176 | dataout.resize(nmmax+1);
|
|---|
| 177 | dataout2.resize(nmmax+1);
|
|---|
| 178 | for (int i = 1;i <= im_max + 1;i++){
|
|---|
| 179 | int m = ksign*(i-1);
|
|---|
| 180 | complex<double> fuck(cos(m*phi0),sin(m*phi0));
|
|---|
| 181 | dataout[i-1]=data[i-1]*fuck;
|
|---|
| 182 | dataout2[i-1]=data2[i-1]*fuck;
|
|---|
| 183 | }
|
|---|
| 184 | for (int i = im_max + 2;i <= nmmax + 1;i++){
|
|---|
| 185 | dataout[i-1] = 0; dataout2[i-1]=0;
|
|---|
| 186 | }
|
|---|
| 187 | delete[] data;
|
|---|
| 188 | delete[] data2;
|
|---|
| 189 | }
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
