c-------------------------------------------------------------------- c Fit halo profiles, c use simple minimization c-------------------------------------------------------------------- PARAMETER (Nmax = 200000) ! PARAMETER (Nbinmax = 85) ! PARAMETER (Nqual = 13) ! COMMON /HALOS/hals(Nqual,Nmax),Nhbin(Nmax), & Prof(7,Nbinmax,Nmax) COMMON /Binn/ Radout(Nbinmax) COMMON /APPROX/rsh(Nmax),Vmaxh(Nmax) DIMENSION Nin(Nbinmax),Radden(Nbinmax),Vc(Nbinmax),Dens(Nbinmax) DIMENSION VVmax(Nmax),Rs(Nmax),r1(Nmax),rho1(Nmax),x2c(Nmax), & a2c(Nmax),error(Nmax) Character Header*100 Open(1,file='CatalogA.999.DAT') Open(2,file='CatFitAf.999.dat') Open(3,file='CatOut.999.dat') Lines =14 ! number of lines in the header preceeding radial bins lbins =85 ! number of bins along radius Lskip =4 ! number of header lines between bins and catalog Do i=1,Lines read(1,'(a)') Header write (*,'(a)') Header write (2,'(a)') Header EndDo read(1,*) (Radout(i),i=1,lbins) write(2,'(10f8.2)') (Radout(i),i=1,lbins) write (*,'(10f8.2)') (Radout(i),i=1,lbins) Do i=1,Lskip read(1,'(a)') Header write (2,'(a)') Header write (*,'(a)') Header EndDo Nh =0 100 Nh =Nh+1 read(1,*,err=400,end=400) (hals(i,Nh),i=1,Nqual),Nhbin(Nh) c write(*,'(" Halo=",i8," coords=",3f8.3," bins=",i3)') c & Nh,(hals(i,Nh),i=1,3),Nhbin(Nh) do i=1,Nhbin(Nh) read(1,*) (Prof(j,i,Nh),j=1,7) EndDo If(Nh.lt.Nmax)goto100 c If(Nh.lt.1000)goto100 Nh =Nh +1 400 Nh =Nh -1 write (*,*) ' Number of halos read =',Nh factor = 8.23e10 Nchunk = max(Nh/20,100) C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (ih,Umax,Rumax,Rvir,Nbin,Radden) C$OMP+PRIVATE (Nin,Vc,Dens) C$OMP+SCHEDULE(DYNAMIC,Nchunk) Do ih=1,Nh If(mod(ih,Nchunk/2).eq.0)write(*,*) ' halo=',ih Umax = hals(10,ih) Rumax= hals(12,ih) Rvir = hals(8,ih) Nbin = Nhbin(ih) do i=1,Nbin Radden(i) = Prof(1,i,ih) Nin(i) = Prof(2,i,ih) Vc(i) = Prof(6,i,ih) Dens(i) = Prof(7,i,ih)/factor EndDo if(hals(7,ih).lt.1.e12)Then CALL HaloFit(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs(ih),VVmax(ih),Nbin,error(ih)) r1(ih) = 0. rho1(ih) =0. x2c(ih) =0. a2c(ih) =0. Else CALL HaloFitExp(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs(ih),VVmax(ih),r1(ih),rho1(ih),x2c(ih), & a2c(ih),Nbin,error(ih)) EndIf EndDo Do ih=1,Nh Rvir = hals(8,ih) Conc = Rvir/Rs(ih) error(ih) =min(error(ih),1.e0) write(2,20) (hals(i,ih),i=1,9), & VVmax(ih),INT(hals(11,ih)),Rs(ih)*2.15,Conc,Nhbin(ih), & sqrt(error(ih)) & ,r1(ih),rho1(ih),x2c(ih),a2c(ih) 20 format(3f8.3,3F8.1,g11.3,f8.2,2f8.1,i9,f8.1,f8.2,i4,f8.3 & f7.3,g11.4,f7.2,f7.4) If((Conc.lt.4..or.Conc.gt.25.).and.hals(7,ih).gt.1.e12)Then write(3,'(" Halo=",i8," coords=",3f8.3," M=",g11.4, & " C=",f8.3," R1,Rs,Rvir=",3f8.3)') & ih,(hals(i,ih),i=1,3),hals(7,ih), & Conc,r1(ih),Rs(ih),hals(8,ih) write (3,'(5x,"r Nin Dens Dens_4exp", & 6x,"Vc Vc_4exp")') do i=1,Nhbin(ih) Radd = Prof(1,i,ih) Ni = Prof(2,i,ih) Vcc = Prof(6,i,ih) Den = Prof(7,i,ih)/factor write(3,50) Radd,Ni,Den, & RhoExp(Radd,r1(ih),rho1(ih),x2c(ih),a2c(ih)), & Vcc,VcExp(Radd,r1(ih),rho1(ih),x2c(ih),a2c(ih)) EndDo EndIf enddo 50 format(f8.3,i7,2g11.4,2f9.2) stop end c-------------------------------------------------------------------- c Fit halo profile SUBROUTINE HaloFit(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs,VVmax,Nbin,error) c output: Rs = estimate of NFW core radius, c Vvmax = estimate of Vmax c-------------------------------------------------------------------- PARAMETER (Nbinmax = 85) ! PARAMETER (facV = 1.25) ! max variation for Vmax PARAMETER (facR = 3.0) ! max variation for rs c PARAMETER (facRvir = 2.0) ! max radius for fit in units of Rvir PARAMETER (facRvir = 1.0) ! max radius for fit in units of Rvir PARAMETER (nStepR = 50) ! number of bins in radius PARAMETER (nStepV = 20) ! number of bins in Vmax DIMENSION Radden(Nbin),Nin(Nbin),Dens(Nbin),Vc(Nbin), & Radout(Nbinmax) 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 error = 1.e10 istart= 1 Do istart =1,Nbin If(Nin(istart).gt.200)goto20 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)goto30 EndDo 30 imax = i if(imax.le.istart)Then c write (*,*)'not enough particles. Bins=',istart,imax return EndIf Do iv =0,nStepV Vmaxx = Vmin +iv*dV Do ir =0,nStepR rss = Rmin +ir*dR rho_0 = 1.022e3*(Vmaxx/rss)**2 ee =0. ww =0. do i=istart,imax w = (Radout(i)-Radout(i-1))/Radden(i) ee = ee + w*(log10(RhoNFW(Radden(i),rho_0,rss)/Dens(i)))**2 ww = ww + w EndDo c ee =ee/(imax-istart) ee =ee/ww If(ee.lt.error)then error = ee Rs = rss VVmax = Vmaxx ivv = iv irr = ir rho = rho_0 c write (*,'(10x,"V=",f8.3," r=",f8.3," err=",3g11.3)') c & Vmaxx,rss,ee,rho_0 EndIf EndDo EndDo c write (*,'(10x,"V=",f8.3," r=",f8.3," err=",3g11.3)') c & VVmax,Rs,error,rho c write(*,*) ' after ----' c do i=1,Nbin c x = Radden(i)/rs c F = log(1.+x)-x/(1.+x) c write(*,50) i,Radden(i),Dens(i),rho/x/(1.+x)**2, c & Vc(i),VVmax*sqrt(4.62*F/x) c enddo 50 format(i3,' r=',f8.3,' dens=',2g11.4,' v=',3f8.2) Return End c-------------------------------------------------------------------- c NFW density profile FUNCTION RhoNFW(r,rho_0,rs) c output: Rs = estimate of NFW core radius, c Vvmax = estimate of Vmax c-------------------------------------------------------------------- x = r/rs RhoNFW = rho_0/x/(1.+x)**2 Return End c-------------------------------------------------------------------- c Double Exp-Fit for halo profile SUBROUTINE HaloFitExp(Umax,Rumax,Rvir,Radden,Radout,Nin,Dens,Vc, & Rs,VVmax,r1,rho1,x2c,a2c,Nbin,error) c output: Rs = estimate of NFW core radius, c Vvmax = estimate of Vmax c Uses 4 parameters: c rho_1, r_1, a_2, x_2 c-------------------------------------------------------------------- PARAMETER (Nbinmax = 85) ! PARAMETER (facRho = 10.) ! max variation for density PARAMETER (facR = 10.0) ! max variation for r1 c PARAMETER (facRvir = 2.0) ! max radius for fit in units of Rvir PARAMETER (facRvir = 0.5) ! max radius for fit in units of Rvir DIMENSION Radden(Nbin),Nin(Nbin),Dens(Nbin),Vc(Nbin), & Radout(Nbinmax) Logical Converge rs = Rumax/2.15 ! initial guess for rs VVmax = Umax r1 = min(max(Rumax/2.15,Radden(2)),Rvir/5.) ! initial guess for r1 dR = r1/20. rho1 = 1.022e3*(Umax/r1)**2 ! initial guess for rho_1 dRho = rho1/10. x2c = 5. a2c = 1./10. dx2 = x2c/10. da2 = a2c/10. error = 100. c write (*,*) c write(*,'(" Guess: R,rho,a2,x2=",4g11.3)')r1,rho1,a2c,x2c c write(*,'(" steps: R,rho,a2,x2=",4g11.3)')dR,dRho,da2,dx2 istart= 1 Do istart =1,Nbin If(Nin(istart).gt.100)goto20 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)goto30 EndDo 30 imax = i c write (*,*) ' max bin=',imax,' Rout=',Radout(imax),Rmax if(imax.le.istart.or.Nin(imax).lt.1000)Then c write (*,*)'not enough particles. Bins=',istart,imax return EndIf Do iter =1,100 r1_old = r1 rho1_old = rho1 a2c_old = a2c x2c_old = x2c e_old = error c write (*,*) ' ==== iteration=== ',iter,error CALL LocalFit(Radden,Radout,Dens,istart,imax,Nbin, & r1,rho1,a2c,x2c,dR,dRho,dx2,da2,error,iFlag) c write (*,'(10x,"r1=",f8.4," rho1=",7g11.4," err=",g11.4,i2)') c & r1,rho1,a2c,x2c,dR,dRho,dx2,da2,error,iFlag quality = (r1/r1_old-1.)**2+(rho1/rho1_old-1.)**2+ & (a2c/a2c_old-1)**2+(x2c/x2c_old-1)**2 Converge = (quality.lt.1.e-4).and. & (abs(error/e_old-1.).lt.1.e-3).and. & (iter.gt.10).and.(iFlag.eq.0) dRho =max(min(dRho,rho1/10.),rho1/30.) If(Converge)Goto40 If(error.ge.e_old.and.iter.gt.10)Then r1 = r1_old rho1= rho1_old a2c = a2c_old x2c = x2c_old dR =dR/1.21 dRho =dRho/1.21 dx2 =dx2/1.21 da2 =da2/1.15 c write (*,*) ' restore center, reduce steps' EndIf If(iFlag.eq.0)Then dR =min(dR/1.41,r1/20.) dRho =min(dRho/1.41,rho1/20.) dx2 =dx2/1.41 da2 =min(da2/1.25,a2c/20.) c write (*,*) ' changed steps' c Else c dR =dr/1.01 c dRho =drho/1.01 EndIf EndDo 40 Rs =r_1 c write (*,'(10x,"r1=",f8.3," rho1=",7g11.4," err=",g11.3,2i4)') c & r1,rho1,a2c,x2c,dR,dRho,dx2,da2,error,iFlag,iter CALL FindVcMax(rmax,Vcmax,Rvir,r_s,r1,rho1,x2c,a2c,0.05) Rout = Radden(imax) rhs = RhoExp(Rout,r1,rho1,x2c,a2c)/4. CALL FindRs(r_s2,r1,rho1,x2c,a2c,0.5,rhs,Rout) r_s =(r_s+r_s2)/2. c r_s =min(r_s,r_s2) c r_s =r_s2 c If(Rvir/r_s.lt.4.)Then c write (*,'(" rmax,Vmax=",2f8.3," rs=",3g11.3)')rmax,Vcmax, c & r_s,Rout,Dens(imax) c do i=istart,imax c write(*,50) i,Radden(i),Dens(i), c & RhoExp(Radden(i),r1,rho1,x2c,a2c), c & rhs*Rout/Radden(i)*(1.+Rout/Radden(i))**2, c & Vc(i),VcExp(Radden(i),r1,rho1,x2c,a2c) c enddo c EndIf c 50 format(i3,' r= ',f8.3,' dens=',3g11.4,' v=',3f8.2) Rs = r_s VVmax = Vcmax Return End c-------------------------------------------------------------------- c Fit halo profile SUBROUTINE LocalFit(Radden,Radout,Dens,istart,imax,Nbin, & Rc,Rhoc,a2c,x2c,dR,dRho,dx2,da2,error,iFlag) PARAMETER (Nbinmax = 85) ! PARAMETER (nStepR = 6) ! number of bins in radius PARAMETER (nStepRho= 6) ! number of bins in Rho PARAMETER (nStepa2 = 6) ! number of bins in a2 PARAMETER (nStepx2 = 6) ! number of bins in x2 DIMENSION Radden(Nbin),Dens(Nbin),Radout(Nbinmax) Logical on_boundary Rmin = max(Rc -dR*nStepR/2.,Rc/10.) Rhomin = max(Rhoc -dRho*nStepRho/2.,Rhoc/10.) x2min = max(x2c - dx2*nStepx2/2.,1.e-2) a2min = max(a2c - da2*nStepa2/2.,a2c/10.) c write (*,'(5x," Rmin/max =",4g11.4)') c & Rmin,Rmin +nStepR*dR,Rc c write (*,'(5x," Rhomin/max=",4g11.4)') c & Rhomin,Rhomin +nStepRho*dRho,Rhoc c write (*,'(5x," amin/max =",4g11.4)') c & a2min,a2min +nStepa2*da2,a2c c write (*,'(5x," xmin/max =",4g11.4)') c & x2min,x2min +nStepx2*dx2,x2c error =1.e+10 Do ir =0,nStepR r11 = Rmin +ir*dR e_old = error Do irho =0,nStepRho rho11 = Rhomin +irho*dRho Do ix2 = 0,nStepx2 x2 = x2min +ix2*dx2 Do ia2 = 0,nStepa2 a2 = a2min +ia2*da2 c write (*,'(4g11.3,4i4)') r11,rho11,x2,a2,ir,irho,ix2,ia2 ee =0. ww =0. do i=istart,imax w = (Radout(i)-Radout(i-1))/Radden(i) ee = ee + w*(log10(RhoExp( Radden(i),r11,rho11,x2,a2)/ & Dens(i)))**2 ww = ww + w EndDo c ee =ee/(imax-istart) ee =ee/ww If(ee.lt.error)then error = ee r_1 = r11 rho_1 = rho11 a_2 = a2 x_2 = x2 ir_1 = ir irho_1= irho ix_1 = ix2 ia_1 = ia2 EndIf EndDo EndDo EndDo c write (*,'(10x,"r1=",f8.3," rho1=",3g11.3," err=",g11.3,4i4)') c & r_1,rho_1,a_2,x_2,error,ir_1,irho_1,ix_1,ia_1 c If(abs(error-e_old).lt.1.e-5)Goto100 EndDo 100 on_boundary = .false. if(ir_1.eq.0.or.ir_1.eq.nStepR)on_boundary = .true. if(irho_1.eq.0.or.irho_1.eq.nStepRho)on_boundary = .true. if(ix_1.eq.0.or.ix_1.eq.nStepx2)on_boundary = .true. if(ia_1.eq.0.or.ia_1.eq.nStepa2)on_boundary = .true. iFlag = 0 If(on_boundary)iFlag =1 Rc = r_1 Rhoc= rho_1 a2c = a_2 x2c = x_2 c write (*,'(10x,"r1=",f8.3," rho1=",3g11.3," err=",g11.3,4i4)') c & r_1,rho_1,a_2,x_2,error,ir_1,irho_1,ix_1,ia_1 Return End c-------------------------------------------------------------------- c Exp-fit density profile FUNCTION RhoExp(r,r1,rho_1,x2,a2) c input: c Vvmax = estimate of Vmax c-------------------------------------------------------------------- x = r/r1 RhoExp = rho_1/x*(exp(-x)+a2*exp(-x/x2))+1. Return End c-------------------------------------------------------------------- c Exp-fit mass profile: dimensionless c F =1-(1+x)exp(-x) FUNCTION FExp(x) c-------------------------------------------------------------------- real*8 y,Yexp y= x YExp = 1.d0-(1.d0+y)*exp(-y) FExp = YExp Return End c-------------------------------------------------------------------- c Exp-fit velocity profile FUNCTION VcExp(r,r1,rho_1,x2,a2) c output: Vcirc in km/s c c-------------------------------------------------------------------- x = r/r1 VcExp = 6.698e-2*r1*sqrt(rho_1/x*(FExp(x)+a2*x2**2*FExp(x/x2))) Return End c-------------------------------------------------------------------- c Vmax and effective r_s for Exp-fit velocity profile SUBROUTINE FindVcMax(r,Vcmax,Rvir,r_s,r_1,rho_1,x2,a2,dr) c output: Vcirc in km/s c c-------------------------------------------------------------------- Vmax =0. iter =0 r = r_1 Vv0 = VcExp(r,r_1,rho_1,x2,a2) Vv1 = VcExp(r+dr,r_1,rho_1,x2,a2) If(Vv1.gt.Vv0)Then r = r-dr Else r = 0 EndIf 10 iter =iter +1 ! Find Vcmax using the approximation r = r +dr Vv= VcExp(r,r_1,rho_1,x2,a2) c write (*,*) ' VV,r=',Vv,r,Vmax If(Vv.gt.Vmax.and.r.lt.1.e4.and.r.lt.Rvir)Then Vmax = Vv GoTo 10 EndIf Vcmax = Vmax rmax = r c write (*,*) ' Vmax=',Vmax,' r=',r c r0 = r-dr/2. c V0= VcExp(r0,r_1,rho_1,x2,a2) c r1 = r-dr c V1= VcExp(r1,r_1,rho_1,x2,a2) c write (*,*) ' V0=',V0,' r0=',r0 c write (*,*) ' V1=',V1,' r1=',r1 Vmax =0.945*Vcmax V2 = Vmax iter =0 r =0. c write (*,*) ' V(rs)=',V2 20 iter =iter +1 ! Find r_s: Vc(rs)=0.945Vmax r = r +dr Vv= VcExp(r,r_1,rho_1,x2,a2) If(Vv.lt.Vmax.and.r.lt.1.e4)Then GoTo 20 EndIf r_s =r r =rmax Return End c-------------------------------------------------------------------- c Using density profile find Vmax and effective r_s for Exp-fit velocity profile SUBROUTINE FindRs(r_s,r_1,rho_1,x2,a2,dr,rhs,Rout) c output: Vcirc in km/s, r_s (kpc) c solve eq: rho(rs) = rho(Rout)/4*Rout/rs*(1+Rout/rs)**2 c-------------------------------------------------------------------- Vmax =0. iter =0 c write (*,*) ' Rout=',Rout,' rhs=',rhs r = r_1/2. 10 rho = RhoExp(r,r_1,rho_1,x2,a2) rho_a= rhs*Rout/r*(1.+Rout/r)**2 c write (*,*) ' rho=',rho,' rho_a=',rho_a,r if(rho_a.lt.rho)Then r =r/2. If(r.lt.1.e-2)Then write (*,*) 'FindRs: initial guess failed' r_s = r_1/2 return EndIf goto 10 EndIf 20 iter =iter +1 ! Find Vcmax using the approximation r = r +dr rho = RhoExp(r,r_1,rho_1,x2,a2) rho_a= rhs*Rout/r*(1.+Rout/r)**2 c write (*,*) ' rho=',rho,' rho_a=',rho_a,r if(rho_a.gt.rho.and.r.lt.1.e3)goto20 r_s =r Return End