
*               --------------------------------------------
*               MONTE CARLO SIMULATION OF COMPTON SCATTERING
*               --------------------------------------------
*
*  @ GEOMETRY ... sphere with taus
*                 taues = sigma_T r n_e
*                          sigma_T ... Thomson scattering depth
*                          r ... radius
*                          n_e ... electron density
*
*  @ ELECTRONS ... isotropic, non-thermal + thermal
*  @ INPUT PHOTONS ... planckian
*
* + INPUT :
*     
*     amxwl: fraction of thermal electrons
*            (1-amxwl) ... fraction of non-thermal electrons
*     te:    temperature of thermal electrons
*     p:     power-law index of nonthermal electrons
*     gam2:  high-energy cutoff of electron's Lorentz factor
*
*     trad:   temperature of soft photons (planck distribution)
*             ( unit = electron mass )
*     (tradk ... soft photon temperature in K)
*
*     taues:  electron scattering depth
*
*     nphoton: photon number
*
*     epsct:  cut off of weight of photon distribution
*
* + OUTPUT :
*     sinpt(j) : relative value of input energy flux from per x
*     snj(j)   : relative value of emitted energy flux from a sphere
*     dsnj(j)  : error ... one sigma
*     mean photon enregies :
*         x0mean ... input  photon energy
*         x1mean ... output photon energy
*
*
*  SYMBOLS :
*     xhn(j) ... photon energy in units of electron mass
*                ( dimension = nx0 )
*     w(i) ... weight
*     --  --  --  --  --  --  --  --  --  -- --  --  --  --  --
*                             V. 0.1   BY M. K.          1/5/96
*                                1.0                    8/26/96
*                                1.01                   9/11/96
*                                1.02                  12/24/96
*                                1.1                   12/28/96
*     --  --  --  --  --  --  --  --  --  -- --  --  --  --  --
*  REFERENCE :
*     SOVIET SCIENTIFIC REVIEWS SECTION E
*     ASTROPHYSICS AND SPACE PHYSICS REVIEWS, VOL. 2 (1983)
*     PAGES 189 - 331
*      BY L.A. POZDNYAKOV, I.M. SOBOL', AND R. SYUNYAEV
* -----------------------------------------------------------
*
      CHARACTER*12 index, hi_gam, tradch, tauch, nch, epsch,
     &             frac_maxwl, el_temp
      INTEGER NX0, NGDIM
      PARAMETER( NX0 = 1026 )
      PARAMETER( NGDIM = 128 )
*
      DOUBLE PRECISION PI, EMASS, VELC, BOLTZK, SBC
      PARAMETER( PI = 3.1415927D0 )
      PARAMETER( EMASS = 9.109534D-28, VELC = 2.99792458D10,
     &           BOLTZK = 1.3807D-16, SBC = 5.67051D-05 )
*
*                                    ! EMASS ... electron mass
*                                    ! VELC ... light velocity
*                                    ! BOLTZK ... Boltzman constant
*                                    ! SBC ... Stefan-Boltzmann constant
*
      INTEGER nx, noph
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION sinpt
      COMMON /penbin/ sinpt(NX0), noph(NX0)
      DOUBLE PRECISION snj, snjs, dsnj, esbj
      COMMON /statis/ snj(NX0), snjs(NX0), dsnj(NX0), esbj(NX0)
      INTEGER nspace
      DOUBLE PRECISION g_sp, p_sp, gdis, pdis
      COMMON /used_e/ g_sp(1000), p_sp(1000), gdis(1000), pdis(1000),
     &                nspace
      DOUBLE PRECISION th_nrm
      COMMON /maxnrm/ th_nrm
      DOUBLE PRECISION gmin
      COMMON /upmxwl/ gmin

      INTEGER n2
      DOUBLE PRECISION gsp(1000), fdens(1000), y2(1000)
*
      INTEGER ijkl, i, j, ntotin, nphoton, imin
      DOUBLE PRECISION p, gam1, gam2, tradkv, taues,
     &     epsct, tradk, trad, xlow, xhigh, dxl, c_nom, x0,
     &     snjwa, x0mean, x1mean, acomp
      DOUBLE PRECISION amxwl, te, gdnorm, pdnorm, thnrm0
c      DOUBLE PRECISION rk2e
c      DOUBLE PRECISION summax, sumnon



      CHARACTER*12, spfile, elfile
      CHARACTER*24, spfname, elfname

*
*
*
*  ! Input parameters
*
      open(8,file='ip1228',status='old')
        read(8,*) spfile,     spfname
        read(8,*) elfile,     elfname
        read(8,*) frac_maxwl, amxwl
        read(8,*) el_temp,    te
        read(8,*) index,      p
        read(8,*) hi_gam,     gam2
        read(8,*) tradch,     tradkv
        read(8,*) tauch,      taues
        read(8,*) nch,        nphoton
        read(8,*) epsch,      epsct
      close(8)
*
      tradk = tradkv / 510.99906e0 * 5.92986e9
      trad = BOLTZK * tradk / EMASS / VELC / VELC
*
* ------------------------------------------------------------
*        Get gam1
*        gam1:  low-energy cutoff of electron's Lorentz factor
*
      if( abs(amxwl-1.d0) .gt. 1.d-7 ) then
         call f_gmin( te, amxwl, p, gam2, gam1 )
      else
         gam1 = 100.d0*te
      end if
      gmin = gam1
*
*        If you want to set amxwl = 0.0, comment above line
*        and set
cc        gam1 = 1.1d0
*        (the value of the right hand side is arbitrary
*         as long as gam1 > 1.0)      
*
* ------------------------------------------------------------
*
      write(6,905)
      write(6,910) amxwl, 1.d0-amxwl, te, te*511.d0,
     &             p, gam1, gam2,
     &             trad, tradkv, tradk, 
     &             taues, nphoton, epsct
*
* ------------------------------------------------------
*                                    ! mesh of photon energy and angle
      xlow = 1.d-10
      xhigh = 1.d6
      dxl = 0.2d0
      call mesh( xlow, xhigh, dxl )
*
*                             ! normalization of electron distribution
      call el_nml( p, gam1, gam2, c_nom )
*
*                             ! mesh in electron Lorentz factor
      call g_mesh( p, gam1, gam2, c_nom )
*
*                             ! normalization for Maxwellian
*                               mean free path used in meanen
c      th_nrm = 1.d0/(4.d0*PI*rk2e(1.d0/te))
      call thnrml( te, gam1, thnrm0 )
      th_nrm = te/exp(1.d0/te) * thnrm0
*
*                             ! make an interpolation table
*                               to select electron momentum
      call table( amxwl, thnrm0, c_nom, te, gam1, gam2, p,
     &            gsp, fdens, y2, n2 )
*
*                                  ! initialize random number generator
* the input value should be in the range: 0 < = ijkl < = 900 000 000
*
      ijkl = 199615
      call rmarin(ijkl)
*
*  ------------------------------------------------------
      call sumgb( te, gam1 )
      call sumofm
      call initfx
