!---------------------------------------------------------------------------------------- ! Find properties of large distinct halos ! ! Input: Halo Catalog (short list) ! ! Output: Distribution functions ! Parameters: ! ! ! MODULE setHalos PUBLIC Integer*4, Parameter :: NhaloMax = 45000000, & ! maximum number of halos NinBin = 600, & ! threshold to start increasing of bin size NtabM = 10000 ! Max size of Pk table Real, PARAMETER :: & dlog = 0.1, & dVLog = 0.05, & dMlog = 0.10, & factorlog = 1.5, & ! increment for variable binning slog_start= 0.02, & ! starting bin size for Vcirc/Vvir mass selection vlog_start= 0.010, & ! starting bin size for Vcirc/Vvir velocity selection dBuffer = 2.5 Integer*8, DIMENSION(NhaloMax) :: LDistinct,HaloNumber,indR Real, DIMENSION (NhaloMax) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc, & Xoff,VirRatio,CAratio,aLambda, & DD,SD ! buffers Real*4 :: deltac22,deltac,Uklow,Ukup , & xkt(0:NtabM),Pkt(0:NtabM),StepK,alog0 Integer :: Ncenter,iFlag,iSnap,Ntab Real*8 :: Box, Rmax, Aexpn, Om0,Omb, Ovdens, sigma8, hsmall, & OmLOm0 Real*8 :: Grf,rf Character :: Letter*2 Contains !---------------------------------------------------- SUBROUTINE Gr(Growth,a) !---------------------------------------------------- ! use SetHalos Real*8 :: a,ww,x,Grow0,Hgrow,Growth,err External Hgrow x = (OmLOm0)**0.33333 ! Normalize to growth at a=1 Call QAG(Hgrow, 1.d-6, x, 1.d-6,1.d-6,2,ww,err,Neval,ier) Grow0 =sqrt(1.+x**3)/x**1.5*ww x = (OmLOm0)**0.33333*a Call QAG(Hgrow, 1.d-6, x, 1.d-6,1.d-6,2,ww,err,Neval,ier) Growth =sqrt(1.+x**3)/x**1.5*ww/Grow0 End SUBROUTINE Gr !-------------------------------------- rms fluctuations ! for mass aMass ! aMass is in Msun/h units REAL*8 FUNCTION LogSlope(aMass) !--------------------------------------- !Use SetHalos Real*8, PARAMETER :: pi = 3.1415926535d0, xx =0.691d0, dM =1.e-1 REAL*8 :: a,aMass,dd,err,dm1,ds1,M1,M2,s1,s2,s REAL*8, EXTERNAL :: Ptophat rf= (aMass/(1.162e12*Om0))**0.33333333 ! in Mpc/h units Call QAG(Ptophat, 1.d-6, 1.d3, 1.d-6,1.d-6,2,dd,err,Neval,ier) s =GrF*sqrt((dd)/rf**3/(2.*pi**2)) M1 = aMass*(1.+dM) rf= (M1/(1.162e12*Om0))**0.33333333 ! in Mpc/h units Call QAG(Ptophat, 1.d-6, 1.d3, 1.d-6,1.d-6,2,dd,err,Neval,ier) s1 =GrF*sqrt((dd)/rf**3/(2.*pi**2)) M2 = aMass*(1.-dM) rf= (M2/(1.162e12*Om0))**0.33333333 ! in Mpc/h units Call QAG(Ptophat, 1.d-6, 1.d3, 1.d-6,1.d-6,2,dd,err,Neval,ier) s2 =GrF*sqrt((dd)/rf**3/(2.*pi**2)) LogSlope = -aMass/s*(s1-s2)/(M1-M2) END FUNCTION LogSlope !-------------------------------------- rms fluctuations ! for mass aMass ! aMass is in Msun/h units REAL*8 FUNCTION SigMass2(aMass) !--------------------------------------- !Use SetHalos Real*8, PARAMETER :: pi = 3.1415926535d0, xx =0.691d0 REAL*8 :: a,aMass,dd,err REAL*8, EXTERNAL :: Ptophat rf= (aMass/(1.162e12*Om0))**0.33333333 ! in Mpc/h units Call QAG(Ptophat, 1.d-6, 1.d3, 1.d-6,1.d-6,2,dd,err,Neval,ier) SigMass2 =GrF*sqrt((dd)/rf**3/(2.*pi**2)) END FUNCTION SigMass2 !--------------------------------------------------- ! ! Setup for the power sapectrum ! Subroutine SetPk !---------------------------------------------------- Real*8, parameter :: pi=3.141592653 Real*8 :: sig8,a,wk REAL*8 :: Ptophat,Ppk,FF,Growth,aM external Ptophat Call ReadTable !-------- Uchuu Om0 = 0.3089 OmL = 1.-Om0 OmLOm0 = OmL/Om0 Aexpn = 1.00 ! Pkt(:) = Pkt(:)*1.023 !!!! scale up Pk for bigMD sim a = Aexpn Call Gr(Growth,a) GrF = Growth dM =1.05 rf = 8. aM = 1.162e12*Om0*rf**3 Sig8 = SigMass2(aM) Write(*,'(1p,3(a,g14.4))') 'a = ',a,' GrFactor =',GrF Write(*,'(1p,3(a,g14.4))') 'M = ',aM, ' R = ',rf,' Sigma8 = ',Sig8 write(*,'(a,$)') ' Do you want to change Sigma8? [New value or "return" for "no"] ' READ (*,'(g12.4)',IOSTAT=iStat) Sig8new if(iStat.eq.0.and.Sig8new.gt.0.) Then Pkt(:) = (Sig8new/Sig8)**2*Pkt(:) sig8 = SigMass2(aM) write(*,*) ' New sigma_8 = ', sig8,Sig8New Else write(*,*) ' No change in sigma_8' End if end Subroutine SetPk !--------------------------------------- SUBROUTINE ReadTable ! interpolate table with p(k) !--------------------------------------- !Use SetHalos Character*120 :: Name,Line Real*8 :: a,ww,x,Grow0,Hgrow,Growth,Omc,OmL,B Logical :: exst Name = 'PkTable.dat' Open(1,file=TRIM(Name)) Call ParseHeader(Omb,B) Call ParseHeader(Omc,B) Call ParseHeader(OmL,B) Call ParseHeader(Om0,B) Call ParseHeader(A,sigma8) write(*,'(a,4g12.4)') 'Omb =',Omb,Omc,A,sigma8 OmLOm0 = OmL/Om0 j = 0 Do i=1,NtabM Read(1,*,iostat=ierr) ak,Pk If(ierr/=0)exit xkt(i) = ak Pkt(i) = Pk j = j+1 EndDo Ntab = j StepK = log10 (xkt (Ntab) / xkt (1) )/(Ntab-1) alog0 = log10 (xkt (1) ) ! test that spacing is the same Do k = 2, Ntab - 1 ss = log10 (xkt (k + 1) / xkt (k) ) If (abs (ss / StepK - 1.) .gt.2.e-2) Then Write (*,*) ' error in K spacing. k=', k, xkt (k + 1),xkt (k) STOP EndIf EndDo write(*,*) ' Number of lines read =',Ntab END SUBROUTINE ReadTable !--------------------------------------------- ! subroutine ParseHeader(A,B) CHARACTER*150 :: Str1 Real*8 :: A,B read(1,'(a)') Str1 ipos = INDEX(Str1,'=') ! First equal position Str1 = Str1(ipos+1:) ! shift left read(Str1,*) A ipos = INDEX(Str1,'=') ! Second equal position if(ipos==0)Then B = 0. Else Str1 = Str1(ipos+1:) ! shift left read(Str1,*) B End if end subroutine ParseHeader !-------------------------------------------------------------- ! Reads and recognizes files ! SUBROUTINE GetFile(iFile) !-------------------------------------------------------------- Integer*4 , parameter :: Nbin = 1000 Integer*4, SAVE :: iFileout Logical :: exst Character*120 :: file1,file2,catname,listname,HEADER*45,line,line2 Open(20,file='ErrorAnalyzeCatalogs.dat',position='append') If(iFile == 0) Then ! recongize configuration write(*,'(a,$)')' Enter snapshot number, Catalog letter, and Box size => ' read(*,*)iSnap,Letter,Box write(file2,'(a,2a1,i4.4,a)')'HaloStaTot',TRIM(Letter),'.',iSnap,'.DAT' Open( 2,file=TRIM(file2)) ! output iFileout = 0 write(file1,'(a,2a1,i4.4,a1,i2.2,a)')'Catshort',TRIM(Letter),'.',iSnap, & '.',iFileout,'.DAT' write(*,'(a)') TRIM(file1) Inquire(file=TRIM(file1),exist = exst) write(*,*) ' Status =',exst If(exst)Then ! multiple files CatshortX.XXX.XX.DAT iFlag = 0 open(3,file = TRIM(file1)) Else iFlag = 1 !--- single file CatshortX.xxx.DAT write(file1,'(a,2a1,i4.4,a)')'Catshort',TRIM(Letter),'.',iSnap, & '.DAT' Inquire(file=TRIM(file1),exist = exst) write(*,'(a)') TRIM(file1) write(*,*) ' Status =',exst If(exst)Then open(3,file = TRIM(file1)) Else Stop '---- Input Catshort file not found ' End If End If read(3,'(a)')HEADER !read(3,'(2(a7,f8.5))') line(1:7),AEXPN,Line2(1:7),Step read(3,'(2(a8,f8.5))') line(1:8),AEXPN !,Line2(1:7),Step write(*,*) AEXPN rewind(3) Else If(iFlag == 1)Then ! single file iFile = 0 Return Else iFileout = iFileout +1 write(file1,'(a,2a1,i4.4,a1,i2.2,a)')'Catshort',TRIM(Letter),'.',iSnap, & '.',iFileout,'.DAT' Inquire(file=TRIM(file1),exist = exst) If(exst)Then ! multiple files CatshortX.XXX.XX.DAT open(3,file = TRIM(file1)) Else ! no more files to read iFile = 0 Return End If End If End If end SUBROUTINE GetFile !-------------------------------------------------------------- ! Mass function of distinct halos ! SUBROUTINE VelocityFunction !-------------------------------------------------------------- Integer, parameter :: Nbin = 1000 Integer, Dimension(-Nbin:Nbin) :: & mHalosD, iHalosD, vHalosD, & mHalosS, iHalosS, vHalosS Real, Dimension(-Nbin:Nbin) :: massbinD,vbinD,massbinS,vbinS vHalosD = 0 ; vbinD = 0. vHalosS = 0 ; vbinS = 0. nD=0 ; nS =0 ! count all halos Do i=1,Ncenter vv = sqrt(VXc(i)**2+VYc(i)**2+VZc(i)**2) If(LDistinct(i) == 0 )Then nD =nD +1 ind =INT(log10(vv)/dVlog+Nbin)-Nbin If(indNbin)Then write(*,*) ' !!!! Error in Catshort file !!!!' write(*,'(a,2i12,a,G12.4)')' Halo=',i,ind,' Vcirc=',Vcirc(i) stop EndIf If(ind0.or.vHalosS(i)>0)Then istart =i ; exit EndIf EndDo ifinish =Nbin Do i =Nbin,-Nbin,-1 If(vHalosD(i)>0)Then ifinish =i ; exit EndIf EndDo Do i =ifinish-1,istart,-1 vHalosD(i) = vHalosD(i) +vHalosD(i+1) vHalosS(i) = vHalosS(i) +vHalosS(i+1) EndDo write(2,'(/2a,T90,a)') 'Vdrift:V(km/s) N(>V) n(>V)Mpch3 V(km/s) dN VdN/dV', & ' Sub: N(>V) n(>V)Mpch3 V(km/s) dN VdN/dV' Do i=istart,ifinish ! print results aVm =10.**(i*dVlog) dV = 10.**((i+1)*dVlog)-aVm If(aVm>30.)Then write(2,'(3x,1P,G12.5,2(2g11.4,2x,1P,4g11.4))')aVm, vHalosD(i),vHalosD(i)/(Box)**3, & vbinD(i),iHalosD(i),vbinD(i)*iHalosD(i)/dV/(Box)**3, & vHalosS(i),vHalosS(i)/(Box)**3,vbinS(i),iHalosS(i),vbinS(i)*iHalosS(i)/dV/(Box)**3 End If EndDo End SUBROUTINE VelocityFunction !-------------------------------------------------------------- ! ! SUBROUTINE VelocityDistr !-------------------------------------------------------------- Integer, parameter :: Nbin = 1000 Real*4, Dimension(-Nbin:Nbin) :: & mHalosDx, iHalosDx, vHalosDx, & mHalosSx, iHalosSx, vHalosSx, & mHalosDy, iHalosDy, vHalosDy, & mHalosSy, iHalosSy, vHalosSy, & mHalosDz, iHalosDz, vHalosDz, & mHalosSz, iHalosSz, vHalosSz Real, Dimension(-Nbin:Nbin) :: vbinDx,vbinSx,vbinDy,vbinSy,vbinDz,vbinSz vHalosDx = 0 ; vbinDx = 0. vHalosSx = 0 ; vbinSx = 0. vHalosDy = 0 ; vbinDy = 0. vHalosSy = 0 ; vbinSy = 0. vHalosDz = 0 ; vbinDz = 0. vHalosSz = 0 ; vbinSz = 0. nD=0 ; nS =0 ! count all halos Do i=1,Ncenter vvx = abs(VXc(i))+1.e-10 vvy = abs(VYc(i))+1.e-10 vvz = abs(VZc(i))+1.e-10 If(LDistinct(i) == 0 )Then nD =nD +1 ind =INT(log10(vvx)/dVlog+Nbin)-Nbin if(ind<-1000)write(*,'(i9,3f9.3,i12)') i,VXc(i),VYc(i),VZc(i),ind If(indNbin)Then write(*,*) ' !!!! Error in Catshort file !!!!' write(*,'(a,2i12,a,G12.4)')' Halo=',i,ind,' Vcirc=',Vcirc(i) stop EndIf If(ind0.or.mHalosS(i)>0)Then istart =i ; exit EndIf EndDo ifinish =Nbin Do i =Nbin,-Nbin,-1 If(mHalosD(i)>0)Then ifinish =i ; exit EndIf EndDo iHalosD = mHalosD iHalosS = mHalosS Do i =ifinish-1,istart,-1 mHalosD(i) = mHalosD(i) +mHalosD(i+1) mHalosS(i) = mHalosS(i) +mHalosS(i+1) EndDo write(2,'(/3a)') ' Distinct Mass(Msunh) sigma N(>M) n(>M)/Mpch3',& ' Mass sigma Slope dN MdN/dM', & ' Sub: N(>M) n(>M)/Mpch Mass dN MdN/dM sigma' Do i=istart,ifinish ! print results aMm =10.**(i*dMlog) dM = 10.**((i+1)*dMlog)-aMm sigma = SigMass2(aMm) aMd = massbinD(i) sigmaD= SigMass2(aMd) Slope = LogSlope(aMd) write(2,'(9x,1P,2G11.4,6x,2(2g12.5,2x,1P,10g11.4))')aMm,sigma, mHalosD(i),mHalosD(i)/(Box)**3, & massbinD(i),sigmaD,Slope,iHalosD(i), massbinD(i)*iHalosD(i)/dM/(Box)**3, & mHalosS(i),mHalosS(i)/(Box)**3, & massbinS(i),iHalosS(i),massbinS(i)*iHalosS(i)/dM/(Box)**3,sigma EndDo ! ------------------------------- print mass function of Distinct and Subhalos iHalosD = vHalosD iHalosS = vHalosS istart =-Nbin ! find first and last non-empty bin Do i =-Nbin,Nbin If(vHalosD(i)>0.or.vHalosS(i)>0)Then istart =i ; exit EndIf EndDo ifinish =Nbin Do i =Nbin,-Nbin,-1 If(vHalosD(i)>0)Then ifinish =i ; exit EndIf EndDo Do i =ifinish-1,istart,-1 vHalosD(i) = vHalosD(i) +vHalosD(i+1) vHalosS(i) = vHalosS(i) +vHalosS(i+1) EndDo write(2,'(/2a,T90,a)') 'Dist:V(km/s) N(>V) n(>V)Mpch3 V(km/s) dN VdN/dV', & ' Sub: N(>V) n(>V)Mpch3 V(km/s) dN VdN/dV' Do i=istart,ifinish ! print results aVm =10.**(i*dVlog) dV = 10.**((i+1)*dVlog)-aVm write(2,'(3x,1P,G12.5,2(2g11.4,2x,1P,4g11.4))')aVm, vHalosD(i),vHalosD(i)/(Box)**3, & vbinD(i),iHalosD(i),vbinD(i)*iHalosD(i)/dV/(Box)**3, & vHalosS(i),vHalosS(i)/(Box)**3,vbinS(i),iHalosS(i),vbinS(i)*iHalosS(i)/dV/(Box)**3 EndDo End SUBROUTINE MassFunction !-------------------------------------------------------------- ! Vcirc -Vvir relation for Distinct Mass selection ! SUBROUTINE VmaxVvirM !-------------------------------------------------------------- Integer, Parameter :: Nbin =1500 Real*8 :: ss write(2,'(/3a,i10)')'Ndistinct sigma Vmax/Vvir: 10% 32%', & ' Median 68 90% Concentration C(10%) C(32%) C(median) C(68%) C(90%) Mass Selected' slog = slog_start aMstart =1.e10 aMin = aMstart aMax = aMin*10.**(slog) !write(*,*) ' ---------- Mass selected Concentrations ',Aexpn Do ibin =1,Nbin ! ibin If(aMax>5.e15)exit !If(vMax>3000.or.vMin<50.)exit iC =0 ss =0. DD(:) = 0 Do i=1,Ncenter If(aMc(i)>=aMin.and.aMc(i)=vMin.and.Vcirc(i)NhaloMax)Stop ' Not enough space for buffer. Error' Vvir = 2.081e-3*sqrt(aMc(i)/Rmc(i)/Aexpn) ! If(aMc(i)>5.e11)write(*,'(i8,3f10.4,1p,8g12.5)')HaloNumber(i),& ! Xc(i),Yc(i),Zc(i),aMc(i),Rmc(i),Vcirc(i),Vvir DD(iC) = max(Vcirc(i)/Vvir,1.) ss = ss + aMc(i) EndIf EndIf EndDo ! Ncenter If(iC>10)Then CALL R_mrgref (DD,indR,iC) M10 = INT(0.1*iC)+1 M32 = INT(0.32*iC)+1 M50 = INT(0.5*iC)+1 M68 = INT(0.68*iC)+1 M90 = INT(0.9*iC)+1 V10 = DD(indR(M10)) V32 = DD(indR(M32)) V50 = DD(indR(M50)) V68 = DD(indR(M68)) V90 = DD(indR(M90)) EndIf !y = 1.e12/(ss/iC) !sigma = Da*16.9*y**0.41/(1.+1.102*y**0.20+6.22*y**0.333) !write(*,'(a,g12.4,i9,2g12.4)') ' ss=',ss,iC !write(*,*) ss/iC,sigma,sigmaNew If(iC>10)Then sigma = SigMass2(ss/iC) write(2,'(i8,1p,2g12.5,7x,18g12.5)') iC,ss/iC,sigma,V10,V32,V50,V68,V90, & Concentration(V10),Concentration(V32),Concentration(V50), & Concentration(V68),Concentration(V90) EndIf If(iCaMstart)Then slog =factorlog*slog End If aMin = aMax aMax = aMin*10.**slog EndDo ! change selection criteria ------------------- Xoffmax = 0.075 CAmin = 0.50 aLambdaMax = 0.1 VirRatLimit = 1.5 write(2,'(/3(a,f8.4))') ' Vmax/Vvir ratio: for Xoffset <',Xoffmax, & ' and axial ratio c/a> ',CAmin write(2,'(3(a,f8.4))') ' Lambda <',aLambdaMax, & ' and 2Kin/Upot-1 < ',VirRatLimit write(2,'(2a)')'Ndistinct sigma Vmax/Vvir: 10% 32% ', & ' Median 68% 90% Concentration C(10%) C(32%) C(median) C(68%) C(90%) Mass Selected' slog = slog_start aMin = aMStart aMax =1.e9*10.**(slog) Do ibin =1,Nbin If(aMax>3.e15)exit iC =0 ss =0. DD(:) =0 Do i=1,Ncenter If(aMc(i)>=aMin.and.aMc(i)NhaloMax)Stop ' Not enough space for buffer. Increase M' Vvir = 2.081e-3*sqrt(aMc(i)/Rmc(i)/Aexpn) DD(iC) = max(Vcirc(i)/Vvir,1.) ss = ss + aMc(i) EndIf EndIf EndIf EndDo If(iC>10)Then CALL R_mrgref (DD,indR,iC) M10 = INT(0.1*iC)+1 M32 = INT(0.32*iC)+1 M50 = INT(0.5*iC)+1 M68 = INT(0.68*iC)+1 M90 = INT(0.9*iC)+1 V10 = DD(indR(M10)) V32 = DD(indR(M32)) V50 = DD(indR(M50)) V68 = DD(indR(M68)) V90 = DD(indR(M90)) EndIf !y = 1.e12/(ss/iC) !sigma = Da*16.9*y**0.41/(1.+1.102*y**0.20+6.22*y**0.333) If(iC>10)Then sigma = SigMass2(ss/iC) write(2,'(i8,1p,2g12.5,7x,18g12.5)') iC,ss/iC,sigma,V10,V32,V50,V68,V90, & Concentration(V10),Concentration(V32),Concentration(V50), & Concentration(V68),Concentration(V90) EndIf If(iCaMstart)slog =factorlog*slog aMin = aMax aMax = aMin*10.**slog EndDo End SUBROUTINE VmaxVvirM !-------------------------------------------------------------- ! Vcirc -Vvir relation for Distinct Vcirc selection ! SUBROUTINE VmaxVvir !-------------------------------------------------------------- Integer, Parameter :: Nbin =1500 Real*8 :: ss write(2,'(/3a,i10)')'Ndistinct sigma Vmax Vmax/Vvir: 10% 32% ', & ' Median 68% 90% Concentration C(10%) C(32%) C(median) C(68%) C(90%) Vcirc Selected' slog = vlog_start vStart = 50. vMin = vStart vMax = vMin*10.**(slog) write(*,*) ' a=',AEXPN Do ibin =1,Nbin ! write(*,'(i4,1p,8g12.4)') ibin,vMin,vMax,aMin,aMax,slog If(vMax>3000.or.vMin<50.)exit iC =0 ss =0. DD(:)=0 Do i=1,Ncenter If(Vcirc(i)>=vMin.and.Vcirc(i)NhaloMax)Stop ' Not enough space for buffer. Error' Vvir = 2.081e-3*sqrt(aMc(i)/Rmc(i)/Aexpn) ! If(aMc(i)>5.e11)write(*,'(i8,3f10.4,1p,8g12.5)')HaloNumber(i),& ! Xc(i),Yc(i),Zc(i),aMc(i),Rmc(i),Vcirc(i),Vvir DD(iC) = max(Vcirc(i)/Vvir,1.) ss = ss + aMc(i) EndIf EndIf EndDo If(iC>10)Then CALL R_mrgref (DD,indR,iC) M10 = INT(0.1*iC)+1 M32 = INT(0.32*iC)+1 M50 = INT(0.5*iC)+1 M68 = INT(0.68*iC)+1 M90 = INT(0.9*iC)+1 V10 = DD(indR(M10)) V32 = DD(indR(M32)) V50 = DD(indR(M50)) V68 = DD(indR(M68)) V90 = DD(indR(M90)) EndIf !y = 1.e12/(ss/iC) !sigma = Da*16.9*y**0.41/(1.+1.102*y**0.20+6.22*y**0.333) If(iC>10)Then sigma = SigMass2(ss/iC) write(2,'(i8,1p,3g12.5,7x,12g12.5)') iC,ss/iC,sigma,(vMin+Vmax)/2.,V10,V32,V50,V68,V90, & Concentration(V10),Concentration(V32),Concentration(V50), & Concentration(V68),Concentration(V90) EndIf If(iC1.e12)slog =factorlog*slog vMin = vMax vMax = vMin*10.**slog EndDo ! change selection criteria ------------------- Xoffmax = 0.075 CAmin = 0.50 aLambdaMax = 0.1 VirRatLimit = 1.5 write(2,'(/3(a,f8.4))') ' Vmax/Vvir ratio: for Xoffset <',Xoffmax, & ' and axial ratio c/a> ',CAmin write(2,'(3(a,f8.4))') ' Lambda <',aLambdaMax, & ' and 2Kin/Upot-1 < ',VirRatLimit write(2,'(2a)')'Ndistinct sigma Vmax/Vvir: 10% 32% ', & ' Median 68% 90% Concentration C(10%) C(32%) C(median) C(68%) C(90%) Vcirc Selected' slog = vlog_start vMin = vStart vMax = vMin*10.**(slog) Do ibin =1,Nbin If(vMax>3000)exit iC =0 ss =0. Do i=1,Ncenter If(Vcirc(i)>=vMin.and.Vcirc(i)5.e11)write(*,'(i8,3f10.4,1p,8g12.5)')HaloNumber(i),& ! Xc(i),Yc(i),Zc(i),aMc(i),Rmc(i),Vcirc(i),Xoff(i),CAratio(i) If(LDistinct(i) == 0.and.Xoff(i)NhaloMax)Stop ' Not enough space for buffer. Increase M' Vvir = 2.081e-3*sqrt(aMc(i)/Rmc(i)/Aexpn) DD(iC) = max(Vcirc(i)/Vvir,1.) ss = ss + aMc(i) EndIf EndIf EndIf EndDo If(iC>10)Then CALL R_mrgref (DD,indR,iC) M10 = INT(0.1*iC)+1 M32 = INT(0.32*iC)+1 M50 = INT(0.5*iC)+1 M68 = INT(0.68*iC)+1 M90 = INT(0.9*iC)+1 V10 = DD(indR(M10)) V32 = DD(indR(M32)) V50 = DD(indR(M50)) V68 = DD(indR(M68)) V90 = DD(indR(M90)) EndIf !y = 1.e12/(ss/iC) !sigma = Da*16.9*y**0.41/(1.+1.102*y**0.20+6.22*y**0.333) If(iC>10)Then sigma = SigMass2(ss/iC) write(2,'(i8,1p,2g12.5,7x,12g12.5)') iC,ss/iC,sigma,V10,V32,V50,V68,V90, & Concentration(V10),Concentration(V32),Concentration(V50), & Concentration(V68),Concentration(V90) EndIf If(iC1.e12)slog =factorlog*slog vMin = vMax vMax = vMin*10.**slog EndDo End SUBROUTINE VmaxVvir !-------------------------------------------------------------- ! Vcirc -Mir relation for Distinct and Subs ! SUBROUTINE VcircMvir !-------------------------------------------------------------- Integer, Parameter :: Nbin =500 write(2,'(/3a,i10)')'Ndistinct V10 Vmedian V90', & ' Nsubs V10 Vmedian V90', & ' Nhalos=',Ncenter Do ibin =1,Nbin aMin = 1.e9*10.**((ibin-1)*dlog) aMax =1.e9*10.**(ibin*dlog) If(aMax>3.e15)exit iC =0 ss =0. iCb =0 sb =0. Do i=1,Ncenter If(aMc(i)>=aMin.and.aMc(i)NhaloMax)Stop ' Not enough space for buffer. Error' DD(iC) = Vcirc(i) ss = ss + aMc(i) Else if(LDistinct(i) /= 0)Then iCb =iCb +1 If(iCb>NhaloMax)Stop ' Not enough space for buffer. Increase M' SD(iCb) = Vcirc(i) sb = sb + aMc(i) EndIf EndIf EndIf EndDo If(iC>10)Then CALL R_mrgref (DD,indR,iC) M10 = INT(0.1*iC)+1 M50 = INT(0.5*iC)+1 M90 = INT(0.9*iC)+1 V10 = DD(indR(M10)) V50 = DD(indR(M50)) V90 = DD(indR(M90)) EndIf If(iCb>10)Then CALL R_mrgref (SD,indR,iCb) M10 = INT(0.1*iCb)+1 M50 = INT(0.5*iCb)+1 M90 = INT(0.9*iCb)+1 Vs10 = SD(indR(M10)) Vs50 = SD(indR(M50)) Vs90 = SD(indR(M90)) Else Vs10 =0.; Vs50 =0.; Vs90 =0. EndIf If(iC>10)Then !write(*,'(a,2i8,a,1p,6g12.5)') ' Selected halos=',iC,iCb,' Masses=',aMin,aMax,V50,Vs50 write(2,'(2(i8,1p,4g12.5),6i7)') iC,ss/iC,V10,V50,V90, iCb,sb/max(iCb,1),Vs10,Vs50,Vs90 EndIf EndDo End SUBROUTINE VcircMvir !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalos !-------------------------------------------------------------- Character :: file1*80,Line*120,txt*80 CHARACTER*45 :: HEADER Logical :: inside Ncenter =0 Rdmax =0. factor =1.162e12*Om0*Ovdens ! --- read header of Catshort catalog read(3,'(a)')HEADER write(*,'(a)')HEADER write(2,'(a)')HEADER ! read(3,'(2(a7,f8.5))') txt(1:7),AEXPN,Line(1:7),Step ! write(*,'(2(a7,f8.5))') txt(1:7),AEXPN,Line(1:7),Step ! write(2,'(2(a7,f8.5))') txt(1:7),AEXPN,Line(1:7),Step ! read(3,'(a7,i5,a)') txt(1:7),ISTEP,Line ! write(*,'(a7,i5,a)') txt(1:7),ISTEP,Trim(Line) ! write(2,'(a7,i5,a)') txt(1:7),ISTEP,Trim(Line) ! Do i=1, 14 ! Read(3,'(a)')Line ! if(i .le.13)write(2,'(a)') Line ! write(*,'(a)') Line ! EndDo Do i=1, 16 Read(3,'(a)')Line if(i .le.13)write(2,'(a)') Line write(*,'(a)') Line EndDo Vclimit = 0. iFile = 1 Ndistinct = 0 Do ! i =1,1800000 ! -------- read halos !read (3,*,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ, & ! aM,rr,vrmss,vcrc,vacc,aMacc,aAcc,iBDM,Conc, & ! iDist,Offset,VirRat,aLam,Rin, & ! ba,ca,xax,yax,zax read (3,*,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ, & aM,aMtot,rr,vrmss,vcrc,iHalo,Conc,aNpart, & iDist,Offset,VirRat,aLam,Rin, & ba,ca,xax,yax,zax !read (3,*,iostat=iStat) ii,X0, Y0 , Z0, & ! aM0,aM,vnow,vcrc !rr = aM0**0.33333/4.82e1 ! virial radius in kpch for Om=0.27 delta=360.8 If(mod(Ncenter,100000)==0)& write (*,'(a,i9,3f10.4,2g12.3,3f9.2)') ' Reading: ',Ncenter, X0, Y0 , Z0,aM,aMtot,vcrc,Conc !,aM,rr,vnow,vcrc 45 Format(F9.4,2F9.4,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4) If(iStat.ne.0)Call GetFile(iFile) If(iFile ==0)Exit !inside = (X0>dBuffer.and.X0dBuffer.and.Y0dBuffer.and.Z0 vcLimit .and. rr > 0.95*rvr)Then !!! selection conditions If(vcrc> vcLimit)Then !!! selection conditions Ncenter = Ncenter +1 iHalo = Ncenter !Npart = INT(aNpart) !If(Ncenter /= iHalo) & ! Call HandleError(' --- Catshort :',' Error in sequence of halos:',& ! iHalo,X0, Y0 , Z0, VvX,VvY,VvZ,aM,rr,vrmss,vcrc,Ncenter,Conc) If(X0 > Box .or. Y0 > Box .or. Z0 > Box) & Call HandleError(' --- Catshort :',' Error in upper limit of coordinates:',& iHalo,X0, Y0 , Z0, VvX,VvY,VvZ,aM,rr,vrmss,vcrc,Npart,Conc) If(X0<0. .or. Y0 < 0. .or. Z0<0.) & Call HandleError(' --- Catshort :',' Error in low limit of coordinates:',& iHalo,X0, Y0 , Z0, VvX,VvY,VvZ,aM,rr,vrmss,vcrc,Npart,Conc) If(aM>1.e16.or.aM<1.e7) & Call HandleError(' --- Catshort :',' Error in halo mass :',& iHalo,X0, Y0 , Z0, VvX,VvY,VvZ,aM,rr,vrmss,vcrc,Npart,Conc) If(vcrc<0.or.vcrc>10000.) & Call HandleError(' --- Catshort :',' Error in halo Vcirc :',& iHalo,X0, Y0 , Z0, VvX,VvY,VvZ,aM,rr,vrmss,vcrc,Npart,Conc) !If(aNpart<0.or.aNpart>1.e10) & ! Call HandleError(' --- Catshort :',' Error in number of particles:',& ! iHalo,X0, Y0 , Z0, VvX,VvY,VvZ,aM,rr,vrmss,vcrc,Npart,Conc) If(iDist<0 .or.iDist > 2000000000) & Call HandleError(' --- HaloList Error :',' Error in distinct index:',& iHalo,X0, Y0 , Z0, VvX,VvY,VvZ,aM,rr,vrmss,vcrc,iDist,Conc) If(Ncenter .le. NhaloMax)Then Xc(Ncenter) = X0 Yc(Ncenter) = Y0 Zc(Ncenter) = Z0 VXc(Ncenter) = VvX VYc(Ncenter) = VvY VZc(Ncenter) = VvZ Vcirc(Ncenter) = vcrc rvr = (aMtot/factor)**0.3333*1000. !Vvir = 2.081e-3*sqrt(aMtot/rvr/Aexpn) !if(aMtot>3.5e13)write(*,'(1p,2g12.4,3x,8g12.4)') aM,aMtot,rr,rvr,rvr/rr,Offset,vcrc,Vvir Rmc(Ncenter) = rr !rvr Amc(Ncenter) = aMtot ! aM ! HaloNumber(Ncenter) = iHalo Xoff(Ncenter) = Offset !CAratio(Ncenter) = ca VirRatio(Ncenter) = VirRat aLambda(Ncenter) = aLam !Vrms(Ncenter) = vrmss Rdmax = MAX(rr/1000.,Rdmax) LDistinct(Ncenter) = max(iDist,0) If(iDist == 0)Ndistinct = Ndistinct +1 EndIf EndIf ! test mass Enddo 40 Format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3) close (3) write(*,'(2(a,T39,a,i10/),a,T39,a,1p,g12.4,a)') ' Read Halos','= ',Ncenter,' Ndistinct','= ',Ndistinct, & ' Halo max Radius','= ',Rdmax,'Mpch' write(2,'(2(a,T39,a,i10/),a,T39,a,1p,g12.4,a)') ' Number of Halos','= ',Ncenter,' Ndistinct','= ',Ndistinct, & ' Halo max Radius','= ',Rdmax,'Mpch' If(Ncenter >= NhaloMax)Stop ' Not enough space for halos. Increase NhaloMax.' Return 10 write(*,*) ' Could not read radial bins' stop 20 write(*,*) ' Could not read text R(kpc/h)' stop 30 write(*,*) ' Unexpected End of file' stop End SUBROUTINE ReadHalos !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadROCK(iFlag) !-------------------------------------------------------------- !--- original outXXX.list --- !1 ID !2 DescID !3 Mvir !4 Vmax !5 Vrms 6 Rvir 7 Rs 8 Np 9-X 10-Y 11-Z !12- VX 13-VY 14-VZ 15-JX 16-JY 17-JZ !18- Spin 19-rs_klypin 20-Mvir_all !21- M200b 22-M200c 23-M500c 24-M2500c !25- Xoff 26-Voff 27-spin_bullock 28-b_to_a 29-c_to_a !30- A[x] 31-A[y] 32- A[z] !33- b_to_a(500c) 34- c_to_a(500c) !35- A[x](500c) 36-A[y](500c) 37-A[z](500c) !38- TUratio 39- Mpe_Behr 40- Mpe_Diemer 41- Halmass_Radius 42- rvmax 43-PID !--- Consistent Trees hlistXXX ------ !scale(0) id(1) desc_scale(2) desc_id(3) num_prog(4) pid(5) upid(6) desc_pid(7) phantom(8) !sam_mvir(9) mvir(10) rvir(11) rs(12) vrms(13) mmp?(14) !scale_of_last_MM(15) vmax(16) x(17) y(18) z(19) vx(20) vy(21) vz(22) ! Jx(23) Jy(24) Jz(25) Spin(26) Breadth_first_ID(27) Depth_first_ID(28) ! Rs_Klypin(34) Mmvir_all(35) M200b(36) M200c(37) M500c(38) M2500c(39) ! Xoff(40) Voff(41) Spin_Bullock(42) ! T/|U|(53) M_pe_Behroozi(54) M_pe_Diemer(55) Macc(56) Mpeak(57) ! Vacc(58) Vpeak(59) Halfmass_Scale(60) ! pid(5) -- ID of least massive host halo (-1 if distinct halo) ! Upid(6) -- ID of most massive host halo Character :: file1*80,Line*120,txt*80,cc*1 integer*8 :: id, dum(6), PID real*8 :: dd(60),a32, a33, a34, a35,a36, a37, a38, a39, a40 Ncenter =0 Rdmax =0. factor =1.162e12*Om0*Ovdens Open(3,file='out_126.list') !Open(3,file='hlist_1.00000.list') If(iFlag==1)Then Aexpn =1.0 Box =400. !Open(2,file='HaloStROCK.newMD.z0.Vpeak.DAT') Open(2,file='HaloStROCK.400.126.DAT') EndIf Do Read(3,'(a1)')cc if(cc/='#')exit !write(*,'(a)') cc EndDo Ndistinct = 0 Do ! k =1,800000 ! -------- read halos !read (3,*,iostat=iStat) aexp,id,dscale,dum,(dd(i),i=1,51) ! d1,dum,d2,aM,rvir,rs,vrms, X0, Y0 , Z0, VvX,VvY,VvZ, & ! aM,aMtot,rr,vrmss,vcrc,iHalo,Conc,aNpart, & ! iDist,Offset,VirRat,aLam,Rin, & ! ba,ca,xax,yax,zax read (3,*,iostat=iStat) ih, iDes, d1, vcrc, d2, rvir, rs,aNpart, X0, Y0 , Z0, VvX,VvY,VvZ, & ! 1 2 3 4 5 6 7 8 9 10 11 12 13 14 sx, sy, sz, aLam, rs2, aMtot, a21, a22, a23, a24, Offset, & a26, a27, a28, a29, a30, a31, a32, a33, a34, a35,a36, a37, a38, a39, a40, & a41, a42, PID ! d1,dum,d2,aM,rvir,rs,vrms, X0, Y0 , Z0, VvX,VvY,VvZ, & ! aM,aMtot,rr,vrmss,vcrc,iHalo,Conc,aNpart, & ! iDist,Offset,VirRat,aLam,Rin, & ! ba,ca,xax,yax,zax !read (3,*,iostat=iStat) ii,X0, Y0 , Z0, & ! aM0,aM,vnow,vcrc !rr = aM0**0.33333/4.82e1 ! virial radius in kpch for Om=0.27 delta=360.8 If(mod(Ncenter,100000)==0)& !write (*,'(a,3i12,f8.3,4g13.4,1p,3g12.4)') ' Reading: ',Ncenter,Ndistinct,dum(4),aexp,dd(9),dd(10),dd(11),dd(51),dd(27) write (*,'(a,2i12,g13.4,4f12.4,2es12.3,i12)') ' Reading: ',Ncenter, Ndistinct,aMtot,X0,Y0,Z0,Offset,a41, a42,PID If(iStat.ne.0)then ! backspace(3) ! backspace(3) ! Do m =1,5 ! read(3,*) aexp,id,dscale,dum,(dd(in),in=1,51) ! write(*,'(a,f8.3,i12,f8.3,6i12,51g13.4)') ' ',aexp,id,dscale,dum,(dd(in),in=1,51) ! End Do exit EndIf Ncenter = Ncenter +1 If(Ncenter .le. NhaloMax.and.iFlag==1)Then Xc(Ncenter) = X0 Yc(Ncenter) = Y0 Zc(Ncenter) = Z0 VXc(Ncenter) = Vvx VYc(Ncenter) = Vvy VZc(Ncenter) = Vvz Vcirc(Ncenter) = vcrc ! -Vpeak Amc(Ncenter) = aMtot !aMtot Xoff(Ncenter) = Offset aLambda(Ncenter) = aLam Rmc(Ncenter) = rvir !rvr ! Xc(Ncenter) = dd(9) ! Yc(Ncenter) = dd(10) ! Zc(Ncenter) = dd(11) ! VXc(Ncenter) = dd(12) ! VYc(Ncenter) = dd(13) ! VZc(Ncenter) = dd(14) !Vcirc(Ncenter) = dd(51) ! -Vpeak ! Vcirc(Ncenter) = dd(8) !-Vmax !rvr = (aMtot/factor)**0.3333*1000. !Vvir = 2.081e-3*sqrt(aMtot/rvr/Aexpn) !if(aMtot>3.5e13)write(*,'(1p,2g12.4,3x,8g12.4)') aM,aMtot,rr,rvr,rvr/rr,Offset,vcrc,Vvir ! Rmc(Ncenter) = dd(3) !rvr ! Amc(Ncenter) = dd(27) !aMtot ! HaloNumber(Ncenter) = id ! Xoff(Ncenter) = dd(32) !CAratio(Ncenter) = ca ! VirRatio(Ncenter) = dd(45) ! aLambda(Ncenter) = dd(18) !Vrms(Ncenter) = vrmss Rdmax = MAX(Rmc(Ncenter),Rdmax) If(PID<0)Then LDistinct(Ncenter) = 0 Else LDistinct(Ncenter) = 1 EndIf If(LDistinct(Ncenter) == 0)Ndistinct = Ndistinct +1 !write(*,'(5x,3i9,15g12.4)') Ncenter,Ndistinct,LDistinct(Ncenter),Amc(Ncenter),Vcirc(Ncenter),PID EndIf Enddo close (3) write(*,*) ih, iDes, d1 write(*,*) vcrc, d2, rvir write(*,*) rs,aNpart write(*,'(2(a,T39,a,i10/),a,T39,a,1p,g12.4,a)') ' Read Halos','= ',Ncenter,' Ndistinct','= ',Ndistinct, & ' Halo max Radius','= ',Rdmax,'Mpch' write(2,'(2(a,T39,a,i10/),a,T39,a,1p,g12.4,a)') ' Number of Halos','= ',Ncenter,' Ndistinct','= ',Ndistinct, & ' Halo max Radius','= ',Rdmax,'Mpch' If(Ncenter >= NhaloMax)Stop ' Not enough space for halos. Increase NhaloMax.' Return End SUBROUTINE ReadROCK !--------------------------------------------------------------- ! ! SUBROUTINE HandleError(Text1,Text2,iHalo,X0, Y0 , Z0, & VvX,VvY,VvZ,aM,rr,vrmss,vcrc,Nhalo,Conc) !--------------------------------------------------------------- character(*) :: Text1, Text2 integer :: ListFiles(2) =(/0,20/) Do j=1,2 write (ListFiles(j),'(2a,/3x,a,i9,a,1p,3g12.4,a,3g12.4, & /3x,2(a,g12.4),g12.4,a,0p,i9,a,1p,g12.4)') & Trim(Text1),TRIM(Text2), & ' Line in file=',iHalo,' Coords=', X0, Y0 , Z0, & ' Veloc=',VvX,VvY,VvZ, & ' M=',aM,' R=',rr, & ' Vrms=',vrmss,vcrc, ' Nhalo=',Nhalo,& ' C =',Conc EndDo stop End SUBROUTINE HandleError ! ------------------------ real function Concentration(Vratio) ! ------------------------ ! Real*8 :: A,Aa,C,dC,x,f If(Vratio<1.)Then Concentration = 1. Else C= 2.162 dC =3.d-5 A = Vratio**2 Aa = 1.d0 Do while(Aa tiny ( resabs ) / (5.0e+01* epsilon ( resabs ) ) ) then abserr = max (( epsilon ( resabs ) *5.0e+01)*resabs,abserr) end if return end subroutine qk15 subroutine qk21 ( f, a, b, result, abserr, resabs, resasc ) !*****************************************************************************80 ! !! QK21 carries out a 21 point Gauss-Kronrod quadrature rule. ! ! Discussion: ! ! This routine approximates ! I = integral ( A <= X <= B ) F(X) dx ! with an error estimate, and ! J = integral ( A <= X <= B ) | F(X) | dx ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, external real*8 F, the name of the function routine, of the form ! function f ( x ) ! real*8 f ! real*8 x ! which evaluates the integrand function. ! ! Input, real*8 A, B, the limits of integration. ! ! Output, real*8 RESULT, the estimated value of the integral. ! RESULT is computed by applying the 21-point Kronrod rule (resk) ! obtained by optimal addition of abscissae to the 10-point Gauss ! rule (resg). ! ! Output, real*8 ABSERR, an estimate of | I - RESULT |. ! ! Output, real*8 RESABS, approximation to the integral of the absolute ! value of F. ! ! Output, real*8 RESASC, approximation to the integral | F-I/(B-A) | ! over [A,B]. ! implicit none real*8 a real*8 absc real*8 abserr real*8 b real*8 centr real*8 dhlgth real*8, external :: f real*8 fc real*8 fsum real*8 fval1 real*8 fval2 real*8 fv1(10) real*8 fv2(10) real*8 hlgth integer j integer jtw integer jtwm1 real*8 resabs real*8 resasc real*8 resg real*8 resk real*8 reskh real*8 result real*8 wg(5) real*8 wgk(11) real*8 xgk(11) ! ! the abscissae and weights are given for the interval (-1,1). ! because of symmetry only the positive abscissae and their ! corresponding weights are given. ! ! xgk - abscissae of the 21-point Kronrod rule ! xgk(2), xgk(4), ... abscissae of the 10-point ! Gauss rule ! xgk(1), xgk(3), ... abscissae which are optimally ! added to the 10-point Gauss rule ! ! wgk - weights of the 21-point Kronrod rule ! ! wg - weights of the 10-point Gauss rule ! data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7),xgk(8), & xgk(9),xgk(10),xgk(11)/ & 9.956571630258081e-01, 9.739065285171717e-01, & 9.301574913557082e-01, 8.650633666889845e-01, & 7.808177265864169e-01, 6.794095682990244e-01, & 5.627571346686047e-01, 4.333953941292472e-01, & 2.943928627014602e-01, 1.488743389816312e-01, & 0.000000000000000e+00/ ! data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7),wgk(8), & wgk(9),wgk(10),wgk(11)/ & 1.169463886737187e-02, 3.255816230796473e-02, & 5.475589657435200e-02, 7.503967481091995e-02, & 9.312545458369761e-02, 1.093871588022976e-01, & 1.234919762620659e-01, 1.347092173114733e-01, & 1.427759385770601e-01, 1.477391049013385e-01, & 1.494455540029169e-01/ ! data wg(1),wg(2),wg(3),wg(4),wg(5)/ & 6.667134430868814e-02, 1.494513491505806e-01, & 2.190863625159820e-01, 2.692667193099964e-01, & 2.955242247147529e-01/ ! ! ! list of major variables ! ! centr - mid point of the interval ! hlgth - half-length of the interval ! absc - abscissa ! fval* - function value ! resg - result of the 10-point Gauss formula ! resk - result of the 21-point Kronrod formula ! reskh - approximation to the mean value of f over (a,b), ! i.e. to i/(b-a) ! centr = 5.0e-01*(a+b) hlgth = 5.0e-01*(b-a) dhlgth = abs(hlgth) ! ! Compute the 21-point Kronrod approximation to the ! integral, and estimate the absolute error. ! resg = 0.0e+00 fc = f(centr) resk = wgk(11)*fc resabs = abs(resk) do j = 1, 5 jtw = 2*j absc = hlgth*xgk(jtw) fval1 = f(centr-absc) fval2 = f(centr+absc) fv1(jtw) = fval1 fv2(jtw) = fval2 fsum = fval1+fval2 resg = resg+wg(j)*fsum resk = resk+wgk(jtw)*fsum resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2)) end do do j = 1, 5 jtwm1 = 2*j-1 absc = hlgth*xgk(jtwm1) fval1 = f(centr-absc) fval2 = f(centr+absc) fv1(jtwm1) = fval1 fv2(jtwm1) = fval2 fsum = fval1+fval2 resk = resk+wgk(jtwm1)*fsum resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2)) end do reskh = resk*5.0e-01 resasc = wgk(11)*abs(fc-reskh) do j = 1, 10 resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) end do result = resk*hlgth resabs = resabs*dhlgth resasc = resasc*dhlgth abserr = abs((resk-resg)*hlgth) if ( resasc /= 0.0e+00.and.abserr /= 0.0e+00) then abserr = resasc*min ( 1.0e+00,(2.0e+02*abserr/resasc)**1.5e+00) end if if ( resabs > tiny ( resabs ) /(5.0e+01* epsilon ( resabs ) )) then abserr = max (( epsilon ( resabs ) *5.0e+01)*resabs,abserr) end if return end subroutine qk21 subroutine qsort ( limit, last, maxerr, ermax, elist, iord, nrmax ) !*****************************************************************************80 ! !! QSORT maintains the order of a list of local error estimates. ! ! Discussion: ! ! This routine maintains the descending ordering in the list of the ! local error estimates resulting from the interval subdivision process. ! At each call two error estimates are inserted using the sequential ! search top-down for the largest error estimate and bottom-up for the ! smallest error estimate. ! ! Author: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner ! ! Reference: ! ! Robert Piessens, Elise de Doncker-Kapenger, ! Christian Ueberhuber, David Kahaner, ! QUADPACK, a Subroutine Package for Automatic Integration, ! Springer Verlag, 1983 ! ! Parameters: ! ! Input, integer LIMIT, the maximum number of error estimates the list can ! contain. ! ! Input, integer LAST, the current number of error estimates. ! ! Input/output, integer MAXERR, the index in the list of the NRMAX-th ! largest error. ! ! Output, real*8 ERMAX, the NRMAX-th largest error = ELIST(MAXERR). ! ! Input, real*8 ELIST(LIMIT), contains the error estimates. ! ! Input/output, integer IORD(LAST). The first K elements contain ! pointers to the error estimates such that ELIST(IORD(1)) through ! ELIST(IORD(K)) form a decreasing sequence, with ! K = LAST ! if ! LAST <= (LIMIT/2+2), ! and otherwise ! K = LIMIT+1-LAST. ! ! Input/output, integer NRMAX. ! implicit none integer last real*8 elist(last) real*8 ermax real*8 errmax real*8 errmin integer i integer ibeg integer iord(last) integer isucc integer j integer jbnd integer jupbn integer k integer limit integer maxerr integer nrmax ! ! Check whether the list contains more than two error estimates. ! if ( last <= 2 ) then iord(1) = 1 iord(2) = 2 go to 90 end if ! ! This part of the routine is only executed if, due to a ! difficult integrand, subdivision increased the error ! estimate. in the normal case the insert procedure should ! start after the nrmax-th largest error estimate. ! errmax = elist(maxerr) do i = 1, nrmax-1 isucc = iord(nrmax-1) if ( errmax <= elist(isucc) ) then exit end if iord(nrmax) = isucc nrmax = nrmax-1 end do ! ! Compute the number of elements in the list to be maintained ! in descending order. This number depends on the number of ! subdivisions still allowed. ! jupbn = last if ( (limit/2+2) < last ) then jupbn = limit+3-last end if errmin = elist(last) ! ! Insert errmax by traversing the list top-down, starting ! comparison from the element elist(iord(nrmax+1)). ! jbnd = jupbn-1 ibeg = nrmax+1 do i = ibeg, jbnd isucc = iord(i) if ( elist(isucc) <= errmax ) then go to 60 end if iord(i-1) = isucc end do iord(jbnd) = maxerr iord(jupbn) = last go to 90 ! ! Insert errmin by traversing the list bottom-up. ! 60 continue iord(i-1) = maxerr k = jbnd do j = i, jbnd isucc = iord(k) if ( errmin < elist(isucc) ) then go to 80 end if iord(k+1) = isucc k = k-1 end do iord(i) = last go to 90 80 continue iord(k+1) = last ! ! Set maxerr and ermax. ! 90 continue maxerr = iord(nrmax) ermax = elist(maxerr) return end subroutine qsort !---------------------------------------------------- Real*8 FUNCTION Hgrow(x) Real*8 :: x Hgrow =(sqrt(x/(1.+x**3)))**3 End FUNCTION Hgrow !-------------------------------------- P*k^2*Top-Hat Filter REAL*8 FUNCTION Ptophat(x) !--------------------------------------- use setHalos Real*8 :: wk,x,TOPHAT,Ppk wk = x/rf TOPHAT =( (SIN(x)-x*COS(x))*3.d0/x**3 )**2 ! write(*,*) ' ',wk,TOPHAT Ptophat=Ppk(wk)*x**2*TOPHAT END FUNCTION Ptophat