Program comb_csim

program comb_csim

        ! Uses
    use s2_types_mod, only: S2_STRING_LEN, s2_sp, pi
    use s2_error_mod
    use s2_distn_mod
    use s2_cmb_mod
    use s2_pl_mod
    use s2_vect_mod, only: s2_vect_arcmin_to_rad
    use s2_wnoise_mod
    use comb_error_mod
    use comb_tmpl_mod
    use comb_obj_mod
    use comb_csky_mod
    use extension, only: getArgument, nArguments
    use paramfile_io, only: paramfile_handle, parse_init, parse_int,     parse_real, parse_lgt, parse_string, concatnl
    use pix_tools, only: nside2npix

        ! Calls
    call comb_csim_param_gen
    call comb_csim_param_read
    call comb_csim_print_header
    call comb_csky_free
    call comb_csky_write_param
    call comb_csky_write_sky_full
    call comb_csky_write_sky_obj
    call comb_obj_free
    call comb_obj_write_sky
    call s2_cmb_free
    call s2_cmb_write_sky
    call s2_error
    call s2_pl_free
    call s2_wnoise_downsample
    call s2_wnoise_free
    call s2_wnoise_write_nobs_file
    call s2_wnoise_write_sky_file
    call s2_wnoise_write_std_file

        ! Variables
    character (len=*), parameter :: CODE_NAME = 'COMB_CSIM'
    character (len=*), parameter :: CODE_DESCRIPTION = 'Compact object embedded CMB simulation'
    character (len=*), parameter :: CODE_VERSION = '1.0'
    character (len=*), parameter :: CODE_AUTHOR = 'Jason McEwen'
    character (len=S2_STRING_LEN) :: filename_param
    character (len=S2_STRING_LEN) :: description
    character (len=S2_STRING_LEN) :: line
    type (paramfile_handle) :: handle
    integer, parameter :: NUM_COMMENT_LINES_CLT_FILE = 29
    integer, parameter :: NUM_COMMENT_LINES_BEAM_FILE = 8
    real (kind=s2_sp) :: TOL = 1e-5
    integer :: nside
    integer :: lmax
    integer :: mmax
    integer, parameter :: lmin = 2
    integer, parameter :: NSIDE_DEFAULT = 64
    integer, parameter :: LMAX_DEFAULT = 128
    integer, parameter :: MMAX_DEFAULT = 128
    character (len=S2_STRING_LEN) :: comb_type
    integer, parameter :: COMB_TYPE_OPT_NUM = 6
    character (len=S2_STRING_LEN), parameter, dimension (COMB_TYPE_OPT_NUM) :: COMB_TYPE_OPT = (/ 'butterfly   ', 'gaussian    ', 'mexhat      ', 'morlet      ', 'point       ', 'cos_thetaon2' /)
    character (len=S2_STRING_LEN), parameter :: COMB_TYPE_DEFAULT = COMB_TYPE_OPT (1)
    integer :: comb_num
    integer, parameter, dimension (2) :: COMB_NUM_OPT = (/ 0, 50 /)
    integer, parameter :: COMB_NUM_DEFAULT = 10
    integer :: comb_seed
    real (kind=s2_sp) :: comb_amp_lower
    real (kind=s2_sp), parameter, dimension (2) :: COMB_AMP_LOWER_OPT = (/ 0.01e0, 5.0e0 /)
    real (kind=s2_sp), parameter :: COMB_AMP_LOWER_DEFAULT = 0.1e0
    real (kind=s2_sp) :: comb_amp_upper
    real (kind=s2_sp), parameter, dimension (2) :: COMB_AMP_UPPER_OPT = (/ 0.01e0, 5.0e0 /)
    real (kind=s2_sp), parameter :: COMB_AMP_UPPER_DEFAULT = 2.0e0
    real (kind=s2_sp) :: comb_dil_lower
    real (kind=s2_sp), parameter, dimension (2) :: COMB_DIL_LOWER_OPT = (/ 0.01e0, 1.0e0 /)
    real (kind=s2_sp), parameter :: COMB_DIL_LOWER_DEFAULT = 0.05e0
    real (kind=s2_sp) :: comb_dil_upper
    real (kind=s2_sp), parameter, dimension (2) :: COMB_DIL_UPPER_OPT = (/ 0.01e0, 1.0e0 /)
    real (kind=s2_sp), parameter :: COMB_DIL_UPPER_DEFAULT = 0.1e0
    integer :: comb_dil_grid_num
    integer, parameter :: COMB_DIL_GRID_NUM_DEFAULT = 5
    integer :: comb_gamma_grid_num
    integer, parameter :: COMB_GAMMA_GRID_NUM_DEFAULT = 5
    logical :: cmb_status
    logical, parameter :: CMB_STATUS_DEFAULT = .true.
    integer :: cmb_seed
    character (len=S2_STRING_LEN) :: cmb_cl_file
    character (len=S2_STRING_LEN), parameter :: CMB_CL_FILE_DEFAULT = 'data_in/wmap_lcdm_pl_model_yr1_v1.txt'
    logical :: beam_status
    logical, parameter :: BEAM_STATUS_DEFAULT = .false.
    logical :: beam_gaussian
    logical, parameter :: BEAM_GAUSSIAN_DEFAULT = .true.
    real (kind=s2_sp) :: beam_fwhm
    real (kind=s2_sp), parameter :: BEAM_FWHM_DEFAULT = 10.0e0
    character (len=S2_STRING_LEN) :: beam_cl_file
    character (len=S2_STRING_LEN), parameter :: BEAM_CL_FILE_DEFAULT = 'data_in/map_w1_ampl_bl_yr1_v1.txt'
    logical :: noise_status
    logical, parameter :: NOISE_STATUS_DEFAULT = .false.
    logical :: noise_full_sky
    logical, parameter :: NOISE_FULL_SKY_DEFAULT = .false.
    integer :: noise_seed
    real (kind=s2_sp) :: noise_std
    real (kind=s2_sp), parameter :: NOISE_STD_DEFAULT = 0.1e0
    real (kind=s2_sp) :: noise_sigma0
    real (kind=s2_sp), parameter :: NOISE_SIGMA0_DEFAULT = 1.0e0
    character (len=S2_STRING_LEN) :: noise_nobs_file
    character (len=S2_STRING_LEN), parameter :: NOISE_NOBS_FILE_DEFAULT = 'data_in/wmap_w1_cleanimap_yr1_v1.fits'
    integer :: noise_nobs_extension
    integer, parameter :: NOISE_NOBS_EXTENSION_DEFAULT = 2
    character (len=S2_STRING_LEN) :: output_file_csky
    character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_CSKY_DEFAULT = 'data_out/csky.fits'
    character (len=S2_STRING_LEN) :: output_file_obj
    character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_OBJ_DEFAULT = 'data_out/obj.fits'
    character (len=S2_STRING_LEN) :: output_file_obj_mother
    character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_OBJ_MOTHER_DEFAULT = 'data_out/obj_mother.fits'
    character (len=S2_STRING_LEN) :: output_file_cmb
    character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_CMB_DEFAULT = 'data_out/cmb.fits'
    character (len=S2_STRING_LEN) :: output_file_wnoise
    character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_WNOISE_DEFAULT = 'data_out/wnoise.fits'
    character (len=S2_STRING_LEN) :: output_file_wnoise_nobs
    character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_WNOISE_NOBS_DEFAULT = 'data_out/wnoise_nobs.fits'
    character (len=S2_STRING_LEN) :: output_file_wnoise_std
    character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_WNOISE_STD_DEFAULT = 'data_out/wnoise_std.fits'
    character (len=S2_STRING_LEN) :: output_file_param
    character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_PARAM_DEFAULT = 'data_out/param_out.par'
    real (kind=s2_sp) :: comb_gamma_lower
    real (kind=s2_sp) :: comb_gamma_upper
    integer :: iopt
    integer :: iobj
    integer :: fail = 0
    real (kind=s2_sp), allocatable, dimension (:) :: amplitude
    real (kind=s2_sp), allocatable, dimension (:) :: dilation
    real (kind=s2_sp), allocatable, dimension (:) :: alpha
    real (kind=s2_sp), allocatable, dimension (:) :: beta
    real (kind=s2_sp), allocatable, dimension (:) :: gamma
    type (comb_obj) :: obj_mother
    type (comb_csky) :: csky
    type (s2_pl) :: beam
    type (s2_cmb) :: cmb
    type (s2_wnoise) :: wnoise

        ! Interfaces
    interface comb_csim_param_check_type
    interface comb_csim_snap_to_grid

        ! Subroutines and functions
    subroutine comb_csim_print_header ()
    subroutine comb_csim_param_read ()
    subroutine comb_csim_param_write ()
    subroutine comb_csim_param_gen ()