*
*  -------------- BEGINNING OF SCATTERING ---------------
*
      do 51 i = 1, nphoton
         if( mod(i,1000) .eq. 0 ) then
            open(25,file='number',status='unknown')
            write(25,*) 'number =',i
            close(25)
         end if
         call phnxyz
         call phengy( trad, x0 )
         call scattr( amxwl, te, p, gam1, gam2, taues, x0, epsct,
     &                gsp, fdens, y2, n2 )
         call sums
 51   continue
*
*  ---------------- END OF SCATTERING --------------------
*
*                                                     ! normalization
      call normal( nphoton, ntotin, snjwa )
*
*                                              ! mean photon energies
      call meanen( trad, x0mean, x1mean, acomp )
*
*  ---------------------------------------------
*
      open(12,file=spfname,status='unknown')
      write(12,941) ( xhn(j), sinpt(j), snj(j), dsnj(j),
     &                j = 1, nx )
 941  format(1h ,1pe11.4,2x,e11.4,2x,e11.4,2x,e11.4)
      close(12)
*
* =-------------------------------------------------------=
      do 121 i = 1, nspace
         if( g_sp(i) .le. gmin ) imin = i
 121  continue
      write(6,*) 'imin =',imin
      summax = 0.d0
      do 122 i = 1, imin
         summax = summax + gdis(i)
 122  continue
      summax = summax*log( g_sp(2)/g_sp(1) )
      sumnon = 0.d0
      do 123 i = imin+1, nspace
         sumnon = sumnon + gdis(i)
 123  continue
      sumnon = sumnon*log( g_sp(2)/g_sp(1) )
      write(6,*) 'summax =',summax,'  sumnon =',sumnon
*
      gdnorm = 0.d0
      pdnorm = 0.d0
      do 151 i = 1, nspace
         gdnorm = gdnorm + gdis(i)
         pdnorm = pdnorm + pdis(i)
 151  continue
      gdnorm = gdnorm * log( g_sp(2)/g_sp(1) )
      pdnorm = pdnorm * log( p_sp(2)/p_sp(1) )
*
      do 155 i = 1, nspace
         gdis(i) = gdis(i) / g_sp(i) / gdnorm 
         pdis(i) = pdis(i) / p_sp(i) / pdnorm 
 155  continue
*
      open(14,file=elfname,status='unknown')
      write(14,950) (g_sp(i), gdis(i), p_sp(i), pdis(i), i = 1, nspace)
 950  format(1h ,1pe12.5,2x,e12.5,2x,e12.5,2x,e12.5)
      close(14)
*  ---------------------------------------------
*
 905  format(1h ,5x,'COMPTON SCATTERING BY RELATIVISTIC ELECTRONS'/
     &    1x,'   + Geometry ... Sphere with tau_es'/
     &    1x,'   + Electrons ... Thermal + Nonthermal'/
     &    1x,'   + Input photons ... Planckian with T_rad'/ )
 910  format(1h ,
     &  1x,'Thermal     ... density fraction =',1pe10.3/
     &  2x,'Non-thermal ... density fraction =',1pe10.3/
     &  2x,'Te/mc^2 =',e11.4,2x,'Te(keV) =',e11.4/
     &  2x,'p =',1pe11.4,2x,'gamma_low =',e11.4,
     &  2x,'gamma_high =',e11.4/
     &  2x,'kT_rad/mc^2 =',e11.4,2x,'T_rad(keV) =',e11.4,2x,
     &  'T_rad(K) =',e11.4/
     &  2x,'tau_es =',e11.4/
     &  2x,'Photon Number =',i7,2x,'eps(cutoff) =',e10.3/)
 925  format(1h ,'Radius of a sphere R (cm) =',1pe11.4/
     &    1x,'compton cooling rate per unit volume =',e11.4/)
*
      STOP
      END
*
*
      SUBROUTINE mesh( xlow, xhigh, dxl )
      DOUBLE PRECISION xlow, xhigh, dxl
* ----------------------------------------
*          mesh in photon energy
* ----------------------------------------
      INTEGER NX0
      PARAMETER( NX0 = 1026 )
      INTEGER nx
      DOUBLE PRECISION xhn
      COMMON /ENRBIN/ xhn(0:NX0), nx
      INTEGER j
      DOUBLE PRECISION xde
*
      nx = int( log( xhigh / xlow ) / dxl ) + 2
      if( nx .ge. nx0 ) then
         write(6,900) nx, NX0
 900     format(1h0,'error in energy mesh points !!!'/
     &     'condition nx < NX0 is violated, nx =',
     &     i5,2x,'NX0 =',i5/)
         pause 'stop in mesh'
      end if
*
      xhn(1) = xlow
      xhn(0) = 0.d0
      xhn(nx+1) = 1.d20
      xde = exp( dxl )
      do 10 j = 2, nx
         xhn(j) = xhn(j-1) * xde
 10   continue
*
      write(6,910) xhn(1), xhn(nx), nx
  910 format(1h ,'x(1) =',1pe10.3, 2x,'x(nx) =',e10.3,2x,'nx =',i4/)
*
      return
      END
*
*
      SUBROUTINE el_nml( p, gam1, gam2, c_nom )
      DOUBLE PRECISION p, gam1, gam2, c_nom
* --------------------------------------------------
*        Normalization of electron distribution
* --------------------------------------------------
      DOUBLE PRECISION PI
      PARAMETER( PI = 3.1415927D0 )
      DOUBLE PRECISION al1
*
      al1 = 1.d0 - p
      if( abs(al1) .gt. 1.d-8) then
         c_nom = al1/(gam2**al1-gam1**al1)/(4.d0*PI)
      else
         c_nom = 1.d0/log(gam2/gam1)/(4.d0*PI)
      end if
      return
      END
*
*
      SUBROUTINE g_mesh( p, gam1, gam2, c_nom )
      DOUBLE PRECISION p, gam1, gam2, c_nom
      INTEGER NGDIM
      PARAMETER( NGDIM = 128 )
      INTEGER ng_nth
      DOUBLE PRECISION gi, gbm, gbp, disfc, delg
      COMMON /gint/gi(NGDIM), gbm(NGDIM), gbp(NGDIM), disfc(NGDIM),
     &             delg, ng_nth
      INTEGER nspace
      DOUBLE PRECISION g_sp, p_sp, gdis, pdis
      COMMON /used_e/ g_sp(1000), p_sp(1000), gdis(1000), pdis(1000),
     &                nspace
*
      INTEGER i
      DOUBLE PRECISION exdel, dgx, dpx
*
      ng_nth = NGDIM
      gi(1) = gam1
      gi(ng_nth) = gam2
      delg = log( gi(ng_nth)/gi(1) ) / float( ng_nth-1 )
      exdel = exp( delg )
      do 11 i = 2, ng_nth-1
         gi(i) = gi(i-1) * exdel
 11   continue
      do 21 i = 1, ng_nth
         gbp(i) = gi(i) + sqrt( (gi(i)-1.d0)*(gi(i)+1.d0) )
         gbm(i) = 1.d0/gbp(i)
 21   continue

      do 31 i = 1, ng_nth
         disfc(i) = c_nom / gi(i)**p 
     &             / sqrt( (gi(i)-1.d0)*(gi(i)+1.d0) )
 31   continue

