ccc gadget2psis_v1.1.f Time-stamp: <2016-06-07 16:30:08 hamana> ! ! This is based on gadget2psis_v1.0.f with some modifications made. ! This assumes that input Nbody data is in the form of snapshot files of Gadget2. ! !* required library: ! Healpix (gadget2psis_v1.1.f was developed and tested with Healpix_3.11) ! cfitsio ! !* compile with gfortran: ! gfortran gadget2psis_v1.1.f -o gadget2psis.exe -O3 -I/usr/local/Healpix_3.11/include -DGFORTRAN -fno-second-underscore -fopenmp -fPIC -L/usr/local/Healpix_3.11/lib -L/usr/local/lib -lhealpix -lhpxgif -lcfitsio -Wl,-R/usr/local/lib ! !!! usage: ! gadget2psis.exe [options] ! [options] are: ! -path : PATH to input Gadget2 snapshot files ! -snap : Gadget snapshot number (integer >= 0) ! -np : [number of particles of Nbody data]^{1/3} ! -nres : Healpix's parameter "nres" ! -d : path to output dir ! -f : output file name ! -th : thickness of shell [Mpc/h] ! -kp : lens-plane number (integer >= 1 ***closest is 1***) ! -wde : w_DE ! -nt : number of openmp threads ! !!! example: ! ./gadget2psis.exe -path DATA/450Mpc/snapdir_002 -snap 2 -np 2048 -nres 12 -d DATA -f b450r004p001 -th 150.0 -kp 1 -wde -1.0 -nt 8 ! ./gadget2psis.exe -path DATA/450Mpc/snapdir_001 -snap 1 -np 2048 -nres 12 -d DATA -f b450r004p002 -th 150.0 -kp 2 -wde -1.0 -nt 8 ! ./gadget2psis.exe -path DATA/450Mpc/snapdir_000 -snap 0 -np 2048 -nres 12 -d DATA -f b450r004p003 -th 150.0 -kp 3 -wde -1.0 -nt 8 ! !!! copyright --- Takashi Hamana ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MODULE cosmoparam real(8) omegam real(8) wde END MODULE cosmoparam !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! program gadget2psis use alm_tools use healpix_types use pix_tools use udgrade_nr use fitstools, only : output_map use head_fits use cosmoparam implicit real (a-h,o-z) implicit integer (i-n) !!! valiables for read_snapshot integer(4) :: SNAPSHOT ! number of dump character(len=256) :: path character(len=256) :: filename, dirname, snumber, fnumber integer(4) :: fn,i,nstart integer(4) :: j,k,nlow integer(4) :: nsnap real(4) :: pos(1:3) !!! header structure in c !struct io_header_1 !{ ! int npart[6]; ! double mass[6]; ! double time; ! double redshift; ! int flag_sfr; ! int flag_feedback; ! int npartTotal[6]; ! int flag_cooling; ! int num_files; ! double BoxSize; ! double Omega0; ! double OmegaLambda; ! double HubbleParam; ! char fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8]; /* fills to 256 Bytes */ !} header1; !*** transforming to Fortran integer(4) :: npart(0:5) real(8) :: massarr(0:5) real(8) :: a real(8) :: redshift integer(4) :: flag_sfr integer(4) :: flag_feedback integer(4) :: nall(0:5) integer(4) :: flag_cooling integer(4) :: FILES real(8) :: boxsize ! shoube Mpc/h unit ***check this*** ! real(8) :: omegam ! defined in MODULE cosmoparam real(8) :: OmegaLambda real(8) :: HubbleParam integer(4) :: unused(24) !!! variables for Healpix integer(i4b) :: nres,nside,lmax,hporder,hpspin integer(i8b) :: npix,ipring,ipnest,ipix,np real(dp) :: hpvector(3) real(sp),allocatable :: hpmap_phi(:),hpmap_alp(:,:) real(sp),allocatable :: hpmap_mag(:,:) real(sp),allocatable :: hpmap(:,:) character(len=80), dimension(1:64) :: mapheader real(dp), allocatable, dimension(:,:) :: dw8 real(dp), dimension(2) :: hpz complex(spc), allocatable, dimension(:,:,:) :: hpalm real(sp),allocatable :: hpcl(:,:) !!! variables for Ray-tracing real(sp),allocatable :: sdens(:) real(4) :: rmin2,rmax2,pixpart real(8) :: sdensnorm real(4) :: alens character(len=256) :: pnumber character(len=256) :: pfname character(len=256) :: fname,indir,outfile,outdir,fin,filein,outlog character*256 opt,arg include 'omp_lib.h' ! !!!!!! definitions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! chubble0=2.9979e3 ! c/H0 in Mpc/h (H0=100h km/sex/Mpc) !!! Healpix related hporder=1 ! integer =1 for ring, =2 for nested hpz(1)=-1.0_dp hpz(2)=1.0_dp ! !!!!!!! set up parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! narg=iargc() c write(*,*)narg !!! cosmological parameters *** FLAT COSMOLOGY ONLY *** wde=-1.0 !!! Gadget related path='/data1/DATA_HSC/SNAP/L450/004/snapdir_000' SNAPSHOT=0 np=1024 !!! Healpix related nres=12 !!! output file names outdir='/home/hamana/T2K_LENS/PSIS' fname='test' !!! RayTrace related thick=150.0 !Mpc/h kplane=2 !!! number of threds nopenmp=8 do i=1,narg call getarg(i,opt) call getarg(i+1,arg) select case (opt) case ('-path') path=arg case ('-snap') read(arg,*)SNAPSHOT ! gadget snapshot number case ('-np') read(arg,*)np ! [number of particles]^{1/3} = ll^3 case ('-nres') read(arg,*)nres ! Healpix parameter case ('-d') outdir=arg ! output dir case ('-f') fname=arg ! output file name case ('-th') read(arg,*)thick ! RayTrace parameter case ('-kp') read(arg,*)kplane ! RayTrace parameter *** >0 *** case ('-wde') read(arg,*)wde case ('-nt') read(arg,*)nopenmp end select enddo !!! set parameters for openmp !$ call omp_set_num_threads(nopenmp) if (kplane.le.0) then write(*,*)'*** ERROR ***' write(*,*)'*** kp MUST be > 0 ***' stop endif !!! Healpix related nside=2**nres npix=12_i8b*nside*nside write(*,*)'#' write(*,*)'### HEALPIX parameters' write(*,*)'#' write(*,*)'nres =',nres write(*,*)'nside =',nside write(*,*)'npix =',npix allocate(sdens(0:npix-1)) !$omp parallel private(l) !$omp do do ipix=0,npix-1 sdens(ipix)=0.0 enddo !$omp end do !$omp end parallel !!! Gadget related file path path=path(1:len_trim(path)) !!! RayTrace related rmin=(kplane-1)*thick rmax=kplane*thick rmin2=rmin*rmin rmax2=rmax*rmax ! !!!!!! read Gadget snapshot !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! first we just see how many particles there are of each type write(snumber,'(I3)')SNAPSHOT write(fnumber,'(I3)')0 do i=1,3 if (snumber(i:i) .eq. ' ') then snumber(i:i)='0' end if end do dirname=path(1:len_trim(path)) if (FILES.eq.1) then filename= dirname(1:len_trim(dirname)) // '/snapshot_' // $ snumber(1:len_trim(snumber)) else filename= dirname(1:len_trim(dirname)) // '/snapshot_' // $ snumber(1:len_trim(snumber)) // '.' // $ fnumber(verify(fnumber,' '):3) endif write(*,*)'#' write(*,*)'### reading GADGET header from snapshot :', $ filename(1:len_trim(filename)) write(*,*)'#' ! now, read in the header open (1, file=filename, form='unformatted') read(1)npart,massarr,a,redshift,flag_sfr,flag_feedback,nall $ ,flag_cooling,FILES,boxsize,omegam,OmegaLambda,HubbleParam $ ,unused close(1) write(*,*)'npart =',npart write(*,*)'massarr =',massarr write(*,*)'a =',a write(*,*)'redshift =',redshift write(*,*)'flag_sfr =',flag_sfr write(*,*)'flag_feedback =',flag_feedback write(*,*)'nall =',nall write(*,*)'flag_cooling =',flag_cooling write(*,*)'FILES =',FILES write(*,*)'boxsize =',boxsize write(*,*)'omegam =',omegam write(*,*)'OmegaLambda =',OmegaLambda write(*,*)'HubbleParam =',HubbleParam write(*,*)'unused =',unused !!! setup Gadget related parameters boxsize=boxsize/1000.0 !kpc/h => Mpc/h boxinv=1.0/boxsize boxvol=boxsize**3 ! in (Mpc/h)^3 !!! setup raytracing related parameters shellvol=4.0*PI*(rmax**3-rmin**3)/(8.0*3.0) !!! 1/8 because of shell-volume in Octant !!! shellpart=np**3*shellvol/boxvol pixpart=shellpart/(npix/8.0) dcent=(kplane-0.5)*thick ! Mpc/h dlplane=0.75*(rmax**4-rmin**4)/(rmax**3-rmin**3) ! particle number weighted d_lensplane in Mpc/h ! change made in v1.1 write(*,*)'#' write(*,*)'### RayTracing parameters' write(*,*)'#' write(*,*)'kplane =',kplane write(*,*)'shell center =',dcent write(*,*)'weighted dlens =',dlplane write(*,*)'pixpart =',pixpart ain=a dlens=dlplane/chubble0 call newlap(dlens,aout,ain) alens=aout ! lens redshift zlens=1.0/alens-1.0 write(*,*)'z_plane =',zlens sdensnorm=1.5*omegam*boxvol*npix/ $ (alens*dlplane*chubble0**2*np**3*4.0*PI) write(*,*)'sdensnorm =',sdensnorm ! Now we just read in the coordinates, and only keep those of type=1 write(*,*)'#' do fn=0,FILES-1 write(snumber,'(I3)') SNAPSHOT write(fnumber,'(I4)') fn do i=1,3 if (snumber(i:i) .eq. ' ') then snumber(i:i)='0' end if end do if (FILES.eq.1) then filename= dirname(1:len_trim(dirname)) // '/snapshot_' // $ snumber(1:len_trim(snumber)) else filename= dirname(1:len_trim(dirname)) // '/snapshot_' // $ snumber(1:len_trim(snumber)) // '.' // $ fnumber(verify(fnumber,' '):4) endif ! now, read in the header open(1,file=filename,form='unformatted') read(1)npart,massarr,a,redshift,flag_sfr,flag_feedback,nall $ ,flag_cooling,FILES,boxsize,omegam,OmegaLambda,HubbleParam $ ,unused close(1) boxsize=boxsize/1000.0 !kpc/h => Mpc/h ! now, read in the snapshot data open(1,file=filename,form='unformatted', $ access='direct',recl=4) nsnap=sum(npart) write(*,*)'### reading Gadget snapshot :', $ filename(1:len_trim(filename)) write(*,*)'Num of particles in this file =',nsnap do i=1,nsnap read(1,rec=3*i+65)pos(1) read(1,rec=3*i+66)pos(2) read(1,rec=3*i+67)pos(3) c write(*,*)pos forall (ii=1:3) pos(ii)=1.e-3*pos(ii) ! kpc => Mpc !!!!!! place particles onto 8 octant-spheres do kz=1,2 if (kz.eq.1) then p3=pos(3) else p3=pos(3)-boxsize endif p33=p3*p3 if (p33.ge.rmax2) cycle do ky=1,2 if (ky.eq.1) then p2=pos(2) else p2=pos(2)-boxsize endif p23=p2*p2+p33 if (p23.ge.rmax2) cycle do kx=1,2 if (kx.eq.1) then p1=pos(1) else p1=pos(1)-boxsize endif rr=p1*p1+p23 if (rr.ge.rmin2.and.rr.lt.rmax2) then hpvector(1)=p1 hpvector(2)=p2 hpvector(3)=p3 call vec2pix_ring(nside,hpvector,ipring) sdens(ipring)=sdens(ipring)+1.0 endif enddo enddo enddo enddo c read(1)idamy close(1) enddo write(*,*)'#' !!! setup output filename write(pnumber,'(I3)')kplane do i=1,3 if (pnumber(i:i).eq.' ') then pnumber(i:i)='0' end if end do pfname=fname(1:len_trim(fname))// $ 'p'//pnumber(1:len_trim(pnumber)) $ write(*,*)'output file name =',pfname(1:len_trim(pfname)) write(*,*)'pixpart =',pixpart ! ! write log ! outlog=ADJUSTL(ADJUSTR(pfname)//'.log') outlog=ADJUSTL(ADJUSTR(outdir)//'/'//outlog) open(10,file=outlog,status='unknown') write(10,*)'pfname = ', $ pfname(1:len_trim(pfname)) write(10,*)'path = ',path(1:len_trim(path)) write(10,*)'SNAPSHOT =',SNAPSHOT write(10,*)'np = ',np write(10,*)'boxsize = ',boxsize write(10,*)'nres =',nres write(10,*)'nside =',nside write(10,*)'npix =',npix write(10,*)'kplane =',kplane write(10,*)'thick =',thick write(10,*)'dlens =',dlplane write(10,*)'alens =',alens write(10,*)'zlens =',1.0/alens-1.0 write(10,*)'dcent =',dcent write(10,*)'pixpart =',pixpart write(10,*)'omegam =',omegam write(10,*)'wde =',wde close(10) ! !!!!!!!!!!!! spherical harmonics transform --- sdens => alm_sdens !!!!!!!!!!!!!!!! ! lmax=3*nside allocate(hpalm(1:1,0:lmax,0:lmax)) allocate(dw8(1:2*nside,1:1)) allocate(hpcl(0:lmax,1:1)) write(*,*)'*** spherical harmonics transformation for ', $ pfname(1:len_trim(pfname)) !!! particle number density => delta Sigma forall (ipix=0:npix-1) sdens(ipix)=sdensnorm*(sdens(ipix)-pixpart) outfile=ADJUSTL(ADJUSTR(pfname)//'.spkap.dat') outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile) open(10,file=outfile,form='unformatted',STATUS='unknown') write(10)nside,npix write(10)sdens close(10) dw8=1.0_dp call map2alm(nside,lmax,lmax,sdens,hpalm,hpz,dw8) deallocate(sdens) deallocate(dw8) outfile=ADJUSTL(ADJUSTR(pfname)//'.spkap.alm') outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile) open(10,file=outfile,form='unformatted',STATUS='unknown') write(10)nside,npix,lmax write(10)hpalm close(10) !!! alm_kappa => cl call alm2cl(lmax,lmax,hpalm,hpcl) outfile=ADJUSTL(ADJUSTR(pfname)//'.spkap.cl') outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile) open(10,file=outfile,status='unknown') do i=1,lmax write(10,*)i,hpcl(i,1) enddo close(10) deallocate(hpcl) !!!!!!!!! alm_kappa => alm_phi (alm_phi=2*alm_kappa/l(l+1)) !!!!!!!!!!!! write(*,*)'*** spkap_lm => psi,i, psi,ij: processing...' !$omp parallel private(l) !$omp do do m=0,lmax hpalm(1,0,m)=0.0 !psi_lm=0 for l=0 do l=1,lmax hpalm(1,l,m)=hpalm(1,l,m)/(0.5*l*(l+1)) ! psi_lm = 2 kappa_lm/l/(1+1) for l not= 0 enddo enddo !$omp end do !$omp end parallel allocate(hpmap_phi(0:npix-1)) allocate(hpmap_alp(0:npix-1,1:2)) allocate(hpmap_mag(0:npix-1,1:3)) call alm2map_der(nside,lmax,lmax,hpalm,hpmap_phi, $ hpmap_alp,hpmap_mag) deallocate(hpmap_phi) deallocate(hpalm) ccc changing hpmap_alp(0:npix-1,1:2) => hpmap(1:2,0:npix-1) allocate(hpmap(1:2,0:npix-1)) !$omp parallel private(j) !$omp do do ipix=0,npix-1 do j=1,2 hpmap(j,ipix)=hpmap_alp(ipix,j) enddo enddo !$omp end do !$omp end parallel deallocate(hpmap_alp) !!!!!!!! output PSIs outfile=ADJUSTL(ADJUSTR(pfname)//'.psis') outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile) write(*,*)'#' write(*,*)'### output' write(*,*)'#' write(*,*)'output PSIs file = ',outfile(1:len_trim(outfile)) open(10,file=outfile,status='unknown',form='unformatted') write(10)nside,npix write(10)hpmap deallocate(hpmap) ccc changing hpmap_mag(0:npix-1,1:3) => hpmap(1:3,0:npix-1) allocate(hpmap(1:3,0:npix-1)) !$omp parallel private(j) !$omp do do ipix=0,npix-1 do j=1,3 hpmap(j,ipix)=hpmap_mag(ipix,j) enddo enddo !$omp end do !$omp end parallel deallocate(hpmap_mag) write(10)hpmap close(10) write(*,*)'*** gadget2psis done' end program gadget2psis cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc Newton-Raphson method cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine newlap(alm,aout,ain) implicit real (a-h,o-z) implicit integer (i-n) parameter (eps=1.0e-5) external dldz rhs=alm a=ain 356 call qromb(dldz,a,1.0,solv) err=abs(solv/rhs-1.0) c write(*,*)solv,err if (err.gt.eps) then dela=-(rhs-solv)/dldz(a) a2=a+dela if (a2.le.0.0) a2=a+0.5*dela a=a2 c write(*,*)a goto 356 endif aout=a return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccc package of solving the integral equation of ccc affine parameter - redshit (scale factor) relation cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function dldz(a) use cosmoparam implicit real (a-h,o-z) implicit integer (i-n) dldz=1.0/sqrt(a*omegam+(1.0-omegam)*a**(1.0-3.0*wde)) return end ccc SUBROUTINE qromb(func,a,b,ss) INTEGER JMAX,JMAXP,K,KM REAL a,b,func,ss,EPS EXTERNAL func PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1) CU USES polint,trapzd INTEGER j REAL dss,h(JMAXP),s(JMAXP) h(1)=1. do 11 j=1,JMAX call trapzd(func,a,b,s(j),j) if (j.ge.K) then call polint(h(j-KM),s(j-KM),K,0.,ss,dss) if (abs(dss).le.EPS*abs(ss)) return endif s(j+1)=s(j) h(j+1)=0.25*h(j) 11 continue write(*,*)'too many steps in qromb' stop END ccc SUBROUTINE trapzd(func,a,b,s,n) INTEGER n REAL a,b,s,func EXTERNAL func INTEGER it,j REAL del,sum,tnm,x if (n.eq.1) then s=0.5*(b-a)*(func(a)+func(b)) else it=2**(n-2) tnm=it del=(b-a)/tnm x=a+0.5*del sum=0. do 11 j=1,it sum=sum+func(x) x=x+del 11 continue s=0.5*(s+(b-a)*sum/tnm) endif return END ccc SUBROUTINE polint(xa,ya,n,x,y,dy) INTEGER n,NMAX REAL dy,x,y,xa(n),ya(n) PARAMETER (NMAX=10) INTEGER i,m,ns REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX) ns=1 dif=abs(x-xa(1)) do 11 i=1,n dift=abs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i)=ya(i) d(i)=ya(i) 11 continue y=ya(ns) ns=ns-1 do 13 m=1,n-1 do 12 i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if(den.eq.0.)then write(*,*)'failure in polint' stop endif den=w/den d(i)=hp*den c(i)=ho*den 12 continue if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 endif y=y+dy 13 continue return END