!-------------------------------------------------------------------- ! Fit halo profiles, ! use simple minimization procedure !-------------------------------------------------------------------- Module Structures IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4, PARAMETER :: Nmax = 150000, & ! Nbinmax = 30, & ! Nqual = 23 ! Real*8, PARAMETER :: factor = 8.23d10 Integer*4 :: Nhalos Real*8 :: Aexpn,Amfit,slope Integer*4, DIMENSION(Nmax) :: Nhbin Real*8, DIMENSION(Nmax) :: rsh,Vmaxh, & rho1,errorh,slopeh Real*8, DIMENSION(Nqual,Nmax) :: hals Real*8, DIMENSION(17,Nbinmax,Nmax) :: Prof !$OMP THREADPRIVATE(slope) Contains !-------------------------------------------------------------------- ! SUBROUTINE ReadHalos ! !-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Character :: Header*100,Txt1*2,Txt2*2 Lines =15 ! number of lines in the header preceeding radial bins lbins =100 ! number of bins along radius Lskip =5 ! number of header lines between bins and catalog aMfit = 1.e12 ! limit for mass for fitting NFW or Sersic profiles read(1,'(a)') Header write (*,'(a)') Header write (2,'(a)') Header read(1,*) Txt1,Txt2,Aexpn write (*,'(2a,f8.4)') Txt1,Txt2,Aexpn read(1,*)Txt1,Txt2,ISTEP write (*,'(2a,i4)') Txt1,Txt2,ISTEP Do i=1,Lines read(1,'(a)') Header write (*,'(a)') Header if(i.le.Lines-2)write (2,'(a)') Header EndDo !stop write(2,'(T10,a,T36,a,T51,a,T62,a,T69,a,T78,a,T88,a,T98,a,T108,a,T116,a)') & 'XYZ(Mpch)','Vxyz','Mvir','Rvir', & 'V(fit)','Vcirc','C(fit)','C','Error', & 'Xoff 2K/P-1 Lambda' Nh =0 Do ! js =1,50001 ! read catalog Nh = Nh +1 If(Nh.gt.Nmax)Stop 'Not enough space for Halos' read(1,*,iostat=ierr) (hals(i,Nh),i=1,Nqual),Nhbin(Nh) If(ierr /= 0)exit !write(*,'(" Halo=",i8," coords=",3f8.3," bins=",i3)') & ! Nh,(hals(i,Nh),i=1,3),Nhbin(Nh) do i=1,Nhbin(Nh) read(1,*) (Prof(j,i,Nh),j=1,17) EndDo End Do Nhalos =Nh -1 write(*,*) ' Number of halos = ',Nhalos end SUBROUTINE ReadHalos !-------------------------------------------------------------------- ! Fit halo Einast=Sersic profile SUBROUTINE HaloFitSer(ih) ! output: Rs = estimate of NFW core radius, ! Vvmax = estimate of Vmax ! slope = 1/n slope for Sersic fit !-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Integer*4, parameter :: nStepS = 50 ! number of bins in slope Integer*4, parameter :: nStepR = 100 ! number of bins in radius Real*8, parameter :: facR = 2.0, & ! max variation for rs facRvir = 1.0, & ! max radius for fit in units of Rvir fourpi = 12.56637d0 Real*8, DIMENSION(Nbinmax) :: Radden, Dens,Vc Real*8, DIMENSION(Nbinmax) :: Radout Integer*4, DIMENSION(Nbinmax) :: Nin Real*8 :: Mvir,aMass,rss,Rvir DATA iFlag/0/ SAVE iFlag Mvir = hals(7,ih) Umax = hals(10,ih) Rumax= hals(8,ih)/hals(12,ih)*2.15 Rvir = hals(8,ih) Nbin = Nhbin(ih) do i=1,Nbin Radden(i) = Prof(3,i,ih)*Rvir Nin(i) = Prof(2,i,ih) Vc(i) = Prof(5,i,ih) Dens(i) = Prof(6,i,ih) Radout(i) = Prof(1,i,ih) EndDo Rs = Rumax/2.15 ! initial guess for rs VVmax = Umax Vcmax = Umax Rmin = rs/facR dR = rs/nStepR*(facR -1./facR) ! step in rs Smin = 0.07 Smax = 0.250 dS = (Smax-Smin)/nStepS ! step in slope iFlag = iFlag +1 !write(*,'(1p,6g12.5)') Mvir,Rvir,Rs,Umax slope = Smin error = 1.e10 istart= 1 Do istart =1,Nbin If(Radout(istart)>50.)exit EndDo ! write (*,*)'not enough particles(first,last=)',Nin(1),Nin(Nbin) !return ! not enough particles; get back 20 Rmax = min(Rvir/facRvir,Radout(Nbin)) ! max radius for the fit do i=Nbin,1,-1 if(Radout(i).lt.Rmax)goto 30 EndDo 30 imax = i if(imax.le.istart)Then write (*,*)'not enough particles. Bins=',istart,imax return EndIf Do is =0,nStepS slope = Smin +is*dS Do ir =0,nStepR rss = Rmin +ir*dR aa = aMass(Rvir/rss) rho_0 = Mvir/(fourpi*rss**3*aa) ! write(*,'(a,1p,4g12.5)') 'Slope,Rs=',slope,rss,rho_0 ! write(*,'(a,1p,4g12.5)') ' =',Mvir,fourpi,aa ee =0. ww =0. do i=istart,imax w = 1. !(Radout(i)-Radout(i-1))/Radden(i) ee=ee+w*(log10(RhoSER(Radden(i),rho_0,rss)/Dens(i)))**2 ww = ww + w EndDo ! ee =ee/(imax-istart) ee =ee/ww If(ee.lt.error)then error = ee Rs = rss Slp = slope irr = ir rho = rho_0 !write (*,'(10x," r/slope=",2f8.3," err=",1p,3g11.3)') & ! rss,Slp,ee,rho_0 EndIf EndDo EndDo slope = slp rho_0 = rho xmax = 2.603/exp(1.2*slope -1.3*slope**2+0.4*slope**4) Conc = Rvir/Rs VVmax = 2.076e-3*sqrt(Mvir/Rs* & xmax**2 /exp(2./slope*(xmax**slope-1.)) /aMass(Conc)) !write (*,'(10x,"V=",f8.3," r/Slp=",2f8.3," err=",1p,7g12.4)') & ! VVmax,Rs,Slp,error,rho,Rvir,Rvir/Rs,Conc,xmax !write (3,'(10x,"V=",f8.3," r/Slp=",2f8.3," err=",1p,7g12.4)') & ! VVmax,Rs,Slp,error,rho,Rvir,Rvir/Rs,Conc,xmax slope = Slp 50 format(i3,' r=',f8.3,' dens=',2g11.4,' v=',3f8.2) !Do i=1,Nbin ! write (3,'(1p,4g11.4,0p,i9)')Radout(i),Dens(i),Vc(i), & ! RhoSER(Radden(i),rho_0,Rs),Nin(i) !EndDo slopeh(ih) = slope rsh(ih) = Rs errorh(ih) = error rho1(ih) = rho Vmaxh(ih) = VVmax ! if(iFlag.ge.10)stop End SUBROUTINE HaloFitSer !--------------------------------------------------------- ! Sersic density profile Real*8 FUNCTION RhoSER(r,rho_0,rs) ! IMPLICIT REAL*8 (A-H,O-Z) x = r/rs RhoSER = rho_0/exp(2./slope*(x**slope-1.)) End FUNCTION RhoSER !---------------------------------------------------------- ! mass profile for Einasto density Real*8 Function aMass(x) IMPLICIT REAL*8 (A-H,O-Z) Real*8, parameter :: x0 =1.d-6 Real*8 :: x,xup,ANS,ER,INTG External Faux xup = max(x,x0) !ANS = INTG(Faux,x0,xup) Call DGAUS8 (Faux, x0, xup, ER, ANS, IERR) aMass = ANS !write(*,'(a,1p,2g12.4)') ' aMass out =',aMass,ANS End Function aMass !-------------------------------------------------------------------- ! Fit NFW halo profile SUBROUTINE HaloFit(ih,iFlag) ! output: Rs = estimate of NFW core radius, ! Vvmax = estimate of Vmax ! iFlag = 0 - for two-parameter fit ! = 1 - fixed Mvir !-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) Real*4, parameter :: facV = 1.25 ! max variation for Vmax Real*4, parameter :: facR = 4.0 ! max variation for rs Real*8, PARAMETER :: facRvir = 1.0 ! max radius for fit in units of Rvir Integer*4, parameter :: nStepR = 450 ! number of bins in radius Integer*4, parameter :: nStepV = 250 ! number of bins in Vmax Real*8, DIMENSION(Nbinmax) :: Radden, Dens,Vc Real*8, DIMENSION(Nbinmax) :: Radout Integer*4, DIMENSION(Nbinmax) :: Nin Real*4 :: Mvir Umax = hals(10,ih) Rumax= hals(8,ih)/hals(12,ih)*2.15 Rvir = hals(8,ih) Mvir = hals(7,ih) Nbin = Nhbin(ih) do i=1,Nbin Radden(i) = Prof(11,i,ih)*Rvir !Prof(3,i,ih)*Rvir Nin(i) = Prof(10,i,ih) !Prof(2,i,ih) Vc(i) = Prof(13,i,ih) !Prof(5,i,ih) Dens(i) = Prof(14,i,ih) !Prof(6,i,ih) Radout(i) = Prof(1,i,ih) !Prof(1,i,ih) EndDo rs = Rumax/2.15 ! initial guess for rs VVmax = Umax Vcmax = Umax Vmin = Umax/facV dV = Umax/nStepV*(facV -1./facV) ! step in Vmax Rmin = rs/facR dR = rs/nStepR*(facR -1./facR) ! step in rs write(*,'(2i7,1p,6g12.4)') ih,Nbin,Umax,Rumax,Rvir,dV,dR error = 1.e10 istart= 1 Do istart =1,Nbin If(Radout(istart).gt.20.)exit EndDo !return ! not enough particles; get back 20 Rmax = min(Rvir/facRvir,Radout(Nbin)) ! max radius for the fit do i=Nbin,1,-1 if(Radout(i).lt.Rmax)goto 30 EndDo 30 imax = i write(*,*) ' Rmax=',Rmax,imax if(imax.le.istart)Then write (*,*)'not enough particles. Bins=',istart,imax return EndIf If(iFlag == 0) Then ! -------------- two-parameter fit Do iv =0,nStepV Vmaxx = Vmin +iv*dV Do ir =0,nStepR rss = Rmin +ir*dR rho_0 = 8.539e4*(Vmaxx/rss)**2*Aexpn ee =0. ww =0. do i=istart,imax w = Nin(i) !(Radout(i)-Radout(i-1))/Radden(i) ee = ee + w*(log10(RhoNFW(Radden(i),rho_0,rss)/Dens(i)))**2 ww = ww + w EndDo ee =ee/ww/(imax-istart+1) If(ee.lt.error)then error = ee Rs = rss VVmax = Vmaxx rho = rho_0 !write (*,'(10x,"New: V=",f8.3," r=",f8.3," err=",3g11.3)') & ! Vmaxx,rss,ee,rho_0 EndIf EndDo !write (*,'(10x,"V=",f8.3," r=",f8.3," err=",3g11.3)') & ! Vmaxx,rss,ee,rho_0 EndDo Else ! ----------------- one-parameter. Fix Virial mass Do ir =0,nStepR rss = Rmin +ir*dR Conc = Rvir/rss Fc = log(1.+Conc) -Conc/(1.+Conc) Vmaxx = 2.076e-3*sqrt(Mvir/(rss*Aexpn)*0.216/Fc) rho_0 = 8.539e4*(Vmaxx/rss)**2*Aexpn ee =0. ww =0. do i=istart,imax w = Nin(i) !1. !(Radout(i)-Radout(i-1))/Radden(i) ee = ee + w*(log10(RhoNFW(Radden(i),rho_0,rss)/Dens(i)))**2 ww = ww + w EndDo ee =ee/ww/(imax-istart+1) If(ee.lt.error)then error = ee Rs = rss VVmax = Vmaxx rho = rho_0 !write (*,'(10x,"New: V=",f8.3," r=",f8.3," err=",3g11.3)') & ! Vmaxx,rss,ee,rho_0 EndIf EndDo end If ! -------------- end iFlag 50 format(i3,' r=',f8.3,' dens=',2g11.4,' v=',3f8.2) rss = hals(8,ih)/hals(12,ih) r0 = 8.539e4*(Umax/rss)**2*Aexpn Do i=1,Nbin xx = Radout(i)/rss rho_vratio = r0/xx/(1.+xx)**2 write (*,'(1p,5g11.4,0p,i9)')Radout(i),Dens(i),Vc(i), & RhoNFW(Radden(i),rho,Rs),rho_vratio,Nin(i) EndDo rsh(ih) = Rs errorh(ih) = error rho1(ih) = rho Vmaxh(ih) = VVmax Conc = Rvir/Rs Fc = log(1.+Conc) -Conc/(1.+Conc) aMvir = 1.073e6* VVmax**2*Rs*Aexpn*Fc Conc = Rvir/rss Fc = log(1.+Conc) -Conc/(1.+Conc) aMvir_vrat = 1.073e6* Umax**2*rss*Aexpn*Fc write(*,'(2(a,2f8.3),a,1p,3g12.4)') 'Rs=',Rs,hals(8,ih)/hals(12,ih),' Vmax=',VVmax,hals(10,ih), & ' Mvir =',hals(7,ih),aMvir, aMvir_vrat End SUBROUTINE HaloFit !-------------------------------------------------------------------- ! NFW density profile Real*8 FUNCTION RhoNFW(r,rho_0,rs) ! output: Rs = estimate of NFW core radius, ! Vvmax = estimate of Vmax !-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) x = r/rs RhoNFW = rho_0/x/(1.+x)**2 End FUNCTION RhoNFW end Module Structures !-------------------------------------------------------------------- ! Program HaloProfiles use Structures IMPLICIT REAL*8 (A-H,O-Z) Open(1,file='CatalogASCIIW.0120.DAT') Open(2,file='CatFitNFW.0120.dat') Call ReadHalos Nchunk = max(Nhalos/20,100) !.... Open_MP !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ih,Umax,Rumax,Rvir,Nbin,Radden) & !$OMP PRIVATE (Nin,Vc,Dens) Do ih=1,Nhalos If(mod(ih,100).eq.0)write(*,*) ' halo=',ih ! CALL HaloFitSER(ih) CALL HaloFit(ih,1) EndDo Do ih=1,Nhalos Rvir = hals(8,ih) Conc = Rvir/Rsh(ih) errorh(ih) =min(errorh(ih),1.d0) Do i=1,Nhbin(ih) If(Prof(1,i,ih).ge.Rvir)Then nn =INT(Prof(2,i,ih)) nL =INT(Prof(2,i,ih)) goto 18 EndIf EndDo 18 write(2,20) (hals(i,ih),i=1,8), & Vmaxh(ih),hals(10,ih),Conc,hals(12,ih), & sqrt(errorh(ih)), (hals(i,ih),i=15,17) ! & ,r1(ih),rho1(ih),x2c(ih),a2c(ih) 20 format(3f8.3,3F8.1,g11.3,f8.2,2f8.1,3x,2f8.3, & 3x,6f8.4) enddo 50 format(f8.3,i7,2g11.4,2f9.2) stop end Program HaloProfiles !---------------------------------------------------------- ! integrant for Einasto=Sersic mass profile Real*8 Function Faux(x) use Structures IMPLICIT REAL*8 (A-H,O-Z) real*8 :: x Faux = exp(-2.d0/slope*(x**slope -1.))*x**2 !Faux = exp(-x) !write(*,*) x,Faux,slope End Function Faux !-------------------------------------------- Simpson integration REAL*8 FUNCTION INTG(FUNC,A,B) !--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (EPS=1.0d-4, JMAX=22) EXTERNAL FUNC OST=-1.D30 OS= -1.D30 ST =0. DO J=1,JMAX CALL TRAPZD(FUNC,A,B,ST,J) INTG=(4.0d0*ST-OST)/3.0d0 IF (ABS(INTG-OS).Le.EPS*ABS(OS)) RETURN OS=INTG OST=ST end DO WRITE (*,*)'Integration did not converge' RETURN END FUNCTION INTG !---------------------------------------------- SUBROUTINE TRAPZD(FUNCC,A,B,S,N) !--------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) SAVE IT EXTERNAL FUNCC IF (N.EQ.1) THEN S=0.5d0*(B-A)*(FUNCC(A)+FUNCC(B)) IT=1 ELSE TNM=IT DEL=(B-A)/TNM X=A+0.5D0*DEL SUM=0.0D0 DO J=1,IT SUM=SUM+FUNCC(X) X=X+DEL END DO S=0.5D0*(S+(B-A)*SUM/TNM) IT=2*IT ENDIF RETURN END SUBROUTINE TRAPZD