*
*                    gamma and momentum meshes for output
      nspace = 128
      g_sp(1) = 1.d0
      g_sp(nspace) = 1.1d0*gi(ng_nth)
      dgx = log( g_sp(nspace)/g_sp(1) ) / float( nspace-1 )
      dgx = exp(dgx)
      p_sp(1) = 1.d-8
      p_sp(nspace) = 1.1d0*sqrt( g_sp(nspace)*g_sp(nspace)-1.d0 )
      dpx = log( p_sp(nspace)/p_sp(1) ) / float( nspace-1 )
      dpx = exp(dpx)
      do 101 i = 2, nspace-1
         g_sp(i) = g_sp(i-1) * dgx
         p_sp(i) = p_sp(i-1) * dpx
 101  continue

      do 111 i = 1, nspace
         gdis(i) = 0.d0
         pdis(i) = 0.d0
 111  continue


      return
      END
*
*
      SUBROUTINE sumgb( te, gam1 )
      DOUBLE PRECISION te, gam1
* ------------------------------------------------
*   Mesh points for
*   integration over gamma factor for Maxwellian
* ------------------------------------------------
      INTEGER NG
      PARAMETER( NG = 40 )
*
      INTEGER ng_th
      DOUBLE PRECISION gb_th, gbm_th, gbp_th, du
      COMMON /gbmesh/ gb_th(NG), gbm_th(NG), gbp_th(NG), du, ng_th
*      
      INTEGER j
      DOUBLE PRECISION up1, upn, u(NG)
*
      ng_th = NG

      up1 = 1.d0/exp( (gam1-1.d0)/te )
      upn = 1.d0
      du = (upn-up1)/float(ng_th)
      u(1) = up1 + 0.5d0*du 
      do 11 j = 2, ng_th
         u(j) = u(j-1) + du
 11   continue

      do 21 j = 1, ng_th
         gb_th(j)  = 1.d0 - te*log( u(j) )
         gbp_th(j) = gb_th(j) 
     &              + sqrt( (gb_th(j)-1.d0)*(gb_th(j)+1.d0) )
         gbm_th(j) = 1.d0 / gbp_th(j)
 21   continue
*
      return
      END
*
*
      SUBROUTINE initfx
* ----------------------------------------
*       Initialization of flux, etc.
* ----------------------------------------
      INTEGER NX0
      PARAMETER( NX0 = 1026 )
      INTEGER nx
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      INTEGER noph
      DOUBLE PRECISION sinpt
      COMMON /penbin/ sinpt(NX0), noph(NX0)
      DOUBLE PRECISION snj, snjs, dsnj, esbj
      COMMON /statis/ snj(NX0), snjs(NX0), dsnj(NX0), esbj(NX0)
      INTEGER j
*      
      do 12 j = 1, nx
         esbj(j) = 0.d0
         snj(j)  = 0.d0
         snjs(j) = 0.d0
         dsnj(j) = 0.d0
         noph(j) = 0
 12   continue
*
      return
      END
*
*
      SUBROUTINE sums
* -----------------------
*       get flux, etc.
* -----------------------
      INTEGER NX0
      PARAMETER( NX0 = 1026 )
      DOUBLE PRECISION PI
      PARAMETER( PI = 3.1415927D0 )
      INTEGER nx
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION snj, snjs, dsnj, esbj
      COMMON /statis/ snj(NX0), snjs(NX0), dsnj(NX0), esbj(NX0)
      INTEGER j
*
      do 10 j = 1, nx
         snj(j)  = snj(j)  + esbj(j)
         snjs(j) = snjs(j) + esbj(j)*esbj(j)
 10   continue
*
*                                      ! initialize for next step
      do 21 j = 1, nx
         esbj(j) = 0.d0
 21   continue
*
      return
      END
*
*
      SUBROUTINE normal( nphoton, ntotin, snjwa )
      INTEGER nphoton, ntotin
      DOUBLE PRECISION snjwa
* ---------------------------------------
*      normalization of flux, etc.
* ---------------------------------------
      INTEGER NX0
      PARAMETER( NX0 = 1026 )
      DOUBLE PRECISION PI
      PARAMETER( PI = 3.1415927D0 )
      INTEGER nx, noph
      DOUBLE PRECISION xhn, snj, snjs, dsnj, esbj, sinpt
      COMMON /enrbin/ xhn(0:NX0), nx
      COMMON /statis/ snj(NX0), snjs(NX0), dsnj(NX0), esbj(NX0)
      COMMON /penbin/ sinpt(NX0), noph(NX0)
      INTEGER j
      DOUBLE PRECISION flnp, devi
*
*                                                ! input photons
      flnp = float( nphoton )
      do 11 j = 1, nx
         sinpt(j) = float(noph(j)) / flnp
 11   continue
*
      devi = 0.675d0 / sqrt( flnp )
      do 15 j = 1, nx
         snj(j)  = snj(j) / flnp
         snjs(j) = snjs(j) / flnp
         dsnj(j) = devi*sqrt( snjs(j) - snj(j)*snj(j) )
 15   continue
*
      ntotin = 0
      snjwa  = 0.d0
      do 21 j = 1, nx
         ntotin = ntotin + noph(j)
         snjwa  = snjwa  + snj(j)
 21   continue
*
      return
      END
*
*
      SUBROUTINE meanen( trad, x0mean, x1mean, acomp )
      DOUBLE PRECISION trad, x0mean, x1mean, acomp
* ------------------------------------------------
*   Mean photon energies:
*     x0mean ... input photons
*     x1mean ... output photons
* ------------------------------------------------
      INTEGER NX0
      PARAMETER( NX0 = 1026 )
      DOUBLE PRECISION x(NX0)
      INTEGER nx, noph
      DOUBLE PRECISION xhn, snj, snjs, dsnj, esbj, sinpt
      COMMON /enrbin/ xhn(0:NX0), nx
      COMMON /statis/ snj(NX0), snjs(NX0), dsnj(NX0), esbj(NX0)
      COMMON /penbin/ sinpt(NX0), noph(NX0)
      INTEGER j
      DOUBLE PRECISION hh, winp, winpn, wtrn, wtrnn, x0thry
*
      do 11 j = 1, nx-1
         hh = 0.5d0*( log( xhn(j) ) + log( xhn(j+1) ) )
         x(j) = exp( hh )
 11   continue
      x(nx) = xhn(nx)
*
      winp = 0.d0
      winpn = 0.d0
      do 21 j = 1, nx
         winp = winp + sinpt(j)*x(j)
         winpn = winpn + sinpt(j)
 21   continue
      x0mean = winp / winpn
