1!----------------------------------------------------------------------------- 2! 3! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon, 4! Benjamin D. Wandelt, Anthony J. Banday, 5! Matthias Bartelmann, Hans K. Eriksen, 6! Frode K. Hansen, Martin Reinecke 7! 8! 9! This file is part of HEALPix. 10! 11! HEALPix is free software; you can redistribute it and/or modify 12! it under the terms of the GNU General Public License as published by 13! the Free Software Foundation; either version 2 of the License, or 14! (at your option) any later version. 15! 16! HEALPix is distributed in the hope that it will be useful, 17! but WITHOUT ANY WARRANTY; without even the implied warranty of 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19! GNU General Public License for more details. 20! 21! You should have received a copy of the GNU General Public License 22! along with HEALPix; if not, write to the Free Software 23! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 24! 25! For more information about HEALPix see http://healpix.sourceforge.net 26! 27!----------------------------------------------------------------------------- 28Module sky_sub 29 30 !Subroutines used by sky_ng_sim and sky_ng_sim_bin 31 32Contains 33 34 !-------------------------------------------------------- 35 36 Subroutine add_subscript(text, number) 37 38 Implicit none 39 40 Character(len = 20), Intent(InOut) :: text 41 Integer, Intent(In) :: number 42 43 Integer i 44 Character(len = 10) :: form 45 46 i = Ceiling(log10(number + 0.5)) 47 If (i .ge. 10) Stop "A ridiculously large number of bins!" 48 49 Write (form, "(a,i1,a)") "(a,a,i",i,")" 50 51 Write (text, form) TRIM(text), "_", number 52 53 End Subroutine add_subscript 54 55 !-------------------------------------------------------- 56 57 Subroutine read_powerspec(infile, nsmax, nlmax, cls, header_PS, fwhm, units_power, winfile, beam_file) 58 59 Use healpix_types 60 USE fitstools, ONLY : fits2cl 61 Use alm_tools, Only : generate_beam, pixel_window 62 use utilities, only : die_alloc 63 Use head_fits, only : add_card, get_card, merge_headers 64 65 Implicit none 66 67 Character(len = *), Intent(In) :: infile 68 Integer(I4B), Intent(In) :: nsmax 69 Integer(I4B), Intent(InOut) :: nlmax 70 Real(SP), dimension(0:nlmax,1:1), Intent(Out) :: cls 71 Character(Len=80), dimension(1:), Intent(Out) :: header_PS 72 Real(DP) :: fwhm 73 Character(len=80), dimension(1:1), Intent(Out) :: units_power 74 Character(len=*), Intent(In), optional :: winfile 75 Character(len=*), Intent(In), optional :: beam_file 76 77 78 Logical fitscl 79 Integer(I4B) :: nlheader 80 Integer(I4B) :: l_max, npw 81 Real(DP) :: lbeam 82 Integer i 83 Integer(I4B) :: status, count_tt = 0 84 Real(SP) :: quadrupole 85 REAL(DP), DIMENSION(:,:), ALLOCATABLE :: pixlw, beamw 86 CHARACTER(LEN=*), PARAMETER :: code = 'READ_POWERSPEC' 87 CHARACTER(LEN=80), DIMENSION(1:180) :: header_file 88 CHARACTER(LEN=20) :: string_quad 89 ! CHARACTER(LEN=80), DIMENSION(:), allocatable :: units_power 90 CHARACTER(LEN=80) :: temptype, com_tt 91 92 ! set maximum multipole (depends on user choice and pixel window function) 93 ! if using pixel window function l_max can't be greater than 4*nsmax 94 l_max = nlmax 95 npw = 4*nsmax + 1 96 if (present(winfile)) then 97 l_max = min(l_max, npw - 1) 98 else 99 npw = l_max + 1 100 endif 101 if (l_max < nlmax) then 102 print*,'C_l are only non-zero for 0 < l <= ',l_max 103 endif 104 nlmax = l_max 105 106 !----------------------------------------------------------------------- 107 108 cls = 0.0 109 nlheader = SIZE(header_file) 110 header_file = '' 111 ! Find out whether input cl file is a fits file 112 fitscl = (index(infile, '.fits') /= 0) 113 If (fitscl) Then 114!!! call read_asctab(infile, cls, nlmax, 1, header_file, nlheader, units=units_power) 115 call fits2cl(infile, cls, nlmax, 1, header_file, units=units_power) 116 call get_card(header_file,"TEMPTYPE",temptype,com_tt,count=count_tt) 117 Else 118 call read_camb(infile, cls, nlmax) 119 units_power = 'microKelvin^2' 120 !Check to see if enough values were present in file 121 If (nlmax .LT. l_max) then 122 Write (*,*) "Maximum l-value in cl file is ",nlmax 123 Write (*,*) "Power spectrum above this value set to zero" 124 l_max = nlmax 125 end if 126 End If 127 128 if (count_tt < 1) then 129 temptype = "THERMO" 130 com_tt = " temperature : either THERMO or ANTENNA" 131 print*,' Will assume '//temptype 132 endif 133 134 quadrupole = cls(2,1) 135 ! --- creates the header relative to power spectrum --- 136 header_PS = '' 137 call add_card(header_PS) 138 call add_card(header_PS,'HISTORY',' alm generated from following power spectrum') 139 call add_card(header_PS) 140 call add_card(header_PS,'COMMENT','----------------------------------------------------') 141 call add_card(header_PS,'COMMENT','Planck Power Spectrum Description Specific Keywords') 142 call add_card(header_PS,'COMMENT','----------------------------------------------------') 143 call add_card(header_PS,'COMMENT','Input power spectrum in : ') 144 call add_card(header_PS,'COMMENT',TRIM(infile)) 145 call add_card(header_PS) 146 call add_card(header_PS,'COMMENT','Quadrupole') 147 write(string_quad,'(1pe15.6)') quadrupole 148 call add_card(header_PS,'COMMENT',' C(2) = '//string_quad//' '//trim(units_power(1))) 149 call add_card(header_PS) 150 call add_card(header_PS,"TEMPTYPE",temptype,com_tt) 151 call add_card(header_PS) 152 ! ------- insert header read in power spectrum file ------- 153 call merge_headers(header_file, header_PS) 154 ! deallocate(units_power) 155 156 ! define beam 157 allocate(beamw(0:l_max,1:1), stat = status) 158 if (status /= 0) call die_alloc(code,'beamw') 159 If ((fwhm .NE. 0d0) .OR. present(beam_file)) Then 160 call generate_beam(real(fwhm,kind=dp), l_max, beamw, beam_file) 161 Else 162 beamw(:,:) = 1.0_dp 163 End If 164 165 ! get the pixel window function 166 allocate(pixlw(0:npw-1,1:1), stat = status) 167 if (status /= 0) call die_alloc(code,'pixlw') 168 if(present(winfile)) then 169 call pixel_window(pixlw, windowfile=winfile) 170 else 171 pixlw(:,:) = 1.0_dp 172 endif 173 174 ! Multiply by beam and pixel window functions 175! cls(0:l_max, 1) = cls(0:l_max,1)*beamw(0:l_max,1)*pixlw(0:l_max,1) !EH-2008-03-05 176 cls(0:l_max, 1) = cls(0:l_max,1)* (beamw(0:l_max,1)*pixlw(0:l_max,1))**2 177 178 deallocate(beamw) 179 deallocate(pixlw) 180 181 End Subroutine read_powerspec 182 183 !-------------------------------------------------------- 184 185 Subroutine read_camb(infile, cl_array, lmax) 186 187 Use healpix_types 188 189 Implicit none 190 191 Character (len = *), intent(IN) :: infile 192 Integer(I4B), intent(INOUT) :: lmax 193 Real(SP), Intent(OUT), Dimension(0:lmax,1:1) :: cl_array 194 195 Real(DP) :: clval 196 Integer :: lval = 0 197 Real(DP) :: cmbT = 2.726d0 198 199 ! Read in the cl values and convert to units of uK^2 200 Open (Unit = 23, file = infile, status = 'old', action = 'read', err = 9980) 201 Do while (lval .LT. lmax) 202 Read (23, *, end = 800) lval, clval 203 cl_array(lval, 1) = TWOPI*clval/(lval*(lval+1))*cmbT*cmbT*1.d12 204 End Do 205800 lmax = lval 206 207 Close (23) 208 Return 209 2109980 Stop 'Error opening cl file to read' 211 212 End Subroutine read_camb 213 214 215End Module sky_sub 216