bianchi2_sim.f90
Go to the documentation of this file.
1 !------------------------------------------------------------------------------
2 ! bianchi2_sim -- BIANCHI2 simulation program
3 !
17 !------------------------------------------------------------------------------
18 
19 program bianchi2_sim
20 
21  use s2_types_mod
22  use s2_sky_mod, only: s2_sky_file_type_map, s2_sky_file_type_sky
23  use bianchi2_sky_mod
24  use bianchi2_lut_mod
25  use bianchi2_error_mod
26  use pix_tools, only: nside2npix
27 
28  use extension, only: getargument, narguments
29  use paramfile_io, only: paramfile_handle, parse_init, parse_int, &
30  parse_real, parse_double, parse_lgt, parse_string, concatnl
31 
32  implicit none
33 
34  character(len=S2_STRING_LEN) :: filename_param
35  character(len=S2_STRING_LEN) :: description
36  character(len=S2_STRING_LEN) :: line
37  type(paramfile_handle) :: handle
38 
39  real(s2_dp), parameter :: omega_matter_lower = 0d0
40  real(s2_dp), parameter :: omega_matter_upper = 1d0
41  real(s2_dp), parameter :: omega_matter_default = 0.3d0
42  real(s2_dp), parameter :: omega_lambda_lower = 0d0
43  real(s2_dp), parameter :: omega_lambda_upper = 1d0
44  real(s2_dp), parameter :: omega_lambda_default = 0.6d0
45  real(s2_dp), parameter :: h_lower = 0d0
46  real(s2_dp), parameter :: h_upper = 10d0
47  real(s2_dp), parameter :: h_default = 1d-2
48  real(s2_dp), parameter :: ze_lower = 1d2
49  real(s2_dp), parameter :: ze_upper = 1d4
50  real(s2_dp), parameter :: ze_default = 1d3
51  real(s2_dp), parameter :: wh_default = 1d-10
52  real(s2_sp), parameter :: fwhm_default = 330d0
53  logical, parameter :: rhand_default = .true.
54  logical, parameter :: harmonic_space_default = .true.
55  logical, parameter :: out_cl_default = .false.
56  logical, parameter :: rescale_cl_default = .true.
57  logical, parameter :: read_lut_default = .false.
58  integer, parameter :: ntheta_default = 100
59  integer, parameter :: lmax_default = 64
60  integer, parameter :: nside_default = 128
61  character(len=*), parameter :: file_type_map_str = 'map'
62  character(len=*), parameter :: file_type_sky_str = 'sky'
63 
64  character(len=S2_STRING_LEN) :: filename_out, filename_out_cl
65  character(len=S2_STRING_LEN) :: filetype_str = file_type_map_str
66  integer :: filetype = s2_sky_file_type_map
67 
68  type(bianchi2_sky) :: b
69  real(s2_sp), allocatable :: cl(:)
70  real(s2_dp) :: omega_matter, omega_lambda, h, ze, wh
71  integer :: nside, lmax
72  ! Default angles in degrees for user input, but converted to radians later.
73  real(s2_sp) :: alpha=0e0, beta=-90e0, gamma=0e0
74  logical :: rhand = .true.
75  logical :: apply_beam = .false.
76  real(s2_sp) :: fwhm = fwhm_default
77  logical :: harmonic_space
78  integer :: ntheta
79  logical :: rotation_alm = .false. ! Choose true only for
80  ! simulation in real space but
81  ! rotation in harmonic space.
82  logical :: out_cl
83  logical :: rescale_cl
84  logical :: read_lut
85  character(len=S2_STRING_LEN) :: filename_lut
86  type(bianchi2_lut) :: lut
87 
88  ! Set default parameter values.
89  filename_out = 'sky.fits'
90  filename_out_cl = 'out_Cl.txt'
91  filename_lut = 'lut.dat'
92  omega_matter = omega_matter_default
93  omega_lambda = omega_lambda_default
94 
95  h = h_default
96  ze = ze_default
97  wh = wh_default
98  nside = nside_default
99  lmax = lmax_default
100  rhand = rhand_default
101 
102  harmonic_space = harmonic_space_default
103  ntheta = ntheta_default
104 
105  write(*,'(a)') '***********************************************'
106  write(*,'(a)') 'BIANCHI2 VII_h rotating universe CMB simulation'
107  write(*,'(a)') 'Anthony Lasenby & Jason McEwen October 2005'
108  write(*,'(a)') 'mcewen@mrao.cam.ac.uk '
109  write(*,'(a)') '***********************************************'
110 
111 
112  !---------------------------------------
113  ! Parse parameters
114  !---------------------------------------
115 
116  ! Initialise file parser.
117  if(narguments() == 0) then
118  filename_param = ''
119  else
120  if(narguments() /= 1) then
121  call bianchi2_error(bianchi2_error_sim_narg, 'bianchi2_sim', &
122  comment_add='Usage: bianchi2_sim [input parameter filename]')
123  end if
124  call getargument(1, filename_param)
125  end if
126  handle = parse_init(trim(filename_param))
127 
128 99 continue
129 
130  ! Get omega_matter.
131  write(line,'(a,f4.1,a,f4.1,a)') '(In the range ', &
132  omega_matter_lower, ' <= omega_matter < ', omega_matter_upper, ')'
133  description = concatnl('', &
134  'Enter omega_matter: ', &
135  line)
136 1 continue
137  omega_matter = parse_double(handle, 'omega_matter', &
138  default=omega_matter_default, descr=description)
139  if(omega_matter < omega_matter_lower .or. &
140  omega_matter >= omega_matter_upper) then
141  if(handle%interactive) goto 1
142  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
143  comment_add='omega_matter invalid')
144  end if
145 
146  ! Get omega_lambda.
147  write(line,'(a,f4.1,a,f4.1,a)') '(In the range ', &
148  omega_lambda_lower, ' <= omega_lambda < ', omega_lambda_upper, ')'
149  description = concatnl('', &
150  'Enter omega_lambda: ', &
151  line)
152 2 continue
153  omega_lambda = parse_double(handle, 'omega_lambda', &
154  default=omega_lambda_default, descr=description)
155  if(omega_lambda < omega_lambda_lower .or. &
156  omega_lambda >= omega_lambda_upper) then
157  if(handle%interactive) goto 1
158  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
159  comment_add='omega_lambda invalid')
160  end if
161 
162  ! Check model closed.
163  if((omega_matter + omega_lambda) >= 1d0) then
164  if(handle%interactive) then
165  write(*,*)
166  write(*,'(a)') ' * Only open models are allowed *'
167  goto 99
168  end if
169  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
170  comment_add='Only open models are allowed')
171  end if
172 
173  ! Get h.
174  write(line,'(a,f4.1,a,f4.1,a)') '(In the range ', &
175  h_lower, ' <= h <= ', h_upper, ')'
176  description = concatnl('', &
177  'Enter h: ', &
178  line)
179 3 continue
180  h = parse_double(handle, 'h', &
181  default=h_default, descr=description)
182  if(h < h_lower .or. h > h_upper) then
183  if(handle%interactive) goto 3
184  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
185  comment_add='h invalid')
186 
187  end if
188 
189  ! Get zE.
190  write(line,'(a,e9.2,a,e9.2,a)') '(In the range ', &
191  ze_lower, ' <= zE <= ', ze_upper, ')'
192  description = concatnl('', &
193  'Enter zE: ', &
194  line)
195 4 continue
196  ze = parse_double(handle, 'zE', &
197  default=ze_default, descr=description)
198  if(ze < ze_lower .or. ze > ze_upper) then
199  if(handle%interactive) goto 4
200  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
201  comment_add='zE invalid')
202 
203  end if
204 
205  ! Get wH.
206  description = concatnl('', &
207  'Enter wH: ')
208 5 continue
209  wh = parse_double(handle, 'wH', &
210  default=wh_default, descr=description)
211  if(wh < 0d0) then
212  if(handle%interactive) goto 5
213  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
214  comment_add='wH invalid')
215  end if
216 
217  ! Get right-handedness status.
218  description = concatnl('', &
219  'Enter right-handedness status (logical): ')
220  rhand = parse_lgt(handle, 'rhand', &
221  default=rhand_default, descr=description)
222 
223  ! Get alpha.
224  description = concatnl('', &
225  'Enter alpha (degrees): ')
226 6 continue
227  alpha = parse_real(handle, 'alpha', &
228  default=alpha, descr=description)
229  if(alpha < -360d0 .or. alpha > 360) then
230  if(handle%interactive) goto 6
231  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
232  comment_add='alpha invalid')
233  end if
234  alpha = alpha / 180e0 * pi
235 
236  ! Get beta.
237  description = concatnl('', &
238  'Enter beta (degrees): ')
239 7 continue
240  beta = parse_real(handle, 'beta', &
241  default=beta, descr=description)
242  if(beta < -180d0 .or. beta > 180) then
243  if(handle%interactive) goto 7
244  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
245  comment_add='beta invalid')
246  end if
247  beta = beta / 180e0 * pi
248 
249  ! Get gamma.
250  description = concatnl('', &
251  'Enter gamma (degrees): ')
252 8 continue
253  gamma = parse_real(handle, 'gamma', &
254  default=gamma, descr=description)
255  if(gamma < -360d0 .or. gamma > 360) then
256  if(handle%interactive) goto 8
257  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
258  comment_add='gamma invalid')
259  end if
260  gamma = gamma / 180e0 * pi
261 
262  ! Get apply_beam.
263  description = concatnl('', &
264  'Enter apply_beam status (logical): ')
265  apply_beam = parse_lgt(handle, 'apply_beam', &
266  default=apply_beam, descr=description)
267 
268  if(apply_beam) then
269  ! Get beam fwhm.
270  description = concatnl('', &
271  'Enter beam fwhm (arcmin): ')
272 9 continue
273  fwhm = parse_real(handle, 'fwhm', &
274  default=fwhm_default, descr=description)
275  if(fwhm <= 0d0) then
276  if(handle%interactive) goto 9
277  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
278  comment_add='fwhm invalid')
279  end if
280  end if
281 
282  ! Get lmax.
283  description = concatnl("", &
284  "Enter the maximum harmonic l (lmax) for the simulated sky: ")
285 10 continue
286  lmax = parse_int(handle, 'lmax', &
287  default=lmax_default, descr=description)
288  if(lmax < 0) then
289  if(handle%interactive) goto 10
290  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
291  comment_add='lmax invalid')
292  endif
293 
294 
295  ! Get nside.
296  description = concatnl("", &
297  "Enter the resolution parameter (nside) for the simulated sky: ", &
298  "(npix = 12*nside**2, where nside must be a power of 2)")
299 11 continue
300  nside = parse_int(handle, 'nside', &
301  default=nside_default, descr=description)
302  if(nside2npix(nside) < 0) then
303  if(handle%interactive) goto 11
304  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
305  comment_add='nside invalid')
306  endif
307 
308  ! Get harmonic_space.
309  description = concatnl('', &
310  'Enter harmonic_space status (logical): ')
311  harmonic_space = parse_lgt(handle, 'harmonic_space', &
312  default=harmonic_space_default, descr=description)
313  if (harmonic_space == .true.) then
314  ! Get Ntheta.
315  description = concatnl("", &
316  "Enter the number of terms (Ntheta) used for the integrations: ")
317 12 continue
318  ntheta = parse_int(handle, 'Ntheta', &
319  default=ntheta_default, descr=description)
320  if (ntheta <9) then
321  if(handle%interactive) goto 12
322  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
323  comment_add='Ntheta invalid')
324  end if
325  end if
326 
327  ! Get filename_out.
328  description = concatnl('', &
329  'Enter filename_out: ')
330  filename_out = parse_string(handle, 'filename_out', &
331  default=trim(filename_out), descr=description)
332 
333  ! Get output file type: map or sky.
334  description = concatnl('', &
335  'Enter output file type (filetype={map; sky}): ')
336 13 continue
337  filetype_str = parse_string(handle, 'filetype', &
338  default=trim(filetype_str), descr=description)
339 
340  ! Set filetype integer status.
341  select case (filetype_str)
342 
343  case (file_type_map_str)
344  filetype = s2_sky_file_type_map
345 
346  case (file_type_sky_str)
347  filetype = s2_sky_file_type_sky
348 
349  case default
350  if(handle%interactive) goto 13
351  call bianchi2_error(bianchi2_error_sim_param_invalid, 'bianchi2_sim', &
352  comment_add='Invalid output file type')
353 
354  end select
355 
356  ! Get out_Cl.
357  description = concatnl('', &
358  'Enter out_Cl (logical): ')
359  out_cl = parse_lgt(handle, 'out_Cl', &
360  default=out_cl_default, descr=description)
361 
362  if (out_cl == .true.) then
363  ! Get rescale_Cl.
364  description = concatnl('', &
365  'Enter rescale_Cl (logical): ')
366  rescale_cl = parse_lgt(handle, 'rescale_Cl', &
367  default=rescale_cl_default, descr=description)
368 
369  ! Get filename_out_Cl.
370  description = concatnl('', &
371  'Enter filename_out_Cl: ')
372  filename_out_cl = parse_string(handle, 'filename_out_Cl', &
373  default=trim(filename_out_cl), descr=description)
374  end if
375 
376  ! Get read_LUT.
377  description = concatnl('', &
378  'Enter read_LUT (logical): ')
379  read_lut = parse_lgt(handle, 'read_LUT', &
380  default=read_lut_default, descr=description)
381 
382  if (read_lut==.true.) then
383  ! Get filename_LUT.
384  description = concatnl('', &
385  'Enter filename_LUT: ')
386  filename_lut = parse_string(handle, 'filename_LUT', &
387  default=trim(filename_lut), descr=description)
388  end if
389 
390  !---------------------------------------
391  ! Run simulation and save sky
392  !---------------------------------------
393 
394  ! Simulated bianchi2 sky.
395  write(*,'(a)')
396 
397 
398  ! Initialise bianchi2 object.
399  if (harmonic_space) then
400 
401  ! If using a Look-Up-Table.
402  if (read_lut==.true.) then
403 
404  lut = bianchi2_lut_read(lmax,ntheta,filename_lut)
405  write(*,'(a)') 'Computing BIANCHI2 simulation in harmonic space, using a lut...'
406  b = bianchi2_sky_init_alm(omega_matter, omega_lambda, h, ze, wh, rhand, &
407  nside, lmax, ntheta, alpha, beta, gamma,lut)
408  call bianchi2_lut_free(lut)
409  else
410  write(*,'(a)') 'Computing BIANCHI2 simulation in harmonic space, without using a lut...'
411  b = bianchi2_sky_init_alm(omega_matter, omega_lambda, h, ze, wh, rhand, &
412  nside, lmax, ntheta, alpha, beta, gamma)
413  end if
414 
415  call bianchi2_sky_compute_map(b,nside)
416 
417  write(*,'(a)') 'Simulation complete'
418  write(*,'(a)')
419 
420  else
421 
422  write(*,'(a)') 'Computing BIANCHI2 simulation in real space...'
423  b = bianchi2_sky_init(omega_matter, omega_lambda, h, ze, wh, rhand, &
424  nside)
425 
426  write(*,'(a)') 'Simulation complete'
427  write(*,'(a)')
428 
429  ! Perform rotation
430  call bianchi2_sky_rotate(b, alpha, beta, gamma, lmax, nside,rotation_alm)
431  ! Compute the alm (needed for Cl).
432  call bianchi2_sky_compute_alm(b,lmax,lmax)
433 
434  end if
435 
436  ! Apply beam if required.
437  if(apply_beam) then
438 
439  ! Recompute alms if necessary (if not necessary does nothing).
440  ! Note, if compute alms directly but then perform
441  ! bianchi2_sky_rotate, alms are removed and must be recomputed.
442  call bianchi2_sky_compute_alm(b, lmax, lmax)
443 
444  ! If map already present then apply beam will recompute the map.
445  call bianchi2_sky_apply_beam(b, fwhm, lmax)
446 
447  end if
448 
449  ! Save sky.
450  call bianchi2_sky_write(b, filename_out, filetype)
451  write(*,'(a,a)') 'Simulated map written to ', trim(filename_out)
452  write(*,'(a)')
453 
454  if (out_cl == .true.) then
455  ! Compute and write the power spectrum.
456  allocate(cl(0:lmax))
457  call bianchi2_sky_compute_cl(b,lmax,cl)
458  call bianchi2_sky_write_cl(cl,lmax,rescale_cl,filename_out_cl)
459  deallocate(cl)
460  write(*,'(a,a)') 'Simulated power spectrum written to ', &
461  trim(filename_out_cl)
462  write(*,'(a)')
463  end if
464 
465  ! Write bianchi2 simulation variables to standard output.
466  write(*,'(a)') 'Simulation parameters:'
467  call bianchi2_sky_param_write(b)
468 
469  ! Free bianchi2 object.
470  call bianchi2_sky_free(b)
471 
472 end program bianchi2_sim
Generated on Tue Jul 28 2015 20:36:54 by   doxygen 1.8.3.1