*
      wtrn = 0.
      wtrnn = 0.
      do 22 j = 1, nx
         wtrn = wtrn + snj(j)*x(j)
         wtrnn = wtrnn + snj(j)
 22   continue
      x1mean = wtrn / wtrnn
*
      acomp = x1mean / x0mean
      x0thry = 2.701178d0 * trad
*
* ------------------------------------------------------------------
      write(6,910) x0thry, x0mean, x1mean, acomp
 910  format(1h ,'mean of initial photon energy (theory) =',1pe11.4/
     &       1x,'++ simulation results ++'/
     &       1x ,'mean of initial     photon energy =', 1pe11.4/
     &       1x, 'mean of emission    photon energy =', e11.4/
     &       1x, 'energy amplification factor =',e11.4/)
* ------------------------------------------------------------------
*
      return
      END
*
*
      SUBROUTINE phnxyz
* -------------------------------------------------
*  initial values of photon position and direction
* -------------------------------------------------
      DOUBLE PRECISION PI, TPI, THIRD
      PARAMETER( PI = 3.1415927D0, TPI = 2.D0*PI, THIRD = 1.D0/3.D0 )
*
      REAL rvec(5)
      DOUBLE PRECISION xph, yph, zph, omg01, omg02, omg03
      COMMON /positn/ xph, yph, zph, omg01, omg02, omg03
      DOUBLE PRECISION r, cs, sn, phi, dir, trv
*
* ------------------------------------------------------------
*  Central source
*                               + ( x, y, z )
*                                  unit = radius of a sphere
c      xph = 0.d0
c      yph = 0.d0
c      zph = 0.d0
*                                + propagation direction
c      omg03 = 1.d0
c      omg02 = 0.d0
c      omg01 = 0.d0
*
* ------------------------------------------------------------
*
*  Uniform source
*
      call ranmar( rvec, 5 )
*                               + ( x, y, z )
*                                  unit = radius of a sphere
      r = (dble(rvec(1)))**THIRD
      cs = dble( rvec(2)+rvec(2) ) - 1.d0
      sn = sqrt( (1.d0-cs)*(1.d0+cs) )
      phi = TPI*dble(rvec(3))
      xph = r*sn*cos(phi)
      yph = r*sn*sin(phi)
      zph = r*cs
*                                + propagation direction
*  Random direction
*
      omg03 = dble( rvec(4)+rvec(4) ) - 1.d0
      dir = sqrt( (1.d0-omg03)*(1.d0+omg03) )
      trv = TPI*dble(rvec(5))
      omg02 = dir*sin(trv)
      omg01 = dir*cos(trv)
*
      return
      END
*
*
      SUBROUTINE scattr( amxwl, te, p, gam1, gam2, taues, x0, epsct,
     &                   gsp, fdens, y2, n2 )
      INTEGER n2
      DOUBLE PRECISION amxwl, te, p, gam1, gam2, taues, x0, epsct
      DOUBLE PRECISION gsp(1000), fdens(1000), y2(1000) 
* -------------------------------------------------------------------
*            propagation and scattering
* -------------------------------------------------------------------
      INTEGER NX0
      PARAMETER( NX0 = 1026 )
      INTEGER NG
      PARAMETER( NG = 40 )
*
      DOUBLE PRECISION xph, yph, zph, omg01, omg02, omg03
      COMMON /positn/ xph, yph, zph, omg01, omg02, omg03
      INTEGER nx, noph
      DOUBLE PRECISION xhn, snj, snjs, dsnj, esbj, sinpt, w, escwti
      COMMON /enrbin/ xhn(0:NX0), nx
      COMMON /statis/ snj(NX0), snjs(NX0), dsnj(NX0), esbj(NX0)
      COMMON /penbin/ sinpt(NX0), noph(NX0)
      COMMON /weight/ w(2), escwti
      INTEGER ng_th
      DOUBLE PRECISION gb_th, gbm_th, gbp_th, du
      COMMON /gbmesh/ gb_th(NG), gbm_th(NG), gbp_th(NG), du, ng_th

      REAL gxi
      INTEGER j
      DOUBLE PRECISION tauinv, hph, rlmean, rovecs, rleng2,
     &     freefr, fracli, rlamph, x0new, xob, pathli,
     &     omgp1, omgp2, omgp3
      DOUBLE PRECISION gmaxa
*
      tauinv = 1.d0 / taues
      w(1) = 1.d0
*
      do 5 j = 1, nx
         if( x0 .lt. xhn(j+1) ) then
            noph(j) = noph(j) + 1
            go to 6
         end if
 5    continue
 6    continue
*
 50   continue
        hph = x0 + x0
*                                         ! mean free path
        gmaxa = max( 2.d0*gam2*hph, 2.d0*gb_th(1)*hph )
        if( gmaxa .le. 0.01d0 ) then
            rlmean = tauinv
          else
            call meanfr( amxwl, te, taues, hph, rlmean )
        end if
*                                         ! distance to boundary

        rovecs = xph*omg01 + yph*omg02 + zph*omg03
        rleng2 = xph*xph + yph*yph + zph*zph
        pathli = -rovecs + sqrt( 1.d0 - rleng2 + rovecs*rovecs )
        freefr = pathli / rlmean
        fracli = exp( -freefr )
        escwti = w(1)*fracli
*
        if( x0 .ge. xhn(nx) ) go to 1100
        do 111 j = 1, nx-1
           if( x0 .lt. xhn( j + 1 ) ) then
              esbj(j) = esbj(j) + escwti
              go to 112
	   end if   
 111    continue
 112    continue
*
        w(2) = w(1) - escwti
*
        if( w(2) .lt. epsct ) go to 1100
*                                       ! proceed to scattering point
        call ranmar( gxi, 1 )
*
        rlamph = -rlmean*log( 1.d0 - dble(gxi)*( 1.d0 - fracli ) )
        xph = xph + rlamph*omg01
        yph = yph + rlamph*omg02
        zph = zph + rlamph*omg03
*                                                       ! scattering
*
*                                       ! select electron momentum
        call slcmom( hph, xob, gsp, fdens, y2, n2 )
*
*                                       ! get a new photon momentum
        call newphn( x0, xob, omgp1, omgp2, omgp3, x0new )
        x0    = x0new
        omg01 = omgp1
        omg02 = omgp2
        omg03 = omgp3
        w(1)  = w(2)
*
      go to 50
*
 1100 continue
*
      return
      END
*
*
      SUBROUTINE phengy( trad, x0 )
      DOUBLE PRECISION trad, x0
* -------------------------------------
*  selection of initial photon energy
*  + Planckian with temperature of trad
* -------------------------------------
      INTEGER MDIM
      PARAMETER( MDIM = 200 )
*
      DOUBLE PRECISION sm
      COMMON /SUMJ3I/ sm(MDIM)
      REAL rvec(4)
      INTEGER i
      DOUBLE PRECISION z3xi
*
 1    call ranmar( rvec, 4 )
*
      z3xi = 1.202d0 * dble(rvec(1))
      if( z3xi .lt. 1.d0 ) then
        x0 = -trad*log( dble(rvec(2)*rvec(3)*rvec(4)) )
        return
      end if
