CosmicFish  Reference documentation for version 1.0
Looking into future Cosmology
007_FUTURCMB_lensnoise.f90
Go to the documentation of this file.
1 !----------------------------------------------------------------------------------------
2 !
3 ! This file is part of CosmicFish.
4 !
5 ! Copyright (C) 2015-2016 by the CosmicFish authors
6 !
7 ! The CosmicFish code is free software;
8 ! You can use it, redistribute it, and/or modify it under the terms
9 ! of the GNU General Public License as published by the Free Software Foundation;
10 ! either version 3 of the License, or (at your option) any later version.
11 ! The full text of the license can be found in the file LICENSE at
12 ! the top level of the CosmicFish distribution.
13 !
14 !----------------------------------------------------------------------------------------
15 
21 
22 !----------------------------------------------------------------------------------------
24 
25 module lensnoise
26  !+
27  !PURPOSE:
28  ! calcul of the deflection field noise power spectrum
29  !AUTHOR:
30  ! J. Lesgourgues
31  !REF:
32  ! Okamoto & Hu's algorithm
33  ! for the calculation of the quadratic estimator variance
34  !VERSION:
35  ! june 2006: made an indep. module by LP
36  ! june 2005: written by JL
37  !-
38 
39  use amlutils
40  use precision
41 
42  implicit none
43 
44  private
45 
46  public :: calc_nldd
47 
48 contains
49 
50  ! ---------------------------------------------------------------------------------------------
52  subroutine calc_nldd (incls, inlcls, clnoise, nldd, jjmaxTmax, lmaxM)
53  ! file_cl_plus(1)='../fidO/fidOa5_scalCls.dat'
54  ! les Cl lus dans CAMB sont a convertir en muK^2
55 
56  ! declaration des arguments :
57  REAL(dl), INTENT(IN), DIMENSION (:,:) :: incls
63 
64  REAL(dl), INTENT(IN), DIMENSION (:,:) :: inlcls
69 
70  REAL(dl), INTENT(IN), DIMENSION (:,:) :: clnoise
73 
74  real(dl), INTENT(OUT), DIMENSION (:) :: nldd
75 
76  integer , INTENT(INOUT) :: jjmaxTmax
78 
79  integer , INTENT(IN) :: lmaxM
80 
81  ! les variables locales
82  integer jjmaxT, jj, ll, l1, l2max, nn, l2, l, err
83  integer, dimension(:), allocatable :: llens
84  double precision all, gall, al1, an, al2
85  double precision aAlTT, aAlEE, aAlTE, aAlTB, aAlEB
86  double precision aBlTTEE, aBlTTTE, aBlEETE, aBlTBEB
87  double precision Al1L0, Al1L2, wigner0, wigner2
88  double precision parity_test, parity, prefactor
89  double precision fTT12, fEE12, fTE12, fTE21
90  double precision fTB12, fTB21, fEB12, fEB21
91  double precision gTE12_den, gTE12_num, gTE21_num
92  double precision Wl20_sup, Wl22_sup, Wl20_mid, Wl22_mid, Wl20_inf, Wl22_inf
93  double precision aa0, bb0, cc0, aa2, bb2, cc2
94  double precision, dimension(5,5) :: LensingNoise, InvLensingNoise
95  double precision GaussVecLens(5,1)
96  double precision, dimension(:), allocatable :: MinVarNoise_L
97  double precision, dimension(:), allocatable :: ddMinVarNoise_L
98  double precision, dimension(:), allocatable :: alogllens
99  integer bonalloc
100  CHARACTER(LEN=100), PARAMETER :: code = 'calc_NLDD'
101 
102 
103  ! debut !
104  ! compute noise for lensing
105 
106  ! define values for interpolation
107  jjmaxt=int((dble(lmaxm)+1700.d0)/50.d0)
108 
109  if (jjmaxt.gt.jjmaxtmax) then
110  jjmaxtmax = jjmaxt+1
111  write(*,*)'WARNING: jjmaxTmax increased to: ',jjmaxt
112  end if
113 
114  allocate(minvarnoise_l(1:jjmaxt), &
115  & alogllens(1:jjmaxt), &
116  & llens(1:jjmaxt), stat = bonalloc)
117  if (bonalloc > 0) then
118  print*,code,'> pb allocation memoire', bonalloc
119  stop
120  endif
121  !allocate(ddMinVarNoise_L(1::jjmaxT))
122 
123  llens(1)=2
124  llens(2)=3
125  llens(3)=4
126  llens(4)=5
127  llens(5)=6
128  llens(6)=7
129  llens(7)=8
130  llens(8)=9
131  do jj=9,17
132  llens(jj)=jj*10-80
133  end do
134  do jj=18,53
135  llens(jj)=jj*25-350
136  end do
137  do jj=54,jjmaxt
138  llens(jj)=jj*50-1700
139  end do
140 
141  ! sum over L
142 
143 ! open(unit=37,file="test_Nls.dat",form="formatted", &
144 ! & status='replace', iostat=err)
145 ! if (err /= 0) print *, "Probleme d'ecriture du fichier test..."
146 
147  do jj=1,jjmaxt !! do 1
148 
149  ll=llens(jj)
150  all=dble(ll)
151 
152  !! test !!
153  !print *, code, incls(1,ll), incls(2,ll), incls(3,ll)
154  !print *, code, inlcls(1,ll),inlcls(2,ll),inlcls(3,ll),inlcls(4,ll)
155  !print *, code, clnoise(1,ll),clnoise(2,ll)
156 
157  ! determine where to cut the sum over l1 in two parts
158 
159  gall=all/2.d0
160  if (dint(gall).ne.gall) gall=gall+0.5d0
161  if (gall.le.2.d0) gall=2.d0
162 
163  ! reset the variable accounting for the various sums over l1,l2
164 
165  aaltt=0.d0
166  aalee=0.d0
167  aalte=0.d0
168  aaltb=0.d0
169  aaleb=0.d0
170  ablttee=0.d0
171  ablttte=0.d0
172  ableete=0.d0
173  abltbeb=0.d0
174 
175  ! start loop on l1 (first part)
176 
177  do l1=2,int(gall)-1 !! do 2
178  al1=dble(l1)
179 
180  ! determine maximum value of l2 in the sum over l2
181 
182  l2max=ll+l1
183  ! if (l2max.gt.lmaxTP) l2max=lmaxTP
184 
185  !compute A(l1,L,0) (useful for first values of W_l2)
186 
187  al1l0=1.d0
188 
189  do nn=1,ll
190  an=dble(nn)
191  al1l0=al1l0*(an-0.5d0)*(al1+an)/an/(al1+an-0.5d0)
192  end do
193 
194  al1l0=dsqrt(al1l0)
195 
196  al1l2=al1l0 &
197  & *dsqrt(al1*(al1-1.d0)/(al1+1.d0)/(al1+2.d0)) &
198  & *dsqrt((al1+all+1.d0)*(al1+all+2.d0) &
199  & /(al1+all)/(al1+all-1.d0))
200 
201  ! case l2=l2max
202 
203  l2=l2max
204  al2=dble(l2)
205 
206  wl20_sup=al1l0/dsqrt(2.d0*(all+al1)+1.d0)
207  wl22_sup=al1l2/dsqrt(2.d0*(all+al1)+1.d0)
208 
209  parity_test=(all+al1)/2.d0
210 
211  if (dint(parity_test).ne.parity_test) then
212  wl20_sup=-wl20_sup
213  wl22_sup=-wl22_sup
214  end if
215 
216  parity=(al1+all+al2)/2.d0
217 
218  if (l2.le.lmaxm) then
219 
220  wigner0=wl20_sup
221  wigner2=wl22_sup
222 
223  prefactor=dsqrt((2.d0*all+1.d0)*(2.d0*al1+1.d0) &
224  & *(2.d0*al2+1.d0)/(16.d0*pi))
225 
226  if (dint(parity).eq.parity) then
227 
228  ftt12=prefactor*wigner0*(incls(1,l1) &
229  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
230  & +incls(1,l2) &
231  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
232  fee12=prefactor*wigner2*(incls(2,l1) &
233  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
234  & +incls(2,l2) &
235  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
236  fte12=prefactor*(wigner2*incls(3,l1) &
237  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
238  & +wigner0*incls(3,l2) &
239  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
240  fte21=prefactor*(wigner2*incls(2,l2) &
241  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)) &
242  & +wigner0*incls(3,l1) &
243  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)))
244 
245  gte12_den=(inlcls(1,l1)+clnoise(1,l1)) &
246  & *(inlcls(1,l2)+clnoise(1,l2)) &
247  & *(inlcls(2,l1)+clnoise(2,l1)) &
248  & *(inlcls(2,l2)+clnoise(2,l2)) &
249  & -(inlcls(3,l1)*inlcls(3,l2))**2
250  gte12_num=(inlcls(1,l2)+clnoise(1,l2)) &
251  & *(inlcls(2,l1)+clnoise(2,l1)) &
252  & *fte12-(inlcls(3,l1)) &
253  & *(inlcls(3,l2))*fte21
254  gte21_num=(inlcls(1,l1)+clnoise(1,l1)) &
255  & *(inlcls(2,l2)+clnoise(2,l2)) &
256  & *fte21-(inlcls(3,l1)) &
257  & *(inlcls(3,l2))*fte12
258 
259  aaltt = aaltt +ftt12**2 &
260  & /(inlcls(1,l1)+clnoise(1,l1))/(inlcls(1,l2)+clnoise(1,l2))
261  aalee = aalee +fee12**2 &
262  & /(inlcls(2,l1)+clnoise(2,l1))/(inlcls(2,l2)+clnoise(2,l2))
263  aalte = aalte + (fte12*gte12_num+fte21*gte21_num) &
264  & /gte12_den
265 
266  ablttee= ablttee+ftt12*fee12 &
267  & *(inlcls(3,l1))*(inlcls(3,l2)) &
268  & /(inlcls(1,l1)+clnoise(1,l1))/(inlcls(1,l2)+clnoise(1,l2)) &
269  & /(inlcls(2,l1)+clnoise(2,l1))/(inlcls(2,l2)+clnoise(2,l2))
270  ablttte= ablttte+ftt12/gte12_den &
271  & /(inlcls(1,l1)+clnoise(1,l1))/(inlcls(1,l2)+clnoise(1,l2)) &
272  & *((inlcls(1,l1)+clnoise(1,l1))*(inlcls(3,l2)) &
273  & *gte12_num &
274  & +(inlcls(1,l2)+clnoise(1,l2))*(inlcls(3,l1)) &
275  & *gte21_num)
276  ableete= ableete+fee12/gte12_den &
277  & /(inlcls(2,l1)+clnoise(2,l1))/(inlcls(2,l2)+clnoise(2,l2)) &
278  & *((inlcls(2,l2)+clnoise(2,l2))*(inlcls(3,l1)) &
279  & *gte12_num &
280  & +(inlcls(2,l1)+clnoise(2,l1))*(inlcls(3,l2)) &
281  & *gte21_num)
282 
283  else
284 
285  ftb12=prefactor*wigner2*incls(3,l1) &
286  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
287  ftb21=prefactor*wigner2*incls(3,l2) &
288  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
289  feb12=prefactor*wigner2*incls(2,l1) &
290  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
291  feb21=prefactor*wigner2*incls(2,l2) &
292  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
293 
294  aaltb = aaltb &
295  & + ftb12**2/(inlcls(1,l1)+clnoise(1,l1)) &
296  & /(inlcls(4,l2)+clnoise(2,l2)) &
297  & + ftb21**2/(inlcls(1,l2)+clnoise(1,l2)) &
298  & /(inlcls(4,l1)+clnoise(2,l1))
299  aaleb = aaleb &
300  & + feb12**2/(inlcls(2,l1)+clnoise(2,l1)) &
301  & /(inlcls(4,l2)+clnoise(2,l2)) &
302  & + feb21**2/(inlcls(2,l2)+clnoise(2,l2)) &
303  & /(inlcls(4,l1)+clnoise(2,l1))
304  abltbeb = abltbeb &
305  & + ftb12*feb12*(inlcls(3,l1)) &
306  & /(inlcls(4,l2)+clnoise(2,l2)) &
307  & /(inlcls(1,l1)+clnoise(1,l1)) &
308  & /(inlcls(2,l1)+clnoise(2,l1)) &
309  & + ftb21*feb21*(inlcls(3,l2)) &
310  & /(inlcls(4,l1)+clnoise(2,l1)) &
311  & /(inlcls(1,l2)+clnoise(1,l2))/(inlcls(2,l2)+clnoise(2,l2))
312 
313  end if
314  end if
315 
316  ! case l2=l2max-1
317 
318  l2=l2max-1
319  al2=dble(l2)
320 
321  wl20_mid=0.
322  wl22_mid=al1l2*2.d0* &
323  & dsqrt(all/al1/(al1+all-2.d0)/(al1+all+2.d0))
324 
325  if (dint(parity_test).ne.parity_test) then
326  wl22_mid=-wl22_mid
327  end if
328 
329  parity=(al1+all+al2)/2.d0
330 
331  if (l2.le.lmaxm) then
332 
333  wigner0=wl20_mid
334  wigner2=wl22_mid
335 
336  prefactor=dsqrt((2.d0*all+1.d0)*(2.d0*al1+1.d0) &
337  & *(2.d0*al2+1.d0)/(16.d0*pi))
338 
339  if (dint(parity).eq.parity) then
340  ftt12=prefactor*wigner0*(incls(1,l1) &
341  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
342  & +incls(1,l2) &
343  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
344  fee12=prefactor*wigner2*(incls(2,l1) &
345  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
346  & +incls(2,l2) &
347  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
348  fte12=prefactor*(wigner2*incls(3,l1) &
349  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
350  & +wigner0*incls(3,l2) &
351  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
352  fte21=prefactor*(wigner2*incls(3,l2) &
353  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)) &
354  & +wigner0*incls(3,l1) &
355  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)))
356 
357  gte12_den=(inlcls(1,l1)+clnoise(1,l1)) &
358  & *(inlcls(1,l2)+clnoise(1,l2)) &
359  & *(inlcls(2,l1)+clnoise(2,l1)) &
360  & *(inlcls(2,l2)+clnoise(2,l2)) &
361  & -((inlcls(3,l1)) &
362  & *(inlcls(3,l2)))**2
363  gte12_num=(inlcls(1,l2)+clnoise(1,l2)) &
364  & *(inlcls(2,l1)+clnoise(2,l1))*fte12-(inlcls(3,l1)) &
365  & *(inlcls(3,l2))*fte21
366  gte21_num=(inlcls(1,l1)+clnoise(1,l1)) &
367  & *(inlcls(2,l2)+clnoise(2,l2))*fte21-(inlcls(3,l1)) &
368  & *(inlcls(3,l2))*fte12
369 
370  aaltt = aaltt +ftt12**2 &
371  & /(inlcls(1,l1)+clnoise(1,l1))/(inlcls(1,l2)+clnoise(1,l2))
372  aalee = aalee +fee12**2 &
373  & /(inlcls(2,l1)+clnoise(2,l1))/(inlcls(2,l2)+clnoise(2,l2))
374  aalte = aalte + (fte12*gte12_num+fte21*gte21_num) &
375  & /gte12_den
376 
377  ablttee= ablttee+ftt12*fee12 &
378  & *(inlcls(3,l1))*(inlcls(3,l2)) &
379  & /(inlcls(1,l1)+clnoise(1,l1))/(inlcls(1,l2)+clnoise(1,l2)) &
380  & /(inlcls(2,l1)+clnoise(2,l1))/(inlcls(2,l2)+clnoise(2,l2))
381  ablttte= ablttte+ftt12/gte12_den &
382  & /(inlcls(1,l1)+clnoise(1,l1))/(inlcls(1,l2)+clnoise(1,l2)) &
383  & *((inlcls(1,l1)+clnoise(1,l1))*(inlcls(3,l2)) &
384  & *gte12_num &
385  & +(inlcls(1,l2)+clnoise(1,l2))*(inlcls(3,l1)) &
386  & *gte21_num)
387  ableete= ableete+fee12/gte12_den &
388  & /(inlcls(2,l1)+clnoise(2,l1))/(inlcls(2,l2)+clnoise(2,l2)) &
389  & *((inlcls(2,l2)+clnoise(2,l2))*(inlcls(3,l1)) &
390  & *gte12_num &
391  & +(inlcls(2,l1)+clnoise(2,l1))*(inlcls(3,l2)) &
392  & *gte21_num)
393 
394  else
395 
396  ftb12=prefactor*wigner2*incls(3,l1) &
397  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
398  ftb21=prefactor*wigner2*incls(3,l2) &
399  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
400  feb12=prefactor*wigner2*inlcls(2,l1) &
401  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
402  feb21=prefactor*wigner2*inlcls(2,l2) &
403  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
404 
405  aaltb = aaltb &
406  & + ftb12**2/(inlcls(1,l1)+clnoise(1,l1)) &
407  & /(inlcls(4,l2)+clnoise(2,l2)) &
408  & + ftb21**2/(inlcls(1,l2)+clnoise(1,l2)) &
409  & /(inlcls(4,l1)+clnoise(2,l1))
410  aaleb = aaleb &
411  & + feb12**2/(inlcls(2,l1)+clnoise(2,l1)) &
412  & /(inlcls(4,l2)+clnoise(2,l2)) &
413  & + feb21**2/(inlcls(2,l2)+clnoise(2,l2)) &
414  & /(inlcls(4,l1)+clnoise(2,l1))
415 
416  abltbeb = abltbeb &
417  & + ftb12*feb12*(inlcls(3,l1)) &
418  & /(inlcls(4,l2)+clnoise(2,l2)) &
419  & /(inlcls(1,l1)+clnoise(1,l1)) &
420  & /(inlcls(2,l1)+clnoise(2,l1)) &
421  & + ftb21*feb21*(inlcls(3,l2)) &
422  & /(inlcls(4,l1)+clnoise(2,l1)) &
423  & /(inlcls(1,l2)+clnoise(1,l2)) &
424  & /(inlcls(2,l2)+clnoise(2,l2))
425 
426  end if
427  end if
428 
429  ! sum over other cases, down to l2=L-l1
430 
431  do nn=1,l2max-ll+l1-1 !! do 3
432 
433  l2=l2max-1-nn
434  al2=dble(l2)
435 
436  ! recursion formula giving W_l2
437 
438  aa0=(al2+2.d0)*(al2+1.d0)* &
439  & dsqrt((-al2+all+al1)*(al2-all+al1+1.d0) &
440  & *(al2+all-al1+1.d0)*(al2+all+al1+2.d0))
441  bb0=0.d0
442  cc0=-(al2+1.d0)*(al2+2.d0)* &
443  & dsqrt((-al2+all+al1-1.d0)*(al2-all+al1+2.d0) &
444  & *(al2+all-al1+2.d0)*(al2+all+al1+3.d0))
445 
446  wl20_inf=(bb0*wl20_mid+cc0*wl20_sup)/aa0
447 
448  aa2=(al2+2.d0)*dsqrt(((al2+1.d0)**2-4.d0) &
449  & *(-al2+all+al1)*(al2-all+al1+1.d0) &
450  & *(al2+all-al1+1.d0)*(al2+all+al1+2.d0))
451  bb2=2.d0*(2.d0*al2+3.d0)*((al2+1.d0)*(al2+2.d0) &
452  & +all*(all+1.d0)-al1*(al1+1.d0))
453  cc2=-(al2+1.d0)*dsqrt(((al2+2.d0)**2-4.d0) &
454  & *(-al2+all+al1-1.d0)*(al2-all+al1+2.d0) &
455  & *(al2+all-al1+2.d0)*(al2+all+al1+3.d0))
456 
457  wl22_inf=(bb2*wl22_mid+cc2*wl22_sup)/aa2
458 
459  parity=(al1+all+al2)/2.d0
460 
461  if (l2.le.lmaxm) then
462 
463  wigner0=wl20_inf
464  wigner2=wl22_inf
465 
466  prefactor=dsqrt((2.d0*all+1.d0)*(2.d0*al1+1.d0) &
467  & *(2.d0*al2+1.d0)/(16.d0*3.14159d0))
468 
469  if (dint(parity).eq.parity) then
470 
471  ftt12=prefactor*wigner0*(incls(1,l1) &
472  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
473  & +incls(1,l2) &
474  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
475  fee12=prefactor*wigner2*(incls(2,l1) &
476  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
477  & +incls(2,l2) &
478  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
479  fte12=prefactor*(wigner2*incls(3,l1) &
480  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
481  & +wigner0*incls(3,l2) &
482  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
483  fte21=prefactor*(wigner2*incls(3,l2) &
484  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)) &
485  & +wigner0*incls(3,l1) &
486  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)))
487 
488  gte12_den=(inlcls(1,l1)+clnoise(1,l1)) &
489  & *(inlcls(1,l2)+clnoise(1,l2)) &
490  & *(inlcls(2,l1)+clnoise(2,l1)) &
491  & *(inlcls(2,l2)+clnoise(2,l2)) &
492  & -((inlcls(3,l1)) &
493  & *(inlcls(3,l2)))**2
494  gte12_num=(inlcls(1,l2)+clnoise(1,l2)) &
495  & *(inlcls(2,l1)+clnoise(2,l1)) &
496  & *fte12-(inlcls(3,l1)) &
497  & *(inlcls(3,l2))*fte21
498  gte21_num=(inlcls(1,l1)+clnoise(1,l1)) &
499  & *(inlcls(2,l2)+clnoise(2,l2)) &
500  & *fte21-(inlcls(3,l1)) &
501  & *(inlcls(3,l2))*fte12
502 
503  aaltt = aaltt +ftt12**2 &
504  & /(inlcls(1,l1)+clnoise(1,l1)) &
505  & /(inlcls(1,l2)+clnoise(1,l2))
506  aalee = aalee +fee12**2 &
507  & /(inlcls(2,l1)+clnoise(2,l1)) &
508  & /(inlcls(2,l2)+clnoise(2,l2))
509  aalte = aalte + (fte12*gte12_num+fte21*gte21_num) &
510  & /gte12_den
511 
512  ablttee= ablttee+ftt12*fee12 &
513  & *(inlcls(3,l1))*(inlcls(3,l2)) &
514  & /(inlcls(1,l1)+clnoise(1,l1)) &
515  & /(inlcls(1,l2)+clnoise(1,l2)) &
516  & /(inlcls(2,l1)+clnoise(2,l1)) &
517  & /(inlcls(2,l2)+clnoise(2,l2))
518  ablttte= ablttte+ftt12/gte12_den &
519  & /(inlcls(1,l1)+clnoise(1,l1)) &
520  & /(inlcls(1,l2)+clnoise(1,l2)) &
521  & *((inlcls(1,l1)+clnoise(1,l1))*(inlcls(3,l2)) &
522  & *gte12_num &
523  & +(inlcls(1,l2)+clnoise(1,l2))*(inlcls(3,l1)) &
524  & *gte21_num)
525  ableete= ableete+fee12/gte12_den &
526  & /(inlcls(2,l1)+clnoise(2,l1)) &
527  & /(inlcls(2,l2)+clnoise(2,l2)) &
528  & *((inlcls(2,l2)+clnoise(2,l2))*(inlcls(3,l1)) &
529  & *gte12_num &
530  & +(inlcls(2,l1)+clnoise(2,l1))*(inlcls(3,l2)) &
531  & *gte21_num)
532 
533  else
534 
535  ftb12=prefactor*wigner2*incls(3,l1) &
536  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
537  ftb21=prefactor*wigner2*incls(3,l2) &
538  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
539  feb12=prefactor*wigner2*incls(2,l1) &
540  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
541  feb21=prefactor*wigner2*incls(2,l2) &
542  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
543 
544  aaltb = aaltb &
545  & + ftb12**2/(inlcls(1,l1)+clnoise(1,l1)) &
546  & /(inlcls(4,l2)+clnoise(2,l2)) &
547  & + ftb21**2/(inlcls(1,l2)+clnoise(1,l2)) &
548  & /(inlcls(4,l1)+clnoise(2,l1))
549  aaleb = aaleb &
550  & + feb12**2/(inlcls(2,l1)+clnoise(2,l1)) &
551  & /(inlcls(4,l2)+clnoise(2,l2)) &
552  & + feb21**2/(inlcls(2,l2)+clnoise(2,l2)) &
553  & /(inlcls(4,l1)+clnoise(2,l1))
554 
555  abltbeb = abltbeb &
556  & + ftb12*feb12*(inlcls(3,l1)) &
557  & /(inlcls(4,l2)+clnoise(2,l2)) &
558  & /(inlcls(1,l1)+clnoise(1,l1)) &
559  & /(inlcls(2,l1)+clnoise(2,l1)) &
560  & + ftb21*feb21*(inlcls(3,l2)) &
561  & /(inlcls(4,l1)+clnoise(2,l1)) &
562  & /(inlcls(1,l2)+clnoise(1,l2)) &
563  & /(inlcls(2,l2)+clnoise(2,l2))
564 
565  end if
566  end if
567 
568  ! shift for the next recursion
569 
570  wl20_sup=wl20_mid
571  wl20_mid=wl20_inf
572 
573  wl22_sup=wl22_mid
574  wl22_mid=wl22_inf
575 
576  end do !! do 3
577 
578  end do !! do 2
579 
580  ! start loop on l1 (second part)
581 
582  do l1=int(gall),lmaxm !! do 2
583  al1=dble(l1)
584 
585  ! determine maximum value of l2 in the sum over l2
586 
587  l2max=ll+l1
588  ! if (l2max.gt.lmaxTP) l2max=lmaxTP
589 
590  ! compute A(l1,L,0) (useful for first values of W_l2)
591 
592  al1l0=1.d0
593 
594  do nn=1,ll
595  an=dble(nn)
596  al1l0=al1l0*(an-0.5d0)*(al1+an)/an/(al1+an-0.5d0)
597  end do
598 
599  al1l0=dsqrt(al1l0)
600 
601  al1l2=al1l0 &
602  & *dsqrt(al1*(al1-1.d0)/(al1+1.d0)/(al1+2.d0)) &
603  & *dsqrt((al1+all+1.d0)*(al1+all+2.d0) &
604  & /(al1+all)/(al1+all-1.d0))
605 
606  ! case l2=l2max
607 
608  l2=l2max
609  al2=dble(l2)
610 
611  wl20_sup=al1l0/dsqrt(2.d0*(all+al1)+1.d0)
612  wl22_sup=al1l2/dsqrt(2.d0*(all+al1)+1.d0)
613 
614  parity_test=(all+al1)/2.d0
615 
616  if (dint(parity_test).ne.parity_test) then
617  wl20_sup=-wl20_sup
618  wl22_sup=-wl22_sup
619  end if
620 
621  parity=(al1+all+al2)/2.d0
622 
623  if (l2.le.lmaxm) then
624 
625  wigner0=wl20_sup
626  wigner2=wl22_sup
627 
628  if (dint(parity).eq.parity) then
629 
630  prefactor=dsqrt((2.d0*all+1.d0)*(2.d0*al1+1.d0) &
631  & *(2.d0*al2+1.d0)/(16.d0*pi))
632 
633  ftt12=prefactor*wigner0*(incls(1,l1) &
634  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
635  & +incls(1,l2) &
636  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
637  fee12=prefactor*wigner2*(incls(2,l1) &
638  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
639  & +incls(2,l2) &
640  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
641  fte12=prefactor*(wigner2*incls(3,l1) &
642  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
643  & +wigner0*incls(3,l2) &
644  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
645  fte21=prefactor*(wigner2*incls(3,l2) &
646  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)) &
647  & +wigner0*incls(3,l1) &
648  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)))
649 
650  gte12_den=(inlcls(1,l1)+clnoise(1,l1)) &
651  & *(inlcls(1,l2)+clnoise(1,l2)) &
652  & *(inlcls(2,l1)+clnoise(2,l1)) &
653  & *(inlcls(2,l2)+clnoise(2,l2)) &
654  & -((inlcls(3,l1)) &
655  & *(inlcls(3,l2)))**2
656  gte12_num=(inlcls(1,l2)+clnoise(1,l2)) &
657  & *(inlcls(2,l1)+clnoise(2,l1))*fte12-(inlcls(3,l1)) &
658  & *(inlcls(3,l2))*fte21
659  gte21_num=(inlcls(1,l1)+clnoise(1,l1)) &
660  & *(inlcls(2,l2)+clnoise(2,l2))*fte21-(inlcls(3,l1)) &
661  & *(inlcls(3,l2))*fte12
662 
663  aaltt = aaltt +ftt12**2 &
664  & /(inlcls(1,l1)+clnoise(1,l1)) &
665  & /(inlcls(1,l2)+clnoise(1,l2))
666  aalee = aalee +fee12**2 &
667  & /(inlcls(2,l1)+clnoise(2,l1)) &
668  & /(inlcls(2,l2)+clnoise(2,l2))
669  aalte = aalte + (fte12*gte12_num+fte21*gte21_num) &
670  & /gte12_den
671 
672  ablttee= ablttee+ftt12*fee12 &
673  & *(inlcls(3,l1))*(inlcls(3,l2)) &
674  & /(inlcls(1,l1)+clnoise(1,l1)) &
675  & /(inlcls(1,l2)+clnoise(1,l2)) &
676  & /(inlcls(2,l1)+clnoise(2,l1)) &
677  & /(inlcls(2,l2)+clnoise(2,l2))
678  ablttte= ablttte+ftt12/gte12_den &
679  & /(inlcls(1,l1)+clnoise(1,l1)) &
680  & /(inlcls(1,l2)+clnoise(1,l2)) &
681  & *((inlcls(1,l1)+clnoise(1,l1))*(inlcls(3,l2)) &
682  & *gte12_num &
683  & +(inlcls(1,l2)+clnoise(1,l2))*(inlcls(3,l1)) &
684  & *gte21_num)
685  ableete= ableete+fee12/gte12_den &
686  & /(inlcls(2,l1)+clnoise(2,l1)) &
687  & /(inlcls(2,l2)+clnoise(2,l2)) &
688  & *((inlcls(2,l2)+clnoise(2,l2))*(inlcls(3,l1)) &
689  & *gte12_num &
690  & +(inlcls(2,l1)+clnoise(2,l1))*(inlcls(3,l2)) &
691  & *gte21_num)
692 
693  else
694 
695  ftb12=prefactor*wigner2*incls(3,l1) &
696  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
697  ftb21=prefactor*wigner2*incls(3,l2) &
698  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
699  feb12=prefactor*wigner2*incls(2,l1) &
700  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
701  feb21=prefactor*wigner2*incls(2,l2) &
702  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
703 
704  aaltb = aaltb &
705  & + ftb12**2/(inlcls(1,l1)+clnoise(1,l1)) &
706  & /(inlcls(4,l2)+clnoise(2,l2)) &
707  & + ftb21**2/(inlcls(1,l2)+clnoise(1,l2)) &
708  & /(inlcls(4,l1)+clnoise(2,l1))
709  aaleb = aaleb &
710  & + feb12**2/(inlcls(2,l1)+clnoise(2,l1)) &
711  & /(inlcls(4,l2)+clnoise(2,l2)) &
712  & + feb21**2/(inlcls(2,l2)+clnoise(2,l2)) &
713  & /(inlcls(4,l1)+clnoise(2,l1))
714 
715  abltbeb = abltbeb &
716  & + ftb12*feb12*(inlcls(3,l1)) &
717  & /(inlcls(4,l2)+clnoise(2,l2)) &
718  & /(inlcls(1,l1)+clnoise(1,l1)) &
719  & /(inlcls(2,l1)+clnoise(2,l1)) &
720  & + ftb21*feb21*(inlcls(3,l2)) &
721  & /(inlcls(4,l1)+clnoise(2,l1)) &
722  & /(inlcls(1,l2)+clnoise(1,l2)) &
723  & /(inlcls(2,l2)+clnoise(2,l2))
724 
725  end if
726  end if
727 
728  ! case l2=l2max-1
729 
730  l2=l2max-1
731  al2=dble(l2)
732 
733  wl20_mid=0.
734  wl22_mid=al1l2*2.d0* &
735  & dsqrt(all/al1/(al1+all-2.d0)/(al1+all+2.d0))
736 
737  if (dint(parity_test).ne.parity_test) then
738  wl22_mid=-wl22_mid
739  end if
740 
741  parity=(al1+all+al2)/2.d0
742 
743  if (l2.le.lmaxm) then
744 
745  wigner0=wl20_mid
746  wigner2=wl22_mid
747 
748  prefactor=dsqrt((2.d0*all+1.d0)*(2.d0*al1+1.d0) &
749  & *(2.d0*al2+1.d0)/(16.d0*3.14159d0))
750 
751  if (dint(parity).eq.parity) then
752 
753  ftt12=prefactor*wigner0*(incls(1,l1) &
754  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
755  & +incls(1,l2) &
756  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
757  fee12=prefactor*wigner2*(incls(2,l1) &
758  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
759  & +incls(2,l2) &
760  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
761  fte12=prefactor*(wigner2*incls(3,l1) &
762  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
763  & +wigner0*incls(3,l2) &
764  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
765  fte21=prefactor*(wigner2*incls(3,l2) &
766  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)) &
767  & +wigner0*incls(3,l1) &
768  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)))
769 
770  gte12_den=(inlcls(1,l1)+clnoise(1,l1)) &
771  & *(inlcls(1,l2)+clnoise(1,l2)) &
772  & *(inlcls(2,l1)+clnoise(2,l1)) &
773  & *(inlcls(2,l2)+clnoise(2,l2)) &
774  & -((inlcls(3,l1)) &
775  & *(inlcls(3,l2)))**2
776  gte12_num=(inlcls(1,l2)+clnoise(1,l2)) &
777  & *(inlcls(2,l1)+clnoise(2,l1)) &
778  & *fte12-(inlcls(3,l1)) &
779  & *(inlcls(3,l2))*fte21
780  gte21_num=(inlcls(1,l1)+clnoise(1,l1)) &
781  & *(inlcls(2,l2)+clnoise(2,l2)) &
782  & *fte21-(inlcls(3,l1)) &
783  & *(inlcls(3,l2))*fte12
784 
785  aaltt = aaltt +ftt12**2 &
786  & /(inlcls(1,l1)+clnoise(1,l1)) &
787  & /(inlcls(1,l2)+clnoise(1,l2))
788  aalee = aalee +fee12**2 &
789  & /(inlcls(2,l1)+clnoise(2,l1)) &
790  & /(inlcls(2,l2)+clnoise(2,l2))
791  aalte = aalte + (fte12*gte12_num+fte21*gte21_num) &
792  & /gte12_den
793 
794  ablttee= ablttee+ftt12*fee12 &
795  & *(inlcls(3,l1))*(inlcls(3,l2)) &
796  & /(inlcls(1,l1)+clnoise(1,l1)) &
797  & /(inlcls(1,l2)+clnoise(1,l2)) &
798  & /(inlcls(2,l1)+clnoise(2,l1)) &
799  & /(inlcls(2,l2)+clnoise(2,l2))
800  ablttte= ablttte+ftt12/gte12_den &
801  & /(inlcls(1,l1)+clnoise(1,l1)) &
802  & /(inlcls(1,l2)+clnoise(1,l2)) &
803  & *((inlcls(1,l1)+clnoise(1,l1)) &
804  & *(inlcls(3,l2))*gte12_num &
805  & +(inlcls(1,l2)+clnoise(1,l2)) &
806  & *(inlcls(3,l1))*gte21_num)
807  ableete= ableete+fee12/gte12_den &
808  & /(inlcls(2,l1)+clnoise(2,l1)) &
809  & /(inlcls(2,l2)+clnoise(2,l2)) &
810  & *((inlcls(2,l2)+clnoise(2,l2)) &
811  & *(inlcls(3,l1))*gte12_num &
812  & +(inlcls(2,l1)+clnoise(2,l1)) &
813  & *(inlcls(3,l2))*gte21_num)
814 
815  else
816 
817  ftb12=prefactor*wigner2*incls(3,l1) &
818  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
819  ftb21=prefactor*wigner2*incls(3,l2) &
820  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
821  feb12=prefactor*wigner2*incls(2,l1) &
822  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
823  feb21=prefactor*wigner2*incls(2,l2) &
824  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
825 
826  aaltb = aaltb &
827  & + ftb12**2/(inlcls(1,l1)+clnoise(1,l1)) &
828  & /(inlcls(4,l2)+clnoise(2,l2)) &
829  & + ftb21**2/(inlcls(1,l2)+clnoise(1,l2)) &
830  & /(inlcls(4,l1)+clnoise(2,l1))
831  aaleb = aaleb &
832  & + feb12**2/(inlcls(2,l1)+clnoise(2,l1)) &
833  & /(inlcls(4,l2)+clnoise(2,l2)) &
834  & + feb21**2/(inlcls(2,l2)+clnoise(2,l2)) &
835  & /(inlcls(4,l1)+clnoise(2,l1))
836 
837  abltbeb = abltbeb &
838  & + ftb12*feb12*(inlcls(3,l1)) &
839  & /(inlcls(4,l2)+clnoise(2,l2)) &
840  & /(inlcls(1,l1)+clnoise(1,l1)) &
841  & /(inlcls(2,l1)+clnoise(2,l1)) &
842  & + ftb21*feb21*(inlcls(3,l2)) &
843  & /(inlcls(4,l1)+clnoise(2,l1)) &
844  & /(inlcls(1,l2)+clnoise(1,l2)) &
845  & /(inlcls(2,l2)+clnoise(2,l2))
846 
847  end if
848  end if
849 
850  ! sum over other cases, down to l2=l1+1
851 
852  do nn=1,l2max-1-l1-1 !! do 3
853 
854  l2=l2max-1-nn
855  al2=dble(l2)
856 
857  ! recursion formula giving W_l2
858 
859  aa0=(al2+2.d0)*(al2+1.d0) &
860  & *dsqrt((-al2+all+al1)*(al2-all+al1+1.d0) &
861  & *(al2+all-al1+1.d0)*(al2+all+al1+2.d0))
862  bb0=0.d0
863  cc0=-(al2+1.d0)*(al2+2.d0) &
864  & *dsqrt((-al2+all+al1-1.d0)*(al2-all+al1+2.d0) &
865  & *(al2+all-al1+2.d0)*(al2+all+al1+3.d0))
866 
867  wl20_inf=(bb0*wl20_mid+cc0*wl20_sup)/aa0
868 
869  aa2=(al2+2.d0)*dsqrt(((al2+1.d0)**2-4.d0) &
870  & *(-al2+all+al1)*(al2-all+al1+1.d0) &
871  & *(al2+all-al1+1.d0)*(al2+all+al1+2.d0))
872  bb2=2.d0*(2.d0*al2+3.d0)*((al2+1.d0)*(al2+2.d0) &
873  & +all*(all+1.d0)-al1*(al1+1.d0))
874  cc2=-(al2+1.d0)*dsqrt(((al2+2.d0)**2-4.d0) &
875  & *(-al2+all+al1-1.d0)*(al2-all+al1+2.d0) &
876  & *(al2+all-al1+2.d0)*(al2+all+al1+3.d0))
877 
878  wl22_inf=(bb2*wl22_mid+cc2*wl22_sup)/aa2
879 
880  parity=(al1+all+al2)/2.d0
881 
882  if (l2.le.lmaxm) then
883 
884  wigner0=wl20_inf
885  wigner2=wl22_inf
886 
887  prefactor=dsqrt((2.d0*all+1.d0)*(2.d0*al1+1.d0) &
888  & *(2.d0*al2+1.d0)/(16.d0*3.14159d0))
889 
890  if (dint(parity).eq.parity) then
891 
892  ftt12=prefactor*wigner0*(incls(1,l1) &
893  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
894  & +incls(1,l2) &
895  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
896  fee12=prefactor*wigner2*(incls(2,l1) &
897  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
898  & +incls(2,l2) &
899  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
900  fte12=prefactor*(wigner2*incls(3,l1) &
901  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)) &
902  & +wigner0*incls(3,l2) &
903  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)))
904  fte21=prefactor*(wigner2*incls(3,l2) &
905  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0)) &
906  & +wigner0*incls(3,l1) &
907  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0)))
908 
909  gte12_den=(inlcls(1,l1)+clnoise(1,l1)) &
910  & *(inlcls(1,l2)+clnoise(1,l2)) &
911  & *(inlcls(2,l1)+clnoise(2,l1)) &
912  & *(inlcls(2,l2)+clnoise(2,l2)) &
913  & -((inlcls(3,l1)) &
914  & *(inlcls(3,l2)))**2
915  gte12_num=(inlcls(1,l2)+clnoise(1,l2)) &
916  & *(inlcls(2,l1)+clnoise(2,l1))*fte12 &
917  & -(inlcls(3,l1))*(inlcls(3,l2))*fte21
918  gte21_num=(inlcls(1,l1)+clnoise(1,l1)) &
919  & *(inlcls(2,l2)+clnoise(2,l2))*fte21 &
920  & -(inlcls(3,l1))*(inlcls(3,l2))*fte12
921 
922  aaltt = aaltt +ftt12**2 &
923  & /(inlcls(1,l1)+clnoise(1,l1)) &
924  & /(inlcls(1,l2)+clnoise(1,l2))
925  aalee = aalee +fee12**2 &
926  & /(inlcls(2,l1)+clnoise(2,l1)) &
927  & /(inlcls(2,l2)+clnoise(2,l2))
928  aalte = aalte + (fte12*gte12_num+fte21*gte21_num) &
929  & /gte12_den
930 
931  ablttee= ablttee+ftt12*fee12 &
932  & *(inlcls(3,l1))*(inlcls(3,l2)) &
933  & /(inlcls(1,l1)+clnoise(1,l1)) &
934  & /(inlcls(1,l2)+clnoise(1,l2)) &
935  & /(inlcls(2,l1)+clnoise(2,l1)) &
936  & /(inlcls(2,l2)+clnoise(2,l2))
937  ablttte= ablttte+ftt12/gte12_den &
938  & /(inlcls(1,l1)+clnoise(1,l1)) &
939  & /(inlcls(1,l2)+clnoise(1,l2)) &
940  & *((inlcls(1,l1)+clnoise(1,l1)) &
941  & *(inlcls(3,l2))*gte12_num &
942  & +(inlcls(1,l2)+clnoise(1,l2)) &
943  & *(inlcls(3,l1))*gte21_num)
944  ableete= ableete+fee12/gte12_den &
945  & /(inlcls(2,l1)+clnoise(2,l1)) &
946  & /(inlcls(2,l2)+clnoise(2,l2)) &
947  & *((inlcls(2,l2)+clnoise(2,l2)) &
948  & *(inlcls(3,l1))*gte12_num &
949  & +(inlcls(2,l1)+clnoise(2,l1)) &
950  & *(inlcls(3,l2))*gte21_num)
951 
952  else
953 
954  ftb12=prefactor*wigner2*incls(3,l1) &
955  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
956  ftb21=prefactor*wigner2*incls(3,l2) &
957  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
958  feb12=prefactor*wigner2*incls(2,l1) &
959  & *(all*(all+1.d0)+al1*(al1+1.d0)-al2*(al2+1.d0))
960  feb21=prefactor*wigner2*incls(2,l2) &
961  & *(all*(all+1.d0)+al2*(al2+1.d0)-al1*(al1+1.d0))
962 
963  aaltb = aaltb &
964  & + ftb12**2/(inlcls(1,l1)+clnoise(1,l1)) &
965  & /(inlcls(4,l2)+clnoise(2,l2)) &
966  & + ftb21**2/(inlcls(1,l2)+clnoise(1,l2)) &
967  & /(inlcls(4,l1)+clnoise(2,l1))
968  aaleb = aaleb &
969  & + feb12**2/(inlcls(2,l1)+clnoise(2,l1)) &
970  & /(inlcls(4,l2)+clnoise(2,l2)) &
971  & + feb21**2/(inlcls(2,l2)+clnoise(2,l2)) &
972  & /(inlcls(4,l1)+clnoise(2,l1))
973 
974  abltbeb = abltbeb &
975  & + ftb12*feb12*(inlcls(3,l1)) &
976  & /(inlcls(4,l2)+clnoise(2,l2)) &
977  & /(inlcls(1,l1)+clnoise(1,l1)) &
978  & /(inlcls(2,l1)+clnoise(2,l1)) &
979  & + ftb21*feb21*inlcls(3,l2) &
980  & /(inlcls(4,l1)+clnoise(2,l1)) &
981  & /(inlcls(1,l2)+clnoise(1,l2)) &
982  & /(inlcls(2,l2)+clnoise(2,l2))
983 
984  end if
985  end if
986 
987  ! shift for the next recursion
988 
989  wl20_sup=wl20_mid
990  wl20_mid=wl20_inf
991 
992  wl22_sup=wl22_mid
993  wl22_mid=wl22_inf
994 
995  end do !! do 3
996 
997  ! last case, for l2=l1
998 
999  l2=l1
1000  al2=al1
1001 
1002  ! recursion formula giving W_l2
1003 
1004  aa0=(al2+2.d0)*(al2+1.d0)* &
1005  & dsqrt((-al2+all+al1)*(al2-all+al1+1.d0) &
1006  & *(al2+all-al1+1.d0)*(al2+all+al1+2.d0))
1007  bb0=0.d0
1008  cc0=-(al2+1.d0)*(al2+2.d0)* &
1009  & dsqrt((-al2+all+al1-1.d0)*(al2-all+al1+2.d0) &
1010  & *(al2+all-al1+2.d0)*(al2+all+al1+3.d0))
1011 
1012  wl20_inf=(bb0*wl20_mid+cc0*wl20_sup)/aa0
1013 
1014  aa2=(al2+2.d0)*dsqrt(((al2+1.d0)**2-4.d0) &
1015  & *(-al2+all+al1)*(al2-all+al1+1.d0) &
1016  & *(al2+all-al1+1.d0)*(al2+all+al1+2.d0))
1017  bb2=2.d0*(2.d0*al2+3.d0)*((al2+1.d0)*(al2+2.d0) &
1018  & +all*(all+1.d0)-al1*(al1+1.d0))
1019  cc2=-(al2+1.d0)*dsqrt(((al2+2.d0)**2-4.d0) &
1020  & *(-al2+all+al1-1.d0)*(al2-all+al1+2.d0) &
1021  & *(al2+all-al1+2.d0)*(al2+all+al1+3.d0))
1022 
1023  wl22_inf=(bb2*wl22_mid+cc2*wl22_sup)/aa2
1024 
1025  parity=(al1+all+al2)/2.d0
1026 
1027  if (l2.le.lmaxm) then
1028 
1029  wigner0=wl20_inf
1030  wigner2=wl22_inf
1031 
1032  prefactor=dsqrt((2.d0*all+1.d0)*(2.d0*al1+1.d0) &
1033  & *(2.d0*al2+1.d0)/(16.d0*pi))
1034 
1035  if (dint(parity).eq.parity) then
1036 
1037  ftt12=prefactor*wigner0*2.d0*incls(1,l1)*all*(all+1.d0)
1038  fee12=prefactor*wigner2*2.d0*incls(2,l1)*all*(all+1.d0)
1039  fte12=prefactor*(wigner2+wigner0)*incls(3,l1) &
1040  & *all*(all+1.d0)
1041 
1042  gte12_den=(inlcls(1,l1)+clnoise(1,l1))**2 &
1043  & *(inlcls(2,l1)+clnoise(2,l1))**2 &
1044  & -(inlcls(3,l1))**4
1045  gte12_num=((inlcls(1,l1)+clnoise(1,l1)) &
1046  & *(inlcls(2,l1)+clnoise(2,l1)) &
1047  & -(inlcls(3,l1))**2)*fte21
1048 
1049  aaltt = aaltt + ftt12**2 &
1050  & /(2.d0*(inlcls(1,l1)+clnoise(1,l1))**2)
1051  aalee = aalee + fee12**2 &
1052  & /(2.d0*(inlcls(2,l1)+clnoise(2,l1))**2)
1053  aalte = aalte + fte12*gte12_num &
1054  & /gte12_den
1055 
1056  ablttee= ablttee+ftt12*fee12 &
1057  & /2.d0*(inlcls(3,l1))**2 &
1058  & /(inlcls(1,l1)+clnoise(1,l1))**2 &
1059  & /(inlcls(2,l1)+clnoise(2,l1))**2
1060  ablttte= ablttte+ftt12/gte12_den &
1061  & /(inlcls(1,l1)+clnoise(1,l1)) &
1062  & *(inlcls(3,l1))*gte12_num
1063  ableete= ableete+fee12/gte12_den &
1064  & /(inlcls(2,l1)+clnoise(2,l1)) &
1065  & *(inlcls(3,l1))*gte12_num
1066 
1067  else
1068 
1069  ftb12=prefactor*wigner2*incls(3,l1)*all*(all+1.d0)
1070  feb12=prefactor*wigner2*incls(2,l1)*all*(all+1.d0)
1071 
1072  aaltb = aaltb &
1073  & + ftb12**2/(inlcls(1,l1)+clnoise(1,l1)) &
1074  & /(inlcls(4,l1)+clnoise(2,l1))
1075  aaleb = aaleb &
1076  & + feb12**2/(inlcls(2,l1)+clnoise(2,l1)) &
1077  & /(inlcls(4,l1)+clnoise(2,l1))
1078 
1079  abltbeb = abltbeb &
1080  & + ftb12*feb12*(inlcls(3,l1)) &
1081  & /(inlcls(4,l1)+clnoise(2,l1)) &
1082  & /(inlcls(1,l1)+clnoise(1,l1)) &
1083  & /(inlcls(2,l1)+clnoise(2,l1))
1084 
1085  end if
1086  end if
1087 
1088  end do !! do 2
1089 
1090  ! finally computing the noise
1091 
1092  aaltt = all*(all+1.d0)*(2.d0*all+1.d0)/aaltt
1093  aalee = all*(all+1.d0)*(2.d0*all+1.d0)/aalee
1094  aalte = all*(all+1.d0)*(2.d0*all+1.d0)/aalte
1095  aaltb = all*(all+1.d0)*(2.d0*all+1.d0)/aaltb
1096  aaleb = all*(all+1.d0)*(2.d0*all+1.d0)/aaleb
1097  ablttee= all*(all+1.d0)*(2.d0*all+1.d0)/ablttee
1098  ablttte= all*(all+1.d0)*(2.d0*all+1.d0)/ablttte
1099  ableete= all*(all+1.d0)*(2.d0*all+1.d0)/ableete
1100  abltbeb= all*(all+1.d0)*(2.d0*all+1.d0)/abltbeb
1101 
1102  lensingnoise(1,1) = aaltt
1103  lensingnoise(2,2) = aalee
1104  lensingnoise(3,3) = aalte
1105  lensingnoise(4,4) = aaltb
1106  lensingnoise(5,5) = aaleb
1107 
1108  lensingnoise(1,2) = aaltt*aalee/ablttee
1109  lensingnoise(1,3) = aaltt*aalte/ablttte
1110  lensingnoise(1,4) = 0.
1111  lensingnoise(1,5) = 0.
1112  lensingnoise(2,3) = aalee*aalte/ableete
1113  lensingnoise(2,4) = 0.
1114  lensingnoise(2,5) = 0.
1115  lensingnoise(3,4) = 0.
1116  lensingnoise(3,5) = 0.
1117  lensingnoise(4,5) = aaltb*aaleb/abltbeb
1118 
1119  lensingnoise(2,1) = lensingnoise(1,2)
1120  lensingnoise(3,1) = lensingnoise(1,3)
1121  lensingnoise(4,1) = lensingnoise(1,4)
1122  lensingnoise(5,1) = lensingnoise(1,5)
1123  lensingnoise(3,2) = lensingnoise(2,3)
1124  lensingnoise(4,2) = lensingnoise(2,4)
1125  lensingnoise(5,2) = lensingnoise(2,5)
1126  lensingnoise(4,3) = lensingnoise(3,4)
1127  lensingnoise(5,3) = lensingnoise(3,5)
1128  lensingnoise(5,4) = lensingnoise(4,5)
1129 
1130  ! compute minimum variance estimator based on TT,EE,TE,TB,EB
1131 
1132  do l1=1,5
1133  gaussveclens(l1,1)=1.
1134  do l2=1,5
1135  invlensingnoise(l1,l2)=lensingnoise(l1,l2)
1136  end do
1137  end do
1138 
1139  call gaussj(invlensingnoise,5,5,gaussveclens,1,1)
1140 
1141  minvarnoise_l(jj)=0.
1142 
1143  do l1=1,5
1144  do l2=1,5
1145  minvarnoise_l(jj)=minvarnoise_l(jj)+invlensingnoise(l1,l2)
1146  end do
1147  end do
1148 
1149  minvarnoise_l(jj)=1.d0/minvarnoise_l(jj)
1150 
1151  ! compute minimum variance estimator based only on EE,EB
1152  !MinVarNoisePol_L(jj)=1.d0/
1153  !& (1.d0/LensingNoise(2,2)+1.d0/LensingNoise(5,5))
1154  ! compute minimum variance estimator based only on TT
1155  !MinVarNoiseTemp_L(jj)=LensingNoise(1,1)
1156 
1157  ! print the results
1158 
1159 ! write(37,'(1I5,8E13.4)')ll, &
1160 ! & LensingNoise(1,1)*all*(all+1.d0)/2.d0/pi, &
1161 ! & LensingNoise(2,2)*all*(all+1.d0)/2.d0/pi, &
1162 ! & LensingNoise(3,3)*all*(all+1.d0)/2.d0/pi, &
1163 ! & LensingNoise(4,4)*all*(all+1.d0)/2.d0/pi, &
1164 ! & LensingNoise(5,5)*all*(all+1.d0)/2.d0/pi, &
1165 ! & MinVarNoise_L(jj)*all*(all+1.d0)/2.d0/pi, &
1166 ! & incls(4,int(all))*all*(all+1.d0)/2.d0/pi
1167 
1168  alogllens(jj)=log(all)
1169  end do
1170 
1171  !close(37)
1172  ! interpolate to any value of l
1173  ! (interpolation performed in log(l) space)
1174  allocate(ddminvarnoise_l(1:jjmaxt))
1175  call spline(alogllens,minvarnoise_l,jjmaxt, &
1176  & 0.d0,0.d0,ddminvarnoise_l)
1177  !call spline(alogllens,MinVarNoisePol_L,jjmaxT,
1178  !& 0.d0,0.d0,ddMinVarNoisePol_L)
1179  !call spline(alogllens,MinVarNoiseTemp_L,jjmaxT,
1180  !& 0.d0,0.d0,ddMinVarNoiseTemp_L)
1181 ! open(unit=34,file="test_Nldd.dat",form="formatted", &
1182 ! & status='replace', iostat=err)
1183 ! if (err /= 0) print *, "Probleme d'ecriture du fichier test..."
1184  do l=2,lmaxm
1185  call splint(alogllens,minvarnoise_l,ddminvarnoise_l, &
1186  & jjmaxt,log(dble(l)),nldd(l))
1187  !call splint(alogllens,MinVarNoisePol_L,ddMinVarNoisePol_L,
1188  !& jjmaxT,log(dble(l)),omBPol_D(l))
1189  !call splint(alogllens,MinVarNoiseTemp_L,ddMinVarNoiseTemp_L,
1190  !& jjmaxT,log(dble(l)),omBTemp_D(l))
1191  !write(34,'(1I5,1E13.4)') l, nldd(l)*l*(l+1.d0)/2.d0/pi
1192  end do
1193  !close(34)
1194 
1195  !print*,code,"FIN"
1196  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1197  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1198  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1199  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1200 
1201  end subroutine calc_nldd
1202 
1203 
1204  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1205  SUBROUTINE gaussj(a,n,np,b,m,mp)
1206  !c implicit double precision (a-h,o-z)
1207  INTEGER m,mp,n,np,NMAX
1208  REAL*8 a(np,np),b(np,mp)
1209  parameter(nmax=50)
1210  !c Linear equation solution by Gauss-Jordan elimination,
1211  !c equation (2.1.1) above. a(1:n,1:n) is an input matrix stored in an
1212  !c array of physical dimensions np by np.
1213  !c b(1:n,1:m) is an in-put matrix containing the m right-hand side
1214  !c vectors, stored in an array of physical dimensions np by mp.
1215  !c On output, a(1:n,1:n) is replaced by its matrix inverse,
1216  !c and b(1:n,1:m) is replaced by the corresponding set of solution
1217  !c vectors. Parameter: NMAX is the largest anticipated value of n.
1218  INTEGER i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
1219  !c The integer arrays ipiv, indxr,andindxc are used for bookkeeping
1220  !c on the pivoting.
1221  REAL*8 big,dum,pivinv
1222 
1223  do 11 j=1,n
1224  ipiv(j)=0
1225 11 enddo
1226  do 22 i=1,n
1227  !c This is the main loop over the columns to be re-duced.
1228  big=0.
1229  do 13 j=1,n
1230  !c This is the outer loop of the search for a pivot ele-ment.
1231  if(ipiv(j).ne.1)then
1232  do 12 k=1,n
1233  if (ipiv(k).eq.0) then
1234  if (abs(a(j,k)).ge.big)then
1235  big=abs(a(j,k))
1236  irow=j
1237  icol=k
1238  endif
1239  endif
1240 12 enddo
1241  endif
1242 13 enddo
1243  ipiv(icol)=ipiv(icol)+1
1244  if (irow.ne.icol) then
1245  do 14 l=1,n
1246  dum=a(irow,l)
1247  a(irow,l)=a(icol,l)
1248  a(icol,l)=dum
1249 14 enddo
1250  do 15 l=1,m
1251  dum=b(irow,l)
1252  b(irow,l)=b(icol,l)
1253  b(icol,l)=dum
1254 15 enddo
1255  endif
1256  indxr(i)=irow
1257  !c We are now ready to divide the pivot row by the pivot element,
1258  !c located at irow and icol.
1259  indxc(i)=icol
1260  if (a(icol,icol).eq.0.) stop 'singular matrix in gaussj'
1261  !c singular matrix in gaussj
1262  pivinv=1./a(icol,icol)
1263  a(icol,icol)=1.
1264  do 16 l=1,n
1265  a(icol,l)=a(icol,l)*pivinv
1266 16 enddo
1267  do 17 l=1,m
1268  b(icol,l)=b(icol,l)*pivinv
1269 17 enddo
1270  do 21 ll=1,n
1271  !c Next, we reduce the rows...
1272  if(ll.ne.icol)then
1273  !c ...except for the pivot one, of course.
1274  dum=a(ll,icol)
1275  a(ll,icol)=0.
1276  do 18 l=1,n
1277  a(ll,l)=a(ll,l)-a(icol,l)*dum
1278 18 enddo
1279  do 19 l=1,m
1280  b(ll,l)=b(ll,l)-b(icol,l)*dum
1281 19 enddo
1282  endif
1283 21 enddo
1284 22 enddo
1285  !c This is the end of the main loop over columns of the reduction.
1286  do 24 l=n,1,-1
1287  !c It only remains to unscramble the solution in view of the
1288  !c column interchanges. We do this by in-terchanging pairs of
1289  !c columns in the reverse order that the permutation was built up.
1290  if(indxr(l).ne.indxc(l))then
1291  do 23 k=1,n
1292  dum=a(k,indxr(l))
1293  a(k,indxr(l))=a(k,indxc(l))
1294  a(k,indxc(l))=dum
1295 23 enddo
1296  endif
1297 24 enddo
1298  return
1299  !c And we are done.
1300  END SUBROUTINE gaussj
1301  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1302 
1303  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1304 
1305  SUBROUTINE spline(x,y,n,yp1,ypn,y2)
1306  INTEGER n,NMAX
1307  double precision yp1,ypn,x(n),y(n),y2(n)
1308  parameter(nmax=500)
1309  INTEGER i,k
1310  double precision p,qn,sig,un,u(nmax)
1311  if (yp1.gt..99d30) then
1312  y2(1)=0.d0
1313  u(1)=0.d0
1314  else
1315  y2(1)=-0.5d0
1316  u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
1317  endif
1318  do 11 i=2,n-1
1319  sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
1320  p=sig*y2(i-1)+2.d0
1321  y2(i)=(sig-1.d0)/p
1322  u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
1323  & /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
1324 11 continue
1325  if (ypn.gt..99d30) then
1326  qn=0.d0
1327  un=0.d0
1328  else
1329  qn=0.5d0
1330  un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
1331  endif
1332  y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
1333  do 12 k=n-1,1,-1
1334  y2(k)=y2(k)*y2(k+1)+u(k)
1335 12 continue
1336  return
1337  !C (C) Copr. 1986-92 Numerical Recipes Software =$j*m,).
1338  END subroutine spline
1339 
1340  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1341 
1342  subroutine splint(xa,ya,y2a,n,x,y)
1343  !c
1344  integer n
1345  double precision x,y,xa(n),y2a(n),ya(n)
1346  !c
1347  integer k,khi,klo
1348  double precision a,b,h
1349  !c
1350  klo=1
1351  khi=n
1352 1 if (khi-klo.gt.1) then
1353  k=(khi+klo)/2
1354  if(xa(k).gt.x)then
1355  khi=k
1356  else
1357  klo=k
1358  end if
1359  goto 1
1360  end if
1361  h=xa(khi)-xa(klo)
1362  if (h.eq.0.) stop 'bad xa input in splint'
1363  a=(xa(khi)-x)/h
1364  b=(x-xa(klo))/h
1365  y=a*ya(klo)+b*ya(khi) + &
1366  & ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
1367  return
1368  end subroutine splint
1369  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1370 
1371 end module lensnoise
lensnoise
This module contains the code to compute the noise for CMB lensing.
Definition: 007_FUTURCMB_lensnoise.f90:25
lensnoise::calc_nldd
subroutine, public calc_nldd(incls, inlcls, clnoise, nldd, jjmaxTmax, lmaxM)
This subroutine computes the noise for the Okamoto & Hu's estimator of the CMB lensing.
Definition: 007_FUTURCMB_lensnoise.f90:53
Generated on Tue Jun 21 2016 03:18:43 for CosmicFish by   doxygen 1.8.9.1