c----------------------------------------------------------------
c                                Find <r>, <M>, period  for
c                                trajectories in  NFW potential
c
      Include 'NFWr.h'
      EXTERNAL F, EUeff, AvRadius,AvMass,Period,AvMR

      CALL InitValues
      Open(1,file='Result.dat')
c
      write(*,'("x",T12,"<x>",T24,"ratio",
     &           T36,"Mass(x)",T48,"Mass(<x>)",
     &           T60,"ratio",T72,"V_c",T84,"RMS V   <r_min/r_max>")')
c                                                     Loop over radii
      do i=0,150
         x=0.001*10.**(i/10.)
         If(x.ge.1.d-3.and.x.lt.min(x_out,10.))Then
            Sigma =sqrt(Sigma_r(x))
            rRad =0.
            rMR  =0.
            eccen =0.
            aLs   = 0.
         Do j=1,norb
 1          CALL SetOrbit2(x)
            CALL Rminmax(r_min,r_max,EUeff)
            If(r_max.ge.x_out)GoTo1       ! take only orbit inside x_out
            If(r_max.gt.r_min*1.01)Then   ! non circular orbit
               r_max = r_max/1.001           ! slightly change limits 
               r_min = r_min*1.001             ! to avoid singularities

            CALL DGAUS8 (Period, r_min,r_max, eps2, Tperiod, IERR)
                  if(IERR.ne.1)write (*,*)  ' ERROR Period=',IERR
            CALL DGAUS8 (AvRadius, r_min,r_max, eps2, RadAver, IERR)
                  if(IERR.ne.1)write (*,*)  ' ERROR AvRadius=',IERR
c            CALL DGAUS8 (AvMR, r_min,r_max, eps2, aMRAver, IERR)
c                  if(IERR.ne.1)write (*,*)  ' ERROR AvMR=',IERR
c               rMR       = rMR  + aMRAver/Tperiod
                  RadAver = RadAver/Tperiod
            Else                ! circular orbit
                RadAver = r_min
                r_max      = r_min
            EndIf
               aLs        = aLs + aL2
               rRad      = rRad  + RadAver
               eccen = eccen + r_min/r_max
         EndDo
         rRad      = rRad/norb
         eccen    = eccen/norb
         aLs         = sqrt(aLs/norb)
               write(*,100) x,rRad,rRad/x,F(x),F(rRad),
     &               F(rRad)/F(x),sqrt(F(x)/x),sqrt(3.)*Sigma,eccen,
     &               aLs/sqrt(F(x)*x)
 100     format(g11.4,T12,g11.4,T24,f7.3,T36,g11.4,T48,g11.4,T60,g11.4,
     &                        T72,g11.4,T84,3g11.4)
         EndIf
      EndDo
      Stop
      End

C---------------------------------------
       SUBROUTINE SetOrbit
c                    set energy and angular momentum
c                    for orbit with given xmax and alpha
c                      L =alpha*L_circular(xmax)
C---------------------------------------
      Include 'NFWr.h'
         En  = Phi(xmax) + 0.5*alpha**2 *F(xmax)/xmax
         aL2 = alpha**2 *F(xmax)*xmax
      RETURN
      END

C---------------------------------------
       SUBROUTINE SetOrbit2(x)
c                    set energy and angular momentum
c                    for orbit with random rad and perp velocities
c                      L =alpha*L_circular(xmax)
C---------------------------------------
      Include 'NFWr.h'
 1      Vrad  = Gauss(Nseed)*Sigma
         Vperp= Gauss(Nseed)*Sigma*sq2
         aL2     = (Vperp*x)**2
         En       = Phi(x) + 0.5d0*(Vperp**2+Vrad**2)
         If(En.ge.Phimax)goto1
      RETURN
      END

C---------------------------------------
       SUBROUTINE InitValues
c                    Initialize parameters
c                    for NFW and truncated Isothermal
C---------------------------------------
      Include 'NFWr.h'

      x_out   =35.             ! truncation for NFW
      x1out = 1./(1.+x_out)
      
      norb = 50000
      Nseed = 121071

      eps1 = 1.d-7           ! errors of intigration
      eps2 = 1.d-5
      RETURN
      END

C---------------------------------------
       SUBROUTINE Rminmax(r_min,r_max,Func)
c                     Find r_min and r_max for NFW
C---------------------------------------
      Include 'NFWr.h'
      External Func
      r0 =3.e-5                  ! find radius r0 such that Func(r0) <0
      n = 0
 1    r0 =r0/1.5
      n = n+1
c      write (*,*)  '  Rminmax: r0=',r0,n,Func(r0)
         if(n.gt.100)Then
            write (*,*) '             r0=',r0,Func(r0)       
         Do i=1,30
            r= 0.001*10.**(i/15.)
            write (*,300) r,Func(r)
 300        format(' r=',g11.3,' F=',g11.3)
         EndDo
         Stop
         EndIf
      If(Func(r0).gt.0.)goto1

 6    fact =1.1
      r1 =r0
      n = 0
      ff0= -1.e+15 
 2    r1 =r1*fact
      ff = ff0
      n = n+1
      ff0 = Func(r1)
