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