| 1 | #include <math.h>
|
|---|
| 2 | #include "syn2fast.h"
|
|---|
| 3 | #include "lambuilder.h"
|
|---|
| 4 | #include "fftserver.h"
|
|---|
| 5 |
|
|---|
| 6 | extern "C" {
|
|---|
| 7 | // void fft_gpd_(long double* ,int& ,int& ,int& ,int& ,long double*);
|
|---|
| 8 | float gasdev2_(int& idum);
|
|---|
| 9 | }
|
|---|
| 10 |
|
|---|
| 11 | void a2lm2map(const int nsmax,const int nlmax,const int nmmax,
|
|---|
| 12 | const vector< vector< complex<double> > >&alme,
|
|---|
| 13 | const vector< vector< complex<double> > >&almb,
|
|---|
| 14 | vector<float>& mapq,
|
|---|
| 15 | vector<float>& mapu,
|
|---|
| 16 | vector< complex<long double> > b_northp,
|
|---|
| 17 | vector< complex<long double> > b_southp,
|
|---|
| 18 | vector< complex<long double> > bwp,
|
|---|
| 19 | vector< complex<long double> > b_northm,
|
|---|
| 20 | vector< complex<long double> > b_southm,
|
|---|
| 21 | vector< complex<long double> > bwm){
|
|---|
| 22 | /*=======================================================================
|
|---|
| 23 | computes a map form its alm for the HEALPIX pixelisation
|
|---|
| 24 | map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
|
|---|
| 25 | = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
|
|---|
| 26 |
|
|---|
| 27 | where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
|
|---|
| 28 |
|
|---|
| 29 | * the recurrence of Ylm is the standard one (cf Num Rec)
|
|---|
| 30 | * the sum over m is done by FFT
|
|---|
| 31 |
|
|---|
| 32 | =======================================================================*/
|
|---|
| 33 | vector< complex<float> > mapp;
|
|---|
| 34 | vector< complex<float> > mapm;
|
|---|
| 35 |
|
|---|
| 36 | //convert the alm from e/b to +/-
|
|---|
| 37 | complex<double> im(0,1);
|
|---|
| 38 | vector< vector< complex<double> > >almp=alme;
|
|---|
| 39 | vector< vector< complex<double> > >almm=almb;
|
|---|
| 40 | for (int i=0;i<(signed) almp.size();i++){
|
|---|
| 41 | for (int j=0;j<(signed) almp[i].size();j++){
|
|---|
| 42 | almp[i][j]=-(alme[i][j]+im*almb[i][j]);
|
|---|
| 43 | almm[i][j]=-(alme[i][j]-im*almb[i][j]);
|
|---|
| 44 | }
|
|---|
| 45 | }
|
|---|
| 46 | //for (int i=0;i<(signed) almp.size();i++){cout<<alme[i][0]<<" "<<almb[i][0]<<" "<<almp[i][0]<<" "<<almm[i][0]<<endl;}
|
|---|
| 47 | mapp.resize(12*nsmax*nsmax);
|
|---|
| 48 | b_northp.resize(2*nmmax+1); //index m corresponds to nmmax+m
|
|---|
| 49 | b_southp.resize(2*nmmax+1);
|
|---|
| 50 | bwp.resize(4*nsmax);
|
|---|
| 51 | mapm.resize(12*nsmax*nsmax);
|
|---|
| 52 | b_northm.resize(2*nmmax+1); //index m corresponds to nmmax+m
|
|---|
| 53 | b_southm.resize(2*nmmax+1);
|
|---|
| 54 | bwm.resize(4*nsmax);
|
|---|
| 55 |
|
|---|
| 56 |
|
|---|
| 57 | /* INTEGER l, indl */
|
|---|
| 58 |
|
|---|
| 59 | int istart_north = 0;
|
|---|
| 60 | int istart_south = 12*nsmax*nsmax;
|
|---|
| 61 |
|
|---|
| 62 | double dth1 = 1. / (3.*nsmax*nsmax);
|
|---|
| 63 | double dth2 = 2. / (3.*nsmax);
|
|---|
| 64 | double dst1 = 1. / (sqrt(6.) * nsmax);
|
|---|
| 65 |
|
|---|
| 66 | for (int ith = 1; ith <= 2*nsmax;ith++){
|
|---|
| 67 | int nph, kphi0;
|
|---|
| 68 | double cth, sth, sth2;
|
|---|
| 69 | if (ith <= nsmax-1){ /* north polar cap */
|
|---|
| 70 | nph = 4*ith;
|
|---|
| 71 | kphi0 = 1;
|
|---|
| 72 | cth = 1. - dth1*ith*ith; /* cos(theta) */
|
|---|
| 73 | sth = sin( 2. * asin( ith * dst1 ) ) ; /* sin(theta) */
|
|---|
| 74 | sth2 = sth*sth;
|
|---|
| 75 | } else { /* tropical band + equat. */
|
|---|
| 76 | nph = 4*nsmax;
|
|---|
| 77 | kphi0 = (ith+1-nsmax) % 2;
|
|---|
| 78 | cth = (2.*nsmax-ith) * dth2;
|
|---|
| 79 | sth = sqrt((1.-cth)*(1.+cth)); /* ! sin(theta)*/
|
|---|
| 80 | sth2=(1.-cth)*(1.+cth);
|
|---|
| 81 | }
|
|---|
| 82 |
|
|---|
| 83 | /* -----------------------------------------------------
|
|---|
| 84 | for each theta, and each m, computes
|
|---|
| 85 | b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
|
|---|
| 86 | ------------------------------------------------------
|
|---|
| 87 | lambda_mm tends to go down when m increases (risk of underflow)
|
|---|
| 88 | lambda_lm tends to go up when l increases (risk of overflow)*/
|
|---|
| 89 | Lambda2Builder l2b(acos(cth),nlmax,nmmax);
|
|---|
| 90 | for (int m = 0; m <= nmmax; m++){
|
|---|
| 91 | complex<double> b_np,b_sp,b_nm,b_sm;
|
|---|
| 92 | // cout <<l2b.lam2lmp(m,m)<<endl;
|
|---|
| 93 | b_np = l2b.lam2lmp(m,m)* almp[m][m];
|
|---|
| 94 | b_sp = l2b.lam2lmp(m,m,-1)* almp[m][m];
|
|---|
| 95 | b_nm = l2b.lam2lmm(m,m)* almm[m][m];
|
|---|
| 96 | b_sm = l2b.lam2lmm(m,m,-1)* almm[m][m];
|
|---|
| 97 | for (int l = m+1; l<= nlmax; l++){
|
|---|
| 98 | //cout<<l2b.lam2lmp(l,m)<<endl;
|
|---|
| 99 | b_np += l2b.lam2lmp(l,m)*almp[l][m];
|
|---|
| 100 | b_sp += l2b.lam2lmp(l,m,-1)*almp[l][m];
|
|---|
| 101 | b_nm += l2b.lam2lmm(l,m)*almm[l][m];
|
|---|
| 102 | b_sm += l2b.lam2lmm(l,m,-1)*almm[l][m];
|
|---|
| 103 | }
|
|---|
| 104 |
|
|---|
| 105 | b_northp[m+nmmax] = b_np;
|
|---|
| 106 | b_southp[m+nmmax] = b_sp;
|
|---|
| 107 | b_northm[m+nmmax] = b_nm;
|
|---|
| 108 | b_southm[m+nmmax] = b_sm;
|
|---|
| 109 | }
|
|---|
| 110 |
|
|---|
| 111 | // obtains the negative m of b(m,theta) (= complex conjugate)
|
|---|
| 112 | for (int m=-nmmax;m<=-1;m++){
|
|---|
| 113 | int fac = 1;//(int) pow(-1,m);// ! either 1 or -1
|
|---|
| 114 | complex<long double>
|
|---|
| 115 | shit(b_northp[-m+nmmax].real(),-b_northp[-m+nmmax].imag());
|
|---|
| 116 | complex<long double>
|
|---|
| 117 | shit2(b_southp[-m+nmmax].real(),-b_southp[-m+nmmax].imag());
|
|---|
| 118 | shit*=fac; shit2*=fac;
|
|---|
| 119 | b_northm[m+nmmax] = shit;
|
|---|
| 120 | b_southm[m+nmmax] = shit2;
|
|---|
| 121 | shit=complex<long double>(b_northm[-m+nmmax].real(),-b_northm[-m+nmmax].imag());
|
|---|
| 122 | shit2=complex<long double>(b_southm[-m+nmmax].real(),-b_southm[-m+nmmax].imag());
|
|---|
| 123 | shit*=fac; shit2*=fac;
|
|---|
| 124 | b_northp[m+nmmax] = shit;
|
|---|
| 125 | b_southp[m+nmmax] = shit2;
|
|---|
| 126 | }
|
|---|
| 127 |
|
|---|
| 128 | // for (int i=0;i<2*nmmax+1;i++){ cout<< b_northp[i]<<" "<<b_southp[i]<<" "<<b_northm[i]
|
|---|
| 129 | // <<" "<<b_southm[i]<<endl;}
|
|---|
| 130 | /* ---------------------------------------------------------------
|
|---|
| 131 | sum_m b(m,theta)*exp(i*m*phi) -> f(phi,theta)
|
|---|
| 132 | ---------------------------------------------------------------*/
|
|---|
| 133 | syn2_phas(nsmax,nlmax,nmmax,b_northp,nph,mapp,istart_north,kphi0,bwp); // north hemisph. + equator
|
|---|
| 134 | syn2_phas(nsmax,nlmax,nmmax,b_northm,nph,mapm,istart_north,kphi0,bwm);
|
|---|
| 135 | istart_north = istart_north + nph;
|
|---|
| 136 | if (ith < 2*nsmax){
|
|---|
| 137 | istart_south=istart_south-nph;
|
|---|
| 138 | syn2_phas(nsmax,nlmax,nmmax,b_southp,nph,mapp,istart_south,kphi0,bwp); // south hemisph. w/o equat
|
|---|
| 139 | syn2_phas(nsmax,nlmax,nmmax,b_southm,nph,mapm,istart_south,kphi0,bwm);
|
|---|
| 140 | }
|
|---|
| 141 | }
|
|---|
| 142 | for (int i=0;i< (signed)mapp.size();i++){
|
|---|
| 143 | //cout << mapp[i]<<" "<<mapm[i]<<endl;
|
|---|
| 144 | mapq[i]=(mapp[i]+mapm[i]).real()/2;
|
|---|
| 145 | mapu[i]=(mapp[i]-mapm[i]).imag()/2;
|
|---|
| 146 | //cout << mapq[i]<<" "<<mapu[i]<<endl;
|
|---|
| 147 | }
|
|---|
| 148 | }
|
|---|
| 149 |
|
|---|
| 150 | void syn2_phas(const int nsmax,const int nlmax,const int nmmax,
|
|---|
| 151 | const vector< complex<long double> >& datain,
|
|---|
| 152 | int nph,vector< complex<float> >& dataout, const int start,
|
|---|
| 153 | int kphi0,
|
|---|
| 154 | vector< complex<long double> >& bw){
|
|---|
| 155 |
|
|---|
| 156 | /*=======================================================================
|
|---|
| 157 | dataout(j) = sum_m datain(m) * exp(i*m*phi(j))
|
|---|
| 158 | with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1
|
|---|
| 159 |
|
|---|
| 160 | as the set of frequencies {m} is larger than nph,
|
|---|
| 161 | we wrap frequencies within {0..nph-1}
|
|---|
| 162 | ie m = k*nph + m' with m' in {0..nph-1}
|
|---|
| 163 | then
|
|---|
| 164 | noting bw(m') = exp(i*m'*phi0)
|
|---|
| 165 | * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
|
|---|
| 166 | with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
|
|---|
| 167 | dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
|
|---|
| 168 | = Fourier Transform of bw
|
|---|
| 169 | is real
|
|---|
| 170 |
|
|---|
| 171 | NB nph is not necessarily a power of 2
|
|---|
| 172 |
|
|---|
| 173 | =======================================================================*/
|
|---|
| 174 | int ksign = 1;
|
|---|
| 175 | complex<double>* data= new complex<double>[4*nsmax];
|
|---|
| 176 |
|
|---|
| 177 | for (int iw=0;iw<nph;iw++){bw[iw]=0;}
|
|---|
| 178 |
|
|---|
| 179 | double phi0 = kphi0*M_PI/nph;
|
|---|
| 180 | int kshift = (int) pow(-1,kphi0);// ! either 1 or -1
|
|---|
| 181 | // all frequencies [-m,m] are wrapped in [0,nph-1]
|
|---|
| 182 | for (int i=1;i<=2*nmmax+1;i++){
|
|---|
| 183 | int m=i-nmmax-1; //in -nmmax, nmmax
|
|---|
| 184 | int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
|
|---|
| 185 | int k=(m-iw)/nph; //number of 'turns'
|
|---|
| 186 | complex<long double> shit(pow(kshift,k));
|
|---|
| 187 | bw[iw]+=datain[i-1]*shit;//complex number
|
|---|
| 188 | }
|
|---|
| 189 |
|
|---|
| 190 | // kshift**k = 1 for even turn numbers
|
|---|
| 191 | // = 1 or -1 for odd turn numbers : results from the shift in space
|
|---|
| 192 |
|
|---|
| 193 | // applies the shift in position <-> phase factor in Fourier space
|
|---|
| 194 | for (int iw=1;iw<=nph;iw++){
|
|---|
| 195 | int m=ksign*(iw-1);
|
|---|
| 196 | complex<long double> shit(cos(m*phi0),sin(m*phi0));
|
|---|
| 197 | data[iw-1]=bw[iw-1]*shit;
|
|---|
| 198 | }
|
|---|
| 199 | for (int i = nph; i< 4*nsmax;i++){
|
|---|
| 200 | data[i] = 0;
|
|---|
| 201 | }
|
|---|
| 202 | //cout<<endl;
|
|---|
| 203 | //fft_gpd_(data,nph,dum0,ksign,dum1,work); // complex to complex
|
|---|
| 204 | FFTServer fft;
|
|---|
| 205 | fft.fftf(nph,data);
|
|---|
| 206 |
|
|---|
| 207 | for (int iw=0;iw<nph;iw++){
|
|---|
| 208 | dataout[iw+start]=data[iw];
|
|---|
| 209 | }
|
|---|
| 210 | delete[] data;
|
|---|
| 211 | }
|
|---|
| 212 |
|
|---|
| 213 |
|
|---|
| 214 | void create_a2lm(const int nsmax,const int nlmax,const int nmmax,
|
|---|
| 215 | vector< vector< complex<double> > >& a2lme,
|
|---|
| 216 | vector< vector< complex<double> > >& a2lmb,
|
|---|
| 217 | vector<float>& cls_e,
|
|---|
| 218 | vector<float>& cls_b,
|
|---|
| 219 | int& iseed,const float fwhm,
|
|---|
| 220 | vector<float> lread){
|
|---|
| 221 |
|
|---|
| 222 |
|
|---|
| 223 | /*=======================================================================
|
|---|
| 224 | creates the a_lm from the power spectrum,
|
|---|
| 225 | assuming they are gaussian and complex
|
|---|
| 226 | with a variance given by C(l)
|
|---|
| 227 |
|
|---|
| 228 | the input file should contain : l and C(l) with *consecutive* l's
|
|---|
| 229 | (missing C(l) are put to 0.)
|
|---|
| 230 |
|
|---|
| 231 | because the map is real we have : a_l-m = (-)^m conjug(a_lm)
|
|---|
| 232 | so we actually compute them only for m >= 0
|
|---|
| 233 |
|
|---|
| 234 | modifie G. Le Meur (13/01/98) :
|
|---|
| 235 | on ne lit plus les C(l) sur un fichier. Le tableau est entre en argument
|
|---|
| 236 | (cls_tt). Ce tableau doit contenir les valeur de C(l) par ordre
|
|---|
| 237 | SEQUENTIEL de l (de l=0 a l=nlmax)
|
|---|
| 238 | =======================================================================*/
|
|---|
| 239 | /* CHARACTER*128 filename
|
|---|
| 240 |
|
|---|
| 241 | integer unit,n_l,j,il,im
|
|---|
| 242 | real fs,quadrupole,xxxtmp,correct
|
|---|
| 243 | real fs_tt
|
|---|
| 244 | real zeta1_r, zeta1_i
|
|---|
| 245 | LOGICAL ok
|
|---|
| 246 | character*20 string_quad
|
|---|
| 247 |
|
|---|
| 248 | real*4 prefact(0:2)
|
|---|
| 249 | logical prefact_bool
|
|---|
| 250 |
|
|---|
| 251 | integer lneffct
|
|---|
| 252 | real gasdev2
|
|---|
| 253 | external lneffct, gasdev2
|
|---|
| 254 |
|
|---|
| 255 | c-----------------------------------------------------------------------*/
|
|---|
| 256 | a2lme.resize(nlmax+1);
|
|---|
| 257 | a2lmb.resize(nlmax+1);
|
|---|
| 258 | for (int i=0; i< (signed) a2lme.size();i++){
|
|---|
| 259 | a2lme[i].resize(nmmax+1);
|
|---|
| 260 | a2lmb[i].resize(nmmax+1);
|
|---|
| 261 | }
|
|---|
| 262 |
|
|---|
| 263 |
|
|---|
| 264 | iseed = -abs(iseed);
|
|---|
| 265 | if (iseed == 0){iseed=-1;}
|
|---|
| 266 | int idum = iseed;
|
|---|
| 267 |
|
|---|
| 268 | float sig_smooth = fwhm/sqrt(8.*log(2.))/(60.*180.)* M_PI;
|
|---|
| 269 |
|
|---|
| 270 | for (int i=0;i<nlmax+1;i++){lread[i]=0;}
|
|---|
| 271 |
|
|---|
| 272 | for (int i=0;i <= nlmax;i++){lread[i] = i;}
|
|---|
| 273 |
|
|---|
| 274 | int n_l = nlmax+1;
|
|---|
| 275 |
|
|---|
| 276 | cout<<lread[0]<<" <= l <= "<<lread[n_l-1]<<endl;
|
|---|
| 277 |
|
|---|
| 278 | // --- smoothes the initial power spectrum ---
|
|---|
| 279 |
|
|---|
| 280 | for (int i=0;i<n_l;i++){
|
|---|
| 281 | int l= (int) lread[i];
|
|---|
| 282 | float gauss=exp(-l*(l+1.)*sig_smooth*sig_smooth);
|
|---|
| 283 | cls_e[i]*=gauss;
|
|---|
| 284 | cls_b[i]*=gauss;
|
|---|
| 285 | }
|
|---|
| 286 |
|
|---|
| 287 |
|
|---|
| 288 | // --- generates randomly the alm according to their power spectrum ---
|
|---|
| 289 | float hsqrt2 = 1.0 / sqrt(2.0);
|
|---|
| 290 |
|
|---|
| 291 | for (int i=0;i<n_l;i++){
|
|---|
| 292 | int l=(int) lread[i];
|
|---|
| 293 | float rms_e=sqrt(cls_e[i]);
|
|---|
| 294 | float rms_b=sqrt(cls_b[i]);
|
|---|
| 295 | // ------ m = 0 ------
|
|---|
| 296 | complex<float> zeta1(gasdev2_(idum));
|
|---|
| 297 | a2lme[l][0] = zeta1 * rms_e;
|
|---|
| 298 | zeta1=complex<float>(gasdev2_(idum));
|
|---|
| 299 | a2lmb[l][0] = zeta1 * rms_b;
|
|---|
| 300 |
|
|---|
| 301 | //------ m > 0 ------
|
|---|
| 302 | for (int m=1;m<=l;m++){
|
|---|
| 303 | complex<float> shit1(hsqrt2);
|
|---|
| 304 | complex<float> shit2(gasdev2_(idum),gasdev2_(idum));
|
|---|
| 305 | zeta1=shit1*shit2;
|
|---|
| 306 | a2lme[l][m]=rms_e*zeta1;
|
|---|
| 307 | shit2= complex<float>(gasdev2_(idum),gasdev2_(idum));
|
|---|
| 308 | zeta1=shit1*shit2;
|
|---|
| 309 | a2lmb[l][m]=rms_b*zeta1;
|
|---|
| 310 | }
|
|---|
| 311 | }
|
|---|
| 312 | }
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