*
      i = 1
 10   i = i + 1
      if( i    .ge. MDIM  ) go to 1
      if( z3xi .ge. sm(i) ) go to 10
      x0 = -trad / float(i) * log( dble(rvec(2)*rvec(3)*rvec(4)) )
*
      return
      END
*
*
      SUBROUTINE meanfr( amxwl, te, taues, hph, rlmean )
      DOUBLE PRECISION amxwl, te, taues, hph, rlmean
* --------------------------------------------------
*   rlmean ... mean free path of a photon
*              for Thermal and Non-Thermal Electrons
* --------------------------------------------------
      DOUBLE PRECISION PI, SIXI, GCOF
      PARAMETER( PI = 3.1415927D0, SIXI = 1.D0/6.D0 )
      PARAMETER( GCOF = 2.D0/3.D0/PI )
      INTEGER NGDIM
      PARAMETER( NGDIM = 128 )
      INTEGER NG
      PARAMETER( NG = 40 )

      INTEGER ng_nth
      DOUBLE PRECISION gi, gbm, gbp, disfc, delg
      COMMON /gint/ gi(NGDIM), gbm(NGDIM), gbp(NGDIM), disfc(NGDIM),
     &              delg, ng_nth
      INTEGER ng_th
      DOUBLE PRECISION gb_th, gbm_th, gbp_th, du
      COMMON /gbmesh/ gb_th(NG), gbm_th(NG), gbp_th(NG), du, ng_th

      DOUBLE PRECISION th_nrm
      COMMON /maxnrm/ th_nrm

      INTEGER j
      DOUBLE PRECISION gghh, sphi, gx, gx2, phil, phiu, gxo
      DOUBLE PRECISION ri_th, ri_nth
*
      gghh = GCOF*hph*hph/taues
*
      if( amxwl .lt. 1.d-8 ) then
         ri_th = 0.d0
         go to 666
      end if
*                                                    Thermal part
*
      sphi = 0.d0
      DO 11 j = 1, ng_th
         gx = hph * gbm_th(j)
         if( gx .le. 0.5d0 ) then
            gx2  = gx * gx
            phil = ( SIXI + 0.047d0*gx - 0.03d0*gx2
     &             + 0.5d0/(1.d0+gx) )*gx2
           else if( gx .gt. 3.5d0 ) then
            gxo  = 1.d0 + gx
            phil = gxo*log( gxo ) - 0.5d0*gx
     &            - 13.16d0*log( 2.d0+0.076d0*gx ) + 9.214d0
           else
            gxo  = 1.d0 + gx
            phil = gxo*log(gxo) - 0.94d0*gx - 0.00925d0
         end if
*
         gx = hph * gbp_th(j)
         if( gx .le. 0.5d0 ) then
            gx2  = gx * gx
            phiu = ( SIXI + 0.047d0 * gx - 0.03d0*gx2
     &             + 0.5d0/(1.d0+gx) )*gx2
           else if( gx .gt. 3.5d0 ) then
            gxo  = 1.d0 + gx
            phiu = gxo*log( gxo ) - 0.5d0*gx
     &             - 13.16d0*log( 2.d0+0.076d0*gx ) + 9.214d0
           else
            gxo  = 1.d0 + gx
            phiu = gxo*log( gxo ) - 0.94d0*gx - 0.00925d0
         end if
*
         sphi = sphi + phiu - phil
 11   continue
      sphi = sphi * du
*
      ri_th = th_nrm * sphi
*
*
*                                                    Non-thermal part
      if( abs(amxwl-1.d0) .lt. 1.d-8 ) then
         ri_nth = 0.d0
         go to 888
      end if

 666  continue
      sphi = 0.d0
      do 100 j = 1, ng_nth
         gx = hph * gbm(j)
         if( gx .le. 0.5d0 ) then
            gx2  = gx*gx
            phil = ( SIXI + 0.047d0*gx - 0.03d0*gx2
     &             + 0.5d0/( 1.d0 + gx ) )*gx2
           else if( gx .gt. 3.5d0 ) then
            gxo  = 1.d0 + gx
            phil = gxo*log( gxo ) - 0.5d0*gx
     &             - 13.16d0*log( 2.d0 + 0.076d0*gx ) + 9.214d0
           else
            gxo  = 1.d0 + gx
            phil = gxo*log( gxo ) - 0.94d0*gx - 0.00925d0
         end if
*
         gx = hph * gbp(j)
         if( gx .le. 0.5d0 ) then
            gx2  = gx*gx
            phiu = ( SIXI + 0.047d0*gx - 0.03d0*gx2
     &             + 0.5d0/( 1.d0 + gx ) )*gx2
          else if( gx .gt. 3.5d0 ) then
            gxo  = 1.d0 + gx
            phiu = gxo*log( gxo ) - 0.5d0*gx
     &             - 13.16d0*log( 2.d0 + 0.076d0*gx ) + 9.214d0
          else
            gxo  = 1.d0 + gx
            phiu = gxo*log( gxo ) - 0.94d0*gx - 0.00925d0
         end if
*
         sphi = sphi + ( phiu - phil )*disfc(j)
 100  continue
      sphi = sphi * delg
      ri_nth = sphi
      
 888  rlmean = gghh / (amxwl*ri_th + (1.d0-amxwl)*ri_nth)
*
      return
      END
*
*
      SUBROUTINE slcmom( hph, xob, gsp, fdens, y2, n2 )
      INTEGER n2
      DOUBLE PRECISION hph, xob
      DOUBLE PRECISION gsp(1000), fdens(1000), y2(1000) 
* ---------------------------------------------------------
*   selection of electron momentum
*   isotropic, relativistic MAXWELLIAN with temperature te
*   isotropic, relativistic non-thermal distribution
* ---------------------------------------------------------
      DOUBLE PRECISION PI, THRI, TPI
      PARAMETER ( PI = 3.14159265359D0, THRI = 1.D0/3.D0, 
     &            TPI = 2.D0*PI )
*
      REAL rvec(4)
      DOUBLE PRECISION xph, yph, zph, omg01, omg02, omg03
      COMMON /POSITN/ xph, yph, zph, omg01, omg02, omg03
      DOUBLE PRECISION gmma, voc, rmu, v01, v02, v03
      COMMON /NEWENG/ gmma, voc, rmu, v01, v02, v03

      INTEGER nspace
      DOUBLE PRECISION g_sp, p_sp, gdis, pdis
      COMMON /used_e/ g_sp(1000), p_sp(1000), gdis(1000), pdis(1000),
     &                nspace

      INTEGER i
      DOUBLE PRECISION v3s, tx2, rmv, xob2, sigmh
      DOUBLE PRECISION p_mom
*
 1    call ranmar( rvec, 4 )
      v03 = dble( rvec(1) + rvec(1) ) - 1.d0
      v3s = sqrt( (1.d0-v03)*(1.d0+v03) )
      tx2 = TPI * dble(rvec(2))
      v02 = v3s * sin( tx2 )
      v01 = v3s * cos( tx2 )
