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