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