c-------------------------------------------------------------------- c Velocity dispersion of NFW c c-------------------------------------------------------------------- Include 'NFWdisp.h' character*80 Name REAL*4 INPUT external Fproj x_out =20. ! truncation for NFW x1out = 1./(1.+x_out) Conc = 12. aMvir = 1.d12 deltacr= 340. Om0 = 0.3 aMvir = INPUT('Enter Virial Mass =') Conc = INPUT('Enter Concentration=') fcons = F(Conc)/Conc**3/deltacr rvir = 0.09534*(aMvir/Om0/deltacr)**(1./3.) rs = rvir/Conc write(Name,'(a,f4.1,a,i2.2,a)')'NFWdispMvir',(log10(amvir)), & 'C',INT(Conc),'.dat' open(1,file=Name) do j=0,1 write (j,'(a,g12.4,a,f6.1,2(a,f8.2,2x))') & '# Mvir=',aMvir,' C=',Conc, & ' Rvir(kpch)=',rvir,' Rs(kpch)=',rs write (j,'(a,T12,a,T24,a,T36,a,a)')'# R(kpch)',' Mass', & ' Vcirc',' Vrms(1d kms)',' overdenisty' EndDo vcons = 2.081d-3*sqrt(aMvir/rs/F(Conc)) eps1 = 1.d-7 ! errors of intigration eps2 = 1.d-5 c Loop over radii do i=0,150 r=1.*10.**(i/15.) x= r/rs If(r.ge.1.d0.and.r.lt.2000.)Then Sigma =vcons*sqrt(Sigma_r(x)) aMass = F(x)+fcons*x**3 xp = x*0.99990 CALL DGAUS8 (Fproj, x,x_out, eps2, ProjSurf, IERR) do j=0,1 write(j,100) x*rs,aMvir*aMass/F(Conc), & vcons*sqrt(aMass/x),Sigma, & rho(x)/(3.*fcons),ProjSurf EndDo 100 format(g11.4,T12,g11.4,T24,f7.3,T36,g11.4,T48,g11.4,T60,g11.4, & T72,g11.4,T84,4g11.4) EndIf EndDo Stop end C---------------------------------- Read in variables REAL FUNCTION INPUT(text) C------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End C--------------------------------------- DOUBLE PRECISION FUNCTION Sigma_r(x) c radial velocity dispersion for truncated NFW C--------------------------------------- Include 'NFWdisp.h' External Fsigma if(x.ge.x_out)Then Sigma_r =0. Else CALL DGAUS8 (Fsigma, x,x_out, eps2, Sigma_r, IERR) Sigma_r =Sigma_r/rho(x) EndIf RETURN END C--------------------------------------- DOUBLE PRECISION FUNCTION Fproj(x) c Integrant for ProjSurf rho(x)x/sqrt(x**2-Rp**2) C--------------------------------------- Include 'NFWdisp.h' Fproj =x*rho(x)/sqrt(x**2-xp**2) RETURN END C--------------------------------------- DOUBLE PRECISION FUNCTION Fsigma(x) c Integrant for Sigma_r f(x)/x^3/(1+x)^2 C--------------------------------------- Include 'NFWdisp.h' Fsigma =(F(x)+fcons*x**3)*rho(x)/x**2 RETURN END C--------------------------------------- DOUBLE PRECISION FUNCTION rho(x) c Integrant for Sigma_r f(x)/x^3/(1+x)^2 C--------------------------------------- Include 'NFWdisp.h' rho =1.d0/x/(1.d0+x)**2+3.*fcons RETURN END C--------------------------------------- DOUBLE PRECISION FUNCTION F(x) c Mass for NFW c x must be smaller than xout C--------------------------------------- Include 'NFWdisp.h' If(x.le.x_out)Then F = Log(1.+x) -x/(1.+x) Else F = Log(1.+x_out) -x_out/(1.+x_out) EndIf RETURN END c-------------------------------------------------------------------- SUBROUTINE DGAUS8 (FUN, A, B, ERR, ANS, IERR) C***BEGIN PROLOGUE DGAUS8 C***PURPOSE Integrate a real function of one variable over a finite C interval using an adaptive 8-point Legendre-Gauss C algorithm. Intended primarily for high accuracy C integration or integration of smooth functions. C***LIBRARY SLATEC C***CATEGORY H2A1A1 C***TYPE DOUBLE PRECISION (GAUS8-S, DGAUS8-D) C***KEYWORDS ADAPTIVE QUADRATURE, AUTOMATIC INTEGRATOR, C GAUSS QUADRATURE, NUMERICAL INTEGRATION C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C Abstract *** a DOUBLE PRECISION routine *** C DGAUS8 integrates real functions of one variable over finite C intervals using an adaptive 8-point Legendre-Gauss algorithm. C DGAUS8 is intended primarily for high accuracy integration C or integration of smooth functions. C C The maximum number of significant digits obtainable in ANS C is the smaller of 18 and the number of digits carried in C double precision arithmetic. C C Description of Arguments C C Input--* FUN, A, B, ERR are DOUBLE PRECISION * C FUN - name of external function to be integrated. This name C must be in an EXTERNAL statement in the calling program. C FUN must be a DOUBLE PRECISION function of one DOUBLE C PRECISION argument. The value of the argument to FUN C is the variable of integration which ranges from A to B. C A - lower limit of integration C B - upper limit of integration (may be less than A) C ERR - is a requested pseudorelative error tolerance. Normally C pick a value of ABS(ERR) so that DTOL .LT. ABS(ERR) .LE. C 1.0D-3 where DTOL is the larger of 1.0D-18 and the C double precision unit roundoff D1MACH(4). ANS will C normally have no more error than ABS(ERR) times the C integral of the absolute value of FUN(X). Usually, C smaller values of ERR yield more accuracy and require C more function evaluations. C C A negative value for ERR causes an estimate of the C absolute error in ANS to be returned in ERR. Note that C ERR must be a variable (not a constant) in this case. C Note also that the user must reset the value of ERR C before making any more calls that use the variable ERR. C C Output--* ERR,ANS are double precision * C ERR - will be an estimate of the absolute error in ANS if the C input value of ERR was negative. (ERR is unchanged if C the input value of ERR was non-negative.) The estimated C error is solely for information to the user and should C not be used as a correction to the computed integral. C ANS - computed value of integral C IERR- a status code C --Normal codes C 1 ANS most likely meets requested error tolerance, C or A=B. C -1 A and B are too nearly equal to allow normal C integration. ANS is set to zero. C --Abnormal code C 2 ANS probably does not meet requested error tolerance. C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, I1MACH, XERMSG C***REVISION HISTORY (YYMMDD) C 810223 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890911 Removed unnecessary intrinsics. (WRB) C 890911 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900326 Removed duplicate information from DESCRIPTION section. C (WRB) C***END PROLOGUE DGAUS8 INTEGER IERR, K, KML, KMX, L, LMN, LMX, LR, MXL, NBITS, 1 NIB, NLMN, NLMX INTEGER I1MACH DOUBLE PRECISION A,AA,AE,ANIB,ANS,AREA,B,C,CE,EE,EF, 1 EPS, ERR, EST, GL, GLR, GR, HH, SQ2, TOL, VL, VR, W1, W2, W3, 2 W4, X1, X2, X3, X4, X, H, 3 D1MACH5, D1MACH4 DOUBLE PRECISION D1MACH, G8, FUN DIMENSION AA(60), HH(60), LR(60), VL(60), GR(60) SAVE X1, X2, X3, X4, W1, W2, W3, W4, SQ2, 1 NLMN, KMX, KML DATA X1, X2, X3, X4/ 1 1.83434642495649805D-01, 5.25532409916328986D-01, 2 7.96666477413626740D-01, 9.60289856497536232D-01/ DATA W1, W2, W3, W4/ 1 3.62683783378361983D-01, 3.13706645877887287D-01, 2 2.22381034453374471D-01, 1.01228536290376259D-01/ DATA SQ2/1.41421356D0/ DATA NLMN/1/,KMX/5000/,KML/6/ G8(X,H)=H*((W1*(FUN(X-X1*H) + FUN(X+X1*H)) 1 +W2*(FUN(X-X2*H) + FUN(X+X2*H))) 2 +(W3*(FUN(X-X3*H) + FUN(X+X3*H)) 3 +W4*(FUN(X-X4*H) + FUN(X+X4*H)))) C***FIRST EXECUTABLE STATEMENT DGAUS8 C C Initialize C D1MACH4 = 1.d-17 c K = I1MACH(14) c ANIB = D1MACH(5)*K/0.30102000D0 c NBITS = ANIB c NLMX = MIN(60,(NBITS*5)/8) NBITS = 16 NLMX = MIN(60,(NBITS*5)/8) ANS = 0.0D0 IERR = 1 CE = 0.0D0 IF (A .EQ. B) GO TO 140 LMX = NLMX LMN = NLMN IF (B .EQ. 0.0D0) GO TO 10 IF (SIGN(1.0D0,B)*A .LE. 0.0D0) GO TO 10 C = ABS(1.0D0-A/B) IF (C .GT. 0.1D0) GO TO 10 IF (C .LE. 0.0D0) GO TO 140 ANIB = 0.5D0 - LOG(C)/0.69314718D0 NIB = ANIB LMX = MIN(NLMX,NBITS-NIB-7) IF (LMX .LT. 1) GO TO 130 LMN = MIN(LMN,LMX) 10 TOL = MAX(ABS(ERR),2.0D0**(5-NBITS))/2.0D0 IF (ERR .EQ. 0.0D0) TOL = SQRT(D1MACH4) EPS = TOL HH(1) = (B-A)/4.0D0 AA(1) = A LR(1) = 1 L = 1 EST = G8(AA(L)+2.0D0*HH(L),2.0D0*HH(L)) K = 8 AREA = ABS(EST) EF = 0.5D0 MXL = 0 C C Compute refined estimates, estimate the error, etc. C 20 GL = G8(AA(L)+HH(L),HH(L)) GR(L) = G8(AA(L)+3.0D0*HH(L),HH(L)) K = K + 16 AREA = AREA + (ABS(GL)+ABS(GR(L))-ABS(EST)) C IF (L .LT .LMN) GO TO 11 GLR = GL + GR(L) EE = ABS(EST-GLR)*EF AE = MAX(EPS*AREA,TOL*ABS(GLR)) IF (EE-AE) 40, 40, 50 30 MXL = 1 40 CE = CE + (EST-GLR) IF (LR(L)) 60, 60, 80 C C Consider the left half of this level C 50 IF (K .GT. KMX) LMX = KML IF (L .GE. LMX) GO TO 30 L = L + 1 EPS = EPS*0.5D0 EF = EF/SQ2 HH(L) = HH(L-1)*0.5D0 LR(L) = -1 AA(L) = AA(L-1) EST = GL GO TO 20 C C Proceed to right half at this level C 60 VL(L) = GLR 70 EST = GR(L-1) LR(L) = 1 AA(L) = AA(L) + 4.0D0*HH(L) GO TO 20 C C Return one level C 80 VR = GLR 90 IF (L .LE. 1) GO TO 120 L = L - 1 EPS = EPS*2.0D0 EF = EF*SQ2 IF (LR(L)) 100, 100, 110 100 VL(L) = VL(L+1) + VR GO TO 70 110 VR = VL(L+1) + VR GO TO 90 C C Exit C 120 ANS = VR IF ((MXL.EQ.0) .OR. (ABS(CE).LE.2.0D0*TOL*AREA)) GO TO 140 IERR = 2 c write (*,*) 'DGAUS8: ', c + 'ANS is probably insufficiently accurate.' GO TO 140 130 IERR = -1 write (*,*) 'DGAUS8: ', + 'A and B are too nearly equal to allow normal integration.' 140 IF (ERR .LT. 0.0D0) ERR = CE RETURN END C------------------------------------------------ C random number generator FUNCTION RANDd(M) C------------------------------------------------ DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4 + /453815927,2147483648.,2147483647,536870912,131072,256, + 16777216,4,16384,8388608,128/ ML=M/K1*K1 M1=(M-ML)*L1 ML=M/K2*K2 M2=(M-ML)*L2 ML=M/K3*K3 M3=(M-ML)*L3 ML=M/K4*K4 M4=(M-ML)*L4 M5=KI-M IF(M1.GE.M5)M1=M1-KI-1 ML=M+M1 M5=KI-ML IF(M2.GE.M5)M2=M2-KI-1 ML=ML+M2 M5=KI-ML IF(M3.GE.M5)M3=M3-KI-1 ML=ML+M3 M5=KI-ML IF(M4.GE.M5)M4=M4-KI-1 ML=ML+M4 M5=KI-ML IF(LC.GE.M5)ML=ML-KI-1 M=ML+LC RANDd=M/AM RETURN END C-------------------------------------- C normal random numbers DOUBLE PRECISION FUNCTION GAUSS(M) C-------------------------------------- DOUBLE PRECISION X,X2 X=0. DO I=1,5 X=X+RANDd(M) EndDo X2 =1.5491933*(X-2.5) GAUSS=X2*(1.-0.01*(3.-X2*X2)) RETURN END