bianchi2_sky_mod.f90
Go to the documentation of this file.
1 !==============================================================================
16 !==============================================================================
17 
18 
19 !------------------------------------------------------------------------------
20 ! bianchi2_sky_mod -- BIANCHI2 library sky class
21 !
29 !------------------------------------------------------------------------------
30 
31 module bianchi2_sky_mod
32 
33  use s2_types_mod
34  use s2_sky_mod
35  use s2_error_mod
36  use bianchi2_error_mod
37  use bianchi2_lut_mod
38  use omp_lib
39 
40  implicit none
41 
42  private
43 
44 
45  !---------------------------------------
46  ! Subroutine and function scope
47  !---------------------------------------
48 
49  public :: &
50  bianchi2_sky_init, &
51  bianchi2_sky_init_alm, &
52  bianchi2_sky_free, &
53  bianchi2_sky_write, &
54  bianchi2_sky_rotate, &
55  bianchi2_sky_param_write, &
56  bianchi2_sky_compute_map, &
57  bianchi2_sky_compute_alm, &
58  bianchi2_sky_compute_cl, &
59  bianchi2_sky_get_sky, &
60  bianchi2_sky_get_alm, &
61  bianchi2_sky_write_cl, &
62  bianchi2_sky_apply_beam
63 
64 
65 
66  !---------------------------------------
67  ! Interfaces
68  !---------------------------------------
69 
70  ! None.
71 
72 
73  !---------------------------------------
74  ! Global variables
75  !---------------------------------------
76 
78 #ifdef MILLIK
79 
80  real(s2_dp), parameter :: BIANCHI2_CMB_T = 2.725d3
81 #else
82 
83  real(s2_dp), parameter :: BIANCHI2_CMB_T = 1d0
84 #endif
85 
86  !---------------------------------------
87  ! Data types
88  !---------------------------------------
89 
91  type, public :: bianchi2_sky
92  private
94  logical :: init = .false.
95 
96  real(s2_dp) :: omega_matter = 0.0d0
97 
98  real(s2_dp) :: omega_lambda = 0.0d0
99 
100  ! characteristic wavelength over which basis vectors change
101  ! orientation) (input parameter).
102  real(s2_dp) :: h = 0.0d0
103 
104  real(s2_dp) :: zE = 0.0d0
105 
106  ! parameter).
107  real(s2_dp) :: wH = 0.0d0
108 
109  logical :: rhand = .true.
110 
111  type(s2_sky) :: sky
112  end type bianchi2_sky
113 
114 
115  !----------------------------------------------------------------------------
116 
117  contains
118 
119 
120  !--------------------------------------------------------------------------
121  ! bianchi2_sky_init
122  !
140  !--------------------------------------------------------------------------
141 
142  function bianchi2_sky_init(omega_matter_in, omega_lambda_in, h, zE_in, wH, rhand, &
143  nside) result(b)
144 
145  use bianchi2_globaldata_mod
146  use pix_tools, only: pix2ang_ring, nside2npix, in_ring
147 
148  real(s2_dp), intent(in) :: omega_matter_in, omega_lambda_in, h, ze_in, wh
149  logical, intent(in) :: rhand
150  integer, intent(in) :: nside
151  type(bianchi2_sky) :: b
152 
153  real(s2_dp) :: thetaob, phiob ! Observing angles.
154  real(s2_dp) :: theta0, phi0 ! Barrow's angle convention.
155 
156  integer :: nt_use
157  real(s2_dp) :: tstop_use
158  real(s2_dp) :: tau_needed
159  real(s2_dp) :: fact, u10_req, phi1
160  real(s2_dp) :: r_final, rh_final, cos_bit_final, sin_bit_final, other_bit
161  real(s2_dp) :: lambda, final_dens, t0
162 
163  integer :: npix, ipix, fail
164  real(s2_dp), allocatable :: map(:)
165  real(s2_sp) :: handedness_sign = +1d0
166 
167  real(s2_dp), parameter :: phi_centre = 0.0d0 ! Extract extire rings.
168  real(s2_dp), parameter :: phi_diff = pi
169  integer :: iring, nring, max_pix_per_ring, npixring
170  integer, allocatable :: ipixring(:)
171 
172  ! Check object not already initialised.
173  if(b%init) then
174  call bianchi2_error(bianchi2_error_init, 'bianchi2_sky_init')
175  return
176  end if
177 
179  b%omega_matter = omega_matter_in
180  b%omega_lambda = omega_lambda_in
181  b%h = h
182  b%zE = ze_in
183  b%wH = wh
184  b%rhand = rhand
185 
186  ! Initialise healpix map settings.
187  npix = nside2npix(nside)
188  allocate(map(0:npix-1), stat=fail)
189  map = 0d0
190  if(fail /= 0) then
191  call bianchi2_error(bianchi2_error_mem_alloc_fail, &
192  'bianchi2_sky_init')
193  end if
194 
196  if(b%rhand) then
197  handedness_sign = 1d0
198  else
199  handedness_sign = -1d0
200  end if
201 
202 ! ! Allocate global data.
203 ! allocate(b2gd_treal(1:b2gd_ntmax), stat=fail)
204 ! allocate(b2gd_xreal(1:b2gd_ntmax), stat=fail)
205 ! allocate(b2gd_tarr(1:b2gd_ntmax), stat=fail)
206 ! allocate(b2gd_xarr(1:b2gd_nvars, 1:b2gd_ntmax), stat=fail)
207 ! if(fail /= 0) then
208 ! call bianchi2_error(BIANCHI2_ERROR_MEM_ALLOC_FAIL, &
209 ! 'bianchi2_sky_init')
210 ! end if
211 
213  b2gd_bianchi_h = b%h
214  b2gd_omega_matter = b%omega_matter
215  b2gd_omega_lambda = b%omega_lambda
216  b2gd_ze = b%zE
217 
218  ! Initialise other variables.
219  nt_use = 3
220  b2gd_alpha = sqrt(b%h)
221  b2gd_rh_start = sqrt(-b2gd_alpha**2/(b%omega_matter+b%omega_lambda-1d0))
222  call get_tau(tau_needed)
223  tstop_use=-tau_needed
224  lambda = 3*b2gd_rh_start**2*b%omega_lambda
225 
226  ! Set ring parameter values.
227  nring = 4*nside - 1
228  max_pix_per_ring = 4*nside
229 
230  ! Allocate space for ring pixel indices.
231  allocate(ipixring(0:max_pix_per_ring-1), stat=fail)
232  if(fail /= 0) then
233  call bianchi2_error(bianchi2_error_mem_alloc_fail, &
234  'bianchi2_sky_init')
235  end if
236  ipixring = 0
237 
238  ! Compute map value for each pixel in each ring.
239  ! Do ring at a time since A and B functions of theta so once need
240  ! to compute once for each ring.
241 
242  do iring = 1, nring
243 
244 ! if(mod(iring,nring/4) == 0) then
245 ! write(*,'(a,f5.1,a)') ' Percent complete: ', &
246 ! iring/real(nring,s2_sp)*100e0, '%'
247 ! end if
248 
249  ! Get ring_pix in order to determine theta of ring.
250  call in_ring(nside, iring, phi_centre, phi_diff, ipixring, &
251  & npixring, nest=0) ! Ring scheme!!!!!!!!!
252 
253  ! Calculate thetaOB for current ring.
254  call pix2ang_ring(nside, ipixring(0), thetaob, phiob)
255 
256  ! Set photon theta angle from observation angle.
257  theta0 = pi - thetaob
258 
259  ! Set global data.
260  ! (Compute theta_0 for ANL's convention.)
261  b2gd_theta_0 = theta0 - pi/2d0
262 
263  ! Compute terms for constant theta.
264  call get_results(tstop_use,r_final,rh_final,cos_bit_final,sin_bit_final)
265  fact = sqrt(9.d0*b2gd_alpha**2+1.d0)*sqrt(1.d0+b2gd_alpha**2)* &
266  sqrt(1.d0-b%omega_matter-b%omega_lambda)**3.d0/b2gd_alpha**3 &
267  /b%omega_matter/6.d0
268  u10_req = b%wH/fact
269 
270  ! Compute map valus for all pixels in current ring.
271  do ipix = 0,npixring-1
272 
273  ! Calculate thetaOB for current ring.
274  call pix2ang_ring(nside, ipixring(ipix), thetaob, phiob)
275 
276  ! Set photon phi angle from observation angle.
277  phi0 = phiob - pi
278 
279  ! Account for handedness
280  ! (also component outside loop to account for handedness).
281  phi0 = phi0 * handedness_sign
282 
283  ! Compute phi_0 for ANL's convention.
284  b2gd_phi_0 = phi0 + 3d0/4d0*pi
285 
286  ! Get first contribution.
287  phi1 = u10_req/r_final*(sin_bit_final*sin(b2gd_phi_0)+cos_bit_final*cos(b2gd_phi_0))
288  map(ipixring(ipix)) = -phi1/(1+b2gd_ze)*2d0 ! JDM: Need factor 2 to make agree
289 
290  ! Now get end contribution.
291  final_dens = (3d0*rh_final**2-lambda*r_final**2-3d0*b2gd_alpha**2)/(r_final**2)
292  t0 = 3.d0/final_dens/r_final**5*cos((b2gd_phi_0*b2gd_alpha-b2gd_alpha*tstop_use- &
293  dlog(cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0+exp(-2.d0*b2gd_alpha*tstop_use)- &
294  exp(-2.d0*b2gd_alpha*tstop_use)*cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0))/b2gd_alpha)* &
295  exp(b2gd_alpha*tstop_use)*cos(b2gd_theta_0)/(-1.d0-sin(b2gd_theta_0)+ &
296  exp(2.d0*b2gd_alpha*tstop_use)*sin(b2gd_theta_0)-exp(2.d0*b2gd_alpha*tstop_use))*b2gd_alpha+ &
297  1.d0/final_dens/r_final**5*sin((b2gd_phi_0*b2gd_alpha-b2gd_alpha*tstop_use- &
298  dlog(cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0+exp(-2.d0*b2gd_alpha*tstop_use)- &
299  exp(-2.d0*b2gd_alpha*tstop_use)*cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0))/b2gd_alpha)* &
300  exp(b2gd_alpha*tstop_use)*cos(b2gd_theta_0)/(-1.d0-sin(b2gd_theta_0)+ &
301  exp(2.d0*b2gd_alpha*tstop_use)*sin(b2gd_theta_0)-exp(2.d0*b2gd_alpha*tstop_use))
302  other_bit = u10_req*t0
303 
304  map(ipixring(ipix)) = map(ipixring(ipix)) - other_bit/(1+b2gd_ze)*2d0
305 
306  end do
307 
308  end do
309 
310 ! write(*,'(a)') ' Percent complete: 100.0%'
311 
312  ! Account for handedness.
313  map = handedness_sign * map
314 
315  ! Convert Delta_T/T map computed to Delta_T map.
316  map = bianchi2_cmb_t * map
317 
318  ! Initialise sky object with map.
319  b%sky = s2_sky_init(real(map,s2_sp), nside, s2_sky_ring)
320 
321  ! Set initialised status.
322  b%init = .true.
323 
324 ! ! Free global data.
325 ! deallocate(b2gd_treal)
326 ! deallocate(b2gd_xreal)
327 ! deallocate(b2gd_tarr)
328 ! deallocate(b2gd_xarr)
329 
330  ! Free memory.
331  deallocate(map)
332  deallocate(ipixring)
333 
334  end function bianchi2_sky_init
335 
336 
337 
338  !--------------------------------------------------------------------------
339  ! bianchi2_sky_init_alm
340  !
364  !--------------------------------------------------------------------------
365 
366  function bianchi2_sky_init_alm(omega_matter_in, omega_lambda_in, h, zE_in, wH, rhand, &
367  nside, lmax, ntheta, alpha, beta, gamma, lut) result(b)
368 
369  use bianchi2_globaldata_mod
370  use s2_dl_mod, only : s2_dl_beta_operator
371 
372  real(s2_dp), intent(in) :: omega_matter_in, omega_lambda_in, h, ze_in, wh
373  logical, intent(in) :: rhand
374  integer, intent(in) :: nside, lmax, ntheta
375  type(bianchi2_lut), intent(in), optional :: lut
376  type(bianchi2_sky) :: b
377 
378  real(s2_dp) :: tstop_use
379  real(s2_dp) :: tau_needed
380  real(s2_dp) :: fact, u10_req
381  real(s2_dp) :: r_final, rh_final, cos_bit_final, sin_bit_final
382  real(s2_dp) :: lambda, final_dens
383 
384  integer :: npix, ipix, fail
385 
386  ! Parameters for computing alm.
387  complex(s2_spc), allocatable :: alm(:,:)
388  real(s2_dp) :: handedness_sign = +1d0, lsign = 1d0
389  integer :: l, itheta
390  real(s2_dp), allocatable :: a_grid(:), b_grid(:)
391  real(s2_dp) :: theta_inc
392  real(s2_dp) :: c_sin, c_cos
393  real(s2_dp) :: ia, ib
394 
395  ! Parameters for the rotation.
396  real(s2_sp), intent(in), optional :: alpha, beta, gamma
397  real(s2_sp), parameter :: zero_tol = 1d-4
398  real(s2_dp), pointer :: dl(:,:) => null()
399  integer :: m
400  complex(s2_spc), allocatable :: alm_rotate(:,:)
401  complex(s2_dpc) :: dm_p1, dm_m1
402  complex(s2_dpc) :: icmpx
403 
404  ! Set icmpx=sqrt(-1).
405  icmpx = cmplx(0d0,1d0)
406 
407  ! Check object not already initialised.
408  if(b%init) then
409  call bianchi2_error(bianchi2_error_init, 'bianchi2_sky_init_alm')
410  return
411  end if
412 
413  ! Initialise parameters passed as arguments.
414  b%omega_matter = omega_matter_in
415  b%omega_lambda = omega_lambda_in
416  b%h = h
417  b%zE = ze_in
418  b%wH = wh
419  b%rhand = rhand
420 
421  ! Initialise healpix alm settings.
422  allocate(alm(0:lmax,0:lmax), stat=fail)
423  alm = cmplx(0d0, 0d0)
424  if(fail /= 0) then
425  call bianchi2_error(bianchi2_error_mem_alloc_fail, &
426  'bianchi2_sky_init_alm')
427  end if
428 
429  ! Set handedness sign.
430  if(b%rhand) then
431  handedness_sign = 1d0
432  else
433  handedness_sign = -1d0
434  end if
435 
436  ! Set global data.
437  b2gd_bianchi_h = b%h
438  b2gd_omega_matter = b%omega_matter
439  b2gd_omega_lambda = b%omega_lambda
440  b2gd_ze = b%zE
441 
442  ! Initialise other variables.
443  b2gd_alpha = sqrt(b%h)
444  b2gd_rh_start = sqrt(-b2gd_alpha**2/(b%omega_matter+b%omega_lambda-1d0))
445  call get_tau(tau_needed)
446  tstop_use=-tau_needed
447 
448  ! Calculate A(theta) and B(theta) terms.
449  allocate(a_grid(0:ntheta-1), stat=fail)
450  allocate(b_grid(0:ntheta-1), stat=fail)
451 
452  if(fail /= 0) then
453  call bianchi2_error(bianchi2_error_mem_alloc_fail, &
454  'bianchi2_sky_init_alm')
455  end if
456 
457  a_grid = 0d0
458  b_grid = 0d0
459  b2gd_theta_0 = -pi/2d0
460  theta_inc = (pi - 0d0) / real(ntheta, s2_dp)
461 
462  !$OMP PARALLEL DEFAULT(none), COPYIN(b2gd_it, b2gd_nt, b2_gd_treal,b2gd_xreal,b2gd_theta_0,b2gd_xarr,b2gd_tarr) &
463  !$OMP FIRSTPRIVATE(itheta,R_final,RH_final,cos_bit_final,sin_bit_final,final_dens,C_sin,C_cos) &
464  !$OMP SHARED(Ntheta,A_grid,B_grid,tstop_use,theta_inc,Lambda,fact,U10_req,b,b2gd_RH_start,b2gd_deltat,b2gd_ze,b2gd_tstop,b2gd_alpha,b2gd_Omega_matter,b2gd_Omega_Lambda)
465  !$OMP DO SCHEDULE (static)
466  do itheta = 0, ntheta-1
467 
468  b2gd_theta_0=-pi/2d0+itheta*theta_inc
469 
470  call get_results(tstop_use,r_final,rh_final,cos_bit_final,sin_bit_final)
471  fact = sqrt(9.d0*b2gd_alpha**2+1.d0)*sqrt(1.d0+b2gd_alpha**2)* &
472  sqrt(1.d0-b%omega_matter-b%omega_lambda)**3.d0/b2gd_alpha**3 &
473  /b%omega_matter/6.d0
474  u10_req = b%wH/fact
475  lambda = 3*b2gd_rh_start**2*b%omega_lambda
476  final_dens = (3d0*rh_final**2-lambda*r_final**2-3d0*b2gd_alpha**2)/(r_final**2)
477 
478  c_sin=-u10_req/(1+b2gd_ze)*2d0*( &
479  1d0/r_final*(sin_bit_final*cos(3d0/4d0*pi)-cos_bit_final*sin(3d0/4d0*pi)) &
480  - 3.d0/final_dens/r_final**5 &
481  * exp(b2gd_alpha*tstop_use)*cos(b2gd_theta_0) &
482  / (-1.d0-sin(b2gd_theta_0) &
483  + exp(2.d0*b2gd_alpha*tstop_use)*sin(b2gd_theta_0) &
484  - exp(2.d0*b2gd_alpha*tstop_use)) &
485  * b2gd_alpha &
486  * sin((-b2gd_alpha*tstop_use- &
487  dlog(cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0 &
488  + exp(-2.d0*b2gd_alpha*tstop_use) &
489  - exp(-2.d0*b2gd_alpha*tstop_use)*cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0) &
490  )/b2gd_alpha+3d0/4d0*pi) &
491  + 1.d0/final_dens/r_final**5 &
492  * exp(b2gd_alpha*tstop_use)*cos(b2gd_theta_0) &
493  / (-1.d0-sin(b2gd_theta_0) &
494  + exp(2.d0*b2gd_alpha*tstop_use)*sin(b2gd_theta_0) &
495  -exp(2.d0*b2gd_alpha*tstop_use)) &
496  * cos((-b2gd_alpha*tstop_use- &
497  dlog(cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0 &
498  + exp(-2.d0*b2gd_alpha*tstop_use) &
499  - exp(-2.d0*b2gd_alpha*tstop_use)*cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0) &
500  )/b2gd_alpha+3d0/4d0*pi) &
501  )
502 
503  c_cos=-u10_req/(1+b2gd_ze)*2d0*( &
504  1d0/r_final*(sin_bit_final*sin(3d0/4d0*pi)+cos_bit_final*cos(3d0/4d0*pi)) &
505  + 3.d0/final_dens/r_final**5 &
506  * exp(b2gd_alpha*tstop_use)*cos(b2gd_theta_0) &
507  / (-1.d0-sin(b2gd_theta_0) &
508  + exp(2.d0*b2gd_alpha*tstop_use)*sin(b2gd_theta_0) &
509  - exp(2.d0*b2gd_alpha*tstop_use)) &
510  * b2gd_alpha &
511  * cos((-b2gd_alpha*tstop_use- &
512  dlog(cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0 &
513  + exp(-2.d0*b2gd_alpha*tstop_use) &
514  - exp(-2.d0*b2gd_alpha*tstop_use)*cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0) &
515  )/b2gd_alpha+3d0/4d0*pi) &
516  + 1.d0/final_dens/r_final**5 &
517  * exp(b2gd_alpha*tstop_use)*cos(b2gd_theta_0) &
518  / (-1.d0-sin(b2gd_theta_0) &
519  + exp(2.d0*b2gd_alpha*tstop_use)*sin(b2gd_theta_0) &
520  -exp(2.d0*b2gd_alpha*tstop_use)) &
521  * sin((-b2gd_alpha*tstop_use- &
522  dlog(cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0 &
523  + exp(-2.d0*b2gd_alpha*tstop_use) &
524  - exp(-2.d0*b2gd_alpha*tstop_use)*cos(pi/4.d0+b2gd_theta_0/2.d0)**2.d0) &
525  )/b2gd_alpha+3d0/4d0*pi) &
526  )
527 
528  a_grid(itheta) = (c_sin - c_cos) / 2d0
529  b_grid(itheta) = (c_sin + c_cos) / 2d0
530 
531  end do
532  !$OMP END DO
533  !$OMP END PARALLEL
534 
535  ! Compute alms.
536  lsign = +1d0
537 
538  !$OMP PARALLEL DEFAULT(none) &
539  !$OMP FIRSTPRIVATE(l,lsign,IA,IB) &
540  !$OMP SHARED(lmax,alm,Ntheta,A_grid,B_grid,handedness_sign,lut)
541  !$OMP DO SCHEDULE(static)
542  do l=1, lmax
543 
544  ! Invert lsign.
545  if (l/2==l/2.d0) then
546  lsign=+1d0
547  else
548  lsign=-1d0
549  end if
550 
551  ! Compute integrals.
552  ! Use precomputed A(theta) and B(theta).
553  ia = bianchi2_sky_comp_ix(l,lmax,ntheta,a_grid,lut)
554  ib = bianchi2_sky_comp_ix(l,lmax,ntheta,b_grid,lut)
555 
556  ! Compute alm for a given 1. Only m=1 is non-zero.
557  alm(l,1) = -lsign * pi * cmplx(- handedness_sign * (ib - ia), (ia + ib))
558 
559  end do
560  !$OMP END DO
561  !$OMP END PARALLEL
562 
563 ! write(*,'(a)') ' Percent complete: 100.0%'
564 
565  ! Convert Delta_T/T alm computed to Delta_T alm.
566  alm = bianchi2_cmb_t * alm
567 
568  ! Rotate alms if Euler angles present and ar least one angle non-zero.
569  if(present(alpha) .and. present(beta) .and. present(gamma) &
570  .and. (abs(alpha)+abs(beta)+abs(gamma) > zero_tol) ) then
571 
572 ! write(*,*) ' Rotation of the alm '
573  allocate(dl(-lmax:lmax,-lmax:lmax),stat=fail)
574  allocate(alm_rotate(0:lmax,0:lmax), stat=fail)
575  if(fail/=0) then
576  call bianchi2_error(bianchi2_error_mem_alloc_fail, &
577  'bianchi_sky_rotate')
578  endif
579  alm_rotate = (0d0,0d0)
580 
581 
582  do l = 1, lmax
583 
584  ! Computing the Wigner functions is the bottleneck for
585  ! computation time. If necessary, we could optimise this
586  ! by precomputing dlmn(pi/2) and then using FFTs to
587  ! compute dlmn(beta).
588  call s2_dl_beta_operator(dl, real(beta,s2_dp), l)
589 
590  do m = 0,l
591 
592  ! Calculation of Dl,m,+-1.
593  dm_p1 = exp(-icmpx*m*alpha) * dl(m, 1) * exp(-icmpx*gamma)
594  dm_m1 = exp(-icmpx*m*alpha) * dl(m,-1) * exp( icmpx*gamma)
595 
596  ! Rotation of the alm.
597  alm_rotate(l,m) = - dm_m1*conjg(alm(l,1)) + dm_p1*alm(l,1)
598 
599  end do
600  end do
601 
602  ! Initialise sky object with alm_rotate.
603  b%sky = s2_sky_init(alm_rotate, lmax, lmax, nside)
604  deallocate(alm_rotate)
605  deallocate(dl)
606 
607  else
608 
609 ! write(*,'(a)') ' No rotation needed '
610 
611  ! Initialise sky object with alm.
612  b%sky = s2_sky_init(alm, lmax, lmax, nside)
613  endif
614 
615  ! Set initialised status.
616  b%init = .true.
617 
618 
619  ! call bianchi2_sky_compute_map(b,nside)
620 
621  ! Free memory.
622  deallocate(alm)
623  deallocate(a_grid)
624  deallocate(b_grid)
625 
626  end function bianchi2_sky_init_alm
627 
628 
629 
630  !--------------------------------------------------------------------------
631  ! bianchi2_sky_free
632  !
638  !--------------------------------------------------------------------------
639 
640  subroutine bianchi2_sky_free(b)
641 
642  type(bianchi2_sky), intent(inout) :: b
643 
644  ! Check object initialised.
645  if(.not. b%init) then
646  call bianchi2_error(bianchi2_error_not_init, 'bianchi2_sky_free')
647  end if
648 
649  ! Free sky.
650  call s2_sky_free(b%sky)
651 
652  ! Reset attributes.
653  b%omega_matter = 0.0d0
654  b%omega_lambda = 0.0d0
655  b%h = 0.0d0
656  b%zE = 0.0d0
657  b%wH = 0.0d0
658  b%rhand = .true.
659  b%init = .false.
660 
661  end subroutine bianchi2_sky_free
662 
663 
664  !--------------------------------------------------------------------------
665  ! bianchi2_sky_param_write
666  !
673  !--------------------------------------------------------------------------
674 
675  subroutine bianchi2_sky_param_write(b)
676 
677  type(bianchi2_sky), intent(in) :: b
678 
679  write(*,'(a,e12.5)') ' b%omega_matter: ', b%omega_matter
680  write(*,'(a,e12.5)') ' b%omega_lambda: ', b%omega_lambda
681  write(*,'(a,e23.5)') ' b%h: ', b%h
682  write(*,'(a,e22.5)') ' b%zE: ', b%zE
683  write(*,'(a,e22.5)') ' b%wH: ', b%wH
684  write(*,'(a,l19)') ' b%rhand: ', b%rhand
685 
686  end subroutine bianchi2_sky_param_write
687 
688 
689  !--------------------------------------------------------------------------
690  ! bianchi2_sky_compute_map
691  !
699  !--------------------------------------------------------------------------
700 
701  subroutine bianchi2_sky_compute_map(b, nside)
702 
703  type(bianchi2_sky), intent(inout) :: b
704  integer, intent(in) :: nside
705 
706  ! Check object initialised.
707  if(.not. b%init) then
708  call bianchi2_error(bianchi2_error_not_init, &
709  'bianchi2_sky_compute_map')
710  end if
711 
712  ! Compute bianchi2 sky alms.
713  call s2_sky_compute_map(b%sky, nside)
714 
715  end subroutine bianchi2_sky_compute_map
716 
717 
718  !--------------------------------------------------------------------------
719  ! bianchi2_sky_compute_alm
720  !
729  !--------------------------------------------------------------------------
730 
731  subroutine bianchi2_sky_compute_alm(b, lmax, mmax)
732 
733  type(bianchi2_sky), intent(inout) :: b
734  integer, intent(in) :: lmax, mmax
735 
736  ! Check object initialised.
737  if(.not. b%init) then
738  call bianchi2_error(bianchi2_error_not_init, &
739  'bianchi2_sky_compute_alm')
740  end if
741 
742  ! Compute bianchi2 sky alms.
743  call s2_sky_compute_alm(b%sky, lmax, mmax)
744 
745  end subroutine bianchi2_sky_compute_alm
746 
747 
748  !--------------------------------------------------------------------------
749  ! bianchi2_sky_compute_Cl
750  !
761  !--------------------------------------------------------------------------
762 
763  subroutine bianchi2_sky_compute_cl(b, lmax, Cl)
764 
765  type(bianchi2_sky), intent(in) :: b
766  integer, intent(in) :: lmax
767  real(s2_sp), dimension(0:lmax), intent(out) :: cl
768 
769  complex(s2_spc), dimension(0:lmax,0:lmax) :: alm
770  integer :: l,m
771 
772  ! Check object initialised.
773  if(.not. b%init) then
774  call bianchi2_error(bianchi2_error_not_init, 'bianchi2_sky_get_Cl')
775  end if
776 
777  ! Get alms from sky.
778  call s2_sky_get_alm(b%sky, alm)
779 
780  ! Compute the Cls.
781  do l=0, lmax
782  cl(l)=real(alm(l,0)*conjg(alm(l,0)))
783  do m=1, l
784  cl(l)=cl(l)+2*real(alm(l,m)*conjg(alm(l,m)))
785  end do
786  cl(l)=cl(l)/(2*l+1)
787  end do
788 
789  end subroutine bianchi2_sky_compute_cl
790 
791 
792  !--------------------------------------------------------------------------
793  ! bianchi2_sky_get_sky
794  !
805  !--------------------------------------------------------------------------
806 
807  function bianchi2_sky_get_sky(b) result(sky)
808 
809  type(bianchi2_sky), intent(in) :: b
810  type(s2_sky) :: sky
811 
812  ! Check object initialised.
813  if(.not. b%init) then
814  call s2_error(bianchi2_error_not_init, 'bianchi2_sky_get_sky')
815  end if
816 
817  ! Make a copy for the returned sky.
818  sky = s2_sky_init(b%sky)
819 
820  end function bianchi2_sky_get_sky
821 
822 
823  !--------------------------------------------------------------------------
824  ! bianchi2_sky_get_alm
825  !
836  !--------------------------------------------------------------------------
837 
838  subroutine bianchi2_sky_get_alm(b, alm)
839 
840  type(bianchi2_sky), intent(in) :: b
841  complex(s2_spc), intent(out) :: alm(:,:)
842 
843  ! Check object initialised.
844  if(.not. b%init) then
845  call bianchi2_error(bianchi2_error_not_init, 'bianchi2_sky_get_alm')
846  end if
847 
848  ! Get alms from sky.
849  call s2_sky_get_alm(b%sky, alm)
850 
851  end subroutine bianchi2_sky_get_alm
852 
853 
854  !--------------------------------------------------------------------------
855  ! bianchi2_sky_apply_beam
856  !
868  !--------------------------------------------------------------------------
869 
870  subroutine bianchi2_sky_apply_beam(b, fwhm, lmax)
871 
872  use s2_pl_mod
873  use s2_vect_mod, only: s2_vect_arcmin_to_rad
874 
875  type(bianchi2_sky), intent(inout) :: b
876  real(s2_sp), intent(in) :: fwhm
877  integer, intent(in) :: lmax
878 
879  real(s2_sp) :: fwhm_rad
880  type(s2_pl) :: beam
881 
882  ! Convert beam_fwhm to radians.
883  fwhm_rad = s2_vect_arcmin_to_rad(fwhm)
884 
885  ! Create beam.
886  beam = s2_pl_init_guassian(fwhm_rad, lmax)
887 
888  ! Apply beam.
889  call s2_sky_conv(b%sky, beam)
890 
891  ! Free beam.
892  call s2_pl_free(beam)
893 
894  end subroutine bianchi2_sky_apply_beam
895 
896 
897  !--------------------------------------------------------------------------
898  ! bianchi2_sky_write
899  !
910  !--------------------------------------------------------------------------
911 
912  subroutine bianchi2_sky_write(b, filename, file_type, comment)
913 
914  type(bianchi2_sky), intent(inout) :: b
915  character(len=*), intent(in) :: filename
916  integer, intent(in) :: file_type
917  character(len=*), intent(in), optional :: comment
918 
919  ! Check object initialised.
920  if(.not. b%init) then
921  call bianchi2_error(bianchi2_error_not_init, 'bianchi2_sky_write')
922  end if
923 
924  select case(file_type)
925 
926  case(s2_sky_file_type_map)
927  call s2_sky_write_map_file(b%sky, filename, comment)
928 
929  case(s2_sky_file_type_alm)
930  call s2_sky_write_alm_file(b%sky, filename, comment)
931 
932  case(s2_sky_file_type_sky)
933  call s2_sky_io_fits_write(filename, b%sky, comment)
934 
935  case default
936  call s2_error(s2_error_sky_file_invalid, 'bianchi2_sky_write', &
937  comment_add='Invalid file type specifier')
938 
939  end select
940 
941  end subroutine bianchi2_sky_write
942 
943 
944  !--------------------------------------------------------------------------
945  ! bianchi2_sky_write_Cl
946  !
955  !--------------------------------------------------------------------------
956 
957  subroutine bianchi2_sky_write_cl (Cl, lmax, rescale_Cl, filename)
958 
959  integer, intent(in) :: lmax
960  real(s2_sp), dimension(0:lmax), intent(inout) :: cl
961  logical, intent(in) :: rescale_cl
962  character(len=*), intent(in) :: filename
963 
964  integer :: l
965 
966  open(11,file=filename)
967  if (rescale_cl) then
968  do l=0, lmax
969  write(11,*) l, cl(l)*l*(l+1)/ (2d0*pi)
970  end do
971  else
972  do l=0, lmax
973  write(11,*) l, cl(l)
974  end do
975  end if
976  close(11)
977 
978  end subroutine bianchi2_sky_write_cl
979 
980 
981  !--------------------------------------------------------------------------
982  ! bianchi2_sky_rotate
983  !
997  !--------------------------------------------------------------------------
998 
999  subroutine bianchi2_sky_rotate(b, alpha, beta, gamma, lmax, nside, &
1000  rotation_alm)
1001 
1002  use s2_dl_mod, only: s2_dl_beta_operator
1003 
1004  type(bianchi2_sky), intent(inout) :: b
1005  real(s2_sp), intent(in), optional :: alpha, beta, gamma
1006 
1007  logical, intent(in) :: rotation_alm
1008  integer, intent(in) :: lmax, nside
1009  integer :: l, m
1010  integer :: fail=10
1011  complex(s2_spc), allocatable :: alm(:,:)
1012  complex(s2_spc), allocatable :: alm_rotate(:,:)
1013  real(s2_dp), pointer :: dl(:,:) => null()
1014  complex(s2_dpc) :: dm_p1, dm_m1
1015  real(s2_sp), parameter :: zero_tol = 1d-4
1016  complex(s2_dpc) :: icmpx
1017 
1018  ! Set icmpx=sqrt(-1).
1019  icmpx = cmplx(0d0,1d0)
1020 
1021  ! Check object initialised.
1022  if(.not. b%init) then
1023  call bianchi2_error(bianchi2_error_not_init, 'bianchi2_sky_rotate')
1024  end if
1025 
1026 
1027  ! Rotate if Euler angles present and at least one angle non-zero.
1028  if(present(alpha) .and. present(beta) .and. present(gamma) &
1029  .and. (abs(alpha)+abs(beta)+abs(gamma) > zero_tol) ) then
1030 
1031  ! Rotation pixel by pixel or rotation of the alm.
1032  if (.not. rotation_alm) then
1033 
1034  write(*,'(a)') ' Rotation pixel by pixel with s2_sky_rotate '
1035  call s2_sky_rotate(b%sky, alpha, beta, gamma)
1036 
1037  else
1038 
1039  write(*,'(a)') ' Rotation of the alm '
1040  allocate(dl(-lmax:lmax,-lmax:lmax),stat=fail)
1041  allocate(alm(0:lmax,0:1), stat=fail)
1042  allocate(alm_rotate(0:lmax,0:lmax), stat=fail)
1043  if(fail/=0) then
1044  call bianchi2_error(bianchi2_error_mem_alloc_fail, &
1045  'bianchi_sky_rotate')
1046  endif
1047  alm_rotate = 0e0
1048 
1049  ! Get alm.
1050  call bianchi2_sky_compute_alm(b,lmax,1)
1051  call bianchi2_sky_get_alm(b, alm)
1052 
1053  ! Perform rotation in harmonic space.
1054  ! Noting Bianchi alms only non zero for m=+-1).
1055  do l = 1,lmax
1056 
1057  call s2_dl_beta_operator(dl, real(beta,s2_dp), l)
1058 
1059  do m = 0,l
1060 
1061 
1062  ! Calculation of Dl,m,+-1.
1063  dm_p1 = exp(-icmpx*m*alpha) * dl(m, 1) * exp(-icmpx*gamma)
1064  dm_m1 = exp(-icmpx*m*alpha) * dl(m,-1) * exp( icmpx*gamma)
1065 
1066  ! Rotation of the alm.
1067  alm_rotate(l,m) = - dm_m1*conjg(alm(l,1)) + dm_p1*alm(l,1)
1068 
1069  end do
1070 
1071  end do
1072 
1073  ! Make free b%sky.
1074  call s2_sky_free(b%sky)
1075 
1076  ! Compute sky object with rotated alms.
1077  b%sky = s2_sky_init(alm_rotate, lmax, lmax, nside)
1078 
1079  call bianchi2_sky_compute_map(b,nside)
1080 
1081  ! Dellocate memory used for rotating alms.
1082  deallocate(alm_rotate)
1083  deallocate(alm)
1084  deallocate(dl)
1085 
1086  end if
1087 
1088  else
1089  write(*,'(a)') ' No rotation needed '
1090  endif
1091 
1092  end subroutine bianchi2_sky_rotate
1093 
1094 
1095  !--------------------------------------------------------------------------
1096  ! Routines to compute Bianchi induced fluctuations
1097  !--------------------------------------------------------------------------
1098 
1099  !--------------------------------------------------------------------------
1100  ! get_results
1101  !
1113  !--------------------------------------------------------------------------
1114 
1115  subroutine get_results(tstop_use, R_final, RH_final, &
1116  cos_bit_final, sin_bit_final)
1117 
1118  use s2_types_mod
1119  use bianchi2_globaldata_mod
1120 
1121  implicit none
1122 
1123  real(s2_dp), intent(in) :: tstop_use
1124  real(s2_dp), intent(out) :: r_final,rh_final,cos_bit_final,sin_bit_final
1125 
1126  integer, parameter :: neq = 4
1127 
1128  real(s2_dp) :: vals(b2gd_nvars)
1129  real(s2_dp) :: warr(21*neq+28)
1130  real(s2_dp) :: tstart,tol
1131  real(s2_dp) :: tmin,tmax,t,tend
1132  real(s2_dp) :: r,rh
1133 #ifdef NAGI8
1134  integer(kind=8) :: neqn,irelab,ifail
1135 #else
1136  integer :: neqn,irelab,ifail=1
1137 #endif
1138 
1139  character(len=1) :: relabs
1140  external d02cjw
1141 
1142  b2gd_tstop=tstop_use
1143  tstart=0
1144 
1145  b2gd_nt=3
1146 
1147  b2gd_deltat=(b2gd_tstop-tstart)/(b2gd_nt-1)
1148 
1149  tmin=tstart
1150  tmax=b2gd_tstop
1151 
1152  tol=1d-12
1153 
1154  t=tstart
1155 
1156  ! Next bit sets up initial value for R.
1157 
1158  vals(1)=1
1159 
1160  ! Next bit sets up initial value for RH.
1161 
1162  vals(2)=b2gd_rh_start
1163 
1164  vals(3)=0d0
1165  vals(4)=0d0
1166 
1167  t=tstart
1168  tend=b2gd_tstop
1169  neqn=neq
1170  irelab=0
1171  ifail=1
1172  b2gd_it=b2gd_nt-2
1173  relabs='D'
1174 
1175  call d02cjf(t,tend,neqn,vals,fcn,tol,relabs,out,d02cjw,warr,ifail)
1176 
1177  ! Check success.
1178  if((ifail/=0).and.(ifail/=1)) then
1179  call bianchi2_error(bianchi2_error_sky_num_fail, 'get_results', &
1180  comment_add='NAG routine d02cjf failed')
1181  end if
1182 
1183  do b2gd_it=1,b2gd_nt
1184  b2gd_treal(b2gd_it)=b2gd_tarr(b2gd_it)
1185  r=b2gd_xarr(1,b2gd_it)
1186  rh=b2gd_xarr(2,b2gd_it)
1187  t=b2gd_tarr(b2gd_it)
1188  b2gd_xreal(b2gd_it)=r
1189  end do
1190 
1191  r_final=b2gd_xarr(1,b2gd_nt)
1192  rh_final=b2gd_xarr(2,b2gd_nt)
1193  cos_bit_final=b2gd_xarr(3,b2gd_nt)
1194  sin_bit_final=b2gd_xarr(4,b2gd_nt)
1195 
1196  end subroutine get_results
1197 
1198 
1199  !--------------------------------------------------------------------------
1200  ! fcn
1209  !--------------------------------------------------------------------------
1210 
1211  subroutine fcn(t,vars,ff)
1212 
1213  use s2_types_mod
1214  use bianchi2_globaldata_mod
1215 
1216  implicit none
1217  integer, parameter :: neq = 4
1218  integer, parameter :: nvars = 4
1219 
1220  real(s2_dp), intent(in) :: t
1221  real(s2_dp), intent(in) :: vars(nvars)
1222  real(s2_dp), intent(out) :: ff(nvars)
1223 
1224  real(s2_dp) ::tau
1225  real(s2_dp) :: r,rh,lambda
1226  real(s2_dp) :: arg,fact
1227 
1228  r=vars(1)
1229  rh=vars(2)
1230  tau=t
1231 
1232  lambda=3*b2gd_rh_start**2*b2gd_omega_lambda
1233 
1234  ff(1)=r*rh
1235 
1236  ff(2)=-rh**2.d0/2.d0+lambda*r**2.d0/2.d0+b2gd_alpha**2/2.d0
1237 
1238  arg=(b2gd_alpha*tau+dlog(cos(0.3141592653589793d1/4.d0+ &
1239  b2gd_theta_0/2.d0)**2.d0+exp(-2.d0*b2gd_alpha*tau)-&
1240  exp(-2.d0*b2gd_alpha*tau)*cos(0.3141592653589793d1/4.d0+ &
1241  b2gd_theta_0/2.d0)**2.d0))/b2gd_alpha
1242 
1243  fact=-2.d0/r**2*exp(b2gd_alpha*tau)*cos(b2gd_theta_0)/(-1.d0-sin(b2gd_theta_0)+ &
1244  exp(2.d0*b2gd_alpha*tau)*sin(b2gd_theta_0)-exp(2.d0*b2gd_alpha*tau))**2.d0* &
1245  (1.d0+sin(b2gd_theta_0)+exp(2.d0*b2gd_alpha*tau)*sin(b2gd_theta_0)- &
1246  exp(2.d0*b2gd_alpha*tau))
1247 
1248  ff(3)=fact*cos(arg)
1249 
1250  ff(4)=fact*sin(arg)
1251 
1252  end subroutine fcn
1253 
1254 
1255  !--------------------------------------------------------------------------
1256  ! out
1257  !
1265  !--------------------------------------------------------------------------
1266 
1267  subroutine out(t,vals)
1268 
1269  use s2_types_mod
1270  use bianchi2_globaldata_mod
1271 
1272  implicit none
1273 
1274  real(s2_dp), intent(inout) :: t
1275  real(s2_dp), intent(in) :: vals(b2gd_nvars)
1276 
1277  real(s2_dp) :: deriv(b2gd_nvars)
1278  integer :: iel, ival
1279 
1280  iel=b2gd_nt-(b2gd_it+1)
1281  b2gd_tarr(iel)=t
1282  do ival=1,b2gd_nvars
1283  b2gd_xarr(ival,iel)=vals(ival)
1284  end do
1285 
1286  call fcn(t,vals,deriv)
1287 
1288  t=b2gd_tstop-b2gd_it*b2gd_deltat
1289  b2gd_it=b2gd_it-1
1290 
1291  end subroutine out
1292 
1293 
1294  !--------------------------------------------------------------------------
1295  ! get_tau
1296  !
1303  !--------------------------------------------------------------------------
1304 
1305  subroutine get_tau(tau_needed)
1306 
1307  use s2_types_mod
1308  use bianchi2_globaldata_mod
1309 
1310  implicit none
1311 
1312  real(s2_dp), intent(out) :: tau_needed
1313 
1314  real(s2_dp) :: aa, abserr, bb, epsabs, epsrel, result_val
1315 
1316 #ifdef NAGI8
1317  integer(kind=8) :: ifail
1318  integer(kind=8), parameter :: liwork = 10000
1319  integer(kind=8), parameter :: lwork = 10000
1320  integer(kind=8) :: iwork(liwork)
1321 #else
1322  integer :: ifail
1323  integer, parameter :: liwork = 10000
1324  integer, parameter :: lwork = 10000
1325  integer :: iwork(liwork)
1326 #endif
1327 
1328  real(s2_dp) :: work(lwork)
1329 
1330  ifail=0
1331  epsabs=1d-6
1332  epsrel=1d-8
1333 
1334  aa=1d0/sqrt(1+b2gd_ze)
1335  bb=1d0
1336  call d01ajf(f1r,aa,bb,epsabs,epsrel,result_val,abserr,work,lwork,&
1337  iwork,liwork,ifail)
1338 
1339  tau_needed=result_val
1340 
1341  end subroutine get_tau
1342 
1343 
1344  !--------------------------------------------------------------------------
1345  ! F1r
1346  !
1352  !--------------------------------------------------------------------------
1353 
1354  function f1r(a)
1355 
1356  use s2_types_mod
1357  use bianchi2_globaldata_mod
1358 
1359  implicit none
1360 
1361  real(s2_dp) :: f1r
1362  real(s2_dp), intent(in) :: a
1363 
1364  real(s2_dp) :: t0
1365 
1366  t0 = -b2gd_alpha**2*(b2gd_omega_lambda*a**6-a**2*b2gd_omega_matter-a**2* &
1367  b2gd_omega_lambda+a**2+b2gd_omega_matter)/(b2gd_omega_matter+b2gd_omega_lambda-1)
1368 
1369  f1r=2d0/sqrt(t0)
1370 
1371  end function f1r
1372 
1373 
1374  !--------------------------------------------------------------------------
1375  ! Bianchi integrals
1376  !--------------------------------------------------------------------------
1377 
1378  !--------------------------------------------------------------------------
1379  ! bianchi2_sky_comp_IX
1380  !
1391  !--------------------------------------------------------------------------
1392 
1393  function bianchi2_sky_comp_ix(l,lmax,Ntheta,X_grid,lut) result(IX)
1394 
1395  integer, intent(in) :: l,lmax,ntheta
1396  real(s2_dp), intent(in) :: x_grid(0:) ! X = A or B.
1397  type(bianchi2_lut), intent(in), optional :: lut
1398  real(s2_dp) :: ix
1399 
1400  real(s2_dp) :: integrand, plm1
1401  integer :: itheta
1402  real(s2_dp) :: dtheta, theta
1403 
1404 
1405  ix = 0d0
1406  dtheta = (pi - 0d0) / Real(ntheta, s2_dp)
1407 
1408  ! Compute the sum of all the terms.
1409  if (present(lut)) then
1410 
1411  do itheta = 0, ntheta-1
1412  theta = itheta*dtheta
1413  plm1 = bianchi2_lut_access(lut,l,itheta) ! Take the value in the table.
1414  integrand = sin(theta) * x_grid(itheta) * plm1
1415  ix = ix + integrand*dtheta
1416  end do
1417 
1418  else
1419 
1420  do itheta = 0, ntheta-1
1421  theta = itheta*dtheta
1422  plm1 = bianchi2_lut_plgndr(l, 1, cos(theta)) ! Re-compute the value.
1423  integrand = sin(theta) * x_grid(itheta) * plm1
1424  ix = ix + integrand*dtheta
1425  end do
1426 
1427  end if
1428 
1429  ! Compute the definitive integral.
1430  ix = ix * sqrt( (2d0*l+1d0) / real(4d0*PI*l*(l+1d0), s2_dp) )
1431  return
1432 
1433  end function bianchi2_sky_comp_ix
1434 
1435 
1436 end module bianchi2_sky_mod
Generated on Tue Jul 28 2015 20:36:53 by   doxygen 1.8.3.1