New Upstream Release - healpix-fortran
Ready changes
Summary
Merged new upstream version: 3.82+ds (was: 3.81).
Resulting package
Built on 2023-07-04T10:54 (took 6m18s)
The resulting binary packages can be installed (if you have the apt repository enabled) by running one of:
apt install -t fresh-releases libhealpix-devapt install -t fresh-releases libhealpix0-dbgsymapt install -t fresh-releases libhealpix0
Lintian Result
Diff
diff --git a/debian/changelog b/debian/changelog
index 7ae0ae3..1764a66 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,10 @@
+healpix-fortran (3.82+ds-1) UNRELEASED; urgency=low
+
+ * New upstream release.
+ * New upstream release.
+
+ -- Debian Janitor <janitor@jelmer.uk> Tue, 04 Jul 2023 10:49:25 -0000
+
healpix-fortran (3.60+ds-1) unstable; urgency=medium
* Initial release. (Closes: #963092)
diff --git a/debian/patches/0001-Link-against-libsharp-and-libcfitsio.patch b/debian/patches/0001-Link-against-libsharp-and-libcfitsio.patch
index a9ba311..3568b42 100644
--- a/debian/patches/0001-Link-against-libsharp-and-libcfitsio.patch
+++ b/debian/patches/0001-Link-against-libsharp-and-libcfitsio.patch
@@ -6,10 +6,10 @@ Subject: Link against libsharp and libcfitsio
src/f90/mod/Makefile | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
-diff --git a/src/f90/mod/Makefile b/src/f90/mod/Makefile
-index d0a45f3..2344ed7 100644
---- a/src/f90/mod/Makefile
-+++ b/src/f90/mod/Makefile
+Index: healpix-fortran.git/src/f90/mod/Makefile
+===================================================================
+--- healpix-fortran.git.orig/src/f90/mod/Makefile
++++ healpix-fortran.git/src/f90/mod/Makefile
@@ -51,7 +51,7 @@ endif
all: $(OD)_mkdir $(LIB)
diff --git a/src/f90/lib/gifmod.f90 b/src/f90/lib/gifmod.f90
index a30276f..a9feeda 100644
--- a/src/f90/lib/gifmod.f90
+++ b/src/f90/lib/gifmod.f90
@@ -222,7 +222,7 @@ module gifmod
subroutine gifstr(x,nx,ny,sx,sy,nc,r,g,b,or_,fn) bind(C, name="gifstr")
use iso_c_binding
implicit none
- integer (c_int), intent(out) :: x(*)
+ integer (c_int), intent(inout) :: x(*) ! 2020-08-25
integer (c_int), intent(in) :: nx, ny, sx, sy, nc, r(*), g(*), b(*), or_
character(c_char), intent(in) :: fn(*)
end subroutine
@@ -485,8 +485,9 @@ contains
end subroutine addttl
! ------------------------------------------------------------------
- ! addstr: used exclusively by addbar to add labels to the ends of
- ! color bar; [a] is the image array, [i] and [j] are the
+ ! addstr: used by addbar to add labels to the ends of color bar
+ ! and by addttl to add title;
+ ! [a] is the (modified) image array, [i] and [j] are the
! array positions where the string [s] should end if [or]
! is negative, or start if [or] is positive
! ------------------------------------------------------------------
@@ -494,10 +495,11 @@ contains
subroutine addstr(a,i,j,or,s)
use healpix_types
implicit none
- integer(i4b) :: nx,ny,i,j
+ integer(i4b), intent(inout), dimension(:,:) :: a
+ integer(i4b), intent(in) :: i,j
integer(i4b), intent(in) :: or
- integer(i4b), dimension(:,:) :: a
character(len=*), intent(in) :: s
+ integer(i4b) :: nx,ny
nx = size(a,1)
ny = size(a,2)
diff --git a/src/f90/mod/Makefile b/src/f90/mod/Makefile
index d0a45f3..b3985e3 100644
--- a/src/f90/mod/Makefile
+++ b/src/f90/mod/Makefile
@@ -13,7 +13,7 @@ incpix = convert_inplace_1d_inc.f90 \
apply_mask_inc.f90 \
pixel_routines.F90
-incfits = fits_s_inc.f90 fits_d_inc.f90
+incfits = fits_s_inc.F90 fits_d_inc.F90
incudgrade = udgrade_s_inc.f90 udgrade_d_inc.f90
@@ -109,7 +109,7 @@ $(OD)/alm_tools.o:alm_tools.F90 $(incalm)
$(OD)/fitstools.o: fitstools.F90 $(incfits)
-$(OD)/udgrade_nr.o: udgrade_nr.f90 $(incudgrade)
+$(OD)/udgrade_nr.o: udgrade_nr.F90 $(incudgrade)
$(OD)/num_rec.o: num_rec.F90 $(incnumrec)
@@ -123,10 +123,10 @@ alm_map_ss_inc.F90: alm_map_template.F90 gen_alm_code
alm_map_dd_inc.F90: alm_map_template.F90 gen_alm_code alm_map_ss_inc.F90
@./gen_alm_code
-fits_s_inc.f90: fits_template.f90 gen_fits_code
+fits_s_inc.F90: fits_template.F90 gen_fits_code
@./gen_fits_code
-fits_d_inc.f90: fits_template.f90 gen_fits_code fits_s_inc.f90
+fits_d_inc.F90: fits_template.F90 gen_fits_code fits_s_inc.F90
@./gen_fits_code
udgrade_s_inc.f90: udgrade_template.f90 gen_udgrade_code
diff --git a/src/f90/mod/alm_map_template.F90 b/src/f90/mod/alm_map_template.F90
index 3fc2921..fadee7c 100644
--- a/src/f90/mod/alm_map_template.F90
+++ b/src/f90/mod/alm_map_template.F90
@@ -327,7 +327,7 @@
real(DP), dimension(1:2) :: zbounds_in
integer(I4B) :: l, m, ith ! alm related
- integer(I8B) :: istart_south, istart_north, npix ! map related
+ integer(I8B) :: istart_south, istart_north, npix, p! map related
integer(I4B) :: nrings, nphmx
real(DP), dimension(1:4) :: b_even, b_odd
@@ -354,13 +354,27 @@
#ifdef USE_SHARP
zbounds_in = (/-1.d0 , 1.d0/)
if (present(zbounds)) zbounds_in = zbounds
-
+! aspin = abs(spin)
+! if ((aspin>0).and.(aspin<=100)) then
+! call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
+! alm(1:2,0:nlmax,0:nmmax),map(0:12*nsmax*nsmax-1,1:2), zbounds_in)
+! return
+! endif
+ ! 2021-06-15: use scalar routine when s=0
aspin = abs(spin)
- if ((aspin>0).and.(aspin<=100)) then
- call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
- alm(1:2,0:nlmax,0:nmmax),map(0:12*nsmax*nsmax-1,1:2), zbounds_in)
- return
+ npix = (12_I8B*nsmax)*nsmax
+ if (aspin > 0) then
+ call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
+ alm(1:2,0:nlmax,0:nmmax),map(0:npix-1,1:2), zbounds_in)
+ else
+ call sharp_hp_alm2map_x_KLOAD(nsmax,nlmax,nmmax,&
+ alm(1:1,0:nlmax,0:nmmax),map(0:npix-1,1), zbounds_in)
+ do p = 0_I8B, npix - 1_I8B
+ map(p,1) = -map(p,1) ! mimic alm sign convention of alm2map_spin
+ map(p,2) = 0.0_KMAP ! unused for s=0
+ enddo
endif
+ return
#endif
!=======================================================================
@@ -2420,7 +2434,6 @@
if (present(zbounds)) zbounds_in = zbounds
w8ring_in = 1.d0
if (present(w8ring)) w8ring_in = w8ring
-
call sharp_hp_map2alm_x_KLOAD(nsmax,nlmax,nmmax,map,alm,zbounds_in,w8ring_in)
#else
@@ -2696,12 +2709,29 @@
if (present(w8ring)) w8ring_in = w8ring
#ifdef USE_SHARP
+! aspin = abs(spin)
+! if ((aspin>0).and.(aspin<=100)) then
+! call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
+! map(0:12*nsmax*nsmax-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
+! return
+! endif
+ ! 2021-06-15: use scalar routine when s=0
aspin = abs(spin)
- if ((aspin>0).and.(aspin<=100)) then
- call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
- map(0:12*nsmax*nsmax-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
- return
+ npix = (12_I8B*nsmax)*nsmax
+ if (aspin > 0) then
+ call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
+ map(0:npix-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
+ else
+ call sharp_hp_map2alm_x_KLOAD(nsmax,nlmax,nmmax, &
+ map(0:npix-1,1),alm(1:1,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
+ do l=0,nlmax
+ do m=l,nmmax
+ alm(1,l,m) = -alm(1,l,m) ! implement alm sign convention of map2alm_spin
+ alm(2,l,m) = 0.0_DPC ! unused for s=0
+ enddo
+ enddo
endif
+ return
#endif
! Healpix definitions
@@ -5112,33 +5142,41 @@
rms_c1 = ZERO
rms_c2 = ZERO
rms_c3 = ZERO
+ ! compute rms_g2 (E->E) and rms_c2 (E->B)
if (cls_tt(l) > ZERO) then
var_g2 = cls_gg(l) - (cls_tg(l)/cls_tt(l))*cls_tg(l) ! to avoid underflow
- ! test for consistency but make sure it is not due to round off error
- if (var_g2 <= ZERO) then
- if (abs(var_g2) > abs(1.e-8*cls_gg(l))) then
- print*,code//'> Inconsistent TT, GG and TG spectra at l=',l
- call fatal_error
- else ! only round off error, keep going
- var_g2 = ZERO
- endif
+ else
+ var_g2 = cls_gg(l)
+ endif
+ ! test for consistency but make sure it is not due to round off error
+ if (var_g2 <= ZERO) then
+ if (abs(var_g2) > abs(1.e-8*cls_gg(l))) then
+ print*,code//'> Inconsistent TT, GG and TG spectra at l=',l
+ call fatal_error
+ else ! only round off error, keep going
+ var_g2 = ZERO
endif
- rms_c1 = cls_tc(l) / sqrt( cls_tt(l) )
- if (var_g2 > ZERO) then
- rms_g2 = sqrt( var_g2 )
+ endif
+ if (var_g2 > ZERO) then
+ rms_g2 = sqrt( var_g2 )
+ if (cls_tt(l) > ZERO) then
rms_c2 = ( cls_gc(l) - cls_tc(l) * (cls_tg(l) / cls_tt(l)) ) / rms_g2
+ else
+ rms_c2 = cls_gc(l) / rms_g2
endif
- var_c3 = cls_cc(l) - rms_c1**2 - rms_c2**2
- if (var_c3 <= ZERO) then
- if (abs(var_c3) > abs(1.e-8*cls_cc(l))) then
- print*,code//'> Inconsistent spectra at l=',l
- call fatal_error
- else ! only round off error, keep going
- var_c3 = ZERO
- endif
+ endif
+ ! compute rms_c3 (B->B)
+ if (cls_tt(l) > ZERO) rms_c1 = cls_tc(l) / sqrt( cls_tt(l) )
+ var_c3 = cls_cc(l) - rms_c1**2 - rms_c2**2
+ if (var_c3 <= ZERO) then
+ if (abs(var_c3) > abs(1.e-8*cls_cc(l))) then
+ print*,code//'> Inconsistent TT, TG, TC, GC spectra at l=',l
+ call fatal_error
+ else ! only round off error, keep going
+ var_c3 = ZERO
endif
- rms_c3 = sqrt( var_c3 )
endif
+ rms_c3 = sqrt( var_c3 )
! ------ m = 0 ------
zeta2_r = rand_gauss(rng_handle)
diff --git a/src/f90/mod/fits_template.f90 b/src/f90/mod/fits_template.F90
similarity index 76%
rename from src/f90/mod/fits_template.f90
rename to src/f90/mod/fits_template.F90
index 337094a..7817a14 100644
--- a/src/f90/mod/fits_template.f90
+++ b/src/f90/mod/fits_template.F90
@@ -39,6 +39,8 @@
! subroutine write_bintabh_KLOAD (8)
! subroutine unfold_weights_KLOAD
! subroutine unfold_weightslist_KLOAD
+! subroutine read_fits_partial_KLOAD (4/8)
+! subroutine write_fits_partial_KLOAD (4)
!
!
! K M A P : map kind either SP or DP
@@ -175,19 +177,23 @@ subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
LOGICAL(LGT) :: anynull, do_polcconv
! Note : read_fits_cut4 still SP and I4B only
- integer(I4B), dimension(:), allocatable :: pixel
-! real(KMAP), dimension(:), allocatable :: signal
- real(SP), dimension(:), allocatable :: signal
+ integer(I4B), dimension(:), allocatable :: pixel4
+ real(SP), dimension(:), allocatable :: signal
+ ! while read_fits_partial is SP/DP and I4B/I8B
+ integer(I8B), dimension(:), allocatable :: pixel8
+ real(KMAP), dimension(:,:), allocatable :: iqu
integer(I4B) :: status
integer(I4B) :: type_fits, nmaps_fits
CHARACTER(LEN=80) :: units1
CHARACTER(LEN=80),dimension(1:30) :: unitsm
character(len=80), dimension(:), allocatable :: hbuffer
! CHARACTER(LEN=80),dimension(:), allocatable :: unitsm
- integer(I4B) :: extno_i, extno_f
+ integer(I4B) :: extno_i, extno_f, n_ext
CHARACTER(LEN=80) :: polcconv
integer(I4B) :: ipolcconv, polar, nside_fits
character(len=*), parameter :: primer_url = 'http://healpix.sf.net/pdf/intro.pdf'
+ character(len=*), parameter :: code = 'input_map'
+
!-----------------------------------------------------------------------
units1 = ' '
@@ -208,6 +214,7 @@ subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
do_polcconv = .not. ignore_polcconv
endif
+ n_ext = getnumext_fits(filename)
obs_npix = getsize_fits(filename, nmaps = nmaps_fits, type=type_fits, extno=extno_i, &
polcconv=ipolcconv, polarisation=polar, nside=nside_fits)
@@ -228,47 +235,71 @@ subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
units(1:size(units)) = unitsm(1:size(units))
endif
call map_bad_pixels(map, fmissing, fmiss_effct, imissing, verbose=.false.)
-! do imap = 1, nmaps
-! anynull = .true.
-! if (anynull) then
-! imissing = 0
-! do i=0,npixtot-1
-! if ( ABS(map(i,imap)/fmissing -1.) < 1.e-5 ) then
-! map(i,imap) = fmiss_effct
-! imissing = imissing + 1
-! endif
-! enddo
-! endif
-! enddo
if (imissing > 0) write(*,'(a,1pe11.4)') 'blank value : ' ,fmissing
- else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3)) then
- ! cut sky FITS file, reading 1 map (I) or 3 maps (I,Q,U), each from a different extension
+ else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3) &
+ & .and. (nmaps_fits > nmaps) .and. n_ext == 1) then
+ ! partial FITS file, reading 1 map (I) or 3 maps (I,Q,U) from the same extension
+ obs_npix = getsize_fits(filename)
+ allocate(pixel8(0:obs_npix-1), stat = status)
+ call assert_alloc(status, code, "pixel8")
+ allocate(iqu (0:obs_npix-1, 1:nmaps), stat = status)
+ call assert_alloc(status, code, "iqu")
+ call read_fits_partial(filename, pixel8, iqu, header=hbuffer)
+ if (present(header)) then
+ do i=1,nlheader
+ header(i) = hbuffer(i)
+ enddo
+ endif
+ if (present(units)) then
+! do imap=1,min(size(units),nmaps-1) ! 2022-06-10
+ do imap=1,min(size(units),nmaps)
+ call get_card(hbuffer,'TUNIT'//string(imap+1,format='(i0.0)'), units(imap))
+ enddo
+ endif
+ maxpix = maxval(pixel8)
+ minpix = maxval(pixel8)
+ if (maxpix > (npixtot-1) .or. minpix < 0) then
+ print*,'map constructed from file '//trim(filename)//', with pixels in ',minpix,maxpix
+ print*,' wont fit in array with ',npixtot,' elements'
+ call fatal_error
+ endif
+ map(:,:) = fmiss_effct
+ do imap=1,nmaps
+ map(pixel8(:),imap) = iqu(:,imap)
+ enddo
+ imissing = npixtot - obs_npix
+ deallocate(iqu)
+ deallocate(pixel8)
+! else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3) .and. n_ext == nmaps) then
+! ! cut sky FITS file, reading 1 map (I) or 3 maps (I,Q,U), each from a different extension
+ else if (type_fits == 3 .and. (nmaps == 3) .and. n_ext == nmaps) then
+ ! cut sky FITS file, reading 3 maps (I,Q,U), each from a different extension
do imap = 1, nmaps
extno_f = extno_i + imap - 1 ! 2016-08-16
obs_npix = getsize_fits(filename, extno = extno_f)
! one partial map (in binary table with explicit pixel indexing)
- allocate(pixel(0:obs_npix-1), stat = status)
+ allocate(pixel4(0:obs_npix-1), stat = status)
allocate(signal(0:obs_npix-1), stat = status)
call read_fits_cut4(filename, int(obs_npix,kind=i4b), &
- & pixel, signal, header=hbuffer, units=units1, extno=extno_f)
+ & pixel4, signal, header=hbuffer, units=units1, extno=extno_f)
if (present(header) .and. imap == 1) then
do i=1,nlheader
header(i) = hbuffer(i)
enddo
endif
if (present(units)) units(imap) = trim(units1)
- maxpix = maxval(pixel)
- minpix = maxval(pixel)
+ maxpix = maxval(pixel4)
+ minpix = maxval(pixel4)
if (maxpix > (npixtot-1) .or. minpix < 0) then
print*,'map constructed from file '//trim(filename)//', with pixels in ',minpix,maxpix
print*,' wont fit in array with ',npixtot,' elements'
call fatal_error
endif
map(:,imap) = fmiss_effct
- map(pixel(:),imap) = signal(:)
+ map(pixel4(:),imap) = signal(:)
imissing = npixtot - obs_npix
deallocate(signal)
- deallocate(pixel)
+ deallocate(pixel4)
enddo
else
print*,'Unable to read the ',nmaps,' required map(s) from file '//trim(filename)
@@ -2357,3 +2388,580 @@ subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
return
end subroutine unfold_weightslist_KLOAD
+ !=======================================================================
+ ! read_fits_partial(filename, pixel, cutmap, [header, extno])
+ ! routine to read FITS file with cut sky data : PIXEL, ?, ?, ?, ?, ...
+ !=======================================================================
+#ifndef NO64BITS
+ !=======================================================================
+ subroutine read_fits_partial4_KLOAD(filename, pixel, cutmap, header, extno)
+ !=======================================================================
+ use long_intrinsic, only: long_size
+ character(len=*), intent(in) :: filename
+ integer(I4B), dimension(0:), intent(out) :: pixel
+ real(KMAP), dimension(0:,1:), intent(out) :: cutmap
+ character(len=*), dimension(1:), intent(out), optional :: header
+ integer(I4B), intent(in), optional :: extno
+ !
+ integer(i8b), dimension(:), allocatable :: pixel8
+ integer(i8b) :: npix
+ character(len=*), parameter :: code = 'read_fits_partial'
+ !---------------------------------------------------------------------
+ npix = long_size(pixel)
+ if (npix > MAX_I4b) then
+ call fatal_error(code//': PIXEL must be of type I8B to read '//trim(filename))
+ endif
+ allocate(pixel8(0:npix-1))
+ call read_fits_partial8_KLOAD(filename, pixel8, cutmap, header, extno)
+ if (maxval(pixel8) > MAX_I4B) then
+ call fatal_error(code//': PIXEL must be of type I8B to read '//trim(filename))
+ endif
+ pixel(0:npix-1) = pixel8(0:npix-1)
+ deallocate(pixel8)
+ return
+
+ end subroutine read_fits_partial4_KLOAD
+#endif
+ !=======================================================================
+ subroutine read_fits_partial8_KLOAD(filename, pixel, cutmap, header, extno)
+ !=======================================================================
+ use long_intrinsic, only: long_size
+ character(len=*), intent(in) :: filename
+ integer(I8B), dimension(0:), intent(out) :: pixel
+ real(KMAP), dimension(0:,1:), intent(out) :: cutmap
+ character(len=*), dimension(1:), intent(out), optional :: header
+ integer(I4B), intent(in), optional :: extno
+
+ integer(I8B) :: obs_npix, npixtot, npix
+
+ integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension
+ integer(I4B) :: blocksize, datacode
+ integer(I4B) :: firstpix, frow, hdutype
+ integer(I4B) :: naxis, nfound, nmove, nrows, nmaps, nmr, i, npix_4
+ integer(I4B) :: readwrite, status, tfields, unit, varidat, width
+ integer(I4B) :: repeat, repeat1, repeat2
+ logical(LGT) :: anynull, extend
+ character(len=20), dimension(MAXDIM) :: ttype, tform, tunit
+ character(len=20) :: extname
+ character(len=80) :: comment
+ real(KMAP) :: nullval
+ !=======================================================================
+
+ npixtot = min(long_size(pixel),long_size(cutmap,1))
+ nmaps = size(cutmap,2)
+ status=0
+
+ unit = 150
+ nfound = -1
+ anynull = .false.
+
+ readwrite=0
+ call ftopen(unit,filename,readwrite,blocksize,status)
+ if (status > 0) call printerror(status)
+
+ ! determines the presence of image
+ call ftgkyj(unit,'NAXIS', naxis, comment, status)
+ if (status > 0) call printerror(status)
+ if (naxis > 0) then ! there is an image
+ print*,'an image was found in the FITS file '//filename
+ print*,'... it is ignored.'
+ endif
+
+ ! determines the presence of an extension
+ call ftgkyl(unit,'EXTEND', extend, comment, status)
+ if (status > 0) then
+ print*,'extension expected and not found in FITS file '//filename
+ print*,'abort code'
+ call fatal_error
+ endif
+
+ nmove = +1
+ if (present(extno)) nmove = +1 + extno
+ call ftmrhd(unit, nmove, hdutype, status)
+
+ call assert (hdutype==2, 'this is not a binary table')
+
+ ! reads all the FITS related keywords
+ call ftghbn(unit, MAXDIM, &
+ & nrows, tfields, ttype, tform, tunit, extname, varidat, &
+ & status)
+
+ if (.not. (trim(ttype(1)) == 'PIXEL' &
+ & .or. trim(ttype(1)) == 'PIX' ) ) call fatal_error('did not find PIXEL column in '//filename)
+
+
+ if (present(header)) then
+ header = ""
+ status = 0
+ !call fitstools_mp_get_clean_header(unit, header, filename, status)
+ call get_clean_header(unit, header, filename, status)
+ endif
+
+! if (present(units)) then
+! ! the second column contains the SIGNAL, for which we
+! ! already read the units from the header
+! units = adjustl(tunit(2))
+! endif
+
+ ! finds the bad data value
+ call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
+ if (status == 202) then ! bad_data not found
+ if (KMAP == SP) nullval = s_bad_value ! default value
+ if (KMAP == DP) nullval = d_bad_value ! default value
+ status = 0
+ endif
+
+ frow = 1
+ firstpix = 1
+ ! parse TFORM keyword to find out the length of the column vector
+ repeat1 = 1
+ repeat2 = 1
+ call ftbnfm(tform(1), datacode, repeat1, width, status)
+ if (tfields > 1) call ftbnfm(tform(2), datacode, repeat2, width, status)
+ repeat = max(repeat1, repeat2)
+
+ !call ftgkyj(unit,'OBS_NPIX',obs_npix,comment,status) ! I4B
+ !call ftgkyk(unit,'OBS_NPIX',obs_npix,comment,status) ! I8B
+ call f90ftgky_(unit,'OBS_NPIX',obs_npix,comment,status) ! generic
+ if (status == 202) then ! obs_npix not found
+ obs_npix = int(nrows, kind=i8b) * repeat
+ status = 0
+ print*,'Keyword OBS_NPIX not found in '//trim(filename)
+ print*,'Will assume OBS_NPIX='//trim(string(obs_npix,format='(i0.0)'))//' instead.' ! necessary to do a printout with gfortran (on M1), 2022-07-06
+ endif
+
+ npix = min(npixtot, obs_npix)
+ npix_4 = npix
+ nmr = min(nmaps, tfields-1)
+ !call ftgcvj (unit, 1_i4b, frow, firstpix, npix_4, i_bad_value, pixel(0), anynull, status) ! I4B
+ !call ftgcvk (unit, 1_i4b, frow, firstpix, npix_4, l_bad_value, pixel(0), anynull, status) ! I8B
+ call f90ftgcv_(unit, 1_i4b, frow, firstpix, npix_4, l_bad_value, pixel, anynull, status) ! generic
+ do i=1,nmr
+ call f90ftgcv_(unit, i+1_i4b, frow, firstpix, npix_4, nullval, cutmap(0:npix-1,i), anynull, status)
+ enddo
+
+ ! close the file
+ call ftclos(unit, status)
+
+ ! check for any error, and if so print out error messages
+ if (status > 0) call printerror(status)
+
+ return
+ end subroutine read_fits_partial8_KLOAD
+
+
+ !=======================================================================
+ ! writes a fits file for (polarized) sky data set with columns:
+ ! PIXEL, TEMPERATURE, [Q_POLARISATION, U_POLARISATION]
+ !
+ ! write_fits_partial(filename, pixel, cutmap, &
+ ! [, header, coord, nside, order, units, extno])
+ !=======================================================================
+ subroutine write_fits_partial4_KLOAD(filename, pixel, cutmap, &
+ & header, coord, nside, order, units, extno)
+ !=======================================================================
+ character(len=*), intent(in) :: filename
+ integer(I4B), dimension(0:), intent(in) :: pixel
+ real(KMAP), dimension(0:,1:), intent(in) :: cutmap
+ character(len=*), dimension(1:), intent(in), optional :: header
+ integer(I4B), intent(in), optional :: nside, order
+ character(len=*), intent(in), optional :: coord, units
+ integer(I4B), intent(in), optional :: extno !, polarisation
+ character(len=*), parameter :: routine = 'write_fits_partial'
+ call sub_write_fits_partial_KLOAD(filename, cutmap, pixel4=pixel, &
+ header=header, coord=coord, nside=nside, order=order, units=units, extno=extno)
+ end subroutine write_fits_partial4_KLOAD
+ !=======================================================================
+#ifndef NO64BITS
+ subroutine write_fits_partial8_KLOAD(filename, pixel, cutmap, &
+ & header, coord, nside, order, units, extno)
+ !=======================================================================
+ character(len=*), intent(in) :: filename
+ integer(I8B), dimension(0:), intent(in) :: pixel
+ real(KMAP), dimension(0:,1:), intent(in) :: cutmap
+ character(len=*), dimension(1:), intent(in), optional :: header
+ integer(I4B), intent(in), optional :: nside, order
+ character(len=*), intent(in), optional :: coord, units
+ integer(I4B), intent(in), optional :: extno !, polarisation
+ character(len=*), parameter :: routine = 'write_fits_partial'
+ call sub_write_fits_partial_KLOAD(filename, cutmap, pixel8=pixel, &
+ header=header, coord=coord, nside=nside, order=order, units=units, extno=extno)
+ end subroutine write_fits_partial8_KLOAD
+#endif
+ !=======================================================================
+ subroutine sub_write_fits_partial_KLOAD(filename, cutmap, &
+ & pixel4, pixel8, &
+ & header, coord, nside, order, units, extno)
+ !=======================================================================
+ use long_intrinsic, only: long_size
+ use pix_tools, only: nside2npix
+ use misc_utils, only: assert
+
+ character(len=*), intent(in) :: filename
+ real(KMAP), dimension(0:,1:), intent(in) :: cutmap
+ integer(I4B), dimension(0:), intent(in), optional :: pixel4
+ integer(I8B), dimension(0:), intent(in), optional :: pixel8
+ character(len=*), dimension(1:), intent(in), optional :: header
+ integer(I4B), intent(in), optional :: nside, order
+ character(len=*), intent(in), optional :: coord, units
+ integer(I4B), intent(in), optional :: extno !, polarisation
+
+ character(len=*), parameter :: routine = 'write_fits_partial'
+ ! --- healpix/fits related variables ----
+ integer(I4B) :: ncol, grain, nd2
+ integer(I4B) :: npix_hd, nside_hd, i
+ integer(I4B) :: maxpix, minpix
+ character(len=1) :: char1, coord_usr
+! character(len=8) :: char8
+ character(len=20) :: units_usr
+ logical(LGT) :: done_nside, done_order, done_coord, done_polar, polar_flag
+ integer(I4B) :: nlheader, extno_i
+
+ ! --- cfitsio related variables ----
+ integer(I4B) :: status, unit, blocksize, bitpix, naxis, naxes(1)
+ logical(LGT) :: simple, extend
+ character(LEN=80) :: svalue, comment
+
+ integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension
+ integer(I4B) :: nrows, tfields, varidat
+ integer(I4B) :: frow, felem, repeat, repeatg, hdutype
+ character(len=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
+ character(len=20), dimension(1:3) :: extnames
+ character(len=80) :: card
+ character(len=4) :: srepeat, srepeatg
+ character(len=1) :: pform1, pform2
+ character(len=filenamelen) sfilename
+ integer(I8B) :: obs_npix8, npix_final, obs_npp
+ integer(I4B) :: obs_npix4, kpix, nside_final
+
+ integer(I4B), save :: nside_old
+ !=======================================================================
+ if (KMAP == SP) pform2='E'
+ if (KMAP == DP) pform2='D'
+ if (present(pixel4) .eqv. present(pixel8)) then
+ call fatal_error(routine//': choose either PIXEL4 or PIXEL8')
+ endif
+ if (present(pixel4)) then
+ kpix = I4B
+ pform1='J'
+ obs_npp = long_size(pixel4)
+ elseif (present(pixel8)) then
+ kpix = I8B
+ pform1='K'
+ obs_npp = long_size(pixel8)
+ endif
+ obs_npix8 = long_size(cutmap, 1)
+ nd2 = size(cutmap, 2)
+ obs_npix4 = int(obs_npix8, kind=I4B)
+ call assert( obs_npix8 == obs_npp, routine//': mismatched size for PIXEL and DATA')
+ ncol = 1 + nd2
+ grain = 1
+ status=0
+ unit = 149
+ blocksize=1
+
+ nlheader = 0
+ if (present(header)) nlheader = size(header)
+ units_usr = ' '
+ if (present(units)) units_usr = units
+ extno_i = 0
+ if (present(extno)) extno_i = extno
+ polar_flag = .false.
+ done_polar = .false.
+! if (present(polarisation)) then
+! polar_flag = (polarisation /= 0)
+! done_polar = .true. ! will ignore the header provided value of POLAR
+! endif
+
+ if (extno_i == 0) then
+ !*************************************
+ ! create the new empty FITS file
+ !*************************************
+ call ftinit(unit,filename,blocksize,status)
+ if (status > 0) call fatal_error("Error while creating file " &
+ & //trim(filename) &
+ & //". Check path and/or access rights.")
+
+ ! -----------------------------------------------------
+ ! initialize parameters about the FITS image
+ simple=.true.
+ bitpix=32 ! integer*4
+ naxis=0 ! no image
+ naxes(1)=0
+ extend=.true. ! there is an extension
+
+ ! ----------------------
+ ! primary header
+ ! ----------------------
+ ! write the required header keywords
+ call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
+
+ ! write the current date
+ call ftpdat(unit,status) ! format ccyy-mm-dd
+
+ else
+ !*********************************************
+ ! reopen an existing file and go to the end
+ !*********************************************
+ ! remove the leading '!' (if any) when reopening the same file
+ sfilename = adjustl(filename)
+ if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
+ call ftopen(unit,sfilename,1_i4b,blocksize,status)
+ ! move to last extension written
+ call ftmahd(unit, 1_i4b+ extno_i, hdutype, status)
+
+ endif
+
+
+ ! ----------------------
+ ! new extension
+ ! ----------------------
+ call ftcrhd(unit, status)
+
+ ! writes required keywords
+ ! repeat = 1024
+ repeat = 1
+ if (obs_npix8 < repeat) repeat = 1
+ nrows = (obs_npix8 + repeat - 1)/ repeat ! naxis1
+ write(srepeat,'(i4)') repeat
+ srepeat = adjustl(srepeat)
+
+ repeatg = repeat * grain
+ write(srepeatg,'(i4)') repeatg
+ srepeatg = adjustl(srepeatg)
+
+ tfields = ncol
+ ttype(1) = 'PIXEL'
+ if (nd2 == 1) then
+ done_polar = .true.
+ polar_flag = .false.
+ ttype(2) = 'TEMPERATURE '
+ elseif (nd2 == 3) then
+ done_polar = .true.
+ polar_flag = .true.
+ ttype(2:4) = (/'TEMPERATURE ' , &
+ & 'Q_POLARISATION' , &
+ & 'U_POLARISATION' /)
+ else
+ do i=2,ncol
+ ttype(i) = 'C'//string(i-1, format='(i2.2)') ! C01, C02, ...
+ enddo
+ endif
+
+ tform(1) = trim(srepeat) //pform1
+ tform(2:ncol) = trim(srepeatg)//pform2
+
+ tunit = ' ' ! optional, will not appear
+ tunit(2:ncol) = units_usr
+
+ extname = 'SKY_OBSERVATION' ! default, will be overridden by user provided one if any
+ ! if (polar_flag) extname = extnames(1+extno_i)
+ varidat = 0
+ call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
+ & extname, varidat, status)
+
+ call ftpcom(unit,'------------------------------------------',status)
+ call ftpcom(unit,' Pixelisation Specific Keywords ',status)
+ call ftpcom(unit,'------------------------------------------',status)
+ call ftpkys(unit,'PIXTYPE','HEALPIX ',' HEALPIX Pixelisation',status)
+ call ftpkyu(unit,'NSIDE', ' ',status) ! place holder, will be updated later on
+ call ftpkyu(unit,'ORDERING',' ',status)
+ call ftpkys(unit,'COORDSYS',' ',' ',status)
+ call ftpcom(unit,' G = Galactic, E = ecliptic, C = celestial = equatorial', status)
+ call ftpkyl(unit,'POLAR',polar_flag,'Polarisation included in file (T/F)',status)
+ call ftpcom(unit,'------------------------------------------',status)
+ call ftpcom(unit,' Data Specific Keywords ',status)
+ call ftpcom(unit,'------------------------------------------',status)
+ call ftpkys(unit,'INDXSCHM','EXPLICIT',' Indexing : IMPLICIT or EXPLICIT', status)
+! call ftpkyj(unit,'GRAIN', grain, ' Grain of pixel indexing',status)
+! call ftpcom(unit,'GRAIN=0 : no indexing of pixel data (IMPLICIT)',status)
+! call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)',status)
+! call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
+ call ftpkys(unit,'OBJECT','PARTIAL ',' Sky coverage represented by data',status)
+ if (KPIX == I4B) then
+ call ftpkyj(unit,'OBS_NPIX',obs_npix4, ' Number of pixels observed and recorded',status)
+ else
+ call ftpkyk(unit,'OBS_NPIX',obs_npix8, ' Number of pixels observed and recorded',status)
+ endif
+
+ ! add required Healpix keywords (NSIDE, ORDER) if provided by user
+ done_order = .false.
+ if (present(order)) then
+ ! char8 = order
+ ! call ftupch(char8)
+ ! if (char8(1:4) == 'RING') then
+ if (order == 1) then
+ call ftukys(unit, 'ORDERING','RING',' Pixel ordering scheme, either RING or NESTED',status)
+ done_order = .true.
+ ! elseif (char8(1:4) == 'NEST') then
+ elseif (order == 2) then
+ call ftukys(unit, 'ORDERING','NESTED',' Pixel ordering scheme, either RING or NESTED ',status)
+ done_order = .true.
+ else
+ print*,'Invalid ORDER given : ',order, ' instead of 1 (RING) or 2 (NESTED)'
+ endif
+ endif
+
+ done_nside = .false.
+ if (present(nside)) then
+ if (nside2npix(nside) > 0) then ! valid nside
+ call ftukyj(unit,'NSIDE',nside,' Resolution parameter for HEALPIX', status)
+ done_nside = .true.
+ nside_final = nside
+ else
+ print*,'Invalid NSIDE given : ',nside
+ endif
+ endif
+
+ ! add non required Healpix keyword (COORD)
+ done_coord = .false.
+ if (present(coord)) then
+ coord_usr = adjustl(coord)
+ char1 = coord_usr(1:1)
+ call ftupch(char1) ! uppercase
+ if (char1 == 'C' .or. char1 == 'Q') then
+ coord_usr = 'C'
+ done_coord = .true.
+ elseif (char1 == 'G') then
+ coord_usr = 'G'
+ done_coord = .true.
+ elseif (char1 == 'E' ) then
+ coord_usr='E'
+ done_coord = .true.
+ else
+ print*,'Unrecognised COORD given : ',coord,' instead of C, G, or E'
+ print*,'Proceed at you own risks '
+ coord_usr = char1
+ done_coord = .true.
+ endif
+ if (done_coord) then
+ call ftukys(unit, 'COORDSYS',coord_usr,' Pixelisation coordinate system',status)
+ endif
+ endif
+
+
+ ! write the user provided header literally, except for PIXTYPE, TFORM*, TTYPE*, TUNIT*, INDXSCHM and GRAIN
+ ! copy NSIDE, ORDERING and COORDSYS and POLAR if they are valid and not already given
+ do i=1,nlheader
+ card = header(i)
+ if (card(1:5) == 'TTYPE' .or. card(1:5) == 'TFORM' .or. card(1:7) == 'PIXTYPE') then
+ continue ! don't keep them
+ else if (card(1:8) == 'INDXSCHM') then
+ continue
+ else if (card(1:5) == 'GRAIN') then ! already written above
+ continue
+ else if (card(1:13) == 'COMMENT GRAIN' .or. card(1:14) == 'COMMENT GRAIN') then ! already written above
+ continue
+ else if (card(1:5) == 'TUNIT') then
+ if (trim(units_usr) == '') then
+ call ftucrd(unit,'TUNIT'//card(6:7),card, status) !update TUNIT2 and above
+ endif
+ else if (card(1:5) == 'NSIDE') then
+ call ftpsvc(card, svalue, comment, status)
+ read(svalue,*) nside_hd
+ npix_hd = nside2npix(nside_hd)
+ if (.not. done_nside .and. npix_hd > 0) then
+ call ftucrd(unit,'NSIDE',card, status) !update NSIDE
+ done_nside = .true.
+ nside_final = nside_hd
+ endif
+ else if (card(1:8) == 'ORDERING') then
+ call ftpsvc(card, svalue, comment, status)
+ svalue = adjustl(svalue)
+ svalue = svalue(2:8) ! remove leading quote
+ call ftupch(svalue)
+ if (.not. done_order .and. (svalue(1:4)=='RING' .or. svalue(1:6) == 'NESTED')) then
+ call ftucrd(unit,'ORDERING',card, status) !update ORDERING
+ done_order = .true.
+ endif
+ else if (card(1:8) == 'COORDSYS') then
+ if (.not. done_coord) call putrec(unit,card, status)
+ done_coord = .true.
+ else if (card(1:5) == 'POLAR') then
+ if (.not. done_polar) then
+ call ftucrd(unit,'POLAR', card, status) ! update POLAR
+ done_polar = .true.
+ endif
+ else if (card(1:7) == 'EXTNAME') then
+ call ftucrd(unit, 'EXTNAME', card, status)
+ else if (card(1:1) /= ' ') then ! if none of the above, copy to FITS file
+ call putrec(unit, card, status)
+ endif
+10 continue
+ enddo
+
+
+ ! test that required keywords have been provided in some way
+ if (.not. done_nside) then
+ print*,routine//'> NSIDE is a Healpix required keyword, '
+ print*,routine//'> it was NOT provided either as routine argument or in the input header'
+ print*,routine//'> abort execution, file not written'
+ call fatal_error
+ endif
+ if (.not. done_order) then
+ print*,routine//'> ORDER is a Healpix required keyword, '
+ print*,routine//'> it was NOT provided either as routine argument or in the input header'
+ print*,routine//'> abort execution, file not written'
+ call fatal_error
+ endif
+! if ((.not. done_polar) .and. extno_i >=2) then
+! print*,routine//'> Warning: POLAR keyword not set while 3 extensions have been written'
+! endif
+
+ ! check that NSIDE is the same for all extensions
+ if (extno_i == 0) then
+ nside_old = nside_final
+ else
+ if (nside_final /= nside_old) then
+ print*,routine//'> Inconsistent NSIDE: ',nside_final, nside_old
+ print*,routine//'> Should use same NSIDE for all extensions'
+ endif
+ endif
+
+ ! check validity of PIXEL
+ npix_final = nside2npix(nside_final)
+ if (kpix == I4B) then
+ minpix = minval(pixel4(0:obs_npix8-1))
+ maxpix = maxval(pixel4(0:obs_npix8-1))
+ else
+ minpix = minval(pixel8(0:obs_npix8-1))
+ maxpix = maxval(pixel8(0:obs_npix8-1))
+ endif
+ if (minpix < 0 .or. maxpix > npix_final-1) then
+ print*,routine//'> Actual pixel range = ',minpix,maxpix
+ print*,routine//'> expected range (for Nside =',nside_final,') : 0, ',npix_final-1
+ print*,routine//'> ABORT execution, file not written.'
+ call fatal_error
+ endif
+ if (obs_npix8 > npix_final) then
+ print*,routine//'> The actual number of pixels ',obs_npix8
+ print*,routine//'> is larger than the number of pixels over the whole sky : ',npix_final
+ print*,routine//'> for Nside = ',nside_final
+ print*,routine//'> ABORT execution, file not written.'
+ call fatal_error
+ endif
+
+ ! write the extension one column by one column
+ frow = 1 ! starting position (row)
+ felem = 1 ! starting position (element)
+ if (kpix == I4B) then
+ call f90ftpcl_(unit, 1_i4b, frow, felem, obs_npix4, pixel4, status)
+ else
+ call f90ftpcl_(unit, 1_i4b, frow, felem, obs_npix4, pixel8, status)
+ endif
+ do i=1, nd2
+ call f90ftpcl_(unit, i+1_i4b, frow, felem, obs_npix4, cutmap(0:,i), status)
+ enddo
+
+ ! ----------------------
+ ! close and exit
+ ! ----------------------
+
+ ! close the file and free the unit number
+ call ftclos(unit, status)
+
+ ! check for any error, and if so print out error messages
+ if (status > 0) call printerror(status)
+
+ return
+ end subroutine sub_write_fits_partial_KLOAD
diff --git a/src/f90/mod/fitstools.F90 b/src/f90/mod/fitstools.F90
index 2f3d3ce..6705185 100644
--- a/src/f90/mod/fitstools.F90
+++ b/src/f90/mod/fitstools.F90
@@ -84,6 +84,8 @@ module fitstools
! subroutine read_bintod
! subroutine write_bintabh
! subroutine unfold_weights
+ ! subroutine read_fits_partial
+ ! subroutine write_fits_partial
! ----------------------------------
!
! subroutine read_fits_cut4 ?
@@ -119,6 +121,7 @@ module fitstools
real(kind=SP), private, parameter :: s_bad_value = HPX_SBADVAL
real(kind=DP), private, parameter :: d_bad_value = HPX_DBADVAL
integer(kind=I4B), private, parameter :: i_bad_value = -1637500000
+ integer(kind=I8B), private, parameter :: l_bad_value = -1637500000_I8B
integer(I4B) , private, parameter :: nchunk_max = 12000
integer(I4B), private, parameter :: MAXDIM_TOP = 199 ! < 999
@@ -167,6 +170,24 @@ module fitstools
#endif
end interface
+ interface read_fits_partial
+#ifdef NO64BITS
+ module procedure read_fits_partial8_s, read_fits_partial8_d
+#else
+ module procedure read_fits_partial8_s, read_fits_partial8_d, &
+ & read_fits_partial4_s, read_fits_partial4_d
+#endif
+ end interface
+
+ interface write_fits_partial
+#ifdef NO64BITS
+ module procedure write_fits_partial4_s, write_fits_partial4_d
+#else
+ module procedure write_fits_partial8_s, write_fits_partial8_d, &
+ & write_fits_partial4_s, write_fits_partial4_d
+#endif
+ end interface
+
!
interface input_tod
module procedure input_tod_s,input_tod_d
@@ -217,11 +238,19 @@ module fitstools
end interface
interface f90ftpcl_
- module procedure f90ftpcle, f90ftpcld
+#ifdef NO64BITS
+ module procedure f90ftpcle, f90ftpcld, f90ftpclj
+#else
+ module procedure f90ftpcle, f90ftpcld, f90ftpclj, f90ftpclk
+#endif
end interface
interface f90ftgcv_
- module procedure f90ftgcve, f90ftgcvd
+#ifdef NO64BITS
+ module procedure f90ftgcve, f90ftgcvd, f90ftgcvj
+#else
+ module procedure f90ftgcve, f90ftgcvd, f90ftgcvj, f90ftgcvk
+#endif
end interface
interface f90ftgpv_
@@ -229,7 +258,11 @@ module fitstools
end interface
interface f90ftgky_
- module procedure f90ftgkye, f90ftgkyd
+#ifdef NO64BITS
+ module procedure f90ftgkye, f90ftgkyd, f90ftgkyj
+#else
+ module procedure f90ftgkye, f90ftgkyd, f90ftgkyj, f90ftgkyk
+#endif
end interface
interface map_bad_pixels
@@ -249,6 +282,7 @@ module fitstools
public :: read_fits_cut4, write_fits_cut4, &
& input_map, read_bintab, &
& output_map, write_bintab
+ public :: read_fits_partial, write_fits_partial
public :: write_plm
public :: fits2cl, read_asctab, write_asctab
public :: read_dbintab, write_dbintab
@@ -262,7 +296,7 @@ module fitstools
contains
!-------------------------------------------------------------------------------
- ! generic interface F90FTPCL_ for FITSIO's FTPCLE and FTPCLD
+ ! generic interface F90FTPCL_ for FITSIO's FTPCL[E,D,J,K]
! writes data in ASCTAB or BINTAB
subroutine f90ftpcle(unit, colnum, frow, felem, np, data, status)
integer(I4B), intent(in) :: unit, colnum, frow, felem, np
@@ -278,8 +312,24 @@ contains
call ftpcld(unit, colnum, frow, felem, np, data, status)
return
end subroutine f90ftpcld
+ subroutine f90ftpclj(unit, colnum, frow, felem, np, data, status)
+ integer(I4B), intent(in) :: unit, colnum, frow, felem, np
+ integer(I4B), intent(out) :: status
+ integer(I4B), intent(in), dimension(0:) :: data
+ call ftpclj(unit, colnum, frow, felem, np, data, status)
+ return
+ end subroutine f90ftpclj
+#ifndef NO64BITS
+ subroutine f90ftpclk(unit, colnum, frow, felem, np, data, status)
+ integer(I4B), intent(in) :: unit, colnum, frow, felem, np
+ integer(I4B), intent(out) :: status
+ integer(I8B), intent(in), dimension(0:) :: data
+ call ftpclk(unit, colnum, frow, felem, np, data, status)
+ return
+ end subroutine f90ftpclk
+#endif
!-------------------------------------------------------------------------------
- ! generic interface F90FTGCV_ for FITSIO's FTGCVE and FTGCVD
+ ! generic interface F90FTGCV_ for FITSIO's FTGCV[E,D,J,K]
! reads data from BINTAB
subroutine f90ftgcve(unit, colnum, frow, felem, np, nullval, data, anynull, status)
integer(I4B), intent(in) :: unit, colnum, frow, felem, np
@@ -299,6 +349,26 @@ contains
call ftgcvd(unit, colnum, frow, felem, np, nullval, data, anynull, status)
return
end subroutine f90ftgcvd
+ subroutine f90ftgcvj(unit, colnum, frow, felem, np, nullval, data, anynull, status)
+ integer(I4B), intent(in) :: unit, colnum, frow, felem, np
+ integer(I4B), intent(out) :: status
+ logical(LGT), intent(out) :: anynull
+ integer(I4B), intent(out), dimension(0:) :: data
+ integer(I4B), intent(in) :: nullval
+ call ftgcvj(unit, colnum, frow, felem, np, nullval, data, anynull, status)
+ return
+ end subroutine f90ftgcvj
+#ifndef NO64BITS
+ subroutine f90ftgcvk(unit, colnum, frow, felem, np, nullval, data, anynull, status)
+ integer(I4B), intent(in) :: unit, colnum, frow, felem, np
+ integer(I4B), intent(out) :: status
+ logical(LGT), intent(out) :: anynull
+ integer(I8B), intent(out), dimension(0:) :: data
+ integer(I8B), intent(in) :: nullval
+ call ftgcvk(unit, colnum, frow, felem, np, nullval, data, anynull, status)
+ return
+ end subroutine f90ftgcvk
+#endif
!-------------------------------------------------------------------------------
! generic interface F90FTGPV_ for FITSIO's FTGPVE and FTGPVD
! reads data from IMAGE
@@ -321,7 +391,7 @@ contains
return
end subroutine f90ftgpvd
!-------------------------------------------------------------------------------
- ! generic interface F90FTGKY_ for FITSIO's FTGKYE and FTGKYD
+ ! generic interface F90FTGKY_ for FITSIO's FTGKY[E,D,J,K]
! reads a keyword
subroutine f90ftgkye(unit, keyword, value, comment, status)
integer(I4B), intent(in) :: unit
@@ -341,14 +411,34 @@ contains
call ftgkyd(unit, keyword, value, comment, status)
return
end subroutine f90ftgkyd
+ subroutine f90ftgkyj(unit, keyword, value, comment, status)
+ integer(I4B), intent(in) :: unit
+ character(len=*), intent(in) :: keyword
+ integer(I4B), intent(out) :: status
+ character(len=*), intent(out) :: comment
+ integer(I4B), intent(out) :: value
+ call ftgkyj(unit, keyword, value, comment, status)
+ return
+ end subroutine f90ftgkyj
+#ifndef NO64BITS
+ subroutine f90ftgkyk(unit, keyword, value, comment, status)
+ integer(I4B), intent(in) :: unit
+ character(len=*), intent(in) :: keyword
+ integer(I4B), intent(out) :: status
+ character(len=*), intent(out) :: comment
+ integer(I8B), intent(out) :: value
+ call ftgkyk(unit, keyword, value, comment, status)
+ return
+ end subroutine f90ftgkyk
+#endif
!-------------------------------------------------------------------------------
! define routine with SP I/O
-#include "fits_s_inc.f90"
+#include "fits_s_inc.F90"
! define routine with DP I/O
-#include "fits_d_inc.f90"
+#include "fits_d_inc.F90"
@@ -402,14 +492,14 @@ contains
call ftgkyj(unit,'NAXIS', naxis, comment, status)
if (status > 0) call printerror(status)
if (naxis > 0) then ! there is an image
- print*,'an image was found in the FITS file '//filename
+ print*,'an image was found in the FITS file '//trim(filename)
print*,'... it is ignored.'
endif
! determines the presence of an extension
call ftgkyl(unit,'EXTEND', extend, comment, status)
if (status > 0) then
- print*,'extension expected and not found in FITS file '//filename
+ print*,'extension expected and not found in FITS file '//trim(filename)
print*,'abort code'
call fatal_error
endif
@@ -425,7 +515,7 @@ contains
& nrows, tfields, ttype, tform, tunit, extname, varidat, &
& status)
if (tfields < 4) then
- print*,'Expected 4 columns in FITS file '//filename
+ print*,'Expected 4 columns in FITS file '//trim(filename)
print*,'found ',tfields
if (tfields < 2) call fatal_error
if (.not. (trim(ttype(1)) == 'PIXEL' &
@@ -614,8 +704,8 @@ contains
! writes required keywords
! repeat = 1024
repeat = 1
- nrows = (obs_npix + repeat - 1)/ repeat ! naxis1
if (obs_npix < repeat) repeat = 1
+ nrows = (obs_npix + repeat - 1)/ repeat ! naxis1
write(srepeat,'(i4)') repeat
srepeat = adjustl(srepeat)
@@ -633,7 +723,7 @@ contains
tunit = ' ' ! optional, will not appear
tunit(2) = units_usr
tunit(4) = units_usr
- extname = 'SKY_OBSERVATION' ! default, will be overide by user provided one if any
+ extname = 'SKY_OBSERVATION' ! default, will be overridden by user provided one if any
if (polar_flag) extname = extnames(1+extno_i)
varidat = 0
call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
@@ -652,10 +742,10 @@ contains
call ftpcom(unit,' Data Specific Keywords ',status)
call ftpcom(unit,'------------------------------------------',status)
call ftpkys(unit,'INDXSCHM','EXPLICIT',' Indexing : IMPLICIT or EXPLICIT', status)
- call ftpkyj(unit,'GRAIN', grain, ' Grain of pixel indexing',status)
- call ftpcom(unit,'GRAIN=0 : no indexing of pixel data (IMPLICIT) ',status)
- call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)',status)
- call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
+! call ftpkyj(unit,'GRAIN', grain, ' Grain of pixel indexing',status)
+! call ftpcom(unit,'GRAIN=0 : no indexing of pixel data (IMPLICIT)',status)
+! call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)',status)
+! call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
call ftpkys(unit,'OBJECT','PARTIAL ',' Sky coverage represented by data',status)
call ftpkyj(unit,'OBS_NPIX',obs_npix, ' Number of pixels observed and recorded',status)
@@ -716,7 +806,7 @@ contains
endif
- ! write the user provided header literally, except for PIXTYPE, TFORM*, TTYPE*, TUNIT* and INDXSCHM
+ ! write the user provided header literally, except for PIXTYPE, TFORM*, TTYPE*, TUNIT*, INDXSCHM and GRAIN
! copy NSIDE, ORDERING and COORDSYS and POLAR if they are valid and not already given
do i=1,nlheader
card = header(i)
@@ -724,6 +814,10 @@ contains
continue ! don't keep them
else if (card(1:8) == 'INDXSCHM') then
continue
+ else if (card(1:5) == 'GRAIN') then ! already written above
+ continue
+ else if (card(1:13) == 'COMMENT GRAIN' .or. card(1:14) == 'COMMENT GRAIN') then ! already written above
+ continue
else if (card(1:5) == 'TUNIT') then
if ((card(6:6) == '2' .or. card(6:6) == '4') .and. trim(units_usr) == '') then
call ftucrd(unit,'TUNIT'//card(6:6),card, status) !update TUNIT2 and TUNIT4
@@ -1968,7 +2062,8 @@ contains
INTEGER(I4B) :: status,unit,readwrite,blocksize,naxis
CHARACTER(LEN=80) :: comment, ttype1
LOGICAL(LGT) :: extend, anyf
- INTEGER(I4B):: nmove, hdutype, idmax, nrows
+ INTEGER(I4B):: nmove, hdutype, nrows
+ integer(I4B), dimension(1) :: idmax ! because ftgcvj expects an array (2020-08-24)
!-----------------------------------------------------------------------
status=0
@@ -2015,9 +2110,9 @@ contains
ttype1 = trim(strupcase(adjustl(ttype1)))
if (trim(ttype1(1:5)) == 'INDEX') then
call ftgkyj(unit, 'NAXIS2', nrows, comment, status) ! find number of rows
- call ftgcvj(unit, 1_i4b, nrows, 1_i4b, 1_i4b, 0_i4b, idmax, anyf, status) ! read element on last row of first column
+ call ftgcvj(unit, 1_i4b, nrows, 1_i4b, 1_i4b, 0_i4b, idmax(1), anyf, status) ! read element on last row of first column
if (status == 0) then
- lmax = int(sqrt( real(idmax-1, kind = DP) ) )
+ lmax = int(sqrt( real(idmax(1)-1, kind = DP) ) )
if (lmax > 0) goto 1000
endif
endif
@@ -2759,28 +2854,29 @@ contains
character(len=*), parameter :: primer_url = 'http://healpix.sf.net/pdf/intro.pdf'
!====================================================================
+ !n_ext = getnumext_fits(mapfile)
npixtot = getsize_fits(mapfile, nmaps = nmaps, ordering=order_map, nside=nsmax,&
& mlpol=mlpol, type = type, polarisation = polar_fits, &
& coordsys=coordsys, polcconv=polcconv)
if (nsmax<=0) then
- print*,"Keyword NSIDE not found in FITS header!"
+ print*,"Keyword NSIDE not found in FITS header of "//trim(mapfile)//" !"
call fatal_error(code)
endif
if (type == 3) npixtot = nside2npix(nsmax) ! cut sky input data set
if (nsmax/=npix2nside(npixtot)) then
print 9000,"FITS header keyword NSIDE does not correspond"
- print 9000,"to the size of the map!"
+ print 9000,"to the size of the map in "//trim(mapfile)//" !"
call fatal_error(code)
endif
if (polarisation .and. (nmaps >=3) .and. polar_fits == -1) then
- print 9000,"The input fits file MAY NOT contain polarisation data."
+ print 9000,"The input FITS file "//trim(mapfile)//" MAY NOT contain polarisation data."
print 9000,"Proceed at your own risk"
endif
if (polarisation .and. (nmaps<3 .or. polar_fits ==0)) then
- print 9000,"The file does NOT contain polarisation maps"
+ print 9000,"The FITS file "//trim(mapfile)//" does NOT contain polarisation maps"
print 9000,"only the temperature field will be analyzed"
polarisation = .false.
endif
@@ -2809,7 +2905,7 @@ contains
! --- check ordering scheme ---
if ((order_map/=1).and.(order_map/=2)) then
print 9000,"The ordering scheme of the map must be RING or NESTED."
- print 9000,"No ordering specification is given in the FITS-header!"
+ print 9000,"No ordering specification is given in the FITS-header of "//trim(mapfile)//" !"
call fatal_error(code)
endif
diff --git a/src/f90/mod/gen_fits_code b/src/f90/mod/gen_fits_code
index a17ab38..bc5fed4 100755
--- a/src/f90/mod/gen_fits_code
+++ b/src/f90/mod/gen_fits_code
@@ -3,9 +3,9 @@
# EH, Dec 2004 -- Jan 2005
SetDef(){
-template_file=fits_template.f90
-outfile_s=fits_s_inc.f90
-outfile_d=fits_d_inc.f90
+template_file=fits_template.F90
+outfile_s=fits_s_inc.F90
+outfile_d=fits_d_inc.F90
}
#--------
diff --git a/src/f90/mod/head_fits.F90 b/src/f90/mod/head_fits.F90
index 477ff6f..b804850 100644
--- a/src/f90/mod/head_fits.F90
+++ b/src/f90/mod/head_fits.F90
@@ -524,7 +524,7 @@ contains
endif
endif
enddo
- if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+ if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
if (present(comment)) comment = ' '
if (present(count)) count = count_in
return
@@ -563,7 +563,7 @@ contains
endif
endif
enddo
- if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+ if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
if (present(comment)) comment = ' '
if (present(count)) count = count_in
return
@@ -673,7 +673,7 @@ contains
return
endif
enddo
- if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+ if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
if (present(comment)) comment = ' '
if (present(count)) count = count_in
return
@@ -712,7 +712,7 @@ contains
endif
endif
enddo
- if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+ if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
if (present(comment)) comment = ' '
if (present(count)) count = count_in
return
@@ -753,7 +753,7 @@ contains
endif
endif
enddo
- if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+ if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
if (present(comment)) comment = ' '
if (present(count)) count = count_in
return
@@ -792,7 +792,7 @@ contains
endif
endif
enddo
- if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
+ if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
if (present(comment)) comment = ' '
if (present(count)) count = count_in
return
@@ -878,14 +878,15 @@ contains
character(len=80), dimension(1:20) :: my_units
character(len=8) :: my_ordering
character(len=FILENAMELEN) :: udtype, my_beam_leg
- logical(LGT) :: do_full, do_alm, do_cl, do_cut, do_polar, &
- & do_bcross, do_reset, do_asym
+ logical(LGT) :: do_full, do_alm, do_cl, do_cut, do_partial, &
+ & do_polar, do_bcross, do_reset, do_asym, do_map
character(len=*), parameter :: code = 'write_minimal_header'
!--------------------------------------------------------------------------------
! find out type of data to be contained in FITS file
udtype = strupcase(dtype)
- do_alm = .false. ; do_cl = .false. ; do_full = .false. ; do_cut = .false.
+ do_alm = .false. ; do_cl = .false. ; do_full = .false.
+ do_cut = .false. ; do_partial = .false.
select case (trim(udtype))
case ('ALM')
do_alm = .true.
@@ -895,11 +896,14 @@ contains
do_full = .true.
case ('CUTMAP')
do_cut = .true.
+ ! case ('PARTIAL')
+ do_partial = .true.
case default
print*,'Invalid choice of datatype:'//trim(dtype)
print*,'Should be one of: ALM, alm, CL, cl, CUTMAP, cutmap, MAP, map'
call fatal_error(code)
end select
+ do_map = (do_cut .or. do_full .or. do_partial)
!
if (present(order) .and. present(ordering)) then
@@ -908,7 +912,7 @@ contains
endif
! check for required keywords for map
- if (do_cut .or. do_full) then
+ if (do_map) then
if (.not. present(nside)) then
print*,'NSIDE required for maps in '//code
call fatal_error(code)
@@ -1108,7 +1112,7 @@ contains
endif ! cl
- if (do_full .or. do_cut) then
+ if (do_map) then
!--------------------------------------------------------------------------------------------
! minimal header for FULL SKY MAP and CUT SKY MAP
!--------------------------------------------------------------------------------------------
@@ -1156,27 +1160,54 @@ contains
call add_card(header,"COMMENT","Full sky data")
call add_card(header,"OBJECT","FULLSKY")
call add_card(header,"INDXSCHM","IMPLICIT"," Indexing : IMPLICIT or EXPLICIT")
- call add_card(header,"GRAIN", 0, " Grain of pixel indexing")
+! call add_card(header,"GRAIN", 0, " Grain of pixel indexing")
endif
- if (do_cut) then
+ if (do_cut .or. do_partial) then
call add_card(header,"EXTNAME","'CUT SKY MAP'",update=.true.)
call add_card(header,"COMMENT","Cut sky data")
call add_card(header,"OBJECT","PARTIAL")
call add_card(header,"INDXSCHM","EXPLICIT"," Indexing : IMPLICIT or EXPLICIT")
- call add_card(header,"GRAIN", 1, " Grain of pixel indexing")
+! call add_card(header,"GRAIN", 1, " Grain of pixel indexing")
endif
- call add_card(header,"COMMENT","GRAIN=0 : no indexing of pixel data (IMPLICIT)")
- call add_card(header,"COMMENT","GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)")
- call add_card(header,"COMMENT","GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)")
+! call add_card(header,"COMMENT","GRAIN=0 : no indexing of pixel data (IMPLICIT)")
+! call add_card(header,"COMMENT","GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)")
+! call add_card(header,"COMMENT","GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)")
call add_card(header) ! blank line
call add_card(header,"POLAR",do_polar," Polarisation included (True/False)")
if (do_cut) then
+ call add_card(header) ! blank line
+ call add_card(header,"TTYPE1", "PIXEL","pixel index")
+ call add_card(header,"TUNIT1", "", "")
+ call add_card(header)
if (present(units)) then
my_units = units
- call add_card(header,"TUNIT2", my_units(1),"physical unit of signal map")
- call add_card(header,"TUNIT4", my_units(1),"physical unit of error map")
+ call add_card(header,"TUNIT2", my_units(1),"physical units of map")
+ call add_card(header,"TUNIT4", my_units(1),"physical units of map")
endif
endif
+! if (do_partial) then
+! if (present(units)) my_units = units
+! it = 1 ; iq = 1 ; iu = 1
+! call add_card(header) ! blank line
+! call add_card(header,"TTYPE1", "PIXEL","pixel index")
+! call add_card(header,"TUNIT1", "", "")
+! call add_card(header)
+
+! call add_card(header) ! blank line
+! call add_card(header,"TTYPE2", "TEMPERATURE","Temperature map")
+! call add_card(header,"TUNIT2", my_units(it), "physical units of map")
+! call add_card(header)
+
+! if (do_polar) then
+! call add_card(header,"TTYPE3", "Q_POLARISATION","Q Polarisation map")
+! call add_card(header,"TUNIT3", my_units(iq), "physical units of map")
+! call add_card(header)
+
+! call add_card(header,"TTYPE4", "U_POLARISATION","U Polarisation map")
+! call add_card(header,"TUNIT4", my_units(iu), "physical units unit")
+! call add_card(header)
+! endif
+! endif
if (do_full) then
call add_card(header,"DERIV",my_deriv," Derivative included (0, 1 or 2)")
@@ -1199,12 +1230,10 @@ contains
call add_card(header)
if (do_polar) then
- !call add_card(header,"TTYPE2", "Q-POLARISATION","Q Polarisation map")
call add_card(header,"TTYPE2", "Q_POLARISATION","Q Polarisation map")
call add_card(header,"TUNIT2", my_units(iq),"map unit")
call add_card(header)
- !call add_card(header,"TTYPE3", "U-POLARISATION","U Polarisation map")
call add_card(header,"TTYPE3", "U_POLARISATION","U Polarisation map")
call add_card(header,"TUNIT3", my_units(iu),"map unit")
call add_card(header)
diff --git a/src/f90/mod/healpix_types.F90 b/src/f90/mod/healpix_types.F90
index fc4dff8..f9cd184 100644
--- a/src/f90/mod/healpix_types.F90
+++ b/src/f90/mod/healpix_types.F90
@@ -39,7 +39,7 @@ MODULE healpix_types
! Mars 2008: i8b same as i4b on machines not supporting 64 bits (NO64BITS flag set)
! Feb 2009: introduce healpix_version
!
- character(len=*), PARAMETER, public :: healpix_version = '3.60'
+ character(len=*), PARAMETER, public :: healpix_version = '3.82'
INTEGER, PARAMETER, public :: i4b = SELECTED_INT_KIND(9)
#ifdef NO64BITS
INTEGER, PARAMETER, public :: i8b = i4b
diff --git a/src/f90/mod/pix_tools.F90 b/src/f90/mod/pix_tools.F90
index f9ec2d0..0df526a 100644
--- a/src/f90/mod/pix_tools.F90
+++ b/src/f90/mod/pix_tools.F90
@@ -871,9 +871,11 @@ contains
integer(i4b), dimension(1:,1:), intent(inout) :: ringphi
integer(i4b), intent(inout) :: ngr
- integer(i4b) :: i, j, k, kk, ngr_out, diff, iphi, i0, nrh
+ !integer(i4b) :: i, j, k, kk, ngr_out, diff, iphi, i0, nrh
+ integer(i4b) :: i, j, k, kk, ngr_out, iphi, i0, nrh, nshrink, np
real(dp), dimension(1:2*nsboost+1) :: phiw, phie
real(dp) :: dd, dph, phic
+ ! 2020-08-31: bug correction
!=======================================================================
if (nsboost <= 1) return
@@ -882,11 +884,13 @@ contains
do i=1, ngr ! loop on low-res rings
i0 = ringphi(1,i) * nsboost - nsboost - irmin
do k=-1,1,2 ! West and East side of disc
+ nshrink = 0
kk = (k+5)/2 ! 2 or 3
222 continue
iphi = ringphi(kk, i)
- if (ringphi(2,i) <= ringphi(3,i) .and. iphi >= 0) then
- call find_pixel_bounds(nside, nsboost, ringphi(1,i), iphi, phiw, phie)
+ !if (ringphi(2,i) <= ringphi(3,i) .and. iphi >= 0) then
+ if (iphi >= 0) then ! corrected 2020-08-31
+ call find_pixel_bounds(nside, nsboost, ringphi(1,i), iphi, phiw, phie, np=np)
do j=1, 2*nsboost+1
if (i0+j >= 1 .and. i0+j <= nrh) then
phic = (phie(j)+phiw(j))*0.5_dp ! pixel center
@@ -896,18 +900,25 @@ contains
if (dd <= dph) goto 1000 ! pixel touched by disc, move to next one
endif
enddo
- ringphi(kk, i)= iphi - k ! pixel not in disc, move edge pixel inwards
- goto 222 ! try next pixel inward
+ !ringphi(kk, i)= iphi - k ! pixel not in disc, move edge pixel inwards
+ ringphi(kk, i) = mod(iphi + np - k, np) ! pixel not in disc, move edge pixel inwards in [0,np-1]
+ nshrink = nshrink + 1
+ if (nshrink <= 2) then
+ goto 222 ! try next pixel inward
+ else
+ ringphi(kk,i) = -999
+ endif
1000 continue
endif
- enddo ! loop on side
+ enddo ! loop on k East-West sides
enddo ! loop on low-res rings
! remove empty rings
ngr_out = 0
do i=1,ngr
- diff = ringphi(3,i) - ringphi(2,i)
- if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0 .and. diff /= -2 .and. diff /= -1) then
+! diff = ringphi(3,i) - ringphi(2,i)
+! if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0 .and. diff /= -2 .and. diff /= -1) then
+ if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0) then ! 2020-08-31
ngr_out = ngr_out + 1
ringphi(1:3, ngr_out) = ringphi(1:3, i)
endif
@@ -921,10 +932,11 @@ contains
return
end subroutine check_edge_pixels
!=======================================================================
- subroutine find_pixel_bounds (nside, nsboost, iring, iphi, phiw, phie)
+ subroutine find_pixel_bounds (nside, nsboost, iring, iphi, phiw, phie, np)
!=======================================================================
integer(i4b), intent(in) :: nside, nsboost, iring, iphi
real(dp), dimension(1:2*nsboost+1), intent(out) :: phiw, phie
+ integer(i4b), optional, intent(out) :: np
real(dp), dimension(1:2*nsboost+1) :: f, f1, phiw_t, phie_t
real(dp) :: c0, quad, phie1, phie2, phiw1, phiw2, cv
@@ -936,6 +948,7 @@ contains
!f = ((/ (i,i=0,2*nsboost) /) - nsboost) / nsboost
f = ((/ (i*1.d0,i=0,2*nsboost) /) - nsboost*1.d0) / nsboost
+ if (present(np)) np = npr
nq = npr/4 ! number of pixels on current ring in [0,Pi/2] (quadrant)
transition = (iring == nside .or. iring == nside*3)
diff --git a/src/f90/mod/pixel_routines.F90 b/src/f90/mod/pixel_routines.F90
index 26a6b4b..b5d7fac 100644
--- a/src/f90/mod/pixel_routines.F90
+++ b/src/f90/mod/pixel_routines.F90
@@ -1920,6 +1920,7 @@
! New algorithm for Inclusive case: test boundary of edge pixels on each ring
! 2013-04-02: bug correction in query_disc in inclusive mode
! 2014-07-07: bug correction in discedge2fulldisc
+! 2020-08-31: bug correction in check_edge_pixels (reported by R. de Belsunce)
!=======================================================================
#ifdef DOI8B
subroutine query_disc_8( nside, vector0, radius, listpix, nlist, nest, inclusive)
@@ -2016,8 +2017,9 @@
call pixels_on_edge(nside, irmin, irmax, phi0, dphitab, ringphi, ngr)
if (do_inclusive) then
! sample edge pixels at larger Nside
- nsboost = 16
- nsideh = min(NS_MAX8, nside * int(nsboost,i8b))
+ nsboost = 2**ceiling( log(1.1_DP/(radius*nside))/log(2.0_DP) ) ! make pixel smaller than disc (2020-08-31)
+ nsboost = max(nsboost,16)
+ nsideh = min(NS_MAX8, nside * int(nsboost,i8b))
radiush = fudge_query_radius(nsideh, radius, quadratic=.true.)
irmin = ring_num(nsideh, zmax, shift=+1) ! shifted South
@@ -2028,6 +2030,7 @@
zlist(iz-irmin+1) = ring2z(nsideh, iz)
enddo
call discphirange_at_z(vector0, radiush, zlist, nrh, dphilist, phi0)
+ phi0 = mod(phi0 + TWOPI, TWOPI) ! [-Pi,Pi] -> [0,2Pi]
call check_edge_pixels(nside, nsboost, irmin, irmax, phi0, dphilist, ringphi, ngr)
deallocate(zlist, dphilist)
endif
diff --git a/src/f90/mod/sharp_healpix_f_inc.c b/src/f90/mod/sharp_healpix_f_inc.c
index ff7b876..901efb7 100644
--- a/src/f90/mod/sharp_healpix_f_inc.c
+++ b/src/f90/mod/sharp_healpix_f_inc.c
@@ -44,6 +44,7 @@ void X(map2alm_spin) (int nside, int lmax, int mmax, int spin,
void *mapptr[2], *almptr[2];
CHECK_STACK_ALIGN(8);
+ UTIL_ASSERT(spin>0, "spin must not be zero");
sharp_make_healpix_geom_info_2 (nside, wgt, zbounds[0], zbounds[1], &ginfo);
sharp_make_rectangular_alm_info (lmax,mmax,2,&ainfo);
@@ -109,6 +110,7 @@ void X(alm2map_spin) (int nside, int lmax, int mmax, int spin,
void *mapptr[2], *almptr[2];
CHECK_STACK_ALIGN(8);
+ UTIL_ASSERT(spin>0, "spin must not be zero");
// sharp_make_healpix_geom_info (nside, 1, &ginfo);
int i,nrings=4*nside-1;
double *wgt=RALLOC(double,nrings);
diff --git a/src/f90/mod/udgrade_nr.f90 b/src/f90/mod/udgrade_nr.F90
similarity index 97%
rename from src/f90/mod/udgrade_nr.f90
rename to src/f90/mod/udgrade_nr.F90
index b0d6178..757115d 100644
--- a/src/f90/mod/udgrade_nr.f90
+++ b/src/f90/mod/udgrade_nr.F90
@@ -59,10 +59,10 @@ module udgrade_nr
contains
! include SP implementation of routines
- include 'udgrade_s_inc.f90'
+#include "udgrade_s_inc.f90"
! include DP implementation of routines
- include 'udgrade_d_inc.f90'
+#include "udgrade_d_inc.f90"
end module udgrade_nr
diff --git a/src/f90/synfast/syn_sub_inc.f90 b/src/f90/synfast/syn_sub_inc.f90
index f777355..72db731 100644
--- a/src/f90/synfast/syn_sub_inc.f90
+++ b/src/f90/synfast/syn_sub_inc.f90
@@ -121,7 +121,7 @@
character(len=100) :: chline
type(planck_rng) :: rng_handle
- integer(i4b) :: deriv, n_pols, n_maps
+ integer(i4b) :: deriv, n_pols, n_maps, n_clin
character(len=20) :: coordsys = ' '
real(DP), dimension(1:2) :: zbounds
! real :: ttime0, ttime1, ttime2
@@ -397,6 +397,8 @@
if (input_cl) then
polar_cl = polar ! 0 or 1: input 1 or 4 field spectrum
!polar_cl = 2 * polar ! 0 or 2; input 1 or 6 field spectrum
+ !junk = getsize_fits(infile, nmaps=n_clin)
+ !if (n_clin == 6) polar_cl = 2*polar ! input 1 or 6 spectra
call create_alm(nsmax, nlmax, nmmax, polar_cl, &
& infile, rng_handle, fwhm_arcmin, &
& alm_TGC, header_PS, windowfile, units=units_alm, beam_file=beam_file)
Debdiff
[The following lists of changes regard files as different if they have different names, permissions or owners.]
Files in second set of .debs but not in first
-rw-r--r-- root/root /usr/lib/debug/.build-id/58/06212d3f53c18752a1f4420f21d150ec670a50.debug
Files in first set of .debs but not in second
-rw-r--r-- root/root /usr/lib/debug/.build-id/02/6d1a16ac8e14c87e79f9f89ca3c44a97499090.debug
Control files of package libhealpix-dev: lines which differ (wdiff format)
Depends: gfortran-11 gfortran-13 | gfortran-mod-15, libhealpix0 (=
No differences were encountered between the control files of package libhealpix0
Control files of package libhealpix0-dbgsym: lines which differ (wdiff format)
Build-Ids: 026d1a16ac8e14c87e79f9f89ca3c44a97499090 5806212d3f53c18752a1f4420f21d150ec670a50