end program comb_csim

Generates simulations of compact objects embedded in primordial cmb (convolved with beam if desired) and noise.

Author: J. D. McEwen (mcewen[AT]mrao.cam.ac.uk)

Version: 0.1 September 2004

Description of Variables

CODE_NAME

character (len=*), parameter :: CODE_NAME = 'COMB_CSIM'

CODE_DESCRIPTION

character (len=*), parameter :: CODE_DESCRIPTION = 'Compact object embedded CMB simulation'

CODE_VERSION

character (len=*), parameter :: CODE_VERSION = '1.0'

CODE_AUTHOR

character (len=*), parameter :: CODE_AUTHOR = 'Jason McEwen'

filename_param

character (len=S2_STRING_LEN) :: filename_param

description

character (len=S2_STRING_LEN) :: description

line

character (len=S2_STRING_LEN) :: line

handle

type (paramfile_handle) :: handle

NUM_COMMENT_LINES_CLT_FILE

integer, parameter :: NUM_COMMENT_LINES_CLT_FILE = 29

NUM_COMMENT_LINES_BEAM_FILE

integer, parameter :: NUM_COMMENT_LINES_BEAM_FILE = 8

TOL

real (kind=s2_sp) :: TOL = 1e-5

nside

integer :: nside

lmax

integer :: lmax

