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
