Commit f71e89ba56acb602abd783125b07188bad96aa51 - healpix-fortran

Import upstream version 3.81 Debian Janitor 2 years ago

15 changed file(s) with 3293 addition(s) and 2513 deletion(s). Raw diff Collapse all Expand all

src/._healpy less more

Binary diff not shown

-5

src/f90/lib/gifmod.f90 less more

221	221	subroutine gifstr(x,nx,ny,sx,sy,nc,r,g,b,or_,fn) bind(C, name="gifstr")
222	222	use iso_c_binding
223	223	implicit none
224		integer (c_int), intent(out) :: x(*)
	224	integer (c_int), intent(inout) :: x(*) ! 2020-08-25
225	225	integer (c_int), intent(in) :: nx, ny, sx, sy, nc, r(), g(), b(*), or_
226	226	character(c_char), intent(in) :: fn(*)
227	227	end subroutine

484	484	end subroutine addttl
485	485
486	486	! ------------------------------------------------------------------
487		! addstr: used exclusively by addbar to add labels to the ends of
488		! color bar; [a] is the image array, [i] and [j] are the
	487	! addstr: used by addbar to add labels to the ends of color bar
	488	! and by addttl to add title;
	489	! [a] is the (modified) image array, [i] and [j] are the
489	490	! array positions where the string [s] should end if [or]
490	491	! is negative, or start if [or] is positive
491	492	! ------------------------------------------------------------------

493	494	subroutine addstr(a,i,j,or,s)
494	495	use healpix_types
495	496	implicit none
496		integer(i4b) :: nx,ny,i,j
	497	integer(i4b), intent(inout), dimension(:,:) :: a
	498	integer(i4b), intent(in) :: i,j
497	499	integer(i4b), intent(in) :: or
498		integer(i4b), dimension(:,:) :: a
499	500	character(len=*), intent(in) :: s
	501	integer(i4b) :: nx,ny
500	502
501	503	nx = size(a,1)
502	504	ny = size(a,2)

-4

src/f90/mod/Makefile less more

12	12	apply_mask_inc.f90 \
13	13	pixel_routines.F90
14	14
15		incfits = fits_s_inc.f90 fits_d_inc.f90
	15	incfits = fits_s_inc.F90 fits_d_inc.F90
16	16
17	17	incudgrade = udgrade_s_inc.f90 udgrade_d_inc.f90
18	18

108	108
109	109	$(OD)/fitstools.o: fitstools.F90 $(incfits)
110	110
111		$(OD)/udgrade_nr.o: udgrade_nr.f90 $(incudgrade)
	111	$(OD)/udgrade_nr.o: udgrade_nr.F90 $(incudgrade)
112	112
113	113	$(OD)/num_rec.o: num_rec.F90 $(incnumrec)
114	114

122	122	alm_map_dd_inc.F90: alm_map_template.F90 gen_alm_code alm_map_ss_inc.F90
123	123	@./gen_alm_code
124	124
125		fits_s_inc.f90: fits_template.f90 gen_fits_code
	125	fits_s_inc.F90: fits_template.F90 gen_fits_code
126	126	@./gen_fits_code
127	127
128		fits_d_inc.f90: fits_template.f90 gen_fits_code fits_s_inc.f90
	128	fits_d_inc.F90: fits_template.F90 gen_fits_code fits_s_inc.F90
129	129	@./gen_fits_code
130	130
131	131	udgrade_s_inc.f90: udgrade_template.f90 gen_udgrade_code

+43

-13

src/f90/mod/alm_map_template.F90 less more

326	326
327	327	real(DP), dimension(1:2) :: zbounds_in
328	328	integer(I4B) :: l, m, ith ! alm related
329		integer(I8B) :: istart_south, istart_north, npix ! map related
	329	integer(I8B) :: istart_south, istart_north, npix, p! map related
330	330	integer(I4B) :: nrings, nphmx
331	331
332	332	real(DP), dimension(1:4) :: b_even, b_odd

353	353	#ifdef USE_SHARP
354	354	zbounds_in = (/-1.d0 , 1.d0/)
355	355	if (present(zbounds)) zbounds_in = zbounds
356
	356	! aspin = abs(spin)
	357	! if ((aspin>0).and.(aspin<=100)) then
	358	! call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
	359	! alm(1:2,0:nlmax,0:nmmax),map(0:12nsmaxnsmax-1,1:2), zbounds_in)
	360	! return
	361	! endif
	362	! 2021-06-15: use scalar routine when s=0
357	363	aspin = abs(spin)
358		if ((aspin>0).and.(aspin<=100)) then
359		call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
360		alm(1:2,0:nlmax,0:nmmax),map(0:12nsmaxnsmax-1,1:2), zbounds_in)
361		return
362		endif
	364	npix = (12_I8Bnsmax)nsmax
	365	if (aspin > 0) then
	366	call sharp_hp_alm2map_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
	367	alm(1:2,0:nlmax,0:nmmax),map(0:npix-1,1:2), zbounds_in)
	368	else
	369	call sharp_hp_alm2map_x_KLOAD(nsmax,nlmax,nmmax,&
	370	alm(1:1,0:nlmax,0:nmmax),map(0:npix-1,1), zbounds_in)
	371	do p = 0_I8B, npix - 1_I8B
	372	map(p,1) = -map(p,1) ! mimic alm sign convention of alm2map_spin
	373	map(p,2) = 0.0_KMAP ! unused for s=0
	374	enddo
	375	endif
	376	return
363	377	#endif
364	378	!=======================================================================
365	379

2419	2433	if (present(zbounds)) zbounds_in = zbounds
2420	2434	w8ring_in = 1.d0
2421	2435	if (present(w8ring)) w8ring_in = w8ring
2422
2423	2436	call sharp_hp_map2alm_x_KLOAD(nsmax,nlmax,nmmax,map,alm,zbounds_in,w8ring_in)
2424	2437	#else
2425	2438