mmax

integer :: mmax

lmin

integer, parameter :: lmin = 2

NSIDE_DEFAULT

integer, parameter :: NSIDE_DEFAULT = 64

LMAX_DEFAULT

integer, parameter :: LMAX_DEFAULT = 128

MMAX_DEFAULT

integer, parameter :: MMAX_DEFAULT = 128

comb_type

character (len=S2_STRING_LEN) :: comb_type

COMB_TYPE_OPT_NUM

integer, parameter :: COMB_TYPE_OPT_NUM = 6

COMB_TYPE_OPT

character (len=S2_STRING_LEN), parameter, dimension (COMB_TYPE_OPT_NUM) :: COMB_TYPE_OPT = (/ 'butterfly   ', 'gaussian    ', 'mexhat      ', 'morlet      ', 'point       ', 'cos_thetaon2' /)

COMB_TYPE_DEFAULT

character (len=S2_STRING_LEN), parameter :: COMB_TYPE_DEFAULT = COMB_TYPE_OPT (1)

comb_num

integer :: comb_num

COMB_NUM_OPT

integer, parameter, dimension (2) :: COMB_NUM_OPT = (/ 0, 50 /)

COMB_NUM_DEFAULT

integer, parameter :: COMB_NUM_DEFAULT = 10

comb_seed

integer :: comb_seed

comb_amp_lower

real (kind=s2_sp) :: comb_amp_lower

COMB_AMP_LOWER_OPT

real (kind=s2_sp), parameter, dimension (2) :: COMB_AMP_LOWER_OPT = (/ 0.01e0, 5.0e0 /)

COMB_AMP_LOWER_DEFAULT

real (kind=s2_sp), parameter :: COMB_AMP_LOWER_DEFAULT = 0.1e0

comb_amp_upper

real (kind=s2_sp) :: comb_amp_upper

COMB_AMP_UPPER_OPT

real (kind=s2_sp), parameter, dimension (2) :: COMB_AMP_UPPER_OPT = (/ 0.01e0, 5.0e0 /)

COMB_AMP_UPPER_DEFAULT

real (kind=s2_sp), parameter :: COMB_AMP_UPPER_DEFAULT = 2.0e0

comb_dil_lower

real (kind=s2_sp) :: comb_dil_lower

COMB_DIL_LOWER_OPT

real (kind=s2_sp), parameter, dimension (2) :: COMB_DIL_LOWER_OPT = (/ 0.01e0, 1.0e0 /)

COMB_DIL_LOWER_DEFAULT

real (kind=s2_sp), parameter :: COMB_DIL_LOWER_DEFAULT = 0.05e0

comb_dil_upper

real (kind=s2_sp) :: comb_dil_upper

COMB_DIL_UPPER_OPT

real (kind=s2_sp), parameter, dimension (2) :: COMB_DIL_UPPER_OPT = (/ 0.01e0, 1.0e0 /)

