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!----------------------------------------------------------------------------- 28!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 29! Fast Smoothing of a full sky map in HEALPIX pixelisation. 30! Written and developed by E. Hivon (hivon@iap.fr) and K. Gorski 31! (krzysztof.m.gorski@jpl.nasa.gov) based on HEALPIX pixelisation developed 32! by K. Gorski 33! 34! Copyright 1997 by Eric Hivon and Krzysztof M. Gorski. 35! All rights reserved. 36! 37!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 38!======================================================================= 39 !======================================================================= 40 ! EXTERNAL LIBRARY: 41 ! this code uses the FITSIO library that can be found at 42 ! http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html 43 ! 44 ! RELATED LITTERATURE: 45 ! about HEALPIX : see Gorski et al, 2005, ApJ, 622, 759 46 ! about this code : see Hivon & Gorski, 1997, in preparation 47 ! 48 ! HISTORY: 49 ! April-October 1997, Eric Hivon, TAC (f77) 50 ! Jan 1998 : translation in Fortran 90 51 ! Dec 2001 : version 1.2, EH, IPAC 52 ! implementation of non-interactive input 53 ! extension to polarised maps 54 ! Sep 2002 : implement new parser 55 ! Feb 2003 : output correct ORDERING keyword 56 ! May 2007 : calls map2alm_iterative 57 ! Jun 2010 : supports large maps 58 ! 59 ! FEEDBACK: 60 ! for any questions : hivon@iap.fr 61 ! 62 !======================================================================= 63 ! version 2.2 64 !======================================================================= 65 ! this file can not be compiled on its own. 66 ! It must be inserted into the file smoothing.f90 by the command include 67 ! 68 integer(I4B), parameter :: KALMC = KIND((1.0_KMAP, 1.0_KMAP)) ! precision of alm arrays 69 !------------------------------------------------------------------------- 70 71 integer(I4B) :: nsmax, nlmax, nmmax 72 73 complex(kind=KALMC), DIMENSION(:,:,:), ALLOCATABLE :: alm_TGC 74 complex(kind=KALMC), DIMENSION(:,:,:), ALLOCATABLE :: alm_TGC2 75 real(kind=KMAP), DIMENSION(:,:), ALLOCATABLE :: map_TQU 76 real(kind=KMAP), DIMENSION(:,:), ALLOCATABLE :: map_TQU2 77 real(kind=DP), DIMENSION(:,:), ALLOCATABLE :: w8ring_TQU 78 real(kind=DP), DIMENSION(:,:), ALLOCATABLE :: dw8 79 real(kind=DP), DIMENSION(:,:), ALLOCATABLE :: plm 80! real(kind=KMAP), dimension(:,:), ALLOCATABLE :: pw8list 81 real(kind=KMAP), dimension(:,:), target, ALLOCATABLE :: pw8map 82 real(kind=KMAP), dimension(:,:), pointer :: pmask 83 integer(kind=I4B) :: status 84 85 integer(kind=I8B) :: npixtot, n_plm, ipix, npixtot_tmp !, n_pw8_read 86 integer(kind=I4B) :: mlpol 87 integer(kind=I4B) :: plm_nside, plm_lmax, plm_pol, plm_mmax 88 !integer(kind=I4B) :: nmaps, type, polar_fits 89 integer(kind=I4B) :: i, j 90 integer(kind=I4B) :: nlheader 91 integer(kind=I4B) :: ordering, order_type 92 integer(kind=I4B) :: iter, iter_order, polar, n_pols 93 integer(kind=I4B) :: simul_type 94 integer(kind=I4B) :: n_rings, n_rings_read, won, nmw8 95 96 real(kind=KMAP) :: fwhm_arcmin, fwhm_deg 97 real(kind=KMAP) :: fmissval 98 real(kind=SP), parameter :: fbad_value = -1.6375e30_sp 99 real(kind=SP) :: pix_size_arcmin 100 real(kind=DP) :: theta_cut_deg, cos_theta_cut 101 real(kind=DP) :: nullval 102 real(kind=DP), dimension(1:2) :: zbounds 103 real(kind=DP) :: fsky 104 real(kind=SP) :: clock_time, time0, time1 105 real(kind=SP) :: ptime, ptime0, ptime1 106 107 character(LEN=filenamelen) :: infile 108 character(LEN=filenamelen) :: outfile 109! character(LEN=filenamelen) :: parafile = '' 110 character(LEN=filenamelen) :: infile_plm 111 character(LEN=filenamelen) :: infile_w8 112 character(LEN=filenamelen) :: w8name 113 character(LEN=filenamelen) :: def_dir, def_file 114 character(LEN=filenamelen) :: usr_dir, usr_file 115 character(LEN=filenamelen) :: final_file 116 character(LEN=filenamelen) :: healpixtestdir 117 character(LEN=filenamelen) :: beam_file 118 character(len=filenamelen) :: description 119 character(len=100) :: chline, chline1, strw8file 120 character(LEN=80), DIMENSION(1:180) :: header 121 LOGICAL(kind=LGT) :: bad, ok, polarisation 122! character(LEN=*), PARAMETER :: code = 'SMOOTHING' 123 character(len=*), parameter :: VERSION = HEALPIX_VERSION 124 integer(kind=i4b), parameter :: nm_max = 5 125 character(len=80), dimension(1:nm_max) :: units_map 126 character(len=20) :: coordsys 127 128 type(paramfile_handle) :: handle 129 130 real(kind=DP), dimension(0:3) :: mono_dip 131 real(kind=DP) :: theta_dip, phi_dip, long_dip, lat_dip 132 integer(kind=I4B) :: lowlreg 133 134 !----------------------------------------------------------------------- 135 ! get input parameters and create arrays 136 !----------------------------------------------------------------------- 137 138 call wall_clock_time(time0) 139 call cpu_time(ptime0) 140 ! --- read parameters interactively if no command-line arguments 141 ! --- are given, otherwise interpret first command-line argument 142 ! --- as parameter file name and read parameters from it: 143 144 ! --- announces program, default banner *** 145 print *, " " 146 print *,' '//code//' '//version 147 write(*,'(a)') ' *** Smoothing of a Temperature (+Polarisation) map *** ' 148 if (KMAP == SP) print*,' Single precision outputs' 149 if (KMAP == DP) print*,' Double precision outputs' 150 print *, '' 151 handle = parse_init(parafile) 152 153 ! --- choose temp. only or temp. + pol. --- 154 description = concatnl( & 155 & " Do you want to analyse", & 156 & " 1) a Temperature only map", & 157 & " 2) Temperature + polarisation maps") 158 simul_type = parse_int(handle, 'simul_type', default=1, vmin=1, vmax=2, descr=description) 159 if (simul_type .eq. 1) polarisation = .false. 160 if (simul_type .eq. 2) polarisation = .true. 161 162 ! --- gets the file name for the map --- 163 chline = "''" 164 healpixtestdir = get_healpix_test_dir() 165 if (trim(healpixtestdir)/='') chline = trim(healpixtestdir)//'/map.fits' 166 description = concatnl( & 167 & "", & 168 & " Enter input file name (Map FITS file): ") 169 infile = parse_string(handle, 'infile',default=chline, descr=description, filestatus='old') 170 171 ! --- check that this map is usable for this analysis 172 call check_input_map(code, infile, polarisation) ! this may change polarisation from T to F 173 174 ! --- finds out the pixel number of the map and its ordering --- 175 ! npixtot = getsize_fits(infile, nmaps = nmaps, ordering=ordering, nside=nsmax, & 176 ! & mlpol=mlpol, type = type, polarisation = polar_fits, coordsys = coordsys) 177 178 ! --- get other relevant information --- 179 npixtot_tmp = getsize_fits(infile, ordering=ordering, nside=nsmax, & 180 mlpol=mlpol, coordsys = coordsys) 181 182! if (nsmax<=0) then 183! print*,'Keyword NSIDE not found in FITS header!' 184! call fatal_error(code) 185! endif 186! if (type == 3) npixtot = nside2npix(nsmax) ! cut sky input data set 187! if (nsmax/=npix2nside(npixtot)) then 188! print*,"FITS header keyword NSIDE does not correspond" 189! print*,"to the size of the map!" 190! call fatal_error(code) 191! endif 192 193! if (polarisation .and. (nmaps >=3) .and. polar_fits == -1) then 194! print*,"The input fits file MAY NOT contain polarisation data." 195! print*,"Proceed at your own risk" 196! endif 197 198! if (polarisation .and. (nmaps<3 .or. polar_fits ==0)) then 199! print *,"The file does NOT contain polarisation maps" 200! print *,"only the temperature field will be smoothed" 201! polarisation = .false. 202! endif 203 204! ! --- check ordering scheme --- 205! if ((ordering/=1).and.(ordering/=2)) then 206! print*,'The ordering scheme of the map must be RING or NESTED.' 207! print*,'No ordering specification is given in the FITS-header!' 208! call fatal_error(code) 209! endif 210 211 order_type = ordering 212 213 ! --- gets the L range for the analysis --- 214 WRITE(chline1,"(a,i5)") " The map has Nside = ",nsmax 215 if (mlpol<=0) WRITE(chline,"(a,i5,a)") "We recommend: (0 <= l <= l_max <= ",3*nsmax-1,")" 216 if (mlpol> 0) WRITE(chline,"(a,i5,a)") "We recommend: (0 <= l <= l_max <= ",mlpol,")" 217 description = concatnl(& 218 & "", & 219 & chline1, & 220 & " Enter the maximum l range (l_max) for the simulation. ", & 221 & chline ) 222 nlmax = parse_int(handle, 'nlmax', default=2*nsmax, descr=description) 223 224 ! ---------- default parameters ---------------- 225 nmmax = nlmax 226 npixtot = nside2npix(nsmax) 227 pix_size_arcmin = 360.*60./SQRT(npixtot*PI) 228 229 ! --- ask about monopole/dipole removal --- 230 description = concatnl(& 231 & "", & 232 & " To improve the temperature map analysis (specially on a cut sky) ", & 233 & " you have the option of first regressing out the best fit monopole and dipole", & 234 & " NOTE : the regression is made on valid (unflagged) pixels, ", & 235 & " out of the symmetric cut (if any)", & 236 & " Do you want to : ", & 237 & " 0) Do the analysis on the raw map", & 238 & " 1) Remove the monopole (l=0) first ", & 239 & " 2) Remove the monopole and dipole (l=1) first") 240 lowlreg = parse_int(handle, 'regression', vmin=0, vmax=2, default=0, descr=description) 241 242 243 ! --- gets the cut applied to the map --- 244 description = concatnl(& 245 & "", & 246 & " Enter the symmetric cut around the equator in DEGREES : ", & 247 & " (One ignores data within |b| < b_cut) 0 <= b_cut = ") 248 theta_cut_deg = parse_double(handle, 'theta_cut_deg', & 249 & vmin=0.0_dp, default=0.0_dp, descr=description) 250 cos_theta_cut = SIN(theta_cut_deg/180.d0*PI) !counted from equator instead of from pole 251 252 zbounds = (/ cos_theta_cut , -cos_theta_cut /) 253 if (theta_cut_deg<1e-4) zbounds = (/ -1.0_dp, 1.0_dp /) !keep all sphere if cut not set 254 fsky = (zbounds(2)-zbounds(1))/2.0_dp 255 if (fsky <= 0.0_dp) fsky = 1.0_dp + fsky 256 write(*,"(a,f6.1,a)") "One keeps ",100.*fsky," % of the original map" 257 258 ! --- gets the order of iteration --- 259 description = concatnl(& 260 & "", & 261 & " Do you want : ", & 262 & " 0) a standard analysis", & 263 & " 1,2,3,4....) an iterative analysis", & 264 & " (enter order of iteration, 3rd order is usually optimal)") 265 iter_order=parse_int(handle, 'iter_order', vmin=0, default=0, descr=description) 266 if (.not. handle%interactive) then 267 select case (iter_order) 268 case (0) 269 print*," Standard analysis" 270 case default 271 print*," Iterative analysis" 272 end select 273 end if 274 275 ! ------------------- precomputed plms --------------------- 27640 continue 277 n_plm = 0 278 plm_pol = 0 279 description = concatnl(& 280 & "", & 281 & " Enter name of file with precomputed p_lms :", & 282 & "(eg, plm.fits) (if '' is typed, no file is read)" ) 283 infile_plm = parse_string(handle, 'plmfile', default="''", descr=description, filestatus='old') 284 285 if (trim(infile_plm) /= "") then 286 ! check that the plm file matches the current simulation (for nside,lmax,mmax) 287 call read_par(infile_plm, plm_nside, plm_lmax, plm_pol, mmax=plm_mmax) 288 289 if ((plm_pol /= 1).and.(.not.polarisation)) plm_pol=1 290 291 if (plm_nside /= nsmax .or. plm_lmax /= nlmax .or. plm_mmax /= nmmax) then 292 print *," "//code//"> Plm file does not match map to smooth in Nside, Lmax or Mmax ! " 293 print*,'nside (plm file, map) = ',plm_nside,nsmax 294 print*,'lmax = ',plm_lmax, nlmax 295 print*,'mmax = ',plm_mmax, nmmax 296 if (handle%interactive) goto 40 297 call fatal_error(code) 298 endif 299 n_plm = (nmmax + 1_i8b)*(2*nlmax + 2_i8b -nmmax)*nsmax 300 print *," " 301 endif 302 303 ! --- inputs the FWHM or beam file --- 304 beam_file = '' 305 description = concatnl(& 306 & "", & 307 & " Enter the Gaussian beam FWHM in arcmin >= 0 (eg: 420.) : ") 308 fwhm_arcmin = parse_double(handle, 'fwhm_arcmin', default=420.0_dp, vmin=0.0_dp,descr=description) 309 310 description = concatnl(& 311 & "", & 312 & " Enter an external file name containing ", & 313 & " a symmetric beam Legendre transform (eg: mybeam.fits)", & 314 & " NB: if set to an existing file, it will override the FWHM chosen above", & 315 & " if set to '', the gaussian FWHM will be used.") 316 beam_file = parse_string(handle, 'beam_file', default="''", & 317 & descr=description, filestatus='old') 318 if (beam_file /= '') then 319 fwhm_arcmin = 0. 320 print*,'fwhm_arcmin is now : 0.' 321 endif 322 323 ! --- gets the output sky map filename --- 324 description = concatnl(& 325 & "", & 326 & " Enter Output map file name (eg, test_smooth.fits) :", & 327 & " (or !test_smooth.fits to overwrite an existing file)" ) 328 outfile = parse_string(handle, "outfile", & 329 default="!test_smooth.fits", descr=description, filestatus="new") 330 331 !----------------------------------------------------------------------- 332 ! ask for weights 333 !----------------------------------------------------------------------- 334 description = concatnl(& 335 & "", & 336 & " Do you want to use quadrature weights in the analysis: ", & 337 & " 0: no weight ", & 338 & " 1: ring-based weights ", & 339 & " 2: pixel-based weights ?") 340 won = parse_int(handle, 'won', vmin=0, vmax=2, default=0, descr=description) 341 342 infile_w8="" 343 strw8file = "" 344 345 if (won == 1 .or. won == 2) then 346 347 ! default weight file name 348 w8name = get_healpix_weight_file(nsmax, won) 349 350 def_file = trim(w8name) 351 def_dir = get_healpix_data_dir() 352 353 if (won == 1) strw8file = "ring-based quadrature weight file" 354 if (won == 2) strw8file = "pixel-based quadrature weight file" 355 35622 continue 357 final_file = '' 358 ok = .false. 359 ! if interactive, try default name in default directories first 360 if (handle%interactive) ok = scan_directories(def_dir, def_file, final_file) 361 if (.not. ok) then 362 ! not found, ask the user 363 description = concatnl("",& 364 & " Could not find "//trim(strw8file), & 365 & " Enter the directory where this file can be found:") 366 usr_dir = parse_string(handle,'w8filedir',default='',descr=description) 367 if (trim(usr_dir) == '') usr_dir = trim(def_dir) 368 description = concatnl("",& 369 & " Enter the name of the "//trim(strw8file)//":") 370 usr_file = parse_string(handle,'w8file',default=def_file,descr=description) 371 ! look for new name in user provided or default directories 372 ok = scan_directories(usr_dir, usr_file, final_file) 373 ! if still fails, crash or ask again if interactive 374 if (.not. ok) then 375 print*,' File not found:'//trim(usr_file) 376 if (handle%interactive) goto 22 377 call fatal_error(code) 378 endif 379 endif 380 infile_w8 = final_file 381 endif 382 383 print *," " 384 call parse_check_unused(handle, code=lcode) 385 call parse_summarize(handle,code=lcode,prec=KMAP) 386 call parse_finish(handle) 387! call parse_summarize(handle) 388 call brag_openmp() 389 390 !----------------------------------------------------------------------- 391 ! allocate space for arrays 392 !----------------------------------------------------------------------- 393 polar = 0 394 if (polarisation) polar = 1 395 n_pols = 1 + 2*polar ! either 1 or 3 396 397 ALLOCATE(map_TQU(0:12*nsmax**2-1,1:n_pols),stat = status) 398 call assert_alloc(status,code,"map_TQU") 399 400 ALLOCATE(w8ring_TQU(1:2*nsmax,1:n_pols),stat = status) 401 call assert_alloc(status,code,"w8ring_TQU") 402 403 ALLOCATE(dw8(1:2*nsmax,1:n_pols),stat = status) 404 call assert_alloc(status,code,"dw8") 405 406 if (n_plm/=0) then 407 ALLOCATE(plm(0:n_plm-1,1:plm_pol),stat = status) 408 call assert_alloc(status,code,"plm") 409 endif 410 !----------------------------------------------------------------------- 411 ! reads the map 412 !----------------------------------------------------------------------- 413 print *," "//code//"> Inputting original map " 414 415 fmissval = 0.0 416 if (lowlreg > 0) fmissval = fbad_value 417 call input_map(infile, map_TQU(0:,1:n_pols), & 418 & npixtot, n_pols, fmissval=fmissval, units=units_map(1:n_pols)) 419 do i=1,n_pols 420 units_map(i) = adjustl(units_map(i)) 421 if (trim(units_map(i)) == '') units_map(i) = 'unknown' 422 enddo 423 424 !----------------------------------------------------------------------- 425 ! remove dipole 426 !----------------------------------------------------------------------- 427 428 if (lowlreg > 0) then 429 PRINT *," "//code//"> Remove monopole (and dipole) from Temperature map" 430 call remove_dipole(nsmax, map_TQU(0:npixtot-1,1), ordering, lowlreg, & 431 & mono_dip(0:), zbounds, fmissval=fmissval) 432 433 if (fmissval /= 0.0) then 434 do j=1,n_pols 435 do ipix=0, npixtot-1 436! if (abs(map_TQU(i,j)/fmissval-1.0) < 1.e-6*fmissval) then 437 if (abs(map_TQU(ipix,j)-fmissval) <= abs(1.e-6*fmissval)) then 438 map_TQU(ipix,j) = 0.0_sp 439 endif 440 enddo 441 enddo 442 endif 443 444 write(unit=*,fmt="(a,g13.3,a)") " Monopole = ",& 445 & mono_dip(0)," "//trim(units_map(1)) 446 if (lowlreg > 1) then 447 call vec2ang( mono_dip(1:3), theta_dip, phi_dip) 448 write(unit=*,fmt="(a,g13.3,a)") " Dipole = ",& 449 & sqrt(sum(mono_dip(1:3)**2))," "//trim(units_map(1)) 450 long_dip = phi_dip /PI*180.0 451 lat_dip = 90.0-theta_dip/PI*180.0 452 write(unit=*,fmt="(a,g10.3,', ',g10.3,a)") & 453 & "(long.,lat.) = (",long_dip,lat_dip,") Deg" 454 endif 455 456 endif 457 !----------------------------------------------------------------------- 458 ! input weights 459 !----------------------------------------------------------------------- 460 461 dw8=0.0_dp ! read as DP, even in SP in FITS file 462 if (trim(infile_w8)/="") then 463 if (won == 1) then 464 ! ring based weights: treated as such in map2alm 465 n_rings = 2 * nsmax 466 n_rings_read = getsize_fits(infile_w8, nmaps=nmw8) 467 468 if (n_rings_read /= n_rings) then 469 print *," " 470 print *,"wrong ring weight file:"//trim(infile_w8) 471 call fatal_error(code) 472 endif 473 nmw8 = min(nmw8, n_pols) 474 475 PRINT *," "//code//"> Inputting "//trim(strw8file) 476 call input_map(infile_w8, dw8, n_rings, nmw8, fmissval=0.0_dp) 477 endif 478 479 if (won == 2) then 480 ! pixel based weights: turned into full sky map and treat as mask 481 allocate(pw8map (0:npixtot-1, 1:1), stat=status) 482 call assert_alloc(status,code,"pw8map") 483 484 PRINT *," "//code//"> Inputting and unfolding "//trim(strw8file) 485 486! n_pw8_read = getsize_fits(infile_w8) 487! allocate(pw8list(0:n_pw8_read-1,1:1), stat=status) 488! call assert_alloc(status,code,"pw8list") 489! call input_map(infile_w8, pw8list, n_pw8_read, 1) 490! call unfold_weightslist(nsmax, won, pw8list(:,1), pw8map(:,1)) 491! deallocate(pw8list) 492 493 call unfold_weightsfile(infile_w8, pw8map(:,1)) 494 495 pmask => pw8map(:,1:1) 496 !print*,pmask(1:12,1) 497 !print*,pmask(npixtot-11:npixtot,1) 498 endif 499 endif 500 w8ring_TQU = 1.d0 + dw8 501 if (won /= 2) allocate(pmask(0:1,1)) 502 503 !----------------------------------------------------------------------- 504 ! reorder maps to RING if necessary 505 !----------------------------------------------------------------------- 506 if (order_type == 2) then 507 print *," "//code//"> Convert Nest -> Ring " 508 call convert_nest2ring (nsmax, map_TQU) 509 endif 510 511 !----------------------------------------------------------------------- 512 ! precomputed plms 513 !----------------------------------------------------------------------- 514 515 if (n_plm/=0) then 516 if (plm_pol.eq.1) then 517 print*," "//code//"> Reading precomputed scalar P_lms " 518 else 519 print*," "//code//"> Reading precomputed tensor P_lms " 520 endif 521 !call read_dbintab(infile_plm,plm,n_plm,plm_pol,nullval,anynull=bad) ! replaced 2010-11-24 522 call read_bintab(infile_plm, plm, n_plm, plm_pol, nullval, anynull=bad) 523 if (bad) call fatal_error ("Missing points in P_lm-file!") 524 endif 525 526 !----------------------------------------------------------------------- 527 ! map to alm 528 !----------------------------------------------------------------------- 529 530 ALLOCATE(alm_TGC(1:n_pols, 0:nlmax, 0:nmmax),stat = status) 531 call assert_alloc(status,code,"alm_TGC") 532 533 PRINT *," "//code//"> Analyse map " 534 535 if (n_plm/=0) then 536 call map2alm_iterative(nsmax, nlmax, nmmax, iter_order, map_TQU, alm_TGC, & 537 & zbounds=zbounds, w8ring=w8ring_TQU, plm=plm, mask=pmask) 538 else 539 call map2alm_iterative(nsmax, nlmax, nmmax, iter_order, map_TQU, alm_TGC, & 540 & zbounds=zbounds, w8ring=w8ring_TQU, mask=pmask) 541 endif 542 543 !----------------------------------------------------------------------- 544 ! smoothing of the alm 545 !----------------------------------------------------------------------- 546 print *,' '//code//'> Alter alm ' 547 call alter_alm(nsmax, nlmax, nmmax, fwhm_arcmin, alm_TGC(1:1+2*polar,0:nlmax,0:nmmax), beam_file) 548 549 !----------------------------------------------------------------------- 550 ! a_lm to map 551 !----------------------------------------------------------------------- 552 553 print *,' '//code//'> Generating sky map(s) ' 554 555 select case (plm_pol+polar*10) 556 case(0) 557 call alm2map(nsmax,nlmax,nmmax,alm_TGC,map_TQU(:,1)) 558 case(1) 559 call alm2map(nsmax,nlmax,nmmax,alm_TGC,map_TQU(:,1),plm=plm(:,1)) 560 case(10) 561 call alm2map(nsmax,nlmax,nmmax,alm_TGC,map_TQU) 562 case(11) 563 call alm2map(nsmax,nlmax,nmmax,alm_TGC,map_TQU,plm=plm(:,1)) 564 case(13) 565 call alm2map(nsmax,nlmax,nmmax,alm_TGC,map_TQU,plm=plm) 566 end select 567 568 !----------------------------------------------------------------------- 569 ! convert back to NEST if necessary 570 !----------------------------------------------------------------------- 571 if (order_type == 2) then 572 print *," "//code//"> Convert Ring -> Nest " 573 call convert_ring2nest (nsmax, map_TQU) 574 endif 575 576 !----------------------------------------------------------------------- 577 ! generates header 578 !----------------------------------------------------------------------- 579 nlheader = SIZE(header) 580 fwhm_deg = fwhm_arcmin/60. 581 print *,' '//code//'> Writing smoothed map to FITS file ' 582 583 call write_minimal_header(header, 'map', & 584 order = order_type, nside = nsmax, coordsys = coordsys, & 585 creator = CODE, version = VERSION, nlmax = nlmax, & 586 beam_leg = trim(beam_file), fwhm_degree = fwhm_deg * 1.d0, & 587 polar = polarisation, units = units_map(1) ) 588 589 call add_card(header,'EXTNAME','SMOOTHED DATA', update = .true.) 590 call add_card(header) 591 call add_card(header,'COMMENT','************************* Input data ************************* ') 592 call add_card(header,'COMMENT','Input Map in '//TRIM(infile)) 593 call add_card(header,'COMMENT','*************************************************************** ') 594 call add_card(header,"NITERALM",iter_order, " Number of iterations for a_lms extraction") 595 call add_card(header,"LREGRESS",lowlreg, " Regression of low multipoles (0/1/2)") 596 call add_card(header,"QUADSCHM",won, " Quadrature scheme (0/1/2)") 597 call add_card(header) ! blank line 598 !----------------------------------------------------------------------- 599 ! write the map to FITS file 600 !----------------------------------------------------------------------- 601 !call output_map(map_TQU(0:npixtot-1,1:n_pols), header, outfile) !not on Sun 602 call write_bintab(map_TQU, npixtot, n_pols, header, nlheader, outfile) 603 604 !----------------------------------------------------------------------- 605 ! deallocate memory for arrays 606 !----------------------------------------------------------------------- 607 DEALLOCATE(map_TQU) 608 DEALLOCATE(w8ring_TQU) 609 DEALLOCATE(dw8) 610 611 !----------------------------------------------------------------------- 612 ! output and report card 613 !----------------------------------------------------------------------- 614 call wall_clock_time(time1) 615 call cpu_time(ptime1) 616 clock_time = time1 - time0 617 ptime = ptime1 - ptime0 618 619 write(*,9000) " " 620 write(*,9000) " Report Card for "//code//" run" 621 write(*,9000) " -----------------------------" 622 write(*,9000) " " 623 if (.not. polarisation) then 624 write(*,9000) " Temperature alone" 625 else 626 write(*,9000) " Temperature + Polarisation" 627 endif 628 write(*,9000) " " 629 630 write(*,9000) " Input map : "//TRIM(infile) 631 write(*,9006) " Multipole range : ",lowlreg ," <= l <= ", nlmax 632 write(*,9010) " Number of pixels : ", npixtot 633 write(*,9020) " Pixel size in arcmin : ", pix_size_arcmin 634 if (trim(beam_file) == '') then 635 write(*,9020) " Gauss. FWHM in arcmin : ", fwhm_arcmin 636 else 637 write(*,9000) " Beam file : "//trim(beam_file) 638 endif 639 write(*,9000) " Output map : "//TRIM(outfile) 640 if (won >0) write(*,9000) " Quadrature weight file : "//trim(infile_w8) 641 write(*,9030) " Clock and CPU time [s] : ", clock_time, ptime 642 643 !----------------------------------------------------------------------- 644 ! end of routine 645 !----------------------------------------------------------------------- 646 write(*,9000) " " 647 write(*,9000) " "//code//"> normal completion" 648 6499000 format(a) 6509005 format(a,i8) 6519006 format(a,i2,a,i8) 6529010 format(a,i16) 6539020 format(a,g20.5) 6549030 format(a,f11.2,f11.2) 655