c      write (*,*)  '   r1=',r1,' F=',ff0
      If(ff0.lt.0.and.ff0.gt.ff)goto2
      If(ff0.lt.0.)Then
c         write(*,*) ' r1=',r1,' F=',ff0,' previous F=',ff
         fact =1.+(fact-1.)/15.
c         EnIso =EnIso/1.05

c         write(*,*) ' r1=',r1,' Fnew=',Func(r1)
         fact = 1./fact
c         write (*,*) '  change factor=',fact
         n =0
 5       r1 =r1*fact
         n = n +1
         ff0 = Func(r1)
         if(n.gt.50)Then ! this is circular orbit
                                  ! make it a bit elliptical
c           write (*,*) ' too many iterations for r1=',r1,ff0
           r1 =r0
          En =En-ff0+1.e-10
          fact =1./fact
          ff0 = -1.e+15
          goto2
         endIf
         If(ff0.lt.0.)Goto5
      EndIf
      If(Func(r0)*Func(r1).gt.0.)Then
         write(*,'(" error1: r=",2g11.3," Func=",2g11.3)')
     &   r0,r1, Func(r0),Func(r1)
         STOP
       EndIf
      n = 0
      r2 =r1
 3    r2 =r2*1.1
      n = n+1
      If(n.gt.250)Then
         write (*,*) '             r2=',r2, Func(r2)   
         write (*,*) '             r0=',r0, Func(r0),En
         rr =r0
         Do ii =0,1500
            rr =rr*1.02
            write (*,*) '             rr=',rr, Func(rr)
         EndDo 
         STOP
      EndIf 
      If(Func(r2).gt.0.)goto3
c      write (*,*) ' first guess =',r0,r1,r2,'   Functions=',Func(r0),
c     &        Func(r1), Func(r2)  
      rr = r1
      ru =r0
      r_guess =(r0+r1)/2.
      CALL DFZERO(Func,r0,rr,r_guess,1.d-8,1.d-8,iflag)
      If(iflag.ne.1)write(*,*)' error1 zero not reached',r0,ru,r1
      r_min =r0
      r_guess =(r1+r2)/2.
      CALL DFZERO(Func,r1,r2,r_guess,1.d-8,1.d-8,iflag)
      If(iflag.ne.1)write(*,*)' error2 zero not reached',r0,r1
      r_max =r1
c      write (*,50) r_min,r_max,Func(r_min) ,Func(r_max) 
c 50   format(' rmin/max=',2g11.4,' functions=',2g11.4)
c      write (*,*)  ' Func(xmax)=',Func(xmax),' xmax=',xmax
      RETURN
      END

C---------------------------------------
       DOUBLE PRECISION FUNCTION F(x)
c                              Mass for NFW
c                               x must be smaller than xout
C---------------------------------------
      Include 'NFWr.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---------------------------------------
       DOUBLE PRECISION FUNCTION Period(x)
c                              integrant for radial period
c                              
C---------------------------------------
      Include 'NFWr.h'
         Period = 1./sqrt(ABS(2.*EUeff(x)))
      RETURN
      END

C---------------------------------------
       DOUBLE PRECISION FUNCTION AvRadius(x)
c                              integrant for average radius
c                              
C---------------------------------------
      Include 'NFWr.h'
         AvRadius = x/sqrt(ABS(2.*EUeff(x)))
      RETURN
      END
C---------------------------------------
       DOUBLE PRECISION FUNCTION AvMass(x)
c                              integrant for average mass
c                              
C---------------------------------------
      Include 'NFWr.h'
         AvMass = F(x)/sqrt(ABS(2.*EUeff(x)))
      RETURN
      END

C---------------------------------------
       DOUBLE PRECISION FUNCTION AvMR(x)
c                              integrant for average massxRadius
c                              
C---------------------------------------
      Include 'NFWr.h'
         AvMR = F(x)*x/sqrt(ABS(2.*EUeff(x)))
      RETURN
      END

C---------------------------------------
       DOUBLE PRECISION FUNCTION Phi(x)
c                              Potential for truncated NFW
c                              Phi(x) < 0
C---------------------------------------
      Include 'NFWr.h'
      If(x.le.x_out)Then
         Phi = -Log(1.+x)/x +x1out
      Else
         Phi = -F(x_out)/x
      EndIf
      RETURN
      END

C---------------------------------------
       DOUBLE PRECISION FUNCTION EUeff(x)
c                    E-Ueff(x)
c                          
C---------------------------------------
      Include 'NFWr.h'
         EUeff = En -Phi(x)-aL2/2./x**2
      RETURN
      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 'NFWr.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 =(x*(1.+x)**2) *Sigma_r
      EndIf

      RETURN
      END

C---------------------------------------
       DOUBLE PRECISION FUNCTION Fsigma(x)
c                     Integrant for Sigma_r f(x)/x^3/(1+x)^2
C---------------------------------------
      Include 'NFWr.h'
         z      = 1.+x
         Fsigma =(LOG(z)-x/z)/x**3/z**2 
      RETURN
      END
C------------------------------------------------
C				                                       random number generator
      DOUBLE PRECISION 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  RANDd
      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