*
*                                        ! use spline interpolation
      call splint( fdens, gsp, y2, n2, dble(rvec(3)), gmma )
*
      voc = sqrt( (gmma-1.d0)*(gmma+1.d0) ) / gmma
      rmu = v01*omg01 + v02*omg02 + v03*omg03
      rmv = 1.d0 - rmu*voc
      xob = hph*gmma*rmv
      if( xob .le. 0.5d0 ) then
        xob2  = xob * xob
        sigmh = THRI + 0.141d0*xob - 0.12d0*xob2
     &          + ( 1.d0 + 0.5d0*xob )/( 1.d0 + xob + xob + xob2 )
       else if( xob .gt. 3.5d0 ) then
        sigmh = ( log( 1.d0 + xob ) + 0.5d0
     &         - 1.d0/( 2.d0 + 0.076d0*xob ) )/xob
       else
        sigmh = ( log( 1.d0 + xob ) + 0.06d0 )/xob
      end if
*
      if( dble(rvec(4)) .ge. 0.375d0*sigmh*rmv ) go to 1
*
*               
*                                           Add selected electron
      do 121 i = 1, nspace-1
         if( gmma .ge. g_sp(i) .and. gmma .lt. g_sp(i+1) ) then
            gdis(i+1) = gdis(i+1) + 1.d0
            go to 122
         end if
 121  continue
 122  continue
      do 131 i = 1, nspace-1
         p_mom = sqrt( gmma*gmma-1.d0 )
         if( p_mom .ge. p_sp(i) .and. p_mom .lt. p_sp(i+1) ) then
            pdis(i+1) = pdis(i+1) + 1.d0
            go to 132
         end if
 131  continue
 132  continue
*
      return
      END
*
*
      SUBROUTINE newphn( x0, xob, omgp1, omgp2, omgp3, x0new )
      DOUBLE PRECISION x0, xob, omgp1, omgp2, omgp3, x0new
* ------------------------------------------------------------
*            scattered photon energy and direction
* ------------------------------------------------------------
      DOUBLE PRECISION PI, TPI
      PARAMETER( PI = 3.1415927D0, TPI = 2.D0*PI )
*
      DOUBLE PRECISION xph, yph, zph, omg01, omg02, omg03
      DOUBLE PRECISION gmma, voc, rmu, v01, v02, v03
      COMMON /positn/ xph, yph, zph, omg01, omg02, omg03
      COMMON /neweng/ gmma, voc, rmu, v01, v02, v03
      REAL rvec(3)
      DOUBLE PRECISION tgx1, rmup, rho1, rmus, rmri, phip, cphi, sphi, 
     &     v3sph, omgvc, rmupv, xpoxi, xpox, xp, 
     &     xis, xx, yy
*
 10   call ranmar( rvec, 3 )
      tgx1 = dble( rvec(1)+rvec(1) ) - 1.d0
      rmup = ( voc + tgx1 ) / ( 1.d0 + voc*tgx1 )
      rho1 = sqrt( v01*v01 + v02*v02 )
      rmus = sqrt( 1.d0 - rmup*rmup )
      rmri = rmus / rho1
      phip = TPI  * dble(rvec(2))
      cphi = cos( phip )
      sphi = sin( phip )
      v3sph = v03 * sphi
*
      omgp1 = rmup*v01 + rmri*(  v02*cphi + v01*v3sph )
      omgp2 = rmup*v02 + rmri*( -v01*cphi + v02*v3sph )
      omgp3 = rmup*v03 - rmus*rho1*sphi
*
      omgvc = omg01*omgp1 + omg02*omgp2 + omg03*omgp3
      rmupv = 1.d0 - rmup * voc
      xpoxi = 1.d0 + x0 / gmma * ( 1.d0 - omgvc ) / rmupv
      xpox  = 1.d0 / xpoxi
      xp  = xpox*xob
      xis = 1.d0/xob - 1.d0/xp
      xx = xpox + xpoxi + 4.d0*xis*( 1.d0 + xis )
      yy = xpox * xpox * xx
*
      if( dble(rvec(3)+rvec(3)) .ge. yy ) go to 10
*
      x0new = 0.5d0*xp/gmma/rmupv
*
      return
      END
*
*
      SUBROUTINE sumofm
* ----------------------
*      sum of 1/j**3
* ----------------------
      INTEGER MDIM
      PARAMETER( MDIM = 200 )
*
      DOUBLE PRECISION sm
      COMMON /sumj3i/ sm(MDIM)
      INTEGER j
      DOUBLE PRECISION dj
*
      sm(1) = 1.d0
      do 10 j = 2, MDIM
        dj = float(j)
        sm(j) = sm(j-1) + 1.d0/( dj*dj*dj )
 10   continue
*
      return
      END
*
*
      SUBROUTINE f_gmin( te, amxwl, p, gam2, gam1 )
      DOUBLE PRECISION te, amxwl, p, gam2, gam1
* --------------------------------------------------
*  Obtain gam1 = gamma_min
*  such that Maxwellian = Power-law at gamma = gam1
* --------------------------------------------------
      DOUBLE PRECISION rmaxn, rnon, theta, pind, g2
      COMMON /toequ/ rmaxn, rnon, theta, pind, g2
      DOUBLE PRECISION zbrent, equi0, equi1
      EXTERNAL equi0, equi1

      rmaxn = amxwl/te
      theta = te

      if( abs(p-1.d0) .gt. 1.d-8 ) then
         rnon = (1.d0-amxwl)*(1.d0-p)
         pind = p
         g2 = gam2
         gam1 = zbrent( equi0, 1.d0+te, g2/1.2d0, 1.d-7 )
      else
         rnon = 1.d0 - amxwl
         g2 = gam2
         gam1 = zbrent( equi1, 1.d0+te, g2/1.2d0, 1.d-7 )
      end if
      return
      END

      FUNCTION equi0(g)
      DOUBLE PRECISION equi0, g
      DOUBLE PRECISION rmaxn, rnon, theta, pind, g2
      COMMON /toequ/ rmaxn, rnon, theta, pind, g2
      INTEGER n, i
      DOUBLE PRECISION up1, upn, u(100), du, gi
      DOUBLE PRECISION thnorm, sss, pone
*
      n = 100
      up1 = 1.d0/exp( (g-1.d0)/theta )
      upn = 1.d0
      du = (upn-up1)/float(n)
      u(1) = up1 + 0.5d0*du 
      do 11 i = 2, n
         u(i) = u(i-1) + du
 11   continue
      sss = 0.d0
      do 21 i = 1, n
         gi = 1.d0 - theta*log(u(i))
         sss = sss + gi*sqrt((gi-1.d0)*(gi+1.d0))
 21   continue
      sss = sss * du
      thnorm = rmaxn / sss 

      pone = 1.d0 - pind
      equi0 = thnorm*g*sqrt((g-1.d0)*(g+1.d0))*exp(-(g-1.d0)/theta)
     &       - rnon/(g2**pone-g**pone)/g**pind
      return
      END
