| 1 | #include <math.h>
|
|---|
| 2 | #include "syngen.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 alm2mapgen(const int nsmax,const int nlmax,const int nmmax,
|
|---|
| 12 | const vector< vector< complex<double> > >&alm,
|
|---|
| 13 | SphericalMap& map,
|
|---|
| 14 | vector< complex<long double> > b_north,
|
|---|
| 15 | vector< complex<long double> > b_south,
|
|---|
| 16 | vector< complex<long double> > bw,
|
|---|
| 17 | vector< complex<long double> > data,
|
|---|
| 18 | vector< complex<long double> > work){
|
|---|
| 19 | /*=======================================================================
|
|---|
| 20 | computes a map form its alm for the HEALPIX pixelisation
|
|---|
| 21 | map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
|
|---|
| 22 | = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
|
|---|
| 23 |
|
|---|
| 24 | where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
|
|---|
| 25 |
|
|---|
| 26 | * the recurrence of Ylm is the standard one (cf Num Rec)
|
|---|
| 27 | * the sum over m is done by FFT
|
|---|
| 28 |
|
|---|
| 29 | =======================================================================*/
|
|---|
| 30 |
|
|---|
| 31 | // map.resize(12*nsmax*nsmax);
|
|---|
| 32 | b_north.resize(2*nmmax+1); //index m corresponds to nmmax+m
|
|---|
| 33 | b_south.resize(2*nmmax+1);
|
|---|
| 34 | bw.resize(4*nsmax);
|
|---|
| 35 | data.resize(4*nsmax);
|
|---|
| 36 | work.resize(4*nsmax);
|
|---|
| 37 |
|
|---|
| 38 |
|
|---|
| 39 | /* INTEGER l, indl */
|
|---|
| 40 |
|
|---|
| 41 | for (int ith = 0; ith <= (map.NbThetaSlices()+1)/2-1;ith++){
|
|---|
| 42 | int nph;
|
|---|
| 43 | double phi0;
|
|---|
| 44 | float theta,thetas;
|
|---|
| 45 | OVector phin, datan,phis,datas;
|
|---|
| 46 | map.GetThetaSlice(map.NbThetaSlices()-ith-1,thetas,phis,datas);
|
|---|
| 47 | map.GetThetaSlice(ith,theta,phin,datan);
|
|---|
| 48 | nph = phin.NElts();
|
|---|
| 49 | phi0 = phin(0);
|
|---|
| 50 | double cth = cos(theta);
|
|---|
| 51 |
|
|---|
| 52 | /* -----------------------------------------------------
|
|---|
| 53 | for each theta, and each m, computes
|
|---|
| 54 | b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
|
|---|
| 55 | ------------------------------------------------------
|
|---|
| 56 | lambda_mm tends to go down when m increases (risk of underflow)
|
|---|
| 57 | lambda_lm tends to go up when l increases (risk of overflow)*/
|
|---|
| 58 | // lambda_0_0 * big number --> to avoid underflow
|
|---|
| 59 | LambdaBuilder lb(acos(cth),nlmax,nmmax);
|
|---|
| 60 | for (int m = 0; m <= nmmax; m++){
|
|---|
| 61 | complex<double> b_n = lb.lamlm(m,m) * alm[m][m];
|
|---|
| 62 | complex<double> b_s = lb.lamlm(m,m) * alm[m][m]; //m,m even
|
|---|
| 63 | for (int l = m+1; l<= nlmax; l++){
|
|---|
| 64 | b_n += lb.lamlm(l,m)*alm[l][m];
|
|---|
| 65 | b_s += lb.lamlm(l,m,-1)*alm[l][m];
|
|---|
| 66 | }
|
|---|
| 67 | b_north[m+nmmax] = b_n;
|
|---|
| 68 | b_south[m+nmmax] = b_s;
|
|---|
| 69 | }
|
|---|
| 70 |
|
|---|
| 71 | // obtains the negative m of b(m,theta) (= complex conjugate)
|
|---|
| 72 |
|
|---|
| 73 | for (int m=-nmmax;m<=-1;m++){
|
|---|
| 74 | //compiler doesn't have conj()
|
|---|
| 75 | complex<long double>
|
|---|
| 76 | shit(b_north[-m+nmmax].real(),-b_north[-m+nmmax].imag());
|
|---|
| 77 | complex<long double>
|
|---|
| 78 | shit2(b_south[-m+nmmax].real(),-b_south[-m+nmmax].imag());
|
|---|
| 79 | b_north[m+nmmax] = shit;
|
|---|
| 80 | b_south[m+nmmax] = shit2;
|
|---|
| 81 | }
|
|---|
| 82 | /* ---------------------------------------------------------------
|
|---|
| 83 | sum_m b(m,theta)*exp(i*m*phi) -> f(phi,theta)
|
|---|
| 84 | ---------------------------------------------------------------*/
|
|---|
| 85 | /*cout <<"gen ";
|
|---|
| 86 | for (int i=0;i<2*nmmax+1;i++){
|
|---|
| 87 | cout << b_north[i]<<" ";}
|
|---|
| 88 | cout<<endl;*/
|
|---|
| 89 |
|
|---|
| 90 | syn_phasgen(nsmax,nlmax,nmmax,b_north,nph,datan,phi0,bw); // north hemisph. + equator
|
|---|
| 91 | if (ith != map.NbThetaSlices()/2){
|
|---|
| 92 | syn_phasgen(nsmax,nlmax,nmmax,b_south,nph,datas,phi0,bw); // south hemisph. w/o equat
|
|---|
| 93 | for (int i=0;i< nph;i++){
|
|---|
| 94 | map.PixVal(map.PixIndexSph(thetas,phis(i)))=datas(2*i);
|
|---|
| 95 | }
|
|---|
| 96 | }
|
|---|
| 97 | //cout <<"gen "<<datan<<endl;
|
|---|
| 98 | for (int i=0;i< nph;i++){
|
|---|
| 99 | map.PixVal(map.PixIndexSph(theta,phin(i)))=datan(2*i);
|
|---|
| 100 | }
|
|---|
| 101 | }
|
|---|
| 102 | }
|
|---|
| 103 |
|
|---|
| 104 | void syn_phasgen(const int nsmax,const int nlmax,const int nmmax,
|
|---|
| 105 | const vector< complex<long double> >& datain,
|
|---|
| 106 | int nph, OVector& dataout, double phi0,
|
|---|
| 107 | vector< complex<long double> >& bw){
|
|---|
| 108 |
|
|---|
| 109 | /*=======================================================================
|
|---|
| 110 | dataout(j) = sum_m datain(m) * exp(i*m*phi(j))
|
|---|
| 111 | with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1
|
|---|
| 112 |
|
|---|
| 113 | as the set of frequencies {m} is larger than nph,
|
|---|
| 114 | we wrap frequencies within {0..nph-1}
|
|---|
| 115 | ie m = k*nph + m' with m' in {0..nph-1}
|
|---|
| 116 | then
|
|---|
| 117 | noting bw(m') = exp(i*m'*phi0)
|
|---|
| 118 | * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
|
|---|
| 119 | with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
|
|---|
| 120 | dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
|
|---|
| 121 | = Fourier Transform of bw
|
|---|
| 122 | is real
|
|---|
| 123 |
|
|---|
| 124 | NB nph is not necessarily a power of 2
|
|---|
| 125 |
|
|---|
| 126 | =======================================================================*/
|
|---|
| 127 |
|
|---|
| 128 | int ksign = 1;
|
|---|
| 129 | //long double * data = new long double[nph*2];
|
|---|
| 130 | complex<double>* data = new complex<double>[nph];
|
|---|
| 131 | long double * work = new long double[nph*2];
|
|---|
| 132 |
|
|---|
| 133 | int iw;
|
|---|
| 134 | for (iw=0;iw<nph;iw++){bw[iw]=0;}
|
|---|
| 135 |
|
|---|
| 136 | long double shit=phi0*nph;
|
|---|
| 137 | //try complex<long double> since generally this will not be an int
|
|---|
| 138 | //int kshift = (int) pow(-1,shit);
|
|---|
| 139 | //cout<<shit<<" "<<phi0*nph/M_PI<<" "<<kshift<<endl;
|
|---|
| 140 | // int kshift = (int) pow(-1,kphi0);// ! either 1 or -1
|
|---|
| 141 | // all frequencies [-m,m] are wrapped in [0,nph-1]
|
|---|
| 142 | for (int i=1;i<=2*nmmax+1;i++){
|
|---|
| 143 | int m=i-nmmax-1; //in -nmmax, nmmax
|
|---|
| 144 | iw=((m % nph) +nph) % nph; //between 0 and nph = m'
|
|---|
| 145 | int k=(m-iw)/nph; //number of 'turns'
|
|---|
| 146 | // cout <<(complex<long double>)pow(kshift,k)<<" ";
|
|---|
| 147 | // bw[iw]+=datain[i-1]* (complex<long double>)(pow(kshift,k));//complex number
|
|---|
| 148 | bw[iw]+=datain[i-1]* complex<long double>(cos(shit*k),sin(shit*k));//complex number
|
|---|
| 149 | }
|
|---|
| 150 |
|
|---|
| 151 |
|
|---|
| 152 | // applies the shift in position <-> phase factor in Fourier space
|
|---|
| 153 | for (int iw=1;iw<=nph;iw++){
|
|---|
| 154 | int m=ksign*(iw-1);
|
|---|
| 155 | complex<long double> shit(cos(m*phi0),sin(m*phi0));
|
|---|
| 156 | data[iw-1]=bw[iw-1]*shit;
|
|---|
| 157 | // data[2*(iw-1)]=(bw[iw-1]*shit).real();
|
|---|
| 158 | //data[2*(iw-1)+1]=(bw[iw-1]*shit).imag();
|
|---|
| 159 | }
|
|---|
| 160 |
|
|---|
| 161 | //FFTVector inv((double*) data,nph,true);
|
|---|
| 162 | //FFTVector outv=FFTServer::solve(inv,1);
|
|---|
| 163 | //outv.Vreal(dataout);
|
|---|
| 164 |
|
|---|
| 165 | FFTServer fft;
|
|---|
| 166 | fft.fftf(nph, data);
|
|---|
| 167 | dataout.Realloc(2*nph);
|
|---|
| 168 | for (int i=0;i<nph;i++) {
|
|---|
| 169 | dataout(2*i)=data[i].real();dataout(2*i+1)=data[i].imag();
|
|---|
| 170 | }
|
|---|
| 171 | delete[] data;
|
|---|
| 172 | delete[] work;
|
|---|
| 173 |
|
|---|
| 174 | /*
|
|---|
| 175 | int dum0=1;
|
|---|
| 176 | int dum1=1;
|
|---|
| 177 |
|
|---|
| 178 | fft_gpd_(data,nph,dum0,ksign,dum1,work); // complex to complex
|
|---|
| 179 |
|
|---|
| 180 | for (int iw=0;iw<nph;iw++){dataout(iw) = data[2*iw];}*/
|
|---|
| 181 | }
|
|---|
| 182 |
|
|---|
| 183 | void create_almgen(const int nsmax,const int nlmax,const int nmmax,
|
|---|
| 184 | vector< vector< complex<double> > >& alm,
|
|---|
| 185 | vector<float>& cls_tt,int& iseed,const float fwhm,
|
|---|
| 186 | vector<float> lread){
|
|---|
| 187 |
|
|---|
| 188 |
|
|---|
| 189 | /*=======================================================================
|
|---|
| 190 | creates the a_lm from the power spectrum,
|
|---|
| 191 | assuming they are gaussian and complex
|
|---|
| 192 | with a variance given by C(l)
|
|---|
| 193 |
|
|---|
| 194 | the input file should contain : l and C(l) with *consecutive* l's
|
|---|
| 195 | (missing C(l) are put to 0.)
|
|---|
| 196 |
|
|---|
| 197 | because the map is real we have : a_l-m = (-)^m conjug(a_lm)
|
|---|
| 198 | so we actually compute them only for m >= 0
|
|---|
| 199 |
|
|---|
| 200 | modifie G. Le Meur (13/01/98) :
|
|---|
| 201 | on ne lit plus les C(l) sur un fichier. Le tableau est entre en argument
|
|---|
| 202 | (cls_tt). Ce tableau doit contenir les valeur de C(l) par ordre
|
|---|
| 203 | SEQUENTIEL de l (de l=0 a l=nlmax)
|
|---|
| 204 | =======================================================================*/
|
|---|
| 205 | /* CHARACTER*128 filename
|
|---|
| 206 |
|
|---|
| 207 | integer unit,n_l,j,il,im
|
|---|
| 208 | real fs,quadrupole,xxxtmp,correct
|
|---|
| 209 | real fs_tt
|
|---|
| 210 | real zeta1_r, zeta1_i
|
|---|
| 211 | LOGICAL ok
|
|---|
| 212 | character*20 string_quad
|
|---|
| 213 |
|
|---|
| 214 | real*4 prefact(0:2)
|
|---|
| 215 | logical prefact_bool
|
|---|
| 216 |
|
|---|
| 217 | integer lneffct
|
|---|
| 218 | real gasdev2
|
|---|
| 219 | external lneffct, gasdev2
|
|---|
| 220 |
|
|---|
| 221 | c-----------------------------------------------------------------------*/
|
|---|
| 222 | alm.resize(nlmax+1);
|
|---|
| 223 | for (int i=0; i< (signed) alm.size();i++){
|
|---|
| 224 | alm[i].resize(nmmax+1);
|
|---|
| 225 | }
|
|---|
| 226 | //cout << "cing createalm"<<endl;
|
|---|
| 227 | iseed = -abs(iseed);
|
|---|
| 228 | if (iseed == 0){iseed=-1;}
|
|---|
| 229 | int idum = iseed;
|
|---|
| 230 |
|
|---|
| 231 | float sig_smooth = fwhm/sqrt(8.*log(2.))/(60.*180.)* M_PI;
|
|---|
| 232 |
|
|---|
| 233 | for (int i=0;i<nlmax+1;i++){lread[i]=0;}
|
|---|
| 234 |
|
|---|
| 235 | for (int i=0;i <= nlmax;i++){lread[i] = i;}
|
|---|
| 236 |
|
|---|
| 237 | int n_l = nlmax+1;
|
|---|
| 238 |
|
|---|
| 239 | cout<<lread[0]<<" <= l <= "<<lread[n_l-1]<<endl;
|
|---|
| 240 |
|
|---|
| 241 | // --- smoothes the initial power spectrum ---
|
|---|
| 242 |
|
|---|
| 243 | for (int i=0;i<n_l;i++){
|
|---|
| 244 | int l= (int) lread[i];
|
|---|
| 245 | float gauss=exp(-l*(l+1.)*sig_smooth*sig_smooth);
|
|---|
| 246 | cls_tt[i]*=gauss;
|
|---|
| 247 | }
|
|---|
| 248 |
|
|---|
| 249 |
|
|---|
| 250 | // --- generates randomly the alm according to their power spectrum ---
|
|---|
| 251 | float hsqrt2 = 1.0 / sqrt(2.0);
|
|---|
| 252 |
|
|---|
| 253 | for (int i=0;i<n_l;i++){
|
|---|
| 254 | int l=(int) lread[i];
|
|---|
| 255 | float rms_tt=sqrt(cls_tt[i]);
|
|---|
| 256 | // ------ m = 0 ------
|
|---|
| 257 | complex<float> zeta1(gasdev2_(idum));
|
|---|
| 258 | //cout << "c2 "<<idum << " "<<zeta1 <<endl;
|
|---|
| 259 | alm[l][0] = zeta1 * rms_tt;
|
|---|
| 260 |
|
|---|
| 261 | //------ m > 0 ------
|
|---|
| 262 | for (int m=1;m<=l;m++){
|
|---|
| 263 | complex<float> shit1(hsqrt2);
|
|---|
| 264 | // cout << "c "<<idum << endl;
|
|---|
| 265 | complex<float> shit2(gasdev2_(idum),gasdev2_(idum));
|
|---|
| 266 | zeta1=shit1*shit2;
|
|---|
| 267 | alm[l][m]=rms_tt*zeta1;
|
|---|
| 268 | }
|
|---|
| 269 | }
|
|---|
| 270 | }
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
