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 fitstools 29 ! 30 ! module fitstools 31 ! toward version 2.1 32 ! 33 ! addition of read_fits_cut4 and write_fits_cut4 34 ! improvement of input_map (addition of header output) 35 ! improvement of read_asctab (more flexible on number of columns) 36 ! removal of ask_* routines 37 ! addition of get_clean_header 38 ! addition of optional keyword units to input_map, read_fits_cut4, read_bintab 39 ! 40 ! cleaning of the mapping between (l,m) <-> lm = l*l+l+m+1 41 ! to remove round-off errors 42 ! 43 ! v1.2_p01 : input_map, read_bintab : fix SUN compiler bug 44 ! getsize_fits : correct handling of ASCII table 45 ! v1.2_p02 getsize_fits : correct handling of files with more than 2^31 elements 46 ! read_bintab : correction to array argument in ftgcve 47 ! 2004-Feb: 48 ! modification of dump_alms to deal with large arrays on MacOsX 49 ! modification of npix in dump_alms 50 ! modification of read_conbintab to deal with large arrrays 51 ! 52 ! read_asctab : removal of nl_header from arguments list 53 ! 54 ! 2004-07-21 : overload (SP/DP) output_map and write_bintab 55 ! write_dbintab has been renamed write_plm 56 ! 2004-11-01 : added function getnumext_fits 57 ! added extno keyword to function getsize_fits 58 ! 2004-12-22 : overload (SP/DP) read_asctab, dump_alms 59 ! 2005-03-03 : added COORDSYS optional output to getsize_fits 60 ! 2006-04-04 : edited to avoid concatenation problem in TFORMs (write_plm and write_bintabh with Ifort 9) 61 ! 2007-02-16 : added EXTNO and POLARISATION keywords to Write_fits_cut4 62 ! 2007-04-04 : close file on exit of getnumext_fits 63 ! 2007-05-10 : increased maxdim (max number of columns per extension) from 20 to 40 64 ! 2007-10-01 : getsize_fits returns -1 in case of error 65 ! 2008-08-27 : in dump_alms and write_alms and write_*tab*: 66 ! do not write TTYPE# and TFORM# in excess of # of fields in the file 67 ! 2010-11-23: implemented support for large (>2^31-1 pixel) map IO (ie, Nside > 8192) 68 ! alm related IO is still limited to l<46340 69 ! 2013-12-13 : increased MAXDIM from 40 to MAXDIM_TOP 70 ! 2014-05-02: fixed problem with keywords having a long string value 71 ! 2016-08-16: debugged input_map on cut4 FITS files with multi-HDU polarized data 72 ! 2018-05-22: added unfold_weightsfile and unfold_weightslist 73 ! ------------------------------------------------------------- 74 ! 75 ! --------------------------- from include file (see fits_template.f90) 76 ! subroutine input_map 77 ! subroutine read_bintab 78 ! subroutine read_conbintab 79 ! subroutine write_bintab (write_bintab4 and write_bintab8) 80 ! subroutine write_asctab 81 ! subroutine dump_alms 82 ! subroutine write_alms 83 ! subroutine read_alms 84 ! subroutine read_bintod 85 ! subroutine write_bintabh 86 ! subroutine unfold_weights 87 ! ---------------------------------- 88 ! 89 ! subroutine read_fits_cut4 ? 90 ! subroutine write_fits_cut4 ? 91 ! subroutine fits2cl NA accepts extname addressing (eg myfile.fits[beam]) 92 ! subroutine read_asctab NA 93 ! 94 ! subroutine output_map (4/8) 95 ! subroutine write_dbintab : Obsolete 96 ! subroutine write_plm NA 97 ! 98 ! subroutine alms2fits NotYet 99 ! subroutine fits2alms NotYet 100 ! subroutine input_tod (8) 101 ! subroutine read_dbintab NA 102 ! 103 ! subroutine printerror NA 104 ! 105 ! subroutine read_par NA 106 ! function getnumext_fits NA 107 ! function getsize_fits (8) accepts extname addressing (eg myfile.fits[beam]) 108 ! function number_of_alms NotYet 109 ! 110 ! subroutine putrec NA 111 ! subroutine get_clean_header NA 112 ! 113 ! subroutine check_input_map 114 ! 115 use healpix_types 116 use misc_utils, only : fatal_error, assert, assert_alloc, strupcase, string 117 implicit none 118 119 real(kind=SP), private, parameter :: s_bad_value = HPX_SBADVAL 120 real(kind=DP), private, parameter :: d_bad_value = HPX_DBADVAL 121 integer(kind=I4B), private, parameter :: i_bad_value = -1637500000 122 integer(I4B) , private, parameter :: nchunk_max = 12000 123 integer(I4B), private, parameter :: MAXDIM_TOP = 199 ! < 999 124 125! interface read_fits_cut4 126! module procedure read_fits_cut4_s,read_fits_cut4_d 127! end interface 128 129! interface write_fits_cut4 130! module procedure write_fits_cut4_s,write_fits_cut4_d 131! end interface 132 133 134 interface input_map 135#ifdef NO64BITS 136 module procedure input_map8_s,input_map8_d 137#else 138 module procedure input_map8_s,input_map8_d, input_map4_s,input_map4_d 139#endif 140 end interface 141 142 interface output_map 143 module procedure output_map_s, output_map_d 144 end interface 145 146 interface write_plm 147#ifdef NO64BITS 148 module procedure write_plm8 149#else 150 module procedure write_plm8, write_plm4 151#endif 152 end interface 153 154 interface write_bintab 155#ifdef NO64BITS 156 module procedure write_bintab8_s, write_bintab8_d 157#else 158 module procedure write_bintab8_s, write_bintab8_d, write_bintab4_s, write_bintab4_d 159#endif 160 end interface 161 162 interface read_bintab 163#ifdef NO64BITS 164 module procedure read_bintab8_s, read_bintab8_d 165#else 166 module procedure read_bintab8_s, read_bintab8_d, read_bintab4_s, read_bintab4_d 167#endif 168 end interface 169 170 ! 171 interface input_tod 172 module procedure input_tod_s,input_tod_d 173 end interface 174 175 interface read_bintod 176 module procedure read_bintod_s,read_bintod_d 177 end interface 178 179 interface write_bintabh 180 module procedure write_bintabh_s,write_bintabh_d 181 end interface 182 ! 183 interface write_asctab 184 module procedure write_asctab_s, write_asctab_d 185 end interface 186 187 interface fits2cl 188 module procedure fits2cl_s, fits2cl_d 189 end interface 190 191 interface read_asctab 192 module procedure read_asctab_s, read_asctab_d, read_asctab_v12s, read_asctab_v12d 193 end interface 194 ! 195 interface fits2alms 196 module procedure fits2alms_s, fits2alms_d 197 end interface 198 199 interface alms2fits 200 module procedure alms2fits_s, alms2fits_d 201 end interface 202 203 interface read_conbintab 204 module procedure read_conbintab_s, read_conbintab_d 205 end interface 206 207 interface dump_alms 208 module procedure dump_alms_s, dump_alms_d 209 end interface 210 211 interface read_alms 212 module procedure read_alms_s, read_alms_d 213 end interface 214 215 interface write_alms 216 module procedure write_alms_s, write_alms_d 217 end interface 218 219 interface f90ftpcl_ 220 module procedure f90ftpcle, f90ftpcld 221 end interface 222 223 interface f90ftgcv_ 224 module procedure f90ftgcve, f90ftgcvd 225 end interface 226 227 interface f90ftgpv_ 228 module procedure f90ftgpve, f90ftgpvd 229 end interface 230 231 interface f90ftgky_ 232 module procedure f90ftgkye, f90ftgkyd 233 end interface 234 235 interface map_bad_pixels 236 module procedure map_bad_pixels_s, map_bad_pixels_d 237 end interface 238 239 interface unfold_weightsfile 240 module procedure unfold_weightsfile_s, unfold_weightsfile_d 241 end interface 242 243 interface unfold_weightslist 244 module procedure unfold_weightslist_s, unfold_weightslist_d 245 end interface 246 247 private 248 249 public :: read_fits_cut4, write_fits_cut4, & 250 & input_map, read_bintab, & 251 & output_map, write_bintab 252 public :: write_plm 253 public :: fits2cl, read_asctab, write_asctab 254 public :: read_dbintab, write_dbintab 255 public :: input_tod, write_bintabh 256 public :: dump_alms, alms2fits, fits2alms, write_alms, read_alms, read_conbintab 257 public :: printerror, read_par, getsize_fits, number_of_alms, getnumext_fits 258 public :: putrec 259 public :: check_input_map 260 public :: unfold_weightsfile, unfold_weightslist 261 262contains 263 264 !------------------------------------------------------------------------------- 265 ! generic interface F90FTPCL_ for FITSIO's FTPCLE and FTPCLD 266 ! writes data in ASCTAB or BINTAB 267 subroutine f90ftpcle(unit, colnum, frow, felem, np, data, status) 268 integer(I4B), intent(in) :: unit, colnum, frow, felem, np 269 integer(I4B), intent(out) :: status 270 real(SP), intent(in), dimension(0:) :: data 271 call ftpcle(unit, colnum, frow, felem, np, data, status) 272 return 273 end subroutine f90ftpcle 274 subroutine f90ftpcld(unit, colnum, frow, felem, np, data, status) 275 integer(I4B), intent(in) :: unit, colnum, frow, felem, np 276 integer(I4B), intent(out) :: status 277 real(DP), intent(in), dimension(0:) :: data 278 call ftpcld(unit, colnum, frow, felem, np, data, status) 279 return 280 end subroutine f90ftpcld 281 !------------------------------------------------------------------------------- 282 ! generic interface F90FTGCV_ for FITSIO's FTGCVE and FTGCVD 283 ! reads data from BINTAB 284 subroutine f90ftgcve(unit, colnum, frow, felem, np, nullval, data, anynull, status) 285 integer(I4B), intent(in) :: unit, colnum, frow, felem, np 286 integer(I4B), intent(out) :: status 287 logical(LGT), intent(out) :: anynull 288 real(SP), intent(out), dimension(0:) :: data 289 real(SP), intent(in) :: nullval 290 call ftgcve(unit, colnum, frow, felem, np, nullval, data, anynull, status) 291 return 292 end subroutine f90ftgcve 293 subroutine f90ftgcvd(unit, colnum, frow, felem, np, nullval, data, anynull, status) 294 integer(I4B), intent(in) :: unit, colnum, frow, felem, np 295 integer(I4B), intent(out) :: status 296 logical(LGT), intent(out) :: anynull 297 real(DP), intent(out), dimension(0:) :: data 298 real(DP), intent(in) :: nullval 299 call ftgcvd(unit, colnum, frow, felem, np, nullval, data, anynull, status) 300 return 301 end subroutine f90ftgcvd 302 !------------------------------------------------------------------------------- 303 ! generic interface F90FTGPV_ for FITSIO's FTGPVE and FTGPVD 304 ! reads data from IMAGE 305 subroutine f90ftgpve(unit, group, firstpix, np, nullval, data, anynull, status) 306 integer(I4B), intent(in) :: unit, group, firstpix, np 307 integer(I4B), intent(out) :: status 308 logical(LGT), intent(out) :: anynull 309 real(SP), intent(out), dimension(0:) :: data 310 real(SP), intent(in) :: nullval 311 call ftgpve(unit, group, firstpix, np, nullval, data, anynull, status) 312 return 313 end subroutine f90ftgpve 314 subroutine f90ftgpvd(unit, group, firstpix, np, nullval, data, anynull, status) 315 integer(I4B), intent(in) :: unit, group, firstpix, np 316 integer(I4B), intent(out) :: status 317 logical(LGT), intent(out) :: anynull 318 real(DP), intent(out), dimension(0:) :: data 319 real(DP), intent(in) :: nullval 320 call ftgpvd(unit, group, firstpix, np, nullval, data, anynull, status) 321 return 322 end subroutine f90ftgpvd 323 !------------------------------------------------------------------------------- 324 ! generic interface F90FTGKY_ for FITSIO's FTGKYE and FTGKYD 325 ! reads a keyword 326 subroutine f90ftgkye(unit, keyword, value, comment, status) 327 integer(I4B), intent(in) :: unit 328 character(len=*), intent(in) :: keyword 329 integer(I4B), intent(out) :: status 330 character(len=*), intent(out) :: comment 331 real(SP), intent(out) :: value 332 call ftgkye(unit, keyword, value, comment, status) 333 return 334 end subroutine f90ftgkye 335 subroutine f90ftgkyd(unit, keyword, value, comment, status) 336 integer(I4B), intent(in) :: unit 337 character(len=*), intent(in) :: keyword 338 integer(I4B), intent(out) :: status 339 character(len=*), intent(out) :: comment 340 real(DP), intent(out) :: value 341 call ftgkyd(unit, keyword, value, comment, status) 342 return 343 end subroutine f90ftgkyd 344 !------------------------------------------------------------------------------- 345 346 347 ! define routine with SP I/O 348#include "fits_s_inc.f90" 349 350 ! define routine with DP I/O 351#include "fits_d_inc.f90" 352 353 354 355 !======================================================================= 356 subroutine read_fits_cut4(filename, npixtot, pixel, signal, n_obs, serror, header, units, extno) 357 !======================================================================= 358 ! 359 ! routine to read FITS file with 4 columns : PIXEL, SIGNAL, N_OBS, SERROR 360 ! 361 ! read_fits_cut4(filename, npixtot, pixel, 362 ! [signal, n_obs, serror, header, units, extno]) 363 !======================================================================= 364 integer(I4b), parameter :: KMAP = SP 365 366 character(len=*), intent(in) :: filename 367 integer(I4B), intent(in) :: npixtot 368 integer(I4B), dimension(0:npixtot-1), intent(out) :: pixel 369 real(KMAP), dimension(0:npixtot-1), intent(out), optional :: signal 370 integer(I4B), dimension(0:npixtot-1), intent(out), optional :: n_obs 371 real(KMAP), dimension(0:npixtot-1), intent(out), optional :: serror 372 character(len=*), dimension(1:), intent(out), optional :: header 373 character(len=*), intent(out), optional :: units 374 integer(I4B), intent(in), optional :: extno 375 376 integer(I4B) :: obs_npix 377 378 integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension 379 integer(I4B) :: blocksize, datacode 380 integer(I4B) :: firstpix, frow, hdutype 381 integer(I4B) :: naxis, nfound, nmove, npix, nrows 382 integer(I4B) :: readwrite, status, tfields, unit, varidat, width 383 integer(I4B) :: repeat, repeat1, repeat2 384 logical(LGT) :: anynull, extend 385 character(len=20), dimension(MAXDIM) :: ttype, tform, tunit 386 character(len=20) :: extname 387 character(len=80) :: comment 388 real(KMAP) :: nullval 389 !======================================================================= 390 391 status=0 392 393 unit = 150 394 nfound = -1 395 anynull = .false. 396 397 readwrite=0 398 call ftopen(unit,filename,readwrite,blocksize,status) 399 if (status > 0) call printerror(status) 400 401 ! determines the presence of image 402 call ftgkyj(unit,'NAXIS', naxis, comment, status) 403 if (status > 0) call printerror(status) 404 if (naxis > 0) then ! there is an image 405 print*,'an image was found in the FITS file '//filename 406 print*,'... it is ignored.' 407 endif 408 409 ! determines the presence of an extension 410 call ftgkyl(unit,'EXTEND', extend, comment, status) 411 if (status > 0) then 412 print*,'extension expected and not found in FITS file '//filename 413 print*,'abort code' 414 call fatal_error 415 endif 416 417 nmove = +1 418 if (present(extno)) nmove = +1 + extno 419 call ftmrhd(unit, nmove, hdutype, status) 420 421 call assert (hdutype==2, 'this is not a binary table') 422 423 ! reads all the FITS related keywords 424 call ftghbn(unit, MAXDIM, & 425 & nrows, tfields, ttype, tform, tunit, extname, varidat, & 426 & status) 427 if (tfields < 4) then 428 print*,'Expected 4 columns in FITS file '//filename 429 print*,'found ',tfields 430 if (tfields < 2) call fatal_error 431 if (.not. (trim(ttype(1)) == 'PIXEL' & 432 & .or. trim(ttype(1)) == 'PIX' ) ) call fatal_error 433 endif 434 435 if (present(header)) then 436 header = "" 437 status = 0 438 call get_clean_header(unit, header, filename, status) 439 endif 440 441 if (present(units)) then 442 ! the second column contains the SIGNAL, for which we 443 ! already read the units from the header 444 units = adjustl(tunit(2)) 445 endif 446 447 ! finds the bad data value 448! if (KMAP == SP) CALL ftgkye(unit,'BAD_DATA',nullval,comment,status) 449! if (KMAP == DP) CALL ftgkyd(unit,'BAD_DATA',nullval,comment,status) 450 call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status) 451 if (status == 202) then ! bad_data not found 452 if (KMAP == SP) nullval = s_bad_value ! default value 453 if (KMAP == DP) nullval = d_bad_value ! default value 454 status = 0 455 endif 456 457 frow = 1 458 firstpix = 1 459 ! parse TFORM keyword to find out the length of the column vector 460 repeat1 = 1 461 repeat2 = 1 462 call ftbnfm(tform(1), datacode, repeat1, width, status) 463 if (tfields > 1) call ftbnfm(tform(2), datacode, repeat2, width, status) 464 repeat = max(repeat1, repeat2) 465 466 call ftgkyj(unit,'OBS_NPIX',obs_npix,comment,status) 467 if (status == 202) then ! obs_npix not found 468 obs_npix = nrows * repeat 469 status = 0 470 endif 471 472 npix = min(npixtot, obs_npix) 473 call ftgcvj(unit, 1_i4b, frow, firstpix, npix, i_bad_value, pixel(0), anynull, status) 474 if (present(signal)) then 475 call f90ftgcv_(unit, 2_i4b, frow, firstpix, npix, nullval, signal(0:npix-1), anynull, status) 476 endif 477 if (tfields >= 3 .and. present(n_obs)) & 478 & call ftgcvj(unit, 3_i4b, frow, firstpix, npix, i_bad_value, n_obs(0), anynull, status) 479 if (tfields >= 4 .and. present(serror)) then 480 call f90ftgcv_(unit, 4_i4b, frow, firstpix, npix, nullval, serror(0:npix-1), anynull, status) 481 endif 482 483 ! close the file 484 call ftclos(unit, status) 485 486 ! check for any error, and if so print out error messages 487 if (status > 0) call printerror(status) 488 489 return 490 end subroutine read_fits_cut4 491 492 !======================================================================= 493 subroutine write_fits_cut4(filename, obs_npix, pixel, signal, n_obs, serror, & 494 & header, coord, nside, order, units, extno, polarisation) 495 !======================================================================= 496 ! writes a fits file for cut sky data set with 4 columns: 497 ! PIXEL, SIGNAL, N_OBS, SERROR 498 ! 499 ! write_fits_cut4(filename, obs_npix, pixel, signal, n_obs, serror 500 ! [, header, coord, nside, order, units, extno, polarisation]) 501 !======================================================================= 502 use pix_tools, only: nside2npix 503 integer(I4b), parameter :: KMAP = SP 504 505 character(len=*), intent(in) :: filename 506 integer(I4B) :: obs_npix 507 integer(I4B), dimension(0:obs_npix-1), intent(in) :: pixel 508 real(KMAP), dimension(0:obs_npix-1), intent(in) :: signal 509 integer(I4B), dimension(0:obs_npix-1), intent(in) :: n_obs 510 real(KMAP), dimension(0:obs_npix-1), intent(in) :: serror 511 character(len=*), dimension(1:), intent(in), optional :: header 512 integer(I4B), intent(in), optional :: nside, order 513 character(len=*), intent(in), optional :: coord, units 514 integer(I4B), intent(in), optional :: extno, polarisation 515 516 character(len=*), parameter :: routine = 'write_fits_cut4' 517 ! --- healpix/fits related variables ---- 518 integer(I4B) :: ncol, grain 519 integer(I4B) :: npix_hd, nside_hd, nside_final, npix_final, i 520 integer(I4B) :: maxpix, minpix 521 character(len=1) :: char1, coord_usr 522! character(len=8) :: char8 523 character(len=20) :: units_usr 524 logical(LGT) :: done_nside, done_order, done_coord, done_polar, polar_flag 525 integer(I4B) :: nlheader, extno_i 526 527 ! --- cfitsio related variables ---- 528 integer(I4B) :: status, unit, blocksize, bitpix, naxis, naxes(1) 529 logical(LGT) :: simple, extend 530 character(LEN=80) :: svalue, comment 531 532 integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension 533 integer(I4B) :: nrows, tfields, varidat 534 integer(I4B) :: frow, felem, repeat, repeatg, hdutype 535 character(len=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname 536 character(len=20), dimension(1:3) :: extnames 537 character(len=80) :: card 538 character(len=4) :: srepeat, srepeatg 539 character(len=1) :: pform 540 character(len=filenamelen) sfilename 541 542 integer(I4B), save :: nside_old 543 !======================================================================= 544 if (KMAP == SP) pform='E' 545 if (KMAP == DP) pform='D' 546 547 extnames = (/ 'TEMPERATURE ' , & ! default in case of polarisation 548 & 'Q_POLARISATION' , & 549 & 'U_POLARISATION' /) 550 551 ncol = 4 552 grain = 1 553 status=0 554 unit = 149 555 blocksize=1 556 557 nlheader = 0 558 if (present(header)) nlheader = size(header) 559 units_usr = ' ' 560 if (present(units)) units_usr = units 561 extno_i = 0 562 if (present(extno)) extno_i = extno 563 polar_flag = .false. 564 if (present(polarisation)) then 565 polar_flag = (polarisation /= 0) 566 done_polar = .true. ! will ignore the header provided value of POLAR 567 endif 568 569 if (extno_i == 0) then 570 !************************************* 571 ! create the new empty FITS file 572 !************************************* 573 call ftinit(unit,filename,blocksize,status) 574 if (status > 0) call fatal_error("Error while creating file " & 575 & //trim(filename) & 576 & //". Check path and/or access rights.") 577 578 ! ----------------------------------------------------- 579 ! initialize parameters about the FITS image 580 simple=.true. 581 bitpix=32 ! integer*4 582 naxis=0 ! no image 583 naxes(1)=0 584 extend=.true. ! there is an extension 585 586 ! ---------------------- 587 ! primary header 588 ! ---------------------- 589 ! write the required header keywords 590 call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status) 591 592 ! write the current date 593 call ftpdat(unit,status) ! format ccyy-mm-dd 594 595 else 596 !********************************************* 597 ! reopen an existing file and go to the end 598 !********************************************* 599 ! remove the leading '!' (if any) when reopening the same file 600 sfilename = adjustl(filename) 601 if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen) 602 call ftopen(unit,sfilename,1_i4b,blocksize,status) 603 ! move to last extension written 604 call ftmahd(unit, 1_i4b+ extno_i, hdutype, status) 605 606 endif 607 608 609 ! ---------------------- 610 ! new extension 611 ! ---------------------- 612 call ftcrhd(unit, status) 613 614 ! writes required keywords 615 ! repeat = 1024 616 repeat = 1 617 nrows = (obs_npix + repeat - 1)/ repeat ! naxis1 618 if (obs_npix < repeat) repeat = 1 619 write(srepeat,'(i4)') repeat 620 srepeat = adjustl(srepeat) 621 622 repeatg = repeat * grain 623 write(srepeatg,'(i4)') repeatg 624 srepeatg = adjustl(srepeatg) 625 626 tfields = ncol 627 ttype(1:4) = (/ 'PIXEL ','SIGNAL','N_OBS ','SERROR' /) 628 tform(1) = trim(srepeat)//'J' 629 tform(2) = trim(srepeatg)//pform 630 tform(3) = trim(srepeatg)//'J' 631 tform(4) = trim(srepeatg)//pform 632 633 tunit = ' ' ! optional, will not appear 634 tunit(2) = units_usr 635 tunit(4) = units_usr 636 extname = 'SKY_OBSERVATION' ! default, will be overide by user provided one if any 637 if (polar_flag) extname = extnames(1+extno_i) 638 varidat = 0 639 call ftphbn(unit, nrows, tfields, ttype, tform, tunit, & 640 & extname, varidat, status) 641 642 call ftpcom(unit,'------------------------------------------',status) 643 call ftpcom(unit,' Pixelisation Specific Keywords ',status) 644 call ftpcom(unit,'------------------------------------------',status) 645 call ftpkys(unit,'PIXTYPE','HEALPIX ',' HEALPIX Pixelisation',status) 646 call ftpkyu(unit,'NSIDE', ' ',status) ! place holder, will be updated later on 647 call ftpkyu(unit,'ORDERING',' ',status) 648 call ftpkys(unit,'COORDSYS',' ',' ',status) 649 call ftpcom(unit,' G = Galactic, E = ecliptic, C = celestial = equatorial', status) 650 call ftpkyl(unit,'POLAR',polar_flag,'Polarisation included in file (T/F)',status) 651 call ftpcom(unit,'------------------------------------------',status) 652 call ftpcom(unit,' Data Specific Keywords ',status) 653 call ftpcom(unit,'------------------------------------------',status) 654 call ftpkys(unit,'INDXSCHM','EXPLICIT',' Indexing : IMPLICIT or EXPLICIT', status) 655 call ftpkyj(unit,'GRAIN', grain, ' Grain of pixel indexing',status) 656 call ftpcom(unit,'GRAIN=0 : no indexing of pixel data (IMPLICIT) ',status) 657 call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)',status) 658 call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status) 659 call ftpkys(unit,'OBJECT','PARTIAL ',' Sky coverage represented by data',status) 660 call ftpkyj(unit,'OBS_NPIX',obs_npix, ' Number of pixels observed and recorded',status) 661 662 663 ! add required Healpix keywords (NSIDE, ORDER) if provided by user 664 done_order = .false. 665 if (present(order)) then 666 ! char8 = order 667 ! call ftupch(char8) 668 ! if (char8(1:4) == 'RING') then 669 if (order == 1) then 670 call ftukys(unit, 'ORDERING','RING',' Pixel ordering scheme, either RING or NESTED',status) 671 done_order = .true. 672 ! elseif (char8(1:4) == 'NEST') then 673 elseif (order == 2) then 674 call ftukys(unit, 'ORDERING','NESTED',' Pixel ordering scheme, either RING or NESTED ',status) 675 done_order = .true. 676 else 677 print*,'Invalid ORDER given : ',order, ' instead of 1 (RING) or 2 (NESTED)' 678 endif 679 endif 680 681 done_nside = .false. 682 if (present(nside)) then 683 if (nside2npix(nside) > 0) then ! valid nside 684 call ftukyj(unit,'NSIDE',nside,' Resolution parameter for HEALPIX', status) 685 done_nside = .true. 686 nside_final = nside 687 else 688 print*,'Invalid NSIDE given : ',nside 689 endif 690 endif 691 692 ! add non required Healpix keyword (COORD) 693 done_coord = .false. 694 if (present(coord)) then 695 coord_usr = adjustl(coord) 696 char1 = coord_usr(1:1) 697 call ftupch(char1) ! uppercase 698 if (char1 == 'C' .or. char1 == 'Q') then 699 coord_usr = 'C' 700 done_coord = .true. 701 elseif (char1 == 'G') then 702 coord_usr = 'G' 703 done_coord = .true. 704 elseif (char1 == 'E' ) then 705 coord_usr='E' 706 done_coord = .true. 707 else 708 print*,'Unrecognised COORD given : ',coord,' instead of C, G, or E' 709 print*,'Proceed at you own risks ' 710 coord_usr = char1 711 done_coord = .true. 712 endif 713 if (done_coord) then 714 call ftukys(unit, 'COORDSYS',coord_usr,' Pixelisation coordinate system',status) 715 endif 716 endif 717 718 719 ! write the user provided header literally, except for PIXTYPE, TFORM*, TTYPE*, TUNIT* and INDXSCHM 720 ! copy NSIDE, ORDERING and COORDSYS and POLAR if they are valid and not already given 721 do i=1,nlheader 722 card = header(i) 723 if (card(1:5) == 'TTYPE' .or. card(1:5) == 'TFORM' .or. card(1:7) == 'PIXTYPE') then 724 continue ! don't keep them 725 else if (card(1:8) == 'INDXSCHM') then 726 continue 727 else if (card(1:5) == 'TUNIT') then 728 if ((card(6:6) == '2' .or. card(6:6) == '4') .and. trim(units_usr) == '') then 729 call ftucrd(unit,'TUNIT'//card(6:6),card, status) !update TUNIT2 and TUNIT4 730 endif 731 else if (card(1:5) == 'NSIDE') then 732 call ftpsvc(card, svalue, comment, status) 733 read(svalue,*) nside_hd 734 npix_hd = nside2npix(nside_hd) 735 if (.not. done_nside .and. npix_hd > 0) then 736 call ftucrd(unit,'NSIDE',card, status) !update NSIDE 737 done_nside = .true. 738 nside_final = nside_hd 739 endif 740 else if (card(1:8) == 'ORDERING') then 741 call ftpsvc(card, svalue, comment, status) 742 svalue = adjustl(svalue) 743 svalue = svalue(2:8) ! remove leading quote 744 call ftupch(svalue) 745 if (.not. done_order .and. (svalue(1:4)=='RING' .or. svalue(1:6) == 'NESTED')) then 746 call ftucrd(unit,'ORDERING',card, status) !update ORDERING 747 done_order = .true. 748 endif 749 else if (card(1:8) == 'COORDSYS') then 750 if (.not. done_coord) call putrec(unit,card, status) 751 done_coord = .true. 752 else if (card(1:5) == 'POLAR') then 753 if (.not. done_polar) then 754 call ftucrd(unit,'POLAR', card, status) ! update POLAR 755 done_polar = .true. 756 endif 757 else if (card(1:7) == 'EXTNAME') then 758 call ftucrd(unit, 'EXTNAME', card, status) 759 else if (card(1:1) /= ' ') then ! if none of the above, copy to FITS file 760 call putrec(unit, card, status) 761 endif 76210 continue 763 enddo 764 765 766 ! test that required keywords have been provided in some way 767 if (.not. done_nside) then 768 print*,routine//'> NSIDE is a Healpix required keyword, ' 769 print*,routine//'> it was NOT provided either as routine argument or in the input header' 770 print*,routine//'> abort execution, file not written' 771 call fatal_error 772 endif 773 if (.not. done_order) then 774 print*,routine//'> ORDER is a Healpix required keyword, ' 775 print*,routine//'> it was NOT provided either as routine argument or in the input header' 776 print*,routine//'> abort execution, file not written' 777 call fatal_error 778 endif 779 if ((.not. done_polar) .and. extno_i >=2) then 780 print*,routine//'> Warning: POLAR keyword not set while 3 extensions have been written' 781 endif 782 783 ! check that NSIDE is the same for all extensions 784 if (extno_i == 0) then 785 nside_old = nside_final 786 else 787 if (nside_final /= nside_old) then 788 print*,routine//'> Inconsistent NSIDE: ',nside_final, nside_old 789 print*,routine//'> Should use same NSIDE for all extensions' 790 endif 791 endif 792 793 ! check validity of PIXEL 794 npix_final = nside2npix(nside_final) 795 minpix = minval(pixel(0:obs_npix-1)) 796 maxpix = maxval(pixel(0:obs_npix-1)) 797 if (minpix < 0 .or. maxpix > npix_final-1) then 798 print*,routine//'> Actual pixel range = ',minpix,maxpix 799 print*,routine//'> expected range (for Nside =',nside_final,') : 0, ',npix_final-1 800 print*,routine//'> ABORT execution, file not written.' 801 call fatal_error 802 endif 803 if (obs_npix > npix_final) then 804 print*,routine//'> The actual number of pixels ',obs_npix 805 print*,routine//'> is larger than the number of pixels over the whole sky : ',npix_final 806 print*,routine//'> for Nside = ',nside_final 807 print*,routine//'> ABORT execution, file not written.' 808 call fatal_error 809 endif 810 811 ! write the extension one column by one column 812 frow = 1 ! starting position (row) 813 felem = 1 ! starting position (element) 814 call ftpclj(unit, 1_i4b, frow, felem, obs_npix, pixel , status) 815 call f90ftpcl_(unit, 2_i4b, frow, felem, obs_npix, signal, status) 816 817 if (tfields >= 3) call ftpclj(unit, 3_i4b, frow, felem, obs_npix, n_obs , status) 818 if (tfields >= 4) then 819 call f90ftpcl_(unit, 4_i4b, frow, felem, obs_npix, serror, status) 820 endif 821 822 ! ---------------------- 823 ! close and exit 824 ! ---------------------- 825 826 ! close the file and free the unit number 827 call ftclos(unit, status) 828 829 ! check for any error, and if so print out error messages 830 if (status > 0) call printerror(status) 831 832 return 833 end subroutine write_fits_cut4 834 835 !======================================================================= 836 ! FITS2CL(filename, clin, lmax, ncl, header, units, fmissval) 837 ! Read C_l from a FITS file 838 ! Aug 2000 : modification by EH of the number of columns actually read 839 ! 840 ! This routine is used for reading input power spectra for synfast 841 ! 842 ! Dec 2004: overloading for single and double precision output array (clin) 843 ! Feb 2013: added fmissval (value to be given to NaN or Inf C(l)) 844 ! if absent, those C(l) are left unchanged 845 ! 846 !======================================================================= 847 subroutine fits2cl_s(filename, clin, lmax, ncl, header, units, fmissval) 848 !======================================================================= 849 ! single precision 850 !======================================================================= 851 integer(I4b), parameter :: KCL = SP 852 ! 853 CHARACTER(LEN=*), INTENT(IN) :: filename 854 INTEGER(I4B), INTENT(IN) :: lmax, ncl 855 REAL(KCL), DIMENSION(0:lmax,1:ncl), INTENT(OUT) :: clin 856 CHARACTER(LEN=*), DIMENSION(1:), INTENT(OUT) :: header 857 CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units 858 real(KCL), optional, intent(IN):: fmissval 859 860 real(DP), dimension(:,:), allocatable :: cl_dp 861 real(DP) :: fmvd 862 863 ! since the arrays involved are small 864 ! read in double precision, and then cast in single 865 allocate(cl_dp(0:lmax, 1:ncl)) 866 if (present(fmissval)) then 867 fmvd = fmissval * 1.0_DP 868 call fits2cl_d(filename, cl_dp, lmax, ncl, header, units, fmvd) 869 else 870 call fits2cl_d(filename, cl_dp, lmax, ncl, header, units) 871 endif 872 clin(0:lmax, 1:ncl) = cl_dp(0:lmax, 1:ncl) 873 874 return 875 end subroutine fits2cl_s 876 877 subroutine fits2cl_d(filename, clin, lmax, ncl, header, units, fmissval) 878 !======================================================================= 879 ! double precision 880 !======================================================================= 881 integer(I4b), parameter :: KCL = DP 882 ! 883 CHARACTER(LEN=*), INTENT(IN) :: filename 884 INTEGER(I4B), INTENT(IN) :: lmax, ncl 885 REAL(KCL), DIMENSION(0:lmax, 1:ncl), INTENT(OUT) :: clin 886 CHARACTER(LEN=*), DIMENSION(1:), INTENT(OUT) :: header 887 CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units 888 real(KCL), optional, intent(IN):: fmissval 889 890 INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),ncl_file, naxis 891 INTEGER(I4B) :: firstpix,lmax_file,lmax_min,nelems,datacode,repeat,width 892 CHARACTER(LEN=80) :: comment, pdmstr 893 LOGICAL :: extend 894 INTEGER(I4B) :: nmove, hdutype, hdunum 895 INTEGER(I4B) :: column, frow 896 REAL(KCL), DIMENSION(:), ALLOCATABLE :: clin_file 897 898 INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension 899 INTEGER(I4B) :: rowlen, nrows, varidat 900 INTEGER(I4B), dimension(1:MAXDIM) :: tbcol 901 CHARACTER(LEN=20), dimension(1:MAXDIM) :: ttype, tform, tunit 902 CHARACTER(LEN=20) :: extname 903 LOGICAL :: anynull 904 REAL(KCL) :: nullval 905 logical :: planck_format 906 integer(i8b) :: nbads 907 908 909 !----------------------------------------------------------------------- 910 status=0 911 anynull = .false. 912 ttype='' 913 tform='' 914 tunit='' 915 extname='' 916 comment='' 917 918 unit = 110 919 naxes(1) = 1 920 naxes(2) = 1 921 922 readwrite=0 923 call ftnopn(unit,filename,readwrite, status) !open primary HDU or specified HDU 924 if (status > 0) call printerror(status) 925 ! ----------------------------------------- 926 call ftghdn(unit, hdunum) 927 if (hdunum == 1) then ! in primary HDU: move to next HDU 928 ! determines the presence of image 929 call ftgkyj(unit,'NAXIS', naxis, comment, status) 930 if (status > 0) call printerror(status) 931 932 ! determines the presence of an extension 933 call ftgkyl(unit,'EXTEND', extend, comment, status) 934 if (status > 0) call printerror(status) 935 936 nmove = +1 937 call ftmrhd(unit, nmove, hdutype, status) 938 else ! already in non primary HDU 939 extend = .true. 940 call ftghdt(unit, hdutype, status) 941 endif 942 943 if (extend) then ! there is an extension 944 945 call assert ((hdutype==1).or.(hdutype==2), 'this is not a table') 946 947 ! reads keywords related to table layout 948 if (hdutype==1) then ! ASCII table 949 call ftghtb(unit, MAXDIM, rowlen, & 950 & nrows, ncl_file, ttype, tbcol, & 951 & tform, tunit, extname, status) 952 repeat = 1 953 else ! binary table 954 call ftghbn(unit, MAXDIM, & 955 & nrows, ncl_file, ttype, & 956 & tform, tunit,extname, varidat, status) 957 call ftbnfm(tform(1), datacode, repeat, width, status) 958 endif 959 960 status = 0 961 header = "" 962 call get_clean_header(unit, header, filename, status) 963 964 if (present(units)) then 965 do column = 1, min(ncl, ncl_file) 966 units(column) = adjustl(tunit(column)) 967 enddo 968 endif 969 970 ! reads the columns 971 column = 1 972 frow = 1 973 firstpix = 1 974 lmax_file = nrows*repeat - 1 975 lmax_min = MIN(lmax,lmax_file) 976 nelems = lmax_min + 1 977 nullval = 0.0_KCL ! CFITSIO will leave bad pixels unchanged, and so will this routine 978 if (present(fmissval)) then 979 if (fmissval == 0.0_KCL) then 980 nullval = HPX_DBADVAL ! sentinel value to be given to bad pixels by CFITSIO, will later be mapped to user defined 0 981 else 982 nullval = fmissval ! user defined value to be given to bad pixels by CFITSIO, and by this routine 983 endif 984 endif 985 986! check for the special Planck format (i.e. one additional column) 987 planck_format=.true. 988 call ftgkys(unit,'PDMTYPE',pdmstr,comment,status) 989 if (status==202) then 990 planck_format=.false. 991 status=0 992 endif 993 allocate(clin_file(0:lmax_min),stat=status) 994 clin = 0.0_KCL ! modification by EH 995 if (planck_format) then 996 do column = 1, MIN(ncl, ncl_file-1) ! modification by EH 997 clin_file(:) = 0.0_KCL 998 call ftgcvd(unit, column+1_i4b, frow, firstpix, nelems, nullval, & 999 & clin_file(0), anynull, status) 1000 clin(0:lmax_min, column) = clin_file(0:lmax_min) 1001 enddo 1002 else 1003 do column = 1, MIN(ncl, ncl_file) ! modification by EH 1004 clin_file(:) = 0.0_KCL 1005 call ftgcvd(unit, column, frow, firstpix, nelems, nullval, & 1006 & clin_file(0), anynull, status) 1007 clin(0:lmax_min, column) = clin_file(0:lmax_min) 1008 enddo 1009 endif 1010 deallocate(clin_file) 1011 if (present(fmissval)) then 1012 if (nullval /= fmissval) call map_bad_pixels(clin, nullval, fmissval, nbads) 1013 endif 1014 else ! no image no extension, you are dead, man 1015 call fatal_error(' No image, no extension in '//trim(filename)) 1016 endif 1017 1018 ! close the file 1019 call ftclos(unit, status) 1020 1021 ! check for any error, and if so print out error messages 1022 if (status > 0) call printerror(status) 1023 1024! CHARACTER(LEN=*), INTENT(IN) :: filename 1025! INTEGER(I4B), INTENT(IN) :: lmax, ncl 1026! REAL(KCL), DIMENSION(0:lmax, 1:ncl), INTENT(OUT) :: clin 1027! CHARACTER(LEN=*), DIMENSION(1:), INTENT(OUT) :: header 1028! CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units 1029 1030! INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),ncl_file, naxis 1031! INTEGER(I4B) :: firstpix,lmax_file,lmax_min,nelems 1032! CHARACTER(LEN=80) :: comment, pdmstr 1033! LOGICAL :: extend 1034! INTEGER(I4B) :: nmove, hdutype 1035! INTEGER(I4B) :: column, frow 1036! REAL(KCL), DIMENSION(:), ALLOCATABLE :: clin_file 1037 1038! INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension 1039! INTEGER(I4B) :: rowlen, nrows, varidat 1040! INTEGER(I4B), dimension(1:MAXDIM) :: tbcol 1041! CHARACTER(LEN=20), dimension(1:MAXDIM) :: ttype, tform, tunit 1042! CHARACTER(LEN=20) :: extname 1043! LOGICAL :: anynull 1044! REAL(KCL) :: nullval 1045! logical :: planck_format 1046 1047 1048! !----------------------------------------------------------------------- 1049! status=0 1050! anynull = .false. 1051! ttype='' 1052! tform='' 1053! tunit='' 1054! extname='' 1055! comment='' 1056 1057! unit = 110 1058! naxes(1) = 1 1059! naxes(2) = 1 1060 1061! readwrite=0 1062! call ftopen(unit,filename,readwrite,blocksize,status) 1063! if (status > 0) call printerror(status) 1064! ! ----------------------------------------- 1065 1066! ! determines the presence of image 1067! call ftgkyj(unit,'NAXIS', naxis, comment, status) 1068! if (status > 0) call printerror(status) 1069 1070! ! determines the presence of an extension 1071! call ftgkyl(unit,'EXTEND', extend, comment, status) 1072! if (status > 0) call printerror(status) 1073 1074! if (extend) then ! there is an extension 1075! nmove = +1 1076! call ftmrhd(unit, nmove, hdutype, status) 1077 1078! call assert ((hdutype==1).or.(hdutype==2), 'this is not a table') 1079 1080! ! reads keywords related to table layout 1081! if (hdutype==1) then 1082! call ftghtb(unit, MAXDIM, & 1083! & rowlen, nrows, ncl_file, ttype, tbcol, tform, tunit, & 1084! & extname, status) 1085! else 1086! call ftghbn(unit, MAXDIM, & 1087! & nrows, ncl_file, ttype, tform, tunit,extname, & 1088! & varidat, status) 1089! endif 1090 1091! status = 0 1092! header = "" 1093! call get_clean_header(unit, header, filename, status) 1094 1095! if (present(units)) then 1096! do column = 1, min(ncl, ncl_file) 1097! units(column) = adjustl(tunit(column)) 1098! enddo 1099! endif 1100 1101! ! reads the columns 1102! column = 1 1103! frow = 1 1104! firstpix = 1 1105! lmax_file = nrows - 1 1106! lmax_min = MIN(lmax,lmax_file) 1107! nullval = 0.0_KCL 1108! nelems = lmax_min + 1 1109 1110! ! check for the special Planck format (i.e. one additional column) 1111! planck_format=.true. 1112! call ftgkys(unit,'PDMTYPE',pdmstr,comment,status) 1113! if (status==202) then 1114! planck_format=.false. 1115! status=0 1116! endif 1117! allocate(clin_file(0:lmax_min),stat=status) 1118! clin = 0.0_KCL ! modification by EH 1119! if (planck_format) then 1120! do column = 1, MIN(ncl, ncl_file-1) ! modification by EH 1121! clin_file(:) = 0.0_KCL 1122! call ftgcvd(unit, column+1_i4b, frow, firstpix, nelems, nullval, & 1123! & clin_file(0), anynull, status) 1124! clin(0:lmax_min, column) = clin_file(0:lmax_min) 1125! enddo 1126! else 1127! do column = 1, MIN(ncl, ncl_file) ! modification by EH 1128! clin_file(:) = 0.0_KCL 1129! call ftgcvd(unit, column, frow, firstpix, nelems, nullval, & 1130! & clin_file(0), anynull, status) 1131! clin(0:lmax_min, column) = clin_file(0:lmax_min) 1132! enddo 1133! endif 1134! deallocate(clin_file) 1135! else ! no image no extension, you are dead, man 1136! call fatal_error(' No image, no extension') 1137! endif 1138 1139! ! close the file 1140! call ftclos(unit, status) 1141 1142! ! check for any error, and if so print out error messages 1143! if (status > 0) call printerror(status) 1144 1145 1146 1147 return 1148 end subroutine fits2cl_d 1149 1150 !======================================================================= 1151 ! READ_ASC_TAB(filename, clin, lmax, ncl, header, units) 1152 ! superseded by fits2cl 1153 !======================================================================= 1154 subroutine read_asctab_v12s(filename, clin, lmax, ncl, header, nlheader, units) 1155 !======================================================================= 1156 ! single precision, legacy code 1157 !======================================================================= 1158 integer(I4b), parameter :: KCL = SP 1159 ! 1160 CHARACTER(LEN=*), INTENT(IN) :: filename 1161 INTEGER(I4B), INTENT(IN) :: lmax, ncl,nlheader 1162 REAL(KCL), DIMENSION(0:lmax,1:ncl), INTENT(OUT) :: clin 1163 CHARACTER(LEN=*), DIMENSION(1:), INTENT(OUT) :: header 1164 CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units 1165 1166 print*,"-------------------------------------------------------------" 1167 print*,"WARNING : the routine read_asctab is now obsolete" 1168 print*," Use " 1169 print*," call fits2cl(filename, clin, lmax, ncl, header, units)" 1170 print*," instead (same module)" 1171 print*,"-------------------------------------------------------------" 1172 1173 call fits2cl(filename, clin, lmax, ncl, header, units) 1174 return 1175 end subroutine read_asctab_v12s 1176 subroutine read_asctab_v12d(filename, clin, lmax, ncl, header, nlheader, units) 1177 !======================================================================= 1178 ! double precision, legacy code 1179 !======================================================================= 1180 integer(I4b), parameter :: KCL = DP 1181 ! 1182 CHARACTER(LEN=*), INTENT(IN) :: filename 1183 INTEGER(I4B), INTENT(IN) :: lmax, ncl,nlheader 1184 REAL(KCL), DIMENSION(0:lmax,1:ncl), INTENT(OUT) :: clin 1185 CHARACTER(LEN=*), DIMENSION(1:), INTENT(OUT) :: header 1186 CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units 1187 1188 print*,"-------------------------------------------------------------" 1189 print*,"WARNING : the routine read_asctab is now obsolete" 1190 print*," Use " 1191 print*," call fits2cl(filename, clin, lmax, ncl, header, units)" 1192 print*," instead (same module)" 1193 print*,"-------------------------------------------------------------" 1194 1195 call fits2cl(filename, clin, lmax, ncl, header, units) 1196 return 1197 end subroutine read_asctab_v12d 1198 !======================================================================= 1199 subroutine read_asctab_s(filename, clin, lmax, ncl, header, units) 1200 !======================================================================= 1201 ! single precision 1202 !======================================================================= 1203 integer(I4b), parameter :: KCL = SP 1204 ! 1205 CHARACTER(LEN=*), INTENT(IN) :: filename 1206 INTEGER(I4B), INTENT(IN) :: lmax, ncl 1207 REAL(KCL), DIMENSION(0:lmax,1:ncl), INTENT(OUT) :: clin 1208 CHARACTER(LEN=*), DIMENSION(1:), INTENT(OUT) :: header 1209 CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units 1210 1211 print*,"-------------------------------------------------------------" 1212 print*,"WARNING : the routine read_asctab is now obsolete" 1213 print*," Use " 1214 print*," call fits2cl(filename, clin, lmax, ncl, header, units)" 1215 print*," instead (same module)" 1216 print*,"-------------------------------------------------------------" 1217 1218 call fits2cl(filename, clin, lmax, ncl, header, units) 1219 return 1220 end subroutine read_asctab_s 1221 1222 subroutine read_asctab_d(filename, clin, lmax, ncl, header, units) 1223 !======================================================================= 1224 ! double precision 1225 !======================================================================= 1226 integer(I4b), parameter :: KCL = DP 1227 ! 1228 CHARACTER(LEN=*), INTENT(IN) :: filename 1229 INTEGER(I4B), INTENT(IN) :: lmax, ncl 1230 REAL(KCL), DIMENSION(0:lmax, 1:ncl), INTENT(OUT) :: clin 1231 CHARACTER(LEN=*), DIMENSION(1:), INTENT(OUT) :: header 1232 CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units 1233 1234 print*,"-------------------------------------------------------------" 1235 print*,"WARNING : the routine read_asctab is now obsolete" 1236 print*," Use " 1237 print*," call fits2cl(filename, clin, lmax, ncl, header, units)" 1238 print*," instead (same module)" 1239 print*,"-------------------------------------------------------------" 1240 1241 call fits2cl(filename, clin, lmax, ncl, header, units) 1242 return 1243 end subroutine read_asctab_d 1244 1245 !======================================================================= 1246 subroutine output_map_s(map, header, outfile, extno) 1247 !======================================================================= 1248 ! this is only a wrapper to write_bintab for SP maps 1249 ! this routine needs an explicit interface 1250 !======================================================================= 1251 use long_intrinsic, only: long_size 1252 1253 REAL(SP), INTENT(IN), DIMENSION(0:,1:), target :: map 1254 CHARACTER(LEN=*), INTENT(IN), DIMENSION(1:) :: header 1255 CHARACTER(LEN=*), INTENT(IN) :: outfile 1256 INTEGER(I4B) , INTENT(IN) , optional :: extno 1257 1258 INTEGER(I8B) :: npix 1259 INTEGER(I4B) :: nlheader, nmap, extno_i 1260 real(SP), pointer, dimension(:,:) :: pmap 1261 !----------------------------------------------------------------------- 1262 npix = long_size(map,1) 1263 nmap = long_size(map,2) 1264 nlheader = SIZE(header) 1265 1266 extno_i = 0 1267 if (present(extno)) extno_i = extno 1268 pmap => map(0:npix-1,1:nmap) 1269 call write_bintab(pmap, npix, nmap, & 1270 & header(1:nlheader), nlheader, outfile, extno=extno_i) 1271 return 1272 end subroutine output_map_s 1273 !======================================================================= 1274 subroutine output_map_d(map, header, outfile, extno) 1275 !======================================================================= 1276 ! this is only a wrapper to write_bintab for DP maps 1277 ! this routine needs an explicit interface 1278 !======================================================================= 1279 use long_intrinsic, only: long_size 1280 1281 REAL(DP), INTENT(IN), DIMENSION(0:,1:), target :: map 1282 CHARACTER(LEN=*), INTENT(IN), DIMENSION(1:) :: header 1283 CHARACTER(LEN=*), INTENT(IN) :: outfile 1284 INTEGER(I4B) , INTENT(IN) , optional :: extno 1285 1286 INTEGER(I8B) :: npix 1287 INTEGER(I4B) :: nlheader, nmap, extno_i 1288 real(DP), pointer, dimension(:,:) :: pmap 1289 !----------------------------------------------------------------------- 1290 npix = long_size(map,1) 1291 nmap = long_size(map,2) 1292 nlheader = SIZE(header) 1293 1294 extno_i = 0 1295 if (present(extno)) extno_i = extno 1296 pmap => map(0:npix-1,1:nmap) 1297 call write_bintab(pmap, npix, nmap, & 1298 & header(1:nlheader), nlheader, outfile, extno=extno_i) 1299 return 1300 end subroutine output_map_d 1301 1302 1303 1304 !======================================================================= 1305 subroutine write_dbintab(plm, nplm, nhar, header, nlheader, filename, nsmax, nlmax) 1306 !======================================================================= 1307 ! now obsolete routine 1308 !======================================================================= 1309 INTEGER(I4B), INTENT(IN) :: nplm, nhar, nlheader, nsmax, nlmax 1310 REAL(DP), INTENT(IN), DIMENSION(0:nplm-1,1:nhar) :: plm 1311 CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header 1312 CHARACTER(LEN=*), INTENT(IN) :: filename 1313 !======================================================================= 1314 print*,'WRITE_DBINTAB is obsolete.' 1315 print*,' ' 1316 print*,'To write a Healpix map into a FITS file' 1317 print*,'use WRITE_BINTAB or OUTPUT_MAP, which both accept ' 1318 print*,'Single and Double Precision arguments' 1319 print*,' ' 1320 print*,'To write Plm polynoms into a FITS file,' 1321 print*,'use WRITE_PLM (same syntax)' 1322 call write_plm(plm, nplm, nhar, header, nlheader, filename, nsmax, nlmax) 1323 1324 return 1325 end subroutine write_dbintab 1326 !======================================================================= 1327 ! WRITE_PLM 1328 ! Create a FITS file containing a binary table extension with 1329 ! the temperature map in the first column 1330 ! written by EH from writeimage and writebintable 1331 ! (fitsio cookbook package) 1332 ! 1333 ! slightly modified to deal with vector column 1334 ! in binary table EH/IAP/Jan-98 1335 ! 1336 ! simplified the calling sequence, the header sould be filled in 1337 ! before calling the routine 1338 ! 1339 ! Changed to write double precision plms for const.real. FKH/Apr-99 1340 ! 1341 !======================================================================= 1342#ifndef NO64BITS 1343 subroutine write_plm4(plm, nplm, nhar, header, nlheader, filename, nsmax, nlmax) 1344 !======================================================================= 1345 INTEGER(I4B), INTENT(IN) :: nplm 1346 INTEGER(I4B), INTENT(IN) :: nhar, nlheader, nsmax, nlmax 1347 REAL(DP), INTENT(IN), DIMENSION(0:nplm-1,1:nhar) :: plm 1348 CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header 1349 CHARACTER(LEN=*), INTENT(IN) :: filename 1350 integer(i8b) :: nplm8 1351 nplm8 = nplm 1352 call write_plm8(plm, nplm8, nhar, header, nlheader, filename, nsmax, nlmax) 1353 return 1354 end subroutine write_plm4 1355 !======================================================================= 1356#endif 1357 subroutine write_plm8(plm, nplm, nhar, header, nlheader, filename, nsmax, nlmax) 1358 !======================================================================= 1359 1360 INTEGER(I8B), INTENT(IN) :: nplm 1361 INTEGER(I4B), INTENT(IN) :: nhar, nlheader, nsmax, nlmax 1362 REAL(DP), INTENT(IN), DIMENSION(0:nplm-1,1:nhar) :: plm 1363 CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header 1364 CHARACTER(LEN=*), INTENT(IN) :: filename 1365 ! 1366 INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1) 1367 INTEGER(I4B) :: i 1368 LOGICAL(LGT) :: simple,extend 1369 CHARACTER(LEN=80) :: comment 1370 1371 INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension 1372 INTEGER(I4B) :: nrows, tfields, varidat, nmmax 1373 INTEGER(I8B) :: nlm, nf, i0, i1 1374 INTEGER(I4B) :: frow, felem, colnum 1375 CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname 1376 CHARACTER(LEN=10) :: card, num 1377 CHARACTER(LEN=2) :: stn 1378 INTEGER(I4B) :: itn, irow, np 1379 real(DP) :: dm 1380 !----------------------------------------------------------------------- 1381 1382 status=0 1383 1384 unit = 148 1385 1386 ! create the new empty FITS file 1387 blocksize=1 1388 call ftinit(unit,filename,blocksize,status) 1389 if (status > 0) call fatal_error("Error while creating file " & 1390 & //trim(filename) & 1391 & //". Check path and/or access rights.") 1392 1393 ! ----------------------------------------------------- 1394 ! initialize parameters about the FITS image 1395 simple=.true. 1396 bitpix=32 ! integer*4 1397 naxis=0 ! no image 1398 naxes(1)=0 1399 extend=.true. ! there is an extension 1400 1401 ! ---------------------- 1402 ! primary header 1403 ! ---------------------- 1404 ! write the required header keywords 1405 call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status) 1406 1407 ! writes supplementary keywords : none 1408 1409 ! write the current date 1410 call ftpdat(unit,status) ! (format ccyy-mm-dd) 1411 1412 ! ---------------------- 1413 ! image : none 1414 ! ---------------------- 1415 1416 ! ---------------------- 1417 ! extension 1418 ! ---------------------- 1419 1420 nlm = nplm / (2*nsmax) ! number of Plm per Healpix ring 1421 dm = (2*nlmax + 3.0_dp)**2 - 8.0_dp*nlm 1422 nmmax = nlmax - (nint(sqrt(dm)) - 1)/2 1423 call assert(nplm == nlm*2*nsmax, 'inconsistent array size in write_plm') 1424 call assert(nplm == (nmmax+1_i8b)*(2*nlmax-nmmax+2_i8b)*nsmax, & 1425 & 'inconsistent array size (nmmax) in write_plm') 1426 1427 ! creates an extension 1428 call ftcrhd(unit, status) 1429 1430 ! writes required keywords 1431 nrows = 2*nsmax ! naxis1 1432 tfields = nhar 1433 nf = nplm / nrows 1434 if (nf * nrows /= nplm) then 1435 print*,' problems in write_plm',nf*nrows,nplm 1436 stop 1437 endif 1438 write (num,'(I10)') nf 1439! tform(1:nhar) = trim(adjustl(num))//'D' ! does not work with Ifort, EH, 2006-04-04 1440 num = trim(adjustl(num))//'D' 1441 tform(1:nhar) = num 1442 ttype(1:nhar) = 'harmonics' ! will be updated 1443 tunit(1:nhar) = '' ! optional, will not appear 1444 extname = '' ! optional, will not appear 1445 varidat = 0 1446 call ftphbn(unit, nrows, tfields, ttype, tform, tunit, & 1447 & extname, varidat, status) 1448 1449 ! write the header literally, putting TFORM1 at the desired place 1450 comment = 'data format of field: 8-byte REAL' 1451 do i=1,nlheader 1452 card = header(i) 1453 if (card(1:5) == 'TTYPE') then ! if TTYPE1 is explicitely given 1454 stn = card(6:6) 1455 read(stn,'(i1)') itn 1456 ! discard at their original location: 1457 call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and 1458 status = 0 1459 call ftdkey(unit,'TFORM'//stn,status) ! TFORMi 1460 status = 0 1461 if (itn <= tfields) then ! only put relevant information 2008-08-27 1462 call putrec(unit,header(i), status) ! write new TTYPE1 1463 status = 0 1464 call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after 1465 endif 1466 elseif (header(i)/=' ') then 1467 call putrec(unit,header(i), status) 1468 endif 1469 status = 0 1470 enddo 1471 1472 ! make sure keywords critical for future use are present and correct 1473 call ftukyj(unit,'NSIDE', nsmax, 'HEALPIX resolution parameter', status) 1474 call ftukyj(unit,'MAX-LPOL',nlmax, 'Maximum multipole order l of P_lm', status) 1475 call ftukyj(unit,'MAX-MPOL',nmmax, 'Maximum degree m of P_lm', status) 1476 call ftukyl(unit,'POLAR', (nhar>1),'Polarisation included (T/F)', status) 1477 1478 ! write the extension one column by one column 1479 do irow=1, nrows 1480 frow = irow ! starting position (row) 1481 felem = 1 ! starting position (element) 1482 np = nf 1483 i0 = (irow-1)*nf 1484 i1 = i0+nf-1 1485 do colnum = 1, nhar 1486 !call ftpcld(unit, colnum, frow, felem, nplm, plm(0,colnum), status) 1487 call ftpcld(unit, colnum, frow, felem, np, plm(i0:i1,colnum), status) 1488 enddo 1489 enddo 1490 1491 ! ---------------------- 1492 ! close and exit 1493 ! ---------------------- 1494 1495 ! close the file and free the unit number 1496 call ftclos(unit, status) 1497 1498 ! check for any error, and if so print out error messages 1499 if (status > 0) call printerror(status) 1500 1501 return 1502 end subroutine write_plm8 1503 1504 !======================================================================= 1505 ! ALMS2FITS 1506 ! Writes alms from to binary FITS file, FKH/Apr-99 1507 ! ncl is the number of columns, in the output fits file, 1508 ! either 3 or 5 (with or without errors respectively) 1509 ! 1510 ! input array (real) FITS file 1511 ! alms(:,1) = l ) 1512 ! alms(:,2) = m )---> col 1: l*l+l+m+1 1513 ! alms(:,3) = real(a_lm) ---> col 2 1514 ! alms(:,4) = imag(a_lm) ---> col 3 1515 ! alms(:,5) = real(delta a_lm) ---> col 4 1516 ! alms(:,6) = imag(delta a_lm) ---> col 5 1517 !======================================================================= 1518 subroutine alms2fits_s(filename, nalms, alms, ncl, header, nlheader, next) 1519 !======================================================================= 1520 character(len=*), intent(in) :: filename 1521 integer(I4B), intent(in) :: nalms, nlheader, next, ncl 1522 real(SP), intent(in), dimension(1:nalms,1:ncl+1,1:next), target :: alms 1523 character(len=80), intent(in), dimension(1:nlheader,1:next) :: header 1524 integer(I4B) :: i 1525 real(SP), dimension(:,:), pointer :: palms 1526 1527 do i=1,next 1528 palms => alms(1:nalms,1:(ncl+1),i) ! use pointer to reduce memory usage 1529 call write_alms_s(filename, nalms, palms, ncl, & 1530 & header(1:nlheader,i), nlheader, i) 1531 enddo 1532 1533 end subroutine alms2fits_s 1534 !======================================================================= 1535 subroutine alms2fits_d(filename, nalms, alms, ncl, header, nlheader, next) 1536 !======================================================================= 1537 character(len=*), intent(in) :: filename 1538 integer(I4B), intent(in) :: nalms, nlheader, next, ncl 1539 real(DP), intent(in), dimension(1:nalms,1:ncl+1,1:next), target :: alms 1540 character(len=80), intent(in), dimension(1:nlheader,1:next) :: header 1541 integer(I4B) :: i 1542 real(DP), dimension(:,:), pointer :: palms 1543 1544 do i=1,next 1545 palms => alms(1:nalms,1:(ncl+1),i) ! use pointer to reduce memory usage 1546 call write_alms_d(filename, nalms, palms, ncl, & 1547 & header(1:nlheader,i), nlheader, i) 1548 enddo 1549 1550 end subroutine alms2fits_d 1551 1552 !======================================================================= 1553 subroutine fits2alms_s(filename, nalms, alms, ncl, header, nlheader, next) 1554 !======================================================================= 1555 character(len=*), intent(in) :: filename 1556 integer(I4B), intent(in) :: nalms, nlheader, next, ncl 1557 real(SP), intent(inout), dimension(1:nalms,1:ncl+1,1:next), target :: alms 1558 character(len=80), intent(out), dimension(1:nlheader,1:next) :: header 1559 integer(I4B) :: i 1560 real(SP), dimension(:,:), pointer :: palms 1561 1562 do i=1,next 1563 palms => alms(1:nalms,1:(ncl+1),i) ! use pointer to reduce memory usage 1564 call read_alms_s(filename, nalms, palms, ncl, & 1565 & header(1:nlheader,i), nlheader, i) 1566 enddo 1567 1568 end subroutine fits2alms_s 1569 !======================================================================= 1570 subroutine fits2alms_d(filename, nalms, alms, ncl, header, nlheader, next) 1571 !======================================================================= 1572 character(len=*), intent(in) :: filename 1573 integer(I4B), intent(in) :: nalms, nlheader, next, ncl 1574 real(DP), intent(inout), dimension(1:nalms,1:ncl+1,1:next), target :: alms 1575 character(len=80), intent(out), dimension(1:nlheader,1:next) :: header 1576 integer(I4B) :: i 1577 real(DP), dimension(:,:), pointer :: palms 1578 1579 do i=1,next 1580 palms => alms(1:nalms,1:(ncl+1),i) ! use pointer to reduce memory usage 1581 call read_alms_d(filename, nalms, palms, ncl, & 1582 & header(1:nlheader,i), nlheader, i) 1583 enddo 1584 1585 end subroutine fits2alms_d 1586 1587 SUBROUTINE input_tod_s(filename, tod, npixtot, ntods, & 1588 & header, firstpix, fmissval, extno) 1589 ! ************************************************************************* 1590 1591 ! ************************************************************************* 1592 ! v1.0 : Olivier Dore et Eric Hivon, June-July 2002 1593 ! v1.1 2002-07-08 : bugs correction by E.H. 1594 ! v1.2 2002-08-22 : pixel number starts at 0 for these routine, 1595 ! but starts at 1 for cftisio routines ... 1596 ! ************************************************************************* 1597 ! =================================================================================== 1598 ! reads fits file 1599 ! filename = fits file (input) 1600 ! tod = data rad from the file (ouput) = real*4 array of size (npixtot,ntods) 1601 ! npixtot = length of the tod (input) 1602 ! ntods = number of tods 1603 ! header = OPTIONAL, FITS header 1604 ! fmissval = OPTIONAL argument (input) with the value to be given to missing 1605 ! pixels, its default value is 0 1606 ! firstpix = OPTIONAL first pixel to be read (starting at 0) 1607 ! extno = OPTIONAL, extension number (starting at 0) 1608 ! O.D. & E.H. 06/02 @ IAP 1609 ! 2002-07-08 : bugs correction by E.H. 1610 ! (consistent use of fmiss_effct) 1611 ! ======================================================================= 1612 1613 IMPLICIT NONE 1614 1615 INTEGER(I8B), INTENT(IN) :: npixtot 1616 INTEGER(I8B), INTENT(IN), OPTIONAL :: firstpix 1617 INTEGER(I4B), INTENT(IN) :: ntods 1618 REAL(SP), INTENT(OUT) ,dimension(0:,1:) :: tod 1619 CHARACTER(LEN=*), INTENT(IN) :: filename 1620 REAL(SP), INTENT(IN), OPTIONAL :: fmissval 1621 character(len=*), intent(out),OPTIONAL, dimension(1:) :: header 1622 INTEGER(I4B), INTENT(IN), OPTIONAL :: extno 1623 1624 INTEGER(I4B) :: i,itod 1625 REAL(SP) :: fmissing, fmiss_effct 1626 INTEGER(I8B) :: imissing 1627 integer(i8b) :: fp = 0_i8b 1628 1629 LOGICAL(LGT) :: anynull 1630 1631 !----------------------------------------------------------------------- 1632 fmiss_effct = 0. 1633 IF (PRESENT(fmissval)) fmiss_effct = fmissval 1634 1635 if (present(firstpix)) fp = firstpix 1636 1637 CALL read_bintod_s(filename, tod, npixtot, ntods, fp, fmissing, anynull, & 1638 & header, extno) 1639 1640 call map_bad_pixels(tod, fmissing, fmiss_effct, imissing, verbose=.true.) 1641! DO itod = 1, ntods 1642! anynull = .TRUE. 1643! IF (anynull) THEN 1644! imissing = 0 1645! DO i=0,npixtot-1 1646! IF ( ABS(tod(i,itod)/fmissing -1.) .LT. 1.e-5 ) THEN 1647! tod(i,itod) = fmiss_effct 1648! imissing = imissing + 1 1649! ENDIF 1650! ENDDO 1651! IF (imissing .GT. 0) THEN 1652! WRITE(*,'(a,1pe11.4)') 'blank value : ' ,fmissing 1653! WRITE(*,'(i7,a,f7.3,a,1pe11.4)') & 1654! & imissing,' missing pixels (', & 1655! & (100.*imissing)/npixtot,' %),'// & 1656! & ' have been set to : ',fmiss_effct 1657! ENDIF 1658! ENDIF 1659! ENDDO 1660 RETURN 1661 1662 END SUBROUTINE input_tod_s 1663 !======================================================================= 1664 1665 !************************************************************************************** 1666 SUBROUTINE input_tod_d(filename, tod, npixtot, ntods, & 1667 & header, firstpix, fmissval, extno) 1668 !************************************************************************************** 1669 1670 !======================================================================= 1671 ! reads fits file 1672 ! filename = fits file (input) 1673 ! tod = data rad from the file (ouput) = real*8 array of size (npixtot,ntods) 1674 ! npixtot = length of the tod (input) 1675 ! ntods = number of tods 1676 ! header = OPTIONAL, FITS header 1677 ! fmissval = OPTIONAL argument (input) with the value to be given to missing 1678 ! pixels, its default value is 0 1679 ! firstpix = OPTIONAL first pixel to be read (starting at 0) 1680 ! extno = OPTIONAL, extension number (starting at 0) 1681 ! O.D. & E.H. 06/02 @ IAP 1682 ! 2002-07-08 : bugs correction by E.H. 1683 ! (consistent use of fmiss_effct) 1684 !======================================================================= 1685 1686 IMPLICIT NONE 1687 1688 INTEGER(I8B), INTENT(IN) :: npixtot 1689 INTEGER(I8B), INTENT(IN), OPTIONAL :: firstpix 1690 INTEGER(I4B), INTENT(IN) :: ntods 1691 REAL(DP), INTENT(OUT) ,dimension(0:,1:) :: tod 1692 CHARACTER(LEN=*), INTENT(IN) :: filename 1693 REAL(DP), INTENT(IN), OPTIONAL :: fmissval 1694 character(len=*), intent(out),OPTIONAL, dimension(1:) :: header 1695 INTEGER(I4B), INTENT(IN), OPTIONAL :: extno 1696 1697 INTEGER(I4B) :: i,itod 1698 REAL(DP) :: fmissing, fmiss_effct 1699 INTEGER(I8B) :: imissing 1700 integer(i8b) :: fp = 0_i8b 1701 1702 LOGICAL(LGT) :: anynull 1703 1704 !----------------------------------------------------------------------- 1705 1706 fmiss_effct = 0. 1707 IF (PRESENT(fmissval)) fmiss_effct = fmissval 1708 1709 if (present(firstpix)) fp = firstpix 1710 1711 CALL read_bintod_d(filename, tod, npixtot, ntods, fp, fmissing, anynull, & 1712 & header, extno) 1713 1714 call map_bad_pixels(tod, fmissing, fmiss_effct, imissing, verbose=.true.) 1715! DO itod = 1, ntods 1716! anynull = .TRUE. 1717! IF (anynull) THEN 1718! imissing = 0 1719! DO i=0,npixtot-1 1720! IF ( ABS(tod(i,itod)/fmissing -1.) .LT. 1.e-5 ) THEN 1721! tod(i,itod) = fmiss_effct 1722! imissing = imissing + 1 1723! ENDIF 1724! ENDDO 1725! IF (imissing .GT. 0) THEN 1726! WRITE(*,'(a,1pe11.4)') 'blank value : ' ,fmissing 1727! WRITE(*,'(i7,a,f7.3,a,1pe11.4)') & 1728! & imissing,' missing pixels (', & 1729! & (100.*imissing)/npixtot,' %),'// & 1730! & ' have been set to : ',fmiss_effct 1731! ENDIF 1732! ENDIF 1733! ENDDO 1734 RETURN 1735 1736 END SUBROUTINE input_tod_d 1737 !======================================================================= 1738 !======================================================================= 1739 subroutine read_dbintab(filename,map,npixtot,nmaps,nullval,anynull, units) 1740 !======================================================================= 1741 ! Read a FITS file 1742 ! 1743 ! slightly modified to deal with vector column 1744 ! in binary table EH/IAP/Jan-98 1745 ! 1746 ! Reads a double-precision array with precomputed plms, used by syn/anafast 1747 ! FKH/Apr-99 1748 !======================================================================= 1749 CHARACTER(LEN=*), INTENT(IN) :: filename 1750 INTEGER(I4B), INTENT(IN) :: npixtot, nmaps 1751 ! REAL(DP), DIMENSION(0:npixtot-1,1:nmaps), INTENT(OUT) :: map 1752 REAL(DP), DIMENSION(0:,1:), INTENT(OUT) :: map 1753 REAL(DP), INTENT(OUT) :: nullval 1754 LOGICAL(LGT), INTENT(OUT) :: anynull 1755 ! CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: header 1756 CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units 1757 1758 INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis 1759 INTEGER(I4B) :: group,firstpix,npix,i 1760 REAL(SP) :: blank, testval 1761 REAL(DP) :: bscale,bzero 1762 CHARACTER(LEN=80) :: comment 1763 LOGICAL(LGT) :: extend 1764 INTEGER(I4B) :: nmove, hdutype 1765 INTEGER(I4B) :: column, frow, imap 1766 INTEGER(I4B) :: datacode, repeat, width 1767 1768 INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension 1769 INTEGER(I4B) :: nrows, tfields, varidat 1770 CHARACTER(LEN=20), dimension(1:MAXDIM) :: ttype, tform, tunit 1771 CHARACTER(LEN=20) :: extname 1772 !----------------------------------------------------------------------- 1773 status=0 1774 1775 unit = 147 1776 naxes(1) = 1 1777 naxes(2) = 1 1778 nfound = -1 1779 anynull = .false. 1780 bscale = 1.0d0 1781 bzero = 0.0d0 1782 blank = -2.e25 1783 nullval = bscale*blank + bzero 1784 1785 readwrite=0 1786 call ftopen(unit,filename,readwrite,blocksize,status) 1787 if (status > 0) call printerror(status) 1788 ! ----------------------------------------- 1789 1790 ! determines the presence of image 1791 call ftgkyj(unit,'NAXIS', naxis, comment, status) 1792 if (status > 0) call printerror(status) 1793 1794 ! determines the presence of an extension 1795 call ftgkyl(unit,'EXTEND', extend, comment, status) 1796 if (status > 0) status = 0 ! no extension : 1797 ! to be compatible with first version of the code 1798 1799 if (naxis > 0) then ! there is an image 1800 ! determine the size of the image (look naxis1 and naxis2) 1801 call ftgknj(unit,'NAXIS',1_i4b,2_i4b,naxes,nfound,status) 1802 1803 ! check that it found only NAXIS1 1804 if (nfound == 2 .and. naxes(2) > 1) then 1805 print *,'multi-dimensional image' 1806 print *,'expected 1-D data.' 1807 call fatal_error 1808 end if 1809 1810 if (nfound < 1) then 1811 call printerror(status) 1812 print *,'can not find NAXIS1.' 1813 call fatal_error 1814 endif 1815 1816 npix=naxes(1) 1817 if (npix /= npixtot) then 1818 print *,'found ',npix,' plms' 1819 print *,'expected ',npixtot 1820 call fatal_error 1821 endif 1822 1823 call ftgkyd(unit,'BSCALE',bscale,comment,status) 1824 if (status == 202) then ! BSCALE not found 1825 bscale = 1.0d0 1826 status = 0 1827 endif 1828 call ftgkyd(unit,'BZERO', bzero, comment,status) 1829 if (status == 202) then ! BZERO not found 1830 bzero = 0.0d0 1831 status = 0 1832 endif 1833 call ftgkye(unit,'BLANK', blank, comment,status) 1834 if (status == 202) then ! BLANK not found 1835 ! (according to fitsio BLANK is integer) 1836 blank = -2.e25 1837 status = 0 1838 endif 1839 nullval = bscale*blank + bzero 1840 1841 ! ----------------------------------------- 1842 1843 group=1 1844 firstpix = 1 1845 call ftgpvd(unit,group,firstpix,npix,nullval,map,anynull,status) 1846 ! if there are any NaN pixels, (real data) 1847 ! or BLANK pixels (integer data) they will take nullval value 1848 ! and anynull will switch to .true. 1849 ! otherwise, switch it by hand if necessary 1850 testval = 1.e-6 * ABS(nullval) 1851 do i=0, npix-1 1852 if (ABS(map(i,1)-nullval) < testval) then 1853 anynull = .true. 1854 goto 111 1855 endif 1856 enddo 1857111 continue 1858 1859 else if (extend) then ! there is an extension 1860 nmove = +1 1861 call ftmrhd(unit, nmove, hdutype, status) 1862 !cc write(*,*) hdutype 1863 1864 call assert (hdutype==2, 'this is not a binary table') 1865 1866 ! reads all the keywords 1867 call ftghbn(unit, MAXDIM, & 1868 & nrows, tfields, ttype, tform, tunit, extname, varidat, & 1869 & status) 1870 1871 if (tfields < nmaps) then 1872 print *,'found ',tfields,' maps in the file' 1873 print *,'expected ',nmaps 1874 call fatal_error 1875 endif 1876 1877 ! finds the bad data value 1878 call ftgkyd(unit,'BAD_DATA',nullval,comment,status) 1879 if (status == 202) then ! bad_data not found 1880 nullval = d_bad_value 1881 status = 0 1882 endif 1883 1884 ! if (nlheader > 0) then 1885 ! call get_clean_header(unit, header, filename, status) 1886 ! endif 1887 1888 if (present(units)) then 1889 units = 'unknown' 1890 do imap = 1, nmaps 1891 units(imap) = adjustl(tunit(imap)) 1892 enddo 1893 endif 1894 1895 do imap = 1, nmaps 1896 !parse TFORM keyword to find out the length of the column vector 1897 call ftbnfm(tform(imap), datacode, repeat, width, status) 1898 1899 !reads the columns 1900 column = imap 1901 frow = 1 1902 firstpix = 1 1903 npix = nrows * repeat 1904 if (npix /= npixtot) then 1905 print *,'found ',npix,' plms' 1906 print *,'expected ',npixtot 1907 call fatal_error("read_dbintab "//trim(filename)) 1908 endif 1909 call ftgcvd(unit, column, frow, firstpix, npix, nullval, & 1910 & map(0:npix-1,imap), anynull, status) 1911 enddo 1912 1913 else ! no image no extension, you are dead, man 1914 call fatal_error(' No image, no extension') 1915 endif 1916 1917 ! close the file 1918 call ftclos(unit, status) 1919 1920 ! check for any error, and if so print out error messages 1921 if (status > 0) call printerror(status) 1922 1923 return 1924 end subroutine read_dbintab 1925 1926 1927 !======================================================================= 1928 subroutine printerror(status) 1929 !======================================================================= 1930 ! Print out the FITSIO error messages to the user 1931 !======================================================================= 1932 INTEGER(I4B), INTENT(IN) :: status 1933 CHARACTER :: errtext*30,errmessage*80 1934 !----------------------------------------------------------------------- 1935 ! check if status is OK (no error); if so, simply return 1936 if (status .le. 0)return 1937 1938 ! get the text string which describes the error 1939 call ftgerr(status,errtext) 1940 print *,'FITSIO Error Status =',status,': ',errtext 1941 1942 ! read and print out all the error messages on the FITSIO stack 1943 call ftgmsg(errmessage) 1944 do while (errmessage /= ' ') 1945 print *,errmessage 1946 call ftgmsg(errmessage) 1947 end do 1948 1949 return 1950 end subroutine printerror 1951 1952 1953 !======================================================================= 1954 subroutine read_par(filename,nside,lmax,tfields,mmax) 1955 !======================================================================= 1956 ! Read nside, lmax, tfields and mmax from a FITS file 1957 ! parameters not found take a value of -1 1958 ! 1959 ! Frode K. Hansen, April 1999 1960 ! EH, Dec 2004 1961 ! 1962 !======================================================================= 1963 CHARACTER(LEN=*), INTENT(IN) :: filename 1964 1965 INTEGER(I4B), INTENT(OUT) :: nside, lmax, tfields 1966 integer(i4b), intent(out), optional :: mmax 1967 1968 INTEGER(I4B) :: status,unit,readwrite,blocksize,naxis 1969 CHARACTER(LEN=80) :: comment, ttype1 1970 LOGICAL(LGT) :: extend, anyf 1971 INTEGER(I4B):: nmove, hdutype, idmax, nrows 1972 1973 !----------------------------------------------------------------------- 1974 status=0 1975 unit = 120 1976 comment='' 1977 ttype1='' 1978 1979 readwrite=0 1980 call ftopen(unit,filename,readwrite,blocksize,status) 1981 if (status > 0) call printerror(status) 1982 ! ----------------------------------------- 1983 call ftgkyj(unit,'NAXIS', naxis, comment, status) 1984 1985 ! determines the presence of an extension 1986 call ftgkyl(unit,'EXTEND', extend, comment, status) 1987 call assert (status<=0, 'No Extension in FITS file!') 1988 1989 nmove = +1 1990 call ftmrhd(unit, nmove, hdutype, status) 1991 1992 call assert (hdutype==2, 'This is not a FITS binary-table') 1993 1994 call ftgkyj(unit,'NSIDE',nside,comment,status) 1995 if (status == 202) then 1996 print*,'WARNING: NSIDE keyword not found!' 1997 nside = -1 1998 status = 0 1999 endif 2000 2001 call ftgkyj(unit,'TFIELDS',tfields,comment,status) 2002 if (status == 202) then 2003 print*,'WARNING: TFIELDS keyword not found!' 2004 tfields = -1 2005 status = 0 2006 endif 2007 2008 call ftgkyj(unit,'MAX-LPOL',lmax,comment,status) 2009 if (status == 202) then 2010 status = 0 2011 ! if not found, determines if file contains indexed list of alm 2012 ! if so, find out largest l 2013 if (tfields >= 3 .and. hdutype==2) then ! 3 column binary table 2014 call ftgkys(unit,'TTYPE1',ttype1,comment,status) 2015 ttype1 = trim(strupcase(adjustl(ttype1))) 2016 if (trim(ttype1(1:5)) == 'INDEX') then 2017 call ftgkyj(unit, 'NAXIS2', nrows, comment, status) ! find number of rows 2018 call ftgcvj(unit, 1_i4b, nrows, 1_i4b, 1_i4b, 0_i4b, idmax, anyf, status) ! read element on last row of first column 2019 if (status == 0) then 2020 lmax = int(sqrt( real(idmax-1, kind = DP) ) ) 2021 if (lmax > 0) goto 1000 2022 endif 2023 endif 2024 endif 2025 print*,'WARNING: MAX-LPOL keyword not found!' 2026 lmax = -1 2027 status = 0 2028 endif 20291000 continue 2030 2031 if (present(mmax)) then 2032 call ftgkyj(unit,'MAX-MPOL',mmax,comment,status) 2033 if (status == 202) then 2034 print*,'WARNING: MAX-MPOL keyword not found!' 2035 mmax = -1 2036 status = 0 2037 endif 2038 endif 2039 call ftclos(unit, status) 2040 2041 end subroutine read_par 2042 2043 !======================================================================= 2044 function getnumext_fits(filename) 2045 !======================================================================= 2046 ! result = getnumext_fits(filename) 2047 ! returns the number of extensions present in FITS file 'filename' 2048 ! 2049 ! EH, Nov 2004 2050 ! April 2007: close file on exit 2051 !======================================================================= 2052 character(LEN=*), intent(IN) :: filename 2053 integer(i4b) :: getnumext_fits 2054 ! 2055 integer(i4b) :: status, unit, readwrite, blocksize, nhdu 2056 !----------------------------------------------------------------------- 2057 status = 0 2058 unit = 149 2059 getnumext_fits = 0 2060 readwrite = 0 ! Read only 2061 call ftopen(unit, filename, readwrite, blocksize, status) 2062 if (status > 0) then 2063 call printerror(status) 2064 call ftclos(unit, status) 2065 return 2066 endif 2067 2068 call ftthdu(unit, nhdu, status) 2069 getnumext_fits = nhdu - 1 2070 2071 call ftclos(unit, status) 2072 return 2073 end function getnumext_fits 2074 !======================================================================= 2075 function getsize_fits(filename, nmaps, ordering, obs_npix, & 2076 & nside, mlpol, type, polarisation, fwhm_arcmin, beam_leg, & 2077 & coordsys, polcconv, extno) 2078 !======================================================================= 2079 ! result = getsize_fits(filename, nmaps, ordering, nside, mlpol, type) 2080 ! Get the number of pixels stored in a map FITS file. 2081 ! Each pixel is counted only once 2082 ! (even if several information is stored on each of them, see nmaps). 2083 ! Depending on the data storage format, this may be : 2084 ! - equal or smaller to the number Npix of Healpix pixels available 2085 ! over the sky for the given resolution (Npix = 12*nside*nside) 2086 ! - equal or larger to the number of non blank pixels (obs_npix) 2087 ! 2088 ! filename = (IN) name of the FITS file containing the map 2089 ! nmaps = (OPTIONAL, OUT) number of maps in the file 2090 ! ordering = (OPTIONAL, OUT) Healpix ordering scheme used 2091 ! 0 = unknown 2092 ! 1 = RING 2093 ! 2 = NESTED 2094 ! obs_npix = (OPTIONAL, OUT) number of non blanck pixels. 2095 ! It is set to -1 if it can not be determined from header 2096 ! information alone 2097 ! nside = (OPTIONAL, OUT) Healpix parameter Nside 2098 ! returns a negative value if not found 2099 ! mlpol = (OPTIONAL, OUT) maximum multipole used to generate the map (for simulated map) 2100 ! returns a negative value if not found 2101 ! type = (OPTIONAL, OUT) Healpix/FITS file type 2102 ! <0 : file not found, or not valid 2103 ! 0 : image only fits file, deprecated Healpix format 2104 ! (result = 12 * nside * nside) 2105 ! 1 : ascii table, generally used for C(l) storage 2106 ! 2 : binary table : with implicit pixel indexing (full sky) 2107 ! (result = 12 * nside * nside) 2108 ! 3 : binary table : with explicit pixel indexing (generally cut sky) 2109 ! (result <= 12 * nside * nside) 2110 ! 999 : unable to determine the type 2111 ! polarisation = (OPTIONAL, OUT) presence of polarisation data in the file 2112 ! <0 : can not find out 2113 ! 0 : no polarisation 2114 ! 1 : contains polarisation (Q,U or G,C) 2115 ! fwhm_arcmin = (OPTIONAL, DP, OUT) returns the beam FWHM read from FITS header, 2116 ! translated from Deg (hopefully) to arcmin 2117 ! returns a negative value if not found 2118 ! 2119 ! beam_leg = (OPTIONAL, CHR, OUT) filename of beam or filtering window function applied to data 2120 ! returns a empty string if not found 2121 ! coordsys = (OPTIONAL, CHR, OUT) string describing coordinate system, 2122 ! G = Galactic, E = ecliptic, C = celestial = equatorial. 2123 ! empty if not found. 2124 ! 2125 ! polcconv = (OPTIONAL, I4B, OUT) coordinate convention for polarisation 2126 ! 0: not found 2127 ! 1: COSMO (default for Healpix) 2128 ! 2: IAU 2129 ! 3: neither COSMO nor IAU 2130 ! 2131 ! extno = (OPTIONAL, IN) specify FITS extension to look at (0 based) 2132 ! default = 0 (first extension) 2133 ! 2134 ! Benjamin D. Wandelt January 1998 2135 ! includes Eric Hivon's modification of FITS columns 2136 ! and draws heavily on the read_bintab routine. 2137 ! 2138 ! addition of optional Ordering by E.H. (July 98) 2139 ! addition of optional Nside and Mlpol by ?? (??) 2140 ! addition of optional type by E.H. (Sept 00) 2141 ! improved for ASCII table E.H (Feb 03) 2142 ! addition of extno, E.H (Nov 04) 2143 ! addition of fwhm_arcmin and beam_leg, EH (Jan 05). 2144 ! addition of polcconv, EH (June 05). 2145 ! accept files with NAXIS2 > 2^31 (2015-04) (requires recent cfitsio >= 3.20) 2146 ! accept files with polarisation columns named Q_S*, Q-S*, ... U_S*, U-S* 2147 !======================================================================= 2148 use pix_tools, only: nside2npix 2149 character(LEN=*), intent(IN) :: filename 2150 integer(I4B), intent(out), optional :: nmaps 2151 integer(I4B), intent(out), optional :: ordering 2152 integer(I4B), intent(out), optional :: nside 2153 integer(I4B), intent(out), optional :: mlpol 2154 integer(I4B), intent(out), optional :: obs_npix 2155 integer(I4B), intent(out), optional :: type 2156 integer(I4B), intent(out), optional :: polarisation 2157 real(DP), intent(out), optional :: fwhm_arcmin 2158 character(LEN=*), intent(out), optional :: beam_leg 2159 character(LEN=*), intent(out), optional :: coordsys 2160 integer(I4B), intent(out), optional :: polcconv 2161 integer(I4B), intent(in), optional :: extno 2162 2163 INTEGER(I4B) :: nmaps_in, ordering_in, nside_in 2164 INTEGER(I4B) :: mlpol_in, obs_npix_in, ftype_in 2165 INTEGER(I4B) :: extno_in, polcconv_in 2166 real(DP) :: fwhm_arcmin_in 2167 character(len=FILENAMELEN) :: beam_leg_in 2168 character(len=20) :: coordsys_in 2169 INTEGER(I4B) :: grain 2170 CHARACTER(len=20) :: indxschm 2171 CHARACTER(LEN=20) :: order_val, object_val, ttype_val, polcconv_val 2172! INTEGER(I4B) :: getsize_fits 2173 INTEGER(I8B) :: getsize_fits 2174 LOGICAL(LGT) :: polar_in 2175 character(len=3), dimension(1:16,1:2) :: defpol 2176 logical(kind=LGT), dimension(1:2) :: pf 2177 integer(kind=I4B) :: ndp, j, k 2178 2179 INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis 2180 INTEGER(I4B) :: i 2181! INTEGER(I4B) :: nsize 2182 INTEGER(I8B) :: nsize 2183 CHARACTER(LEN=80) :: comment 2184 LOGICAL(LGT) :: extend 2185 INTEGER(I4B) :: nmove, hdutype, hdunum 2186 INTEGER(I4B) :: datacode, repeat1, repeat2, width 2187 2188 INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension 2189 INTEGER(I4B) :: tfields, varidat, rowlen 2190 CHARACTER(LEN=20), dimension(1:MAXDIM) :: ttype, tform, tunit 2191 INTEGER(I4B), dimension(1:MAXDIM) :: tbcol 2192 CHARACTER(LEN=20) :: extname 2193 integer(I4B) :: nrows32 2194 integer(I8B) :: nrows64 2195 !----------------------------------------------------------------------- 2196 status=0 2197 order_val = '' 2198 nmaps_in = 0 2199 ordering_in = 0 2200 mlpol_in = -1 2201 nside_in = -1 2202 ftype_in = 999 ! 2203 obs_npix_in = -1 2204 unit = 119 2205 naxes(1) = 1 2206 naxes(2) = 1 2207 nfound = -1 2208 extno_in = 0 2209 polcconv_in = 0 2210 comment='' 2211 ttype='' 2212 tform='' 2213 tunit='' 2214 extname='' 2215 if (present(extno)) extno_in = extno 2216 naxis = 0 2217 extend = .false. 2218 2219 readwrite=0 2220 !call ftopen(unit,filename,readwrite,blocksize,status) 2221 call ftnopn(unit,filename,readwrite,status) 2222 if (status > 0) then 2223 ftype_in = -1 2224 getsize_fits = -1 2225 call printerror(status) 2226 call ftclos(unit, status) 2227 return 2228 endif 2229 ! ----------------------------------------- 2230 call ftghdn(unit, hdunum) 2231 if (hdunum == 1) then ! in primary HDU: move to next HDU 2232 ! determines the presence of image 2233 call ftgkyj(unit,'NAXIS', naxis, comment, status) 2234 2235 ! determines the presence of an extension 2236 call ftgkyl(unit,'EXTEND', extend, comment, status) 2237 if (status > 0) then 2238 ftype_in = 0 2239 status = 0 ! no extension : 2240 ! to be compatible with first version of the code 2241 endif 2242 2243 else ! already in non primary HDU 2244 extend = .true. 2245 call ftghdt(unit, hdutype, status) 2246 endif 2247 2248 if (naxis > 0) then 2249 !--------------------------------- 2250 ! there is an image 2251 !--------------------------------- 2252 ! determine the size of the image (look naxis1 and naxis2) 2253 call ftgknj(unit,'NAXIS',1_i4b,2_i4b,naxes,nfound,status) 2254 2255 ! check that it found only NAXIS1 2256 if (nfound == 2 .and. naxes(2) > 1) then 2257 print *,'multi-dimensional image' 2258 print *,'expected 1-D data.' 2259 call ftclos(unit, status) 2260 call fatal_error 2261 end if 2262 2263 if (nfound < 1) then 2264 call printerror(status) 2265 print *,'can not find NAXIS1.' 2266 call ftclos(unit, status) 2267 call fatal_error 2268 endif 2269 2270 nsize=naxes(1) 2271 2272 call ftgkys(unit,'ORDERING',order_val,comment,status) 2273 if (status == 202) then ! Ordering not found 2274 ordering_in = 0 2275 order_val = '' 2276 status = 0 2277 endif 2278 2279 call ftgkyj(unit,'NSIDE',nside_in,comment,status) 2280 if (status == 202) then ! Nside not found 2281 nside_in = -1 2282 status = 0 2283 endif 2284 2285 if (present(mlpol)) then 2286 call ftgkyj(unit,'MAX-LPOL',mlpol_in,comment,status) 2287 if (status == 202) then ! max-lpol not found 2288 mlpol_in = -1 2289 status = 0 2290 endif 2291 endif 2292 2293 if (present(polarisation)) then 2294 polarisation = 0 2295 endif 2296 else if (extend) then 2297 !----------------------------------------------------------------- 2298 ! there is an extension 2299 !----------------------------------------------------------------- 2300 2301 nmove = extno_in 2302 if (hdunum == 1) nmove = extno_in + 1 2303 call ftmrhd(unit, nmove, hdutype, status) 2304 if (status > 0) then ! extension not found 2305 print*,'Extension #',extno_in,' not found in '//trim(filename) 2306 call printerror(status) 2307 call ftclos(unit, status) 2308 call fatal_error 2309 endif 2310 !c write(*,*) hdutype 2311 2312 ! reads all the keywords 2313 if (hdutype == 2) then ! binary table 2314 call ftghbn(unit, MAXDIM, & 2315 & nrows32, tfields, ttype, tform, tunit, & 2316 & extname, varidat, status) 2317 else ! ASCII table (hdutype = 1) 2318 ftype_in = 1 2319 call ftghtb(unit, MAXDIM, & 2320 & rowlen, & 2321 & nrows32, tfields, ttype, tbcol, tform, tunit, & 2322 & extname, status) 2323 endif 2324 ! deal with NAXIS2 > 2^31, 2015-04 2325#ifdef NO64BITS 2326 call ftgnrw (unit, nrows32, status) 2327 nrows64 = nrows32 2328#else 2329 call ftgnrwll(unit, nrows64, status) 2330#endif 2331 2332 ! parse TFORM keyword to find out the length of the column vector 2333 repeat1 = 1 2334 repeat2 = 1 2335 call ftbnfm(tform(1), datacode, repeat1, width, status) 2336 if (tfields > 1) call ftbnfm(tform(2), datacode, repeat2, width, status) 2337 2338 !nsize = int(nrows32, kind=i8b) * max(repeat1,repeat2) ! corrected Oct-03 2339 nsize = int(nrows64, kind=i8b) * max(repeat1,repeat2) ! 2015-04 2340 2341 nmaps_in = tfields 2342 2343 call ftgkys(unit,'ORDERING',order_val,comment,status) 2344 if (status == 202) then ! Ordering not found 2345 ordering_in = 0 2346 order_val = '' 2347 status = 0 2348 endif 2349 2350 call ftgkyj(unit,'NSIDE',nside_in,comment,status) 2351 if (status == 202) then ! Nside not found 2352 nside_in = -1 2353 status = 0 2354 endif 2355 2356 call ftgkyj(unit,'OBS_NPIX',obs_npix_in,comment,status) 2357 if (status == 202) then ! obs_npix not found 2358 obs_npix_in = -1 2359 status = 0 2360 endif 2361 2362 if (present(mlpol)) then 2363 call ftgkyj(unit,'MAX-LPOL',mlpol_in,comment,status) 2364 if (status == 202) then ! max-lpol not found 2365 status = 0 2366 call ftgkyj(unit,'LMAX',mlpol_in,comment,status) 2367 if (status == 202) then 2368 mlpol_in = -1 2369 status = 0 2370 endif 2371 endif 2372 endif 2373 2374 if (present(fwhm_arcmin)) then 2375 call ftgkyd(unit,'FWHM',fwhm_arcmin_in, comment, status) 2376 if (status == 202) then ! fwhm not found 2377 fwhm_arcmin_in = -1. 2378 status = 0 2379 else 2380 fwhm_arcmin_in = 60.0_dp * fwhm_arcmin_in 2381 endif 2382 endif 2383 2384 if (present(beam_leg)) then 2385 call ftgkys(unit,'BEAM_LEG',beam_leg_in, comment, status) 2386 if (status == 202) then ! beam_leg not found 2387 beam_leg_in = ' ' 2388 status = 0 2389 endif 2390 endif 2391 2392 if (present(coordsys)) then 2393 call ftgkys(unit,'COORDSYS',coordsys_in, comment, status) 2394 if (status == 202) then ! coordsys not found 2395 coordsys_in = ' ' 2396 status = 0 2397 endif 2398 endif 2399 2400 ! determines pixel indexing (for binary tables) 2401 if (present(type) .and. ftype_in == 999) then 2402 ! most stringent test 2403 call ftgkys(unit,'INDXSCHM',indxschm,comment,status) 2404 if (status == 0) then ! found 2405 ftype_in = 3 2406 if (trim(indxschm) == 'IMPLICIT') ftype_in = 2 2407 goto 1000 2408 else 2409 status = 0 2410 endif 2411 ! 2nd most stringent test 2412 call ftgkyj(unit,'GRAIN',grain,comment,status) 2413 if (status == 0) then ! found 2414 ftype_in = 3 2415 if (grain == 0) ftype_in = 2 2416 goto 1000 2417 else 2418 status = 0 2419 endif 2420 ! 3rd most stringent test 2421 if (trim(ttype(1)) /= '') then 2422 if (trim(ttype(1)) == 'PIXEL' .or. trim(ttype(1)) == 'PIX') then 2423 ftype_in = 3 2424 else 2425 ftype_in = 2 2426 endif 2427 goto 1000 2428 endif 2429 ! lousy test 2430 call ftgkys(unit,'OBJECT',object_val,comment,status) 2431 if (status == 0) then 2432 if (trim(object_val) == 'PARTIAL') ftype_in = 3 2433 if (trim(object_val) == 'FULLSKY') ftype_in = 2 2434 if (ftype_in /= 999) goto 1000 2435 else 2436 status = 0 2437 endif 2438 ! very lousy test 2439 if (nside_in > 0) then 2440 ftype_in = 3 2441 if (nside2npix(nside_in) == nsize) ftype_in = 2 2442 goto 1000 2443 endif 2444 endif 24451000 continue 2446 2447 ! find out if polarisation data is present 2448 if (present(polarisation)) then 2449 if (tfields < 3) then 2450 polarisation = 0 ! no polarisation 2451 goto 2000 2452 endif 2453 call ftgkyl(unit,'POLAR',polar_in,comment,status) 2454 if (status == 0) then ! polar found 2455 polarisation = 0 2456 if (polar_in) polarisation = 1 2457 goto 2000 2458 else ! polar not found 2459 status = 0 2460 polarisation = -1 2461 if (hdutype <= 0) goto 2000 2462 if (hdutype == 1) then ! ascii table -> power spectra 2463 ndp = 4 2464 defpol(1:ndp,1) = (/ "GRA","E-M","POW","EE " /) 2465 defpol(1:ndp,2) = (/ "CUR","B-M","POW","BB " /) 2466 endif 2467 if (hdutype == 2) then ! binary table -> maps or power spectra 2468 ndp = 10 2469 defpol(1:ndp,1) = (/ "Q-P","Q_P","Q P", "Q-S","Q_S", & 2470 & "Q ", "GRA","E-M","POW","EE " /) 2471 defpol(1:ndp,2) = (/ "U-P","U_P","U P", "U-S","U_S", & 2472 & "U ", "CUR","B-M","POW","BB " /) 2473 endif 2474 pf(:) = .false. 2475 do i = 2, tfields ! do not consider first field (generally temperature) 2476 ttype_val = adjustl(ttype(i)) 2477 call ftupch(ttype_val) ! upper case 2478 do k=1,2 2479 do j=1, ndp 2480 if (index( ttype_val, trim(defpol(j,k)) ) == 1) then 2481 pf(k) = .true. 2482 goto 1500 ! move to next field 2483 endif 2484 enddo 2485 enddo ! loop on k 24861500 continue 2487 enddo ! loop on i 2488 if (pf(1) .and. pf(2)) polarisation = 1 2489 if (.not. (pf(1) .or. pf(2))) polarisation = 0 2490 endif ! polar not found 2491 endif ! present(polarisation) 24922000 continue 2493 2494 call ftgkys(unit,'POLCCONV',polcconv_val,comment,status) 2495 if (status == 0) then 2496 polcconv_in = 3 ! neither COSMO nor IAU 2497 if (trim(polcconv_val) == 'COSMO') polcconv_in = 1 2498 if (trim(polcconv_val) == 'IAU') polcconv_in = 2 2499 else 2500 polcconv_in = 0 ! not found 2501 status = 0 2502 endif 2503 2504 else ! no image no extension, you are dead, man 2505 ftype_in = -1 2506 call fatal_error(' No image, no extension') 2507 endif 2508 2509 ! close the file 2510 call ftclos(unit, status) 2511 2512 ! check for any error, and if so print out error messages 2513 if (status > 0) call printerror(status) 2514 2515 getsize_fits=nsize 2516 2517 call ftupch(order_val) ! convert order_val to upper case 2518 if (order_val(1:4) == 'RING') ordering_in = 1 2519 if (order_val(1:4) == 'NEST') ordering_in = 2 2520 2521 if (present(nmaps)) nmaps = nmaps_in 2522 if (present(mlpol)) mlpol = mlpol_in 2523 if (present(obs_npix)) obs_npix = obs_npix_in 2524 if (present(ordering)) ordering = ordering_in 2525 if (present(nside)) nside = nside_in 2526 if (present(type)) type = ftype_in 2527 if (present(fwhm_arcmin)) fwhm_arcmin = fwhm_arcmin_in 2528 if (present(beam_leg)) beam_leg = adjustl(beam_leg_in) 2529 if (present(coordsys)) coordsys = adjustl(coordsys_in) 2530 if (present(polcconv)) polcconv = polcconv_in 2531 2532 return 2533 end function getsize_fits 2534 !======================================================================= 2535 function number_of_alms(filename, extnum) 2536 !======================================================================= 2537 ! Read the number of alms from a FITS-file containing alms 2538 ! for constrained realisations in synfast. 2539 ! 2540 ! Frode K. Hansen, April 1999 2541 ! EH. Jan 2004: use repeat information 2542 !======================================================================= 2543 character(len=*), intent(in) :: filename 2544 integer(I4B), optional, intent(inout) :: extnum 2545 integer(I4B) :: number_of_alms 2546 !integer(I8B) :: number_of_alms 2547 2548 integer(I4B), dimension(2) :: naxes 2549 integer(I4B) :: status, unit, readwrite, blocksize, naxis, nfound 2550 integer(I4B) :: nmove, hdutype, hdunum 2551 integer(I4B) :: datacode, repeat, width 2552 character(len=80) :: comment 2553 character(len=20) :: tform 2554 logical(LGT) :: extend 2555 2556 !----------------------------------------------------------------------- 2557 status=0 2558 unit = 118 2559 comment='' 2560 2561 readwrite=0 2562 call ftopen(unit, filename, readwrite, blocksize, status) 2563 if (status > 0) call printerror(status) 2564 ! ----------------------------------------- 2565 call ftgkyj(unit,'NAXIS', naxis, comment, status) 2566 2567 ! determines the presence of an extension 2568 call ftgkyl(unit,'EXTEND', extend, comment, status) 2569 call assert (status<=0, 'No Extension in FITS file!') 2570 2571 nmove = +1 2572 call ftmrhd(unit, nmove, hdutype, status) 2573 2574 call assert (hdutype==2, 'This is not a FITS binary-table') 2575 2576 call ftgknj(unit,'NAXIS', 1_i4b, 2_i4b, naxes, nfound, status) 2577 call assert (nfound>=2, 'NAXIS2-keyword not found!') 2578 2579 call ftgkys(unit,'TFORM1', tform, comment, status) 2580 call ftbnfm(tform, datacode, repeat, width, status) 2581 2582 number_of_alms = int(naxes(2),kind=i8b) * repeat 2583 2584 if (present(extnum)) then 2585 call ftthdu(unit,hdunum,status) 2586 extnum = hdunum - 1 ! first HDU : primary array 2587 endif 2588 2589 call ftclos(unit, status) 2590 2591 end function number_of_alms 2592 2593 !====================================================================== 2594 subroutine putrec(unit, card, status) 2595 !====================================================================== 2596 ! append, delete or update a card from the current header 2597 ! (deletion of keyword KWD is described by : - KWD) 2598 ! EH, version 1.0, Dec 2001 2599 !====================================================================== 2600 integer(kind=I4B), intent(IN) :: unit 2601 character(len=*), intent(IN) :: card 2602 integer(kind=I4B), intent(OUT) :: status 2603 2604 character(len=80) :: cardfits,record 2605 character(len=8) :: kwd 2606 integer(kind=I4B) :: hdtype 2607 character(len=80) :: nullarr(0) 2608 !====================================================================== 2609 2610 status = 0 2611 cardfits='' 2612 record='' 2613 call ftgthd(card, cardfits, hdtype, status) 2614 kwd = cardfits(1:8) 2615 status = 0 2616 2617 select case (hdtype) 2618 case (-1) ! delete keyword (starting from the top) 2619 call ftgrec(unit,0_i4b,record,status) 2620 ! call ftdkey(unit, kwd, status) 2621 ! patch around cfitsio bug 2622 ! (ftdkey does not recognize wild cards while ftgnxk does) 2623 do 2624 call ftgnxk(unit, kwd, 1_i4b, nullarr, 0_i4b, record, status) 2625 if (status /= 0) exit 2626 call ftdkey(unit, record(1:8), status) 2627 enddo 2628 case (0) ! append or update 2629 if (kwd == 'CONTINUE') then 2630 call ftprec(unit, trim(card), status) 2631 else 2632 ! if long string, put LongStringWarning 2633 if (index(card, "&'")>0) call ftplsw(unit, status) 2634 ! delete keyword in its current location (if any) 2635 call ftdkey(unit, kwd, status) 2636 status = 0 2637 ! append 2638 call ftprec(unit, trim(cardfits), status) 2639 endif 2640 case (1) ! append (for HISTORY and COMMENT) 2641 call ftprec(unit, trim(cardfits), status) 2642 case default 2643 write(unit=*,fmt=*)" Unexpected card format in fits header :" 2644 write(unit=*,fmt="(a80)") card 2645 write(unit=*,fmt=*)" card not written." 2646 end select 2647 status = 0 2648 2649 return 2650 end subroutine putrec 2651 2652 !==================================================================== 2653 subroutine get_clean_header(unit, header, filename, error, xalso, xonly) 2654 !==================================================================== 2655 ! get_clean_header(unit, header, error [, xalso, xonly]) 2656 ! gets the FITS header from unit, after excluding some default keywords 2657 ! defined in def_excl 2658 ! if header in non void on input, its content will be concatenated with that 2659 ! of the FITS header 2660 ! if xalso is defined as a list of keywords, they are also excluded from the header 2661 ! if xonly is defined as a list of keywords, only those keywords are excluded from 2662 ! the header. 2663 ! xonly and xalso are exclusive 2664 !==================================================================== 2665 INTEGER(I4B), intent(IN) :: unit 2666 CHARACTER(LEN=*), DIMENSION(1:), INTENT(IN OUT) :: header 2667 CHARACTER(LEN=*), INTENT(IN) :: filename 2668 INTEGER(I4B), intent(OUT) :: error 2669 character(len=8), dimension(1:), intent(IN), optional :: xalso 2670 character(len=8), dimension(1:), intent(IN), optional :: xonly 2671 2672 INTEGER(I4B) :: nlheader, status, i, n_excl 2673 CHARACTER(LEN=80) :: record 2674 CHARACTER(len=8), dimension(:), allocatable :: to_excl 2675 2676 CHARACTER(len=8), dimension(1:21) :: def_excl 2677 !==================================================================== 2678 2679 ! keywords to be excluded by default from output header 2680 ! Note that TTYPE# and TUNIT# keywords are not excluded 2681 ! from the output header as they might be useful to the 2682 ! calling routines 2683 def_excl=(/& 2684 & "SIMPLE ","BITPIX ","NAXIS ",& 2685 & "NAXIS# ","PCOUNT ","GCOUNT ",& 2686 & "EXTEND ","ORIGIN ","DATE* ",& 2687 & "TFIELDS ","TFORM# ", & 2688 & "TBCOL# ","EXTNAME ","CTYPE# ",& 2689 & "CRVAL# ","CRPIX# ","CDELT# ",& 2690 & "XTENSION","INSTRUME","TELESCOP",& 2691 & "PDMTYPE "/) 2692 2693 error = 0 2694 record='' 2695 2696 if (present(xonly)) then 2697 n_excl = size(xonly) 2698 allocate(to_excl(1:n_excl)) 2699 to_excl = xonly 2700 2701 else if (present(xalso)) then 2702 n_excl = size(xalso) + size(def_excl) 2703 allocate(to_excl(1:n_excl)) 2704 to_excl(1:size(def_excl)) = def_excl 2705 to_excl(size(def_excl)+1:n_excl) = xalso 2706 2707 else 2708 n_excl = size(def_excl) 2709 allocate(to_excl(1:n_excl)) 2710 to_excl = def_excl 2711 endif 2712 2713 nlheader=size(header) 2714 ! go to end of fortran header 2715 do i = 1, nlheader 2716 if (trim(header(i)) == "") exit 2717 enddo 2718 ! go to top of fits file header 2719 status=0 2720 call ftgrec(unit,0_i4b,record,status) 2721 ! read in all header lines except those excluded 2722 do 2723 call ftgnxk(unit,'*',1_i4b,to_excl,n_excl,record,status) 2724 if (status > 0) exit ! end of header 2725 if (i > nlheader) then 2726 write(unit=*,fmt="(a,i5,a)") & 2727 & " WARNING : The header in "// & 2728 & trim(filename)//" has more than ", & 2729 & nlheader," lines." 2730 print*," It will be truncated." 2731 error = 1 2732 exit 2733 endif 2734 header(i)=record 2735 i=i+1 2736 enddo 2737 status=0 2738 2739 return 2740 end subroutine get_clean_header 2741 2742 2743 !==================================================================== 2744 subroutine check_input_map(code, mapfile, polarisation) 2745 !==================================================================== 2746 ! check out that file contains a valid HEALPIX map 2747 ! on input polarisation is set to T if one wants to get polarisation data from map, 2748 ! and on output it will be set to F if that is not possible 2749 !==================================================================== 2750 use pix_tools, only: nside2npix, npix2nside 2751 character(len=*) :: code 2752 character(len=*) :: mapfile 2753 logical(LGT), intent(inout) :: polarisation 2754 2755 ! 2756 integer(I4B) :: nmaps, order_map, nsmax, mlpol, type, polar_fits, polcconv 2757 integer(I4B) :: npixtot 2758 character(len=1) :: coordsys 2759 character(len=*), parameter :: primer_url = 'http://healpix.sf.net/pdf/intro.pdf' 2760 !==================================================================== 2761 2762 npixtot = getsize_fits(mapfile, nmaps = nmaps, ordering=order_map, nside=nsmax,& 2763 & mlpol=mlpol, type = type, polarisation = polar_fits, & 2764 & coordsys=coordsys, polcconv=polcconv) 2765 2766 if (nsmax<=0) then 2767 print*,"Keyword NSIDE not found in FITS header!" 2768 call fatal_error(code) 2769 endif 2770 if (type == 3) npixtot = nside2npix(nsmax) ! cut sky input data set 2771 if (nsmax/=npix2nside(npixtot)) then 2772 print 9000,"FITS header keyword NSIDE does not correspond" 2773 print 9000,"to the size of the map!" 2774 call fatal_error(code) 2775 endif 2776 2777 if (polarisation .and. (nmaps >=3) .and. polar_fits == -1) then 2778 print 9000,"The input fits file MAY NOT contain polarisation data." 2779 print 9000,"Proceed at your own risk" 2780 endif 2781 2782 if (polarisation .and. (nmaps<3 .or. polar_fits ==0)) then 2783 print 9000,"The file does NOT contain polarisation maps" 2784 print 9000,"only the temperature field will be analyzed" 2785 polarisation = .false. 2786 endif 2787 2788 ! POLCCONV now dealt with in input_map 2789! if (polarisation .and. (polcconv == 3)) then 2790! print 9000,"The polarisation coordinate convention (POLCCONV) is neither COSMO nor IAU." 2791! print 9000,code//" can not proceed with these data" 2792! print 9000,"See Healpix primer ("//primer_url//") for details." 2793! call fatal_error(code) 2794! endif 2795 2796! if (polarisation .and. (polcconv == 2)) then 2797! print 9000,"The input map contains polarized data in the IAU coordinate convention (POLCCONV)" 2798! print 9000,code//" can not proceed with these data" 2799! print 9000,"See Healpix primer ("//primer_url//") for details." 2800! call fatal_error(code) 2801! endif 2802 2803! if (polarisation .and. (polcconv == 0)) then 2804! print 9000,"WARNING: the polarisation coordinate convention (POLCCONV) can not be determined" 2805! print 9000," COSMO will be assumed." 2806! print 9000,"See Healpix primer ("//primer_url//") for details." 2807! endif 2808 2809 ! --- check ordering scheme --- 2810 if ((order_map/=1).and.(order_map/=2)) then 2811 print 9000,"The ordering scheme of the map must be RING or NESTED." 2812 print 9000,"No ordering specification is given in the FITS-header!" 2813 call fatal_error(code) 2814 endif 2815 28169000 format(a) 2817 2818 return 2819 end subroutine check_input_map 2820 2821 2822end module fitstools 2823