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!----------------------------------------------------------------------------- 28 29!****************************************************** 30! Sarah Smith (sjm84@mrao.cam.ac.uk) 31! and Graca Rocha (graca@mrao.cam.ac.uk) 32! June 2004 33! Module sub_ngpdf_sho based on code by Michael 34! Hobson and Graca Rocha 35! Module sub_ngpdf_powergauss based on code by 36! Anthony Challinor 37! Please include an appropriate acknowledgement 38! in any publications based on work that has 39! made use of this package - 'NGsims' 40!****************************************************** 41Program sky_ng_sim 42 ! Makes a full-sky non-Gaussian simulation using 43 ! values drawn from pdf made from eigenstates of 44 ! a simple harmonic oscillator or powers of a Gaussian 45 ! Uses Healpix pixelisation and subroutines 46 ! Some of the program is based on the Healpix synfast program 47 48 USE healpix_types 49 USE alm_tools, ONLY : map2alm_iterative, alm2map, alm2map_der, pow2alm_units 50 USE fitstools, ONLY : read_asctab, write_bintab, dump_alms 51 USE pix_tools, ONLY : nside2npix 52 USE head_fits, ONLY : add_card, merge_headers, get_card, write_minimal_header, del_card 53! USE utilities, ONLY : die_alloc 54 use misc_utils, only : wall_clock_time, assert_alloc, brag_openmp, fatal_error 55 USE extension, ONLY : getArgument, nArguments 56 USE paramfile_io, ONLY : paramfile_handle, parse_int, parse_init, & 57 & parse_string, parse_lgt, parse_double, parse_check_unused, & 58 & concatnl, scan_directories, parse_summarize, parse_finish, & 59 & get_healpix_data_dir, get_healpix_test_dir,get_healpix_pixel_window_file 60 61 USE sub_ngpdf_sho 62 USE sub_ngpdf_powergauss 63 USE sky_sub 64 65 Implicit none 66 67 INTEGER(I4B) :: nsmax ! The value of N_{side} in the pixelisation scheme 68 INTEGER(I4B) :: nlmax, lmax ! The maximum l-value used 69 INTEGER(I4B) :: nmmax ! The maximum m-value used (set equal to nlmax) 70 71 COMPLEX(SPC), DIMENSION(:,:,:), ALLOCATABLE :: alm_T !The a_{lm} values 72 REAL(SP), DIMENSION(:, :), ALLOCATABLE :: map_T !The pixel values in real space 73 REAL(SP), DIMENSION(:,:), ALLOCATABLE :: cl_T !The Cl values 74 REAL(KIND=DP), DIMENSION(:,:), ALLOCATABLE :: w8ring_T 75 REAL(SP), DIMENSION(:,:), ALLOCATABLE :: tmp_2d 76 77! INTEGER(I4B), DIMENSION(8,2) :: values_time 78! REAL(SP) :: clock_time 79 real(SP) :: time0, time1, clock_time, ptime0, ptime1, ptime 80 INTEGER(I4B) :: status 81 82 INTEGER(I4B) npixtot !The total number of pixels 83 CHARACTER(LEN=filenamelen) :: parafile = '' 84 CHARACTER(LEN=filenamelen) :: infile 85 character(len=FILENAMELEN) :: outfile_alms 86 CHARACTER(LEN=filenamelen) :: outfile 87 CHARACTER(LEN=filenamelen) :: windowfile 88 CHARACTER(LEN=filenamelen) :: windowname 89 CHARACTER(LEN=filenamelen) :: def_dir, def_file 90 CHARACTER(LEN=filenamelen) :: usr_dir, usr_file 91 CHARACTER(LEN=filenamelen) :: final_file 92 CHARACTER(LEN=filenamelen) :: healpixtestdir 93 Character(LEN=filenamelen) :: beam_file 94 character(len=filenamelen) :: description 95 character(len=100) :: chline 96 LOGICAL(LGT) :: ok, fitscl, polarisation = .FALSE. 97 logical(LGT) :: do_map, output_alms 98 CHARACTER(LEN=80), DIMENSION(1:180) :: header, header_PS !, header_file 99 CHARACTER(LEN=*), PARAMETER :: code = "sky_ng_sim" 100 character(len=*), parameter :: VERSION = HEALPIX_VERSION 101 character(len=80), dimension(1:1) :: units_power, units_map 102 CHARACTER(LEN=20) :: string_quad 103! CHARACTER(LEN=80) :: temptype, com_tt 104! REAL(SP) :: quadrupole 105! INTEGER(I4B) :: count_tt 106 INTEGER(I4B) nlheader 107 Integer :: i,m, n_pols, n_maps, simul_type, deriv 108 Real(DP) :: sigma0, factor 109 Real(DP),dimension(1:3) :: nu !Added to match with f90 version of shodev_driver 110 !nu(i) is the ith moment of the distribution. 111 !The size of nu is fixed to 3 to prevent problems with passing unallocated arrays 112 Real(DP) :: power, rms_alm, mean 113 Integer :: count, iter_order 114 Integer(I4B) :: pdf_choice !Used to chose type of pdf required 115 logical :: plot 116 ! External shodev_driver 117 Real(SP) :: Tmin, Tmax, convert 118 Integer, parameter :: npts = 1000 119 Real(SP) :: xline(npts), yline(npts), xstep, xval 120 type(paramfile_handle) :: handle 121 Real(DP) :: fwhm_arcmin 122 Real(SP) :: fwhm_deg 123 Character(len=30) :: xlabel 124 !----------------------------------------------------------------------- 125 ! get input parameters and create arrays 126 !----------------------------------------------------------------------- 127 128 call wall_clock_time(time0) 129 call cpu_time(ptime0) 130 ! --- read parameters interactively if no command-line arguments 131 ! --- are given, otherwise interpret first command-line argument 132 ! --- as parameter file name and read parameters from it: 133 if (nArguments() == 0) then 134 parafile='' 135 else 136 if (nArguments() /= 1) then 137 print *, " Usage: "//code//" [parameter file name]" 138 call fatal_error() 139 end if 140 call getArgument(1,parafile) 141 end if 142 handle = parse_init(parafile) 143 144 PRINT *, " " 145 PRINT *," "//code//" "//version 146 Write(*, '(a)') "*** Simulation of a non-Gaussian full-sky temperature map ***" 147 148 ! --- choose temp. only or temp. + deriv. --- 149 description = concatnl( & 150 & " Do you want to simulate", & 151 & " 1) Temperature only", & 152 !& " 2) Temperature + polarisation", & 153 & " 3) Temperature + 1st derivatives", & 154 & " 4) Temperature + 1st & 2nd derivatives") 155 simul_type = parse_int(handle, 'simul_type', default=1, valid=(/1,3,4/), descr=description) 156 deriv = 0 ! no derivatives 157 if (simul_type == 3) deriv = 1 158 if (simul_type == 4) deriv = 2 159 n_pols = 1 160 n_maps = max(1, 3*deriv) 161 162 ! --- gets the effective resolution of the sky map --- 1633 continue 164 description = concatnl( & 165 & "", & 166 & " Enter the resolution parameter (Nside) for the simulated skymap: ",& 167 & " (Npix = 12*Nside**2, where Nside HAS to be a power of 2, eg: 32, 512, ...)" ) 168 nsmax = parse_int(handle, 'nsmax', default=32, descr=description) 169 if (nside2npix(nsmax) < 0) then 170 print *, " Error: nsmax is not a power of two." 171 if (handle%interactive) goto 3 172 call fatal_error(code) 173 endif 174 175 ! --- gets the L range for the simulation --- 176 WRITE(chline,"(a,i5,a)") "We recommend: (0 <= l <= l_max <= ",3*nsmax-1,")" 177 description = concatnl(& 178 & "", & 179 & " Enter the maximum l range (l_max) for the simulation. ", & 180 & chline ) 181 nlmax = parse_int(handle, 'nlmax', default=2*nsmax, descr=description) 182 183 184 ! --- get filename for input power spectrum --- 185 chline = '' 186 healpixtestdir = get_healpix_test_dir() 187 if (trim(healpixtestdir)/='') chline = trim(healpixtestdir)//'/cl.fits' 188 description = concatnl( "", & 189 & " Enter input Power Spectrum filename", & 190 & " Can be either in FITS format or in the .dat format produced by CAMB.") 191 infile = parse_string(handle, 'infile', default=chline, descr=description, filestatus='old') 192 ! Find out whether input cl file is a fits file 193 fitscl = (index(infile, '.fits') /= 0) 194 195 ! --- gets the output sky map filename --- 196 description = concatnl(& 197 & "", & 198 & " Enter Output map file name (eg, test.fits) :", & 199 & " (or !test.fits to overwrite an existing file)" ) 200 outfile = parse_string(handle, "outfile", & 201 default="!test.fits", descr=description, filestatus="new") 202 do_map = (trim(outfile) /= '') 203 204 ! --- gets the output alm-filename --- 205 description = concatnl(& 206 & "", & 207 & " Enter file name in which to write a_lm used to synthesize the map : ", & 208 & " (eg alm.fits or !alms.fits to overwrite an existing file): ", & 209 & " (If '', the alms are not written to a file) ") 210 outfile_alms = parse_string(handle, 'outfile_alms', & 211 & default="''", descr=description, filestatus='new') 212 output_alms = (trim(outfile_alms) /= '') 213 214 ! --- gets the fwhm of beam --- 215 description = concatnl(& 216 & "", & 217 & " Enter FWMH of beam in arcminutes (0 for no beam):") 218 fwhm_arcmin = parse_double(handle, "fwhm_arcmin", & 219 default=0d0, descr=description, vmin = 0d0) 220 221 description = concatnl(& 222 & "", & 223 & " Enter an external file name containing ", & 224 & " a symmetric beam Legendre transform (eg: mybeam.fits)", & 225 & " NB: if set to an existing file, it will override the FWHM chosen above", & 226 & " if set to '', the gaussian FWHM will be used.") 227 beam_file = parse_string(handle, 'beam_file', default="''", & 228 & descr=description, filestatus='old') 229 if (beam_file /= '') then 230 fwhm_arcmin = 0. 231 print*,'fwhm_arcmin is now : 0.' 232 print*,'The beam file '//trim(beam_file)//' will be used instead.' 233 endif 234 235 ! including pixel window function, EH-2008-03-05 236 ! --- check for pixel-window-files --- 237 windowname = get_healpix_pixel_window_file(nsmax) 238 239 def_file = trim(windowname) 240 def_dir = get_healpix_data_dir() 241 24222 continue 243 final_file = '' 244 ok = .false. 245 ! if interactive, try default name in default directories first 246 if (handle%interactive) ok = scan_directories(def_dir, def_file, final_file) 247 if (.not. ok) then 248 ! not found, ask the user 249 description = concatnl("",& 250 & " Could not find window file", & 251 & " Enter the directory where this file can be found:") 252 usr_dir = parse_string(handle,'winfiledir',default='',descr=description) 253 if (trim(usr_dir) == '') usr_dir = trim(def_dir) 254 description = concatnl("",& 255 & " Enter the name of the window file:") 256 usr_file = parse_string(handle,'windowfile',default=def_file,descr=description) 257 ! look for new name in user provided or default directories 258 ok = scan_directories(usr_dir, usr_file, final_file) 259 ! if still fails, crash or ask again if interactive 260 if (.not. ok) then 261 print*,' File not found' 262 if (handle%interactive) goto 22 263 call fatal_error(code) 264 endif 265 endif 266 windowfile = final_file 267 268 PRINT *," " 269 270 !----------------------------------------------------------------------- 271 272 nmmax = nlmax 273 npixtot = nside2npix(nsmax) 274 275 !----------------------------------------------------------------------- 276 ! allocates space for arrays 277 !----------------------------------------------------------------------- 278 279 ! ALLOCATE(units_alm(1:1),units_map(1:1+2*polar),stat = status) 280 ! call assert_alloc(status, code,"units_alm & units_map") 281 282 ALLOCATE(alm_T(1:1,0:nlmax, 0:nmmax),stat = status) 283 call assert_alloc(status, code,"alm_T") 284 285 ALLOCATE(map_T(0:npixtot-1, 1:n_maps),stat = status) 286 call assert_alloc(status, code,"map_T") 287 288 ALLOCATE(w8ring_T(1:2*nsmax,1:1),stat = status) 289 call assert_alloc(status, code,"w8ring_T") 290 291 ! For now, not using ring weights for quadrature correction 292 w8ring_T = 1.d0 293 294 ! single analysis 295 iter_order = 0 296 297 ALLOCATE(cl_T(0:nlmax,1:1),stat = status) 298 call assert_alloc(status, code,"cl_T") 299 300 !------------------------------------------------------------------------ 301 ! Read in the input power spectrum 302 !----------------------------------------------------------------------- 303 304 cl_T = 0.0 305 !New lines added 8th June 2004 306 lmax = nlmax 307 call read_powerspec(infile, nsmax, lmax, cl_T, header_PS, fwhm_arcmin, units_power, beam_file = beam_file, winfile = windowfile) 308 call pow2alm_units(units_power, units_map) 309 call del_card(header_PS, (/ "TUNIT#","TTYPE#"/)) ! remove TUNIT* and TTYPE* from header to avoid confusion later on 310 311 !------------------------------------------------------------------------ 312 ! Draw pixel values in real space from non-Gaussian distribution 313 !------------------------------------------------------------------------ 314 315 Write (*,*) "Creating non-Gaussian map with flat power spectrum" 316 317 ! --- Chose the type of pdf to use --- 318 description = concatnl(& 319 & "", & 320 & " Select non-Gaussian pdf to use: Simple harmonic oscillator (1)", & 321 & "or powers of a Gaussian (2)" ) 322 pdf_choice = parse_int(handle, 'pdf_choice', default=1, vmin = 1, vmax = 2, descr=description) 323 324 !Write (*,*) "Use SHO pdf (0) or powers of a Gaussian pdf (1) ?" 325 !Read (*,*) m 326 If (pdf_choice .eq. 1) Then 327 call shodev_driver(handle, npixtot, sigma0, map_T(:,1), nu) 328 Else 329 ! call powergauss_driver(npixtot, npixtot, sigma0, map_T, nu) 330 call powergauss_driver(handle, npixtot, sigma0, map_T(:,1), nu) 331 End If 332 print*,'NG map MIN and MAX:',minval(map_T(:,1)),maxval(map_T(:,1)) 333 !Normalise the map 334 If (nu(1) .ge. 0) Then !nu (1) is set to -1 if highest non-zero alpha > 3 335 ! ie if we can't calculate moments analytically 336 map_T = map_T / sqrt(nu(2)) 337 338 !Compute the rms value of the pixels 339 power = 0 340 Do i = 0, npixtot-1 341 power = power + map_T(i,1)**2 342 End Do 343 power = Sqrt(power/npixtot) 344 Write (*,*) "rms value of the pixels is ", power 345 Else 346 Write (*,*) "Can't calculate theoretical value of nu(2) to normalise" 347 Write (*,*) "Normalising using rms pixel value instead" 348 Write (*,*) "note - this method may change the statistical properties of the map" 349 !Compute the mean value of the pixels 350 mean = SUM(map_T(:,1)*1.d0)/npixtot 351 Write (*,*) 'Mean pixel value is ',mean,' - adjusting to zero' 352 map_T(:,1) = map_T(:,1) - mean 353 !Compute the rms value of the pixels 354 power = 0 355 Do i = 0, npixtot-1 356 power = power + map_T(i,1)**2 357 End Do 358 power = Sqrt(power/npixtot) 359 Write (*,*) "rms value of the pixels is ", power ," - setting to 1.0" 360 map_T(:,1) = map_T(:,1) / power 361 !Check 362 power = 0 363 Do i = 0, npixtot-1 364 power = power + map_T(i,1)**2 365 End Do 366 power = Sqrt(power/npixtot) 367 Write (*,*) "Now rms value of the pixels is ", power 368 End If 369 !------------------------------------------------------------------------ 370 ! Now compute the a_{lm} values from this distribution 371 !------------------------------------------------------------------------ 372 373 Write (*,*) "Computing the values of the a_{lm}" 374! call map2alm(nsmax, nlmax, nmmax, map_T, alm_T, -1000.d0, w8ring_T) 375! call map2alm(nsmax, nlmax, nmmax, map_T, alm_T, (/ 0.d0, 0.d0 /), w8ring_T) 376 call map2alm_iterative(nsmax, nlmax, nmmax, iter_order, map_T(:,1:1), alm_T, w8ring=w8ring_T) 377 378 379 !Compute the rms value of the a_{lm} 380 power = 0.0 381 count = 0 382 Do i = 1, nlmax 383 power = power + abs(alm_T(1,i,0))**2. 384 count = count + 1 385 Do m = 1, i 386 power = power + 2.0*abs(alm_T(1,i,m))**2. 387 count = count + 2 388 End Do 389 End Do 390 rms_alm = sqrt(power/count) 391 Write (*,*) "rms value of the a_{lm} is ",rms_alm 392 factor = sqrt(FOURPI/npixtot) 393 Write (*,*) "Expected rms value is ", factor 394 395 !------------------------------------------------------------------------ 396 ! Multiply the a_{lm} by the correct power spectrum 397 !------------------------------------------------------------------------ 398 399 Do i = 1, nlmax 400 !Write (*,*) 'cl(',i,') is ', cl_T(i,1) 401 alm_T(1,i,0:i) = alm_T(1,i,0:i)*sqrt(cl_T(i,1))/factor 402 End Do 403 404 !------------------------------------------------------------------------ 405 ! Now transform back to a real space map 406 !------------------------------------------------------------------------ 407 408 select case (deriv) 409 case(0) ! temperature only 410 call alm2map( nsmax,nlmax,nmmax,alm_T,map_T(:,1)) 411 case(1) 412 call alm2map_der(nsmax,nlmax,nmmax,alm_T,map_T(:,1),map_T(:,2:3)) 413 case(2) 414 call alm2map_der(nsmax,nlmax,nmmax,alm_T,map_T(:,1),map_T(:,2:3),map_T(:,4:6)) 415 case default 416 print*,'Not valid case' 417 print*,'Derivatives: ', deriv 418 call fatal_error(code) 419 end select 420 421 !------------------------------------------------------------------------ 422 ! Plot histogram of pixel distribution (if required) 423 !------------------------------------------------------------------------ 424 425#ifdef PGPLOT 426 description = concatnl( "", & 427 & " Plot histogram of pixel values? (True or False)") 428 plot = parse_lgt(handle, 'plot', default=.false., descr=description) 429! Write (*,*) 'Plot histrogram of pixel values (0) or not (1)?' 430! Read (*,*) m 431 If (plot) Then 432 Tmin = 0.0 433 Tmax = 0.0 434 do i=1,npixtot 435 if (map_T(i,1) .lt. Tmin) Tmin=map_T(i,1) 436 if (map_T(i,1) .gt. Tmax) Tmax=map_T(i,1) 437 end do 438 call pgbegin(0,'?',1,1) 439 if (.not. fitscl) then 440 ! If reading in from .dat file Cls are already converted to uK 441 convert = 1.0 442 xlabel = 'Temperature / \gmK' 443 else if (ABS(cl_T(2,1) - 1.0) .lt. 0.05) then 444 ! If values in fits file are normalised to Cl(2) = 1.0 445 ! convert to uK for plot, assuming COBE normalisation 446 ! of Qrms = 18uK 447 convert = sqrt(4.0*PI/5.0)*18.0 448 xlabel = 'Temperature / \gmK' 449 else 450 ! Otherwise, units unknown 451 convert = 1.0 452 xlabel = 'Pixel value' 453 end if 454 Tmin = Tmin * convert 455 Tmax = Tmax * convert 456 call pghist(npixtot,map_T*convert,Tmin,Tmax,200,0) 457 !Pixel values scaled to convert to uK if input fits file of cl values 458 !which are normalised to C_2 = 1.0 459 call pglab(xlabel,'Number of pixels','') 460 !Calculate 2*mean square value of T (ie 2 sigma^2) 461 power = 0 462 Do i = 1, npixtot 463 power = power + map_T(i,1)**2 464 End Do 465 power = 2 * power*(convert**2) / npixtot 466 xstep=(Tmax-Tmin)/real(npts-1) 467 do i=1,npts 468 xval=Tmin+real(i)*xstep 469 xline(i)=xval 470 yline(i)=exp(-xval**2/power)/sqrt(PI*power) 471 yline(i)=yline(i)*real(npixtot)*(Tmax-Tmin)/200. 472 end do 473 call pgsci(13) 474 call pgline(npts, xline, yline) 475 call pgend 476 End If 477#endif 478 479 call parse_check_unused(handle, code=code) 480 call parse_summarize(handle, code=code) 481 call parse_finish(handle) 482 call brag_openmp() 483 484 !----------------------------------------------------------------------- 485 ! write the alm to FITS file 486 !----------------------------------------------------------------------- 487 if (output_alms) then 488 489 PRINT *," "//code//"> outputting alms " 490 491 do i = 1, n_pols 492 header = "" 493 ! put inherited information immediatly, so that keyword values can be updated later on 494 ! by current code values 495 call add_card(header,"COMMENT","****************************************************************") 496 call merge_headers(header_PS, header) ! insert header_PS in header at this point 497 call add_card(header,"COMMENT","****************************************************************") 498 499 ! start putting information relative to this code and run 500 call write_minimal_header(header, 'alm', append=.true., & 501 creator = CODE, version = VERSION, & 502 nlmax = nlmax, nmmax = nmmax, & 503 !randseed = ioriginseed, & 504 units = units_map(i), polar = polarisation) 505 if (i == 1) then 506 call add_card(header,"EXTNAME","'SIMULATED a_lms (TEMPERATURE)'", update = .true.) 507 elseif (i == 2) then 508 call add_card(header,"EXTNAME","'SIMULATED a_lms (GRAD / ELECTRIC component)'", update = .true.) 509 elseif (i == 3) then 510 call add_card(header,"EXTNAME","'SIMULATED a_lms (CURL / MAGNETIC component)'", update = .true.) 511 endif 512 !if (input_cl) then 513 call add_card(header,"HISTORY","File with input C(l):") 514 call add_card(header,"HISTORY",trim(infile)) 515 call add_card(header,"HISTORY","These alms are multiplied by pixel and beam window functions") 516 call add_card(header,"NSIDE" ,nsmax, "Resolution parameter for HEALPIX") 517 if (trim(beam_file) == '') then 518 call add_card(header,"FWHM" ,fwhm_deg ," [deg] FWHM of gaussian symmetric beam") 519 else 520 call add_card(header,"BEAM_LEG",trim(beam_file), & 521 & "File containing Legendre transform of symmetric beam") 522 endif 523 call add_card(header,"PDF_TYPE",pdf_choice,"1: Harmon. Oscill. ;2: Power of Gauss.") 524 !endif 525 call add_card(header,"HISTORY") 526 nlheader = SIZE(header) 527 call dump_alms (outfile_alms,alm_T(i,0:nlmax,0:nmmax),nlmax,header,nlheader,i-1_i4b) 528 enddo 529 endif 530 !----------------------------------------------------------------------- 531 ! write the map to FITS file 532 !----------------------------------------------------------------------- 533 534 nlheader = SIZE(header) 535 header = "" 536 PRINT *," "//code//"> Writing sky map to FITS file " 537 ! put inherited information immediatly, so that keyword values can be updated later on 538 ! by current code values 539 call add_card(header,"COMMENT","****************************************************************") 540 call merge_headers(header_PS, header) ! insert header_PS in header at this point 541 call add_card(header,"COMMENT","****************************************************************") 542 ! start putting information relative to this code and run 543 call write_minimal_header(header, 'map', append=.true., & 544 nside = nsmax, ordering = 'RING', & !, coordsys = coordsys, & 545 fwhm_degree = fwhm_arcmin / 60.d0, & 546 beam_leg = trim(beam_file), & 547 polar = polarisation, & 548 deriv = deriv, & 549 creator = CODE, version = VERSION, & 550 nlmax = lmax, & 551 !randseed = ioriginseed, & 552 units = units_map(1) ) 553 !Note using lmax rather than nlmax - lmax will be smaller than nlmax if not enough values 554 ! in input power spectrum file 555 556 call add_card(header,"PDF_TYPE",pdf_choice,"1: Harmon. Oscill. ;2: Power of Gauss.") 557 call add_card(header,"EXTNAME","'SIMULATED MAP'", update=.true.) 558 call add_card(header,"COMMENT","*************************************") 559 560! allocate(tmp_2d(0:npixtot-1,1:1)) 561! tmp_2d(:,1) = map_T 562! call write_bintab(tmp_2d, npixtot, 1, header, nlheader, outfile) 563! deallocate(tmp_2d) 564 call write_bintab(map_T, npixtot, n_maps, header, nlheader, outfile) 565!!$ endif 566 567 !----------------------------------------------------------------------- 568 ! deallocate memory for arrays 569 !----------------------------------------------------------------------- 570 571 DEALLOCATE( map_T ) 572 DEALLOCATE( alm_T ) 573 ! DEALLOCATE( units_alm, units_map ) 574 575 !----------------------------------------------------------------------- 576 ! output and report card 577 !----------------------------------------------------------------------- 578 call wall_clock_time(time1) 579 call cpu_time(ptime1) 580 clock_time = time1 - time0 581 ptime = ptime1 - ptime0 582 583 WRITE(*,9000) " " 584 WRITE(*,9000) " Report Card for "//code//" simulation run" 585 WRITE(*,9000) "----------------------------------------" 586 WRITE(*,9000) " " 587 WRITE(*,9000) " Input power spectrum : "//TRIM(infile) 588 WRITE(*,9010) " Multipole range : 0 < l <= ", nlmax 589 WRITE(*,9010) " Number of pixels : ", npixtot 590 ! WRITE(*,9020) " Pixel size in arcmin : ", pix_size_arcmin 591 ! WRITE(*,9010) " Initial random # seed: ", ioriginseed 592 if (trim(beam_file) == '') then 593 write(*,9020) " Gauss. FWHM in arcmin: ", fwhm_arcmin 594 else 595 write(*,9000) " Beam file: "//trim(beam_file) 596 endif 597 WRITE(*,9000) " Output map : "//TRIM(outfile) 598 write(*,9030) " Clock and CPU time [s] : ", clock_time, ptime 599 600 !----------------------------------------------------------------------- 601 ! end of routine 602 !----------------------------------------------------------------------- 603 604 WRITE(*,9000) " " 605 WRITE(*,9000) " "//code//"> normal completion" 606 6079000 format(a) 6089010 format(a,i16) 6099020 format(a,g20.5) 6109030 format(a,f11.2,f11.2) 611 612 613 Stop 614 615End Program sky_ng_sim 616