source: Sophya/trunk/Poubelle/DPC:FitsIOServer/Samba/syn2fast.cc@ 3252

Last change on this file since 3252 was 658, checked in by ansari, 26 years ago

no message
File size: 10.4 KB
RevLine 
[658]1#include <math.h>
2#include "syn2fast.h"
3#include "lambuilder.h"
4#include "fftserver.h"
5
6extern "C" {
7 // void fft_gpd_(long double* ,int& ,int& ,int& ,int& ,long double*);
8 float gasdev2_(int& idum);
9 }
10
11void 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
150void 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
214void 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
255c-----------------------------------------------------------------------*/
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}
Note: See TracBrowser for help on using the repository browser.