*
      FUNCTION equi1(g)
      DOUBLE PRECISION equi1, g
      DOUBLE PRECISION rmaxn, rnon, theta, pind, g2
      COMMON /toequ/ rmaxn, rnon, theta, pind, g2
      INTEGER n, i
      DOUBLE PRECISION up1, upn, u(100), du, gi
      DOUBLE PRECISION thnorm, sss
*
      n = 100
      up1 = 1.d0/exp( (g-1.d0)/theta )
      upn = 1.d0
      du = (upn-up1)/float(n)
      u(1) = up1 + 0.5d0*du 
      do 11 i = 2, n
         u(i) = u(i-1) + du
 11   continue
      sss = 0.d0
      do 21 i = 1, n
         gi = 1.d0 - theta*log(u(i))
         sss = sss + gi*sqrt((gi-1.d0)*(gi+1.d0))
 21   continue
      sss = sss * du
      thnorm = rmaxn / sss 

      equi1 = rmaxn*g*sqrt((g-1.d0)*(g+1.d0))*exp(-(g-1.d0)/theta)
     &       - rnon/log(g2/g)/g
      return
      END
*
      SUBROUTINE thnrml( te, gam1, thnrm0 )
      DOUBLE PRECISION te, gam1, thnrm0
* ------------------------------------------
*  Normalization of thermal electrons
* ------------------------------------------
      DOUBLE PRECISION PI, FPI
      PARAMETER( PI = 3.1415927D0 )
      PARAMETER( FPI = 4.D0*PI )
      INTEGER n, i
      DOUBLE PRECISION up1, upn, u(100), du, gi, sss

      n = 100
      up1 = 1.d0/exp( (gam1-1.d0)/te )
      upn = 1.d0
      du = (upn-up1)/float(n)
      u(1) = up1 + 0.5d0*du 
      do 11 i = 2, n
         u(i) = u(i-1) + du
 11   continue
      sss = 0.d0
      do 15 i = 1, n
         gi = 1.d0 - te*log(u(i))
         sss = sss + gi*sqrt((gi-1.d0)*(gi+1.d0))
 15   continue
      sss = sss * du

      thnrm0 = exp(1.d0/te)/(FPI*te*sss)
      return
      END
*
      SUBROUTINE table( amxwl, thnrm0, c_nom, te, gam1, gam2, p,
     &                  gsp, fdens, y2, n2 )
      INTEGER n2
      DOUBLE PRECISION amxwl, thnrm0, c_nom, te, gam1, gam2, p 
      DOUBLE PRECISION gsp(1000), fdens(1000), y2(1000) 
* ----------------------------------------------------------------
*  Make a table for intermpolation of momentum-probability
* ----------------------------------------------------------------
      DOUBLE PRECISION yp1, ypn, dg, delg1, delg2
      DOUBLE PRECISION fth, fnt
      INTEGER i
      INTEGER n1
* 
      n1 = 500
      n2 = 1000
      gsp(1) = 1.d0
      gsp(n1) = gam1
      gsp(n2) = gam2
      delg1 = dlog(gsp(n1)/gsp(1))/float(n1-1)
      dg = dexp(delg1)
      do 3 i = 2, n1-1
         gsp(i) = gsp(i-1)*dg
 3    continue
      delg2 = dlog(gsp(n2)/gsp(n1))/float(n2-n1)
      dg = dexp(delg2)
      do 5 i = n1+1, n2-1
         gsp(i) = gsp(i-1)*dg
 5    continue
*
      fdens(1) = 0.d0
      do 51 i = 2, n1
         fdens(i) = fdens(i-1) 
     &             + amxwl*thnrm0*fth(gsp(i),te)*gsp(i)*delg1
 51   continue
      do 53 i = n1+1, n2
         fdens(i) = fdens(i-1) 
     &       + (1.d0-amxwl)*c_nom*fnt(gsp(i),p,gam1,gam2)*gsp(i)*delg2
 53   continue
      do 55 i = 1, n2
         fdens(i) = fdens(i) / fdens(n2)
 55   continue

      yp1 = (gsp(2)-gsp(1))/(fdens(2)-fdens(1))
      ypn = (gsp(n2)-gsp(n2-1))/(fdens(n2)-fdens(n2-1))
      call spline(fdens, gsp, n2, yp1, ypn, y2 )
*
      return
      END
*
      FUNCTION fth( g, te )
      DOUBLE PRECISION fth, g, te
* ------------------------------------------
*   Maxwellian distribution for a given Te
* ------------------------------------------
      fth = g*sqrt((g-1.d0)*(g+1.d0))*exp(-g/te)
      return
      END
*
      FUNCTION fnt( g, p, gam1, gam2 )
      DOUBLE PRECISION fnt, g, p, gam1, gam2
* ---------------------------------------------------------
*   Nonthermal distribution
* ---------------------------------------------------------
      if( g .ge. gam1 .and. g .le. gam2 ) then
         fnt = 1.d0/g**p
      else
         fnt = 0.d0
      end if
      return
      END
*
      SUBROUTINE ranmar( rvec, len )
      INTEGER len
      REAL rvec(len)
*
* Universal random number generator proposed by MARSAGLIA and ZAMAN
* in REPORT FSU-SCRI-87-50
* Statistics and Probability Letters Vol.9 (1990) 35-39
* NORTH-HOLLAND
*
*  @ Random Number:  0 < RAND NUM < 1
*    (more exactly, NUMBER > = 10**(-48) ) 
*  @ Period = 2**144 = 2 * 10**43
*
*   Slightly modified by F. JAMES, 1988, to generate a vector
*   of pseudorandom numbers rvec of length len
*   and making the common blocks include everything needed to
*   specify completely the sate of the generator.
*
      INTEGER i97, j97
      REAL u, c, cd, cm
      COMMON /raset1/ u(97), c, cd, cm, i97, j97
      REAL uni
      INTEGER ivec
*
      do 100 ivec = 1, len
        uni = u(i97) - u(j97)
        if( uni .lt. 0. ) uni = uni + 1.
        u(i97) = uni
        i97 = i97 - 1
        if( i97 .eq. 0 ) i97 = 97
        j97 = j97 - 1
        if( j97 .eq. 0 ) j97 = 97
        c = c - cd
        if( c .lt. 0. ) c = c + cm
        uni = uni - c
        if( uni .lt. 0. ) uni = uni + 1.
*                                         ! to make uni not equal zero
        if( uni .eq. 0. ) then
           uni = u(j97) * 2.**(-24)
        end if
*
        rvec(ivec) = uni
 100  continue
      return
      END
*
      SUBROUTINE rmarin(ijkl)
      INTEGER ijkl
*
* Initializing routine for ranmar, must be called before 
* generating any pseudorandom numbers with ranmar.
* The input value sould be in the range: 0 < = ijkl < = 900 000 000
*
      INTEGER i97, j97
      REAL u, c, cd, cm
      COMMON /raset1/ u(97), c, cd, cm, i97, j97
