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