2695	2708	if (present(w8ring)) w8ring_in = w8ring
2696	2709
2697	2710	#ifdef USE_SHARP
	2711	! aspin = abs(spin)
	2712	! if ((aspin>0).and.(aspin<=100)) then
	2713	! call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
	2714	! map(0:12nsmaxnsmax-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
	2715	! return
	2716	! endif
	2717	! 2021-06-15: use scalar routine when s=0
2698	2718	aspin = abs(spin)
2699		if ((aspin>0).and.(aspin<=100)) then
2700		call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
2701		map(0:12nsmaxnsmax-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
2702		return
2703		endif
	2719	npix = (12_I8Bnsmax)nsmax
	2720	if (aspin > 0) then
	2721	call sharp_hp_map2alm_spin_x_KLOAD(nsmax,nlmax,nmmax,aspin, &
	2722	map(0:npix-1,1:2),alm(1:2,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
	2723	else
	2724	call sharp_hp_map2alm_x_KLOAD(nsmax,nlmax,nmmax, &
	2725	map(0:npix-1,1),alm(1:1,0:nlmax,0:nmmax),zbounds_in,w8ring_in)
	2726	do l=0,nlmax
	2727	do m=l,nmmax
	2728	alm(1,l,m) = -alm(1,l,m) ! implement alm sign convention of map2alm_spin
	2729	alm(2,l,m) = 0.0_DPC ! unused for s=0
	2730	enddo
	2731	enddo
	2732	endif
	2733	return
2704	2734	#endif
2705	2735
2706	2736	! Healpix definitions

+2964

-0

src/f90/mod/fits_template.F90 less more

	0	!-----------------------------------------------------------------------------
	1	!
	2	! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon,
	3	! Benjamin D. Wandelt, Anthony J. Banday,
	4	! Matthias Bartelmann, Hans K. Eriksen,
	5	! Frode K. Hansen, Martin Reinecke
	6	!
	7	!
	8	! This file is part of HEALPix.
	9	!
	10	! HEALPix is free software; you can redistribute it and/or modify
	11	! it under the terms of the GNU General Public License as published by
	12	! the Free Software Foundation; either version 2 of the License, or
	13	! (at your option) any later version.
	14	!
	15	! HEALPix is distributed in the hope that it will be useful,
	16	! but WITHOUT ANY WARRANTY; without even the implied warranty of
	17	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	18	! GNU General Public License for more details.
	19	!
	20	! You should have received a copy of the GNU General Public License
	21	! along with HEALPix; if not, write to the Free Software
	22	! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
	23	!
	24	! For more information about HEALPix see http://healpix.sourceforge.net
	25	!
	26	!-----------------------------------------------------------------------------
	27	! template for routine SP/DP overloading for module fitstools
	28	!
	29	! subroutine input_map_KLOAD (4/8)
	30	! subroutine read_bintab_KLOAD (4/8)
	31	! subroutine read_conbintab_KLOAD NotYet
	32	! subroutine write_bintab_KLOAD (4/8)
	33	! subroutine write_asctab_KLOAD NA
	34	! subroutine dump_alms_KLOAD NotYet
	35	! subroutine write_alms_KLOAD NotYet
	36	! subroutine read_alms_KLOAD NotYet
	37	! subroutine read_bintod_KLOAD (8)
	38	! subroutine write_bintabh_KLOAD (8)
	39	! subroutine unfold_weights_KLOAD
	40	! subroutine unfold_weightslist_KLOAD
	41	! subroutine read_fits_partial_KLOAD (4/8)
	42	! subroutine write_fits_partial_KLOAD (4)
	43	!
	44	!
	45	! K M A P : map kind either SP or DP
	46	!
	47	! K L O A D : suffixe of routine name, to be replaced by either s (SP) or d (DP)
	48	!
	49	! edited Jan 11, 2006 to deal correctly with polarised alms in write_alms and alms2fits
	50	! edited Apr 04, 2006 to avoid concatenation problem in TFORMs (write_plm and write_bintabh with Ifort 9)
	51	! edited Dec 20, 2006 to accept alm file with less elements than expected (in particular if vanishing ones are not included)
	52	! 2007-05-10 : increased maxdim (max number of columns per extension) from 20 to 40
	53	! 2008-08-27 : in dump_alms and write_alms and write_tab:
	54	! do not write TTYPE# and TFORM# in excess of # of fields in the file
	55	! 2008-10-14: corrected bug introduced in write_asctab
	56	! 2012-02-23: correction of a possible bug with index writing in dump_alms and write_alms
	57	! 2013-12-13: increased MAXDIM from 40 to MAXDIM_TOP
	58	! 2018-05-22: added unfold_weights, unfold_weightslist
	59	! 2019-10-14: can write TTYPE??? (ie with up to 3 digits) in write_* and dump_alms
	60	!
	61
	62	!=======================================================================
	63	! map_bad_pixels(map, fin, fout, nbads, verbose)
	64	! map: input data (2D)
	65	! fin: input value of 'bad' pixels
	66	! fout: output value of same 'bad' pixels
	67	! nbads: number of bad pixels found
	68	! verbose: OPTIONAL, if set, be verbose
	69	!=======================================================================
	70
	71	subroutine map_bad_pixels_KLOAD(map, fin, fout, nbads, verbose)
	72	use long_intrinsic, only: long_size
	73	real(KMAP), dimension(0:,1:), intent(inout) :: map
	74	real(KMAP), intent(in) :: fin, fout
	75	integer(I8B), intent(out) :: nbads
	76	logical(LGT), optional, intent(in) :: verbose
	77	integer(I8B) :: i, npix
	78	integer(I8B), dimension(1:100) :: imiss
	79	integer(I4B) :: imap, nmaps
	80	logical(LGT) :: be_verbose
	81	real(KMAP) :: threshold
	82	!-----------------------------------------
	83
	84	npix = long_size(map, 1)
	85	nmaps = long_size(map, 2)
	86	threshold = abs(fin * 1.e-5_KMAP)
	87
	88	imiss(:) = 0
	89	do imap = 1, nmaps
	90	do i=0,npix-1
	91	if ( abs( map(i,imap) - fin ) < threshold ) then
	92	map(i,imap) = fout
	93	imiss(imap) = imiss(imap)+1
	94	endif
	95	enddo
	96	enddo
	97	nbads = sum(imiss)
	98
	99	be_verbose = .false.
	100	if (present(verbose)) be_verbose=verbose
	101	if (be_verbose) then
	102	write(*,'(a,1pe11.4)') 'blank value : ' ,fin
	103	do imap = 1, nmaps
	104	if (imiss(imap) > 0) then
	105	write(*,'(i12,a,f7.3,a,1pe11.4)') &
	106	& imiss(imap),' missing pixels (', &
	107	& (100.0_KMAP*imiss(imap))/npix,' %),'// &
	108	& ' have been set to : ',fout
	109	endif
	110	enddo
	111	endif
	112
	113	return
	114	end subroutine map_bad_pixels_KLOAD
	115	!=======================================================================
	116	! input_map
	117	! reads fits file
	118	! filename = fits file (input)
	119	! map = data read from the file (ouput) = real*4 array of size (npixtot,nmaps)
	120	! npixtot = length of the map (input)
	121	! nmaps = number of maps
	122	! fmissval = OPTIONAL argument (input) with the value to be given to missing
	123	! pixels, its default value is 0
	124	! header = OPTIONAL (output) contains extension header
	125	! units = OPTIONAL (output) contains the maps units
	126	! extno = OPTIONAL (input) contains the unit number to read from (0 based)
	127	! ignore_polcconv = OPTIONAL (input), LGT, default=.false.
	128	! take into account or not the POLCCONV FITS keyword
	129	!
	130	! modified Feb 03 for units argument to run with Sun compiler
	131	! 2017-02-15: deals with POLCCONV
	132	!=======================================================================
	133	#ifndef NO64BITS
	134	subroutine input_map4_KLOAD(filename, map, npixtot, nmaps, &
	135	& fmissval, header, units, extno, ignore_polcconv)
	136	!=======================================================================
	137
	138	INTEGER(I4B), INTENT(IN) :: npixtot
	139	INTEGER(I4B), INTENT(IN) :: nmaps
	140	REAL(KMAP), INTENT(OUT), dimension(0:,1:) :: map
	141	CHARACTER(LEN=*), INTENT(IN) :: filename
	142	REAL(KMAP), INTENT(IN), OPTIONAL :: fmissval
	143	CHARACTER(LEN=*), INTENT(OUT), dimension(1:), OPTIONAL :: header
	144	CHARACTER(LEN=*), INTENT(OUT), dimension(1:), OPTIONAL :: units
	145	INTEGER(I4B), INTENT(IN) , optional :: extno
	146	logical(LGT), intent(IN) , optional :: ignore_polcconv
	147	integer(i8b) :: npixtot8
	148
	149	npixtot8 = npixtot
	150	call input_map8_KLOAD(filename, map, npixtot8, nmaps, &
	151	fmissval, header, units, extno, ignore_polcconv)
	152
	153	return
	154	end subroutine input_map4_KLOAD
	155	#endif
	156
	157	!=======================================================================
	158	subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
	159	fmissval, header, units, extno, ignore_polcconv)
	160	!=======================================================================
	161	use head_fits, only: get_card, add_card
	162	use pix_tools, only: nside2npix
	163	INTEGER(I8B), INTENT(IN) :: npixtot
	164	INTEGER(I4B), INTENT(IN) :: nmaps
	165	REAL(KMAP), INTENT(OUT), dimension(0:,1:) :: map
	166	CHARACTER(LEN=*), INTENT(IN) :: filename
	167	REAL(KMAP), INTENT(IN), OPTIONAL :: fmissval
	168	CHARACTER(LEN=*), INTENT(OUT), dimension(1:), OPTIONAL :: header
	169	CHARACTER(LEN=*), INTENT(OUT), dimension(1:), OPTIONAL :: units
	170	INTEGER(I4B), INTENT(IN) , optional :: extno
	171	logical(LGT), intent(IN) , optional :: ignore_polcconv
	172
	173	INTEGER(I8B) :: i, imissing, obs_npix, maxpix, minpix
	174	REAL(KMAP) :: fmissing, fmiss_effct
	175	integer(I4B) :: imap, nlheader
	176
	177	LOGICAL(LGT) :: anynull, do_polcconv
	178	! Note : read_fits_cut4 still SP and I4B only
	179	integer(I4B), dimension(:), allocatable :: pixel4
	180	real(SP), dimension(:), allocatable :: signal
	181	! while read_fits_partial is SP/DP and I4B/I8B
	182	integer(I8B), dimension(:), allocatable :: pixel8
	183	real(KMAP), dimension(:,:), allocatable :: iqu
	184	integer(I4B) :: status
	185	integer(I4B) :: type_fits, nmaps_fits
	186	CHARACTER(LEN=80) :: units1
	187	CHARACTER(LEN=80),dimension(1:30) :: unitsm
	188	character(len=80), dimension(:), allocatable :: hbuffer
	189	! CHARACTER(LEN=80),dimension(:), allocatable :: unitsm
	190	integer(I4B) :: extno_i, extno_f, n_ext
	191	CHARACTER(LEN=80) :: polcconv
	192	integer(I4B) :: ipolcconv, polar, nside_fits
	193	character(len=*), parameter :: primer_url = 'http://healpix.sf.net/pdf/intro.pdf'
	194	character(len=*), parameter :: code = 'input_map'
	195
	196	!-----------------------------------------------------------------------
	197
	198	units1 = ' '
	199	unitsm(:) = ' '
	200	fmiss_effct = 0.
	201	imissing = 0
	202	if (PRESENT(fmissval)) fmiss_effct = fmissval
	203	if (PRESENT(header)) then
	204	nlheader = size(header)
	205	else
	206	nlheader = 36*100
	207	endif
	208	extno_i = 0
	209	if (present(extno)) extno_i = extno
	210	allocate(hbuffer(1:nlheader))
	211	do_polcconv = .true.
	212	if (present(ignore_polcconv)) then
	213	do_polcconv = .not. ignore_polcconv
	214	endif
	215
	216	n_ext = getnumext_fits(filename)
	217	obs_npix = getsize_fits(filename, nmaps = nmaps_fits, type=type_fits, extno=extno_i, &
	218	polcconv=ipolcconv, polarisation=polar, nside=nside_fits)
	219
	220	! if (nmaps_fits > nmaps) then
	221	! print*,trim(filename)//' only contains ',nmaps_fits,' maps'
	222	! print*,' You try to read ',nmaps
	223	! endif
	224	if (type_fits == 0 .or. type_fits == 2) then ! full sky map (in image or binary table)
	225	call read_bintab(filename, map(0:,1:), &
	226	& npixtot, nmaps, fmissing, anynull, header=hbuffer(1:), &
	227	& units=unitsm(1:), extno=extno_i)
	228	if (present(header)) then
	229	do i=1,nlheader
	230	header(i) = hbuffer(i)
	231	enddo
	232	endif
	233	if (present(units)) then
	234	units(1:size(units)) = unitsm(1:size(units))
	235	endif
	236	call map_bad_pixels(map, fmissing, fmiss_effct, imissing, verbose=.false.)
	237	if (imissing > 0) write(*,'(a,1pe11.4)') 'blank value : ' ,fmissing
	238	else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3) &
	239	& .and. (nmaps_fits > nmaps) .and. n_ext == 1) then
	240	! partial FITS file, reading 1 map (I) or 3 maps (I,Q,U) from the same extension
	241	obs_npix = getsize_fits(filename)
	242	allocate(pixel8(0:obs_npix-1), stat = status)
	243	call assert_alloc(status, code, "pixel8")
	244	allocate(iqu (0:obs_npix-1, 1:nmaps), stat = status)
	245	call assert_alloc(status, code, "iqu")
	246	call read_fits_partial(filename, pixel8, iqu, header=hbuffer)
	247	if (present(header)) then
	248	do i=1,nlheader
	249	header(i) = hbuffer(i)
	250	enddo
	251	endif
	252	if (present(units)) then
	253	do imap=1,min(size(units),nmaps-1)
	254	call get_card(hbuffer,'TUNIT'//string(imap+1,format='(i0.0)'), units(imap))
	255	enddo
	256	endif
	257	maxpix = maxval(pixel8)
	258	minpix = maxval(pixel8)
	259	if (maxpix > (npixtot-1) .or. minpix < 0) then
	260	print*,'map constructed from file '//trim(filename)//', with pixels in ',minpix,maxpix
	261	print*,' wont fit in array with ',npixtot,' elements'
	262	call fatal_error
	263	endif
	264	map(:,:) = fmiss_effct
	265	do imap=1,nmaps
	266	map(pixel8(:),imap) = iqu(:,imap)
	267	enddo
	268	imissing = npixtot - obs_npix
	269	deallocate(iqu)
	270	deallocate(pixel8)
	271	! else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3) .and. n_ext == nmaps) then
	272	! ! cut sky FITS file, reading 1 map (I) or 3 maps (I,Q,U), each from a different extension
	273	else if (type_fits == 3 .and. (nmaps == 3) .and. n_ext == nmaps) then
	274	! cut sky FITS file, reading 3 maps (I,Q,U), each from a different extension
	275	do imap = 1, nmaps
	276	extno_f = extno_i + imap - 1 ! 2016-08-16
	277	obs_npix = getsize_fits(filename, extno = extno_f)
	278	! one partial map (in binary table with explicit pixel indexing)
	279	allocate(pixel4(0:obs_npix-1), stat = status)
	280	allocate(signal(0:obs_npix-1), stat = status)
	281	call read_fits_cut4(filename, int(obs_npix,kind=i4b), &
	282	& pixel4, signal, header=hbuffer, units=units1, extno=extno_f)
	283	if (present(header) .and. imap == 1) then
	284	do i=1,nlheader
	285	header(i) = hbuffer(i)
	286	enddo
	287	endif
	288	if (present(units)) units(imap) = trim(units1)
	289	maxpix = maxval(pixel4)
	290	minpix = maxval(pixel4)
	291	if (maxpix > (npixtot-1) .or. minpix < 0) then
	292	print*,'map constructed from file '//trim(filename)//', with pixels in ',minpix,maxpix
	293	print*,' wont fit in array with ',npixtot,' elements'
	294	call fatal_error
	295	endif
	296	map(:,imap) = fmiss_effct
	297	map(pixel4(:),imap) = signal(:)
	298	imissing = npixtot - obs_npix
	299	deallocate(signal)
	300	deallocate(pixel4)
	301	enddo
	302	else
	303	print*,'Unable to read the ',nmaps,' required map(s) from file '//trim(filename)
	304	print*,'file type = ',type_fits
	305	call fatal_error
	306	endif
	307	!-----------------------------------------------------------------------
	308	if (imissing > 0) then
	309	write(*,'(i12,a,f7.3,a,1pe11.4)') &
	310	& imissing,' missing pixels (', &
	311	& (100.*imissing)/npixtot,' %),'// &
	312	& ' have been set to : ',fmiss_effct
	313	endif
	314
	315	! deallocate(unitsm)
	316
	317	!-----------------------------------------------------------------------
	318	! deal with polcconv
	319	! print,'****** in input_map ',nmaps, polar, type_fits, nside_fits
	320
	321	if ( do_polcconv .and. ( &
	322	& (nmaps >= 3 .and. type_fits == 3) & ! cut-sky, multiple maps
	323	& .or. &
	324	& (nmaps >=3 .and. polar==1 .and. type_fits ==2 .and. &
	325	& obs_npix==nside2npix(nside_fits)) & ! full-sky, polarized map
	326	& ) ) then
	327	if (ipolcconv == 0) then
	328	print 9000,' The POLCCONV keyword was not found in '//trim(filename)
	329	print 9000,' COSMO (HEALPix default) will be assumed, and map is unchanged.'
	330	print 9000,' See HEALPix primer ('//primer_url//') for details.'
	331	endif
	332	! if (ipolcconv == 1) then
	333	! print 9000,' POLCCONV=COSMO found in '//trim(filename)//'. Map is unchanged.'
	334	! endif
	335	if (ipolcconv == 2) then
	336	print 9000,' POLCCONV=IAU found in '//trim(filename)
	337	map(:,3) = - map(:,3)
	338	if (present(header)) then
	339	print 9000,' The sign of the U polarisation is changed in memory,'&
	340	& //' and the header is updated.'
	341	call add_card(header, 'POLCCONV', 'COSMO', &
	342	comment='Changed by input_map', update=.true.)
	343	else
	344	print 9000,' The sign of the U polarisation is changed in memory.'
	345	endif
	346	print 9000,' See HEALPix primer ('//primer_url//') for details.'
	347	endif
	348	if (ipolcconv == 3) then
	349	call get_card(hbuffer,'POLCCONV',polcconv)
	350	print 9000,' POLCCONV='//trim(polcconv)//' found in '//trim(filename)
	351	print 9000,' It is neither COSMO nor IAU. Aborting!'
	352	print 9000,' See HEALPix primer ('//primer_url//') for details.'
	353	call fatal_error
	354	endif
	355	endif
	356	9000 format(a)
	357	if (allocated(hbuffer)) deallocate(hbuffer)
	358
	359	RETURN
	360	END subroutine input_map8_KLOAD
	361	!=======================================================================
	362	! Read a FITS file
	363	! This routine is used for reading MAPS by anafast.
	364	! modified Feb 03 for units argument to run with Sun compiler
	365	! Jan 2005, EH, improved for faster writing
	366	!=======================================================================
	367	#ifndef NO64BITS
	368	subroutine read_bintab4_KLOAD(filename, map, npixtot, nmaps, nullval, anynull, header, units, extno)
	369	character(len=*), intent(IN) :: filename
	370	integer(I4B), intent(IN) :: npixtot
	371	integer(I4B), intent(IN) :: nmaps
	372	real(KMAP), dimension(0:,1:), intent(OUT) :: map
	373	real(KMAP), intent(OUT) :: nullval
	374	logical(LGT), intent(OUT) :: anynull
	375	character(LEN=*), dimension(1:), optional, intent(OUT) :: header
	376	character(LEN=*), dimension(1:), optional, intent(OUT) :: units
	377	integer(I4B) , optional, intent(IN) :: extno
	378
	379	integer(i8b):: npixtot8
	380
	381	npixtot8 = int(npixtot,kind=i8b)
	382	call read_bintab8_KLOAD(filename, map, npixtot8, nmaps, nullval, anynull, header, units, extno)
	383	return
	384
	385	end subroutine read_bintab4_KLOAD
	386	#endif
	387	subroutine read_bintab8_KLOAD(filename, map, npixtot, nmaps, nullval, anynull, header, units, extno)
	388	!=======================================================================
	389	character(len=*), intent(IN) :: filename
	390	integer(I8B), intent(IN) :: npixtot
	391	integer(I4B), intent(IN) :: nmaps
	392	real(KMAP), dimension(0:,1:), intent(OUT) :: map
	393	real(KMAP), intent(OUT) :: nullval
	394	logical(LGT), intent(OUT) :: anynull
	395	character(LEN=*), dimension(1:), optional, intent(OUT) :: header
	396	character(LEN=*), dimension(1:), optional, intent(OUT) :: units
	397	integer(I4B) , optional, intent(IN) :: extno
	398
	399	integer(I4B) :: nl_header, len_header, nl_units, len_units
	400	integer(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis
	401	integer(I4B) :: group, firstpix, i, npix32
	402	real(KMAP) :: blank, testval
	403	real(DP) :: bscale,bzero
	404	character(len=80) :: comment
	405	logical(LGT) :: extend
	406	integer(I4B) :: nmove, hdutype, hdunum
	407	integer(I4B) :: frow, imap
	408	integer(I4B) :: datacode, width
	409	LOGICAL(LGT) :: anynull_i
	410
	411	integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	412	integer(i8b) :: npix_old
	413	integer(i8b), dimension(1:MAXDIM) :: npix
	414	integer(i8b), dimension(1:MAXDIM) :: i0, i1
	415	integer(i4b), dimension(1:MAXDIM) :: repeat
	416	integer(i4b) :: nrow2read, nelem
	417
	418	integer(I4B) :: nrows, tfields, varidat
	419	character(len=20), dimension(1:MAXDIM) :: ttype, tform, tunit
	420	character(len=20) :: extname
	421	character(len=*), parameter :: code='read_bintab'
	422	!-----------------------------------------------------------------------
	423	status=0
	424
	425	unit = 146
	426	naxes(1) = 1
	427	naxes(2) = 1
	428	nfound = -1
	429	anynull = .false.
	430	bscale = 1.0d0
	431	bzero = 0.0d0
	432	blank = -2.e25
	433	nullval = bscale*blank + bzero
	434	comment=''
	435	ttype=''
	436	tform=''
	437	tunit=''
	438	extname=''
	439
	440	nl_header = 0
	441	if (present(header)) then
	442	nl_header = size(header)
	443	len_header = 80
	444	endif
	445
	446	nl_units = 0
	447	if (present(units)) then
	448	nl_units = size(units)
	449	len_units = min(80,len(units(1))) ! due to SUN compiler bug
	450	endif
	451
	452	readwrite=0
	453	!call ftopen(unit,filename,readwrite,blocksize,status)
	454	call ftnopn(unit,filename,readwrite, status) !open primary HDU or specified HDU
	455	if (status > 0) call printerror(status)
	456	! -----------------------------------------
	457	call ftghdn(unit, hdunum)
	458
	459	if (hdunum == 1) then ! in primary HDU
	460	! determines the presence of image
	461	call ftgkyj(unit,'NAXIS', naxis, comment, status)
	462	if (status > 0) call printerror(status)
	463
	464	! determines the presence of an extension
	465	call ftgkyl(unit,'EXTEND', extend, comment, status)
	466	if (status > 0) then
	467	extend = .false.
	468	status = 0 ! no extension :
	469	! to be compatible with first version of the code
	470	endif
	471	endif
	472
	473	if (naxis > 0 .and. .not.extend .and. hdunum==1) then ! there is an image
	474	! determine the size of the image (look naxis1 and naxis2)
	475	call ftgknj(unit,'NAXIS',1_i4b,2_i4b,naxes,nfound,status)
	476
	477	! check that it found only NAXIS1
	478	if (nfound == 2 .and. naxes(2) > 1) then
	479	print *,'multi-dimensional image'
	480	print *,'expected 1-D data.'
	481	call fatal_error
	482	end if
	483
	484	if (nfound < 1) then
	485	call printerror(status)
	486	print *,'can not find NAXIS1.'
	487	call fatal_error
	488	endif
	489
	490	npix(1)=naxes(1)
	491	if (npix(1) /= npixtot) then
	492	print *,'WARNING: found ',npix(1),' pixels in '//trim(filename)
	493	print *,' expected ',npixtot
	494	npix(1) = min(npix(1), npixtot)
	495	print *,' only ',npix(1),' will be read'
	496	endif
	497
	498	call ftgkyd(unit,'BSCALE',bscale,comment,status)
	499	if (status == 202) then ! BSCALE not found
	500	bscale = 1.0d0
	501	status = 0
	502	endif
	503	call ftgkyd(unit,'BZERO', bzero, comment,status)
	504	if (status == 202) then ! BZERO not found
	505	bzero = 0.0d0
	506	status = 0
	507	endif
	508	call f90ftgky_(unit, 'BLANK', blank, comment, status)
	509	if (status == 202) then ! BLANK not found
	510	! (according to fitsio BLANK is integer)
	511	blank = -2.e25
	512	status = 0
	513	endif
	514	nullval = bscale*blank + bzero
	515
	516	! -----------------------------------------
	517
	518	group=1
	519	firstpix = 1
	520	npix32 = npix(1)
	521	call f90ftgpv_(unit, group, firstpix, npix32, nullval, map(0:npix(1)-1,1), anynull, status)
	522	! if there are any NaN pixels, (real data)
	523	! or BLANK pixels (integer data) they will take nullval value
	524	! and anynull will switch to .true.
	525	! otherwise, switch it by hand if necessary
	526	testval = 1.e-6 * ABS(nullval)
	527	do i=0, npix(1)-1
	528	if (ABS(map(i,1)-nullval) < testval) then
	529	anynull = .true.
	530	goto 111
	531	endif
	532	enddo
	533	111 continue
	534
	535	else if (extend .or. hdunum>1) then
	536	if (hdunum == 1) then
	537	nmove = +1
	538	if (present(extno)) nmove = +1 + extno
	539	call ftmrhd(unit, nmove, hdutype, status)
	540	else
	541	call ftghdt(unit, hdutype, status)
	542	endif
	543
	544	call assert (hdutype==2, 'this is not a binary table')
	545
	546	! reads all the keywords
	547	call ftghbn(unit, MAXDIM, &
	548	& nrows, tfields, ttype, tform, tunit, extname, varidat, &
	549	& status)
	550
	551	if (tfields < nmaps) then
	552	print *,'found ',tfields,' maps in file '//trim(filename)
	553	print *,'expected ',nmaps
	554	call fatal_error
	555	endif
	556
	557	! finds the bad data value
	558	call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
	559	if (status == 202) then ! bad_data not found
	560	if (KMAP == SP) nullval = s_bad_value ! default value
	561	if (KMAP == DP) nullval = d_bad_value ! default value
	562	status = 0
	563	endif
	564
	565	if (nl_header > 0) then
	566	do i=1,nl_header
	567	header(i)(1:len_header) = ""
	568	enddo
	569	call get_clean_header(unit, header, filename, status)
	570	status = 0
	571	endif
	572
	573	if (nl_units > 0) then
	574	do i=1,nl_units
	575	units(i)(1:len_units) = 'unknown' ! default
	576	enddo
	577	do imap = 1, min(nmaps, nl_units)
	578	units(imap)(1:len_units) = adjustl(tunit(imap))
	579	enddo
	580	endif
	581
	582	npix_old = npixtot
	583	do imap = 1, nmaps
	584	!parse TFORM keyword to find out the length of the column vector
	585	call ftbnfm(tform(imap), datacode, repeat(imap), width, status)
	586	npix(imap) = int(nrows,i8b) * repeat(imap)
	587	if (npix(imap) /= npixtot .and. npix_old /= npix(imap)) then
	588	print *,'WARNING: found ',npix(imap),' pixels in '//trim(filename)//', column ',imap
	589	print *,' expected ',npixtot,' or ',npix_old
	590	npix_old = npix(imap)
	591	npix(imap) = min(npix(imap), npixtot)
	592	print *,' only ',npix(imap),' will be read'
	593	endif
	594	enddo
	595
	596
	597	call ftgrsz(unit, nrow2read, status)
	598	nrow2read = max(nrow2read, 1)
	599	firstpix = 1 ! starting position in FITS within row, 1 based
	600	i0(:) = 0_i8b ! starting element in array, 0 based
	601	do frow = 1, nrows, nrow2read
	602	do imap = 1, nmaps
	603	i1(imap) = min(i0(imap) + int(nrow2read,i8b) * repeat(imap), npix(imap)) - 1_i8b
	604	nelem = i1(imap) - i0(imap) + 1
	605	call f90ftgcv_(unit, imap, frow, firstpix, nelem, &
	606	& nullval, map(i0(imap):i1(imap),imap), anynull_i, status)
	607	anynull = anynull .or. anynull_i
	608	i0(imap) = i1(imap) + 1_i8b
	609	enddo
	610	enddo
	611	! sanity check
	612	do imap = 1, nmaps
	613	if (i0(imap) /= npix(imap)) then
	614	call fatal_error('something wrong during piece wise reading')
	615	endif
	616	enddo
	617
	618	else ! no image no extension, you are dead, man
	619	call fatal_error(' No image, no extension in '//trim(filename))
	620	endif
	621	! close the file
	622	call ftclos(unit, status)
	623
	624	! check for any error, and if so print out error messages
	625	if (status > 0) call printerror(status)
	626
	627	return
	628	end subroutine read_bintab8_KLOAD
	629
	630	!=======================================================================
	631	subroutine read_conbintab_KLOAD(filename, alms, nalms, units, extno)
	632	!=======================================================================
	633	! Read a FITS file containing alms values
	634	!
	635	! slightly modified to deal with vector column
	636	! in binary table EH/IAP/Jan-98
	637	!
	638	! Used by synfast when reading a binary file with alms for cons.real.
	639	! FKH/Apr-99
	640	!
	641	! extno : optional, number of extension to be read, default=0: first extension
	642	! Jan 2005, EH, improved for faster reading
	643	!=======================================================================
	644	CHARACTER(LEN=*), INTENT(IN) :: filename
	645	INTEGER(I4B), INTENT(IN) :: nalms !, nlheader
	646	REAL(KMAP), DIMENSION(0:nalms-1,1:6), INTENT(OUT) :: alms
	647	CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units
	648	INTEGER(I4B) , optional, INTENT(IN) :: extno
	649
	650	REAL(KMAP) :: nullval
	651	LOGICAL(LGT) :: anynull
	652
	653	INTEGER(I4B), DIMENSION(:), allocatable :: lm
	654	INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis
	655	INTEGER(I4B) :: npix
	656	CHARACTER(LEN=80) :: comment
	657	LOGICAL(LGT) :: extend
	658	INTEGER(I4B) :: nmove, hdutype
	659	INTEGER(I4B) :: frow, imap
	660	INTEGER(I4B) :: datacode, repeat, width
	661	integer(I4B) :: i, l, m
	662	integer(i4b) :: nrow2read, nelem
	663	integer(i8b) :: i0, i1
	664
	665	INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	666	INTEGER(I4B) :: nrows, tfields, varidat
	667	CHARACTER(LEN=20), dimension(1:MAXDIM) :: ttype, tform, tunit
	668	CHARACTER(LEN=20) :: extname
	669	character(len=*), parameter :: code="read_conbintab"
	670
	671	!-----------------------------------------------------------------------
	672	status=0
	673
	674	unit = 145
	675	naxes(1) = 1
	676	naxes(2) = 1
	677	nfound = -1
	678	anynull = .false.
	679	alms=0. ! set output to 0.
	680	readwrite=0
	681	comment=''
	682	ttype=''
	683	tform=''
	684	tunit=''
	685	extname=''
	686	call ftopen(unit,filename,readwrite,blocksize,status)
	687	if (status > 0) call printerror(status)
	688	! -----------------------------------------
	689
	690	! determines the presence of image
	691	call ftgkyj(unit,'NAXIS', naxis, comment, status)
	692	if (status > 0) call printerror(status)
	693
	694	! determines the presence of an extension
	695	call ftgkyl(unit,'EXTEND', extend, comment, status)
	696	if (status > 0) status = 0 ! no extension :
	697	! to be compatible with first version of the code
	698
	699	if (.not.extend) then
	700	print*,'No extension!'
	701	call fatal_error
	702	endif
	703
	704	! go to assigned extension (if it exists)
	705	nmove = +1
	706	if (present(extno)) nmove = +1 + extno
	707	call ftmrhd(unit, nmove, hdutype, status)
	708	if (status > 0) then
	709	! if required extension not found:
	710	! print a warning, fill with dummy values, return to calling routine
	711	print*,code//' WARNING: the extension ',extno,' was not found in ',trim(filename)
	712	alms(0:nalms-1,1)=0. ! l = 0
	713	alms(0:nalms-1,2)=1. ! m = 1
	714	alms(0:nalms-1,3:6)=0.
	715	status = 0
	716	call ftclos(unit, status) ! close the file
	717	return
	718	endif
	719
	720	if (hdutype /= 2) then ! not a binary table
	721	print*, 'this is not a binary table'
	722	call fatal_error
	723	endif
	724
	725	! reads all the keywords
	726	call ftghbn(unit, MAXDIM, &
	727	& nrows, tfields, ttype, tform, tunit, extname, varidat, &
	728	& status)
	729
	730	if ((tfields/=3).and.(tfields/=5)) then
	731	print *,'found ',tfields,' columns in the file'
	732	print *,'expected 3 or 5'
	733	call fatal_error
	734	endif
	735	! finds the bad data value
	736	! if (KMAP == SP) call ftgkye(unit,'BAD_DATA',nullval,comment,status)
	737	! if (KMAP == DP) call ftgkyd(unit,'BAD_DATA',nullval,comment,status)
	738	call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
	739	if (status == 202) then ! bad_data not found
	740	if (KMAP == SP) nullval = s_bad_value ! default value
	741	if (KMAP == DP) nullval = d_bad_value ! default value
	742	status = 0
	743	endif
	744
	745	! if (nlheader > 0) then
	746	! header = ""
	747	! status = 0
	748	! call get_clean_header(unit, header, filename, status)
	749	! endif
	750
	751	if (present(units)) then
	752	units(1) = tunit(2)
	753	endif
	754
	755	!parse TFORM keyword to find out the length of the column vector
	756	call ftbnfm(tform(1), datacode, repeat, width, status)
	757	npix = nrows * repeat
	758	if (npix /= nalms) then
	759	print *,'found ',npix,' alms'
	760	print *,'expected ',nalms
	761	call fatal_error
	762	endif
	763
	764	call ftgrsz(unit, nrow2read, status)
	765	nrow2read = max(nrow2read, 1)
	766	nelem = nrow2read * repeat
	767	i0 = 0_i8b
	768	allocate(lm(0:nelem-1))
	769	do frow = 1, nrows, nrow2read
	770	i1 = min(i0 + nrow2read * repeat, int(npix,i8b)) - 1_i8b
	771	nelem = i1 - i0 + 1
	772	! first column -> index
	773	call ftgcvj(unit, 1_i4b, frow, 1_i4b, nelem, i_bad_value, &
	774	& lm(0), anynull, status)
	775	call assert(.not. anynull, 'There are undefined values in the table!')
	776	! other columns -> a(l,m)
	777	do imap = 2, tfields
	778	call f90ftgcv_(unit, imap, frow, 1_i4b, nelem, nullval, &
	779	& alms(i0:i1,imap+1), anynull, status)
	780	call assert (.not. anynull, 'There are undefined values in the table!')
	781	enddo
	782	! recoding of the mapping, EH, 2002-08-08
	783	! lm = l*l + l + m + 1
	784	do i = i0, i1
	785	l = int(sqrt( real(lm(i-i0)-1, kind = DP) ) )
	786	m = lm(i-i0) - l*(l+1) - 1
	787	! check that round-off did not screw up mapping
	788	if (abs(m) > l) then
	789	print*,'Inconsistent l^2+l+m+1 -> l,m mapping'
	790	print,l, m, l(l+1)+m+1, lm(i-i0)
	791	call fatal_error
	792	endif
	793	alms(i,1) = real( l, kind=KMAP)
	794	alms(i,2) = real( m, kind=KMAP)
	795	enddo
	796	i0 = i1 + 1_i8b
	797	enddo
	798	deallocate(lm)
	799	! sanity check
	800	if (i0 /= npix) then
	801	print*,'something wrong during piece-wise reading'
	802	call fatal_error
	803	endif
	804
	805	! !reads the columns
	806	! ! first column : index -> (l,m)
	807	! allocate(lm(0:nalms-1))
	808	! column = 1
	809	! frow = 1
	810	! firstpix = 1
	811	! npix = nrows * repeat
	812	! if (npix /= nalms) then
	813	! print *,'found ',npix,' alms'
	814	! print *,'expected ',nalms
	815	! call fatal_error
	816	! endif
	817	! call ftgcvj(unit, column, frow, firstpix, npix, i_bad_value, &
	818	! & lm(0), anynull, status)
	819	! call assert (.not. anynull, 'There are undefined values in the table!')
	820	! ! recoding of the mapping, EH, 2002-08-08
	821	! ! lm = l*l + l + m + 1
	822	! do i = 0, nalms - 1
	823	! l = int(sqrt( real(lm(i)-1, kind = DP) ) )
	824	! m = lm(i) - l*(l+1) - 1
	825	! ! check that round-off did not screw up mapping
	826	! if (abs(m) > l) then
	827	! print*,'Inconsistent l^2+l+m+1 -> l,m mapping'
	828	! print,l, m, l(l+1)+m+1, lm(i)
	829	! call fatal_error
	830	! endif
	831	! alms(i,1) = real( l, kind=KMAP )
	832	! alms(i,2) = real( m, kind=KMAP )
	833	! enddo
	834	! deallocate(lm)
	835	!
	836	! do imap = 2, tfields
	837	! !parse TFORM keyword to find out the length of the column vector
	838	! call ftbnfm(tform(imap), datacode, repeat, width, status)
	839	!
	840	! !reads the columns
	841	! column = imap
	842	! frow = 1
	843	! firstpix = 1
	844	! npix = nrows * repeat
	845	! if (npix /= nalms) then
	846	! print *,'found ',npix,' alms'
	847	! print *,'expected ',nalms
	848	! call fatal_error
	849	! endif
	850	! call f90ftgcv_(unit, column, frow, firstpix, npix, nullval, &
	851	! & alms(0:npix-1,imap+1), anynull, status)
	852	! call assert (.not. anynull, 'There are undefined values in the table!')
	853	! enddo
	854	! close the file
	855	call ftclos(unit, status)
	856
	857
	858	! check for any error, and if so print out error messages
	859	if (status > 0) call printerror(status)
	860	return
	861	end subroutine read_conbintab_KLOAD
	862
	863	#ifndef NO64BITS
	864	subroutine write_bintab4_KLOAD(map, npix, nmap, header, nlheader, filename, extno)
	865	!=======================================================================
	866	INTEGER(I4B), INTENT(IN) :: npix, nmap, nlheader
	867	REAL(KMAP), INTENT(IN), DIMENSION(0:npix-1,1:nmap) :: map
	868	CHARACTER(LEN=*), INTENT(IN), DIMENSION(1:nlheader) :: header
	869	CHARACTER(LEN=*), INTENT(IN) :: filename
	870	INTEGER(I4B) , INTENT(IN), optional :: extno
	871
	872	integer(i8b) :: npix8
	873
	874	npix8 = npix
	875	if (present(extno)) then
	876	call write_bintab8_KLOAD(map, npix8, nmap, header, nlheader, filename, extno)
	877	else
	878	call write_bintab8_KLOAD(map, npix8, nmap, header, nlheader, filename)
	879	endif
	880	return
	881	end subroutine write_bintab4_KLOAD
	882	#endif
	883
	884	subroutine write_bintab8_KLOAD(map, npix, nmap, header, nlheader, filename, extno)
	885	!=======================================================================
	886	! Create a FITS file containing a binary table extension with
	887	! the temperature map in the first column
	888	! written by EH from writeimage and writebintable
	889	! (fitsio cookbook package)
	890	!
	891	! slightly modified to deal with vector column (ie TFORMi = '1024E')
	892	! in binary table EH/IAP/Jan-98
	893	!
	894	! simplified the calling sequence, the header sould be filled in
	895	! before calling the routine
	896	!
	897	! July 21, 2004: SP version
	898	! Jan 2005, EH, improved for faster writing
	899	! 2009-02-25: EH, accepts I4B and I8B npix
	900	!=======================================================================
	901	use pix_tools, only: npix2nside
	902	INTEGER(I8B), INTENT(IN) :: npix
	903	INTEGER(I4B), INTENT(IN) :: nmap, nlheader
	904	! REAL(KMAP), INTENT(IN), DIMENSION(0:npix-1,1:nmap), target :: map
	905	REAL(KMAP), INTENT(IN), DIMENSION(0:npix-1,1:nmap) :: map
	906	CHARACTER(LEN=*), INTENT(IN), DIMENSION(1:nlheader) :: header
	907	CHARACTER(LEN=*), INTENT(IN) :: filename
	908	INTEGER(I4B) , INTENT(IN), optional :: extno
	909
	910	INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1)
	911	INTEGER(I4B) :: i, nside
	912	LOGICAL(LGT) :: simple,extend
	913	CHARACTER(LEN=80) :: comment, ch
	914
	915	INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	916	INTEGER(I4B) :: nrows, tfields, varidat
	917	integer(i4b) :: repeat, nrow2write, nelem
	918	integer(i8b) :: i0, i1
	919	INTEGER(I4B) :: frow, felem, colnum, hdutype
	920	CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
	921	CHARACTER(LEN=10) :: card
	922	CHARACTER(LEN=2) :: stn
	923	INTEGER(I4B) :: itn, extno_i
	924	character(len=filenamelen) sfilename
	925	character(len=1) :: pform
	926	! character(len=80) :: junk_k, junk_v, junk_c
	927	!-----------------------------------------------------------------------
	928
	929	if (KMAP == SP) pform = 'E'
	930	if (KMAP == DP) pform = 'D'
	931	status=0
	932	unit = 144
	933	blocksize=1
	934
	935	extno_i = 0
	936	if (present(extno)) extno_i = extno
	937
	938	if (extno_i == 0) then
	939	!*************************************
	940	! create the new empty FITS file
	941	!*************************************
	942	call ftinit(unit,filename,blocksize,status)
	943	if (status > 0) call fatal_error("Error while creating file " &
	944	& //trim(filename) &
	945	& //". Check path and/or access rights.")
	946
	947	! -----------------------------------------------------
	948	! initialize parameters about the FITS image
	949	simple=.true.
	950	bitpix=32 ! integer*4
	951	naxis=0 ! no image
	952	naxes(1)=0
	953	extend=.true. ! there is an extension
	954
	955	! ----------------------
	956	! primary header
	957	! ----------------------
	958	! write the required header keywords
	959	call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
	960
	961	! writes supplementary keywords : none
	962
	963	! write the current date
	964	call ftpdat(unit,status) ! format (yyyy-mm-dd)
	965
	966	! ----------------------
	967	! image : none
	968	! ----------------------
	969
	970	! ----------------------
	971	! extension
	972	! ----------------------
	973	else
	974
	975	!*********************************************
	976	! reopen an existing file and go to the end
	977	!*********************************************
	978	! remove the leading '!' (if any) when reopening the same file
	979	sfilename = adjustl(filename)
	980	if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
	981	call ftopen(unit,sfilename,1_i4b,blocksize,status)
	982	call ftmahd(unit,1_i4b+extno_i,hdutype,status)
	983
	984	endif
	985
	986	! creates an extension
	987	call ftcrhd(unit, status)
	988
	989	! writes required keywords
	990	tfields = nmap
	991	repeat = 1024
	992	if (npix < 1024) repeat = 1 ! for nside <= 8
	993	! for large npix increase repeat so that nrows < 2^31-1
	994	nside = npix2nside(npix)
	995	if (nside > 1024*256) repeat = nside/256
	996	nrows = npix / repeat ! naxis1
	997	ch = string(repeat, format='(i8)')
	998	ch = trim(adjustl(ch))//pform
	999	tform(1:nmap) = ch
	1000	ttype(1:nmap) = 'simulation' ! will be updated
	1001	tunit(1:nmap) = '' ! optional, will not appear
	1002	extname = '' ! optional, will not appear
	1003	varidat = 0
	1004	call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
	1005	& extname, varidat, status)
	1006
	1007	! write the header literally, putting TFORM1 and TUNIT1 at the desired place
	1008	if (KMAP == SP) comment = 'data format of field: 4-byte REAL'
	1009	if (KMAP == DP) comment = 'data format of field: 8-byte REAL'
	1010	do i=1,nlheader
	1011	card = header(i)
	1012	if (card(1:5) == 'TTYPE') then ! if TTYPEi is explicitely given
	1013	stn = card(6:8)
	1014	read(stn,'(i3)') itn
	1015	! discard at their original location:
	1016	call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and
	1017	status = 0
	1018	call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
	1019	status = 0
	1020	if (itn <= tfields) then ! only put relevant information 2008-08-27
	1021	call putrec(unit,header(i), status) ! write new TTYPE1
	1022	status = 0
	1023	call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after
	1024	endif
	1025	elseif (header(i)/=' ') then
	1026	call putrec(unit,header(i), status)
	1027	endif
	1028	status = 0
	1029	enddo
	1030
	1031	! write the extension buffer by buffer
	1032	call ftgrsz(unit, nrow2write, status)
	1033	nrow2write = max(nrow2write, 1)
	1034	felem = 1 ! starting position in FITS (element), 1 based
	1035	i0 = 0_i8b ! starting element in array, 0 based
	1036	do frow = 1, nrows, nrow2write
	1037	i1 = min(i0 + nrow2write * repeat, npix) - 1_i8b
	1038	nelem = i1 - i0 + 1
	1039	do colnum = 1, nmap
	1040	call f90ftpcl_(unit, colnum, frow, felem, nelem, map(i0:i1, colnum), status)
	1041	enddo
	1042	i0 = i1 + 1_i8b
	1043	enddo
	1044	! sanity check
	1045	if (i0 /= npix) then
	1046	call fatal_error("something wrong during piece wise writing")
	1047	endif
	1048
	1049	! ----------------------
	1050	! close and exit
	1051	! ----------------------
	1052
	1053	! close the file and free the unit number
	1054	call ftclos(unit, status)
	1055
	1056	! check for any error, and if so print out error messages
	1057	if (status > 0) call printerror(status)
	1058
	1059	return
	1060	end subroutine write_bintab8_KLOAD
	1061	!=======================================================================
	1062	subroutine write_asctab_KLOAD &
	1063	& (clout, lmax, ncl, header, nlheader, filename, extno)
	1064	!=======================================================================
	1065	! Create a FITS file containing an ASCII table extension with
	1066	! the measured power spectra
	1067	! written by EH from writeimage and writeascii
	1068	! (fitsio cookbook package)
	1069	!
	1070	! clout = power spectra with l in (0:lmax)
	1071	! ncl = number of spectra
	1072	! header = FITS header to be put on top of the file
	1073	! nlheader = number of lines of the header
	1074	! filename = FITS output file name
	1075	!
	1076	! 2015-09-10: EH, added extno
	1077	!=======================================================================
	1078	!
	1079	INTEGER(I4B), INTENT(IN) :: lmax, ncl,nlheader
	1080	REAL(KMAP), INTENT(IN) :: clout(0:lmax,1:ncl)
	1081	CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header
	1082	CHARACTER(LEN=*), INTENT(IN) :: filename
	1083	INTEGER(I4B) , INTENT(IN), optional :: extno
	1084
	1085	INTEGER(I4B) :: bitpix,naxis,naxes(1)
	1086	LOGICAL(LGT) :: simple,extend
	1087	CHARACTER(LEN=10) :: card
	1088
	1089	INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP
	1090	INTEGER(I4B) :: status,unit,blocksize,tfields,nrows,rowlen
	1091	INTEGER(I4B) :: nspace,tbcol(MAXDIM),colnum,frow,felem,hdutype
	1092	CHARACTER(LEN=16) :: ttype(MAXDIM),tform(MAXDIM),tunit(MAXDIM),extname
	1093	CHARACTER(LEN=80) :: comment, card_tbcol
	1094	CHARACTER(LEN=2) :: stn
	1095	INTEGER(I4B) :: itn, i, extno_i
	1096	character(len=filenamelen) sfilename
	1097	character(len=6) :: form
	1098	!=======================================================================
	1099	if (KMAP == SP) form = 'E15.7'
	1100	if (KMAP == DP) form = 'D24.15'
	1101	status=0
	1102	unit = 109
	1103	blocksize=1
	1104
	1105	extno_i = 0
	1106	if (present(extno)) extno_i = extno
	1107
	1108	if (extno_i == 0) then
	1109	call ftinit(unit,filename,blocksize,status)
	1110	if (status > 0) call fatal_error("Error while creating file " &
	1111	& //trim(filename) &
	1112	& //". Check path and/or access rights.")
	1113
	1114	! -----------------------------------------------------
	1115	! initialize parameters about the FITS image
	1116	simple=.true.
	1117	bitpix=32 ! integer*4
	1118	naxis=0 ! no image
	1119	naxes(1)=0
	1120	extend=.true. ! there is an extension
	1121
	1122	! ----------------------
	1123	! primary header
	1124	! ----------------------
	1125	! write the required header keywords
	1126	call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
	1127
	1128	! writes supplementary keywords : none
	1129
	1130	! write the current date
	1131	call ftpdat(unit,status) ! (format ccyy-mm-dd)
	1132
	1133	! ----------------------
	1134	! image : none
	1135	! ----------------------
	1136
	1137	! ----------------------
	1138	! extension
	1139	! ----------------------
	1140	else
	1141	!*********************************************
	1142	! reopen an existing file and go to the end
	1143	!*********************************************
	1144	! remove the leading '!' (if any) when reopening the same file
	1145	sfilename = adjustl(filename)
	1146	if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
	1147	call ftopen(unit,sfilename,1_i4b,blocksize,status)
	1148	call ftmahd(unit,1_i4b+extno_i,hdutype,status)
	1149
	1150	endif
	1151
	1152	! append a new empty extension onto the end of the primary array
	1153	call ftcrhd(unit,status)
	1154
	1155	! define parameters for the ASCII table
	1156	nrows = lmax+1
	1157	tfields = ncl
	1158	tform(1:ncl) = form
	1159	ttype(1:ncl) = 'power spectrum' ! is updated by the value given in the header
	1160	tunit(1:ncl) = '' ! optional, will not appear
	1161	extname = '' ! optional, will not appear
	1162
	1163	! calculate the starting position of each column, and the total row length
	1164	nspace=1
	1165	call ftgabc(tfields,tform,nspace,rowlen,tbcol,status)
	1166
	1167	! write the required header parameters for the ASCII table
	1168	call ftphtb(unit,rowlen,nrows,tfields,ttype,tbcol,tform,tunit, &
	1169	& extname,status)
	1170
	1171	! write the header literally, putting TFORM1 at the desired place
	1172	comment = ''
	1173	card_tbcol=''
	1174	do i=1,nlheader
	1175	card = header(i)
	1176	if (card(1:5) == 'TTYPE') then ! if TTYPE1 is explicitely given
	1177	stn = card(6:8)
	1178	read(stn,'(i3)') itn
	1179	! discard at their original location:
	1180	!!!!!!!call ftdkey(unit,card(1:6),status) ! old TTYPEi 2019-10-14
	1181	call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi
	1182	status = 0
	1183	call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
	1184	status = 0
	1185	call ftgcrd(unit,'TBCOL'//stn,card_tbcol,status)
	1186	status = 0
	1187	call ftdkey(unit,'TBCOL'//stn,status) ! TBCOLi
	1188	status = 0
	1189	! and rewrite
	1190	if (itn <= tfields) then
	1191	call putrec(unit,card_tbcol,status) ! TBCOLi
	1192	status = 0
	1193	call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! TFORMi right after
	1194	status = 0
	1195	call putrec(unit,header(i), status) ! TTYPEi
	1196	endif
	1197	elseif (header(i)/=' ') then
	1198	call putrec(unit,header(i), status)
	1199	endif
	1200	status = 0
	1201	enddo
	1202
	1203	frow=1
	1204	felem=1
	1205	do colnum = 1, ncl
	1206	call f90ftpcl_(unit, colnum, frow, felem, nrows, clout(0:nrows-1,colnum), status)
	1207	enddo
	1208
	1209	! close the FITS file
	1210	call ftclos(unit, status)
	1211
	1212	! check for any error, and if so print out error messages
	1213	if (status > 0) call printerror(status)
	1214
	1215	return
	1216	end subroutine write_asctab_KLOAD
	1217	!=======================================================================
	1218	! DUMP_ALMS
	1219	!
	1220	! Create/extend a FITS file containing a binary table extension with
	1221	! the a_lm coefficients in each extension.
	1222	!
	1223	! Format of alm FITS file: in each extension, 3 columns:
	1224	! index=l*l+l+m+1, real(a_lm), imag(a_lm)
	1225	! the a_lm are obtained using the complex Y_lm
	1226	! only the modes with m>=0 are stored (because the map is assumed real)
	1227	!
	1228	! The extensions contain, in general, T, E (= G) and B (= C) in that order
	1229	! First extension has extno = 0
	1230	!
	1231	! Adapted from write_bintab, FKH/Apr-99
	1232	! SP/DP overloaded, 2004-12, EH
	1233	! reduced size of intermediate arrays (alms_out, lm)
	1234	!
	1235	subroutine dump_alms_KLOAD(filename, alms, nlmax, header, nlheader, extno)
	1236	!=======================================================================
	1237	INTEGER(I4B), INTENT(IN) :: nlmax, nlheader, extno
	1238	COMPLEX(KMAPC), INTENT(IN), DIMENSION(0:,0:) :: alms
	1239	CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header
	1240	CHARACTER(LEN=*), INTENT(IN) :: filename
	1241	! local variables
	1242	INTEGER(I4B), DIMENSION(:), allocatable :: lm
	1243	REAL(KMAP), DIMENSION(:,:), allocatable :: alms_out
	1244	INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1)
	1245	INTEGER(I4B) :: i,l,m,cnt,hdutype, nmmax
	1246	LOGICAL(LGT) :: simple,extend, found_lmax, found_mmax
	1247	CHARACTER(LEN=80) :: comment
	1248
	1249	integer(i8b) :: npix
	1250	INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	1251	INTEGER(I4B) :: nrows, tfields, varidat
	1252	INTEGER(I4B) :: frow, felem
	1253	CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
	1254	CHARACTER(LEN=10) :: card
	1255	CHARACTER(LEN=2) :: stn
	1256	INTEGER(I4B) :: itn
	1257	character(len=filenamelen) sfilename
	1258	character(len=1) :: pform
	1259	character(len=*), parameter :: code = 'dump_alms'
	1260	!-----------------------------------------------------------------------
	1261
	1262	if (KMAP == SP) pform = 'E'
	1263	if (KMAP == DP) pform = 'D'
	1264
	1265	nmmax = size(alms,2) - 1
	1266	if (nmmax < 0 .or. nmmax > nlmax) call fatal_error('inconsistent Lmax and Mmax in dump_alms')
	1267	!! npix=((nlmax+1)*(nlmax+2))/2
	1268	npix = ((nmmax+1_i8b)(2_i8bnlmax-nmmax+2))/2
	1269	if (npix > MAX_I4B .or. nlmax >= int(sqrt(real(MAX_I4B,kind=dp)))) then
	1270	print*,code//'> Index of a_lm too large for current file format,'
	1271	print*,code//'> '//trim(filename)//' not written.'
	1272	call fatal_error
	1273	endif
	1274
	1275	status=0
	1276	unit = 141
	1277	blocksize=1
	1278
	1279	if (extno==0) then
	1280	!*********************************************
	1281	! create the new empty FITS file
	1282	!*********************************************
	1283	call ftinit(unit,filename,blocksize,status)
	1284	if (status > 0) call fatal_error("Error while creating file " &
	1285	& //trim(filename) &
	1286	& //". Check path and/or access rights.")
	1287
	1288	! -----------------------------------------------------
	1289	! initialize parameters about the FITS image
	1290	simple=.true.
	1291	bitpix=32 ! integer*4
	1292	naxis=0 ! no image
	1293	naxes(1)=0
	1294	extend=.true. ! there is an extension
	1295
	1296	! ----------------------
	1297	! primary header
	1298	! ----------------------
	1299	! write the required header keywords
	1300	call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
	1301
	1302	! writes supplementary keywords : none
	1303
	1304	! write the current date
	1305	call ftpdat(unit,status) ! format (yyyy-mm-dd)
	1306
	1307	! ----------------------
	1308	! image : none
	1309	! extension
	1310	! ----------------------
	1311	else
	1312
	1313	!*********************************************
	1314	! reopen an existing file and go to the end
	1315	!*********************************************
	1316	! remove the leading '!' (if any) when reopening the same file
	1317	sfilename = adjustl(filename)
	1318	if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
	1319	call ftopen(unit,sfilename,1_i4b,blocksize,status)
	1320	call ftmahd(unit,1_i4b+extno,hdutype,status)
	1321
	1322	endif
	1323
	1324	! creates an extension
	1325	call ftcrhd(unit, status)
	1326
	1327	! writes required keywords
	1328	nrows = npix ! naxis1
	1329	tfields = 3
	1330	tform = '' ; ttype = '' ; tunit = ''
	1331	tform(1)(1:2)= '1J' ! necessary with ifort 11.1.076
	1332	tform(2)(1:2) = '1'//pform
	1333	tform(3)(1:2) = '1'//pform
	1334	ttype(1)(1:15) = 'index=l^2+l+m+1'
	1335	ttype(2)(1:15) = 'alm (real) '
	1336	ttype(3)(1:15) = 'alm (imaginary)'
	1337	tunit(1:3) = '' ! optional, will not appear
	1338	extname = '' ! optional, will not appear
	1339	varidat = 0
	1340	call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
	1341	& extname, varidat, status)
	1342
	1343	! write the header literally, putting TFORM1 at the desired place
	1344	do i=1,nlheader
	1345	card = header(i)
	1346	if (card(1:5) == 'TTYPE') then ! if TTYPEi is explicitely given
	1347	stn = card(6:8)
	1348	read(stn,'(i3)') itn
	1349	! discard at their original location:
	1350	call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and ! remove
	1351	status = 0
	1352	call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
	1353	status = 0
	1354	if (itn <= tfields) then ! only put relevant information 2008-08-27
	1355	call putrec(unit,header(i), status) ! write new TTYPEi
	1356	status = 0
	1357	if (itn==1) then
	1358	comment = 'data format of field: 4-byte INTEGER'
	1359	else
	1360	if (KMAP == SP) comment = 'data format of field: 4-byte REAL'
	1361	if (KMAP == DP) comment = 'data format of field: 8-byte REAL'
	1362	endif
	1363	call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after
	1364	endif
	1365	elseif (header(i)/=' ') then
	1366	call putrec(unit,header(i), status)
	1367	endif
	1368	status = 0
	1369	enddo
	1370
	1371	call ftukyj(unit, 'MIN-LPOL', 0_i4b, 'Minimum L multipole order', status)
	1372	call ftukyj(unit, 'MAX-LPOL', nlmax, 'Maximum L multipole order', status)
	1373	call ftukyj(unit, 'MAX-MPOL', nmmax, 'Maximum M multipole degree', status)
	1374
	1375	allocate(lm (0:nmmax))
	1376	allocate(alms_out(0:nmmax,1:2))
	1377	! write the extension one column by one column
	1378	frow = 1 ! starting position (row)
	1379	felem = 1 ! starting position (element)
	1380	do l = 0, nlmax
	1381	cnt = 0
	1382	do m = 0, min(l,nmmax)
	1383	lm(cnt) = l**2 + l + m + 1
	1384	alms_out(cnt,1)=REAL( alms(l,m))
	1385	alms_out(cnt,2)=AIMAG(alms(l,m))
	1386	cnt = cnt + 1
	1387	enddo
	1388	call ftpclj( unit, 1_i4b, frow, felem, cnt, lm, status)
	1389	call f90ftpcl_(unit, 2_i4b, frow, felem, cnt, alms_out(0:cnt-1,1), status)
	1390	call f90ftpcl_(unit, 3_i4b, frow, felem, cnt, alms_out(0:cnt-1,2), status)
	1391	frow = frow + cnt
	1392	enddo
	1393	deallocate(lm)
	1394	deallocate(alms_out)
	1395
	1396	! ----------------------
	1397	! close and exit
	1398	! ----------------------
	1399
	1400	! close the file and free the unit number
	1401	call ftclos(unit, status)
	1402
	1403	! check for any error, and if so print out error messages
	1404	if (status > 0) call printerror(status)
	1405
	1406
	1407	return
	1408	end subroutine dump_alms_KLOAD
	1409	!=======================================================================
	1410	subroutine write_alms_KLOAD(filename, nalms, alms, ncl, header, nlheader, extno)
	1411	!=======================================================================
	1412	! Writes alms from to binary FITS file, FKH/Apr-99
	1413	! ncl is the number of columns, in the output fits file,
	1414	! either 3 or 5 (with or without errors respectively)
	1415	!
	1416	! input array (real) FITS file
	1417	! alms(:,1) = l )
	1418	! alms(:,2) = m )---> col 1: l*l+l+m+1
	1419	! alms(:,3) = real(a_lm) ---> col 2
	1420	! alms(:,4) = imag(a_lm) ---> col 3
	1421	! alms(:,5) = real(delta a_lm) ---> col 4
	1422	! alms(:,6) = imag(delta a_lm) ---> col 5
	1423	!
	1424	!=======================================================================
	1425
	1426	INTEGER(I4B), INTENT(IN) :: nalms, nlheader, ncl, extno
	1427	REAL(KMAP), INTENT(IN), DIMENSION(0:nalms-1,1:(ncl+1)) :: alms
	1428	CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header
	1429	CHARACTER(LEN=*), INTENT(IN) :: filename
	1430
	1431	INTEGER(I4B), DIMENSION(:), allocatable :: lm
	1432	INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1)
	1433	INTEGER(I4B) :: i,hdutype, lmax, mmax, lmin
	1434	LOGICAL(LGT) :: simple,extend
	1435	CHARACTER(LEN=80) :: comment
	1436
	1437	INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	1438	INTEGER(I4B) :: nrows, npix, tfields, varidat, repeat
	1439	INTEGER(I4B) :: frow, felem, colnum, stride, istart, iend, k
	1440	CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
	1441	CHARACTER(LEN=10) :: card
	1442	CHARACTER(LEN=2) :: stn
	1443	INTEGER(I4B) :: itn
	1444	integer(I4B) :: l, m
	1445	character(len=filenamelen) sfilename
	1446	character(len=1) :: pform
	1447	!-----------------------------------------------------------------------
	1448
	1449	if (KMAP == SP) pform = 'E'
	1450	if (KMAP == DP) pform = 'D'
	1451
	1452	status=0
	1453	unit = 140
	1454
	1455	! create the new empty FITS file
	1456	blocksize=1
	1457
	1458	if (extno==1) then
	1459
	1460	call ftinit(unit,filename,blocksize,status)
	1461	if (status > 0) call fatal_error("Error while creating file " &
	1462	& //trim(filename) &
	1463	& //". Check path and/or access rights.")
	1464
	1465	! -----------------------------------------------------
	1466	! initialize parameters about the FITS image
	1467	simple=.true.
	1468	bitpix=32 ! integer*4
	1469	naxis=0 ! no image
	1470	naxes(1)=0
	1471	extend=.true. ! there is an extension
	1472
	1473	! ----------------------
	1474	! primary header
	1475	! ----------------------
	1476	! write the required header keywords
	1477	call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
	1478
	1479	! writes supplementary keywords : none
	1480
	1481	! write the current date
	1482	call ftpdat(unit,status) ! format ccyy-mm-dd
	1483
	1484	! ----------------------
	1485	! image : none
	1486	! ----------------------
	1487
	1488	! ----------------------
	1489	! extension
	1490	! ----------------------
	1491
	1492	else
	1493
	1494	!*********************************************
	1495	! reopen an existing file and go to the end
	1496	!*********************************************
	1497	! remove the leading '!' (if any) when reopening the same file
	1498	sfilename = adjustl(filename)
	1499	if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
	1500	call ftopen(unit,sfilename,1_i4b,blocksize,status)
	1501	call ftmahd(unit,extno,hdutype,status)
	1502
	1503	endif
	1504
	1505	! creates an extension
	1506	call ftcrhd(unit, status)
	1507
	1508	! writes required keywords
	1509	nrows = nalms ! naxis1
	1510	tfields = ncl
	1511	repeat = 1
	1512	tform(1)(1:2)='1J' ! necessary with ifort 11.1.076
	1513	do i=2, ncl
	1514	tform(i)(1:2) = '1'//pform
	1515	enddo
	1516	ttype(1)(1:15) = 'index=l^2+l+m+1'
	1517	ttype(2)(1:15) = 'alm (real) '
	1518	ttype(3)(1:15) = 'alm (imaginary)'
	1519	if (ncl>3) then
	1520	ttype(4)(1:17) = 'error (real) '
	1521	ttype(5)(1:17) = 'error (imaginary)'
	1522	endif
	1523	tunit(1:ncl) = '' ! optional, will not appear
	1524	extname = '' ! optional, will not appear
	1525	varidat = 0
	1526	call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
	1527	& extname, varidat, status)
	1528
	1529	! write the header literally, putting TFORMi at the desired place
	1530	do i=1,nlheader
	1531	card = header(i)
	1532	if (card(1:5) == 'TTYPE') then ! if TTYPEi is explicitely given
	1533	stn = card(6:8)
	1534	read(stn,'(i3)') itn
	1535	! discard at their original location:
	1536	call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and ! remove
	1537	status = 0
	1538	call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
	1539	status = 0
	1540	if (itn <= tfields) then ! only put relevant information 2008-08-27
	1541	call putrec(unit,header(i), status) ! write new TTYPE1
	1542	status = 0
	1543	if (itn==1) then
	1544	comment = 'data format of field: 4-byte INTEGER'
	1545	else
	1546	if (KMAP == SP) comment = 'data format of field: 4-byte REAL'
	1547	if (KMAP == DP) comment = 'data format of field: 8-byte REAL'
	1548	endif
	1549	call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after
	1550	endif
	1551	elseif (header(i)/=' ') then
	1552	call putrec(unit,header(i), status)
	1553	endif
	1554	status = 0
	1555	enddo
	1556
	1557	lmax = nint(maxval(alms(:,1)))
	1558	lmin = nint(minval(alms(:,1)))
	1559	mmax = nint(maxval(alms(:,2)))
	1560	call ftukyj(unit, 'MIN-LPOL', lmin, 'Minimum L multipole order', status)
	1561	call ftukyj(unit, 'MAX-LPOL', lmax, 'Maximum L multipole order', status)
	1562	call ftukyj(unit, 'MAX-MPOL', mmax, 'Maximum M multipole degree', status)
	1563
	1564	! write the extension by blocks of rows ! EH, Dec 2004
	1565	felem = 1 ! starting position (element)
	1566	!!! stride = 1000 ! 1000 rows at a time
	1567	call ftgrsz(unit, stride, status) ! find optimal stride in rows
	1568	stride = max( stride, 1)
	1569	allocate(lm(0:stride-1))
	1570	do k = 0, (nalms-1)/(stride * repeat)
	1571	istart = k * (stride * repeat)
	1572	iend = min(nalms, istart + stride * repeat) - 1
	1573	do i = istart, iend ! recode the (l,m) -> lm mapping, EH, 2002-08-08
	1574	l = nint(alms(i,1))
	1575	m = nint(alms(i,2))
	1576	lm(i-istart) = (l+1)*l + m + 1
	1577	enddo
	1578
	1579	frow = istart/repeat + 1
	1580	npix = iend - istart + 1
	1581	call ftpclj( unit, 1_i4b, frow, felem, npix, lm, status)
	1582	do colnum = 2, ncl
	1583	call f90ftpcl_(unit, colnum, frow, felem, npix, alms(istart:iend,colnum+1), status)
	1584	enddo
	1585	enddo
	1586	deallocate(lm)
	1587
	1588	! ----------------------
	1589	! close and exit
	1590	! ----------------------
	1591
	1592	! close the file and free the unit number
	1593	call ftclos(unit, status)
	1594
	1595	! check for any error, and if so print out error messages
	1596	if (status > 0) call printerror(status)
	1597
	1598	return
	1599	end subroutine write_alms_KLOAD
	1600	!=======================================================================
	1601	subroutine read_alms_KLOAD(filename, nalms, alms, ncl, header, nlheader, extno)
	1602	!=======================================================================
	1603	! Read a FITS file
	1604	!
	1605	! slightly modified to deal with vector column
	1606	! in binary table EH/IAP/Jan-98
	1607	!
	1608	! Used by synfast when reading a binary file with alms for cons.real.
	1609	! FKH/Apr-99
	1610	!
	1611	! called by fits2alms
	1612	! Jan 2005, EH, improved for faster reading
	1613	!=======================================================================
	1614	CHARACTER(LEN=*), INTENT(IN) :: filename
	1615	INTEGER(I4B), INTENT(IN) :: nalms, ncl,nlheader,extno
	1616	REAL(KMAP), DIMENSION(0:nalms-1,1:(ncl+1)), INTENT(OUT) :: alms
	1617	CHARACTER(LEN=80), INTENT(OUT), DIMENSION(1:nlheader) :: header
	1618	REAL(KMAP) :: nullval
	1619	LOGICAL(LGT) :: anynull
	1620
	1621	INTEGER(I4B), DIMENSION(:), allocatable :: lm
	1622	INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis
	1623	INTEGER(I4B) :: npix
	1624	CHARACTER(LEN=80) :: comment ! , record
	1625	LOGICAL(LGT) :: extend
	1626	INTEGER(I4B) :: nmove, hdutype ! , nkeys , nspace
	1627	INTEGER(I4B) :: frow, imap
	1628	INTEGER(I4B) :: datacode, repeat, width
	1629	integer(I4B) :: i, l, m
	1630	integer(i4b) :: nrow2read, nelem
	1631	integer(i8b) :: i0, i1
	1632
	1633	INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	1634	INTEGER(I4B) :: nrows, tfields, varidat
	1635	CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
	1636
	1637	!-----------------------------------------------------------------------
	1638	status=0
	1639	header=''
	1640	comment=''
	1641	ttype=''
	1642	tform=''
	1643	tunit=''
	1644	extname=''
	1645	unit = 139
	1646	naxes(1) = 1
	1647	naxes(2) = 1
	1648	nfound = -1
	1649	anynull = .false.
	1650	alms=0.
	1651	readwrite=0
	1652	call ftopen(unit,filename,readwrite,blocksize,status)
	1653	if (status > 0) then
	1654	call printerror(status)
	1655	call fatal_error("Aborting.")
	1656	endif
	1657	! -----------------------------------------
	1658
	1659	! determines the presence of image
	1660	call ftgkyj(unit,'NAXIS', naxis, comment, status)
	1661	if (status > 0) call printerror(status)
	1662
	1663	! determines the presence of an extension
	1664	call ftgkyl(unit,'EXTEND', extend, comment, status)
	1665	if (status > 0) status = 0 ! no extension :
	1666	! to be compatible with first version of the code
	1667
	1668	call assert (extend, 'No extension!')
	1669	nmove = +extno
	1670	call ftmrhd(unit, nmove, hdutype, status)
	1671	!cc write(,) hdutype
	1672
	1673	call assert(hdutype==2, 'this is not a binary table')
	1674
	1675	header = ""
	1676	call get_clean_header( unit, header, filename, status)
	1677
	1678	! reads all the keywords
	1679	call ftghbn(unit, MAXDIM, &
	1680	& nrows, tfields, ttype, tform, tunit, extname, varidat, &
	1681	& status)
	1682
	1683	if (tfields<ncl) then
	1684	print *,'found ',tfields,' columns in the file'
	1685	print *,'expected ',ncl
	1686	call fatal_error
	1687	endif
	1688	! finds the bad data value
	1689	! if (KMAP == SP) call ftgkye(unit,'BAD_DATA',nullval,comment,status)
	1690	! if (KMAP == DP) call ftgkyd(unit,'BAD_DATA',nullval,comment,status)
	1691	call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
	1692	if (status == 202) then ! bad_data not found
	1693	if (KMAP == SP) nullval = s_bad_value ! default value
	1694	if (KMAP == DP) nullval = d_bad_value ! default value
	1695	status = 0
	1696	endif
	1697	!parse TFORM keyword to find out the length of the column vector
	1698	call ftbnfm(tform(1), datacode, repeat, width, status)
	1699	npix = nrows * repeat
	1700	! if (npix /= nalms) then
	1701	if (npix > nalms) then
	1702	print *,'found ',npix,' alms'
	1703	! print *,'expected ',nalms
	1704	print *,'expected ',nalms,' or less'
	1705	call fatal_error
	1706	endif
	1707
	1708	call ftgrsz(unit, nrow2read, status)
	1709	nrow2read = max(nrow2read, 1)
	1710	nelem = nrow2read * repeat
	1711	i0 = 0_i8b
	1712	allocate(lm(0:nelem-1))
	1713	do frow = 1, nrows, nrow2read
	1714	i1 = min(i0 + nrow2read * repeat, int(npix,i8b)) - 1_i8b
	1715	nelem = i1 - i0 + 1
	1716	! first column -> index
	1717	call ftgcvj(unit, 1_i4b, frow, 1_i4b, nelem, i_bad_value, &
	1718	& lm(0), anynull, status)
	1719	call assert(.not. anynull, 'There are undefined values in the table!')
	1720	! other columns -> a(l,m)
	1721	do imap = 2, ncl
	1722	call f90ftgcv_(unit, imap, frow, 1_i4b, nelem, nullval, &
	1723	& alms(i0:i1,imap+1), anynull, status)
	1724	call assert (.not. anynull, 'There are undefined values in the table!')
	1725	enddo
	1726	! recoding of the mapping, EH, 2002-08-08
	1727	! lm = l*l + l + m + 1
	1728	do i = i0, i1
	1729	l = int(sqrt( real(lm(i-i0)-1, kind = DP) ) )
	1730	m = lm(i-i0) - l*(l+1) - 1
	1731	! check that round-off did not screw up mapping
	1732	if (abs(m) > l) then
	1733	print*,'Inconsistent l^2+l+m+1 -> l,m mapping'
	1734	print,l, m, l(l+1)+m+1, lm(i-i0)
	1735	call fatal_error
	1736	endif
	1737	alms(i,1) = real( l, kind=KMAP)
	1738	alms(i,2) = real( m, kind=KMAP)
	1739	enddo
	1740	i0 = i1 + 1_i8b
	1741	enddo
	1742	deallocate(lm)
	1743	! sanity check
	1744	if (i0 /= npix) then
	1745	print*,'something wrong during piece-wise reading'
	1746	call fatal_error
	1747	endif
	1748
	1749	! close the file
	1750	call ftclos(unit, status)
	1751
	1752
	1753	! check for any error, and if so print out error messages
	1754	if (status > 0) call printerror(status)
	1755	return
	1756	end subroutine read_alms_KLOAD
	1757	!**************************************************************************
	1758	SUBROUTINE read_bintod_KLOAD(filename, tod, npixtot, ntods, firstpix, nullval, anynull, &
	1759	header, extno, units)
	1760	!**************************************************************************
	1761	!=======================================================================
	1762	! Read a FITS file
	1763	!
	1764	! slightly modified to deal with vector column (ie TFORMi = '1024E')
	1765	! in binary table EH/IAP/Jan-98
	1766	!
	1767	! This routine is used for reading TODS by anafast.
	1768	! Modified to start at a given pix numb OD & RT 02/02
	1769	! Modified to handle huge array (npix_tot > 2^32) OD & EH 07/02
	1770	! 2002-07-08 : bugs correction by E.H.
	1771	! 2007-01-31 : added units optional output
	1772	!=======================================================================
	1773
	1774	IMPLICIT NONE
	1775
	1776	CHARACTER(LEN=*), INTENT(IN) :: filename
	1777	INTEGER(I8B) , INTENT(IN) :: npixtot,firstpix
	1778	INTEGER(I4B), INTENT(IN) :: ntods
	1779	REAL(KMAP), DIMENSION(0:,1:), INTENT(OUT) :: tod
	1780	REAL(KMAP), INTENT(OUT) :: nullval
	1781	LOGICAL(LGT), INTENT(OUT) :: anynull
	1782	character(len=*), dimension(1:),intent(out), optional :: header
	1783	INTEGER(I4B), INTENT(IN), OPTIONAL :: extno
	1784	character(LEN=*), dimension(1:),intent(OUT), optional :: units
	1785
	1786
	1787	integer(I4B) :: i, nl_header, len_header, nl_units, len_units
	1788	INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis
	1789	INTEGER(I4B) :: npix_32 !,firstpix_32
	1790	CHARACTER(LEN=80) :: comment
	1791	LOGICAL(LGT) :: extend
	1792	INTEGER(I4B) :: nmove, hdutype
	1793	INTEGER(I4B) :: column, frow, itod
	1794	INTEGER(I4B) :: datacode, repeat, width
	1795
	1796	INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	1797	INTEGER(I4B) :: nrows, tfields, varidat,felem
	1798	CHARACTER(LEN=20), dimension(1:MAXDIM) :: ttype, tform,tunit
	1799	CHARACTER(LEN=20) :: extname
	1800
	1801	INTEGER(I8B) :: q,iq,npix_tmp,firstpix_tmp, i0, i1
	1802
	1803	!-----------------------------------------------------------------------
	1804	status=0
	1805
	1806	unit = 138
	1807	naxes(1) = 1
	1808	naxes(2) = 1
	1809	nfound = -1
	1810	anynull = .FALSE.
	1811
	1812	nl_header = 0
	1813	if (present(header)) then
	1814	nl_header = size(header)
	1815	len_header = 80
	1816	endif
	1817
	1818	nl_units = 0
	1819	if (present(units)) then
	1820	nl_units = size(units)
	1821	len_units = min(80,len(units(1))) ! due to SUN compiler bug
	1822	endif
	1823
	1824	readwrite=0
	1825	CALL ftopen(unit,filename,readwrite,blocksize,status)
	1826	IF (status .GT. 0) CALL printerror(status)
	1827	! -----------------------------------------
	1828
	1829	! determines the presence of image
	1830	CALL ftgkyj(unit,'NAXIS', naxis, comment, status)
	1831	IF (status .GT. 0) CALL printerror(status)
	1832
	1833	! determines the presence of an extension
	1834	CALL ftgkyl(unit,'EXTEND', extend, comment, status)
	1835	IF (status .GT. 0) status = 0 ! no extension :
	1836	! to be compatible with first version of the code
	1837
	1838	IF (naxis .GT. 0) THEN ! there is an image
	1839	print*,'WARNING : Image is ignored in '//trim(filename)
	1840	ENDIF
	1841	IF (extend) THEN ! there is an extension
	1842
	1843	nmove = +1
	1844	if (present(extno)) nmove = +1 + extno
	1845	CALL ftmrhd(unit, nmove, hdutype, status)
	1846	!cc write(,) hdutype
	1847
	1848	call assert(hdutype==2, 'this is not a binary table')
	1849
	1850	! reads all the keywords
	1851	CALL ftghbn(unit, MAXDIM, &
	1852	& nrows, tfields, ttype, tform, tunit, extname, varidat, &
	1853	& status)
	1854
	1855	IF (tfields .LT. ntods) THEN
	1856	PRINT *,'found ',tfields,' tods in the file'
	1857	PRINT *,'expected ',ntods
	1858	call fatal_error
	1859	ENDIF
	1860
	1861	if (nl_header > 0) then
	1862	header = ""
	1863	status = 0
	1864	call get_clean_header(unit, header, filename, status)
	1865	endif
	1866
	1867	if (nl_units > 0) then
	1868	do i=1,nl_units
	1869	units(i)(1:len_units) = 'unknown' ! default
	1870	enddo
	1871	do itod = 1, min(ntods, nl_units)
	1872	units(itod)(1:len_units) = adjustl(tunit(itod))
	1873	enddo
	1874	endif
	1875
	1876	! finds the bad data value
	1877	! if (KMAP == SP) CALL ftgkye(unit,'BAD_DATA',nullval,comment,status)
	1878	! if (KMAP == DP) CALL ftgkyd(unit,'BAD_DATA',nullval,comment,status)
	1879	call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
	1880	IF (status .EQ. 202) THEN ! bad_data not found
	1881	if (KMAP == SP) nullval = s_bad_value ! default value
	1882	if (KMAP == DP) nullval = d_bad_value ! default value
	1883	status = 0
	1884	ENDIF
	1885
	1886	IF (npixtot .LT. nchunk_max) THEN
	1887
	1888	DO itod = 1, ntods
	1889
	1890	!parse TFORM keyword to find out the length of the column vector (repeat)
	1891	CALL ftbnfm(tform(itod), datacode, repeat, width, status)
	1892	frow = (firstpix)/repeat+1 ! 1 based
	1893	felem = firstpix-(frow-1)*repeat+1 ! 1 based
	1894	npix_32 = npixtot
	1895
	1896	!reads the columns
	1897	column = itod
	1898	CALL f90ftgcv_(unit, column, frow, felem, npix_32, nullval, &
	1899	& tod(0:npix_32-1,itod), anynull, status)
	1900	END DO
	1901
	1902	ELSE
	1903
	1904	q = (npixtot-1)/nchunk_max
	1905	DO iq = 0,q
	1906	IF (iq .LT. q) THEN
	1907	npix_tmp = nchunk_max
	1908	ELSE
	1909	npix_tmp = npixtot - iq*nchunk_max
	1910	ENDIF
	1911	firstpix_tmp = firstpix + iq*nchunk_max
	1912	npix_32 = npix_tmp
	1913	i0 = firstpix_tmp-firstpix
	1914	i1 = i0 + npix_tmp - 1_i8b
	1915
	1916	DO itod = 1, ntods
	1917	! parse TFORM keyword to find out the length of the column vector
	1918	CALL ftbnfm(tform(itod), datacode, repeat, width, status)
	1919	frow = (firstpix_tmp)/repeat+1 ! 1 based
	1920	felem = firstpix_tmp-(frow-1)*repeat+1 ! 1 based
	1921	CALL f90ftgcv_(unit, itod, frow, felem, npix_32, nullval, &
	1922	& tod(i0:i1,itod), anynull, status)
	1923	END DO
	1924
	1925	ENDDO
	1926
	1927	ENDIF
	1928
	1929	ELSE ! no image no extension, you are dead, man
	1930	call fatal_error(' No image, no extension')
	1931	ENDIF
	1932
	1933	! close the file
	1934	CALL ftclos(unit, status)
	1935
	1936	! check for any error, and if so print out error messages
	1937	IF (status .GT. 0) CALL printerror(status)
	1938
	1939	RETURN
	1940
	1941	END SUBROUTINE read_bintod_KLOAD
	1942	!=======================================================================
	1943
	1944	!======================================================================================
	1945	SUBROUTINE write_bintabh_KLOAD(tod, npix, ntod, header, nlheader, filename, extno, firstpix, repeat)
	1946	!======================================================================================
	1947
	1948	! =================================================================================
	1949	! Create a FITS file containing a binary table extension in the first extension
	1950	!
	1951	! Designed to deal with Huge file, (n_elements > 2^31)
	1952	!
	1953	! OPTIONAL NEW PARAMETERS:
	1954	! firstpix : position in the file of the first element to be written (starts at 0)
	1955	! default value =0
	1956	! 8-byte integer
	1957	! if NE 0 then suppose that the file already exists
	1958	!
	1959	! repeat : length of vector per unit rows and columns of the first binary extension
	1960	! default value = 12000 (\equiv 1 mn of PLANCK/HFI data)
	1961	! 4-byte integer
	1962	!
	1963	! OTHER PARAMETERS
	1964	! unchanged with respect to the standard write_bintab of the HEALPIX package except
	1965	! npix which is an 8-byte integer
	1966	!
	1967	! Adapted from write_bintab
	1968	! E.H. & O.D. @ IAP 07/02
	1969	!
	1970	! Requires a compilation in 64 bits of the CFITSIO
	1971	! Note that the flag -D_FILE_OFFSETS_BITS=64 has to be added
	1972	! (cf page CFITIO 2.2 User's guide Chap 4, section 4-13)
	1973	!
	1974	! 2002-07-08 : bugs correction by E.H.
	1975	! (uniform use of firstpix_tmp, introduction of firstpix_chunk)
	1976	! 2015-04-28: better handling of repeat (suggested by R. Keskitalo)
	1977	!==========================================================================================
	1978
	1979	USE healpix_types
	1980	IMPLICIT NONE
	1981
	1982	INTEGER(I8B) , INTENT(IN) :: npix
	1983	INTEGER(I8B) , INTENT(IN), OPTIONAL :: firstpix
	1984	INTEGER(I4B) , INTENT(IN), OPTIONAL :: repeat
	1985	INTEGER(I4B) , INTENT(IN) :: ntod,nlheader
	1986	REAL(KMAP) , INTENT(IN), DIMENSION(0:npix-1,1:ntod) :: tod
	1987	CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header
	1988	CHARACTER(LEN=*), INTENT(IN) :: filename
	1989	INTEGER(I4B) , INTENT(IN) , OPTIONAL :: extno
	1990
	1991	INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1),repeat_fits
	1992	INTEGER(I4B) :: i,npix_32
	1993	LOGICAL(LGT) :: simple,extend
	1994	CHARACTER(LEN=80) :: comment, ch
	1995	integer(I8B) :: repeat_tmp
	1996
	1997	INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	1998	INTEGER(I4B) :: nrows,tfields,varidat
	1999	INTEGER(I4B) :: frow,felem,colnum,readwrite,width,datacode,hdutype
	2000	CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
	2001	CHARACTER(LEN=10) :: card
	2002	CHARACTER(LEN=2) :: stn
	2003	INTEGER(I4B) :: itn
	2004
	2005	real(KMAP), dimension(:), allocatable :: padding
	2006	integer(I8B) :: lastpix
	2007	integer(I4B) :: lengap
	2008
	2009	INTEGER(I4B) :: extno_i
	2010	character(len=filenamelen) :: sfilename
	2011	INTEGER(I8B) :: q,iq,npix_tmp,firstpix_tmp,firstpix_chunk, i0, i1
	2012	character(len=1) :: pform
	2013	character(len=*), parameter :: code='write_bintabh'
	2014	!-----------------------------------------------------------------------
	2015
	2016	if (KMAP == SP) pform = 'E'
	2017	if (KMAP == DP) pform = 'D'
	2018
	2019	IF (.NOT. PRESENT(repeat) ) THEN
	2020	repeat_tmp = 1
	2021	if (mod(npix,1024_i8b) == 0) then
	2022	repeat_tmp = 1024
	2023	elseif (npix >= 12000) then
	2024	repeat_tmp = 12000
	2025	endif
	2026	ELSE
	2027	repeat_tmp = repeat
	2028	ENDIF
	2029	IF (.NOT. PRESENT(firstpix) ) THEN
	2030	firstpix_tmp = 0
	2031	ELSE
	2032	firstpix_tmp = firstpix
	2033	ENDIF
	2034
	2035	extno_i = 0
	2036	if (present(extno)) extno_i = extno
	2037
	2038	status=0
	2039	unit = 137
	2040	blocksize=1
	2041
	2042	! remove the leading '!' (if any) when reopening the same file
	2043	sfilename = adjustl(filename)
	2044	if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
	2045
	2046	! create the new empty FITS file
	2047
	2048	IF (firstpix_tmp .EQ. 0) THEN
	2049
	2050	if (extno_i == 0) then
	2051	CALL ftinit(unit,filename,blocksize,status)
	2052	if (status > 0) call fatal_error("Error while creating file " &
	2053	& //trim(filename) &
	2054	& //". Check path and/or access rights.")
	2055
	2056	! -----------------------------------------------------
	2057	! Initialize parameters about the FITS image
	2058	simple=.TRUE.
	2059	bitpix=32 ! integer*4
	2060	naxis=0 ! no image
	2061	naxes(1)=0
	2062	extend=.TRUE. ! there is an extension
	2063
	2064	! ----------------------
	2065	! primary header
	2066	! ----------------------
	2067	! write the required header keywords
	2068	CALL ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
	2069
	2070	! writes supplementary keywords : none
	2071
	2072	! write the current date
	2073	CALL ftpdat(unit,status) ! format ccyy-mm-dd
	2074
	2075	! ----------------------
	2076	! image : none
	2077	! ----------------------
	2078
	2079	! ----------------------
	2080	! extension
	2081	! ----------------------
	2082	else
	2083
	2084	!*********************************************
	2085	! reopen an existing file and go to the end
	2086	!*********************************************
	2087	call ftopen(unit,sfilename,1_i4b,blocksize,status)
	2088	call ftmahd(unit,1_i4b+extno_i,hdutype,status)
	2089
	2090	endif
	2091
	2092	! creates an extension
	2093	CALL ftcrhd(unit, status)
	2094
	2095	! writes required keywords
	2096	! if (npix < repeat_tmp) repeat_tmp = npix ! 2015-04-28
	2097
	2098	! nrows = npix / repeat_tmp + 1 ! naxis1
	2099	nrows = (npix - 1_i8b) / repeat_tmp + 1_i4b ! 2015-04-28
	2100	tfields = ntod
	2101	WRITE(ch,'(i8)') repeat_tmp
	2102	! tform(1:ntod) = TRIM(ADJUSTL(ch))//pform ! does not work with Ifort, EH, 2006-04-04
	2103	ch = TRIM(ADJUSTL(ch))//pform
	2104	tform(1:ntod) = ch
	2105
	2106	! IF (npix .LT. repeat_tmp) THEN ! 2015-04-28
	2107	! nrows = npix
	2108	! ch = '1'//pform
	2109	! tform(1:ntod) = ch
	2110	! ENDIF
	2111	ttype(1:ntod) = 'simulation' ! will be updated
	2112	tunit(1:ntod) = '' ! optional, will not appear
	2113
	2114	extname = '' ! optional, will not appear
	2115	varidat = 0
	2116
	2117	CALL ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
	2118	& extname, varidat, status)
	2119
	2120	! write the header literally, putting TFORM1 at the desired place
	2121	if (KMAP == SP) comment = 'data format of field: 4-byte REAL'
	2122	if (KMAP == DP) comment = 'data format of field: 8-byte REAL'
	2123	DO i=1,nlheader
	2124	card = header(i)
	2125	IF (card(1:5) == 'TTYPE') THEN ! if TTYPE1 is explicitely given
	2126	stn = card(6:8)
	2127	READ(stn,'(i3)') itn
	2128	! discard at their original location:
	2129	call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and
	2130	status = 0
	2131	call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
	2132	status = 0
	2133	if (itn <= tfields) then ! only put relevant information 2008-08-27
	2134	call putrec(unit,header(i), status) ! write new TTYPE1
	2135	status = 0
	2136	call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after
	2137	!CALL ftmcrd(unit,'TTYPE'//stn,'COMMENT',status) ! old TTYPEi and
	2138	!CALL ftmcrd(unit,'TFORM'//stn,'COMMENT',status) ! TFORMi
	2139	!CALL ftprec(unit,header(i), status) ! write new TTYPE1
	2140	!CALL ftpkys(unit,'TFORM'//stn,tform(1),comment,status) ! and write new TFORM1 right after
	2141	endif
	2142	ELSEIF (header(i).NE.' ') THEN
	2143	call putrec(unit,header(i), status)
	2144	! CALL ftprec(unit,header(i), status)
	2145	ENDIF
	2146	status = 0
	2147	ENDDO
	2148
	2149	ELSE
	2150	! The file already exists
	2151	readwrite=1
	2152	CALL ftopen(unit,sfilename,readwrite,blocksize,status)
	2153	CALL ftmahd(unit,2_i4b+extno_i,hdutype,status)
	2154
	2155	CALL ftgkys(unit,'TFORM1',tform(1),comment,status)
	2156	CALL ftbnfm(tform(1),datacode,repeat_fits,width,status)
	2157
	2158	IF (repeat_tmp .NE. repeat_fits) THEN
	2159	if (present(repeat)) then
	2160	write(*,'(a)') code//'> WARNING: In file '//trim(sfilename)
	2161	write(*,'(a)') &
	2162	& code//'> user provided REPEAT value (' // &
	2163	& trim(adjustl(string(repeat_tmp))) // &
	2164	& ') differs from value read from file (' // &
	2165	& trim(adjustl(string(repeat_fits))) // &
	2166	& ').'
	2167	write(*,'(a)') code//'> The latter will be used.'
	2168	! else
	2169	! write(*,'(a,i0.0,a)') &
	2170	! & code//'> WARNING: REPEAT value read from file (', &
	2171	! & repeat_fits, &
	2172	! & ') will be used.'
	2173	endif
	2174	repeat_tmp = repeat_fits
	2175	ENDIF
	2176
	2177	ENDIF
	2178
	2179
	2180	IF (npix .LT. nchunk_max) THEN ! data is small enough to be written in one chunk
	2181
	2182	frow = (firstpix_tmp)/repeat_tmp + 1
	2183	felem = firstpix_tmp-(frow-1)*repeat_tmp + 1
	2184	npix_32 = npix
	2185
	2186	DO colnum = 1, ntod
	2187	call f90ftpcl_(unit, colnum, frow, felem, npix_32, tod(0:npix_32-1,colnum), status)
	2188	END DO
	2189
	2190	ELSE ! data has to be written in several chunks
	2191
	2192	q = (npix-1)/nchunk_max
	2193	DO iq = 0,q
	2194	IF (iq .LT. q) THEN
	2195	npix_tmp = nchunk_max
	2196	ELSE
	2197	npix_tmp = npix - iq*nchunk_max
	2198	ENDIF
	2199	i0 = iq*nchunk_max
	2200	i1 = i0 + npix_tmp - 1_i8b
	2201	firstpix_chunk = firstpix_tmp + i0
	2202	frow = (firstpix_chunk)/repeat_tmp + 1
	2203	felem = firstpix_chunk-(frow-1)*repeat_tmp + 1
	2204	npix_32 = npix_tmp
	2205	DO colnum = 1, ntod
	2206	call f90ftpcl_(unit, colnum, frow, felem, npix_32, &
	2207	& tod(i0:i1, colnum), status)
	2208	END DO
	2209	ENDDO
	2210
	2211	ENDIF
	2212
	2213	! pad entry if necessary ! 2015-04-28
	2214	lastpix = firstpix_tmp + npix ! number of pixels written above
	2215	lengap = modulo(lastpix, repeat_tmp) ! remaining gap
	2216	if (lengap > 0) then
	2217	firstpix_tmp = lastpix
	2218	npix_32 = repeat_tmp - lengap
	2219	frow = (firstpix_tmp)/repeat_tmp + 1
	2220	felem = firstpix_tmp-(frow-1)*repeat_tmp + 1
	2221	allocate(padding(0:npix_32-1))
	2222	if (KMAP == SP) padding(:) = HPX_SBADVAL
	2223	if (KMAP == DP) padding(:) = HPX_DBADVAL
	2224	do colnum = 1, ntod
	2225	call f90ftpcl_(unit, colnum, frow, felem, npix_32, padding, status)
	2226	enddo
	2227	deallocate(padding)
	2228	endif
	2229
	2230	! ----------------------
	2231	! close and exit
	2232	! ----------------------
	2233
	2234	! close the file and free the unit number
	2235	CALL ftclos(unit, status)
	2236
	2237	! check for any error, and if so print out error messages
	2238	IF (status .GT. 0) CALL printerror(status)
	2239
	2240	RETURN
	2241
	2242	END SUBROUTINE write_bintabh_KLOAD
	2243	! ==============================================================================
	2244
	2245	!******************************************************************************
	2246	!---------------------------------------
	2247	! reads a FITS file containing a list of
	2248	! ring-based or pixel-based quadrature weights
	2249	! and turns it into a (ring-ordered) full sky map
	2250	!
	2251	! adapted from unfold_weights.pro
	2252	! 2018-05-25
	2253	!---------------------------------------
	2254	subroutine unfold_weightsfile_KLOAD(w8file, w8map)
	2255
	2256	USE pix_tools, only: nside2npix, nside2npweights
	2257
	2258	character(len=*), intent(in) :: w8file
	2259	real(KMAP), intent(out), dimension(0:) :: w8map
	2260	real(KMAP), allocatable, dimension(:,:) :: w8list
	2261	integer(i4b) :: nside, status, won
	2262	integer(i8b) :: nfw8, npw8, nrw8
	2263	character(len=*), parameter :: code = 'unfold_weightsfile'
	2264
	2265	nfw8 = getsize_fits(w8file, nside=nside)
	2266	npw8 = nside2npweights(nside)
	2267	nrw8 = 2*nside
	2268
	2269	won = 0
	2270	if (nfw8 == nrw8) won = 1
	2271	if (nfw8 == npw8) won = 2
	2272	if (won == 0) then
	2273	print*,'Weights file = '//trim(w8file)
	2274	print*,'contains ',nfw8,' weights'
	2275	print*,'for Nside = ',nside
	2276	print*,'Expected either ',nrw8,' or ',npw8
	2277	call fatal_error
	2278	endif
	2279
	2280	allocate(w8list(0:nfw8-1,1:1), stat=status)
	2281	call assert_alloc(status,code,"w8list")
	2282	call input_map(w8file, w8list, nfw8, 1)
	2283	call unfold_weightslist(nside, won, w8list(:,1), w8map)
	2284	deallocate(w8list)
	2285
	2286	return
	2287	end subroutine unfold_weightsfile_KLOAD
	2288	!---------------------------------------
	2289	! turns a list of ring-based or pixel-based quadrature weights
	2290	! into a (ring-ordered) full sky map
	2291	!
	2292	! adapted from unfold_weights.pro
	2293	! 2018-05-18
	2294	! 2018-09-12: promoted qpix to I8B to avoid consistency problems during compilation with g95
	2295	!---------------------------------------
	2296	subroutine unfold_weightslist_KLOAD(nside, won, w8list, w8map)
	2297
	2298	USE pix_tools, only: nside2npix, nside2npweights
	2299	USE long_intrinsic, only: long_size
	2300
	2301	integer(i4b), intent(in) :: nside, won
	2302	real(KMAP), dimension(0:), intent(in) :: w8list
	2303	real(KMAP), dimension(0:), intent(out) :: w8map
	2304
	2305	integer(i8b) :: np, npix
	2306	integer(i8b) :: nf, nw8, pnorth, psouth, vpix, qp4, p, qpix
	2307	integer(i4b) :: ring, odd, shifted, rpix, wpix, j4, n4, it
	2308
	2309	! test nside
	2310	npix = nside2npix(nside)
	2311	if (npix <= 0) then
	2312	print*,'Unvalid Nside = ',nside
	2313	call fatal_error
	2314	endif
	2315
	2316	! test won
	2317	if (won < 1 .or. won > 2) then
	2318	print*,'Expected either 1 or 2, got ',won
	2319	call fatal_error
	2320	endif
	2321
	2322	! test input list size
	2323	nf = long_size(w8list)
	2324	if (won == 1) then
	2325	nw8 = 2 * nside
	2326	else
	2327	nw8 = nside2npweights(nside)
	2328	endif
	2329	if (nf /= nw8) then
	2330	print*,'Expected in input weight list',nw8
	2331	print*,'got ',nf
	2332	call fatal_error
	2333	endif
	2334
	2335	! test output map size
	2336	np = long_size(w8map)
	2337	if (np /= npix) then
	2338	print*,'Expected in output weight map',npix
	2339	print*,'got ',np
	2340	call fatal_error
	2341	endif
	2342
	2343	! do actual unfolding
	2344	if (won == 1) then
	2345
	2346	pnorth = 0_i8b
	2347	n4 = 4 * nside
	2348	do ring=1,n4-1 ! loop on all rings
	2349	it = min(ring, n4 - ring) ! ring index, counting from closest pole
	2350	qp4= 4_i8b * min(it, nside) ! pixels / ring
	2351	w8map(pnorth:pnorth+qp4-1) = w8list(it-1)
	2352	pnorth = pnorth + qp4
	2353	enddo
	2354
	2355	else
	2356
	2357	pnorth = 0_i8b ! position in expanded weights
	2358	vpix = 0_i8b ! position in compress list
	2359	do ring=0, 2*nside-1
	2360	qpix = min(ring+1, nside) ! 1/4 of number of pixels per ring
	2361	shifted = 0
	2362	if (ring < (nside-1) .or. iand(ring+nside,1_i4b) == 1_i4b) shifted = 1
	2363	odd = iand(qpix,1_i8b)
	2364	qp4 = 4_i8b*qpix ! number of pixels per ring
	2365	! fill the weight map
	2366	do p=0,qp4-1
	2367	j4 = mod(p, qpix) ! g95: p and qpix must be of same type
	2368	rpix = min(j4, qpix - shifted - j4) ! seesaw
	2369	w8map(pnorth+p) = w8list(vpix + rpix)
	2370	enddo
	2371
	2372	if (ring < 2*nside-1) then ! south part (except equator)
	2373	psouth = npix - pnorth - qp4
	2374	w8map(psouth:psouth+qp4-1) = w8map(pnorth:pnorth+qp4-1)
	2375	endif
	2376	! locations on next ring
	2377	pnorth = pnorth + qp4
	2378	wpix = (qpix+1)/2 + 1 - ior(odd, shifted)
	2379	vpix = vpix + wpix
	2380	enddo
	2381	endif
	2382
	2383	! add 1 to get final weight
	2384	w8map = w8map + 1.0_KMAP
	2385
	2386	return
	2387	end subroutine unfold_weightslist_KLOAD
	2388
	2389	!=======================================================================
	2390	! read_fits_partial(filename, pixel, cutmap, [header, extno])
	2391	! routine to read FITS file with cut sky data : PIXEL, ?, ?, ?, ?, ...
	2392	!=======================================================================
	2393	#ifndef NO64BITS
	2394	!=======================================================================
	2395	subroutine read_fits_partial4_KLOAD(filename, pixel, cutmap, header, extno)
	2396	!=======================================================================
	2397	use long_intrinsic, only: long_size
	2398	character(len=*), intent(in) :: filename
	2399	integer(I4B), dimension(0:), intent(out) :: pixel
	2400	real(KMAP), dimension(0:,1:), intent(out) :: cutmap
	2401	character(len=*), dimension(1:), intent(out), optional :: header
	2402	integer(I4B), intent(in), optional :: extno
	2403	!
	2404	integer(i8b), dimension(:), allocatable :: pixel8
	2405	integer(i8b) :: npix
	2406	character(len=*), parameter :: code = 'read_fits_partial'
	2407	!---------------------------------------------------------------------
	2408	npix = long_size(pixel)
	2409	if (npix > MAX_I4b) then
	2410	call fatal_error(code//': PIXEL must be of type I8B to read '//trim(filename))
	2411	endif
	2412	allocate(pixel8(0:npix-1))
	2413	call read_fits_partial8_KLOAD(filename, pixel8, cutmap, header, extno)
	2414	if (maxval(pixel8) > MAX_I4B) then
	2415	call fatal_error(code//': PIXEL must be of type I8B to read '//trim(filename))
	2416	endif
	2417	pixel(0:npix-1) = pixel8(0:npix-1)
	2418	deallocate(pixel8)
	2419	return
	2420
	2421	end subroutine read_fits_partial4_KLOAD
	2422	#endif
	2423	!=======================================================================
	2424	subroutine read_fits_partial8_KLOAD(filename, pixel, cutmap, header, extno)
	2425	!=======================================================================
	2426	use long_intrinsic, only: long_size
	2427	character(len=*), intent(in) :: filename
	2428	integer(I8B), dimension(0:), intent(out) :: pixel
	2429	real(KMAP), dimension(0:,1:), intent(out) :: cutmap
	2430	character(len=*), dimension(1:), intent(out), optional :: header
	2431	integer(I4B), intent(in), optional :: extno
	2432
	2433	integer(I8B) :: obs_npix, npixtot, npix
	2434
	2435	integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	2436	integer(I4B) :: blocksize, datacode
	2437	integer(I4B) :: firstpix, frow, hdutype
	2438	integer(I4B) :: naxis, nfound, nmove, nrows, nmaps, nmr, i, npix_4
	2439	integer(I4B) :: readwrite, status, tfields, unit, varidat, width
	2440	integer(I4B) :: repeat, repeat1, repeat2
	2441	logical(LGT) :: anynull, extend
	2442	character(len=20), dimension(MAXDIM) :: ttype, tform, tunit
	2443	character(len=20) :: extname
	2444	character(len=80) :: comment
	2445	real(KMAP) :: nullval
	2446	!=======================================================================
	2447
	2448	npixtot = min(long_size(pixel),long_size(cutmap,1))
	2449	nmaps = size(cutmap,2)
	2450	status=0
	2451
	2452	unit = 150
	2453	nfound = -1
	2454	anynull = .false.
	2455
	2456	readwrite=0
	2457	call ftopen(unit,filename,readwrite,blocksize,status)
	2458	if (status > 0) call printerror(status)
	2459
	2460	! determines the presence of image
	2461	call ftgkyj(unit,'NAXIS', naxis, comment, status)
	2462	if (status > 0) call printerror(status)
	2463	if (naxis > 0) then ! there is an image
	2464	print*,'an image was found in the FITS file '//filename
	2465	print*,'... it is ignored.'
	2466	endif
	2467
	2468	! determines the presence of an extension
	2469	call ftgkyl(unit,'EXTEND', extend, comment, status)
	2470	if (status > 0) then
	2471	print*,'extension expected and not found in FITS file '//filename
	2472	print*,'abort code'
	2473	call fatal_error
	2474	endif
	2475
	2476	nmove = +1
	2477	if (present(extno)) nmove = +1 + extno
	2478	call ftmrhd(unit, nmove, hdutype, status)
	2479
	2480	call assert (hdutype==2, 'this is not a binary table')
	2481
	2482	! reads all the FITS related keywords
	2483	call ftghbn(unit, MAXDIM, &
	2484	& nrows, tfields, ttype, tform, tunit, extname, varidat, &
	2485	& status)
	2486
	2487	if (.not. (trim(ttype(1)) == 'PIXEL' &
	2488	& .or. trim(ttype(1)) == 'PIX' ) ) call fatal_error('did not find PIXEL column in '//filename)
	2489
	2490
	2491	if (present(header)) then
	2492	header = ""
	2493	status = 0
	2494	!call fitstools_mp_get_clean_header(unit, header, filename, status)
	2495	call get_clean_header(unit, header, filename, status)
	2496	endif
	2497
	2498	! if (present(units)) then
	2499	! ! the second column contains the SIGNAL, for which we
	2500	! ! already read the units from the header
	2501	! units = adjustl(tunit(2))
	2502	! endif
	2503
	2504	! finds the bad data value
	2505	call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
	2506	if (status == 202) then ! bad_data not found
	2507	if (KMAP == SP) nullval = s_bad_value ! default value
	2508	if (KMAP == DP) nullval = d_bad_value ! default value
	2509	status = 0
	2510	endif
	2511
	2512	frow = 1
	2513	firstpix = 1
	2514	! parse TFORM keyword to find out the length of the column vector
	2515	repeat1 = 1
	2516	repeat2 = 1
	2517	call ftbnfm(tform(1), datacode, repeat1, width, status)
	2518	if (tfields > 1) call ftbnfm(tform(2), datacode, repeat2, width, status)
	2519	repeat = max(repeat1, repeat2)
	2520
	2521	!call ftgkyj(unit,'OBS_NPIX',obs_npix,comment,status) ! I4B
	2522	!call ftgkyk(unit,'OBS_NPIX',obs_npix,comment,status) ! I8B
	2523	call f90ftgky_(unit,'OBS_NPIX',obs_npix,comment,status) ! generic
	2524	if (status == 202) then ! obs_npix not found
	2525	obs_npix = int(nrows, kind=i8b) * repeat
	2526	status = 0
	2527	endif
	2528
	2529	npix = min(npixtot, obs_npix)
	2530	npix_4 = npix
	2531	nmr = min(nmaps, tfields-1)
	2532	!call ftgcvj (unit, 1_i4b, frow, firstpix, npix_4, i_bad_value, pixel(0), anynull, status) ! I4B
	2533	!call ftgcvk (unit, 1_i4b, frow, firstpix, npix_4, l_bad_value, pixel(0), anynull, status) ! I8B
	2534	call f90ftgcv_(unit, 1_i4b, frow, firstpix, npix_4, l_bad_value, pixel, anynull, status) ! generic
	2535	do i=1,nmr
	2536	call f90ftgcv_(unit, i+1_i4b, frow, firstpix, npix_4, nullval, cutmap(0:npix-1,i), anynull, status)
	2537	enddo
	2538
	2539	! close the file
	2540	call ftclos(unit, status)
	2541
	2542	! check for any error, and if so print out error messages
	2543	if (status > 0) call printerror(status)
	2544
	2545	return
	2546	end subroutine read_fits_partial8_KLOAD
	2547
	2548
	2549	!=======================================================================
	2550	! writes a fits file for (polarized) sky data set with columns:
	2551	! PIXEL, TEMPERATURE, [Q_POLARISATION, U_POLARISATION]
	2552	!
	2553	! write_fits_partial(filename, pixel, cutmap, &
	2554	! [, header, coord, nside, order, units, extno])
	2555	!=======================================================================
	2556	subroutine write_fits_partial4_KLOAD(filename, pixel, cutmap, &
	2557	& header, coord, nside, order, units, extno)
	2558	!=======================================================================
	2559	character(len=*), intent(in) :: filename
	2560	integer(I4B), dimension(0:), intent(in) :: pixel
	2561	real(KMAP), dimension(0:,1:), intent(in) :: cutmap
	2562	character(len=*), dimension(1:), intent(in), optional :: header
	2563	integer(I4B), intent(in), optional :: nside, order
	2564	character(len=*), intent(in), optional :: coord, units
	2565	integer(I4B), intent(in), optional :: extno !, polarisation
	2566	character(len=*), parameter :: routine = 'write_fits_partial'
	2567	call sub_write_fits_partial_KLOAD(filename, cutmap, pixel4=pixel, &
	2568	header=header, coord=coord, nside=nside, order=order, units=units, extno=extno)
	2569	end subroutine write_fits_partial4_KLOAD
	2570	!=======================================================================
	2571	#ifndef NO64BITS
	2572	subroutine write_fits_partial8_KLOAD(filename, pixel, cutmap, &
	2573	& header, coord, nside, order, units, extno)
	2574	!=======================================================================
	2575	character(len=*), intent(in) :: filename
	2576	integer(I8B), dimension(0:), intent(in) :: pixel
	2577	real(KMAP), dimension(0:,1:), intent(in) :: cutmap
	2578	character(len=*), dimension(1:), intent(in), optional :: header
	2579	integer(I4B), intent(in), optional :: nside, order
	2580	character(len=*), intent(in), optional :: coord, units
	2581	integer(I4B), intent(in), optional :: extno !, polarisation
	2582	character(len=*), parameter :: routine = 'write_fits_partial'
	2583	call sub_write_fits_partial_KLOAD(filename, cutmap, pixel8=pixel, &
	2584	header=header, coord=coord, nside=nside, order=order, units=units, extno=extno)
	2585	end subroutine write_fits_partial8_KLOAD
	2586	#endif
	2587	!=======================================================================
	2588	subroutine sub_write_fits_partial_KLOAD(filename, cutmap, &
	2589	& pixel4, pixel8, &
	2590	& header, coord, nside, order, units, extno)
	2591	!=======================================================================
	2592	use long_intrinsic, only: long_size
	2593	use pix_tools, only: nside2npix
	2594	use misc_utils, only: assert
	2595
	2596	character(len=*), intent(in) :: filename
	2597	real(KMAP), dimension(0:,1:), intent(in) :: cutmap
	2598	integer(I4B), dimension(0:), intent(in), optional :: pixel4
	2599	integer(I8B), dimension(0:), intent(in), optional :: pixel8
	2600	character(len=*), dimension(1:), intent(in), optional :: header
	2601	integer(I4B), intent(in), optional :: nside, order
	2602	character(len=*), intent(in), optional :: coord, units
	2603	integer(I4B), intent(in), optional :: extno !, polarisation
	2604
	2605	character(len=*), parameter :: routine = 'write_fits_partial'
	2606	! --- healpix/fits related variables ----
	2607	integer(I4B) :: ncol, grain, nd2
	2608	integer(I4B) :: npix_hd, nside_hd, i
	2609	integer(I4B) :: maxpix, minpix
	2610	character(len=1) :: char1, coord_usr
	2611	! character(len=8) :: char8
	2612	character(len=20) :: units_usr
	2613	logical(LGT) :: done_nside, done_order, done_coord, done_polar, polar_flag
	2614	integer(I4B) :: nlheader, extno_i
	2615
	2616	! --- cfitsio related variables ----
	2617	integer(I4B) :: status, unit, blocksize, bitpix, naxis, naxes(1)
	2618	logical(LGT) :: simple, extend
	2619	character(LEN=80) :: svalue, comment
	2620
	2621	integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension
	2622	integer(I4B) :: nrows, tfields, varidat
	2623	integer(I4B) :: frow, felem, repeat, repeatg, hdutype
	2624	character(len=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
	2625	character(len=20), dimension(1:3) :: extnames
	2626	character(len=80) :: card
	2627	character(len=4) :: srepeat, srepeatg
	2628	character(len=1) :: pform1, pform2
	2629	character(len=filenamelen) sfilename
	2630	integer(I8B) :: obs_npix8, npix_final, obs_npp
	2631	integer(I4B) :: obs_npix4, kpix, nside_final
	2632
	2633	integer(I4B), save :: nside_old
	2634	!=======================================================================
	2635	if (KMAP == SP) pform2='E'
	2636	if (KMAP == DP) pform2='D'
	2637	if (present(pixel4) .eqv. present(pixel8)) then
	2638	call fatal_error(routine//': choose either PIXEL4 or PIXEL8')
	2639	endif
	2640	if (present(pixel4)) then
	2641	kpix = I4B
	2642	pform1='J'
	2643	obs_npp = long_size(pixel4)
	2644	elseif (present(pixel8)) then
	2645	kpix = I8B
	2646	pform1='K'
	2647	obs_npp = long_size(pixel8)
	2648	endif
	2649	obs_npix8 = long_size(cutmap, 1)
	2650	nd2 = size(cutmap, 2)
	2651	obs_npix4 = int(obs_npix8, kind=I4B)
	2652	call assert( obs_npix8 == obs_npp, routine//': mismatched size for PIXEL and DATA')
	2653	ncol = 1 + nd2
	2654	grain = 1
	2655	status=0
	2656	unit = 149
	2657	blocksize=1
	2658
	2659	nlheader = 0
	2660	if (present(header)) nlheader = size(header)
	2661	units_usr = ' '
	2662	if (present(units)) units_usr = units
	2663	extno_i = 0
	2664	if (present(extno)) extno_i = extno
	2665	polar_flag = .false.
	2666	done_polar = .false.
	2667	! if (present(polarisation)) then
	2668	! polar_flag = (polarisation /= 0)
	2669	! done_polar = .true. ! will ignore the header provided value of POLAR
	2670	! endif
	2671
	2672	if (extno_i == 0) then
	2673	!*************************************
	2674	! create the new empty FITS file
	2675	!*************************************
	2676	call ftinit(unit,filename,blocksize,status)
	2677	if (status > 0) call fatal_error("Error while creating file " &
	2678	& //trim(filename) &
	2679	& //". Check path and/or access rights.")
	2680
	2681	! -----------------------------------------------------
	2682	! initialize parameters about the FITS image
	2683	simple=.true.
	2684	bitpix=32 ! integer*4
	2685	naxis=0 ! no image
	2686	naxes(1)=0
	2687	extend=.true. ! there is an extension
	2688
	2689	! ----------------------
	2690	! primary header
	2691	! ----------------------
	2692	! write the required header keywords
	2693	call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
	2694
	2695	! write the current date
	2696	call ftpdat(unit,status) ! format ccyy-mm-dd
	2697
	2698	else
	2699	!*********************************************
	2700	! reopen an existing file and go to the end
	2701	!*********************************************
	2702	! remove the leading '!' (if any) when reopening the same file
	2703	sfilename = adjustl(filename)
	2704	if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
	2705	call ftopen(unit,sfilename,1_i4b,blocksize,status)
	2706	! move to last extension written
	2707	call ftmahd(unit, 1_i4b+ extno_i, hdutype, status)
	2708
	2709	endif
	2710
	2711
	2712	! ----------------------
	2713	! new extension
	2714	! ----------------------
	2715	call ftcrhd(unit, status)
	2716
	2717	! writes required keywords
	2718	! repeat = 1024
	2719	repeat = 1
	2720	if (obs_npix8 < repeat) repeat = 1
	2721	nrows = (obs_npix8 + repeat - 1)/ repeat ! naxis1
	2722	write(srepeat,'(i4)') repeat
	2723	srepeat = adjustl(srepeat)
	2724
	2725	repeatg = repeat * grain
	2726	write(srepeatg,'(i4)') repeatg
	2727	srepeatg = adjustl(srepeatg)
	2728
	2729	tfields = ncol
	2730	ttype(1) = 'PIXEL'
	2731	if (nd2 == 1) then
	2732	done_polar = .true.
	2733	polar_flag = .false.
	2734	ttype(2) = 'TEMPERATURE '
	2735	elseif (nd2 == 3) then
	2736	done_polar = .true.
	2737	polar_flag = .true.
	2738	ttype(2:4) = (/'TEMPERATURE ' , &
	2739	& 'Q_POLARISATION' , &
	2740	& 'U_POLARISATION' /)
	2741	else
	2742	do i=2,ncol
	2743	ttype(i) = 'C'//string(i-1, format='(i2.2)') ! C01, C02, ...
	2744	enddo
	2745	endif
	2746
	2747	tform(1) = trim(srepeat) //pform1
	2748	tform(2:ncol) = trim(srepeatg)//pform2
	2749
	2750	tunit = ' ' ! optional, will not appear
	2751	tunit(2:ncol) = units_usr
	2752
	2753	extname = 'SKY_OBSERVATION' ! default, will be overridden by user provided one if any
	2754	! if (polar_flag) extname = extnames(1+extno_i)
	2755	varidat = 0
	2756	call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
	2757	& extname, varidat, status)
	2758
	2759	call ftpcom(unit,'------------------------------------------',status)
	2760	call ftpcom(unit,' Pixelisation Specific Keywords ',status)
	2761	call ftpcom(unit,'------------------------------------------',status)
	2762	call ftpkys(unit,'PIXTYPE','HEALPIX ',' HEALPIX Pixelisation',status)
	2763	call ftpkyu(unit,'NSIDE', ' ',status) ! place holder, will be updated later on
	2764	call ftpkyu(unit,'ORDERING',' ',status)
	2765	call ftpkys(unit,'COORDSYS',' ',' ',status)
	2766	call ftpcom(unit,' G = Galactic, E = ecliptic, C = celestial = equatorial', status)
	2767	call ftpkyl(unit,'POLAR',polar_flag,'Polarisation included in file (T/F)',status)
	2768	call ftpcom(unit,'------------------------------------------',status)
	2769	call ftpcom(unit,' Data Specific Keywords ',status)
	2770	call ftpcom(unit,'------------------------------------------',status)
	2771	call ftpkys(unit,'INDXSCHM','EXPLICIT',' Indexing : IMPLICIT or EXPLICIT', status)
	2772	! call ftpkyj(unit,'GRAIN', grain, ' Grain of pixel indexing',status)
	2773	! call ftpcom(unit,'GRAIN=0 : no indexing of pixel data (IMPLICIT)',status)
	2774	! call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)',status)
	2775	! call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
	2776	call ftpkys(unit,'OBJECT','PARTIAL ',' Sky coverage represented by data',status)
	2777	if (KPIX == I4B) then
	2778	call ftpkyj(unit,'OBS_NPIX',obs_npix4, ' Number of pixels observed and recorded',status)
	2779	else
	2780	call ftpkyk(unit,'OBS_NPIX',obs_npix8, ' Number of pixels observed and recorded',status)
	2781	endif
	2782
	2783	! add required Healpix keywords (NSIDE, ORDER) if provided by user
	2784	done_order = .false.
	2785	if (present(order)) then
	2786	! char8 = order
	2787	! call ftupch(char8)
	2788	! if (char8(1:4) == 'RING') then
	2789	if (order == 1) then
	2790	call ftukys(unit, 'ORDERING','RING',' Pixel ordering scheme, either RING or NESTED',status)
	2791	done_order = .true.
	2792	! elseif (char8(1:4) == 'NEST') then
	2793	elseif (order == 2) then
	2794	call ftukys(unit, 'ORDERING','NESTED',' Pixel ordering scheme, either RING or NESTED ',status)
	2795	done_order = .true.
	2796	else
	2797	print*,'Invalid ORDER given : ',order, ' instead of 1 (RING) or 2 (NESTED)'
	2798	endif
	2799	endif
	2800
	2801	done_nside = .false.
	2802	if (present(nside)) then
	2803	if (nside2npix(nside) > 0) then ! valid nside
	2804	call ftukyj(unit,'NSIDE',nside,' Resolution parameter for HEALPIX', status)
	2805	done_nside = .true.
	2806	nside_final = nside
	2807	else
	2808	print*,'Invalid NSIDE given : ',nside
	2809	endif
	2810	endif
	2811
	2812	! add non required Healpix keyword (COORD)
	2813	done_coord = .false.
	2814	if (present(coord)) then
	2815	coord_usr = adjustl(coord)
	2816	char1 = coord_usr(1:1)
	2817	call ftupch(char1) ! uppercase
	2818	if (char1 == 'C' .or. char1 == 'Q') then
	2819	coord_usr = 'C'
	2820	done_coord = .true.
	2821	elseif (char1 == 'G') then
	2822	coord_usr = 'G'
	2823	done_coord = .true.
	2824	elseif (char1 == 'E' ) then
	2825	coord_usr='E'
	2826	done_coord = .true.
	2827	else
	2828	print*,'Unrecognised COORD given : ',coord,' instead of C, G, or E'
	2829	print*,'Proceed at you own risks '
	2830	coord_usr = char1
	2831	done_coord = .true.
	2832	endif
	2833	if (done_coord) then
	2834	call ftukys(unit, 'COORDSYS',coord_usr,' Pixelisation coordinate system',status)
	2835	endif
	2836	endif
	2837
	2838
	2839	! write the user provided header literally, except for PIXTYPE, TFORM, TTYPE, TUNIT*, INDXSCHM and GRAIN
	2840	! copy NSIDE, ORDERING and COORDSYS and POLAR if they are valid and not already given
	2841	do i=1,nlheader
	2842	card = header(i)
	2843	if (card(1:5) == 'TTYPE' .or. card(1:5) == 'TFORM' .or. card(1:7) == 'PIXTYPE') then
	2844	continue ! don't keep them
	2845	else if (card(1:8) == 'INDXSCHM') then
	2846	continue
	2847	else if (card(1:5) == 'GRAIN') then ! already written above
	2848	continue
	2849	else if (card(1:13) == 'COMMENT GRAIN' .or. card(1:14) == 'COMMENT GRAIN') then ! already written above
	2850	continue
	2851	else if (card(1:5) == 'TUNIT') then
	2852	if (trim(units_usr) == '') then
	2853	call ftucrd(unit,'TUNIT'//card(6:7),card, status) !update TUNIT2 and above
	2854	endif
	2855	else if (card(1:5) == 'NSIDE') then
	2856	call ftpsvc(card, svalue, comment, status)
	2857	read(svalue,*) nside_hd
	2858	npix_hd = nside2npix(nside_hd)
	2859	if (.not. done_nside .and. npix_hd > 0) then
	2860	call ftucrd(unit,'NSIDE',card, status) !update NSIDE
	2861	done_nside = .true.
	2862	nside_final = nside_hd
	2863	endif
	2864	else if (card(1:8) == 'ORDERING') then
	2865	call ftpsvc(card, svalue, comment, status)
	2866	svalue = adjustl(svalue)
	2867	svalue = svalue(2:8) ! remove leading quote
	2868	call ftupch(svalue)
	2869	if (.not. done_order .and. (svalue(1:4)=='RING' .or. svalue(1:6) == 'NESTED')) then
	2870	call ftucrd(unit,'ORDERING',card, status) !update ORDERING
	2871	done_order = .true.
	2872	endif
	2873	else if (card(1:8) == 'COORDSYS') then
	2874	if (.not. done_coord) call putrec(unit,card, status)
	2875	done_coord = .true.
	2876	else if (card(1:5) == 'POLAR') then
	2877	if (.not. done_polar) then
	2878	call ftucrd(unit,'POLAR', card, status) ! update POLAR
	2879	done_polar = .true.
	2880	endif
	2881	else if (card(1:7) == 'EXTNAME') then
	2882	call ftucrd(unit, 'EXTNAME', card, status)
	2883	else if (card(1:1) /= ' ') then ! if none of the above, copy to FITS file
	2884	call putrec(unit, card, status)
	2885	endif
	2886	10 continue
	2887	enddo
	2888
	2889
	2890	! test that required keywords have been provided in some way
	2891	if (.not. done_nside) then
	2892	print*,routine//'> NSIDE is a Healpix required keyword, '
	2893	print*,routine//'> it was NOT provided either as routine argument or in the input header'
	2894	print*,routine//'> abort execution, file not written'
	2895	call fatal_error
	2896	endif
	2897	if (.not. done_order) then
	2898	print*,routine//'> ORDER is a Healpix required keyword, '
	2899	print*,routine//'> it was NOT provided either as routine argument or in the input header'
	2900	print*,routine//'> abort execution, file not written'
	2901	call fatal_error
	2902	endif
	2903	! if ((.not. done_polar) .and. extno_i >=2) then
	2904	! print*,routine//'> Warning: POLAR keyword not set while 3 extensions have been written'
	2905	! endif
	2906
	2907	! check that NSIDE is the same for all extensions
	2908	if (extno_i == 0) then
	2909	nside_old = nside_final
	2910	else
	2911	if (nside_final /= nside_old) then
	2912	print*,routine//'> Inconsistent NSIDE: ',nside_final, nside_old
	2913	print*,routine//'> Should use same NSIDE for all extensions'
	2914	endif
	2915	endif
	2916
	2917	! check validity of PIXEL
	2918	npix_final = nside2npix(nside_final)
	2919	if (kpix == I4B) then
	2920	minpix = minval(pixel4(0:obs_npix8-1))
	2921	maxpix = maxval(pixel4(0:obs_npix8-1))
	2922	else
	2923	minpix = minval(pixel8(0:obs_npix8-1))
	2924	maxpix = maxval(pixel8(0:obs_npix8-1))
	2925	endif
	2926	if (minpix < 0 .or. maxpix > npix_final-1) then
	2927	print*,routine//'> Actual pixel range = ',minpix,maxpix
	2928	print*,routine//'> expected range (for Nside =',nside_final,') : 0, ',npix_final-1
	2929	print*,routine//'> ABORT execution, file not written.'
	2930	call fatal_error
	2931	endif
	2932	if (obs_npix8 > npix_final) then
	2933	print*,routine//'> The actual number of pixels ',obs_npix8
	2934	print*,routine//'> is larger than the number of pixels over the whole sky : ',npix_final
	2935	print*,routine//'> for Nside = ',nside_final
	2936	print*,routine//'> ABORT execution, file not written.'
	2937	call fatal_error
	2938	endif
	2939
	2940	! write the extension one column by one column
	2941	frow = 1 ! starting position (row)
	2942	felem = 1 ! starting position (element)
	2943	if (kpix == I4B) then
	2944	call f90ftpcl_(unit, 1_i4b, frow, felem, obs_npix4, pixel4, status)
	2945	else
	2946	call f90ftpcl_(unit, 1_i4b, frow, felem, obs_npix4, pixel8, status)
	2947	endif
	2948	do i=1, nd2
	2949	call f90ftpcl_(unit, i+1_i4b, frow, felem, obs_npix4, cutmap(0:,i), status)
	2950	enddo
	2951
	2952	! ----------------------
	2953	! close and exit
	2954	! ----------------------
	2955
	2956	! close the file and free the unit number
	2957	call ftclos(unit, status)
	2958
	2959	! check for any error, and if so print out error messages
	2960	if (status > 0) call printerror(status)
	2961
	2962	return
	2963	end subroutine sub_write_fits_partial_KLOAD

-2359

~~src/f90/mod/fits_template.f90~~ less more

0		!-----------------------------------------------------------------------------
1		!
2		! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon,
3		! Benjamin D. Wandelt, Anthony J. Banday,
4		! Matthias Bartelmann, Hans K. Eriksen,
5		! Frode K. Hansen, Martin Reinecke
6		!
7		!
8		! This file is part of HEALPix.
9		!
10		! HEALPix is free software; you can redistribute it and/or modify
11		! it under the terms of the GNU General Public License as published by
12		! the Free Software Foundation; either version 2 of the License, or
13		! (at your option) any later version.
14		!
15		! HEALPix is distributed in the hope that it will be useful,
16		! but WITHOUT ANY WARRANTY; without even the implied warranty of
17		! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18		! GNU General Public License for more details.
19		!
20		! You should have received a copy of the GNU General Public License
21		! along with HEALPix; if not, write to the Free Software
22		! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
23		!
24		! For more information about HEALPix see http://healpix.sourceforge.net
25		!
26		!-----------------------------------------------------------------------------
27		! template for routine SP/DP overloading for module fitstools
28		!
29		! subroutine input_map_KLOAD (4/8)
30		! subroutine read_bintab_KLOAD (4/8)
31		! subroutine read_conbintab_KLOAD NotYet
32		! subroutine write_bintab_KLOAD (4/8)
33		! subroutine write_asctab_KLOAD NA
34		! subroutine dump_alms_KLOAD NotYet
35		! subroutine write_alms_KLOAD NotYet
36		! subroutine read_alms_KLOAD NotYet
37		! subroutine read_bintod_KLOAD (8)
38		! subroutine write_bintabh_KLOAD (8)
39		! subroutine unfold_weights_KLOAD
40		! subroutine unfold_weightslist_KLOAD
41		!
42		!
43		! K M A P : map kind either SP or DP
44		!
45		! K L O A D : suffixe of routine name, to be replaced by either s (SP) or d (DP)
46		!
47		! edited Jan 11, 2006 to deal correctly with polarised alms in write_alms and alms2fits
48		! edited Apr 04, 2006 to avoid concatenation problem in TFORMs (write_plm and write_bintabh with Ifort 9)
49		! edited Dec 20, 2006 to accept alm file with less elements than expected (in particular if vanishing ones are not included)
50		! 2007-05-10 : increased maxdim (max number of columns per extension) from 20 to 40
51		! 2008-08-27 : in dump_alms and write_alms and write_tab:
52		! do not write TTYPE# and TFORM# in excess of # of fields in the file
53		! 2008-10-14: corrected bug introduced in write_asctab
54		! 2012-02-23: correction of a possible bug with index writing in dump_alms and write_alms
55		! 2013-12-13: increased MAXDIM from 40 to MAXDIM_TOP
56		! 2018-05-22: added unfold_weights, unfold_weightslist
57		! 2019-10-14: can write TTYPE??? (ie with up to 3 digits) in write_* and dump_alms
58		!
59
60		!=======================================================================
61		! map_bad_pixels(map, fin, fout, nbads, verbose)
62		! map: input data (2D)
63		! fin: input value of 'bad' pixels
64		! fout: output value of same 'bad' pixels
65		! nbads: number of bad pixels found
66		! verbose: OPTIONAL, if set, be verbose
67		!=======================================================================
68
69		subroutine map_bad_pixels_KLOAD(map, fin, fout, nbads, verbose)
70		use long_intrinsic, only: long_size
71		real(KMAP), dimension(0:,1:), intent(inout) :: map
72		real(KMAP), intent(in) :: fin, fout
73		integer(I8B), intent(out) :: nbads
74		logical(LGT), optional, intent(in) :: verbose
75		integer(I8B) :: i, npix
76		integer(I8B), dimension(1:100) :: imiss
77		integer(I4B) :: imap, nmaps
78		logical(LGT) :: be_verbose
79		real(KMAP) :: threshold
80		!-----------------------------------------
81
82		npix = long_size(map, 1)
83		nmaps = long_size(map, 2)
84		threshold = abs(fin * 1.e-5_KMAP)
85
86		imiss(:) = 0
87		do imap = 1, nmaps
88		do i=0,npix-1
89		if ( abs( map(i,imap) - fin ) < threshold ) then
90		map(i,imap) = fout
91		imiss(imap) = imiss(imap)+1
92		endif
93		enddo
94		enddo
95		nbads = sum(imiss)
96
97		be_verbose = .false.
98		if (present(verbose)) be_verbose=verbose
99		if (be_verbose) then
100		write(*,'(a,1pe11.4)') 'blank value : ' ,fin
101		do imap = 1, nmaps
102		if (imiss(imap) > 0) then
103		write(*,'(i12,a,f7.3,a,1pe11.4)') &
104		& imiss(imap),' missing pixels (', &
105		& (100.0_KMAP*imiss(imap))/npix,' %),'// &
106		& ' have been set to : ',fout
107		endif
108		enddo
109		endif
110
111		return
112		end subroutine map_bad_pixels_KLOAD
113		!=======================================================================
114		! input_map
115		! reads fits file
116		! filename = fits file (input)
117		! map = data read from the file (ouput) = real*4 array of size (npixtot,nmaps)
118		! npixtot = length of the map (input)
119		! nmaps = number of maps
120		! fmissval = OPTIONAL argument (input) with the value to be given to missing
121		! pixels, its default value is 0
122		! header = OPTIONAL (output) contains extension header
123		! units = OPTIONAL (output) contains the maps units
124		! extno = OPTIONAL (input) contains the unit number to read from (0 based)
125		! ignore_polcconv = OPTIONAL (input), LGT, default=.false.
126		! take into account or not the POLCCONV FITS keyword
127		!
128		! modified Feb 03 for units argument to run with Sun compiler
129		! 2017-02-15: deals with POLCCONV
130		!=======================================================================
131		#ifndef NO64BITS
132		subroutine input_map4_KLOAD(filename, map, npixtot, nmaps, &
133		& fmissval, header, units, extno, ignore_polcconv)
134		!=======================================================================
135
136		INTEGER(I4B), INTENT(IN) :: npixtot
137		INTEGER(I4B), INTENT(IN) :: nmaps
138		REAL(KMAP), INTENT(OUT), dimension(0:,1:) :: map
139		CHARACTER(LEN=*), INTENT(IN) :: filename
140		REAL(KMAP), INTENT(IN), OPTIONAL :: fmissval
141		CHARACTER(LEN=*), INTENT(OUT), dimension(1:), OPTIONAL :: header
142		CHARACTER(LEN=*), INTENT(OUT), dimension(1:), OPTIONAL :: units
143		INTEGER(I4B), INTENT(IN) , optional :: extno
144		logical(LGT), intent(IN) , optional :: ignore_polcconv
145		integer(i8b) :: npixtot8
146
147		npixtot8 = npixtot
148		call input_map8_KLOAD(filename, map, npixtot8, nmaps, &
149		fmissval, header, units, extno, ignore_polcconv)
150
151		return
152		end subroutine input_map4_KLOAD
153		#endif
154
155		!=======================================================================
156		subroutine input_map8_KLOAD(filename, map, npixtot, nmaps, &
157		fmissval, header, units, extno, ignore_polcconv)
158		!=======================================================================
159		use head_fits, only: get_card, add_card
160		use pix_tools, only: nside2npix
161		INTEGER(I8B), INTENT(IN) :: npixtot
162		INTEGER(I4B), INTENT(IN) :: nmaps
163		REAL(KMAP), INTENT(OUT), dimension(0:,1:) :: map
164		CHARACTER(LEN=*), INTENT(IN) :: filename
165		REAL(KMAP), INTENT(IN), OPTIONAL :: fmissval
166		CHARACTER(LEN=*), INTENT(OUT), dimension(1:), OPTIONAL :: header
167		CHARACTER(LEN=*), INTENT(OUT), dimension(1:), OPTIONAL :: units
168		INTEGER(I4B), INTENT(IN) , optional :: extno
169		logical(LGT), intent(IN) , optional :: ignore_polcconv
170
171		INTEGER(I8B) :: i, imissing, obs_npix, maxpix, minpix
172		REAL(KMAP) :: fmissing, fmiss_effct
173		integer(I4B) :: imap, nlheader
174
175		LOGICAL(LGT) :: anynull, do_polcconv
176		! Note : read_fits_cut4 still SP and I4B only
177		integer(I4B), dimension(:), allocatable :: pixel
178		! real(KMAP), dimension(:), allocatable :: signal
179		real(SP), dimension(:), allocatable :: signal
180		integer(I4B) :: status
181		integer(I4B) :: type_fits, nmaps_fits
182		CHARACTER(LEN=80) :: units1
183		CHARACTER(LEN=80),dimension(1:30) :: unitsm
184		character(len=80), dimension(:), allocatable :: hbuffer
185		! CHARACTER(LEN=80),dimension(:), allocatable :: unitsm
186		integer(I4B) :: extno_i, extno_f
187		CHARACTER(LEN=80) :: polcconv
188		integer(I4B) :: ipolcconv, polar, nside_fits
189		character(len=*), parameter :: primer_url = 'http://healpix.sf.net/pdf/intro.pdf'
190		!-----------------------------------------------------------------------
191
192		units1 = ' '
193		unitsm(:) = ' '
194		fmiss_effct = 0.
195		imissing = 0
196		if (PRESENT(fmissval)) fmiss_effct = fmissval
197		if (PRESENT(header)) then
198		nlheader = size(header)
199		else
200		nlheader = 36*100
201		endif
202		extno_i = 0
203		if (present(extno)) extno_i = extno
204		allocate(hbuffer(1:nlheader))
205		do_polcconv = .true.
206		if (present(ignore_polcconv)) then
207		do_polcconv = .not. ignore_polcconv
208		endif
209
210		obs_npix = getsize_fits(filename, nmaps = nmaps_fits, type=type_fits, extno=extno_i, &
211		polcconv=ipolcconv, polarisation=polar, nside=nside_fits)
212
213		! if (nmaps_fits > nmaps) then
214		! print*,trim(filename)//' only contains ',nmaps_fits,' maps'
215		! print*,' You try to read ',nmaps
216		! endif
217		if (type_fits == 0 .or. type_fits == 2) then ! full sky map (in image or binary table)
218		call read_bintab(filename, map(0:,1:), &
219		& npixtot, nmaps, fmissing, anynull, header=hbuffer(1:), &
220		& units=unitsm(1:), extno=extno_i)
221		if (present(header)) then
222		do i=1,nlheader
223		header(i) = hbuffer(i)
224		enddo
225		endif
226		if (present(units)) then
227		units(1:size(units)) = unitsm(1:size(units))
228		endif
229		call map_bad_pixels(map, fmissing, fmiss_effct, imissing, verbose=.false.)
230		! do imap = 1, nmaps
231		! anynull = .true.
232		! if (anynull) then
233		! imissing = 0
234		! do i=0,npixtot-1
235		! if ( ABS(map(i,imap)/fmissing -1.) < 1.e-5 ) then
236		! map(i,imap) = fmiss_effct
237		! imissing = imissing + 1
238		! endif
239		! enddo
240		! endif
241		! enddo
242		if (imissing > 0) write(*,'(a,1pe11.4)') 'blank value : ' ,fmissing
243		else if (type_fits == 3 .and. (nmaps == 1 .or. nmaps == 3)) then
244		! cut sky FITS file, reading 1 map (I) or 3 maps (I,Q,U), each from a different extension
245		do imap = 1, nmaps
246		extno_f = extno_i + imap - 1 ! 2016-08-16
247		obs_npix = getsize_fits(filename, extno = extno_f)
248		! one partial map (in binary table with explicit pixel indexing)
249		allocate(pixel(0:obs_npix-1), stat = status)
250		allocate(signal(0:obs_npix-1), stat = status)
251		call read_fits_cut4(filename, int(obs_npix,kind=i4b), &
252		& pixel, signal, header=hbuffer, units=units1, extno=extno_f)
253		if (present(header) .and. imap == 1) then
254		do i=1,nlheader
255		header(i) = hbuffer(i)
256		enddo
257		endif
258		if (present(units)) units(imap) = trim(units1)
259		maxpix = maxval(pixel)
260		minpix = maxval(pixel)
261		if (maxpix > (npixtot-1) .or. minpix < 0) then
262		print*,'map constructed from file '//trim(filename)//', with pixels in ',minpix,maxpix
263		print*,' wont fit in array with ',npixtot,' elements'
264		call fatal_error
265		endif
266		map(:,imap) = fmiss_effct
267		map(pixel(:),imap) = signal(:)
268		imissing = npixtot - obs_npix
269		deallocate(signal)
270		deallocate(pixel)
271		enddo
272		else
273		print*,'Unable to read the ',nmaps,' required map(s) from file '//trim(filename)
274		print*,'file type = ',type_fits
275		call fatal_error
276		endif
277		!-----------------------------------------------------------------------
278		if (imissing > 0) then
279		write(*,'(i12,a,f7.3,a,1pe11.4)') &
280		& imissing,' missing pixels (', &
281		& (100.*imissing)/npixtot,' %),'// &
282		& ' have been set to : ',fmiss_effct
283		endif
284
285		! deallocate(unitsm)
286
287		!-----------------------------------------------------------------------
288		! deal with polcconv
289		! print,'****** in input_map ',nmaps, polar, type_fits, nside_fits
290
291		if ( do_polcconv .and. ( &
292		& (nmaps >= 3 .and. type_fits == 3) & ! cut-sky, multiple maps
293		& .or. &
294		& (nmaps >=3 .and. polar==1 .and. type_fits ==2 .and. &
295		& obs_npix==nside2npix(nside_fits)) & ! full-sky, polarized map
296		& ) ) then
297		if (ipolcconv == 0) then
298		print 9000,' The POLCCONV keyword was not found in '//trim(filename)
299		print 9000,' COSMO (HEALPix default) will be assumed, and map is unchanged.'
300		print 9000,' See HEALPix primer ('//primer_url//') for details.'
301		endif
302		! if (ipolcconv == 1) then
303		! print 9000,' POLCCONV=COSMO found in '//trim(filename)//'. Map is unchanged.'
304		! endif
305		if (ipolcconv == 2) then
306		print 9000,' POLCCONV=IAU found in '//trim(filename)
307		map(:,3) = - map(:,3)
308		if (present(header)) then
309		print 9000,' The sign of the U polarisation is changed in memory,'&
310		& //' and the header is updated.'
311		call add_card(header, 'POLCCONV', 'COSMO', &
312		comment='Changed by input_map', update=.true.)
313		else
314		print 9000,' The sign of the U polarisation is changed in memory.'
315		endif
316		print 9000,' See HEALPix primer ('//primer_url//') for details.'
317		endif
318		if (ipolcconv == 3) then
319		call get_card(hbuffer,'POLCCONV',polcconv)
320		print 9000,' POLCCONV='//trim(polcconv)//' found in '//trim(filename)
321		print 9000,' It is neither COSMO nor IAU. Aborting!'
322		print 9000,' See HEALPix primer ('//primer_url//') for details.'
323		call fatal_error
324		endif
325		endif
326		9000 format(a)
327		if (allocated(hbuffer)) deallocate(hbuffer)
328
329		RETURN
330		END subroutine input_map8_KLOAD
331		!=======================================================================
332		! Read a FITS file
333		! This routine is used for reading MAPS by anafast.
334		! modified Feb 03 for units argument to run with Sun compiler
335		! Jan 2005, EH, improved for faster writing
336		!=======================================================================
337		#ifndef NO64BITS
338		subroutine read_bintab4_KLOAD(filename, map, npixtot, nmaps, nullval, anynull, header, units, extno)
339		character(len=*), intent(IN) :: filename
340		integer(I4B), intent(IN) :: npixtot
341		integer(I4B), intent(IN) :: nmaps
342		real(KMAP), dimension(0:,1:), intent(OUT) :: map
343		real(KMAP), intent(OUT) :: nullval
344		logical(LGT), intent(OUT) :: anynull
345		character(LEN=*), dimension(1:), optional, intent(OUT) :: header
346		character(LEN=*), dimension(1:), optional, intent(OUT) :: units
347		integer(I4B) , optional, intent(IN) :: extno
348
349		integer(i8b):: npixtot8
350
351		npixtot8 = int(npixtot,kind=i8b)
352		call read_bintab8_KLOAD(filename, map, npixtot8, nmaps, nullval, anynull, header, units, extno)
353		return
354
355		end subroutine read_bintab4_KLOAD
356		#endif
357		subroutine read_bintab8_KLOAD(filename, map, npixtot, nmaps, nullval, anynull, header, units, extno)
358		!=======================================================================
359		character(len=*), intent(IN) :: filename
360		integer(I8B), intent(IN) :: npixtot
361		integer(I4B), intent(IN) :: nmaps
362		real(KMAP), dimension(0:,1:), intent(OUT) :: map
363		real(KMAP), intent(OUT) :: nullval
364		logical(LGT), intent(OUT) :: anynull
365		character(LEN=*), dimension(1:), optional, intent(OUT) :: header
366		character(LEN=*), dimension(1:), optional, intent(OUT) :: units
367		integer(I4B) , optional, intent(IN) :: extno
368
369		integer(I4B) :: nl_header, len_header, nl_units, len_units
370		integer(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis
371		integer(I4B) :: group, firstpix, i, npix32
372		real(KMAP) :: blank, testval
373		real(DP) :: bscale,bzero
374		character(len=80) :: comment
375		logical(LGT) :: extend
376		integer(I4B) :: nmove, hdutype, hdunum
377		integer(I4B) :: frow, imap
378		integer(I4B) :: datacode, width
379		LOGICAL(LGT) :: anynull_i
380
381		integer(I4B), parameter :: MAXDIM = MAXDIM_TOP !number of columns in the extension
382		integer(i8b) :: npix_old
383		integer(i8b), dimension(1:MAXDIM) :: npix
384		integer(i8b), dimension(1:MAXDIM) :: i0, i1
385		integer(i4b), dimension(1:MAXDIM) :: repeat
386		integer(i4b) :: nrow2read, nelem
387
388		integer(I4B) :: nrows, tfields, varidat
389		character(len=20), dimension(1:MAXDIM) :: ttype, tform, tunit
390		character(len=20) :: extname
391		character(len=*), parameter :: code='read_bintab'
392		!-----------------------------------------------------------------------
393		status=0
394
395		unit = 146
396		naxes(1) = 1
397		naxes(2) = 1
398		nfound = -1
399		anynull = .false.
400		bscale = 1.0d0
401		bzero = 0.0d0
402		blank = -2.e25
403		nullval = bscale*blank + bzero
404		comment=''
405		ttype=''
406		tform=''
407		tunit=''
408		extname=''
409
410		nl_header = 0
411		if (present(header)) then
412		nl_header = size(header)
413		len_header = 80
414		endif
415
416		nl_units = 0
417		if (present(units)) then
418		nl_units = size(units)
419		len_units = min(80,len(units(1))) ! due to SUN compiler bug
420		endif
421
422		readwrite=0
423		!call ftopen(unit,filename,readwrite,blocksize,status)
424		call ftnopn(unit,filename,readwrite, status) !open primary HDU or specified HDU
425		if (status > 0) call printerror(status)
426		! -----------------------------------------
427		call ftghdn(unit, hdunum)
428
429		if (hdunum == 1) then ! in primary HDU
430		! determines the presence of image
431		call ftgkyj(unit,'NAXIS', naxis, comment, status)
432		if (status > 0) call printerror(status)
433
434		! determines the presence of an extension
435		call ftgkyl(unit,'EXTEND', extend, comment, status)
436		if (status > 0) then
437		extend = .false.
438		status = 0 ! no extension :
439		! to be compatible with first version of the code
440		endif
441		endif
442
443		if (naxis > 0 .and. .not.extend .and. hdunum==1) then ! there is an image
444		! determine the size of the image (look naxis1 and naxis2)
445		call ftgknj(unit,'NAXIS',1_i4b,2_i4b,naxes,nfound,status)
446
447		! check that it found only NAXIS1
448		if (nfound == 2 .and. naxes(2) > 1) then
449		print *,'multi-dimensional image'
450		print *,'expected 1-D data.'
451		call fatal_error
452		end if
453
454		if (nfound < 1) then
455		call printerror(status)
456		print *,'can not find NAXIS1.'
457		call fatal_error
458		endif
459
460		npix(1)=naxes(1)
461		if (npix(1) /= npixtot) then
462		print *,'WARNING: found ',npix(1),' pixels in '//trim(filename)
463		print *,' expected ',npixtot
464		npix(1) = min(npix(1), npixtot)
465		print *,' only ',npix(1),' will be read'
466		endif
467
468		call ftgkyd(unit,'BSCALE',bscale,comment,status)
469		if (status == 202) then ! BSCALE not found
470		bscale = 1.0d0
471		status = 0
472		endif
473		call ftgkyd(unit,'BZERO', bzero, comment,status)
474		if (status == 202) then ! BZERO not found
475		bzero = 0.0d0
476		status = 0
477		endif
478		call f90ftgky_(unit, 'BLANK', blank, comment, status)
479		if (status == 202) then ! BLANK not found
480		! (according to fitsio BLANK is integer)
481		blank = -2.e25
482		status = 0
483		endif
484		nullval = bscale*blank + bzero
485
486		! -----------------------------------------
487
488		group=1
489		firstpix = 1
490		npix32 = npix(1)
491		call f90ftgpv_(unit, group, firstpix, npix32, nullval, map(0:npix(1)-1,1), anynull, status)
492		! if there are any NaN pixels, (real data)
493		! or BLANK pixels (integer data) they will take nullval value
494		! and anynull will switch to .true.
495		! otherwise, switch it by hand if necessary
496		testval = 1.e-6 * ABS(nullval)
497		do i=0, npix(1)-1
498		if (ABS(map(i,1)-nullval) < testval) then
499		anynull = .true.
500		goto 111
501		endif
502		enddo
503		111 continue
504
505		else if (extend .or. hdunum>1) then
506		if (hdunum == 1) then
507		nmove = +1
508		if (present(extno)) nmove = +1 + extno
509		call ftmrhd(unit, nmove, hdutype, status)
510		else
511		call ftghdt(unit, hdutype, status)
512		endif
513
514		call assert (hdutype==2, 'this is not a binary table')
515
516		! reads all the keywords
517		call ftghbn(unit, MAXDIM, &
518		& nrows, tfields, ttype, tform, tunit, extname, varidat, &
519		& status)
520
521		if (tfields < nmaps) then
522		print *,'found ',tfields,' maps in file '//trim(filename)
523		print *,'expected ',nmaps
524		call fatal_error
525		endif
526
527		! finds the bad data value
528		call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
529		if (status == 202) then ! bad_data not found
530		if (KMAP == SP) nullval = s_bad_value ! default value
531		if (KMAP == DP) nullval = d_bad_value ! default value
532		status = 0
533		endif
534
535		if (nl_header > 0) then
536		do i=1,nl_header
537		header(i)(1:len_header) = ""
538		enddo
539		call get_clean_header(unit, header, filename, status)
540		status = 0
541		endif
542
543		if (nl_units > 0) then
544		do i=1,nl_units
545		units(i)(1:len_units) = 'unknown' ! default
546		enddo
547		do imap = 1, min(nmaps, nl_units)
548		units(imap)(1:len_units) = adjustl(tunit(imap))
549		enddo
550		endif
551
552		npix_old = npixtot
553		do imap = 1, nmaps
554		!parse TFORM keyword to find out the length of the column vector
555		call ftbnfm(tform(imap), datacode, repeat(imap), width, status)
556		npix(imap) = int(nrows,i8b) * repeat(imap)
557		if (npix(imap) /= npixtot .and. npix_old /= npix(imap)) then
558		print *,'WARNING: found ',npix(imap),' pixels in '//trim(filename)//', column ',imap
559		print *,' expected ',npixtot,' or ',npix_old
560		npix_old = npix(imap)
561		npix(imap) = min(npix(imap), npixtot)
562		print *,' only ',npix(imap),' will be read'
563		endif
564		enddo
565
566
567		call ftgrsz(unit, nrow2read, status)
568		nrow2read = max(nrow2read, 1)
569		firstpix = 1 ! starting position in FITS within row, 1 based
570		i0(:) = 0_i8b ! starting element in array, 0 based
571		do frow = 1, nrows, nrow2read
572		do imap = 1, nmaps
573		i1(imap) = min(i0(imap) + int(nrow2read,i8b) * repeat(imap), npix(imap)) - 1_i8b
574		nelem = i1(imap) - i0(imap) + 1
575		call f90ftgcv_(unit, imap, frow, firstpix, nelem, &
576		& nullval, map(i0(imap):i1(imap),imap), anynull_i, status)
577		anynull = anynull .or. anynull_i
578		i0(imap) = i1(imap) + 1_i8b
579		enddo
580		enddo
581		! sanity check
582		do imap = 1, nmaps
583		if (i0(imap) /= npix(imap)) then
584		call fatal_error('something wrong during piece wise reading')
585		endif
586		enddo
587
588		else ! no image no extension, you are dead, man
589		call fatal_error(' No image, no extension in '//trim(filename))
590		endif
591		! close the file
592		call ftclos(unit, status)
593
594		! check for any error, and if so print out error messages
595		if (status > 0) call printerror(status)
596
597		return
598		end subroutine read_bintab8_KLOAD
599
600		!=======================================================================
601		subroutine read_conbintab_KLOAD(filename, alms, nalms, units, extno)
602		!=======================================================================
603		! Read a FITS file containing alms values
604		!
605		! slightly modified to deal with vector column
606		! in binary table EH/IAP/Jan-98
607		!
608		! Used by synfast when reading a binary file with alms for cons.real.
609		! FKH/Apr-99
610		!
611		! extno : optional, number of extension to be read, default=0: first extension
612		! Jan 2005, EH, improved for faster reading
613		!=======================================================================
614		CHARACTER(LEN=*), INTENT(IN) :: filename
615		INTEGER(I4B), INTENT(IN) :: nalms !, nlheader
616		REAL(KMAP), DIMENSION(0:nalms-1,1:6), INTENT(OUT) :: alms
617		CHARACTER(LEN=*), dimension(1:), optional, INTENT(OUT) :: units
618		INTEGER(I4B) , optional, INTENT(IN) :: extno
619
620		REAL(KMAP) :: nullval
621		LOGICAL(LGT) :: anynull
622
623		INTEGER(I4B), DIMENSION(:), allocatable :: lm
624		INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis
625		INTEGER(I4B) :: npix
626		CHARACTER(LEN=80) :: comment
627		LOGICAL(LGT) :: extend
628		INTEGER(I4B) :: nmove, hdutype
629		INTEGER(I4B) :: frow, imap
630		INTEGER(I4B) :: datacode, repeat, width
631		integer(I4B) :: i, l, m
632		integer(i4b) :: nrow2read, nelem
633		integer(i8b) :: i0, i1
634
635		INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
636		INTEGER(I4B) :: nrows, tfields, varidat
637		CHARACTER(LEN=20), dimension(1:MAXDIM) :: ttype, tform, tunit
638		CHARACTER(LEN=20) :: extname
639		character(len=*), parameter :: code="read_conbintab"
640
641		!-----------------------------------------------------------------------
642		status=0
643
644		unit = 145
645		naxes(1) = 1
646		naxes(2) = 1
647		nfound = -1
648		anynull = .false.
649		alms=0. ! set output to 0.
650		readwrite=0
651		comment=''
652		ttype=''
653		tform=''
654		tunit=''
655		extname=''
656		call ftopen(unit,filename,readwrite,blocksize,status)
657		if (status > 0) call printerror(status)
658		! -----------------------------------------
659
660		! determines the presence of image
661		call ftgkyj(unit,'NAXIS', naxis, comment, status)
662		if (status > 0) call printerror(status)
663
664		! determines the presence of an extension
665		call ftgkyl(unit,'EXTEND', extend, comment, status)
666		if (status > 0) status = 0 ! no extension :
667		! to be compatible with first version of the code
668
669		if (.not.extend) then
670		print*,'No extension!'
671		call fatal_error
672		endif
673
674		! go to assigned extension (if it exists)
675		nmove = +1
676		if (present(extno)) nmove = +1 + extno
677		call ftmrhd(unit, nmove, hdutype, status)
678		if (status > 0) then
679		! if required extension not found:
680		! print a warning, fill with dummy values, return to calling routine
681		print*,code//' WARNING: the extension ',extno,' was not found in ',trim(filename)
682		alms(0:nalms-1,1)=0. ! l = 0
683		alms(0:nalms-1,2)=1. ! m = 1
684		alms(0:nalms-1,3:6)=0.
685		status = 0
686		call ftclos(unit, status) ! close the file
687		return
688		endif
689
690		if (hdutype /= 2) then ! not a binary table
691		print*, 'this is not a binary table'
692		call fatal_error
693		endif
694
695		! reads all the keywords
696		call ftghbn(unit, MAXDIM, &
697		& nrows, tfields, ttype, tform, tunit, extname, varidat, &
698		& status)
699
700		if ((tfields/=3).and.(tfields/=5)) then
701		print *,'found ',tfields,' columns in the file'
702		print *,'expected 3 or 5'
703		call fatal_error
704		endif
705		! finds the bad data value
706		! if (KMAP == SP) call ftgkye(unit,'BAD_DATA',nullval,comment,status)
707		! if (KMAP == DP) call ftgkyd(unit,'BAD_DATA',nullval,comment,status)
708		call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
709		if (status == 202) then ! bad_data not found
710		if (KMAP == SP) nullval = s_bad_value ! default value
711		if (KMAP == DP) nullval = d_bad_value ! default value
712		status = 0
713		endif
714
715		! if (nlheader > 0) then
716		! header = ""
717		! status = 0
718		! call get_clean_header(unit, header, filename, status)
719		! endif
720
721		if (present(units)) then
722		units(1) = tunit(2)
723		endif
724
725		!parse TFORM keyword to find out the length of the column vector
726		call ftbnfm(tform(1), datacode, repeat, width, status)
727		npix = nrows * repeat
728		if (npix /= nalms) then
729		print *,'found ',npix,' alms'
730		print *,'expected ',nalms
731		call fatal_error
732		endif
733
734		call ftgrsz(unit, nrow2read, status)
735		nrow2read = max(nrow2read, 1)
736		nelem = nrow2read * repeat
737		i0 = 0_i8b
738		allocate(lm(0:nelem-1))
739		do frow = 1, nrows, nrow2read
740		i1 = min(i0 + nrow2read * repeat, int(npix,i8b)) - 1_i8b
741		nelem = i1 - i0 + 1
742		! first column -> index
743		call ftgcvj(unit, 1_i4b, frow, 1_i4b, nelem, i_bad_value, &
744		& lm(0), anynull, status)
745		call assert(.not. anynull, 'There are undefined values in the table!')
746		! other columns -> a(l,m)
747		do imap = 2, tfields
748		call f90ftgcv_(unit, imap, frow, 1_i4b, nelem, nullval, &
749		& alms(i0:i1,imap+1), anynull, status)
750		call assert (.not. anynull, 'There are undefined values in the table!')
751		enddo
752		! recoding of the mapping, EH, 2002-08-08
753		! lm = l*l + l + m + 1
754		do i = i0, i1
755		l = int(sqrt( real(lm(i-i0)-1, kind = DP) ) )
756		m = lm(i-i0) - l*(l+1) - 1
757		! check that round-off did not screw up mapping
758		if (abs(m) > l) then
759		print*,'Inconsistent l^2+l+m+1 -> l,m mapping'
760		print,l, m, l(l+1)+m+1, lm(i-i0)
761		call fatal_error
762		endif
763		alms(i,1) = real( l, kind=KMAP)
764		alms(i,2) = real( m, kind=KMAP)
765		enddo
766		i0 = i1 + 1_i8b
767		enddo
768		deallocate(lm)
769		! sanity check
770		if (i0 /= npix) then
771		print*,'something wrong during piece-wise reading'
772		call fatal_error
773		endif
774
775		! !reads the columns
776		! ! first column : index -> (l,m)
777		! allocate(lm(0:nalms-1))
778		! column = 1
779		! frow = 1
780		! firstpix = 1
781		! npix = nrows * repeat
782		! if (npix /= nalms) then
783		! print *,'found ',npix,' alms'
784		! print *,'expected ',nalms
785		! call fatal_error
786		! endif
787		! call ftgcvj(unit, column, frow, firstpix, npix, i_bad_value, &
788		! & lm(0), anynull, status)
789		! call assert (.not. anynull, 'There are undefined values in the table!')
790		! ! recoding of the mapping, EH, 2002-08-08
791		! ! lm = l*l + l + m + 1
792		! do i = 0, nalms - 1
793		! l = int(sqrt( real(lm(i)-1, kind = DP) ) )
794		! m = lm(i) - l*(l+1) - 1
795		! ! check that round-off did not screw up mapping
796		! if (abs(m) > l) then
797		! print*,'Inconsistent l^2+l+m+1 -> l,m mapping'
798		! print,l, m, l(l+1)+m+1, lm(i)
799		! call fatal_error
800		! endif
801		! alms(i,1) = real( l, kind=KMAP )
802		! alms(i,2) = real( m, kind=KMAP )
803		! enddo
804		! deallocate(lm)
805		!
806		! do imap = 2, tfields
807		! !parse TFORM keyword to find out the length of the column vector
808		! call ftbnfm(tform(imap), datacode, repeat, width, status)
809		!
810		! !reads the columns
811		! column = imap
812		! frow = 1
813		! firstpix = 1
814		! npix = nrows * repeat
815		! if (npix /= nalms) then
816		! print *,'found ',npix,' alms'
817		! print *,'expected ',nalms
818		! call fatal_error
819		! endif
820		! call f90ftgcv_(unit, column, frow, firstpix, npix, nullval, &
821		! & alms(0:npix-1,imap+1), anynull, status)
822		! call assert (.not. anynull, 'There are undefined values in the table!')
823		! enddo
824		! close the file
825		call ftclos(unit, status)
826
827
828		! check for any error, and if so print out error messages
829		if (status > 0) call printerror(status)
830		return
831		end subroutine read_conbintab_KLOAD
832
833		#ifndef NO64BITS
834		subroutine write_bintab4_KLOAD(map, npix, nmap, header, nlheader, filename, extno)
835		!=======================================================================
836		INTEGER(I4B), INTENT(IN) :: npix, nmap, nlheader
837		REAL(KMAP), INTENT(IN), DIMENSION(0:npix-1,1:nmap) :: map
838		CHARACTER(LEN=*), INTENT(IN), DIMENSION(1:nlheader) :: header
839		CHARACTER(LEN=*), INTENT(IN) :: filename
840		INTEGER(I4B) , INTENT(IN), optional :: extno
841
842		integer(i8b) :: npix8
843
844		npix8 = npix
845		if (present(extno)) then
846		call write_bintab8_KLOAD(map, npix8, nmap, header, nlheader, filename, extno)
847		else
848		call write_bintab8_KLOAD(map, npix8, nmap, header, nlheader, filename)
849		endif
850		return
851		end subroutine write_bintab4_KLOAD
852		#endif
853
854		subroutine write_bintab8_KLOAD(map, npix, nmap, header, nlheader, filename, extno)
855		!=======================================================================
856		! Create a FITS file containing a binary table extension with
857		! the temperature map in the first column
858		! written by EH from writeimage and writebintable
859		! (fitsio cookbook package)
860		!
861		! slightly modified to deal with vector column (ie TFORMi = '1024E')
862		! in binary table EH/IAP/Jan-98
863		!
864		! simplified the calling sequence, the header sould be filled in
865		! before calling the routine
866		!
867		! July 21, 2004: SP version
868		! Jan 2005, EH, improved for faster writing
869		! 2009-02-25: EH, accepts I4B and I8B npix
870		!=======================================================================
871		use pix_tools, only: npix2nside
872		INTEGER(I8B), INTENT(IN) :: npix
873		INTEGER(I4B), INTENT(IN) :: nmap, nlheader
874		! REAL(KMAP), INTENT(IN), DIMENSION(0:npix-1,1:nmap), target :: map
875		REAL(KMAP), INTENT(IN), DIMENSION(0:npix-1,1:nmap) :: map
876		CHARACTER(LEN=*), INTENT(IN), DIMENSION(1:nlheader) :: header
877		CHARACTER(LEN=*), INTENT(IN) :: filename
878		INTEGER(I4B) , INTENT(IN), optional :: extno
879
880		INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1)
881		INTEGER(I4B) :: i, nside
882		LOGICAL(LGT) :: simple,extend
883		CHARACTER(LEN=80) :: comment, ch
884
885		INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
886		INTEGER(I4B) :: nrows, tfields, varidat
887		integer(i4b) :: repeat, nrow2write, nelem
888		integer(i8b) :: i0, i1
889		INTEGER(I4B) :: frow, felem, colnum, hdutype
890		CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
891		CHARACTER(LEN=10) :: card
892		CHARACTER(LEN=2) :: stn
893		INTEGER(I4B) :: itn, extno_i
894		character(len=filenamelen) sfilename
895		character(len=1) :: pform
896		! character(len=80) :: junk_k, junk_v, junk_c
897		!-----------------------------------------------------------------------
898
899		if (KMAP == SP) pform = 'E'
900		if (KMAP == DP) pform = 'D'
901		status=0
902		unit = 144
903		blocksize=1
904
905		extno_i = 0
906		if (present(extno)) extno_i = extno
907
908		if (extno_i == 0) then
909		!*************************************
910		! create the new empty FITS file
911		!*************************************
912		call ftinit(unit,filename,blocksize,status)
913		if (status > 0) call fatal_error("Error while creating file " &
914		& //trim(filename) &
915		& //". Check path and/or access rights.")
916
917		! -----------------------------------------------------
918		! initialize parameters about the FITS image
919		simple=.true.
920		bitpix=32 ! integer*4
921		naxis=0 ! no image
922		naxes(1)=0
923		extend=.true. ! there is an extension
924
925		! ----------------------
926		! primary header
927		! ----------------------
928		! write the required header keywords
929		call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
930
931		! writes supplementary keywords : none
932
933		! write the current date
934		call ftpdat(unit,status) ! format (yyyy-mm-dd)
935
936		! ----------------------
937		! image : none
938		! ----------------------
939
940		! ----------------------
941		! extension
942		! ----------------------
943		else
944
945		!*********************************************
946		! reopen an existing file and go to the end
947		!*********************************************
948		! remove the leading '!' (if any) when reopening the same file
949		sfilename = adjustl(filename)
950		if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
951		call ftopen(unit,sfilename,1_i4b,blocksize,status)
952		call ftmahd(unit,1_i4b+extno_i,hdutype,status)
953
954		endif
955
956		! creates an extension
957		call ftcrhd(unit, status)
958
959		! writes required keywords
960		tfields = nmap
961		repeat = 1024
962		if (npix < 1024) repeat = 1 ! for nside <= 8
963		! for large npix increase repeat so that nrows < 2^31-1
964		nside = npix2nside(npix)
965		if (nside > 1024*256) repeat = nside/256
966		nrows = npix / repeat ! naxis1
967		ch = string(repeat, format='(i8)')
968		ch = trim(adjustl(ch))//pform
969		tform(1:nmap) = ch
970		ttype(1:nmap) = 'simulation' ! will be updated
971		tunit(1:nmap) = '' ! optional, will not appear
972		extname = '' ! optional, will not appear
973		varidat = 0
974		call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
975		& extname, varidat, status)
976
977		! write the header literally, putting TFORM1 and TUNIT1 at the desired place
978		if (KMAP == SP) comment = 'data format of field: 4-byte REAL'
979		if (KMAP == DP) comment = 'data format of field: 8-byte REAL'
980		do i=1,nlheader
981		card = header(i)
982		if (card(1:5) == 'TTYPE') then ! if TTYPEi is explicitely given
983		stn = card(6:8)
984		read(stn,'(i3)') itn
985		! discard at their original location:
986		call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and
987		status = 0
988		call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
989		status = 0
990		if (itn <= tfields) then ! only put relevant information 2008-08-27
991		call putrec(unit,header(i), status) ! write new TTYPE1
992		status = 0
993		call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after
994		endif
995		elseif (header(i)/=' ') then
996		call putrec(unit,header(i), status)
997		endif
998		status = 0
999		enddo
1000
1001		! write the extension buffer by buffer
1002		call ftgrsz(unit, nrow2write, status)
1003		nrow2write = max(nrow2write, 1)
1004		felem = 1 ! starting position in FITS (element), 1 based
1005		i0 = 0_i8b ! starting element in array, 0 based
1006		do frow = 1, nrows, nrow2write
1007		i1 = min(i0 + nrow2write * repeat, npix) - 1_i8b
1008		nelem = i1 - i0 + 1
1009		do colnum = 1, nmap
1010		call f90ftpcl_(unit, colnum, frow, felem, nelem, map(i0:i1, colnum), status)
1011		enddo
1012		i0 = i1 + 1_i8b
1013		enddo
1014		! sanity check
1015		if (i0 /= npix) then
1016		call fatal_error("something wrong during piece wise writing")
1017		endif
1018
1019		! ----------------------
1020		! close and exit
1021		! ----------------------
1022
1023		! close the file and free the unit number
1024		call ftclos(unit, status)
1025
1026		! check for any error, and if so print out error messages
1027		if (status > 0) call printerror(status)
1028
1029		return
1030		end subroutine write_bintab8_KLOAD
1031		!=======================================================================
1032		subroutine write_asctab_KLOAD &
1033		& (clout, lmax, ncl, header, nlheader, filename, extno)
1034		!=======================================================================
1035		! Create a FITS file containing an ASCII table extension with
1036		! the measured power spectra
1037		! written by EH from writeimage and writeascii
1038		! (fitsio cookbook package)
1039		!
1040		! clout = power spectra with l in (0:lmax)
1041		! ncl = number of spectra
1042		! header = FITS header to be put on top of the file
1043		! nlheader = number of lines of the header
1044		! filename = FITS output file name
1045		!
1046		! 2015-09-10: EH, added extno
1047		!=======================================================================
1048		!
1049		INTEGER(I4B), INTENT(IN) :: lmax, ncl,nlheader
1050		REAL(KMAP), INTENT(IN) :: clout(0:lmax,1:ncl)
1051		CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header
1052		CHARACTER(LEN=*), INTENT(IN) :: filename
1053		INTEGER(I4B) , INTENT(IN), optional :: extno
1054
1055		INTEGER(I4B) :: bitpix,naxis,naxes(1)
1056		LOGICAL(LGT) :: simple,extend
1057		CHARACTER(LEN=10) :: card
1058
1059		INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP
1060		INTEGER(I4B) :: status,unit,blocksize,tfields,nrows,rowlen
1061		INTEGER(I4B) :: nspace,tbcol(MAXDIM),colnum,frow,felem,hdutype
1062		CHARACTER(LEN=16) :: ttype(MAXDIM),tform(MAXDIM),tunit(MAXDIM),extname
1063		CHARACTER(LEN=80) :: comment, card_tbcol
1064		CHARACTER(LEN=2) :: stn
1065		INTEGER(I4B) :: itn, i, extno_i
1066		character(len=filenamelen) sfilename
1067		character(len=6) :: form
1068		!=======================================================================
1069		if (KMAP == SP) form = 'E15.7'
1070		if (KMAP == DP) form = 'D24.15'
1071		status=0
1072		unit = 109
1073		blocksize=1
1074
1075		extno_i = 0
1076		if (present(extno)) extno_i = extno
1077
1078		if (extno_i == 0) then
1079		call ftinit(unit,filename,blocksize,status)
1080		if (status > 0) call fatal_error("Error while creating file " &
1081		& //trim(filename) &
1082		& //". Check path and/or access rights.")
1083
1084		! -----------------------------------------------------
1085		! initialize parameters about the FITS image
1086		simple=.true.
1087		bitpix=32 ! integer*4
1088		naxis=0 ! no image
1089		naxes(1)=0
1090		extend=.true. ! there is an extension
1091
1092		! ----------------------
1093		! primary header
1094		! ----------------------
1095		! write the required header keywords
1096		call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
1097
1098		! writes supplementary keywords : none
1099
1100		! write the current date
1101		call ftpdat(unit,status) ! (format ccyy-mm-dd)
1102
1103		! ----------------------
1104		! image : none
1105		! ----------------------
1106
1107		! ----------------------
1108		! extension
1109		! ----------------------
1110		else
1111		!*********************************************
1112		! reopen an existing file and go to the end
1113		!*********************************************
1114		! remove the leading '!' (if any) when reopening the same file
1115		sfilename = adjustl(filename)
1116		if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
1117		call ftopen(unit,sfilename,1_i4b,blocksize,status)
1118		call ftmahd(unit,1_i4b+extno_i,hdutype,status)
1119
1120		endif
1121
1122		! append a new empty extension onto the end of the primary array
1123		call ftcrhd(unit,status)
1124
1125		! define parameters for the ASCII table
1126		nrows = lmax+1
1127		tfields = ncl
1128		tform(1:ncl) = form
1129		ttype(1:ncl) = 'power spectrum' ! is updated by the value given in the header
1130		tunit(1:ncl) = '' ! optional, will not appear
1131		extname = '' ! optional, will not appear
1132
1133		! calculate the starting position of each column, and the total row length
1134		nspace=1
1135		call ftgabc(tfields,tform,nspace,rowlen,tbcol,status)
1136
1137		! write the required header parameters for the ASCII table
1138		call ftphtb(unit,rowlen,nrows,tfields,ttype,tbcol,tform,tunit, &
1139		& extname,status)
1140
1141		! write the header literally, putting TFORM1 at the desired place
1142		comment = ''
1143		card_tbcol=''
1144		do i=1,nlheader
1145		card = header(i)
1146		if (card(1:5) == 'TTYPE') then ! if TTYPE1 is explicitely given
1147		stn = card(6:8)
1148		read(stn,'(i3)') itn
1149		! discard at their original location:
1150		!!!!!!!call ftdkey(unit,card(1:6),status) ! old TTYPEi 2019-10-14
1151		call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi
1152		status = 0
1153		call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
1154		status = 0
1155		call ftgcrd(unit,'TBCOL'//stn,card_tbcol,status)
1156		status = 0
1157		call ftdkey(unit,'TBCOL'//stn,status) ! TBCOLi
1158		status = 0
1159		! and rewrite
1160		if (itn <= tfields) then
1161		call putrec(unit,card_tbcol,status) ! TBCOLi
1162		status = 0
1163		call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! TFORMi right after
1164		status = 0
1165		call putrec(unit,header(i), status) ! TTYPEi
1166		endif
1167		elseif (header(i)/=' ') then
1168		call putrec(unit,header(i), status)
1169		endif
1170		status = 0
1171		enddo
1172
1173		frow=1
1174		felem=1
1175		do colnum = 1, ncl
1176		call f90ftpcl_(unit, colnum, frow, felem, nrows, clout(0:nrows-1,colnum), status)
1177		enddo
1178
1179		! close the FITS file
1180		call ftclos(unit, status)
1181
1182		! check for any error, and if so print out error messages
1183		if (status > 0) call printerror(status)
1184
1185		return
1186		end subroutine write_asctab_KLOAD
1187		!=======================================================================
1188		! DUMP_ALMS
1189		!
1190		! Create/extend a FITS file containing a binary table extension with
1191		! the a_lm coefficients in each extension.
1192		!
1193		! Format of alm FITS file: in each extension, 3 columns:
1194		! index=l*l+l+m+1, real(a_lm), imag(a_lm)
1195		! the a_lm are obtained using the complex Y_lm
1196		! only the modes with m>=0 are stored (because the map is assumed real)
1197		!
1198		! The extensions contain, in general, T, E (= G) and B (= C) in that order
1199		! First extension has extno = 0
1200		!
1201		! Adapted from write_bintab, FKH/Apr-99
1202		! SP/DP overloaded, 2004-12, EH
1203		! reduced size of intermediate arrays (alms_out, lm)
1204		!
1205		subroutine dump_alms_KLOAD(filename, alms, nlmax, header, nlheader, extno)
1206		!=======================================================================
1207		INTEGER(I4B), INTENT(IN) :: nlmax, nlheader, extno
1208		COMPLEX(KMAPC), INTENT(IN), DIMENSION(0:,0:) :: alms
1209		CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header
1210		CHARACTER(LEN=*), INTENT(IN) :: filename
1211		! local variables
1212		INTEGER(I4B), DIMENSION(:), allocatable :: lm
1213		REAL(KMAP), DIMENSION(:,:), allocatable :: alms_out
1214		INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1)
1215		INTEGER(I4B) :: i,l,m,cnt,hdutype, nmmax
1216		LOGICAL(LGT) :: simple,extend, found_lmax, found_mmax
1217		CHARACTER(LEN=80) :: comment
1218
1219		integer(i8b) :: npix
1220		INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
1221		INTEGER(I4B) :: nrows, tfields, varidat
1222		INTEGER(I4B) :: frow, felem
1223		CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
1224		CHARACTER(LEN=10) :: card
1225		CHARACTER(LEN=2) :: stn
1226		INTEGER(I4B) :: itn
1227		character(len=filenamelen) sfilename
1228		character(len=1) :: pform
1229		character(len=*), parameter :: code = 'dump_alms'
1230		!-----------------------------------------------------------------------
1231
1232		if (KMAP == SP) pform = 'E'
1233		if (KMAP == DP) pform = 'D'
1234
1235		nmmax = size(alms,2) - 1
1236		if (nmmax < 0 .or. nmmax > nlmax) call fatal_error('inconsistent Lmax and Mmax in dump_alms')
1237		!! npix=((nlmax+1)*(nlmax+2))/2
1238		npix = ((nmmax+1_i8b)(2_i8bnlmax-nmmax+2))/2
1239		if (npix > MAX_I4B .or. nlmax >= int(sqrt(real(MAX_I4B,kind=dp)))) then
1240		print*,code//'> Index of a_lm too large for current file format,'
1241		print*,code//'> '//trim(filename)//' not written.'
1242		call fatal_error
1243		endif
1244
1245		status=0
1246		unit = 141
1247		blocksize=1
1248
1249		if (extno==0) then
1250		!*********************************************
1251		! create the new empty FITS file
1252		!*********************************************
1253		call ftinit(unit,filename,blocksize,status)
1254		if (status > 0) call fatal_error("Error while creating file " &
1255		& //trim(filename) &
1256		& //". Check path and/or access rights.")
1257
1258		! -----------------------------------------------------
1259		! initialize parameters about the FITS image
1260		simple=.true.
1261		bitpix=32 ! integer*4
1262		naxis=0 ! no image
1263		naxes(1)=0
1264		extend=.true. ! there is an extension
1265
1266		! ----------------------
1267		! primary header
1268		! ----------------------
1269		! write the required header keywords
1270		call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
1271
1272		! writes supplementary keywords : none
1273
1274		! write the current date
1275		call ftpdat(unit,status) ! format (yyyy-mm-dd)
1276
1277		! ----------------------
1278		! image : none
1279		! extension
1280		! ----------------------
1281		else
1282
1283		!*********************************************
1284		! reopen an existing file and go to the end
1285		!*********************************************
1286		! remove the leading '!' (if any) when reopening the same file
1287		sfilename = adjustl(filename)
1288		if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
1289		call ftopen(unit,sfilename,1_i4b,blocksize,status)
1290		call ftmahd(unit,1_i4b+extno,hdutype,status)
1291
1292		endif
1293
1294		! creates an extension
1295		call ftcrhd(unit, status)
1296
1297		! writes required keywords
1298		nrows = npix ! naxis1
1299		tfields = 3
1300		tform = '' ; ttype = '' ; tunit = ''
1301		tform(1)(1:2)= '1J' ! necessary with ifort 11.1.076
1302		tform(2)(1:2) = '1'//pform
1303		tform(3)(1:2) = '1'//pform
1304		ttype(1)(1:15) = 'index=l^2+l+m+1'
1305		ttype(2)(1:15) = 'alm (real) '
1306		ttype(3)(1:15) = 'alm (imaginary)'
1307		tunit(1:3) = '' ! optional, will not appear
1308		extname = '' ! optional, will not appear
1309		varidat = 0
1310		call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
1311		& extname, varidat, status)
1312
1313		! write the header literally, putting TFORM1 at the desired place
1314		do i=1,nlheader
1315		card = header(i)
1316		if (card(1:5) == 'TTYPE') then ! if TTYPEi is explicitely given
1317		stn = card(6:8)
1318		read(stn,'(i3)') itn
1319		! discard at their original location:
1320		call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and ! remove
1321		status = 0
1322		call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
1323		status = 0
1324		if (itn <= tfields) then ! only put relevant information 2008-08-27
1325		call putrec(unit,header(i), status) ! write new TTYPEi
1326		status = 0
1327		if (itn==1) then
1328		comment = 'data format of field: 4-byte INTEGER'
1329		else
1330		if (KMAP == SP) comment = 'data format of field: 4-byte REAL'
1331		if (KMAP == DP) comment = 'data format of field: 8-byte REAL'
1332		endif
1333		call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after
1334		endif
1335		elseif (header(i)/=' ') then
1336		call putrec(unit,header(i), status)
1337		endif
1338		status = 0
1339		enddo
1340
1341		call ftukyj(unit, 'MIN-LPOL', 0_i4b, 'Minimum L multipole order', status)
1342		call ftukyj(unit, 'MAX-LPOL', nlmax, 'Maximum L multipole order', status)
1343		call ftukyj(unit, 'MAX-MPOL', nmmax, 'Maximum M multipole degree', status)
1344
1345		allocate(lm (0:nmmax))
1346		allocate(alms_out(0:nmmax,1:2))
1347		! write the extension one column by one column
1348		frow = 1 ! starting position (row)
1349		felem = 1 ! starting position (element)
1350		do l = 0, nlmax
1351		cnt = 0
1352		do m = 0, min(l,nmmax)
1353		lm(cnt) = l**2 + l + m + 1
1354		alms_out(cnt,1)=REAL( alms(l,m))
1355		alms_out(cnt,2)=AIMAG(alms(l,m))
1356		cnt = cnt + 1
1357		enddo
1358		call ftpclj( unit, 1_i4b, frow, felem, cnt, lm, status)
1359		call f90ftpcl_(unit, 2_i4b, frow, felem, cnt, alms_out(0:cnt-1,1), status)
1360		call f90ftpcl_(unit, 3_i4b, frow, felem, cnt, alms_out(0:cnt-1,2), status)
1361		frow = frow + cnt
1362		enddo
1363		deallocate(lm)
1364		deallocate(alms_out)
1365
1366		! ----------------------
1367		! close and exit
1368		! ----------------------
1369
1370		! close the file and free the unit number
1371		call ftclos(unit, status)
1372
1373		! check for any error, and if so print out error messages
1374		if (status > 0) call printerror(status)
1375
1376
1377		return
1378		end subroutine dump_alms_KLOAD
1379		!=======================================================================
1380		subroutine write_alms_KLOAD(filename, nalms, alms, ncl, header, nlheader, extno)
1381		!=======================================================================
1382		! Writes alms from to binary FITS file, FKH/Apr-99
1383		! ncl is the number of columns, in the output fits file,
1384		! either 3 or 5 (with or without errors respectively)
1385		!
1386		! input array (real) FITS file
1387		! alms(:,1) = l )
1388		! alms(:,2) = m )---> col 1: l*l+l+m+1
1389		! alms(:,3) = real(a_lm) ---> col 2
1390		! alms(:,4) = imag(a_lm) ---> col 3
1391		! alms(:,5) = real(delta a_lm) ---> col 4
1392		! alms(:,6) = imag(delta a_lm) ---> col 5
1393		!
1394		!=======================================================================
1395
1396		INTEGER(I4B), INTENT(IN) :: nalms, nlheader, ncl, extno
1397		REAL(KMAP), INTENT(IN), DIMENSION(0:nalms-1,1:(ncl+1)) :: alms
1398		CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header
1399		CHARACTER(LEN=*), INTENT(IN) :: filename
1400
1401		INTEGER(I4B), DIMENSION(:), allocatable :: lm
1402		INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1)
1403		INTEGER(I4B) :: i,hdutype, lmax, mmax, lmin
1404		LOGICAL(LGT) :: simple,extend
1405		CHARACTER(LEN=80) :: comment
1406
1407		INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
1408		INTEGER(I4B) :: nrows, npix, tfields, varidat, repeat
1409		INTEGER(I4B) :: frow, felem, colnum, stride, istart, iend, k
1410		CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
1411		CHARACTER(LEN=10) :: card
1412		CHARACTER(LEN=2) :: stn
1413		INTEGER(I4B) :: itn
1414		integer(I4B) :: l, m
1415		character(len=filenamelen) sfilename
1416		character(len=1) :: pform
1417		!-----------------------------------------------------------------------
1418
1419		if (KMAP == SP) pform = 'E'
1420		if (KMAP == DP) pform = 'D'
1421
1422		status=0
1423		unit = 140
1424
1425		! create the new empty FITS file
1426		blocksize=1
1427
1428		if (extno==1) then
1429
1430		call ftinit(unit,filename,blocksize,status)
1431		if (status > 0) call fatal_error("Error while creating file " &
1432		& //trim(filename) &
1433		& //". Check path and/or access rights.")
1434
1435		! -----------------------------------------------------
1436		! initialize parameters about the FITS image
1437		simple=.true.
1438		bitpix=32 ! integer*4
1439		naxis=0 ! no image
1440		naxes(1)=0
1441		extend=.true. ! there is an extension
1442
1443		! ----------------------
1444		! primary header
1445		! ----------------------
1446		! write the required header keywords
1447		call ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
1448
1449		! writes supplementary keywords : none
1450
1451		! write the current date
1452		call ftpdat(unit,status) ! format ccyy-mm-dd
1453
1454		! ----------------------
1455		! image : none
1456		! ----------------------
1457
1458		! ----------------------
1459		! extension
1460		! ----------------------
1461
1462		else
1463
1464		!*********************************************
1465		! reopen an existing file and go to the end
1466		!*********************************************
1467		! remove the leading '!' (if any) when reopening the same file
1468		sfilename = adjustl(filename)
1469		if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
1470		call ftopen(unit,sfilename,1_i4b,blocksize,status)
1471		call ftmahd(unit,extno,hdutype,status)
1472
1473		endif
1474
1475		! creates an extension
1476		call ftcrhd(unit, status)
1477
1478		! writes required keywords
1479		nrows = nalms ! naxis1
1480		tfields = ncl
1481		repeat = 1
1482		tform(1)(1:2)='1J' ! necessary with ifort 11.1.076
1483		do i=2, ncl
1484		tform(i)(1:2) = '1'//pform
1485		enddo
1486		ttype(1)(1:15) = 'index=l^2+l+m+1'
1487		ttype(2)(1:15) = 'alm (real) '
1488		ttype(3)(1:15) = 'alm (imaginary)'
1489		if (ncl>3) then
1490		ttype(4)(1:17) = 'error (real) '
1491		ttype(5)(1:17) = 'error (imaginary)'
1492		endif
1493		tunit(1:ncl) = '' ! optional, will not appear
1494		extname = '' ! optional, will not appear
1495		varidat = 0
1496		call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
1497		& extname, varidat, status)
1498
1499		! write the header literally, putting TFORMi at the desired place
1500		do i=1,nlheader
1501		card = header(i)
1502		if (card(1:5) == 'TTYPE') then ! if TTYPEi is explicitely given
1503		stn = card(6:8)
1504		read(stn,'(i3)') itn
1505		! discard at their original location:
1506		call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and ! remove
1507		status = 0
1508		call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
1509		status = 0
1510		if (itn <= tfields) then ! only put relevant information 2008-08-27
1511		call putrec(unit,header(i), status) ! write new TTYPE1
1512		status = 0
1513		if (itn==1) then
1514		comment = 'data format of field: 4-byte INTEGER'
1515		else
1516		if (KMAP == SP) comment = 'data format of field: 4-byte REAL'
1517		if (KMAP == DP) comment = 'data format of field: 8-byte REAL'
1518		endif
1519		call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after
1520		endif
1521		elseif (header(i)/=' ') then
1522		call putrec(unit,header(i), status)
1523		endif
1524		status = 0
1525		enddo
1526
1527		lmax = nint(maxval(alms(:,1)))
1528		lmin = nint(minval(alms(:,1)))
1529		mmax = nint(maxval(alms(:,2)))
1530		call ftukyj(unit, 'MIN-LPOL', lmin, 'Minimum L multipole order', status)
1531		call ftukyj(unit, 'MAX-LPOL', lmax, 'Maximum L multipole order', status)
1532		call ftukyj(unit, 'MAX-MPOL', mmax, 'Maximum M multipole degree', status)
1533
1534		! write the extension by blocks of rows ! EH, Dec 2004
1535		felem = 1 ! starting position (element)
1536		!!! stride = 1000 ! 1000 rows at a time
1537		call ftgrsz(unit, stride, status) ! find optimal stride in rows
1538		stride = max( stride, 1)
1539		allocate(lm(0:stride-1))
1540		do k = 0, (nalms-1)/(stride * repeat)
1541		istart = k * (stride * repeat)
1542		iend = min(nalms, istart + stride * repeat) - 1
1543		do i = istart, iend ! recode the (l,m) -> lm mapping, EH, 2002-08-08
1544		l = nint(alms(i,1))
1545		m = nint(alms(i,2))
1546		lm(i-istart) = (l+1)*l + m + 1
1547		enddo
1548
1549		frow = istart/repeat + 1
1550		npix = iend - istart + 1
1551		call ftpclj( unit, 1_i4b, frow, felem, npix, lm, status)
1552		do colnum = 2, ncl
1553		call f90ftpcl_(unit, colnum, frow, felem, npix, alms(istart:iend,colnum+1), status)
1554		enddo
1555		enddo
1556		deallocate(lm)
1557
1558		! ----------------------
1559		! close and exit
1560		! ----------------------
1561
1562		! close the file and free the unit number
1563		call ftclos(unit, status)
1564
1565		! check for any error, and if so print out error messages
1566		if (status > 0) call printerror(status)
1567
1568		return
1569		end subroutine write_alms_KLOAD
1570		!=======================================================================
1571		subroutine read_alms_KLOAD(filename, nalms, alms, ncl, header, nlheader, extno)
1572		!=======================================================================
1573		! Read a FITS file
1574		!
1575		! slightly modified to deal with vector column
1576		! in binary table EH/IAP/Jan-98
1577		!
1578		! Used by synfast when reading a binary file with alms for cons.real.
1579		! FKH/Apr-99
1580		!
1581		! called by fits2alms
1582		! Jan 2005, EH, improved for faster reading
1583		!=======================================================================
1584		CHARACTER(LEN=*), INTENT(IN) :: filename
1585		INTEGER(I4B), INTENT(IN) :: nalms, ncl,nlheader,extno
1586		REAL(KMAP), DIMENSION(0:nalms-1,1:(ncl+1)), INTENT(OUT) :: alms
1587		CHARACTER(LEN=80), INTENT(OUT), DIMENSION(1:nlheader) :: header
1588		REAL(KMAP) :: nullval
1589		LOGICAL(LGT) :: anynull
1590
1591		INTEGER(I4B), DIMENSION(:), allocatable :: lm
1592		INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis
1593		INTEGER(I4B) :: npix
1594		CHARACTER(LEN=80) :: comment ! , record
1595		LOGICAL(LGT) :: extend
1596		INTEGER(I4B) :: nmove, hdutype ! , nkeys , nspace
1597		INTEGER(I4B) :: frow, imap
1598		INTEGER(I4B) :: datacode, repeat, width
1599		integer(I4B) :: i, l, m
1600		integer(i4b) :: nrow2read, nelem
1601		integer(i8b) :: i0, i1
1602
1603		INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
1604		INTEGER(I4B) :: nrows, tfields, varidat
1605		CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
1606
1607		!-----------------------------------------------------------------------
1608		status=0
1609		header=''
1610		comment=''
1611		ttype=''
1612		tform=''
1613		tunit=''
1614		extname=''
1615		unit = 139
1616		naxes(1) = 1
1617		naxes(2) = 1
1618		nfound = -1
1619		anynull = .false.
1620		alms=0.
1621		readwrite=0
1622		call ftopen(unit,filename,readwrite,blocksize,status)
1623		if (status > 0) then
1624		call printerror(status)
1625		call fatal_error("Aborting.")
1626		endif
1627		! -----------------------------------------
1628
1629		! determines the presence of image
1630		call ftgkyj(unit,'NAXIS', naxis, comment, status)
1631		if (status > 0) call printerror(status)
1632
1633		! determines the presence of an extension
1634		call ftgkyl(unit,'EXTEND', extend, comment, status)
1635		if (status > 0) status = 0 ! no extension :
1636		! to be compatible with first version of the code
1637
1638		call assert (extend, 'No extension!')
1639		nmove = +extno
1640		call ftmrhd(unit, nmove, hdutype, status)
1641		!cc write(,) hdutype
1642
1643		call assert(hdutype==2, 'this is not a binary table')
1644
1645		header = ""
1646		call get_clean_header( unit, header, filename, status)
1647
1648		! reads all the keywords
1649		call ftghbn(unit, MAXDIM, &
1650		& nrows, tfields, ttype, tform, tunit, extname, varidat, &
1651		& status)
1652
1653		if (tfields<ncl) then
1654		print *,'found ',tfields,' columns in the file'
1655		print *,'expected ',ncl
1656		call fatal_error
1657		endif
1658		! finds the bad data value
1659		! if (KMAP == SP) call ftgkye(unit,'BAD_DATA',nullval,comment,status)
1660		! if (KMAP == DP) call ftgkyd(unit,'BAD_DATA',nullval,comment,status)
1661		call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
1662		if (status == 202) then ! bad_data not found
1663		if (KMAP == SP) nullval = s_bad_value ! default value
1664		if (KMAP == DP) nullval = d_bad_value ! default value
1665		status = 0
1666		endif
1667		!parse TFORM keyword to find out the length of the column vector
1668		call ftbnfm(tform(1), datacode, repeat, width, status)
1669		npix = nrows * repeat
1670		! if (npix /= nalms) then
1671		if (npix > nalms) then
1672		print *,'found ',npix,' alms'
1673		! print *,'expected ',nalms
1674		print *,'expected ',nalms,' or less'
1675		call fatal_error
1676		endif
1677
1678		call ftgrsz(unit, nrow2read, status)
1679		nrow2read = max(nrow2read, 1)
1680		nelem = nrow2read * repeat
1681		i0 = 0_i8b
1682		allocate(lm(0:nelem-1))
1683		do frow = 1, nrows, nrow2read
1684		i1 = min(i0 + nrow2read * repeat, int(npix,i8b)) - 1_i8b
1685		nelem = i1 - i0 + 1
1686		! first column -> index
1687		call ftgcvj(unit, 1_i4b, frow, 1_i4b, nelem, i_bad_value, &
1688		& lm(0), anynull, status)
1689		call assert(.not. anynull, 'There are undefined values in the table!')
1690		! other columns -> a(l,m)
1691		do imap = 2, ncl
1692		call f90ftgcv_(unit, imap, frow, 1_i4b, nelem, nullval, &
1693		& alms(i0:i1,imap+1), anynull, status)
1694		call assert (.not. anynull, 'There are undefined values in the table!')
1695		enddo
1696		! recoding of the mapping, EH, 2002-08-08
1697		! lm = l*l + l + m + 1
1698		do i = i0, i1
1699		l = int(sqrt( real(lm(i-i0)-1, kind = DP) ) )
1700		m = lm(i-i0) - l*(l+1) - 1
1701		! check that round-off did not screw up mapping
1702		if (abs(m) > l) then
1703		print*,'Inconsistent l^2+l+m+1 -> l,m mapping'
1704		print,l, m, l(l+1)+m+1, lm(i-i0)
1705		call fatal_error
1706		endif
1707		alms(i,1) = real( l, kind=KMAP)
1708		alms(i,2) = real( m, kind=KMAP)
1709		enddo
1710		i0 = i1 + 1_i8b
1711		enddo
1712		deallocate(lm)
1713		! sanity check
1714		if (i0 /= npix) then
1715		print*,'something wrong during piece-wise reading'
1716		call fatal_error
1717		endif
1718
1719		! close the file
1720		call ftclos(unit, status)
1721
1722
1723		! check for any error, and if so print out error messages
1724		if (status > 0) call printerror(status)
1725		return
1726		end subroutine read_alms_KLOAD
1727		!**************************************************************************
1728		SUBROUTINE read_bintod_KLOAD(filename, tod, npixtot, ntods, firstpix, nullval, anynull, &
1729		header, extno, units)
1730		!**************************************************************************
1731		!=======================================================================
1732		! Read a FITS file
1733		!
1734		! slightly modified to deal with vector column (ie TFORMi = '1024E')
1735		! in binary table EH/IAP/Jan-98
1736		!
1737		! This routine is used for reading TODS by anafast.
1738		! Modified to start at a given pix numb OD & RT 02/02
1739		! Modified to handle huge array (npix_tot > 2^32) OD & EH 07/02
1740		! 2002-07-08 : bugs correction by E.H.
1741		! 2007-01-31 : added units optional output
1742		!=======================================================================
1743
1744		IMPLICIT NONE
1745
1746		CHARACTER(LEN=*), INTENT(IN) :: filename
1747		INTEGER(I8B) , INTENT(IN) :: npixtot,firstpix
1748		INTEGER(I4B), INTENT(IN) :: ntods
1749		REAL(KMAP), DIMENSION(0:,1:), INTENT(OUT) :: tod
1750		REAL(KMAP), INTENT(OUT) :: nullval
1751		LOGICAL(LGT), INTENT(OUT) :: anynull
1752		character(len=*), dimension(1:),intent(out), optional :: header
1753		INTEGER(I4B), INTENT(IN), OPTIONAL :: extno
1754		character(LEN=*), dimension(1:),intent(OUT), optional :: units
1755
1756
1757		integer(I4B) :: i, nl_header, len_header, nl_units, len_units
1758		INTEGER(I4B) :: status,unit,readwrite,blocksize,naxes(2),nfound, naxis
1759		INTEGER(I4B) :: npix_32 !,firstpix_32
1760		CHARACTER(LEN=80) :: comment
1761		LOGICAL(LGT) :: extend
1762		INTEGER(I4B) :: nmove, hdutype
1763		INTEGER(I4B) :: column, frow, itod
1764		INTEGER(I4B) :: datacode, repeat, width
1765
1766		INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
1767		INTEGER(I4B) :: nrows, tfields, varidat,felem
1768		CHARACTER(LEN=20), dimension(1:MAXDIM) :: ttype, tform,tunit
1769		CHARACTER(LEN=20) :: extname
1770
1771		INTEGER(I8B) :: q,iq,npix_tmp,firstpix_tmp, i0, i1
1772
1773		!-----------------------------------------------------------------------
1774		status=0
1775
1776		unit = 138
1777		naxes(1) = 1
1778		naxes(2) = 1
1779		nfound = -1
1780		anynull = .FALSE.
1781
1782		nl_header = 0
1783		if (present(header)) then
1784		nl_header = size(header)
1785		len_header = 80
1786		endif
1787
1788		nl_units = 0
1789		if (present(units)) then
1790		nl_units = size(units)
1791		len_units = min(80,len(units(1))) ! due to SUN compiler bug
1792		endif
1793
1794		readwrite=0
1795		CALL ftopen(unit,filename,readwrite,blocksize,status)
1796		IF (status .GT. 0) CALL printerror(status)
1797		! -----------------------------------------
1798
1799		! determines the presence of image
1800		CALL ftgkyj(unit,'NAXIS', naxis, comment, status)
1801		IF (status .GT. 0) CALL printerror(status)
1802
1803		! determines the presence of an extension
1804		CALL ftgkyl(unit,'EXTEND', extend, comment, status)
1805		IF (status .GT. 0) status = 0 ! no extension :
1806		! to be compatible with first version of the code
1807
1808		IF (naxis .GT. 0) THEN ! there is an image
1809		print*,'WARNING : Image is ignored in '//trim(filename)
1810		ENDIF
1811		IF (extend) THEN ! there is an extension
1812
1813		nmove = +1
1814		if (present(extno)) nmove = +1 + extno
1815		CALL ftmrhd(unit, nmove, hdutype, status)
1816		!cc write(,) hdutype
1817
1818		call assert(hdutype==2, 'this is not a binary table')
1819
1820		! reads all the keywords
1821		CALL ftghbn(unit, MAXDIM, &
1822		& nrows, tfields, ttype, tform, tunit, extname, varidat, &
1823		& status)
1824
1825		IF (tfields .LT. ntods) THEN
1826		PRINT *,'found ',tfields,' tods in the file'
1827		PRINT *,'expected ',ntods
1828		call fatal_error
1829		ENDIF
1830
1831		if (nl_header > 0) then
1832		header = ""
1833		status = 0
1834		call get_clean_header(unit, header, filename, status)
1835		endif
1836
1837		if (nl_units > 0) then
1838		do i=1,nl_units
1839		units(i)(1:len_units) = 'unknown' ! default
1840		enddo
1841		do itod = 1, min(ntods, nl_units)
1842		units(itod)(1:len_units) = adjustl(tunit(itod))
1843		enddo
1844		endif
1845
1846		! finds the bad data value
1847		! if (KMAP == SP) CALL ftgkye(unit,'BAD_DATA',nullval,comment,status)
1848		! if (KMAP == DP) CALL ftgkyd(unit,'BAD_DATA',nullval,comment,status)
1849		call f90ftgky_(unit, 'BAD_DATA', nullval, comment, status)
1850		IF (status .EQ. 202) THEN ! bad_data not found
1851		if (KMAP == SP) nullval = s_bad_value ! default value
1852		if (KMAP == DP) nullval = d_bad_value ! default value
1853		status = 0
1854		ENDIF
1855
1856		IF (npixtot .LT. nchunk_max) THEN
1857
1858		DO itod = 1, ntods
1859
1860		!parse TFORM keyword to find out the length of the column vector (repeat)
1861		CALL ftbnfm(tform(itod), datacode, repeat, width, status)
1862		frow = (firstpix)/repeat+1 ! 1 based
1863		felem = firstpix-(frow-1)*repeat+1 ! 1 based
1864		npix_32 = npixtot
1865
1866		!reads the columns
1867		column = itod
1868		CALL f90ftgcv_(unit, column, frow, felem, npix_32, nullval, &
1869		& tod(0:npix_32-1,itod), anynull, status)
1870		END DO
1871
1872		ELSE
1873
1874		q = (npixtot-1)/nchunk_max
1875		DO iq = 0,q
1876		IF (iq .LT. q) THEN
1877		npix_tmp = nchunk_max
1878		ELSE
1879		npix_tmp = npixtot - iq*nchunk_max
1880		ENDIF
1881		firstpix_tmp = firstpix + iq*nchunk_max
1882		npix_32 = npix_tmp
1883		i0 = firstpix_tmp-firstpix
1884		i1 = i0 + npix_tmp - 1_i8b
1885
1886		DO itod = 1, ntods
1887		! parse TFORM keyword to find out the length of the column vector
1888		CALL ftbnfm(tform(itod), datacode, repeat, width, status)
1889		frow = (firstpix_tmp)/repeat+1 ! 1 based
1890		felem = firstpix_tmp-(frow-1)*repeat+1 ! 1 based
1891		CALL f90ftgcv_(unit, itod, frow, felem, npix_32, nullval, &
1892		& tod(i0:i1,itod), anynull, status)
1893		END DO
1894
1895		ENDDO
1896
1897		ENDIF
1898
1899		ELSE ! no image no extension, you are dead, man
1900		call fatal_error(' No image, no extension')
1901		ENDIF
1902
1903		! close the file
1904		CALL ftclos(unit, status)
1905
1906		! check for any error, and if so print out error messages
1907		IF (status .GT. 0) CALL printerror(status)
1908
1909		RETURN
1910
1911		END SUBROUTINE read_bintod_KLOAD
1912		!=======================================================================
1913
1914		!======================================================================================
1915		SUBROUTINE write_bintabh_KLOAD(tod, npix, ntod, header, nlheader, filename, extno, firstpix, repeat)
1916		!======================================================================================
1917
1918		! =================================================================================
1919		! Create a FITS file containing a binary table extension in the first extension
1920		!
1921		! Designed to deal with Huge file, (n_elements > 2^31)
1922		!
1923		! OPTIONAL NEW PARAMETERS:
1924		! firstpix : position in the file of the first element to be written (starts at 0)
1925		! default value =0
1926		! 8-byte integer
1927		! if NE 0 then suppose that the file already exists
1928		!
1929		! repeat : length of vector per unit rows and columns of the first binary extension
1930		! default value = 12000 (\equiv 1 mn of PLANCK/HFI data)
1931		! 4-byte integer
1932		!
1933		! OTHER PARAMETERS
1934		! unchanged with respect to the standard write_bintab of the HEALPIX package except
1935		! npix which is an 8-byte integer
1936		!
1937		! Adapted from write_bintab
1938		! E.H. & O.D. @ IAP 07/02
1939		!
1940		! Requires a compilation in 64 bits of the CFITSIO
1941		! Note that the flag -D_FILE_OFFSETS_BITS=64 has to be added
1942		! (cf page CFITIO 2.2 User's guide Chap 4, section 4-13)
1943		!
1944		! 2002-07-08 : bugs correction by E.H.
1945		! (uniform use of firstpix_tmp, introduction of firstpix_chunk)
1946		! 2015-04-28: better handling of repeat (suggested by R. Keskitalo)
1947		!==========================================================================================
1948
1949		USE healpix_types
1950		IMPLICIT NONE
1951
1952		INTEGER(I8B) , INTENT(IN) :: npix
1953		INTEGER(I8B) , INTENT(IN), OPTIONAL :: firstpix
1954		INTEGER(I4B) , INTENT(IN), OPTIONAL :: repeat
1955		INTEGER(I4B) , INTENT(IN) :: ntod,nlheader
1956		REAL(KMAP) , INTENT(IN), DIMENSION(0:npix-1,1:ntod) :: tod
1957		CHARACTER(LEN=80), INTENT(IN), DIMENSION(1:nlheader) :: header
1958		CHARACTER(LEN=*), INTENT(IN) :: filename
1959		INTEGER(I4B) , INTENT(IN) , OPTIONAL :: extno
1960
1961		INTEGER(I4B) :: status,unit,blocksize,bitpix,naxis,naxes(1),repeat_fits
1962		INTEGER(I4B) :: i,npix_32
1963		LOGICAL(LGT) :: simple,extend
1964		CHARACTER(LEN=80) :: comment, ch
1965		integer(I8B) :: repeat_tmp
1966
1967		INTEGER(I4B), PARAMETER :: MAXDIM = MAXDIM_TOP !number of columns in the extension
1968		INTEGER(I4B) :: nrows,tfields,varidat
1969		INTEGER(I4B) :: frow,felem,colnum,readwrite,width,datacode,hdutype
1970		CHARACTER(LEN=20) :: ttype(MAXDIM), tform(MAXDIM), tunit(MAXDIM), extname
1971		CHARACTER(LEN=10) :: card
1972		CHARACTER(LEN=2) :: stn
1973		INTEGER(I4B) :: itn
1974
1975		real(KMAP), dimension(:), allocatable :: padding
1976		integer(I8B) :: lastpix
1977		integer(I4B) :: lengap
1978
1979		INTEGER(I4B) :: extno_i
1980		character(len=filenamelen) :: sfilename
1981		INTEGER(I8B) :: q,iq,npix_tmp,firstpix_tmp,firstpix_chunk, i0, i1
1982		character(len=1) :: pform
1983		character(len=*), parameter :: code='write_bintabh'
1984		!-----------------------------------------------------------------------
1985
1986		if (KMAP == SP) pform = 'E'
1987		if (KMAP == DP) pform = 'D'
1988
1989		IF (.NOT. PRESENT(repeat) ) THEN
1990		repeat_tmp = 1
1991		if (mod(npix,1024_i8b) == 0) then
1992		repeat_tmp = 1024
1993		elseif (npix >= 12000) then
1994		repeat_tmp = 12000
1995		endif
1996		ELSE
1997		repeat_tmp = repeat
1998		ENDIF
1999		IF (.NOT. PRESENT(firstpix) ) THEN
2000		firstpix_tmp = 0
2001		ELSE
2002		firstpix_tmp = firstpix
2003		ENDIF
2004
2005		extno_i = 0
2006		if (present(extno)) extno_i = extno
2007
2008		status=0
2009		unit = 137
2010		blocksize=1
2011
2012		! remove the leading '!' (if any) when reopening the same file
2013		sfilename = adjustl(filename)
2014		if (sfilename(1:1) == '!') sfilename = sfilename(2:filenamelen)
2015
2016		! create the new empty FITS file
2017
2018		IF (firstpix_tmp .EQ. 0) THEN
2019
2020		if (extno_i == 0) then
2021		CALL ftinit(unit,filename,blocksize,status)
2022		if (status > 0) call fatal_error("Error while creating file " &
2023		& //trim(filename) &
2024		& //". Check path and/or access rights.")
2025
2026		! -----------------------------------------------------
2027		! Initialize parameters about the FITS image
2028		simple=.TRUE.
2029		bitpix=32 ! integer*4
2030		naxis=0 ! no image
2031		naxes(1)=0
2032		extend=.TRUE. ! there is an extension
2033
2034		! ----------------------
2035		! primary header
2036		! ----------------------
2037		! write the required header keywords
2038		CALL ftphpr(unit,simple,bitpix,naxis,naxes,0_i4b,1_i4b,extend,status)
2039
2040		! writes supplementary keywords : none
2041
2042		! write the current date
2043		CALL ftpdat(unit,status) ! format ccyy-mm-dd
2044
2045		! ----------------------
2046		! image : none
2047		! ----------------------
2048
2049		! ----------------------
2050		! extension
2051		! ----------------------
2052		else
2053
2054		!*********************************************
2055		! reopen an existing file and go to the end
2056		!*********************************************
2057		call ftopen(unit,sfilename,1_i4b,blocksize,status)
2058		call ftmahd(unit,1_i4b+extno_i,hdutype,status)
2059
2060		endif
2061
2062		! creates an extension
2063		CALL ftcrhd(unit, status)
2064
2065		! writes required keywords
2066		! if (npix < repeat_tmp) repeat_tmp = npix ! 2015-04-28
2067
2068		! nrows = npix / repeat_tmp + 1 ! naxis1
2069		nrows = (npix - 1_i8b) / repeat_tmp + 1_i4b ! 2015-04-28
2070		tfields = ntod
2071		WRITE(ch,'(i8)') repeat_tmp
2072		! tform(1:ntod) = TRIM(ADJUSTL(ch))//pform ! does not work with Ifort, EH, 2006-04-04
2073		ch = TRIM(ADJUSTL(ch))//pform
2074		tform(1:ntod) = ch
2075
2076		! IF (npix .LT. repeat_tmp) THEN ! 2015-04-28
2077		! nrows = npix
2078		! ch = '1'//pform
2079		! tform(1:ntod) = ch
2080		! ENDIF
2081		ttype(1:ntod) = 'simulation' ! will be updated
2082		tunit(1:ntod) = '' ! optional, will not appear
2083
2084		extname = '' ! optional, will not appear
2085		varidat = 0
2086
2087		CALL ftphbn(unit, nrows, tfields, ttype, tform, tunit, &
2088		& extname, varidat, status)
2089
2090		! write the header literally, putting TFORM1 at the desired place
2091		if (KMAP == SP) comment = 'data format of field: 4-byte REAL'
2092		if (KMAP == DP) comment = 'data format of field: 8-byte REAL'
2093		DO i=1,nlheader
2094		card = header(i)
2095		IF (card(1:5) == 'TTYPE') THEN ! if TTYPE1 is explicitely given
2096		stn = card(6:8)
2097		READ(stn,'(i3)') itn
2098		! discard at their original location:
2099		call ftdkey(unit,'TTYPE'//stn,status) ! old TTYPEi and
2100		status = 0
2101		call ftdkey(unit,'TFORM'//stn,status) ! TFORMi
2102		status = 0
2103		if (itn <= tfields) then ! only put relevant information 2008-08-27
2104		call putrec(unit,header(i), status) ! write new TTYPE1
2105		status = 0
2106		call ftpkys(unit,'TFORM'//stn,tform(itn),comment,status) ! and write new TFORM1 right after
2107		!CALL ftmcrd(unit,'TTYPE'//stn,'COMMENT',status) ! old TTYPEi and
2108		!CALL ftmcrd(unit,'TFORM'//stn,'COMMENT',status) ! TFORMi
2109		!CALL ftprec(unit,header(i), status) ! write new TTYPE1
2110		!CALL ftpkys(unit,'TFORM'//stn,tform(1),comment,status) ! and write new TFORM1 right after
2111		endif
2112		ELSEIF (header(i).NE.' ') THEN
2113		call putrec(unit,header(i), status)
2114		! CALL ftprec(unit,header(i), status)
2115		ENDIF
2116		status = 0
2117		ENDDO
2118
2119		ELSE
2120		! The file already exists
2121		readwrite=1
2122		CALL ftopen(unit,sfilename,readwrite,blocksize,status)
2123		CALL ftmahd(unit,2_i4b+extno_i,hdutype,status)
2124
2125		CALL ftgkys(unit,'TFORM1',tform(1),comment,status)
2126		CALL ftbnfm(tform(1),datacode,repeat_fits,width,status)
2127
2128		IF (repeat_tmp .NE. repeat_fits) THEN
2129		if (present(repeat)) then
2130		write(*,'(a)') code//'> WARNING: In file '//trim(sfilename)
2131		write(*,'(a)') &
2132		& code//'> user provided REPEAT value (' // &
2133		& trim(adjustl(string(repeat_tmp))) // &
2134		& ') differs from value read from file (' // &
2135		& trim(adjustl(string(repeat_fits))) // &
2136		& ').'
2137		write(*,'(a)') code//'> The latter will be used.'
2138		! else
2139		! write(*,'(a,i0.0,a)') &
2140		! & code//'> WARNING: REPEAT value read from file (', &
2141		! & repeat_fits, &
2142		! & ') will be used.'
2143		endif
2144		repeat_tmp = repeat_fits
2145		ENDIF
2146
2147		ENDIF
2148
2149
2150		IF (npix .LT. nchunk_max) THEN ! data is small enough to be written in one chunk
2151
2152		frow = (firstpix_tmp)/repeat_tmp + 1
2153		felem = firstpix_tmp-(frow-1)*repeat_tmp + 1
2154		npix_32 = npix
2155
2156		DO colnum = 1, ntod
2157		call f90ftpcl_(unit, colnum, frow, felem, npix_32, tod(0:npix_32-1,colnum), status)
2158		END DO
2159
2160		ELSE ! data has to be written in several chunks
2161
2162		q = (npix-1)/nchunk_max
2163		DO iq = 0,q
2164		IF (iq .LT. q) THEN
2165		npix_tmp = nchunk_max
2166		ELSE
2167		npix_tmp = npix - iq*nchunk_max
2168		ENDIF
2169		i0 = iq*nchunk_max
2170		i1 = i0 + npix_tmp - 1_i8b
2171		firstpix_chunk = firstpix_tmp + i0
2172		frow = (firstpix_chunk)/repeat_tmp + 1
2173		felem = firstpix_chunk-(frow-1)*repeat_tmp + 1
2174		npix_32 = npix_tmp
2175		DO colnum = 1, ntod
2176		call f90ftpcl_(unit, colnum, frow, felem, npix_32, &
2177		& tod(i0:i1, colnum), status)
2178		END DO
2179		ENDDO
2180
2181		ENDIF
2182
2183		! pad entry if necessary ! 2015-04-28
2184		lastpix = firstpix_tmp + npix ! number of pixels written above
2185		lengap = modulo(lastpix, repeat_tmp) ! remaining gap
2186		if (lengap > 0) then
2187		firstpix_tmp = lastpix
2188		npix_32 = repeat_tmp - lengap
2189		frow = (firstpix_tmp)/repeat_tmp + 1
2190		felem = firstpix_tmp-(frow-1)*repeat_tmp + 1
2191		allocate(padding(0:npix_32-1))
2192		if (KMAP == SP) padding(:) = HPX_SBADVAL
2193		if (KMAP == DP) padding(:) = HPX_DBADVAL
2194		do colnum = 1, ntod
2195		call f90ftpcl_(unit, colnum, frow, felem, npix_32, padding, status)
2196		enddo
2197		deallocate(padding)
2198		endif
2199
2200		! ----------------------
2201		! close and exit
2202		! ----------------------
2203
2204		! close the file and free the unit number
2205		CALL ftclos(unit, status)
2206
2207		! check for any error, and if so print out error messages
2208		IF (status .GT. 0) CALL printerror(status)
2209
2210		RETURN
2211
2212		END SUBROUTINE write_bintabh_KLOAD
2213		! ==============================================================================
2214
2215		!******************************************************************************
2216		!---------------------------------------
2217		! reads a FITS file containing a list of
2218		! ring-based or pixel-based quadrature weights
2219		! and turns it into a (ring-ordered) full sky map
2220		!
2221		! adapted from unfold_weights.pro
2222		! 2018-05-25
2223		!---------------------------------------
2224		subroutine unfold_weightsfile_KLOAD(w8file, w8map)
2225
2226		USE pix_tools, only: nside2npix, nside2npweights
2227
2228		character(len=*), intent(in) :: w8file
2229		real(KMAP), intent(out), dimension(0:) :: w8map
2230		real(KMAP), allocatable, dimension(:,:) :: w8list
2231		integer(i4b) :: nside, status, won
2232		integer(i8b) :: nfw8, npw8, nrw8
2233		character(len=*), parameter :: code = 'unfold_weightsfile'
2234
2235		nfw8 = getsize_fits(w8file, nside=nside)
2236		npw8 = nside2npweights(nside)
2237		nrw8 = 2*nside
2238
2239		won = 0
2240		if (nfw8 == nrw8) won = 1
2241		if (nfw8 == npw8) won = 2
2242		if (won == 0) then
2243		print*,'Weights file = '//trim(w8file)
2244		print*,'contains ',nfw8,' weights'
2245		print*,'for Nside = ',nside
2246		print*,'Expected either ',nrw8,' or ',npw8
2247		call fatal_error
2248		endif
2249
2250		allocate(w8list(0:nfw8-1,1:1), stat=status)
2251		call assert_alloc(status,code,"w8list")
2252		call input_map(w8file, w8list, nfw8, 1)
2253		call unfold_weightslist(nside, won, w8list(:,1), w8map)
2254		deallocate(w8list)
2255
2256		return
2257		end subroutine unfold_weightsfile_KLOAD
2258		!---------------------------------------
2259		! turns a list of ring-based or pixel-based quadrature weights
2260		! into a (ring-ordered) full sky map
2261		!
2262		! adapted from unfold_weights.pro
2263		! 2018-05-18
2264		! 2018-09-12: promoted qpix to I8B to avoid consistency problems during compilation with g95
2265		!---------------------------------------
2266		subroutine unfold_weightslist_KLOAD(nside, won, w8list, w8map)
2267
2268		USE pix_tools, only: nside2npix, nside2npweights
2269		USE long_intrinsic, only: long_size
2270
2271		integer(i4b), intent(in) :: nside, won
2272		real(KMAP), dimension(0:), intent(in) :: w8list
2273		real(KMAP), dimension(0:), intent(out) :: w8map
2274
2275		integer(i8b) :: np, npix
2276		integer(i8b) :: nf, nw8, pnorth, psouth, vpix, qp4, p, qpix
2277		integer(i4b) :: ring, odd, shifted, rpix, wpix, j4, n4, it
2278
2279		! test nside
2280		npix = nside2npix(nside)
2281		if (npix <= 0) then
2282		print*,'Unvalid Nside = ',nside
2283		call fatal_error
2284		endif
2285
2286		! test won
2287		if (won < 1 .or. won > 2) then
2288		print*,'Expected either 1 or 2, got ',won
2289		call fatal_error
2290		endif
2291
2292		! test input list size
2293		nf = long_size(w8list)
2294		if (won == 1) then
2295		nw8 = 2 * nside
2296		else
2297		nw8 = nside2npweights(nside)
2298		endif
2299		if (nf /= nw8) then
2300		print*,'Expected in input weight list',nw8
2301		print*,'got ',nf
2302		call fatal_error
2303		endif
2304
2305		! test output map size
2306		np = long_size(w8map)
2307		if (np /= npix) then
2308		print*,'Expected in output weight map',npix
2309		print*,'got ',np
2310		call fatal_error
2311		endif
2312
2313		! do actual unfolding
2314		if (won == 1) then
2315
2316		pnorth = 0_i8b
2317		n4 = 4 * nside
2318		do ring=1,n4-1 ! loop on all rings
2319		it = min(ring, n4 - ring) ! ring index, counting from closest pole
2320		qp4= 4_i8b * min(it, nside) ! pixels / ring
2321		w8map(pnorth:pnorth+qp4-1) = w8list(it-1)
2322		pnorth = pnorth + qp4
2323		enddo
2324
2325		else
2326
2327		pnorth = 0_i8b ! position in expanded weights
2328		vpix = 0_i8b ! position in compress list
2329		do ring=0, 2*nside-1
2330		qpix = min(ring+1, nside) ! 1/4 of number of pixels per ring
2331		shifted = 0
2332		if (ring < (nside-1) .or. iand(ring+nside,1_i4b) == 1_i4b) shifted = 1
2333		odd = iand(qpix,1_i8b)
2334		qp4 = 4_i8b*qpix ! number of pixels per ring
2335		! fill the weight map
2336		do p=0,qp4-1
2337		j4 = mod(p, qpix) ! g95: p and qpix must be of same type
2338		rpix = min(j4, qpix - shifted - j4) ! seesaw
2339		w8map(pnorth+p) = w8list(vpix + rpix)
2340		enddo
2341
2342		if (ring < 2*nside-1) then ! south part (except equator)
2343		psouth = npix - pnorth - qp4
2344		w8map(psouth:psouth+qp4-1) = w8map(pnorth:pnorth+qp4-1)
2345		endif
2346		! locations on next ring
2347		pnorth = pnorth + qp4
2348		wpix = (qpix+1)/2 + 1 - ior(odd, shifted)
2349		vpix = vpix + wpix
2350		enddo
2351		endif
2352
2353		! add 1 to get final weight
2354		w8map = w8map + 1.0_KMAP
2355
2356		return
2357		end subroutine unfold_weightslist_KLOAD
2358

+122

-26

src/f90/mod/fitstools.F90 less more

83	83	! subroutine read_bintod
84	84	! subroutine write_bintabh
85	85	! subroutine unfold_weights
	86	! subroutine read_fits_partial
	87	! subroutine write_fits_partial
86	88	! ----------------------------------
87	89	!
88	90	! subroutine read_fits_cut4 ?

118	120	real(kind=SP), private, parameter :: s_bad_value = HPX_SBADVAL
119	121	real(kind=DP), private, parameter :: d_bad_value = HPX_DBADVAL
120	122	integer(kind=I4B), private, parameter :: i_bad_value = -1637500000
	123	integer(kind=I8B), private, parameter :: l_bad_value = -1637500000_I8B
121	124	integer(I4B) , private, parameter :: nchunk_max = 12000
122	125	integer(I4B), private, parameter :: MAXDIM_TOP = 199 ! < 999
123	126

166	169	#endif
167	170	end interface
168	171
	172	interface read_fits_partial
	173	#ifdef NO64BITS
	174	module procedure read_fits_partial8_s, read_fits_partial8_d
	175	#else
	176	module procedure read_fits_partial8_s, read_fits_partial8_d, &
	177	& read_fits_partial4_s, read_fits_partial4_d
	178	#endif
	179	end interface
	180
	181	interface write_fits_partial
	182	#ifdef NO64BITS
	183	module procedure write_fits_partial4_s, write_fits_partial4_d
	184	#else
	185	module procedure write_fits_partial8_s, write_fits_partial8_d, &
	186	& write_fits_partial4_s, write_fits_partial4_d
	187	#endif
	188	end interface
	189
169	190	!
170	191	interface input_tod
171	192	module procedure input_tod_s,input_tod_d

216	237	end interface
217	238
218	239	interface f90ftpcl_
219		module procedure f90ftpcle, f90ftpcld
	240	#ifdef NO64BITS
	241	module procedure f90ftpcle, f90ftpcld, f90ftpclj
	242	#else
	243	module procedure f90ftpcle, f90ftpcld, f90ftpclj, f90ftpclk
	244	#endif
220	245	end interface
221	246
222	247	interface f90ftgcv_
223		module procedure f90ftgcve, f90ftgcvd
	248	#ifdef NO64BITS
	249	module procedure f90ftgcve, f90ftgcvd, f90ftgcvj
	250	#else
	251	module procedure f90ftgcve, f90ftgcvd, f90ftgcvj, f90ftgcvk
	252	#endif
224	253	end interface
225	254
226	255	interface f90ftgpv_

228	257	end interface
229	258
230	259	interface f90ftgky_
231		module procedure f90ftgkye, f90ftgkyd
	260	#ifdef NO64BITS
	261	module procedure f90ftgkye, f90ftgkyd, f90ftgkyj
	262	#else
	263	module procedure f90ftgkye, f90ftgkyd, f90ftgkyj, f90ftgkyk
	264	#endif
232	265	end interface
233	266
234	267	interface map_bad_pixels

248	281	public :: read_fits_cut4, write_fits_cut4, &
249	282	& input_map, read_bintab, &
250	283	& output_map, write_bintab
	284	public :: read_fits_partial, write_fits_partial
251	285	public :: write_plm
252	286	public :: fits2cl, read_asctab, write_asctab
253	287	public :: read_dbintab, write_dbintab

261	295	contains
262	296
263	297	!-------------------------------------------------------------------------------
264		! generic interface F90FTPCL_ for FITSIO's FTPCLE and FTPCLD
	298	! generic interface F90FTPCL_ for FITSIO's FTPCL[E,D,J,K]
265	299	! writes data in ASCTAB or BINTAB
266	300	subroutine f90ftpcle(unit, colnum, frow, felem, np, data, status)
267	301	integer(I4B), intent(in) :: unit, colnum, frow, felem, np

277	311	call ftpcld(unit, colnum, frow, felem, np, data, status)
278	312	return
279	313	end subroutine f90ftpcld
	314	subroutine f90ftpclj(unit, colnum, frow, felem, np, data, status)
	315	integer(I4B), intent(in) :: unit, colnum, frow, felem, np
	316	integer(I4B), intent(out) :: status
	317	integer(I4B), intent(in), dimension(0:) :: data
	318	call ftpclj(unit, colnum, frow, felem, np, data, status)
	319	return
	320	end subroutine f90ftpclj
	321	#ifndef NO64BITS
	322	subroutine f90ftpclk(unit, colnum, frow, felem, np, data, status)
	323	integer(I4B), intent(in) :: unit, colnum, frow, felem, np
	324	integer(I4B), intent(out) :: status
	325	integer(I8B), intent(in), dimension(0:) :: data
	326	call ftpclk(unit, colnum, frow, felem, np, data, status)
	327	return
	328	end subroutine f90ftpclk
	329	#endif
280	330	!-------------------------------------------------------------------------------
281		! generic interface F90FTGCV_ for FITSIO's FTGCVE and FTGCVD
	331	! generic interface F90FTGCV_ for FITSIO's FTGCV[E,D,J,K]
282	332	! reads data from BINTAB
283	333	subroutine f90ftgcve(unit, colnum, frow, felem, np, nullval, data, anynull, status)
284	334	integer(I4B), intent(in) :: unit, colnum, frow, felem, np

298	348	call ftgcvd(unit, colnum, frow, felem, np, nullval, data, anynull, status)
299	349	return
300	350	end subroutine f90ftgcvd
	351	subroutine f90ftgcvj(unit, colnum, frow, felem, np, nullval, data, anynull, status)
	352	integer(I4B), intent(in) :: unit, colnum, frow, felem, np
	353	integer(I4B), intent(out) :: status
	354	logical(LGT), intent(out) :: anynull
	355	integer(I4B), intent(out), dimension(0:) :: data
	356	integer(I4B), intent(in) :: nullval
	357	call ftgcvj(unit, colnum, frow, felem, np, nullval, data, anynull, status)
	358	return
	359	end subroutine f90ftgcvj
	360	#ifndef NO64BITS
	361	subroutine f90ftgcvk(unit, colnum, frow, felem, np, nullval, data, anynull, status)
	362	integer(I4B), intent(in) :: unit, colnum, frow, felem, np
	363	integer(I4B), intent(out) :: status
	364	logical(LGT), intent(out) :: anynull
	365	integer(I8B), intent(out), dimension(0:) :: data
	366	integer(I8B), intent(in) :: nullval
	367	call ftgcvk(unit, colnum, frow, felem, np, nullval, data, anynull, status)
	368	return
	369	end subroutine f90ftgcvk
	370	#endif
301	371	!-------------------------------------------------------------------------------
302	372	! generic interface F90FTGPV_ for FITSIO's FTGPVE and FTGPVD
303	373	! reads data from IMAGE

320	390	return
321	391	end subroutine f90ftgpvd
322	392	!-------------------------------------------------------------------------------
323		! generic interface F90FTGKY_ for FITSIO's FTGKYE and FTGKYD
	393	! generic interface F90FTGKY_ for FITSIO's FTGKY[E,D,J,K]
324	394	! reads a keyword
325	395	subroutine f90ftgkye(unit, keyword, value, comment, status)
326	396	integer(I4B), intent(in) :: unit

340	410	call ftgkyd(unit, keyword, value, comment, status)
341	411	return
342	412	end subroutine f90ftgkyd
	413	subroutine f90ftgkyj(unit, keyword, value, comment, status)
	414	integer(I4B), intent(in) :: unit
	415	character(len=*), intent(in) :: keyword
	416	integer(I4B), intent(out) :: status
	417	character(len=*), intent(out) :: comment
	418	integer(I4B), intent(out) :: value
	419	call ftgkyj(unit, keyword, value, comment, status)
	420	return
	421	end subroutine f90ftgkyj
	422	#ifndef NO64BITS
	423	subroutine f90ftgkyk(unit, keyword, value, comment, status)
	424	integer(I4B), intent(in) :: unit
	425	character(len=*), intent(in) :: keyword
	426	integer(I4B), intent(out) :: status
	427	character(len=*), intent(out) :: comment
	428	integer(I8B), intent(out) :: value
	429	call ftgkyk(unit, keyword, value, comment, status)
	430	return
	431	end subroutine f90ftgkyk
	432	#endif
343	433	!-------------------------------------------------------------------------------
344	434
345	435
346	436	! define routine with SP I/O
347		#include "fits_s_inc.f90"
	437	#include "fits_s_inc.F90"
348	438
349	439	! define routine with DP I/O
350		#include "fits_d_inc.f90"
	440	#include "fits_d_inc.F90"
351	441
352	442
353	443

401	491	call ftgkyj(unit,'NAXIS', naxis, comment, status)
402	492	if (status > 0) call printerror(status)
403	493	if (naxis > 0) then ! there is an image
404		print*,'an image was found in the FITS file '//filename
	494	print*,'an image was found in the FITS file '//trim(filename)
405	495	print*,'... it is ignored.'
406	496	endif
407	497
408	498	! determines the presence of an extension
409	499	call ftgkyl(unit,'EXTEND', extend, comment, status)
410	500	if (status > 0) then
411		print*,'extension expected and not found in FITS file '//filename
	501	print*,'extension expected and not found in FITS file '//trim(filename)
412	502	print*,'abort code'
413	503	call fatal_error
414	504	endif

424	514	& nrows, tfields, ttype, tform, tunit, extname, varidat, &
425	515	& status)
426	516	if (tfields < 4) then
427		print*,'Expected 4 columns in FITS file '//filename
	517	print*,'Expected 4 columns in FITS file '//trim(filename)
428	518	print*,'found ',tfields
429	519	if (tfields < 2) call fatal_error
430	520	if (.not. (trim(ttype(1)) == 'PIXEL' &

613	703	! writes required keywords
614	704	! repeat = 1024
615	705	repeat = 1
	706	if (obs_npix < repeat) repeat = 1
616	707	nrows = (obs_npix + repeat - 1)/ repeat ! naxis1
617		if (obs_npix < repeat) repeat = 1
618	708	write(srepeat,'(i4)') repeat
619	709	srepeat = adjustl(srepeat)
620	710

632	722	tunit = ' ' ! optional, will not appear
633	723	tunit(2) = units_usr
634	724	tunit(4) = units_usr
635		extname = 'SKY_OBSERVATION' ! default, will be overide by user provided one if any
	725	extname = 'SKY_OBSERVATION' ! default, will be overridden by user provided one if any
636	726	if (polar_flag) extname = extnames(1+extno_i)
637	727	varidat = 0
638	728	call ftphbn(unit, nrows, tfields, ttype, tform, tunit, &

651	741	call ftpcom(unit,' Data Specific Keywords ',status)
652	742	call ftpcom(unit,'------------------------------------------',status)
653	743	call ftpkys(unit,'INDXSCHM','EXPLICIT',' Indexing : IMPLICIT or EXPLICIT', status)
654		call ftpkyj(unit,'GRAIN', grain, ' Grain of pixel indexing',status)
655		call ftpcom(unit,'GRAIN=0 : no indexing of pixel data (IMPLICIT) ',status)
656		call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)',status)
657		call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
	744	! call ftpkyj(unit,'GRAIN', grain, ' Grain of pixel indexing',status)
	745	! call ftpcom(unit,'GRAIN=0 : no indexing of pixel data (IMPLICIT)',status)
	746	! call ftpcom(unit,'GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)',status)
	747	! call ftpcom(unit,'GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)',status)
658	748	call ftpkys(unit,'OBJECT','PARTIAL ',' Sky coverage represented by data',status)
659	749	call ftpkyj(unit,'OBS_NPIX',obs_npix, ' Number of pixels observed and recorded',status)
660	750

715	805	endif
716	806
717	807
718		! write the user provided header literally, except for PIXTYPE, TFORM, TTYPE, TUNIT* and INDXSCHM
	808	! write the user provided header literally, except for PIXTYPE, TFORM, TTYPE, TUNIT*, INDXSCHM and GRAIN
719	809	! copy NSIDE, ORDERING and COORDSYS and POLAR if they are valid and not already given
720	810	do i=1,nlheader
721	811	card = header(i)
722	812	if (card(1:5) == 'TTYPE' .or. card(1:5) == 'TFORM' .or. card(1:7) == 'PIXTYPE') then
723	813	continue ! don't keep them
724	814	else if (card(1:8) == 'INDXSCHM') then
	815	continue
	816	else if (card(1:5) == 'GRAIN') then ! already written above
	817	continue
	818	else if (card(1:13) == 'COMMENT GRAIN' .or. card(1:14) == 'COMMENT GRAIN') then ! already written above
725	819	continue
726	820	else if (card(1:5) == 'TUNIT') then
727	821	if ((card(6:6) == '2' .or. card(6:6) == '4') .and. trim(units_usr) == '') then

1967	2061	INTEGER(I4B) :: status,unit,readwrite,blocksize,naxis
1968	2062	CHARACTER(LEN=80) :: comment, ttype1
1969	2063	LOGICAL(LGT) :: extend, anyf
1970		INTEGER(I4B):: nmove, hdutype, idmax, nrows
	2064	INTEGER(I4B):: nmove, hdutype, nrows
	2065	integer(I4B), dimension(1) :: idmax ! because ftgcvj expects an array (2020-08-24)
1971	2066
1972	2067	!-----------------------------------------------------------------------
1973	2068	status=0

2014	2109	ttype1 = trim(strupcase(adjustl(ttype1)))
2015	2110	if (trim(ttype1(1:5)) == 'INDEX') then
2016	2111	call ftgkyj(unit, 'NAXIS2', nrows, comment, status) ! find number of rows
2017		call ftgcvj(unit, 1_i4b, nrows, 1_i4b, 1_i4b, 0_i4b, idmax, anyf, status) ! read element on last row of first column
	2112	call ftgcvj(unit, 1_i4b, nrows, 1_i4b, 1_i4b, 0_i4b, idmax(1), anyf, status) ! read element on last row of first column
2018	2113	if (status == 0) then
2019		lmax = int(sqrt( real(idmax-1, kind = DP) ) )
	2114	lmax = int(sqrt( real(idmax(1)-1, kind = DP) ) )
2020	2115	if (lmax > 0) goto 1000
2021	2116	endif
2022	2117	endif

2758	2853	character(len=*), parameter :: primer_url = 'http://healpix.sf.net/pdf/intro.pdf'
2759	2854	!====================================================================
2760	2855
	2856	!n_ext = getnumext_fits(mapfile)
2761	2857	npixtot = getsize_fits(mapfile, nmaps = nmaps, ordering=order_map, nside=nsmax,&
2762	2858	& mlpol=mlpol, type = type, polarisation = polar_fits, &
2763	2859	& coordsys=coordsys, polcconv=polcconv)
2764	2860
2765	2861	if (nsmax<=0) then
2766		print*,"Keyword NSIDE not found in FITS header!"
	2862	print*,"Keyword NSIDE not found in FITS header of "//trim(mapfile)//" !"
2767	2863	call fatal_error(code)
2768	2864	endif
2769	2865	if (type == 3) npixtot = nside2npix(nsmax) ! cut sky input data set
2770	2866	if (nsmax/=npix2nside(npixtot)) then
2771	2867	print 9000,"FITS header keyword NSIDE does not correspond"
2772		print 9000,"to the size of the map!"
	2868	print 9000,"to the size of the map in "//trim(mapfile)//" !"
2773	2869	call fatal_error(code)
2774	2870	endif
2775	2871
2776	2872	if (polarisation .and. (nmaps >=3) .and. polar_fits == -1) then
2777		print 9000,"The input fits file MAY NOT contain polarisation data."
	2873	print 9000,"The input FITS file "//trim(mapfile)//" MAY NOT contain polarisation data."
2778	2874	print 9000,"Proceed at your own risk"
2779	2875	endif
2780	2876
2781	2877	if (polarisation .and. (nmaps<3 .or. polar_fits ==0)) then
2782		print 9000,"The file does NOT contain polarisation maps"
	2878	print 9000,"The FITS file "//trim(mapfile)//" does NOT contain polarisation maps"
2783	2879	print 9000,"only the temperature field will be analyzed"
2784	2880	polarisation = .false.
2785	2881	endif

2808	2904	! --- check ordering scheme ---
2809	2905	if ((order_map/=1).and.(order_map/=2)) then
2810	2906	print 9000,"The ordering scheme of the map must be RING or NESTED."
2811		print 9000,"No ordering specification is given in the FITS-header!"
	2907	print 9000,"No ordering specification is given in the FITS-header of "//trim(mapfile)//" !"
2812	2908	call fatal_error(code)
2813	2909	endif
2814	2910

-3

src/f90/mod/gen_fits_code less more

2	2	# EH, Dec 2004 -- Jan 2005
3	3
4	4	SetDef(){
5		template_file=fits_template.f90
6		outfile_s=fits_s_inc.f90
7		outfile_d=fits_d_inc.f90
	5	template_file=fits_template.F90
	6	outfile_s=fits_s_inc.F90
	7	outfile_d=fits_d_inc.F90
8	8	}
9	9
10	10	#--------

+52

-23

src/f90/mod/head_fits.F90 less more

523	523	endif
524	524	endif
525	525	enddo
526		if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
	526	if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
527	527	if (present(comment)) comment = ' '
528	528	if (present(count)) count = count_in
529	529	return

562	562	endif
563	563	endif
564	564	enddo
565		if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
	565	if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
566	566	if (present(comment)) comment = ' '
567	567	if (present(count)) count = count_in
568	568	return

672	672	return
673	673	endif
674	674	enddo
675		if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
	675	if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
676	676	if (present(comment)) comment = ' '
677	677	if (present(count)) count = count_in
678	678	return

711	711	endif
712	712	endif
713	713	enddo
714		if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
	714	if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
715	715	if (present(comment)) comment = ' '
716	716	if (present(count)) count = count_in
717	717	return

752	752	endif
753	753	endif
754	754	enddo
755		if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
	755	if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
756	756	if (present(comment)) comment = ' '
757	757	if (present(count)) count = count_in
758	758	return

791	791	endif
792	792	endif
793	793	enddo
794		if (verbose) print*,' >>>>> keyword '//kwd//' not found <<<<< '
	794	if (verbose) print*,' >>>>> keyword '//trim(kwd)//' not found <<<<< '
795	795	if (present(comment)) comment = ' '
796	796	if (present(count)) count = count_in
797	797	return

877	877	character(len=80), dimension(1:20) :: my_units
878	878	character(len=8) :: my_ordering
879	879	character(len=FILENAMELEN) :: udtype, my_beam_leg
880		logical(LGT) :: do_full, do_alm, do_cl, do_cut, do_polar, &
881		& do_bcross, do_reset, do_asym
	880	logical(LGT) :: do_full, do_alm, do_cl, do_cut, do_partial, &
	881	& do_polar, do_bcross, do_reset, do_asym, do_map
882	882
883	883	character(len=*), parameter :: code = 'write_minimal_header'
884	884	!--------------------------------------------------------------------------------
885	885	! find out type of data to be contained in FITS file
886	886	udtype = strupcase(dtype)
887		do_alm = .false. ; do_cl = .false. ; do_full = .false. ; do_cut = .false.
	887	do_alm = .false. ; do_cl = .false. ; do_full = .false.
	888	do_cut = .false. ; do_partial = .false.
888	889	select case (trim(udtype))
889	890	case ('ALM')
890	891	do_alm = .true.

894	895	do_full = .true.
895	896	case ('CUTMAP')
896	897	do_cut = .true.
	898	! case ('PARTIAL')
	899	do_partial = .true.
897	900	case default
898	901	print*,'Invalid choice of datatype:'//trim(dtype)
899	902	print*,'Should be one of: ALM, alm, CL, cl, CUTMAP, cutmap, MAP, map'
900	903	call fatal_error(code)
901	904	end select
	905	do_map = (do_cut .or. do_full .or. do_partial)
902	906
903	907	!
904	908	if (present(order) .and. present(ordering)) then

907	911	endif
908	912
909	913	! check for required keywords for map
910		if (do_cut .or. do_full) then
	914	if (do_map) then
911	915	if (.not. present(nside)) then
912	916	print*,'NSIDE required for maps in '//code
913	917	call fatal_error(code)

1107	1111
1108	1112	endif ! cl
1109	1113
1110		if (do_full .or. do_cut) then
	1114	if (do_map) then
1111	1115	!--------------------------------------------------------------------------------------------
1112	1116	! minimal header for FULL SKY MAP and CUT SKY MAP
1113	1117	!--------------------------------------------------------------------------------------------

1155	1159	call add_card(header,"COMMENT","Full sky data")
1156	1160	call add_card(header,"OBJECT","FULLSKY")
1157	1161	call add_card(header,"INDXSCHM","IMPLICIT"," Indexing : IMPLICIT or EXPLICIT")
1158		call add_card(header,"GRAIN", 0, " Grain of pixel indexing")
1159		endif
1160		if (do_cut) then
	1162	! call add_card(header,"GRAIN", 0, " Grain of pixel indexing")
	1163	endif
	1164	if (do_cut .or. do_partial) then
1161	1165	call add_card(header,"EXTNAME","'CUT SKY MAP'",update=.true.)
1162	1166	call add_card(header,"COMMENT","Cut sky data")
1163	1167	call add_card(header,"OBJECT","PARTIAL")
1164	1168	call add_card(header,"INDXSCHM","EXPLICIT"," Indexing : IMPLICIT or EXPLICIT")
1165		call add_card(header,"GRAIN", 1, " Grain of pixel indexing")
1166		endif
1167		call add_card(header,"COMMENT","GRAIN=0 : no indexing of pixel data (IMPLICIT)")
1168		call add_card(header,"COMMENT","GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)")
1169		call add_card(header,"COMMENT","GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)")
	1169	! call add_card(header,"GRAIN", 1, " Grain of pixel indexing")
	1170	endif
	1171	! call add_card(header,"COMMENT","GRAIN=0 : no indexing of pixel data (IMPLICIT)")
	1172	! call add_card(header,"COMMENT","GRAIN=1 : 1 pixel index -> 1 pixel data (EXPLICIT)")
	1173	! call add_card(header,"COMMENT","GRAIN>1 : 1 pixel index -> data of GRAIN consecutive pixels (EXPLICIT)")
1170	1174	call add_card(header) ! blank line
1171	1175	call add_card(header,"POLAR",do_polar," Polarisation included (True/False)")
1172	1176	if (do_cut) then
	1177	call add_card(header) ! blank line
	1178	call add_card(header,"TTYPE1", "PIXEL","pixel index")
	1179	call add_card(header,"TUNIT1", "", "")
	1180	call add_card(header)
1173	1181	if (present(units)) then
1174	1182	my_units = units
1175		call add_card(header,"TUNIT2", my_units(1),"physical unit of signal map")
1176		call add_card(header,"TUNIT4", my_units(1),"physical unit of error map")
	1183	call add_card(header,"TUNIT2", my_units(1),"physical units of map")
	1184	call add_card(header,"TUNIT4", my_units(1),"physical units of map")
1177	1185	endif
1178	1186	endif
	1187	! if (do_partial) then
	1188	! if (present(units)) my_units = units
	1189	! it = 1 ; iq = 1 ; iu = 1
	1190	! call add_card(header) ! blank line
	1191	! call add_card(header,"TTYPE1", "PIXEL","pixel index")
	1192	! call add_card(header,"TUNIT1", "", "")
	1193	! call add_card(header)
	1194
	1195	! call add_card(header) ! blank line
	1196	! call add_card(header,"TTYPE2", "TEMPERATURE","Temperature map")
	1197	! call add_card(header,"TUNIT2", my_units(it), "physical units of map")
	1198	! call add_card(header)
	1199
	1200	! if (do_polar) then
	1201	! call add_card(header,"TTYPE3", "Q_POLARISATION","Q Polarisation map")
	1202	! call add_card(header,"TUNIT3", my_units(iq), "physical units of map")
	1203	! call add_card(header)
	1204
	1205	! call add_card(header,"TTYPE4", "U_POLARISATION","U Polarisation map")
	1206	! call add_card(header,"TUNIT4", my_units(iu), "physical units unit")
	1207	! call add_card(header)
	1208	! endif
	1209	! endif
1179	1210	if (do_full) then
1180	1211	call add_card(header,"DERIV",my_deriv," Derivative included (0, 1 or 2)")
1181	1212

1198	1229	call add_card(header)
1199	1230
1200	1231	if (do_polar) then
1201		!call add_card(header,"TTYPE2", "Q-POLARISATION","Q Polarisation map")
1202	1232	call add_card(header,"TTYPE2", "Q_POLARISATION","Q Polarisation map")
1203	1233	call add_card(header,"TUNIT2", my_units(iq),"map unit")
1204	1234	call add_card(header)
1205	1235
1206		!call add_card(header,"TTYPE3", "U-POLARISATION","U Polarisation map")
1207	1236	call add_card(header,"TTYPE3", "U_POLARISATION","U Polarisation map")
1208	1237	call add_card(header,"TUNIT3", my_units(iu),"map unit")
1209	1238	call add_card(header)

-1

src/f90/mod/healpix_types.F90 less more

38	38	! Mars 2008: i8b same as i4b on machines not supporting 64 bits (NO64BITS flag set)
39	39	! Feb 2009: introduce healpix_version
40	40	!
41		character(len=*), PARAMETER, public :: healpix_version = '3.60'
	41	character(len=*), PARAMETER, public :: healpix_version = '3.81'
42	42	INTEGER, PARAMETER, public :: i4b = SELECTED_INT_KIND(9)
43	43	#ifdef NO64BITS
44	44	INTEGER, PARAMETER, public :: i8b = i4b

+22

-9

src/f90/mod/pix_tools.F90 less more

870	870	integer(i4b), dimension(1:,1:), intent(inout) :: ringphi
871	871	integer(i4b), intent(inout) :: ngr
872	872
873		integer(i4b) :: i, j, k, kk, ngr_out, diff, iphi, i0, nrh
	873	!integer(i4b) :: i, j, k, kk, ngr_out, diff, iphi, i0, nrh
	874	integer(i4b) :: i, j, k, kk, ngr_out, iphi, i0, nrh, nshrink, np
874	875	real(dp), dimension(1:2*nsboost+1) :: phiw, phie
875	876	real(dp) :: dd, dph, phic
	877	! 2020-08-31: bug correction
876	878
877	879	!=======================================================================
878	880	if (nsboost <= 1) return

881	883	do i=1, ngr ! loop on low-res rings
882	884	i0 = ringphi(1,i) * nsboost - nsboost - irmin
883	885	do k=-1,1,2 ! West and East side of disc
	886	nshrink = 0
884	887	kk = (k+5)/2 ! 2 or 3
885	888	222 continue
886	889	iphi = ringphi(kk, i)
887		if (ringphi(2,i) <= ringphi(3,i) .and. iphi >= 0) then
888		call find_pixel_bounds(nside, nsboost, ringphi(1,i), iphi, phiw, phie)
	890	!if (ringphi(2,i) <= ringphi(3,i) .and. iphi >= 0) then
	891	if (iphi >= 0) then ! corrected 2020-08-31
	892	call find_pixel_bounds(nside, nsboost, ringphi(1,i), iphi, phiw, phie, np=np)
889	893	do j=1, 2*nsboost+1
890	894	if (i0+j >= 1 .and. i0+j <= nrh) then
891	895	phic = (phie(j)+phiw(j))*0.5_dp ! pixel center

895	899	if (dd <= dph) goto 1000 ! pixel touched by disc, move to next one
896	900	endif
897	901	enddo
898		ringphi(kk, i)= iphi - k ! pixel not in disc, move edge pixel inwards
899		goto 222 ! try next pixel inward
	902	!ringphi(kk, i)= iphi - k ! pixel not in disc, move edge pixel inwards
	903	ringphi(kk, i) = mod(iphi + np - k, np) ! pixel not in disc, move edge pixel inwards in [0,np-1]
	904	nshrink = nshrink + 1
	905	if (nshrink <= 2) then
	906	goto 222 ! try next pixel inward
	907	else
	908	ringphi(kk,i) = -999
	909	endif
900	910	1000 continue
901	911	endif
902		enddo ! loop on side
	912	enddo ! loop on k East-West sides
903	913	enddo ! loop on low-res rings
904	914
905	915	! remove empty rings
906	916	ngr_out = 0
907	917	do i=1,ngr
908		diff = ringphi(3,i) - ringphi(2,i)
909		if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0 .and. diff /= -2 .and. diff /= -1) then
	918	! diff = ringphi(3,i) - ringphi(2,i)
	919	! if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0 .and. diff /= -2 .and. diff /= -1) then
	920	if (ringphi(2,i) >=0 .and. ringphi(3,i) >=0) then ! 2020-08-31
910	921	ngr_out = ngr_out + 1
911	922	ringphi(1:3, ngr_out) = ringphi(1:3, i)
912	923	endif

920	931	return
921	932	end subroutine check_edge_pixels
922	933	!=======================================================================
923		subroutine find_pixel_bounds (nside, nsboost, iring, iphi, phiw, phie)
	934	subroutine find_pixel_bounds (nside, nsboost, iring, iphi, phiw, phie, np)
924	935	!=======================================================================
925	936	integer(i4b), intent(in) :: nside, nsboost, iring, iphi
926	937	real(dp), dimension(1:2*nsboost+1), intent(out) :: phiw, phie
	938	integer(i4b), optional, intent(out) :: np
927	939
928	940	real(dp), dimension(1:2*nsboost+1) :: f, f1, phiw_t, phie_t
929	941	real(dp) :: c0, quad, phie1, phie2, phiw1, phiw2, cv

935	947	!f = ((/ (i,i=0,2*nsboost) /) - nsboost) / nsboost
936	948	f = ((/ (i1.d0,i=0,2nsboost) /) - nsboost*1.d0) / nsboost
937	949
	950	if (present(np)) np = npr
938	951	nq = npr/4 ! number of pixels on current ring in [0,Pi/2] (quadrant)
939	952	transition = (iring == nside .or. iring == nside*3)
940	953

-2

src/f90/mod/pixel_routines.F90 less more

1919	1919	! New algorithm for Inclusive case: test boundary of edge pixels on each ring
1920	1920	! 2013-04-02: bug correction in query_disc in inclusive mode
1921	1921	! 2014-07-07: bug correction in discedge2fulldisc
	1922	! 2020-08-31: bug correction in check_edge_pixels (reported by R. de Belsunce)
1922	1923	!=======================================================================
1923	1924	#ifdef DOI8B
1924	1925	subroutine query_disc_8( nside, vector0, radius, listpix, nlist, nest, inclusive)

2015	2016	call pixels_on_edge(nside, irmin, irmax, phi0, dphitab, ringphi, ngr)
2016	2017	if (do_inclusive) then
2017	2018	! sample edge pixels at larger Nside
2018		nsboost = 16
2019		nsideh = min(NS_MAX8, nside * int(nsboost,i8b))
	2019	nsboost = 2*ceiling( log(1.1_DP/(radiusnside))/log(2.0_DP) ) ! make pixel smaller than disc (2020-08-31)
	2020	nsboost = max(nsboost,16)
	2021	nsideh = min(NS_MAX8, nside * int(nsboost,i8b))
2020	2022	radiush = fudge_query_radius(nsideh, radius, quadratic=.true.)
2021	2023
2022	2024	irmin = ring_num(nsideh, zmax, shift=+1) ! shifted South

2027	2029	zlist(iz-irmin+1) = ring2z(nsideh, iz)
2028	2030	enddo
2029	2031	call discphirange_at_z(vector0, radiush, zlist, nrh, dphilist, phi0)
	2032	phi0 = mod(phi0 + TWOPI, TWOPI) ! [-Pi,Pi] -> [0,2Pi]
2030	2033	call check_edge_pixels(nside, nsboost, irmin, irmax, phi0, dphilist, ringphi, ngr)
2031	2034	deallocate(zlist, dphilist)
2032	2035	endif

-0

src/f90/mod/sharp_healpix_f_inc.c less more

43	43	void mapptr[2], almptr[2];
44	44
45	45	CHECK_STACK_ALIGN(8);
	46	UTIL_ASSERT(spin>0, "spin must not be zero");
46	47	sharp_make_healpix_geom_info_2 (nside, wgt, zbounds[0], zbounds[1], &ginfo);
47	48	sharp_make_rectangular_alm_info (lmax,mmax,2,&ainfo);
48	49

108	109	void mapptr[2], almptr[2];
109	110
110	111	CHECK_STACK_ALIGN(8);
	112	UTIL_ASSERT(spin>0, "spin must not be zero");
111	113	// sharp_make_healpix_geom_info (nside, 1, &ginfo);
112	114	int i,nrings=4*nside-1;
113	115	double *wgt=RALLOC(double,nrings);

+68

-0

src/f90/mod/udgrade_nr.F90 less more

	0	!-----------------------------------------------------------------------------
	1	!
	2	! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon,
	3	! Benjamin D. Wandelt, Anthony J. Banday,
	4	! Matthias Bartelmann, Hans K. Eriksen,
	5	! Frode K. Hansen, Martin Reinecke
	6	!
	7	!
	8	! This file is part of HEALPix.
	9	!
	10	! HEALPix is free software; you can redistribute it and/or modify
	11	! it under the terms of the GNU General Public License as published by
	12	! the Free Software Foundation; either version 2 of the License, or
	13	! (at your option) any later version.
	14	!
	15	! HEALPix is distributed in the hope that it will be useful,
	16	! but WITHOUT ANY WARRANTY; without even the implied warranty of
	17	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	18	! GNU General Public License for more details.
	19	!
	20	! You should have received a copy of the GNU General Public License
	21	! along with HEALPix; if not, write to the Free Software
	22	! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
	23	!
	24	! For more information about HEALPix see http://healpix.sourceforge.net
	25	!
	26	!-----------------------------------------------------------------------------
	27	module udgrade_nr
	28
	29	! --- in udgrade_template.f90 ---
	30	! subroutine sub_udgrade_nest
	31	! subroutine udgrade_ring
	32	! subroutine udgrade_nest
	33
	34	use healpix_types
	35	use misc_utils
	36
	37	IMPLICIT none
	38
	39	private
	40
	41	public :: sub_udgrade_nest, udgrade_ring, udgrade_nest
	42
	43	! SP and DP interface
	44	interface sub_udgrade_nest
	45	module procedure sub_udgrade_nest_s, sub_udgrade_nest_d
	46	end interface
	47
	48	! SP and DP, 1-dim and N-dim interface
	49	interface udgrade_ring
	50	module procedure udgrade_ring_1d_s, udgrade_ring_1d_d, udgrade_ring_nd_s, udgrade_ring_nd_d
	51	end interface
	52
	53	interface udgrade_nest
	54	module procedure udgrade_nest_1d_s, udgrade_nest_1d_d, udgrade_nest_nd_s, udgrade_nest_nd_d
	55	end interface
	56	!---------------------
	57
	58	contains
	59
	60	! include SP implementation of routines
	61	#include "udgrade_s_inc.f90"
	62
	63	! include DP implementation of routines
	64	#include "udgrade_d_inc.f90"
	65
	66
	67	end module udgrade_nr

-68

~~src/f90/mod/udgrade_nr.f90~~ less more

0		!-----------------------------------------------------------------------------
1		!
2		! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon,
3		! Benjamin D. Wandelt, Anthony J. Banday,
4		! Matthias Bartelmann, Hans K. Eriksen,
5		! Frode K. Hansen, Martin Reinecke
6		!
7		!
8		! This file is part of HEALPix.
9		!
10		! HEALPix is free software; you can redistribute it and/or modify
11		! it under the terms of the GNU General Public License as published by
12		! the Free Software Foundation; either version 2 of the License, or
13		! (at your option) any later version.
14		!
15		! HEALPix is distributed in the hope that it will be useful,
16		! but WITHOUT ANY WARRANTY; without even the implied warranty of
17		! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18		! GNU General Public License for more details.
19		!
20		! You should have received a copy of the GNU General Public License
21		! along with HEALPix; if not, write to the Free Software
22		! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
23		!
24		! For more information about HEALPix see http://healpix.sourceforge.net
25		!
26		!-----------------------------------------------------------------------------
27		module udgrade_nr
28
29		! --- in udgrade_template.f90 ---
30		! subroutine sub_udgrade_nest
31		! subroutine udgrade_ring
32		! subroutine udgrade_nest
33
34		use healpix_types
35		use misc_utils
36
37		IMPLICIT none
38
39		private
40
41		public :: sub_udgrade_nest, udgrade_ring, udgrade_nest
42
43		! SP and DP interface
44		interface sub_udgrade_nest
45		module procedure sub_udgrade_nest_s, sub_udgrade_nest_d
46		end interface
47
48		! SP and DP, 1-dim and N-dim interface
49		interface udgrade_ring
50		module procedure udgrade_ring_1d_s, udgrade_ring_1d_d, udgrade_ring_nd_s, udgrade_ring_nd_d
51		end interface
52
53		interface udgrade_nest
54		module procedure udgrade_nest_1d_s, udgrade_nest_1d_d, udgrade_nest_nd_s, udgrade_nest_nd_d
55		end interface
56		!---------------------
57
58		contains
59
60		! include SP implementation of routines
61		include 'udgrade_s_inc.f90'
62
63		! include DP implementation of routines
64		include 'udgrade_d_inc.f90'
65
66
67		end module udgrade_nr