*
* This shows correspondence between the simplified input seed ijkl
* and the original marsaglia-zaman paper,
* ( i = 12, j = 34, k = 56, l 78 ) put ijkl = 5421713
*
      INTEGER ij, kl, i, j, k, l, m, ii, jj
      REAL s, t

      ij = ijkl / 30082
      kl = ijkl - 30082 * ij
      i = mod( ij / 177, 177 ) + 2
      j = mod( ij, 177 ) + 2
      k = mod( kl / 169, 178 ) + 1
      l = mod( kl, 169 )
      write(6,12) ijkl, i,j,k,l
   12 format(1h ,'ranmar initialized:', i15, 2x, 4i4//)
      do 2 ii = 1, 97
         s = 0.
         t = 0.5
         do 3 jj = 1, 24
            m = mod( mod( i * j, 179 ) * k, 179 )
            i = j
            j = k
            k = m
            l = mod( 53 * l + 1, 169 )
            if( mod( l * m ,64 ) .ge. 32 ) s = s + t
            t = 0.5 * t
 3       continue
         u(ii) = s
 2    continue
      c  = 362436.   / 16777216.
      cd = 7654321.  / 16777216.
      cm = 16777213. / 16777216.
      i97 = 97
      j97 = 33
      return
      END
*
*
      FUNCTION rk2e(x)
      DOUBLE PRECISION rk2e, x
* -------------------------------------------------------
*  Modified Bessel function of order 2, scaled by exp(x)
*
*   rk2e(x) = exp(x) K_2(x)
* -------------------------------------------------------     
      DOUBLE PRECISION  beske0, beske1
      rk2e = beske0(x) + 2.d0/x * beske1(x)
      return
      END
*
      FUNCTION beske0(x)
      DOUBLE PRECISION beske0,x
CU    USES bessi0
      DOUBLE PRECISION bessi0
      DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
      DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0,
     *0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/
      DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1,
     *-0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/
      if (x.le.2.0d0) then
        y=x*x/4.0d0
        beske0=(-log(x/2.0d0)*bessi0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*
     *(p6+y*p7)))))) * exp(x)
      else
        y=(2.0d0/x)
        beske0=(1.d0/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     *q7))))))
      endif
      return
      END
*
      FUNCTION beske1(x)
      DOUBLE PRECISION beske1,x
CU    USES bessi1
      DOUBLE PRECISION bessi1
      DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
      DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,0.15443144d0,-0.67278579d0,
     *-0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/
      DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1,
     *0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/
      if (x.le.2.0d0) then
        y=x*x/4.0d0
        beske1=(log(x/2.0d0)*bessi1(x))+(1.0d0/x)*(p1+y*(p2+y*(p3+y*(p4
     *       +y*
     *       (p5+y*(p6+y*p7)))))) * exp(x)
      else
        y=2.0d0/x
        beske1=(1.d0/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     *       q7))))))
      endif
      return
      END
*
      FUNCTION bessi0(x)
      DOUBLE PRECISION bessi0,x
      DOUBLE PRECISION ax
      DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
      DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,
     *     1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/
      DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1,
     * 0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,0.2635537d-1,
     * -0.1647633d-1,0.392377d-2/
      if (abs(x).lt.3.75d0) then
        y=(x/3.75d0)**2
        bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
      else
        ax=abs(x)
        y=3.75d0/ax
        bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     *       (q7+y*(q8+y*q9))))))))
      endif
      return
      END
*
      FUNCTION bessi1(x)
      DOUBLE PRECISION bessi1,x
      DOUBLE PRECISION ax
      DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
      DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,
     *0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
      DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1,
     *-0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,-0.2895312d-1,
     *0.1787654d-1,-0.420059d-2/
      if (abs(x).lt.3.75d0) then
        y=(x/3.75d0)**2
        bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
      else
        ax=abs(x)
        y=3.75d0/ax
        bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     *       (q7+y*(q8+y*q9))))))))
        if(x.lt.0.d0)bessi1=-bessi1
      endif
      return
      END
*      
      FUNCTION zbrent(func,x1,x2,tol)
      INTEGER ITMAX
      DOUBLE PRECISION zbrent,tol,x1,x2,func,EPS
      EXTERNAL func
      PARAMETER (ITMAX=100,EPS=3.d-16)
      INTEGER iter
      DOUBLE PRECISION a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
      a=x1
      b=x2
      fa=func(a)
      fb=func(b)
      if((fa.gt.0.d0.and.fb.gt.0.d0).or.(fa.lt.0.d0.and.fb.lt.0.d0))
     * then
         write(*,*) 'root must be bracketed for zbrent: no gam1 found'
         write(*,*) 'amxwl should be changed'
         pause
      end if
      c=b
      fc=fb
      do 11 iter=1,ITMAX
        if((fb.gt.0.d0.and.fc.gt.0.d0).or.(fb.lt.0.d0.and.fc.lt.0.d0)
     *)then
          c=a
          fc=fa
          d=b-a
          e=d
        endif
        if(abs(fc).lt.abs(fb)) then
          a=b
          b=c
          c=a
          fa=fb
          fb=fc
          fc=fa
        endif
        tol1=2.d0*EPS*abs(b)+0.5d0*tol
        xm=.5d0*(c-b)
        if(abs(xm).le.tol1 .or. fb.eq.0.d0)then
          zbrent=b
          return
        endif
        if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
          s=fb/fa
          if(a.eq.c) then
            p=2.d0*xm*s
            q=1.d0-s
          else
            q=fa/fc
            r=fb/fc
            p=s*(2.d0*xm*q*(q-r)-(b-a)*(r-1.d0))
            q=(q-1.d0)*(r-1.d0)*(s-1.d0)
          endif
          if(p.gt.0.d0) q=-q
          p=abs(p)
          if(2.d0*p .lt. min(3.d0*xm*q-abs(tol1*q),abs(e*q))) then
            e=d
            d=p/q
          else
            d=xm
            e=d
          endif
        else
          d=xm
          e=d
        endif
        a=b
        fa=fb
        if(abs(d) .gt. tol1) then
          b=b+d
        else
          b=b+sign(tol1,xm)
        endif
        fb=func(b)
11    continue
      pause 'zbrent exceeding maximum iterations'
      zbrent=b
      return
      END
*
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      DOUBLE PRECISION yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=5000)
      INTEGER i,k
      DOUBLE PRECISION p,qn,sig,un,u(NMAX)
      if (yp1.gt..99d30) then
        y2(1)=0.d0
        u(1)=0.d0
      else
        y2(1)=-0.5d0
        u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.d0
        y2(i)=(sig-1.d0)/p
        u(i)=(6.d0*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99d30) then
        qn=0.d0
        un=0.d0
      else
        qn=0.5d0
        un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END
*
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      DOUBLE PRECISION x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      DOUBLE PRECISION a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.d0) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.d0
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software ]m!+)+535.1n.