COMB_DIL_UPPER_DEFAULT

real (kind=s2_sp), parameter :: COMB_DIL_UPPER_DEFAULT = 0.1e0

comb_dil_grid_num

integer :: comb_dil_grid_num

COMB_DIL_GRID_NUM_DEFAULT

integer, parameter :: COMB_DIL_GRID_NUM_DEFAULT = 5

comb_gamma_grid_num

integer :: comb_gamma_grid_num

COMB_GAMMA_GRID_NUM_DEFAULT

integer, parameter :: COMB_GAMMA_GRID_NUM_DEFAULT = 5

cmb_status

logical :: cmb_status

CMB_STATUS_DEFAULT

logical, parameter :: CMB_STATUS_DEFAULT = .true.

cmb_seed

integer :: cmb_seed

cmb_cl_file

character (len=S2_STRING_LEN) :: cmb_cl_file

CMB_CL_FILE_DEFAULT

character (len=S2_STRING_LEN), parameter :: CMB_CL_FILE_DEFAULT = 'data_in/wmap_lcdm_pl_model_yr1_v1.txt'

beam_status

logical :: beam_status

BEAM_STATUS_DEFAULT

logical, parameter :: BEAM_STATUS_DEFAULT = .false.

beam_gaussian

logical :: beam_gaussian

BEAM_GAUSSIAN_DEFAULT

logical, parameter :: BEAM_GAUSSIAN_DEFAULT = .true.

beam_fwhm

real (kind=s2_sp) :: beam_fwhm

BEAM_FWHM_DEFAULT

real (kind=s2_sp), parameter :: BEAM_FWHM_DEFAULT = 10.0e0

beam_cl_file

character (len=S2_STRING_LEN) :: beam_cl_file

BEAM_CL_FILE_DEFAULT

character (len=S2_STRING_LEN), parameter :: BEAM_CL_FILE_DEFAULT = 'data_in/map_w1_ampl_bl_yr1_v1.txt'

noise_status

logical :: noise_status

NOISE_STATUS_DEFAULT

logical, parameter :: NOISE_STATUS_DEFAULT = .false.

noise_full_sky

logical :: noise_full_sky

NOISE_FULL_SKY_DEFAULT

logical, parameter :: NOISE_FULL_SKY_DEFAULT = .false.

noise_seed

integer :: noise_seed

noise_std

real (kind=s2_sp) :: noise_std

NOISE_STD_DEFAULT

real (kind=s2_sp), parameter :: NOISE_STD_DEFAULT = 0.1e0

noise_sigma0

real (kind=s2_sp) :: noise_sigma0

NOISE_SIGMA0_DEFAULT

real (kind=s2_sp), parameter :: NOISE_SIGMA0_DEFAULT = 1.0e0

noise_nobs_file

character (len=S2_STRING_LEN) :: noise_nobs_file

NOISE_NOBS_FILE_DEFAULT

character (len=S2_STRING_LEN), parameter :: NOISE_NOBS_FILE_DEFAULT = 'data_in/wmap_w1_cleanimap_yr1_v1.fits'

noise_nobs_extension

integer :: noise_nobs_extension

NOISE_NOBS_EXTENSION_DEFAULT

integer, parameter :: NOISE_NOBS_EXTENSION_DEFAULT = 2

output_file_csky

character (len=S2_STRING_LEN) :: output_file_csky

OUTPUT_FILE_CSKY_DEFAULT

character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_CSKY_DEFAULT = 'data_out/csky.fits'

output_file_obj

character (len=S2_STRING_LEN) :: output_file_obj

OUTPUT_FILE_OBJ_DEFAULT

character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_OBJ_DEFAULT = 'data_out/obj.fits'

output_file_obj_mother

character (len=S2_STRING_LEN) :: output_file_obj_mother

OUTPUT_FILE_OBJ_MOTHER_DEFAULT

character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_OBJ_MOTHER_DEFAULT = 'data_out/obj_mother.fits'

output_file_cmb

character (len=S2_STRING_LEN) :: output_file_cmb

OUTPUT_FILE_CMB_DEFAULT

character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_CMB_DEFAULT = 'data_out/cmb.fits'

output_file_wnoise

character (len=S2_STRING_LEN) :: output_file_wnoise

OUTPUT_FILE_WNOISE_DEFAULT

character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_WNOISE_DEFAULT = 'data_out/wnoise.fits'

