1!----------------------------------------------------------------------------- 2! 3! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon, 4! Benjamin D. Wandelt, Anthony J. Banday, 5! Matthias Bartelmann, Hans K. Eriksen, 6! Frode K. Hansen, Martin Reinecke 7! 8! 9! This file is part of HEALPix. 10! 11! HEALPix is free software; you can redistribute it and/or modify 12! it under the terms of the GNU General Public License as published by 13! the Free Software Foundation; either version 2 of the License, or 14! (at your option) any later version. 15! 16! HEALPix is distributed in the hope that it will be useful, 17! but WITHOUT ANY WARRANTY; without even the implied warranty of 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19! GNU General Public License for more details. 20! 21! You should have received a copy of the GNU General Public License 22! along with HEALPix; if not, write to the Free Software 23! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 24! 25! For more information about HEALPix see http://healpix.sourceforge.net 26! 27!----------------------------------------------------------------------------- 28module alm_tools 29 ! Scalar+Open_MP implementation 30 ! 31 ! function do_opemp 32 ! subroutine init_rescale 33 ! subroutine get_pixel_layout 34 ! subroutine select_rings 35 ! subroutine gen_recfac 36 ! subroutine gen_recfac_spin 37 ! subroutine gen_lamfac 38 ! subroutine gen_lamfac_der 39 ! subroutine gen_mfac 40 ! subroutine gen_mfac_spin 41 ! subroutine compute_lam_mm 42 ! subroutine do_lam_lm 43 ! subroutine do_lam_lm_spin 44 ! subroutine do_lam_lm_pol 45 ! subroutine gen_normpol 46 ! function l_min_ylm 47 ! subroutine warning_oldbounds 48 ! 49 ! subroutine ring_synthesis 50 ! subroutine ring_analysis 51 ! 52 ! -------------------- in include files (see alm_map_template.f90) --------- 53 ! subroutine alm2map_sc 54 ! subroutine alm2map_spin 55 ! subroutine alm2map_sc_pre 56 ! subroutine alm2map_pol 57 ! subroutine alm2map_pol_pre1 58 ! subroutine alm2map_pol_pre2 59 ! subroutine alm2map_sc_der 60 ! subroutine alm2map_pol_der 61 ! 62 ! subroutine map2alm_sc 63 ! subroutine map2alm_spin 64 ! subroutine map2alm_sc_pre 65 ! subroutine map2alm_pol 66 ! subroutine map2alm_pol_pre1 67 ! subroutine map2alm_pol_pre2 68 ! 69 ! subroutine map2alm_iterative 70 ! 71 ! subroutine alter_alm 72 ! subroutine create_alm 73 ! subroutine alm2cl 74 ! 75 ! subroutine rotate_alm 76 ! ------------------------------------------------------------------ 77 ! 78 ! subroutine plm_gen 79 ! 80 ! subroutine pow2alm_units 81 ! subroutine generate_beam 82 ! subroutine gaussbeam 83 ! subroutine pixel_window 84 ! 85 use misc_utils, only: assert_alloc, assert_present, assert, fatal_error, string, strupcase !, wall_clock_time 86 use healpix_types 87 use healpix_fft, only: real_fft2, planck_fft2_plan, make_fft2_plan, destroy_fft2_plan, & 88 & fft2_forward, fft2_backward 89 IMPLICIT none 90 91 ! keep everything private unless stated otherwise 92 private 93 !-------------------------------------------------------------------------- 94 ! define large and small numbers used to renormalise the recursion on the Legendre Polynomials 95 integer(I4B), private, parameter :: LOG2LG = 100 96 real(KIND=DP), private, parameter :: FL_LARGE = 2.0_dp ** LOG2LG 97 real(KIND=DP), private, parameter :: FL_SMALL = 2.0_dp ** (-LOG2LG) 98 ! declare array for dynamic rescaling of the Ylm 99 integer(kind=i4b), private, parameter :: RSMAX = 20, RSMIN = -20 100 real(dp), private, dimension(RSMIN:RSMAX) :: rescale_tab 101 real(DP), private, parameter :: ALN2_INV = 1.4426950408889634073599246810_dp ! 1/log(2) 102 ! misc 103 integer(i4b), private, parameter :: SMAXCHK = 50 ! maximum size of chunk (in number of ring pairs) 104 ! parameters of Ylm short-cut 105 integer(kind=i4b), private, parameter :: HPX_MXL0 = 40 ! minimum used, choose <=0 to do full loop 106 real (kind=dp), private, parameter :: HPX_MXL1 = 1.35_dp 107 !-------------------------------------------------------------------------- 108 109 ! make (front end) routines public 110 public :: alm2map, map2alm, alm2map_der, alm2map_spin, map2alm_spin 111 public :: map2alm_iterative 112 public :: map2alm_iterative_old 113 public :: alter_alm, create_alm, alm2cl, rotate_alm 114 public :: create_alm_old ! for tests only 115 public :: plm_gen 116 public :: ring_synthesis, ring_analysis 117 public :: generate_beam, gaussbeam, pixel_window, pow2alm_units 118 119 interface alm2cl 120 module procedure alm2cl2_s, alm2cl2_d, alm2cl1_s, alm2cl1_d 121 end interface 122 123 interface sub_alm2cl 124 module procedure sub_alm2cl_s, sub_alm2cl_d 125 end interface 126 127 interface rotate_alm 128 module procedure rotate_alm_s, rotate_alm_d 129 end interface 130 131 interface alter_alm 132 module procedure alter_alm_s, alter_alm_d 133 end interface 134 135 interface create_alm 136 module procedure create_alm_s, create_alm_d, create_alm_v12_s, create_alm_v12_d 137 end interface 138 139 interface create_alm_old ! for tests only 140 module procedure create_alm_old_s, create_alm_old_d 141 end interface 142 143 interface alm2map_der 144 module procedure alm2map_sc_der_s, alm2map_sc_der_d, alm2map_pol_der_s, alm2map_pol_der_d 145 end interface 146 147 interface alm2map_spin 148 module procedure alm2map_spin_s, alm2map_spin_d 149 end interface 150 151 interface map2alm_spin 152 module procedure map2alm_spin_s, map2alm_spin_d 153 end interface 154 155 interface alm2map 156 ! Call the correct alm2map routine, depending on whether polarisation is included 157 ! (determined from the rank of the map_TQU array) or precomputed plms 158 ! (scalar or tensor, determined from the rank of the plm array) are included, 159 ! or whether the map and alm are single or double precision 160 module procedure alm2map_sc_wrapper_s, & 161 & alm2map_pol_wrapper_s, alm2map_pol_pre2_s, & 162 & alm2map_sc_wrapper_d, & 163 & alm2map_pol_wrapper_d, alm2map_pol_pre2_d 164!!! module procedure alm2map_sc_s, alm2map_sc_pre_s, & 165!!! & alm2map_pol_s, alm2map_pol_pre1_s, alm2map_pol_pre2_s, & 166!!! & alm2map_sc_d, alm2map_sc_pre_d, & 167!!! & alm2map_pol_d, alm2map_pol_pre1_d, alm2map_pol_pre2_d 168 end interface 169 170 interface map2alm 171 ! Call the correct map2alm routine, depending on whether polarisation is included 172 ! (determined from the rank of the map_TQU array) or precomputed plms 173 ! (scalar or tensor, determined from the rank of the plm array) are included, 174 ! or whether the map and alm are single or double precision 175 module procedure map2alm_old_sc_s, map2alm_old_pol_s, map2alm_old_pol2_s, & 176 & map2alm_sc_s, map2alm_sc_pre_s, & 177 & map2alm_pol_s, map2alm_pol_pre1_s, map2alm_pol_pre2_s, & 178 & map2alm_old_sc_d, map2alm_old_pol_d, map2alm_old_pol2_d, & 179 & map2alm_sc_d, map2alm_sc_pre_d, & 180 & map2alm_pol_d, map2alm_pol_pre1_d, map2alm_pol_pre2_d 181 end interface 182 183 interface map2alm_iterative 184 module procedure map2alm_iterative_s, map2alm_iterative_d 185 end interface 186 interface map2alm_iterative_old 187 module procedure map2alm_iterative_old_s, map2alm_iterative_old_d 188 end interface 189 190 interface sub_map2ring 191 module procedure sub_map2ring_1d_s, sub_map2ring_1d_d, sub_map2ring_nd_s, sub_map2ring_nd_d 192 end interface 193 interface sub_ring2map 194 module procedure sub_ring2map_1d_s, sub_ring2map_1d_d, sub_ring2map_nd_s, sub_ring2map_nd_d 195 end interface 196 ! make routines public as most of them are called by mpi_alm* routines 197 public :: do_openmp 198 public :: init_rescale 199 public :: do_lam_lm, do_lam_lm_pol 200 ! needed by mpi_alm* 201 public :: l_min_ylm 202 public :: get_pixel_layout 203 public :: gen_recfac, gen_lamfac, gen_lamfac_der, gen_mfac, compute_lam_mm, gen_normpol 204 public :: select_rings ! added Feb 2006 205 ! 206 public :: gen_recfac_spin, gen_mfac_spin, do_lam_lm_spin ! July 2007 207 208 !========================================================= 209 210 ! 211 ! Aug 14, 2000, EH, Caltech, switch to CMBFAST normalisation of polarisation power spectra 212 ! (see Zaldarriaga, astro-ph/9709271; ApJ, 503, 1 (1998)) 213 ! the factor normal_l (used for synthesis and analysis) is reduced by sqrt(2) 214 ! 215 ! Oct 17, 2001 216 ! polarised alm expression is larger by a factor 2 217 ! introduced FL_LARGE and FL_SMALL parameters 218 ! 219 ! Sept 2001, EH, Caltech: map2alm* : skip calculations for cut out pixels 220 ! : use explicit loops for multiplication of alm by omega_pix 221 ! 222 ! Nov 2001, EH, Caltech : added KIND=DP in all CMPLX commands that missed it 223 ! turned 1.d0 in 1.0_dp and so on 224 ! moved 'use healpix_types' and 'implicit none' to module top level 225 ! 226 ! Dec 2001, EH, Caltech : corrected error on inner loop upper bound of lam_fact initialisation 227 ! in alm2map_pol_pre1, alm2map_pol, map2alm_pol_pre1, map2alm_pol 228 ! pointed out by M. Ashdown to Level S 229 ! 230 ! added pow2alm_units, gaussbeam, generate_beam and alter_alm (moved for smoothing/alter_alm.f90) 231 ! 232 ! 233 ! Nov 2002, EH, reordered routines to simplify maintenance of scalar vs parallel routines 234 ! ------post 1.21 235 ! Dec 2003-Jan 2004, EH, remove 'trig' array from ring_analysis and ring_synthesis 236 ! 2004-05-28, EH, edition of pixel_window 237 ! Aug 2004 : put ring in DP in ring_synthesis 238 ! add alm2map_sc_der for calculation of derivatives 239 ! moved plm_gen here from plmgen main code 240 ! Oct 2004: reduced size of recfac and lam_fact 241 ! : introduced m_max_syn and m_max_ana 242 ! Nov 2004: streamlining of *alm* routines to gain CPU time 243 ! Dec 2004: editions to alter_alm (new argument), pixel_window 244 ! April 2005: separate Lambda_lm calculation (with do_lam_lm[_pol]) from scalar 245 ! product in alm2map_sc and alm2map_pol for speed up. 246 ! does not improve the speed of map2alm though, so stick to former code for those 247 ! May 2005: pixel window returns 1. if nside = 0 (interpreted as infinitely small pixel) 248 ! Aug 2005: added alm2map_pol_der 249 ! ------post 2.00 250 ! Sep-Oct 2005: made internal subroutines public so they can be called by mpi_alm, 251 ! corrected OpenMP+SGI bug in alm2map_pol_der 252 ! Feb 2006: added select_rings to public routines (for MPI use) 253 ! Sep 2006: corrected resetting of lam_fact in gen_lamfac 254 ! July, Sep-Oct 2007: added alm2map_spin, map2alm_spin 255 ! Oct 2007: zbounds and w8rings optional in map2alm (without plm) 256 ! Nov 2008: correction of a crash bug in map2alm_iterative 257 ! Jan 2009: correction of a PGF90 specific bug in create_alm 258 ! Jan 2010: correction of a bug in alm2map_pol_der affecting Q and U derivatives 259 ! Jan 2011: use get_healpix_pixel_window_file in pixel_window() 260 ! ===================================================== 261 ! about the F90 compilers 'features' 262 ! - Intel compiler (ifc) (and maybe the other compilers as well) 263 ! is apparently very slow at allocating arrays 264 ! therefore as much as possible the 'allocate' should be done outside the loops. 265 ! This what is done below, expect when OpenMP require thread safe arrays that much 266 ! be created each time we enter the loop. 267 ! 268 ! - NagF95 (f95) is very slow at creating vectors with the operator (/ ... /) 269 ! therefore, in a loop, the operation 270 ! x(0:4) = x(0:4) + (/ a, b, c, d /) should be replaced by the equivalent 271 ! but MUCH faster, explicit form 272 ! x(0) = x(0) + a 273 ! x(1) = x(1) + b 274 ! ... 275 ! unless (/ a,b,c,d /) can be constructed before entering the loop 276 ! 277contains 278 279 !************************************************************************** 280 ! 281 ! ALM2MAP/MAP2ALM SIDEKICKS 282 ! 283 !************************************************************************** 284 !================================================ 285 function do_openmp() 286 !================================================ 287 ! returns .true. if code was compiled with OpenMP 288 !================================================ 289 logical(LGT) :: do_openmp 290 !------------------------------------------------ 291 292 do_openmp = .false. 293! DO NOT REMOVE FOLLOWING LINES 294!$ do_openmp = .true. ! Intel f90 295!IBMP do_openmp = .true. ! IBM xlf90 296! ----------------------------- 297 298 return 299 end function do_openmp 300 301 !================================================ 302 subroutine init_rescale() 303 !================================================ 304 ! local variables 305 integer(i4b) :: s, smax 306 real(dp) :: logOVFLOW 307 character(len=*), parameter :: code = 'gen_rescale' 308 !------------------------------------------------ 309 logOVFLOW=log(FL_LARGE) 310 smax = INT( log(MAX_DP) / logOVFLOW ) 311 312 if (smax > (RSMAX-1)) then 313 print*,'Array rescale_tab too small in '//code 314 print*,smax ,'>', RSMAX 315 stop 316 endif 317 318 rescale_tab(RSMIN:RSMAX) = 0.0_dp 319 do s = -smax, smax 320 rescale_tab(s) = FL_LARGE ** s 321 enddo 322 rescale_tab(0) = 1.0_dp 323 324 return 325 end subroutine init_rescale 326 327 !======================================================================= 328 subroutine get_pixel_layout(nside, ith, cth, sth, nphi, startpix, kphi0) 329 !======================================================================= 330 ! output Healpix pixel layout for the ring ith in [0,2*nside] 331 !======================================================================= 332 integer(I4B), intent(IN) :: nside, ith 333 real(DP) , intent(OUT) :: cth, sth 334 integer(I4B), intent(OUT) :: nphi, kphi0 335 integer(I8B), intent(OUT) :: startpix 336 ! 337 integer(I4B) :: nrings 338 real(DP) :: dth1, dth2, dst1 339 !======================================================================= 340 341 nrings = 2*nside 342 if (ith < 1 .or. ith> nrings) then 343 print*,'ith out of bounds ',ith,1,nrings 344 call fatal_error 345 endif 346 347 dth1 = 1.0_dp / (3.0_dp*DBLE(nside)**2) 348 dth2 = 2.0_dp / (3.0_dp*DBLE(nside)) 349 dst1 = 1.0_dp / (SQRT(6.0_dp) * DBLE(nside) ) 350 351 if (ith < nside) then ! polar cap (north) 352 cth = 1.0_dp - DBLE(ith)**2 * dth1 353 nphi = 4*ith 354 kphi0 = 1 355 sth = SIN( 2.0_dp * ASIN( ith * dst1 ) ) ! sin(theta) 356 startpix = 2*ith*(ith-1_I8B) 357 else ! tropical band (north) + equator 358 cth = DBLE(2*nside-ith) * dth2 359 nphi = 4*nside 360 kphi0 = MOD(ith+1-nside,2) 361 sth = DSQRT((1.0_dp-cth)*(1.0_dp+cth)) ! sin(theta) 362 startpix = 2*nside*(nside-1_I8B) + (ith-nside)*int(nphi,kind=I8B) 363 endif 364 365 return 366 end subroutine get_pixel_layout 367 !======================================================================= 368 subroutine select_rings_old(z, zbounds, keep_north, keep_south, keep_either) 369 !======================================================================= 370 ! select rings laying within zbounds 371 ! if zbounds(1) < zbounds(2) : keep zbounds(1) < z < zbounds(2) 372 ! if zbounds(2) < zbounds(1) : keep z < zbounds(2) Union zbounds(1) < z 373 ! if zbounds(1)=zbounds(2) : keep everything 374 ! input z should be >= 0 375 !======================================================================= 376 real(DP) , intent(in) :: z 377 real(DP) , intent(in), dimension(1:2) :: zbounds 378 logical(LGT), intent(OUT) :: keep_north, keep_south, keep_either 379 ! 380 real(DP) :: zn, zs 381 !======================================================================= 382 383 ! if (zbounds(1) = zbounds(2)) keep everything 384 if (abs(zbounds(1)-zbounds(2)) < 1.e-6) then 385 keep_north = .true. 386 keep_south = .true. 387 keep_either = .true. 388 return 389 endif 390 391 zn = abs(z) 392 zs = -zn 393 394 if (zbounds(1) < zbounds(2)) then 395 ! inner ring 396 keep_north = (zn >= zbounds(1) .and. zn <= zbounds(2)) 397 keep_south = (zs >= zbounds(1) .and. zs <= zbounds(2)) 398 399 else 400 ! outter ring 401 keep_north = (zn > zbounds(1) .or. zn < zbounds(2)) 402 keep_south = (zs > zbounds(1) .or. zs < zbounds(2)) 403 endif 404 keep_either = keep_north .or. keep_south 405 406 407 return 408 end subroutine select_rings_old 409 410 !======================================================================= 411 subroutine select_rings(z, zbounds, keep_north, keep_south, keep_either) 412 !======================================================================= 413 ! select rings laying within zbounds 414 ! if zbounds(1) < zbounds(2) : keep zbounds(1) < z < zbounds(2) 415 ! if zbounds(2) <= zbounds(1) : keep z < zbounds(2) Union zbounds(1) < z 416 ! input z should be >= 0 417 ! edited 2018-11-09 to be identical to remove_dipole, apply_mask 418! and libsharp's hpsharp_make_healpix_geom_info_2 419 !======================================================================= 420 real(DP) , intent(in) :: z 421 real(DP) , intent(in), dimension(1:2) :: zbounds 422 logical(LGT), intent(OUT) :: keep_north, keep_south, keep_either 423 ! 424 real(DP) :: zn, zs 425 !======================================================================= 426 427 428 zn = abs(z) 429 zs = -zn 430 431 if (zbounds(1) < zbounds(2)) then 432 ! inner ring 433 keep_north = (zn > zbounds(1) .and. zn < zbounds(2)) 434 keep_south = (zs > zbounds(1) .and. zs < zbounds(2)) 435 436 else 437 ! outter ring 438 keep_north = (zn > zbounds(1) .or. zn < zbounds(2)) 439 keep_south = (zs > zbounds(1) .or. zs < zbounds(2)) 440 endif 441 keep_either = keep_north .or. keep_south 442 443 444 return 445 end subroutine select_rings 446 447 !======================================================================= 448 subroutine gen_recfac( l_max, m, recfac) 449 !======================================================================= 450 ! generates recursion factors used to computes the Ylm of degree m 451 ! for all l in m<=l<=l_max 452 !======================================================================= 453 integer(I4B), intent(IN) :: l_max, m 454 real(DP), intent(OUT), dimension(0:1, 0:l_max) :: recfac 455 ! 456 real(DP) :: fm2, fl2 457 integer(I4B) :: l 458 459 recfac(0:1,0:m) = 0.0_dp 460 fm2 = DBLE(m) **2 461 do l = m, l_max 462 fl2 = DBLE(l+1) **2 463 recfac(0,l) = DSQRT( (4.0_dp * fl2 - 1.0_dp) / (fl2-fm2) ) 464 enddo 465 ! put outside the loop because of problem on some compilers 466 recfac(1,m:l_max) = 1.0_dp / recfac(0,m:l_max) 467 468 469 return 470 end subroutine gen_recfac 471 472 !======================================================================= 473 subroutine gen_recfac_spin( l_max, m, spin, recfac_spin) 474 !======================================================================= 475 ! generates recursion factors used to computes the Ylms of degree m and spin s 476 ! for all l in max(m,s)<=l<=l_max 477 !======================================================================= 478 integer(I4B), intent(IN) :: l_max, m, spin 479 real(DP), intent(OUT), dimension(0:2, 0:l_max) :: recfac_spin 480 ! 481 real(DP) :: fm2, fl2, fs2 482 integer(I4B) :: l, l_low, aspin 483 484 485 aspin = abs(spin) 486 l_low = max(m, aspin) 487 488 recfac_spin(0:2, 0:l_max) = 0.0_dp 489 fm2 = DBLE(m) **2 490 fs2 = DBLE(spin) **2 491 do l = l_low, l_max 492 fl2 = DBLE(l+1) **2 493 recfac_spin(0,l) = DSQRT( (4.0_dp * fl2 - 1.0_dp) / (fl2-fm2) /(1.0_dp-fs2/fl2) ) 494 enddo 495 do l = max(l_low, 1), l_max 496 recfac_spin(2,l) = aspin * m / dble( l * (l+1) ) 497 enddo 498 ! put outside the loop because of problem on some compilers 499 recfac_spin(1,l_low:l_max) = 1.0_dp / recfac_spin(0,l_low:l_max) 500 501 502 return 503 end subroutine gen_recfac_spin 504 505 !======================================================================= 506 subroutine gen_lamfac( l_max, m, lam_fact) 507 !======================================================================= 508 ! generates factor relating scalar Ylm to polar Ylm 509 ! for all l in m<=l<=l_max 510 !======================================================================= 511 integer(I4B), intent(IN) :: l_max, m 512 real(DP), intent(OUT), dimension(0:l_max) :: lam_fact 513 ! 514 real(DP) :: fm2, fl, fl2 515 integer(I4B) :: l 516 517! lam_fact(0:m) = 0. 518 lam_fact(0:max(1,m)) = 0.0_dp 519 fm2 = real(m * m, kind=DP) 520 do l = max(2,m+1), l_max 521 fl = real(l, kind=dp) 522 fl2 = fl * fl 523 lam_fact(l) = 2.0_dp * SQRT( (2.0_dp * fl + 1.0_dp) / (2.0_dp * fl - 1.0_dp) * (fl2-fm2) ) 524 enddo 525 526 return 527 end subroutine gen_lamfac 528 529 !======================================================================= 530 subroutine gen_lamfac_der(l_max, m, lam_fact) 531 !======================================================================= 532 ! generates factor relating scalar Ylm to its derivatives 533 ! for all l in m<=l<=l_max 534 !======================================================================= 535 integer(I4B), intent(IN) :: l_max, m 536 real(DP), intent(OUT), dimension(0:l_max) :: lam_fact 537 ! 538 real(DP) :: fm2, fl, fl2 539 integer(I4B) :: l 540 541 lam_fact(0:m) = 0. 542 fm2 = real(m * m, kind=DP) 543 do l = max(1,m+1), l_max ! different lower bound than pol. factor 544 fl = real(l, kind=dp) 545 fl2 = fl * fl 546 lam_fact(l) = SQRT( (2.0_dp * fl + 1.0_dp) / (2.0_dp * fl - 1.0_dp) * (fl2-fm2) ) 547 ! different normalization than polarization factor 548 enddo 549 550 return 551 end subroutine gen_lamfac_der 552 553 !======================================================================= 554 subroutine gen_mfac( m_max, m_fact) 555 !======================================================================= 556 ! generates factor used in lam_mm calculation 557 ! for all m in 0<=m<=m_max 558 !======================================================================= 559 integer(I4B), intent(IN) :: m_max 560 real(DP), intent(OUT), dimension(0:m_max) :: m_fact 561 ! 562 integer(I4B) :: m 563 564 ! fact(m) = fact(m-1) * sqrt( (2m+1) / (2m) ) 565 m_fact(0) = 1.0_dp 566 do m=1,m_max 567 m_fact(m) = m_fact(m-1)*sqrt(dble(2*m+1)/dble(2*m)) 568 end do 569 570 ! Log_2 ( fact(m) / sqrt(4 Pi) ) 571 do m=0,m_max 572 m_fact(m) = log(SQ4PI_INV * m_fact(m)) * ALN2_INV 573 enddo 574 575 return 576 end subroutine gen_mfac 577 !======================================================================= 578 subroutine gen_mfac_spin( m_max, spin, m_fact) 579 !======================================================================= 580 ! generates factor used in lam_mm,s calculation 581 ! for all m in 0<=m<=m_max 582 !======================================================================= 583 integer(I4B), intent(IN) :: m_max, spin 584 real(DP), intent(OUT), dimension(0:m_max) :: m_fact 585 ! 586 integer(I4B) :: m, aspin 587 real(DP) :: tmp 588 589 aspin = abs(spin) 590 m_fact(:) = -1.e30 591 if (aspin <= m_max) m_fact(aspin) = 1.0_dp 592 ! fact(m) = fact(m+1) * sqrt((s+m+1)/(s-m) ) 593 if (aspin > 0) then 594 tmp = 1.d0 595 do m = aspin-1, 0, -1 596 tmp = tmp * sqrt( dble(aspin+m+1)/dble(aspin-m) ) 597 if (m <= m_max) m_fact(m) = tmp 598 enddo 599 endif 600 ! fact(m) = fact(m-1) * sqrt(m*(2m+1)/(m-s)/(m+s)/2 ) 601 do m = aspin+1, m_max 602 m_fact(m) = m_fact(m-1)*sqrt( m*dble(2*m+1)/dble(2*(m-aspin)*(m+aspin)) ) 603 end do 604 605 ! Log_2 ( fact(m) / sqrt(4 Pi)) 606 do m=0,m_max 607 m_fact(m) = log(SQ4PI_INV * m_fact(m)) * ALN2_INV 608 enddo 609 610 return 611 end subroutine gen_mfac_spin 612 !======================================================================= 613 subroutine compute_lam_mm(mfac, m, sth, lam_mm, corfac, scalem) 614 !======================================================================= 615 ! computes lam_mm 616 ! the true lam_mm is lam_mm * corfac 617 !======================================================================= 618 integer(I4B), intent(in) :: m 619 real(DP), intent(in) :: sth, mfac 620 real(DP), intent(out) :: lam_mm, corfac 621 integer(I4B), intent(out) :: scalem 622 ! 623 real(DP) :: log2val, dlog2lg 624 625 dlog2lg = real(LOG2LG, kind=DP) 626 627 log2val = mfac + m*log(sth) * ALN2_INV ! log_2(lam_mm) 628 scalem = int (log2val / dlog2lg) 629 corfac = rescale_tab(max(scalem,RSMIN)) 630 lam_mm = 2.0_dp **(log2val - scalem * dlog2lg) ! rescaled lam_mm 631 if (IAND(m,1)>0) lam_mm = -lam_mm ! negative for odd m 632 633 return 634 end subroutine compute_lam_mm 635 636 !======================================================================= 637 subroutine do_lam_lm(lmax, m, cth, sth, mfac, recfac, lam_lm) 638 !======================================================================= 639 ! computes scalar lambda_lm(theta) for all l in [m,lmax] for a given m, and given theta 640 ! input: lmax, m, cos(theta), sin(theta) 641 ! mfac: precomputed (by gen_mfac) quantity useful for lambda_mm calculation 642 ! recfac: precomputed (by gen_recfac) quantities useful for 643 ! lambda_lm recursion for a given m 644 ! output: lam_lm 645 ! the routine also needs the array rescale_tac initialized by init_rescale 646 !======================================================================= 647 integer(I4B), intent(in) :: lmax, m 648 real(DP), intent(in) :: cth, sth, mfac 649 real(DP), dimension(0:1,0:lmax), intent(in) :: recfac 650 real(DP), dimension( 0:lmax), intent(out) :: lam_lm 651 ! 652 real(DP) :: log2val, dlog2lg 653 real(DP) :: ovflow, unflow, corfac 654 real(DP) :: lam_mm, lam_0, lam_1, lam_2 655 integer(I4B) :: scalel, l, l_min 656 !--------------------------------------------------------------- 657 658 ! define constants 659 ovflow = rescale_tab(1) 660 unflow = rescale_tab(-1) 661 l_min = l_min_ylm(m, sth) 662 dlog2lg = real(LOG2LG, kind=DP) 663 664 ! computes lamba_mm 665 log2val = mfac + m*log(sth) * ALN2_INV ! log_2(lam_mm) 666 scalel = int (log2val / dlog2lg) 667 corfac = rescale_tab(max(scalel,RSMIN)) 668 lam_mm = 2.0_dp **(log2val - scalel * dlog2lg) ! rescaled lam_mm 669 if (IAND(m,1)>0) lam_mm = -lam_mm ! negative for odd m 670 671 lam_lm(0:lmax) = 0.0_dp 672 ! --- l = m --- 673 lam_lm(m) = lam_mm * corfac 674 675 ! --- l > m --- 676 lam_0 = 0.0_dp 677 lam_1 = 1.0_dp 678 lam_2 = cth * lam_1 * recfac(0,m) 679 do l = m+1, lmax 680 ! do recursion 681 if (l >= l_min) then 682 lam_lm(l) = lam_2 * corfac * lam_mm 683 endif 684 lam_0 = lam_1 * recfac(1,l-1) 685 lam_1 = lam_2 686 lam_2 = (cth * lam_1 - lam_0) * recfac(0,l) 687 688 ! do dynamic rescaling 689 if (abs(lam_2) > ovflow) then 690 lam_1 = lam_1*unflow 691 lam_2 = lam_2*unflow 692 scalel= scalel + 1 693 corfac = rescale_tab(max(scalel,RSMIN)) 694 elseif (abs(lam_2) < unflow .and. abs(lam_2) /= 0.0_dp) then 695 lam_1 = lam_1*ovflow 696 lam_2 = lam_2*ovflow 697 scalel= scalel - 1 698 corfac = rescale_tab(max(scalel,RSMIN)) 699 endif 700 701 enddo ! loop on l 702 end subroutine do_lam_lm 703 !======================================================================= 704 subroutine do_lam_lm_spin(lmax, m, spin, cth, sth, mfac, mfac_spin, recfac_spin, lam_lm_spin) 705 !======================================================================= 706 ! computes spin lambda_lm(theta) for all l in [m,lmax] for a given m, spin, and theta 707 ! input: lmax, m, spin, cos(theta), sin(theta) 708 ! mfac: precomputed (by gen_mfac) quantity useful for lambda_mm calculation 709 ! mfac_spin: precomputed (by gen_mfac_spin) quantity useful for lambda_mm calculation 710 ! recfac_spin: precomputed (by gen_recfac_spin) quantities useful for 711 ! lambda_lm recursion for a given m 712 ! output: lam_lm 713 ! the routine also needs the array rescale_tab initialized by init_rescale 714 !======================================================================= 715 integer(I4B), intent(in) :: lmax, m, spin 716 real(DP), intent(in) :: cth, sth, mfac, mfac_spin 717 real(DP), dimension(0:2,0:lmax), intent(in) :: recfac_spin 718 real(DP), dimension(1:2,0:lmax), intent(out) :: lam_lm_spin 719 ! 720 real(DP) :: dlog2lg 721 real(DP) :: ovflow, unflow 722 real(DP) :: lam_0, lam_1, lam_2 723 real(DP) :: thetao2, ttho2, stho2, ctho2, kss, tmp1, tmp2 724 real(DP), dimension(1:2) :: log10sss, corfac, log2val, lam_mm 725 integer(I4B) :: l, l_min, l_low, aspin, iss 726 integer(I4B), dimension(1:2) :: scalel 727 !--------------------------------------------------------------- 728 729 lam_lm_spin(1:2,0:lmax) = 0.0_dp 730 aspin = abs(spin) 731 l_low = max(m, aspin) 732 if (l_low > lmax) return 733 734 ! define constants 735 ovflow = rescale_tab(1) 736 unflow = rescale_tab(-1) 737 l_min = l_min_ylm(m, sth) 738 dlog2lg = real(LOG2LG, kind=DP) 739 740 thetao2 = atan2(sth, cth)/2.d0 741 ttho2 = tan(thetao2) 742 stho2 = sin(thetao2) ! sqrt(1-cth)/sqrt(2) 743 ctho2 = cos(thetao2) ! sqrt(1+cth)/sqrt(2) 744 ! lambda(s,s, s) = (-)^s * sqrt(2s+1/4Pi) sin(theta/2)^(2s) 745 ! lambda(s,s,-s) = (-)^s * sqrt(2s+1/4Pi) cos(theta/2)^(2s) 746 log10sss(1) = (2*aspin *log(stho2) + 0.5d0*log(2*aspin+1.0_dp)) ! log_10(abs(lam_sss)*sqrt(4Pi)) 747 log10sss(2) = (2*aspin *log(ctho2) + 0.5d0*log(2*aspin+1.0_dp)) ! log_10(abs(lam_ss-s)*sqrt(4Pi)) 748 ! 749 if (m >= aspin) then 750 ! computes directly lambda(m,m,s) 751 ! lambda(m,m,s) = - lambda(m-1,m-1,s) * sin(theta) * sqrt(m(2m+1)/2/(m^2-s^2)) 752 log2val(1:2) = mfac_spin + ((m-aspin)*log(sth) + log10sss(1:2)) * ALN2_INV ! log_2(lam_mms) 753 scalel(1:2) = int (log2val(1:2) / dlog2lg) 754 corfac(1) = rescale_tab(max(scalel(1),RSMIN)) 755 corfac(2) = rescale_tab(max(scalel(2),RSMIN)) 756 lam_mm(1:2) = 2.0_dp **(log2val(1:2) - scalel(1:2) * dlog2lg) ! rescaled lam_mm 757 if (IAND(m,1)>0) lam_mm(1:2) = -lam_mm(1:2) ! negative for odd m 758 else 759 ! computes lambda(s,m,s) 760 ! lambda(s,m,s) = - lambda(s,m+1, s) / tan(theta/2) * sqrt((s+m+1)/(s-m)) 761 ! lambda(s,m,-s) = + lambda(s,m+1,-s) * tan(theta/2) * sqrt((s+m+1)/(s-m)) 762 log2val(1) = mfac_spin + ( (m-aspin)*log(ttho2) + log10sss(1)) * ALN2_INV ! log_2(lam_mms) 763 log2val(2) = mfac_spin + (-(m-aspin)*log(ttho2) + log10sss(2)) * ALN2_INV ! log_2(lam_mms) 764 scalel(1:2) = int (log2val(1:2) / dlog2lg) 765 corfac(1) = rescale_tab(max(scalel(1),RSMIN)) 766 corfac(2) = rescale_tab(max(scalel(2),RSMIN)) 767 lam_mm(1:2) = 2.0_dp **(log2val(1:2) - scalel(1:2) * dlog2lg) ! rescaled lam_mm 768 if (IAND(m,1)>0) lam_mm(1) = -lam_mm(1) ! negative for odd m, and s>0 769 if (IAND(aspin,1)>0) lam_mm(2) = -lam_mm(2) ! negative for odd s, and s<0 770 endif 771 772 ! --- l = m --- 773 lam_lm_spin(1:2,l_low) = lam_mm(1:2) * corfac(1:2) 774 775 do iss=1,2 ! +spin and then -spin 776 kss = 1.0_dp 777 if (iss == 2) kss = -1.0_dp 778 ! --- l > m --- 779 lam_0 = 0.0_dp 780 lam_1 = 1.0_dp 781 lam_2 = (cth + kss*recfac_spin(2,l_low)) * lam_1 * recfac_spin(0,l_low) 782 do l = l_low+1, lmax 783 ! do recursion 784 if (l >= l_min) then 785 lam_lm_spin(iss,l) = lam_2 * corfac(iss) * lam_mm(iss) 786 endif 787 lam_0 = lam_1 * recfac_spin(1,l-1) 788 lam_1 = lam_2 789 lam_2 = ((cth + kss*recfac_spin(2,l)) * lam_1 - lam_0) * recfac_spin(0,l) 790 791 ! do dynamic rescaling 792 if (abs(lam_2) > ovflow) then 793 lam_1 = lam_1*unflow 794 lam_2 = lam_2*unflow 795 scalel(iss)= scalel(iss) + 1 796 corfac(iss) = rescale_tab(max(scalel(iss),RSMIN)) 797 elseif (abs(lam_2) < unflow .and. abs(lam_2) /= 0.0_dp) then 798 lam_1 = lam_1*ovflow 799 lam_2 = lam_2*ovflow 800 scalel(iss)= scalel(iss) - 1 801 corfac(iss) = rescale_tab(max(scalel(iss),RSMIN)) 802 endif 803 enddo ! loop on l 804 enddo ! loop on spin sign (iss) 805 806 ! compute functions with fixed parity 807 ! beware: the first one is always the half sum 808 ! the second one is always the half difference, independently of spin 809 do l=0, lmax 810 tmp1 = lam_lm_spin(1,l) * 0.5_dp 811 tmp2 = lam_lm_spin(2,l) * 0.5_dp 812 lam_lm_spin(1,l) = tmp1 + tmp2 813 lam_lm_spin(2,l) = tmp1 - tmp2 814 enddo 815 816 end subroutine do_lam_lm_spin 817 !======================================================================= 818 subroutine do_lam_lm_pol(lmax, m, cth, sth, mfac, recfac, lam_fact, lam_lm) 819 !======================================================================= 820 ! computes temperature&polar lambda_lm(p,theta) for all l in [m,lmax] for a given m, and given theta 821 ! input: lmax, m, cos(theta), sin(theta) 822 ! mfac: precomputed (by gen_mfac) quantity useful for lambda_mm calculation 823 ! recfac: precomputed (by gen_recfac) quantities useful for 824 ! lambda_lm recursion for a given m 825 ! lam_fact: precomputed (by gen_lamfac) factor useful for polarised lambda recursion 826 ! output: lam_lm for T and P 827 ! the routine also needs the array rescale_tac initialized by init_rescale 828 !======================================================================= 829 integer(I4B), intent(in) :: lmax, m 830 real(DP), intent(in) :: cth, sth, mfac 831 real(DP), dimension(0:1,0:lmax), intent(in) :: recfac 832 real(DP), dimension( 0:lmax), intent(in) :: lam_fact 833 real(DP), dimension(1:3,0:lmax), intent(out) :: lam_lm 834 ! 835 real(DP) :: log2val, dlog2lg 836 real(DP) :: ovflow, unflow, corfac 837 real(DP) :: lam_mm, lam_0, lam_1, lam_2, lam_lm1m 838 integer(I4B) :: scalel, l, l_min 839 real(DP) :: normal_m, fm2, fl, flm1 840 real(DP) :: two_cth, one_on_s2, fm_on_s2, two_on_s2, c_on_s2 841 real(DP) :: a_w, a_x, b_w 842 !--------------------------------------------------------------- 843 844 ! define constants 845 ovflow = rescale_tab(1) 846 unflow = rescale_tab(-1) 847 l_min = l_min_ylm(m, sth) 848 dlog2lg = real(LOG2LG, kind=DP) 849 850 fm2 = real(m * m, kind = DP) 851 normal_m = (2.0_dp * m) * (1 - m) 852 two_cth = 2.0_dp * cth 853 one_on_s2 = 1.0_dp / (sth * sth) 854 fm_on_s2 = m * one_on_s2 855 two_on_s2 = 2.0_dp * one_on_s2 856 c_on_s2 = cth * one_on_s2 857 b_w = c_on_s2 858 859 860 ! computes lamba_mm 861 log2val = mfac + m*log(sth) * ALN2_INV ! log_2(lam_mm) 862 scalel = int (log2val / dlog2lg) 863 corfac = rescale_tab(max(scalel,RSMIN)) 864 lam_mm = 2.0_dp **(log2val - scalel * dlog2lg) ! rescaled lam_mm 865 if (IAND(m,1)>0) lam_mm = -lam_mm ! negative for odd m 866 867 lam_lm(1:3,m:lmax) = 0.0_dp 868 ! --- l = m --- 869 lam_lm(1,m) = corfac * lam_mm !Actual lam_mm 870 871 if (m >= l_min) then ! skip Ymm if too small 872 lam_lm(2,m) = (normal_m * lam_lm(1,m)) * ( 0.5_dp - one_on_s2 ) 873 lam_lm(3,m) = (normal_m * lam_lm(1,m)) * c_on_s2 874 endif 875 876 ! --- l > m --- 877 lam_0 = 0.0_dp 878 lam_1 = 1.0_dp 879 lam_2 = cth * lam_1 * recfac(0,m) 880 881 do l = m+1, lmax 882 ! do recursion 883 lam_lm1m = lam_lm(1,l-1) * lam_fact(l) ! must be incremented even if not used 884 lam_lm(1,l) = lam_2 * corfac * lam_mm 885 if (l >= l_min) then 886 fl = real(l, kind = DP) 887 flm1 = fl - 1.0_dp 888 a_w = two_on_s2 * (fl - fm2) + flm1 * fl 889 a_x = two_cth * flm1 890 lam_lm(2,l) = b_w * lam_lm1m - a_w * lam_lm(1,l) 891 lam_lm(3,l) = fm_on_s2 * ( lam_lm1m - a_x * lam_lm(1,l)) 892 endif 893 lam_0 = lam_1 * recfac(1,l-1) 894 lam_1 = lam_2 895 lam_2 = (cth * lam_1 - lam_0) * recfac(0,l) 896 897 ! do dynamic rescaling 898 if (abs(lam_2) > ovflow) then 899 lam_1 = lam_1*unflow 900 lam_2 = lam_2*unflow 901 scalel= scalel + 1 902 corfac = rescale_tab(max(scalel,RSMIN)) 903 elseif (abs(lam_2) < unflow .and. abs(lam_2) /= 0.0_dp) then 904 lam_1 = lam_1*ovflow 905 lam_2 = lam_2*ovflow 906 scalel= scalel - 1 907 corfac = rescale_tab(max(scalel,RSMIN)) 908 endif 909 910 enddo ! loop on l 911 end subroutine do_lam_lm_pol 912 !======================================================================= 913 subroutine gen_normpol(l_max, normal_l) 914 !======================================================================= 915 ! generates normalisaton factors for polarisation basis functions 916 ! for all l 917 !======================================================================= 918 integer(I4B), intent(IN) :: l_max 919 real(DP), intent(OUT), dimension(0:l_max) :: normal_l 920 ! 921 integer(I4B) :: l 922 real(DP) :: fl, xx 923 924 normal_l(0:1) = 0.0_dp 925 do l = 2, l_max 926 fl = DBLE(l) 927 xx = (fl+2.0_dp) * (fl+1.0_dp) * fl * (fl-1.0_dp) 928 normal_l(l) = SQRT ( KvS / xx) 929 ! either CMBFAST (KvS=1) or Kamionkowski et al. (KvS=2) definition 930 enddo 931 932 return 933 end subroutine gen_normpol 934 935 !================================================================ 936 function l_min_ylm(m, sth) result(lmin) 937 !================================================================ 938 ! returns minimal order l at which to keep Ylm 939 ! |Ylm| < eps * Y00 ==> 940 ! m_cut(theta, l) = theta * l * e / 2 + | ln(eps)| + ln(l)/2 941 ! if eps = 1.e-15 and l < 1.e4 942 ! m_cut(theta, l) = theta * l * 1.35 + 40 943 ! the choice of 1.35 (or larger) 944 ! also insures that the equatorial rings will have all their Ylm's computed 945 ! default parameters are HPX_MXL0 = 40 and HPX_MXL1 = 1.35_DP 946 !====================================================== 947 ! parameters of short-cut: defined in module header 948 ! dummy variables 949 integer(I4B) :: lmin 950 integer(I4B), intent(IN) :: m 951 real(DP), intent(IN) :: sth 952 953 lmin = m ! default 954 if (HPX_MXL0 > 0) lmin = max(lmin, int((m - HPX_MXL0)/(HPX_MXL1 * sth))) 955 956 return 957 end function l_min_ylm 958 !========================================================= 959 subroutine warning_oldbounds(cos_theta_cut, zbounds) 960 !========================================================= 961 real(DP), intent(in) :: cos_theta_cut 962 real(DP), intent(out), dimension(1:2) :: zbounds 963 964 if (cos_theta_cut <= 0.0_dp) then ! no cut 965 zbounds(1:2) = 0.0_dp 966 else 967 zbounds(1) = cos_theta_cut 968 zbounds(2) = - cos_theta_cut 969 endif 970 print*,' -------------------------------------' 971 print*,'WARNING: obsolete interface to MAP2ALM: ' 972 print*,' cos_theta_cut (6th argument) currently a DP scalar with value' 973 write(*,9000) ' ',cos_theta_cut 974 print*,' should now be replaced with a 2-element vector with values:' 975 write(*,9001) ' ',zbounds(1),zbounds(2) 976 print*,' See documentation for details.' 977 print*,' -------------------------------------' 9789000 format (a,g13.6) 9799001 format (a,g13.6,g13.6) 980 981 return 982 end subroutine warning_oldbounds 983 984 !************************************************************************** 985 ! 986 ! FOURIER TRANSFORM ON RINGS 987 ! 988 !************************************************************************** 989 !======================================================================= 990 subroutine ring_synthesis(nsmax,nlmax,nmmax,datain,nph,dataout,kphi0) 991 !======================================================================= 992 ! RING_SYNTHESIS 993 ! called by alm2map 994 ! calls real_fft 995 ! 996 ! dataout(j) = sum_m datain(m) * exp(i*m*phi(j)) 997 ! with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1 998 ! 999 ! as the set of frequencies {m} is larger than nph, 1000 ! we wrap frequencies within {0..nph-1} 1001 ! ie m = k*nph + m' with m' in {0..nph-1} 1002 ! then 1003 ! noting bw(m') = exp(i*m'*phi0) 1004 ! * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0)) 1005 ! with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m))) 1006 ! dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ] 1007 ! = Fourier Transform of bw 1008 ! is real 1009 ! 1010 ! NB nph is not necessarily a power of 2 1011 ! 1012 !======================================================================= 1013 1014 1015 INTEGER(I4B), INTENT(IN) :: nsmax, nlmax, nmmax 1016 INTEGER(I4B), INTENT(IN) :: nph, kphi0 1017 COMPLEX(DPC), DIMENSION(0:nmmax), INTENT(IN) :: datain 1018 REAL(DP), DIMENSION(0:nph-1), INTENT(OUT) :: dataout 1019 1020 INTEGER(I4B) :: iw,ksign,m,k,kshift 1021 COMPLEX(DPC), DIMENSION(0:nph-1) :: bw 1022 type(planck_fft2_plan) :: plan 1023 COMPLEX(DPC) :: dat 1024 real(DP) :: arg 1025 !======================================================================= 1026 1027 ksign = + 1 1028 kshift = (-1)**kphi0 ! either 1 or -1 1029 bw(0:nph-1) = CMPLX(0.0_dp, 0.0_dp, KIND=DP) 1030 1031 ! all frequencies [-m,m] are wrapped in [0,nph-1] 1032 bw(0)=datain(0) 1033 do m = 1, nmmax ! in -nmmax, nmmax 1034 iw = MODULO(m, nph) ! between 0 and nph-1 = m', F90 intrisic 1035 k = (m - iw) / nph ! number of 'turns' 1036 bw(iw) = bw(iw) + datain(m)*(kshift**k) ! complex number 1037 iw = MODULO(-m, nph) ! between 0 and nph-1 = m', F90 intrisic 1038 k = (-m - iw) / nph ! number of 'turns' 1039 bw(iw) = bw(iw) + CONJG(datain(m))*(kshift**k) ! complex number 1040 enddo 1041 1042 ! kshift**k = 1 for even turn numbers 1043 ! = 1 or -1 for odd turn numbers : results from the shift in space 1044 1045 ! applies the shift in position <-> phase factor in Fourier space 1046 dataout(0)=REAL(bw(0), kind=DP) 1047 do iw = 1, nph/2-1 1048 m = ksign*(iw) 1049 if(kphi0==1) then 1050 arg = m * PI / dble(nph) 1051 dat =bw(iw) * CMPLX( DCOS(arg), DSIN(arg), kind=DP) 1052 else 1053 dat =bw(iw) 1054 endif 1055 dataout(iw*2-1) = REAL(dat, kind=DP) 1056 dataout(iw*2 ) = AIMAG(dat) 1057 enddo 1058 iw=nph/2 1059 m = ksign*(iw) 1060 if(kphi0==1) then 1061 arg = m * PI / dble(nph) 1062 dat =bw(iw) * CMPLX( DCOS(arg), DSIN(arg), kind=DP) 1063 else 1064 dat =bw(iw) 1065 endif 1066 dataout(iw*2-1) = REAL(dat, kind=DP) 1067 1068 call make_fft2_plan(plan,length=nph,direction=fft2_backward) 1069 call real_fft2 (plan, dataout) 1070 ! ^^^^^^^^^^^^ 1071 call destroy_fft2_plan (plan) 1072 RETURN 1073 END subroutine ring_synthesis 1074 1075 !======================================================================= 1076 subroutine ring_analysis(nsmax,nlmax,nmmax,datain,nph,dataout,kphi0) 1077 !======================================================================= 1078 ! ring_analysis 1079 ! called by map2alm 1080 ! calls real_fft 1081 ! 1082 ! integrates (data * phi-dependence-of-Ylm) over phi 1083 ! --> function of m can be computed by FFT 1084 ! with 0<= m <= npoints/2 (: Nyquist) 1085 ! because the data is real the negative m are the conjugate of the 1086 ! positive ones 1087 !======================================================================= 1088 1089 INTEGER(I4B), INTENT(IN) :: nsmax, nlmax, nmmax 1090 INTEGER(I4B), INTENT(IN) :: nph, kphi0 1091 REAL(DP), DIMENSION(0:nph-1), INTENT(IN) :: datain 1092 COMPLEX(DPC), DIMENSION(0:nmmax), INTENT(OUT) :: dataout 1093 1094 INTEGER(I4B) :: i,m,im_max,ksign 1095! REAL(DP) :: phi0 1096 REAL(DP), DIMENSION(0:nph-1) :: data 1097 type(planck_fft2_plan) :: plan 1098 real(DP) :: arg 1099 1100 !======================================================================= 1101 1102 ksign = - 1 1103 data=0. 1104 data(0:nph-1)=datain(0:nph-1) 1105 1106 call make_fft2_plan(plan,length=nph,direction=fft2_forward) 1107 call real_fft2(plan,data) 1108 call destroy_fft2_plan(plan) 1109 1110 im_max = MIN(nph/2,nmmax) 1111 1112 ! m = 0, i=0 1113 dataout(0)=CMPLX(data(0),0.0_dp,kind=DP) 1114 1115 ! 1 <= m <= im_max, --> i=1,2,3,..., im_max 1116 do i = 1, im_max*2-3, 2 1117 dataout((i+1)/2) = CMPLX( data(i), data(i+1),kind= DP) 1118 enddo 1119 1120 if(im_max==nph/2) then 1121 dataout(im_max)= CMPLX( data(nph-1),0.0_dp,kind=DP) ! m = Nyquist 1122 else 1123 dataout(im_max)= CMPLX( data(2*im_max-1),data(2*im_max),kind=DP) 1124 endif 1125 1126 if(im_max==nmmax) goto 1000 ! m_max <= Nyquist 1127 1128 ! if (m > Nyquist) 1129 do i = im_max+1,min(nph,nmmax) 1130 dataout(i) = conjg(dataout(2*im_max-i) ) 1131 end do 1132 1133 if(min(nph,nmmax)==nmmax) goto 1000 ! nph > nmmax 1134 1135 do i = 2*im_max+1,nmmax 1136 dataout(i) = dataout(mod(i,2*im_max)) 1137 end do 1138 11391000 continue 1140 1141 if(kphi0==1)then 1142 do i = 0,nmmax 1143 m = ksign*i 1144! dataout(i)=dataout(i)* CONJG(trig(-m,nph/4)) 1145 arg = m * PI / dble(nph) 1146 dataout(i)=dataout(i)* CMPLX( DCOS(arg), DSIN(arg), kind=DP) 1147 enddo 1148 end if 1149 1150 RETURN 1151 END subroutine ring_analysis 1152 1153 !************************************************************************** 1154 ! ALM2MAP 1155 ! MAP2ALM 1156 ! ALM creation and alteration 1157 ! import overloaded routines 1158 !************************************************************************** 1159 1160 ! single precision routines 1161#include "alm_map_ss_inc.F90" 1162 1163 ! double precision routines 1164#include "alm_map_dd_inc.F90" 1165 1166 1167 !************************************************************************** 1168 ! 1169 ! PLM GENERATION 1170 ! 1171 !************************************************************************** 1172 !======================================================== 1173 subroutine plm_gen(nsmax, nlmax, nmmax, plm) 1174 !======================================================== 1175 use long_intrinsic, only: long_size 1176 1177 integer(i4b), intent(IN) :: nsmax, nlmax, nmmax 1178 real(dp), dimension(0:,1:), intent(OUT):: plm 1179 1180 INTEGER(I4B) :: l, m, ith, scalem, scalel, nd2, nrings 1181 integer(I8B) :: nd1, n_lm, n_plm, i_mm, i_up 1182 real(DP) :: lam_mm, lam_lm, lam_0, lam_1, lam_2, corfac, cth_ring 1183!! real(DP) :: lambda_w, lambda_x 1184 real(DP) :: normal_m, lam_lm1m 1185 real(DP) :: fm2, fl, flm1, fpol 1186 real(DP) :: a_w, b_w, a_x, fm_on_s2, two_on_s2, two_cth_ring 1187 real(DP) :: ovflow, unflow 1188 real(DP), dimension(:,:,:), allocatable :: plm_sub 1189 real(DP), dimension(:,:), allocatable :: recfac 1190 real(DP), dimension(:), allocatable :: lam_fact 1191 real(DP), dimension(:), allocatable :: mfac 1192 1193 real(DP), dimension(:), allocatable :: normal_l 1194 integer(i4b) :: nchunks, chunksize, ichunk, lchk, uchk, ithl 1195 integer(i4b) :: nph, kphi0, i 1196 integer(i8b) :: startpix 1197 real(DP), dimension(0:SMAXCHK) :: cth, sth 1198 real(DP), dimension(0:SMAXCHK) :: one_on_s2, c_on_s2 1199 1200 INTEGER(I4B) :: status 1201 LOGICAL(LGT) :: polarisation 1202 character(len=*), parameter :: code = 'PLM_GEN' 1203 1204 !================================================================= 1205 1206 ! Healpix definitions 1207 nrings = 2*nsmax ! number of isolatitude rings on N. hemisphere + equat 1208 1209 n_lm = ((nmmax+1)*(2*nlmax-nmmax+2))/2 !number of (l,m) with m in[0,M] and l in [m,L] 1210 n_plm = n_lm * nrings 1211 nd1 = long_size(plm, 1) 1212 nd2 = long_size(plm, 2) 1213 1214 if (nd1 < n_plm) then 1215 print*,code//' > Plm array too small:', nd1, n_plm 1216 stop 1217 endif 1218 if (nd2 /= 1 .and. nd2 /= 3) then 1219 print*,code//' > Plm array should have dimension 1 or 3:',nd2 1220 stop 1221 endif 1222 polarisation = (nd2 == 3) 1223 1224 ! --- allocates space for arrays --- 1225 nchunks = nrings/SMAXCHK + 1 ! number of chunks 1226 chunksize = (nrings+nchunks-1)/nchunks ! <= SMAXCHK 1227 1228 allocate(mfac(0:nmmax),stat = status) 1229 call assert_alloc(status,code,'mfac') 1230 if (polarisation) then 1231 allocate(normal_l(0:nlmax),stat = status) 1232 call assert_alloc(status,code,'normal_l') 1233 endif 1234 1235 if (.not. do_openmp()) then 1236 allocate(recfac(0:1,0:nlmax), plm_sub(1:nd2,0:nlmax,0:chunksize-1), stat = status) 1237 call assert_alloc(status,code,'recfac & plm_sub') 1238 if (polarisation) then 1239 allocate(lam_fact(0:nlmax),stat = status) 1240 call assert_alloc(status,code,'lam_fact') 1241 endif 1242 endif 1243 1244 ! ------------ initiate variables and arrays ---------------- 1245 1246 call gen_mfac(nmmax, mfac) 1247 ! generate Polarization normalisation 1248 if (polarisation) call gen_normpol(nlmax, normal_l) 1249 call init_rescale() 1250 ovflow = rescale_tab(1) 1251 unflow = rescale_tab(-1) 1252 plm = 0.0_dp 1253 1254 do ichunk = 0, nchunks-1 1255 lchk = ichunk * chunksize + 1 1256 uchk = min(lchk+chunksize - 1, nrings) 1257 1258 do ith = lchk, uchk 1259 ithl = ith - lchk !local index 1260 ! get pixel location information 1261 call get_pixel_layout(nsmax, ith, cth(ithl), sth(ithl), nph, startpix, kphi0) 1262 one_on_s2(ithl) = 1.0_dp / sth(ithl)**2 1263 c_on_s2(ithl) = cth(ithl) / sth(ithl)**2 1264 enddo 1265 1266!$OMP parallel default(none) & 1267!$OMP shared(nlmax, nmmax, lchk, uchk, nd2, chunksize, n_lm, & 1268!$OMP rescale_tab, ovflow, unflow, polarisation, & 1269!$OMP plm, cth, sth, mfac, normal_l, one_on_s2, c_on_s2) & 1270!$OMP private(plm_sub, recfac, lam_fact, status, & 1271!$OMP m, fm2, normal_m, ithl, ith, i_mm, i_up, & 1272!$OMP scalem, scalel, corfac, fpol, & 1273!$OMP lam_mm, lam_lm, lam_lm1m, lam_0, lam_1, lam_2, & 1274!$OMP cth_ring, fm_on_s2, two_on_s2, two_cth_ring, a_w, a_x, b_w, & 1275!$OMP l, fl, flm1, i) 1276 1277 if (do_openmp()) then 1278 ! allocate thread safe arrays 1279 allocate(recfac(0:1,0:nlmax), plm_sub(1:nd2,0:nlmax,0:chunksize-1), stat = status) 1280 call assert_alloc(status,code,'recfac & plm_sub') 1281 if (polarisation) then 1282 allocate(lam_fact(0:nlmax),stat = status) 1283 call assert_alloc(status,code,'lam_fact') 1284 endif 1285 endif 1286 1287!$OMP do schedule(dynamic,1) 1288 do m = 0, nmmax 1289 ! generate recursion factors (recfac) for Ylm of degree m 1290 call gen_recfac(nlmax, m, recfac) 1291 if (polarisation) then 1292 ! generate Ylm relation factor for degree m 1293 call gen_lamfac(nlmax, m, lam_fact) 1294 fm2 = real(m * m, kind = DP) 1295 normal_m = (2.0_dp * m) * (1 - m) 1296 endif 1297 1298 do ithl = 0, uchk - lchk 1299 ! determine lam_mm 1300 call compute_lam_mm(mfac(m), m, sth(ithl), lam_mm, corfac, scalem) 1301 ! ---------- l = m ---------- 1302 ! temperature 1303 lam_lm = corfac * lam_mm !Actual lam_mm 1304 plm_sub(1, m, ithl) = lam_lm 1305 1306 lam_0 = 0.0_dp 1307 lam_1 = 1.0_dp 1308 scalel=0 1309 cth_ring = cth(ithl) 1310 lam_2 = cth_ring * lam_1 * recfac(0,m) 1311 1312 if (polarisation) then 1313 fpol = normal_m * normal_l(m) * lam_lm 1314 plm_sub(2, m, ithl) = fpol * ( 0.5_dp - one_on_s2(ithl) ) 1315 plm_sub(3, m, ithl) = fpol * c_on_s2(ithl) 1316 ! 1317 fm_on_s2 = m * one_on_s2(ithl) 1318 two_on_s2 = 2.0_dp * one_on_s2(ithl) 1319 two_cth_ring = 2.0_dp * cth_ring 1320 b_w = c_on_s2(ithl) 1321 endif 1322 ! ---------- l > m ---------- 1323 do l = m+1, nlmax 1324 if (polarisation) lam_lm1m = lam_lm * lam_fact(l) 1325 lam_lm = lam_2 * corfac * lam_mm 1326 plm_sub(1, l, ithl) = lam_lm 1327 1328 if (polarisation) then 1329 fl = real(l, kind = DP) 1330 flm1 = fl - 1.0_dp 1331 a_w = two_on_s2 * (fl - fm2) + flm1 * fl 1332 a_x = two_cth_ring * flm1 1333 plm_sub(2, l, ithl) = ( b_w * lam_lm1m - a_w * lam_lm) * normal_l(l) 1334 plm_sub(3, l, ithl) = fm_on_s2 * ( lam_lm1m - a_x * lam_lm) * normal_l(l) 1335 endif 1336 1337 lam_0 = lam_1 * recfac(1,l-1) 1338 lam_1 = lam_2 1339 lam_2 = (cth_ring * lam_1 - lam_0) * recfac(0,l) 1340 if (abs(lam_2) > OVFLOW) then 1341 lam_1 = lam_1*UNFLOW 1342 lam_2 = lam_2*UNFLOW 1343 scalel= scalel + 1 1344 corfac = rescale_tab(max(scalel+scalem,RSMIN)) 1345 elseif (abs(lam_2) < UNFLOW) then 1346 lam_1 = lam_1*OVFLOW 1347 lam_2 = lam_2*OVFLOW 1348 scalel= scalel - 1 1349 corfac = rescale_tab(max(scalel+scalem,RSMIN)) 1350 endif 1351 enddo ! loop on l 1352 enddo ! loop on rings (ithl) 1353 1354 ! do memory skipping operations outside inner loops 1355 do ith = lchk, uchk 1356 i_mm = n_lm * (ith-1) + ((2*nlmax+3-m)*m)/2 ! location of Ym,m for ring ith 1357 i_up = i_mm + nlmax - m ! location of Ynlmax,m for ring ith 1358 ithl = ith - lchk 1359 do i=1, nd2 1360 plm(i_mm:i_up, i) = plm_sub(i, m:nlmax, ithl) 1361 enddo 1362 enddo 1363!!!!!!!! 1364! i_mm = n_lm*chunksize*ichunk + ((2*nlmax+3-m)*m)/2 1365! do ithl = 0, uchk-lchk 1366! i_up = i_mm+(nlmax-m+1)*ithl 1367! do i=1,nd2 1368! plm(i_up:i_up+nlmax-m,i) = plm_sub(i, m:nlmax, ithl) 1369! enddo 1370! enddo 1371 !------------------------------------------------ 1372 enddo ! loop on m 1373!$OMP end do 1374 if (do_openmp()) then 1375 deallocate (recfac, plm_sub) 1376 if (polarisation) deallocate(lam_fact) 1377 endif 1378!$OMP end parallel 1379 1380 1381 enddo ! loop on chunks 1382 1383 ! -------------------- 1384 ! free memory and exit 1385 ! -------------------- 1386 if (.not. do_openmp()) then 1387 deallocate (recfac, plm_sub) 1388 if (polarisation) deallocate(lam_fact) 1389 endif 1390 deallocate(mfac) 1391 if (polarisation) deallocate(normal_l) 1392 1393 return 1394 end subroutine plm_gen 1395 !************************************************************************** 1396 ! 1397 ! M I S C 1398 ! 1399 !************************************************************************** 1400 !======================================================================= 1401 subroutine pow2alm_units(units_pow, units_alm) 1402 !======================================================================= 1403 ! does the unit conversion between power and alm (and map) 1404 !======================================================================= 1405 character(len=*), dimension(1:), intent(in) :: units_pow 1406 character(len=*), dimension(1:), intent(out) :: units_alm 1407 1408 integer(kind=i4b) :: i, nu, j, ntemplate, idx 1409 character(len=80) :: uu, ui, uo 1410 character(len=5),dimension(1:7) :: template = & 1411 & (/ "_SQUA","-SQUA","SQUA ","^2 ","^ 2 ","**2 ","** 2 " /) 1412 !======================================================================= 1413 nu = min( size(units_pow), size(units_alm)) 1414 ntemplate = size(template) 1415 1416 units_alm = "" 1417 do i = 1, nu 1418 ui = adjustl(units_pow(i)) 1419 uu = trim(strupcase(ui)) 1420 uo = "unknown" ! default 1421 do j = 1, ntemplate 1422 idx = index(uu, template(j)) 1423 if (idx > 0) then 1424 uo = ui(1:idx-1) 1425 exit 1426 endif 1427 enddo 1428 units_alm(i) = uo 1429 ! print*,i,trim(units_pow(i))," -> ",trim(units_alm(i)) 1430 enddo 1431 1432 return 1433 end subroutine pow2alm_units 1434 !**************************************************************************** 1435 subroutine generate_beam(fwhm_arcmin, lmax, gb, beam_file) 1436 !========================================================================== 1437 ! 1438 ! generate_beam(fwhm_arcmin, lmax, gb [, beam_file]) 1439 ! generate the beam window function gb up to lmax 1440 ! either a gaussian of FHWM : fwhm_arcmin (in arcmin) 1441 ! or the function read from beam_file 1442 ! 1443 ! July 2003, EH : replicate temperature beam for polarization when reading 1444 ! standard ASCII file 1445 ! Sep 2014, EH: accept ASCII files with 1,2 or 3 columns (on top of multipoles) 1446 ! read all requested columns from FITS file 1447 ! FITSIO extended file syntax (eg file.fits[col 1]) must be used to select row(s) 1448 ! missing columns will be copied from existing ones 1449 !========================================================================== 1450 use fitstools, only : getsize_fits, fits2cl 1451 real(kind=DP), intent(in) :: fwhm_arcmin 1452 real(kind=DP), dimension(0:,1:), intent(out) :: gb 1453 integer(kind=I4B), intent(in) :: lmax 1454 character(len=*), optional, intent(in) :: beam_file 1455 1456 character(len=FILENAMELEN) :: new_beam_file 1457 logical(kind=lgt) :: extfile, found_unix 1458 integer(kind=i4b) :: type, nl, nd, lunit, il, i, l100, iline, junk, nc 1459 character(len=80), dimension(1:180) :: header 1460 character(len=1600) :: str 1461 character(len=80) :: card 1462 real(dp) :: bref 1463 real(dp), dimension(1:3) :: buffer 1464 character(len=*), parameter :: code = 'generate_beam' 1465 !========================================================================== 1466 ! test if name of external is given and valid 1467 extfile = .false. 1468 if (present(beam_file)) then 1469 extfile = (trim(beam_file) /= "") 1470 endif 1471 1472 lunit = 15 1473 nl = size(gb, 1) 1474 nd = size(gb, 2) 1475 gb = 0.0_dp 1476 1477 if (nl <= lmax) then 1478 print 9000,' Generate_beam: beam array only available up to ',nl 1479 endif 1480 nl = min(nl, lmax+1) 1481 1482 if (extfile) then 1483 call assert_present(beam_file) 1484 new_beam_file = beam_file 1485 print 9000,' Reading beam information from '//trim(new_beam_file) 1486 1487 ! find out file nature 1488 type = 1 ! FITS by default 1489 inquire(file=new_beam_file, exist=found_unix) 1490 if (found_unix) then 1491 open(unit=lunit,file=new_beam_file,status='old', & 1492 & form='formatted',action='read') 1493 read(lunit,'(a)') card 1494 close(lunit) 1495 card = adjustl(card) 1496 if (card(1:8) /= 'SIMPLE ' .AND. card(1:8) /= 'XTENSION') type = -1 1497 endif 1498 1499 ! read file according to its type 1500 if (type < 0) then 1501 ! ordinary ascii file ? 1502 lunit = 32 1503 iline = 0 ! line index 1504 open(unit=lunit,file=new_beam_file,status='old', & 1505 & form='formatted',action='read') 1506 do 1507 iline = iline+1 1508 read(lunit,'(a)', end=100, err=100) str 1509 if (len_trim(str)>0 .and. str(1:1) /= '#') then 1510 buffer = 0_dp 1511 1512 ! try 3 columns 15133 continue 1514 read(str,*, err=2, end=2) il, buffer(1:3) 1515 nc = 3 1516 goto 99 1517 1518 ! try 2 columns 15192 continue 1520 read(str,*, err=1, end=1) il, buffer(1:2) 1521 nc = 2 1522 goto 99 1523 1524 ! try 1 column 15251 continue 1526 read(str,*, err=666, end=666) il, buffer(1) 1527 nc = 1 1528 goto 99 1529 1530 ! failure 1531666 continue 1532 print 9000, 'ERROR in '//code 1533 print 9000, 'Unable to parse ASCII file '//trim(new_beam_file) 1534 print 9000, 'at line '//trim(adjustl(string(iline)))//' containing' 1535 print 9000, trim(str) 1536 print 9000, 'Expected format :' 1537 print 9000, ' Multipole, column_1 [, column_2, column_3]' 1538 print 9000, ' with OPTIONAL colum_2 and column_3' 1539 call fatal_error 1540 1541 ! record data, keep going until end of file 154299 continue 1543 gb(il, 1:nd) = buffer(1:nd) 1544 if (il == nl-1) exit 1545 endif 1546 enddo 1547100 continue 1548 close(lunit) 1549 if (il < (nl-1)) then 1550 print 9001,' WARNING: Beam transfer function only available up to l= ',il ! 1551 print 9000,' The larger multipoles will be set to 0' 1552 endif 1553 1554 1555 else if (type == 1) then 1556 ! FITS file with ascii or binary table 1557 call fits2cl(new_beam_file, gb, nl-1_i4b, nd, header, fmissval=0.0_dp) 1558 junk = getsize_fits(new_beam_file, nmaps=nc) 1559 else 1560 print 9000,' the file '//trim(new_beam_file) & 1561 & //' is of unknown type.' 1562 call fatal_error 1563 endif 1564 1565 ! if Grad absent, replicate Temperature; if Curl absent, replicate Grad 1566 if (nc < nd) then 1567 print 9000,'Not enough columns found in '//trim(new_beam_file) 1568 print 9000,'Expected '//trim(adjustl(string(nd)))& 1569 //', found '//trim(adjustl(string(nc))) 1570 do i=nc+1, nd 1571 print 9000,' column #'//trim(adjustl(string(i))) & 1572 & //' empty, fill in with column #' & 1573 & //trim(adjustl(string(i-1))) 1574 gb(:,i) = gb(:,i-1) 1575 enddo 1576 else 1577 print 9000,'Read '//trim(adjustl(string(nd)))& 1578 //' columns from '//trim(new_beam_file) 1579 print 9000,'(out of '//trim(adjustl(string(nc)))//')' 1580 endif 1581 1582! ! if Grad absent, replicate Temperature; if Curl absent, replicate Grad 1583! l100 = min(100, nl-1) 1584! bref = 1.e-4*sum(abs(gb(0:l100,1))) 1585! do i=2,nd 1586! if ( sum(abs(gb(0:l100,i))) < bref ) then 1587! print 9002,' column #',i,' empty, fill in with column #',i-1 1588! gb(:,i) = gb(:,i-1) 1589! endif 1590! enddo 1591 1592 else 1593 ! gaussian beam 1594 print*,'Generating gaussian beam of FHWM [arcmin] = ',fwhm_arcmin 1595 call gaussbeam(fwhm_arcmin, nl-1, gb) 1596 endif 1597 15989000 format(a) 15999001 format(a,i6) 16009002 format(a,i3,a,i3) 1601 1602 return 1603 end subroutine generate_beam 1604 !************************************************************* 1605 subroutine gaussbeam(fwhm_arcmin, lmax, gb) 1606 !=========================================================== 1607 ! gaussbeam(fwhm_arcmin, gb) 1608 ! returns in gb the window function on [0,lmax] corresponding 1609 ! to the gaussian beam of FWHM = fwhm_arcmin 1610 ! The returned beam function has up to 3 components, 1611 ! gb(*,1) = bt : temperature 1612 ! gb(*,2) = bt * exp(2*sigma^2) : grad 1613 ! gb(*,3) = bt * exp(2*sigma^2) : curl 1614 ! with sigma = gaussian rms in radian 1615 ! and bt = exp( l(l+1) * sigma^2 / 2) 1616 !=========================================================== 1617 real(kind=DP), intent(in) :: fwhm_arcmin 1618 real(kind=DP), dimension(0:,1:), intent(out) :: gb 1619 integer(kind=I4B), intent(in) :: lmax 1620 1621 integer(kind=I4B) :: l, nd 1622 real(kind=DP) :: sigma, arcmin2rad, sigma2fwhm, fact_pol 1623 !=========================================================== 1624 1625 call assert (fwhm_arcmin>=0.0_dp,'fwhm of gaussian beam should be positive') 1626 1627 nd = size(gb,2) 1628 1629 arcmin2rad = PI / (180.0_dp * 60.0_dp) 1630 sigma2fwhm = sqrt(8.0_dp * log(2.0_dp)) 1631 1632 sigma = fwhm_arcmin * arcmin2rad / sigma2fwhm ! in radians 1633 1634 fact_pol = exp(2.0_dp*sigma**2) ! correction for polarised fields 1635 1636 ! temperature 1637 do l=0,lmax 1638 gb(l,1) = exp(-.5_dp * l*(l+1.0_dp) * sigma**2) 1639 enddo 1640 ! electric or gradient 1641 if (nd > 1) gb(0:lmax,2) = gb(0:lmax,1) * fact_pol 1642 ! magnetic or curl 1643 if (nd > 2) gb(0:lmax,3) = gb(0:lmax,1) * fact_pol 1644 1645 return 1646 end subroutine gaussbeam 1647 1648 !======================================================================= 1649 subroutine pixel_window(pixlw, nside, windowfile) 1650 !======================================================================= 1651 !======================================================================= 1652 1653 use fitstools, only : getsize_fits, read_bintab 1654 use extension, only : getEnvironment 1655 use pix_tools, only : nside2npix 1656 use paramfile_io, only: get_healpix_data_dir, scan_directories, & 1657 & get_healpix_pixel_window_file 1658 1659 real(DP), dimension(0:,1:), intent(OUT) :: pixlw 1660 character(len=*), intent(IN), optional :: windowfile 1661 integer(i4b), intent(IN), optional :: nside 1662 1663 character(len=filenamelen) :: wfile, wfile_def, healpixdatadir 1664 character(len=4) :: sstr 1665 integer(I4B) :: npw_file, ncolw_file 1666 integer(I8B) :: npix 1667 integer(I4B) :: npw, ncolw, i 1668 logical(LGT) :: good, bad 1669 REAL(DP) :: nullval 1670 character(len=*), parameter :: code = 'Pixel_Window > ' 1671 real(DP), dimension(:,:), allocatable :: pixtmp 1672 !======================================================================= 1673 1674 call assert (present(nside) .or. present(windowfile), & 1675 code//'Nside or windowfile have to be present') 1676 1677 npw = size(pixlw, dim=1) 1678 ncolw = size(pixlw, dim=2) 1679 1680 if (present(windowfile)) then 1681 wfile = trim(windowfile) 1682 inquire(file=wfile, exist = good) 1683 else 1684 if (nside == 0) then 1685 pixlw = 1.0_dp 1686 return 1687 endif 1688 npix = nside2npix(nside) 1689 if (npix < 0) then 1690 print*,code//'invalid Nside = ',nside 1691 call fatal_error 1692 endif 1693 healpixdatadir = get_healpix_data_dir() 1694 wfile_def = get_healpix_pixel_window_file(nside) 1695 good = scan_directories(healpixdatadir, wfile_def, wfile) 1696 endif 1697 if (.not.good) then 1698 print*,code//trim(wfile)//' not found' 1699 call fatal_error 1700 endif 170110 continue 1702 1703 npw_file = getsize_fits(wfile, nmaps = ncolw_file) 1704 if (ncolw_file > 3) then 1705 print*,code//' Too many columns in '//trim(wfile) 1706 call fatal_error 1707 endif 1708 ! allocate tmp array large enough to read the whole file 1709 allocate(pixtmp(0:npw_file-1, 1:ncolw_file)) 1710 1711 if (npw_file < npw) then 1712 print*,'Only ',npw_file,' multipoles available in '//trim(wfile)//', expected ',npw 1713 endif 1714 1715! call read_dbintab(windowfile, pixlw, npw, ncolw_file, nullval, bad) 1716! edited 2004-05-28 1717! call read_dbintab(wfile, pixlw, npw, ncolw_file, nullval, bad) 1718 1719! call read_dbintab(wfile, pixtmp, npw_file, ncolw_file, nullval, bad) 1720 call read_bintab(wfile, pixtmp, npw_file, ncolw_file, nullval, bad) 1721 1722 npw = min(npw_file, npw) 1723 ncolw_file = min(ncolw_file, ncolw) 1724 1725 do i=1, ncolw_file 1726 pixlw(0:npw-1, i) = pixtmp(0:npw-1, i) 1727 enddo 1728 if (ncolw_file == 1) then ! T only in file 1729 if (ncolw > 1) pixlw(:,2) = pixlw(:,1) ! w_G = w_T 1730 if (ncolw > 2) pixlw(:,3) = pixlw(:,1) ! w_C = w_T 1731 else if (ncolw_file == 2) then ! T and G in file 1732 if (ncolw > 2) pixlw(:,3) = pixlw(:,2) ! w_C = w_G 1733 endif 1734 1735 deallocate(pixtmp) 1736 1737 return 1738 end subroutine pixel_window 1739 1740 1741!*********************************************************************************** 1742 1743end module alm_tools 1744