output_file_wnoise_nobs

character (len=S2_STRING_LEN) :: output_file_wnoise_nobs

OUTPUT_FILE_WNOISE_NOBS_DEFAULT

character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_WNOISE_NOBS_DEFAULT = 'data_out/wnoise_nobs.fits'

output_file_wnoise_std

character (len=S2_STRING_LEN) :: output_file_wnoise_std

OUTPUT_FILE_WNOISE_STD_DEFAULT

character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_WNOISE_STD_DEFAULT = 'data_out/wnoise_std.fits'

output_file_param

character (len=S2_STRING_LEN) :: output_file_param

OUTPUT_FILE_PARAM_DEFAULT

character (len=S2_STRING_LEN), parameter :: OUTPUT_FILE_PARAM_DEFAULT = 'data_out/param_out.par'

comb_gamma_lower

real (kind=s2_sp) :: comb_gamma_lower

comb_gamma_upper

real (kind=s2_sp) :: comb_gamma_upper

iopt

integer :: iopt

iobj

integer :: iobj

fail

integer :: fail = 0

amplitude

real (kind=s2_sp), allocatable, dimension (:) :: amplitude

dilation

real (kind=s2_sp), allocatable, dimension (:) :: dilation

alpha

real (kind=s2_sp), allocatable, dimension (:) :: alpha

beta

real (kind=s2_sp), allocatable, dimension (:) :: beta

gamma

real (kind=s2_sp), allocatable, dimension (:) :: gamma

obj_mother

type (comb_obj) :: obj_mother

csky

type (comb_csky) :: csky

beam

type (s2_pl) :: beam

cmb

type (s2_cmb) :: cmb

wnoise

type (s2_wnoise) :: wnoise

Description of Interfaces

comb_csim_param_check_type

interface comb_csim_param_check_type
    function comb_csim_param_check_type_char (param, types) result (pass)
        character (len=*), intent(in) :: param
        character (len=*), intent(in), dimension (:) :: types
        logical :: pass
    end function comb_csim_param_check_type_char
    function comb_csim_param_check_type_int (param, range) result (pass)
        integer, intent(in) :: param
        integer, intent(in), dimension (2) :: range
        logical :: pass
    end function comb_csim_param_check_type_int
    function comb_csim_param_check_type_real (param, range) result (pass)
        real, intent(in) :: param
        real, intent(in), dimension (2) :: range
        logical :: pass
    end function comb_csim_param_check_type_real
end interface comb_csim_param_check_type

comb_csim_snap_to_grid

interface comb_csim_snap_to_grid
    subroutine comb_csim_snap_to_grid (vals, lower, upper, num)
        real (kind=s2_sp), intent(inout), dimension (:) :: vals
        real (kind=s2_sp), intent(in) :: lower
        real (kind=s2_sp), intent(in) :: upper
        integer, intent(in) :: num
    end subroutine comb_csim_snap_to_grid
end interface comb_csim_snap_to_grid

Description of Subroutines and Functions

comb_csim_print_header

subroutine comb_csim_print_header ()
end subroutine comb_csim_print_header

Print comb_csim header to standard out.

Variables:

Has local scope so accesses variables defined in comb_csim program.

Author: J. D. McEwen

Version: 0.1 September 2004

comb_csim_param_read

subroutine comb_csim_param_read ()
    ! Calls: comb_error, getArgument, system_clock
end subroutine comb_csim_param_read

Parse comb_csim input parameters. Either read from input parameter file or read interactively from standard input.

Variables:

Has local scope so accesses variables defined in comb_csim program.

Author: J. D. McEwen

Version: 0.1 September 2004

comb_csim_param_write

subroutine comb_csim_param_write ()
end subroutine comb_csim_param_write

Write parameter values to screen.

Variables:

Has local scope so accesses variables defined in comb_csim program.

Author: J. D. McEwen

Version: 0.1 September 2004

comb_csim_param_gen

subroutine comb_csim_param_gen ()
    ! Calls: comb_csim_snap_to_grid
end subroutine comb_csim_param_gen

Generate samples compact object parameters in specified range.

Notes:

For now gamma angle and dilation are snapped to an equispaced grid so that corresponds to angles considered in a usual spherical wavelet analysis.

Variables:

Has local scope so accesses variables defined in comb_csim program.

Author: J. D. McEwen

Version: 0.1 September 2004