Subject: disk code

From: [email-masked]
Date: Thu, 6 Jun 1996 16:09:21 +0900

*
*               --------------------------------------------
*               MONTE CARLO SIMULATION OF COMPTON SCATTERING
*               --------------------------------------------
*
*  @ GEOMETRY ... DISK ( SLAB )
*                 THICKNESS OF A DISK ... TAUES
*
*  @ ELECTRONS ... Isotropic, Nonthermal
*  @ INPUT PHOTONS ... Planckian
*
* + INPUT :
*
*     alpha: power-law index of nonthermal electrons ... gamma**(-alpha)
*     gam1:  low-energy cutoff of electron's Lorentz factor
*     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
*                from one side of a disk
*     fftrn(j) : relative value of transmitted energy flux from a disk
*     ffref(j) : relative value of reflected energy flux from a disk
*
*     mean photon enregies :
*         x0mean ... input  photon energy
*         x1mean ... output photon energy
*         x1trnm ... transmitted photon energy
*         x1refm ... reflected photon energy
*
*  SYMBOLS :
*     xhn(j) ... photon energy in units of electron mass
*                (dimension = NX0)
*     rmu(m) ... mesh of cos (theta) (dimension 2*mu+1)
*     w(i) ... weight
*     --  --  --  --  --  --  --  --  --  -- --  --  --  --  --
*
*      V. 1.0   BY M. KUSUNOSE   1/1/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, low_gam, hi_gam, tradch, tauch, nch, epsch
      INTEGER NX0, MU, NGDIM
      PARAMETER( NX0 = 1026, MU = 10 )
      PARAMETER( NGDIM = 128 )
*
      DOUBLE PRECISION EMASS, VELC, BOLTZK, SBC
      PARAMETER( EMASS = 9.109534D-28, VELC = 2.99792458D10,
     &     BOLTZK = 1.3807D-16, SBC = 5.67051D-05 )
*
*                                    ! EMASS ... electron mass
*                                    ! VELC ... light speed
*                                    ! BOLTZK ... Boltzmann constant
*                                    ! SBC ... Stefan-Boltzmann constant
*
      INTEGER nx
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      INTEGER noph
      DOUBLE PRECISION esbb, sinpt
      COMMON /penbin/ esbb(NX0,-MU:2*MU+1), sinpt(NX0), noph(NX0)
      DOUBLE PRECISION fftrn, ffref
      COMMON /flaptn/ fftrn(NX0), ffref(NX0)
      INTEGER ng
      DOUBLE PRECISION gi, gbm, gbp, disfc, delg
      COMMON /gint/ gi(NGDIM), gbm(NGDIM), gbp(NGDIM), disfc(NGDIM),
     &              delg, ng
      DOUBLE PRECISION sum_g
      COMMON /sum_el/ sum_g(NGDIM)

      INTEGER ijkl, i, j, nphoton, ntotin
      REAL power_t, power_r, sigpi_t, sigpi_r
      DOUBLE PRECISION alpha, gam1, gam2, tradkv, taues, epsct, trad,
     &     tradk, xlow, xhigh, dxl, c_nom, x0, ftjwa, frjwa,
     &     x0mean, x1mean, x1trnm, x1refm, acomp,
     &     e_min, e_max, elnorm
*
*
*  ! Input parameters
*
      open(8,file='input1',status='old')
        read(8,*) index,   alpha
        read(8,*) low_gam, gam1
        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
*  ------------------------------------------------
      write(6,905)
      write(6,910) alpha, gam1, gam2,
     &             trad, tradkv, tradk, 
     &             taues, nphoton, epsct
*  ------------------------------------------------
*
*
*                                  ! mesh of photon energy and angle
      xlow = 1.d-13
      xhigh = 1.d+04
      dxl = 0.2d0
      call mesh( xlow, xhigh, dxl )
*
*                             ! normalization of electron distribution
      call el_nml( alpha, gam1, gam2, c_nom )
*
*                             ! mesh in electron Lorentz factor
      call g_mesh( alpha, gam1, gam2, c_nom )
      do 21 i = 1, ng
         sum_g(i) = 0.d0
 21   continue
*
*                                  ! initialize random number generator
* The input value should be in the range: 0 <= ijkl <= 900 000 000
*
      ijkl = 199609
      call rmarin(ijkl)
*
*  -------------- beginning of scattering ---------------
*
      call sumofm
      call initfx
      do 51 i = 1, nphoton
         if( mod(i,1000) .eq. 0 ) then
            open(22,file='number',status='unknown')
            write(22,*) 'number =',i
            close(22)
         end if
         call phnxyz
         call phengy( trad, x0 )
         call scattr( alpha, gam1, gam2, taues, x0, epsct )
         call sums
 51   continue
*
*  ---------------- end of scattering --------------------
*
*                                                     ! normalization
      call normal( nphoton, ntotin, ftjwa, frjwa )
*
*                                              ! mean photon energies
      call meanen( trad, x0mean, x1mean, x1trnm, x1refm, acomp )
*
*
*                                              ! angular distribution
c      call angular
*
*                                              ! power-law index
*                                                of photon spectrum
      e_min = 5.d0
      e_max = 20.d0
      call sp_index( e_min, e_max, power_t, power_r, 
     &               sigpi_t, sigpi_r )
*
*  ---------------------------------------------
*
      open(12,file='spect1.dat',status='unknown')
        write(12,941) ( xhn(j), sinpt(j), fftrn(j), ffref(j),
     &                  j = 1, nx )
 941    format(1h ,1pe11.4,2x,e11.4,2x,e11.4,2x,e11.4)
      close(12)

      elnorm = 0.d0
      do 111 i = 1, ng
         elnorm = elnorm + sum_g(i)
 111  continue
      do 121 i= 1, ng
         sum_g(i) = sum_g(i) / elnorm
 121  continue
      open(14,file='electron.dat',status='unknown')
        write(14,950) (gi(i), sum_g(i), i = 1, ng)
 950    format(1h ,1pe12.5,2x,e12.5)
      close(14)
*
*  ---------------------------------------------
*
 905  format(1h ,5x,'COMPTON SCATTERING BY RELATIVISTIC ELECTRONS'/
     &    1x,'   + Geometry ... Disk with tau_es'/
     &    1x,'   + Electrons ... Nonthermal'/
     &    1x,'   + Input photons ... Planckian with T_rad'// )
 910  format(1h ,'alpha =',1pe11.4,2x,'gam_low =',e11.4,
     &  2x,'gam_high =',e11.4/
     &  1x,'kT_rad/mc^2 =',e11.4,2x,'T_rad(keV) =',e11.4,2x,
     &  'T_rad(K) =',e11.4/
     &  1x,'tau_es =',e11.4/
     &  1x,'Photon Number =',i7,2x,'eps(cutoff) =',e10.3/)
 915  format(1h ,'compton cooling rate per unit column =',1pe11.4/)
*
      stop
      END
*
*
      SUBROUTINE el_nml( alpha, gam1, gam2, c_nom )
      DOUBLE PRECISION alpha, gam1, gam2, c_nom
* --------------------------------------------------
*        Normalization of electron distribution
* --------------------------------------------------
      DOUBLE PRECISION PI
      PARAMETER( PI = 3.1415927D0 )
      DOUBLE PRECISION al1
*
      al1 = 1.d0 - alpha
      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( alpha, gam1, gam2, c_nom )
      DOUBLE PRECISION alpha, gam1, gam2, c_nom
      INTEGER NGDIM
      PARAMETER( NGDIM = 128 )
      INTEGER ng
      DOUBLE PRECISION gi, gbm, gbp, disfc, delg
      COMMON /gint/ gi(NGDIM), gbm(NGDIM), gbp(NGDIM), disfc(NGDIM),
     &              delg, ng
      INTEGER i
      DOUBLE PRECISION exdel, temp
*
      ng = NGDIM
      gi(1) = gam1
      gi(ng) = gam2
      delg = log( gi(ng)/gi(1) ) / float( ng-1 )
      exdel = exp( delg )
      do 11 i = 2, ng-1
         gi(i) = gi(i-1) * exdel
 11   continue
      do 21 i = 1, ng
         temp = sqrt( gi(i)*gi(i) - 1.d0 )
         gbp(i) = gi(i) + temp
         gbm(i) = gi(i) - temp
 21   continue

      do 31 i = 1, ng
         disfc(i) = c_nom / gi(i)**alpha 
     &             / sqrt( gi(i)*gi(i) - 1.d0 )
 31   continue
      return
      END
*
*
      SUBROUTINE mesh( xlow, xhigh, dxl )
      DOUBLE PRECISION xlow, xhigh, dxl 
* ----------------------------------------------
*  Mesh in photon energy
*  Mesh in angle
* ----------------------------------------------
      INTEGER NX0, MU
      PARAMETER( NX0 = 1026, MU = 10 )
      INTEGER nx, mc
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION rmu, rmh, drmu
      COMMON /anglms/ rmu(-MU:2*MU+1), rmh(-MU:2*MU+1), drmu, mc
      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
*
*                                    ! mesh of mu = cos theta
      mc = MU
      rmu(-mc) = -1.d0
      rmu( mc) = 1.d0
      drmu = 1.0 / float( mc )
      do 50 j = -mc+1, mc-1
         rmu(j) = rmu(j-1) + drmu
 50   continue
      rmu(0) = 0.d0
*
      do 60 j = -mc, -1
         rmh(j) = 0.5d0*( rmu(j+1) + rmu(j) )
 60   continue
      do 61 j = 1, mc
         rmh(j) = 0.5d0*( rmu(j-1) + rmu(j) )
 61   continue
      rmh(0) = 0.d0
*
      return
      END
*
*
      SUBROUTINE initfx
* ----------------------------------------
*       Initialization of flux, etc.
* ----------------------------------------
      INTEGER NX0, MU
      PARAMETER( NX0 = 1026, MU = 10 )
      INTEGER nx, mc
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION rmu, rmh, drmu
      COMMON /anglms/ rmu(-MU:2*MU+1), rmh(-MU:2*MU+1), drmu, mc
      INTEGER noph
      DOUBLE PRECISION esbb, sinpt
      COMMON /penbin/ esbb(NX0,-MU:2*MU+1), sinpt(NX0), noph(NX0)
      DOUBLE PRECISION esj, fftrn, ffref, eout
      COMMON /epcptn/ esj(NX0,-MU:2*MU+1)
      COMMON /flaptn/ fftrn(NX0), ffref(NX0)
      COMMON /angdis/ eout(NX0,-MU:2*MU+1)
      INTEGER m, j
*      
      do 11 m = -mc, mc
         do 10 j = 1, nx
            esbb(j,m) = 0.d0
            esj(j,m)  = 0.d0
            eout(j,m) = 0.d0
 10      continue
 11   continue
      do 12 j = 1, nx
         fftrn(j) = 0.d0
         ffref(j) = 0.d0
         noph(j)  = 0
 12   continue
*
      return
      END
*
*
      SUBROUTINE sums
* ------------------------------------
*       Get flux, etc.
* ------------------------------------
      INTEGER NX0, MU
      PARAMETER( NX0 = 1026, MU = 10 )
      DOUBLE PRECISION PI
      PARAMETER( PI = 3.1415927D0 )
      INTEGER nx, mc
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION rmu, rmh, drmu
      COMMON /anglms/ rmu(-MU:2*MU+1), rmh(-MU:2*MU+1), drmu, mc
      DOUBLE PRECISION trn, ref
      COMMON /sums00/trn(NX0), ref(NX0)
      INTEGER noph
      DOUBLE PRECISION esbb, sinpt      
      COMMON /penbin/ esbb(NX0,-MU:2*MU+1), sinpt(NX0), noph(NX0)
      DOUBLE PRECISION esj, fftrn, ffref, eout
      COMMON /epcptn/ esj(NX0,-MU:2*MU+1)
      COMMON /flaptn/ fftrn(NX0), ffref(NX0)
      COMMON /angdis/ eout(NX0,-MU:2*MU+1)
      INTEGER j, m
      DOUBLE PRECISION cirmu
*
*                    ! fftrn(j): normalized transmitted energy flux
*                    ! ffref(j): normalized reflected energy flux
*
      cirmu = 2.d0*PI*drmu
      do 18 j = 1, nx
         trn(j) = 0.5d0*( esj(j,1)*rmh(1) + esj(j,mc)*rmh(mc) )
	 do 16 m = 2, mc-1
   	    trn(j) = trn(j) + esj(j,m)*rmh(m)
 16      continue
         trn(j) = cirmu*trn(j)
         ref(j) = 0.5d0*( esj(j,-mc)*abs(rmh(-mc))
     &	                + esj(j,-1)*abs(rmh(-1)) )
	 do 17 m = -mc+1, -2
   	    ref(j) = ref(j) + esj(j,m)*abs(rmh(m))
 17      continue
         ref(j) = cirmu*ref(j)
 18   continue
*
      do 22 m = -mc, mc
         do 21 j = 1, nx
	    eout(j,m) = eout(j,m) + esj(j,m)
 21      continue
 22   continue	    
*  
      do 50 j = 1, nx
         fftrn(j) = fftrn(j) + trn(j)
         ffref(j) = ffref(j) + ref(j)
 50   continue
*
*                                      ! initialize for next step
      do 68 m = -mc, mc
         do 65 j = 1, nx
            esj(j,m) = 0.d0
 65      continue
 68   continue
*
      return
      END
*
*
      SUBROUTINE normal( nphoton, ntotin, ftjwa, frjwa )
      INTEGER nphoton, ntotin
      DOUBLE PRECISION ftjwa, frjwa
* ------------------------------------------------------
*      normalization of flux, etc.
*      normalization:  sum f(i) = int f(x) dln(x) = 1
* ------------------------------------------------------
      INTEGER NX0, MU
      PARAMETER( NX0 = 1026, MU = 10 )
      DOUBLE PRECISION PI
      PARAMETER( PI = 3.1415927D0 )
      INTEGER nx, mc
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION rmu, rmh, drmu
      COMMON /anglms/ rmu(-MU:2*MU+1), rmh(-MU:2*MU+1), drmu, mc
      DOUBLE PRECISION  fftrn, ffref
      COMMON /flaptn/ fftrn(NX0), ffref(NX0)
      INTEGER noph
      DOUBLE PRECISION esbb, sinpt
      COMMON /penbin/ esbb(NX0,-MU:2*MU+1), sinpt(NX0), noph(NX0)
      INTEGER j, m
      DOUBLE PRECISION cirmu, was, fwa, flnp, devi
*
*                                                ! input photons
      cirmu = 2.d0*PI*drmu
      do 11 j = 1, nx
   	 sinpt(j) = 0.5d0*( esbb(j,1)*rmh(1) + esbb(j,mc)*rmh(mc) )
         do 10 m = 2, mc-1
   	    sinpt(j) = sinpt(j) + esbb(j,m)*rmh(m)
 10      continue
         sinpt(j) = cirmu*sinpt(j)
 11   continue
      was = 0.d0
      do 12 j = 1, nx
         was = was + sinpt(j)
 12   continue
      do 13 j = 1, nx
         sinpt(j) = sinpt(j)/was
 13   continue
*
*                                                ! output photons
      ftjwa = 0.d0
      frjwa = 0.d0
      do 25 j = 1, nx
         ftjwa = ftjwa + fftrn(j)
         frjwa = frjwa + ffref(j)
 25   continue
      fwa = ftjwa + frjwa
      ftjwa = ftjwa/fwa
      frjwa = frjwa/fwa
*
      flnp = float( nphoton )
      devi = 0.675d0/sqrt( flnp )
      do 31 j = 1, nx
         fftrn(j) = fftrn(j)/fwa
         ffref(j) = ffref(j)/fwa
 31   continue
*
      ntotin = 0
      do 41 j = 1, nx
         ntotin = ntotin + noph(j)
 41   continue
*
* ------------------------------------------------------------
      write(6,911) ntotin, ftjwa, frjwa
  911 format(1h ,'output photon number =',i8/
     &       1x, 'transmitted energy flux fraction =',1pe11.4/
     &       1x, 'reflected   energy flux fraction =',e11.4/)
* ------------------------------------------------------------
      return
      END
*
*
      SUBROUTINE phnxyz
* ------------------------------------------------
*  Initial values of photon position and direction
* ------------------------------------------------
      DOUBLE PRECISION PI, TPI
      PARAMETER( PI = 3.1415927D0, TPI = 2.D0*PI )
*
      REAL rvec1(1), rvec2(2)
      DOUBLE PRECISION xph, yph, zph, omg01, omg02, omg03
      COMMON /positn/xph, yph, zph, omg01, omg02, omg03
      DOUBLE PRECISION sqom3, tpxxi
*
*                               + ( x, y, z )
*                                  unit = thickness of a disk (h)
* Injection at the bottom
c      xph = 0.d0
c      yph = 0.d0
c      zph = 0.d0
*
* Uniform Injection
      call ranmar( rvec1, 1 )
      xph = 0.d0
      yph = 0.d0
      zph = rvec1(1)
*                               + propagation direction
      call ranmar( rvec2, 2 )
      omg03 = rvec2(1)
      sqom3 = sqrt( 1.d0 - omg03*omg03 )
      tpxxi = TPI * rvec2(2)
      omg02 = sqom3*sin( tpxxi )
      omg01 = sqom3*cos( tpxxi )

* Single direction (perpendicular injection)
*
c      omg01 = 0.d0
c      omg02 = 0.d0
c      omg03 = 1.d0
*
      return
      END
*
*
      SUBROUTINE scattr( alpha, gam1, gam2, taues, x0, epsct )
      DOUBLE PRECISION alpha, gam1, gam2, taues, x0, epsct
* -------------------------------------------------------------
*              Propagation and scattering
* -------------------------------------------------------------
      INTEGER NX0, MU 
      PARAMETER( NX0 = 1026, MU = 10 )
*
      DOUBLE PRECISION w(2)
      DOUBLE PRECISION xph, yph, zph, omg01, omg02, omg03
      COMMON /positn/ xph, yph, zph, omg01, omg02, omg03
      INTEGER nx, mc
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION rmu, rmh, drmu
      COMMON /anglms/ rmu(-MU:2*MU+1), rmh(-MU:2*MU+1), drmu, mc
      INTEGER noph
      DOUBLE PRECISION esbb, sinpt
      COMMON /penbin/ esbb(NX0,-MU:2*MU+1), sinpt(NX0), noph(NX0)
      DOUBLE PRECISION esj
      COMMON /epcptn/ esj(NX0,-MU:2*MU+1)
      INTEGER j, m
      REAL gxi
      DOUBLE PRECISION tauinv, hph, rlmean, pathli, freefr, fracli,
     &     escwti, rlamph, xob, omgp1, omgp2, omgp3, x0new
*
      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
            do 2 m = 1, mc
               if( omg03 .le. rmu(m) ) then
                  esbb(j,m) = esbb(j,m) + 1.d0
                  go to 6
               end if
 2          continue
         end if
 5    continue
 6    continue
*
 50   continue
        hph = x0 + x0
*                                         ! mean free path
        if( gam2*hph .le. 0.01d0 ) then
           rlmean = tauinv
         else
           call meanfr( taues, hph, rlmean )
        end if
*                                         ! distance to boundary
        if( omg03 .gt. 0.d0 ) then
            pathli = ( 1.d0 - zph ) / omg03
          else if( omg03 .lt. 0.d0 ) then
            pathli = -zph / omg03
        end if
        freefr = pathli / rlmean
        fracli = exp( -freefr )
        if( omg03 .eq. 0.d0 ) fracli = 0.d0
*
*                                         ! escaping photons
        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
              if( omg03 .lt. 0.d0 ) then
                 do 101 m = -mc, -1
                    if( omg03 .lt. rmu(m+1) ) then
                       esj(j,m) = esj(j,m) + escwti
                       go to 112
                    end if
 101             continue
              else
                 do 102 m = 1, mc
                    if( omg03 .le. rmu(m) ) then
                       esj(j,m) = esj(j,m) + escwti
                       go to 112
                    end if
 102             continue	     
              end if
           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 - gxi*( 1.d0 - fracli ) )
        xph = xph + rlamph*omg01
        yph = yph + rlamph*omg02
        zph = zph + rlamph*omg03
*                                                       ! scattering
        call slcpmn( alpha, gam1, gam2, hph, xob )
        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)
      INTEGER i
      REAL rvec(4)
      DOUBLE PRECISION z3xi 
*
 1    call ranmar( rvec, 4 )
*
      z3xi = 1.202d0 * rvec(1)
      if( z3xi .lt. 1.d0 ) then
         x0 = -trad * log( 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( rvec(2)*rvec(3)*rvec(4) )
*
      return
      END
*
*
      SUBROUTINE meanfr( taues, hph, rlmean )
      DOUBLE PRECISION taues, hph, rlmean
* ----------------------------------------
*   rlmean ... mean free path of a photon
* ----------------------------------------
      DOUBLE PRECISION PI, SIXI
      PARAMETER( PI = 3.1415927D0, SIXI = 1.D0/6.D0 )
      INTEGER NGDIM
      PARAMETER( NGDIM = 128 )
      INTEGER ng
      DOUBLE PRECISION gi, gbm, gbp, disfc, delg
      COMMON /gint/ gi(NGDIM), gbm(NGDIM), gbp(NGDIM), disfc(NGDIM),
     &              delg, ng
      INTEGER j
      DOUBLE PRECISION gghh, sphi, gx, gx2, phil, phiu, gxo
*
      gghh = 2.d0/(3.d0*PI)*hph*hph/taues
*
      sphi = 0.d0
      do 100 j = 1, ng
         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.5 ) 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
      rlmean = gghh / sphi
*
      return
      END
*
*
      SUBROUTINE slcpmn( alpha, gam1, gam2, hph, xob )
      DOUBLE PRECISION alpha, gam1, gam2, hph, xob
* ----------------------------------------------------
*   Selection of electron momentum
*   Isotropic, relativistic non-thermal distribution
* ----------------------------------------------------
      INTEGER NGDIM
      PARAMETER( NGDIM = 128 )
      DOUBLE PRECISION PI, THRI, TPI
      PARAMETER( PI = 3.1415927D0, THRI = 1.D0/3.D0, TPI = 2.D0*PI )
*
      DOUBLE PRECISION xph, yph, zph, omg01, omg02, omg03
      COMMON /positn/ xph, yph, zph, omg01, omg02, omg03
      DOUBLE PRECISION gmma, voc, rmu, etap, v01, v02, v03
      COMMON /neweng/ gmma, voc, rmu, etap, v01, v02, v03
      INTEGER ng
      DOUBLE PRECISION gi, gbm, gbp, disfc, delg
      COMMON /gint/ gi(NGDIM), gbm(NGDIM), gbp(NGDIM), disfc(NGDIM),
     &              delg, ng
      DOUBLE PRECISION sum_g
      COMMON /sum_el/ sum_g(NGDIM)
      INTEGER i
      REAL rvec(3), xi
      DOUBLE PRECISION al1, a_nom, tx2, v3s, rmv, xob2, sigmh
*
      al1 = 1.d0 - alpha
      if( abs(al1) .gt. 1.d-8) then
         a_nom = al1/(gam2**al1-gam1**al1)      
      else
         a_nom = 1.d0/log(gam2/gam1)
      end if
*                                                    ! direction
 1    call ranmar( rvec, 3 )
      v03 = rvec(1) + rvec(1) - 1.d0
      v3s = 2.d0*sqrt( rvec(1) )*sqrt( 1.d0 - rvec(1) )
      tx2 = TPI*rvec(2)
      v02 = v3s*sin( tx2 )
      v01 = v3s*cos( tx2 )
*
*                                      ! momentum
*                                             etap = p/mc
*                                             gmma = Lorentz factor
*
      call ranmar( xi, 1 )
      if( abs(al1) .gt. 1.d-8 ) then
         gmma = ( al1/a_nom*xi + gam1**al1 )**(1.d0/al1)
      else
         gmma = gam1*exp( xi/a_nom )
      end if
      etap = sqrt( gmma - 1.d0 ) * sqrt( gmma + 1.d0 )
*
      voc = sqrt(1.d0-1.d0/gmma) * sqrt(1.d0+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( rvec(3) .ge. 0.375d0 * sigmh * rmv ) go to 1
*
      do 101 i = 1, ng
         if( gmma .ge. gi(i) .and. gmma .lt. gi(i+1) ) then
            sum_g(i) = sum_g(i) + 1.d0
            go to 111
         end if
 101  continue
 111  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
      COMMON /positn/ xph, yph, zph, omg01, omg02, omg03
      DOUBLE PRECISION gmma, voc, rmu, etap, v01, v02, v03
      COMMON /neweng/ gmma, voc, rmu, etap, 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 = 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*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
       if(rmupv .lt. 1.d-15 ) then
          write(6,*) 'rmup =',rmup,'  voc =',voc
          write(6,*) 'gmma =',gmma
       end if
      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( 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
*
      sm(1) = 1.d0
      do 10 j = 2, MDIM
        sm(j) = sm(j-1) + 1.d0 / float( j*j*j )
 10   continue
*
      return
      END
*
*
      SUBROUTINE meanen( trad, x0mean, x1mean, x1trnm, x1refm, acomp )
      DOUBLE PRECISION trad, x0mean, x1mean, x1trnm, x1refm, acomp
* --------------------------------------------------------------------
*   Mean photon energies:
*
*   (average wighted by number flux
*    <x> = int F(x) dx / int F(x)/x dx = sum F(i) x(i)
*    normalization  sum F(i) = 1 = int F(x)/x dx = int F(x) dln(x)
*    note: integration is in logarithm)
*
*       x0mean ... input photons
*       x1mean ... output photons
*       x1trnm ... transmitted photons
*       x1refm ... refrected photons
*       acomp ... luminosity enhancement factor
*                 L = acomp * L_0
* --------------------------------------------------------------------
      INTEGER NX0, MU
      PARAMETER( NX0 = 1026, MU = 10 )
      DOUBLE PRECISION x(NX0)
      INTEGER nx
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION fftrn, ffref
      COMMON /flaptn/ fftrn(NX0), ffref(NX0)
      INTEGER noph
      DOUBLE PRECISION esbb, sinpt
      COMMON /penbin/ esbb(NX0,-MU:2*MU+1), sinpt(NX0), noph(NX0)
      INTEGER j
      DOUBLE PRECISION hh, winp, winpn, wtrn, wtrnn, wref, wrefn,
     &     ftr, frf, waf, acompm1, 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.d0
      wtrnn = 0.d0
      wref = 0.d0
      wrefn = 0.d0
      do 22 j = 1, nx
         wtrn = wtrn + fftrn(j)*x(j)
         wtrnn = wtrnn + fftrn(j)
         wref = wref + ffref(j)*x(j)
	 wrefn = wrefn + ffref(j)
 22   continue
      x1trnm = wtrn/wtrnn
      x1refm = wref/wrefn
*
      ftr = wtrnn/winpn
      frf = wrefn/winpn
      waf = ftr + frf
      ftr = ftr/waf
      frf = frf/waf
*         
      x1mean = x1trnm*ftr + x1refm*frf
      acomp = x1mean/x0mean
      acompm1 = (x1mean - x0mean)/x0mean
      x0thry = 2.701178d0*trad
*
* ------------------------------------------------------------------
      write(6,910) x0thry, x0mean, x1mean, x1trnm, x1refm, acomp,
     &     acompm1
 910  format(1h ,'Mean of initial photon energy (theory) =',1pe11.4/
     &   1x,'(average by number flux)'/
     &   1x,'++ Rimulation Results ++'/
     &   1x ,'Mean of initial     photon energy =', 1pe11.4/
     &   1x, 'Mean of emission    photon energy =', e11.4/
     &   1x, 'Mean of transmitted photon energy =', e11.4/
     &   1x, 'Mean of reflected   photon energy =', e11.4/
     &   1x, 'Luminosity enhancement factor a =',e11.4/
     &   1x, 'a - 1 =',e11.4/
     &   1x, '(Comptonized Luminosity L = a L_0)'/)
* ------------------------------------------------------------------
*
      return
      END
*
*
      SUBROUTINE angular
* --------------------------------------------------------
*  Angular distribution and Reflection/Transmission ratio
*   rt1, etc ... Energy flux ratio F-/F+
*                F- ... reflection, F+ ... transmission
*   arf1(m) ... Angular distribution of reflected intensity * cos theta
*               I * mu
*   atr1(m) ... Angular distribution of transmietted intensity
*               I * mu
* in energy range (xa1, xa2), etc.
* --------------------------------------------------------
      INTEGER NX0, MU
      PARAMETER( NX0 = 1026, MU = 10 )
      DOUBLE PRECISION XA1, XA2, XB1, XB2, XC1, XC2, XD1, XD2,
     &                 XE1, XE2
      PARAMETER( XA1 = 0.2D0/511.D0,  XA2 = 0.5D0/511.D0,
     &           XB1 = 1.D0/511.D0,   XB2 = 3.D0/511.D0,
     &           XC1 = 5.D0/511.D0,   XC2 = 15.D0/511.D0,
     &           XD1 = 25.D0/511.D0,  XD2 = 75.D0/511.D0,
     &           XE1 = 100.D0/511.D0, XE2 = 1000.D0/511.D0 )
      DOUBLE PRECISION x(NX0)
      INTEGER nx, mc
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION rmu, rmh, drmu
      COMMON /anglms/ rmu(-MU:2*MU+1), rmh(-MU:2*MU+1), drmu, mc
      DOUBLE PRECISION eout, fftrn, ffref
      COMMON /angdis/ eout(NX0,-MU:2*MU+1)
      COMMON /flaptn/ fftrn(NX0), ffref(NX0)
      DOUBLE PRECISION rt1, rt2, rt3, rt4, rt5, arf1, arf2, arf3,
     &     arf4, arf5, atr1, atr2, atr3, atr4, atr5
      COMMON /tempem/ rt1, rt2, rt3, rt4, rt5,
     &   arf1(MU), atr1(MU), arf2(MU), atr2(MU), arf3(MU), atr3(MU),
     &   arf4(MU), atr4(MU), arf5(MU), atr5(MU)
*
      INTEGER j, k, m
      DOUBLE PRECISION fp1, fp2, fp3, fp4, fp5, fm1, fm2, fm3, fm4, fm5
      DOUBLE PRECISION hh, s1r, s2r, s3r, s4r, s5r, s1t, s2t, s3t, s4t,
     &                 s5t
*
      DATA fp1, fp2, fp3, fp4, fp5, fm1, fm2, fm3, fm4, fm5/
     &      0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0./
      DATA arf1/MU*0./, atr1/MU*0./, arf2/MU*0./, atr2/MU*0./,
     &     arf3/MU*0./, atr3/MU*0./, arf4/MU*0./, atr4/MU*0./,
     &     arf5/MU*0./, atr5/MU*0./
*
      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)
*
      do 111 j = 1, nx-1
         if( xhn(j) .ge. XA1 .and. xhn(j+1) .le. XA2 ) then
            fp1 = fp1 + fftrn(j)
            fm1 = fm1 + ffref(j)
            do 121 m = -1, -mc, -1
               k = -m
               arf1(k) = arf1(k) + eout(j,m)*abs( rmh(m) )
 121        continue
            do 122 m = 1, mc
               atr1(m) = atr1(m) + eout(j,m)*rmh(m)
 122        continue
         else if( xhn(j) .ge. XB1 .and. xhn(j+1) .le. XB2 ) then
            fp2 = fp2 + fftrn(j)
            fm2 = fm2 + ffref(j)
            do 123 m = -1, -mc, -1
               k = -m
               arf2(k) = arf2(k) + eout(j,m)*abs( rmh(m) )
 123        continue
            do 124 m = 1, mc
               atr2(m) = atr2(m) + eout(j,m)*rmh(m)
 124        continue
         else if( xhn(j) .ge. XC1 .and. xhn(j+1) .le. XC2 ) then
            fp3 = fp3 + fftrn(j)
            fm3 = fm3 + ffref(j)
            do 125 m = -1, -mc, -1
               k = -m
               arf3(k) = arf3(k) + eout(j,m)*abs( rmh(m) )
 125        continue
            do 126 m = 1, mc
               atr3(m) = atr3(m) + eout(j,m)*rmh(m)
 126        continue
         else if( xhn(j) .ge. XD1 .and. xhn(j+1) .le. XD2 ) then
            fp4 = fp4 + fftrn(j)
            fm4 = fm4 + ffref(j)
            do 127 m = -1, -mc, -1
               k = -m
               arf4(k) = arf4(k) + eout(j,m)*abs( rmh(m) )
 127        continue
            do 128 m = 1, mc
               atr4(m) = atr4(m) + eout(j,m)*rmh(m)
 128        continue
         else if( xhn(j) .ge. XE1 .and. xhn(j+1) .le. XE2 ) then
            fp5 = fp5 + fftrn(j)
            fm5 = fm5 + ffref(j)
            do 129 m = -1, -mc, -1
               k = -m
               arf5(k) = arf5(k) + eout(j,m)*abs( rmh(m) )
 129        continue
            do 130 m = 1, mc
               atr5(m) = atr5(m) + eout(j,m)*rmh(m)
 130        continue
         end if
 111  continue	 
      rt1 = fm1/fp1
      rt2 = fm2/fp2
      rt3 = fm3/fp3
      rt4 = fm4/fp4
      rt5 = fm5/fp5
*
      s1r = 0.5d0*( arf1(1) + arf1(mc) )
      s1t = 0.5d0*( atr1(1) + atr1(mc) )
      s2r = 0.5d0*( arf2(1) + arf2(mc) )
      s2t = 0.5d0*( atr2(1) + atr2(mc) )
      s3r = 0.5d0*( arf3(1) + arf3(mc) )
      s3t = 0.5d0*( atr3(1) + atr3(mc) )
      s4r = 0.5d0*( arf4(1) + arf4(mc) )
      s4t = 0.5d0*( atr4(1) + atr4(mc) )
      s5r = 0.5d0*( arf5(1) + arf5(mc) )
      s5t = 0.5d0*( atr5(1) + atr5(mc) )
      do 200 m = 2, mc-1
         s1r = s1r + arf1(m)
         s1t = s1t + atr1(m)
         s2r = s2r + arf2(m)
         s2t = s2t + atr2(m)
         s3r = s3r + arf3(m)
         s3t = s3t + atr3(m)
         s4r = s4r + arf4(m)
         s4t = s4t + atr4(m)
         s5r = s5r + arf5(m)
         s5t = s5t + atr5(m)
 200  continue
      do 201 m = 1, mc
         arf1(m) = arf1(m) / s1r / drmu
         atr1(m) = atr1(m) / s1t / drmu
         arf2(m) = arf2(m) / s2r / drmu
         atr2(m) = atr2(m) / s2t / drmu
         arf3(m) = arf3(m) / s3r / drmu
         atr3(m) = atr3(m) / s3t / drmu
         arf4(m) = arf4(m) / s4r / drmu 
         atr4(m) = atr4(m) / s4t / drmu
         arf5(m) = arf5(m) / s5r / drmu
         atr5(m) = atr5(m) / s5t / drmu
 201  continue
*
      write(6,910) rt1, rt2, rt3, rt4, rt5
 910  format(1h ,'f- : reflected energy flux'/
     &       1x, 'f+ : transmitted energy flux'/
     &       1x, '0.2 - 0.5  kev  f-/f+ =',1pe11.4/
     &       1x, '1 - 3      kev  f-/f+ =',e11.4/
     &       1x, '5 - 15     kev  f-/f+ =',e11.4/
     &       1x, '25 - 75    kev  f-/f+ =',e11.4/
     &       1x, '100 - 1000 kev  f-/f+ =',e11.4)
      write(13,920) ( rmh(m), arf1(m), atr1(m), arf2(m), atr2(m),
     &  arf3(m), atr3(m), m = 1, mc )
 920  format(1h ,7(1pe11.4))
      write(14,930) ( rmh(m), arf4(m), atr4(m), arf5(m), atr5(m),
     &   m = 1, mc )
 930  format(1h ,5(1pe11.4))
*
      return
      END
*
*
      SUBROUTINE sp_index( e_min, e_max, power_t, power_r, 
     &                     sigpi_t, sigpi_r )
      DOUBLE PRECISION e_min, e_max
      REAL power_t, power_r, sigpi_t, sigpi_r
* ------------------------------------------------------
*   Power-law fit of emission spectrum
*   in the energy rage (E_MIN, E_MAX)
*     input ... e_min and e_max in keV
*     output ... power_t: power-law index for transmitted photons
*                sigpi_t: probable uncertaity
*                power_r: power-law index for reflected photons
*                sigpi_r: probable uncertaity
* ------------------------------------------------------
      INTEGER NX0
      PARAMETER( NX0 = 1026 )
      REAL sig(NX0), xin(NX0), trn(NX0), ref(NX0)
      INTEGER nx
      DOUBLE PRECISION xhn
      COMMON /enrbin/ xhn(0:NX0), nx
      DOUBLE PRECISION fftrn, ffref
      COMMON /flaptn/ fftrn(NX0), ffref(NX0)
      INTEGER i, ndata, indatr, indarf, mwt
      REAL acons, siga, chi2, q
      DOUBLE PRECISION x_min, x_max
*
      x_min = e_min / 510.99906d0
      x_max = e_max / 510.99906d0
*                                           ! transmitted spectrum
      indatr = 99
      ndata = 0
      do 11 i = 1, nx
         if( xhn(i) .ge. x_min .and. xhn(i) .lt. x_max ) then
            if( fftrn(i) .le. 0.d0 ) then
               indatr = 0
               go to 15
            end if
            ndata = ndata + 1
            xin(ndata) = sngl( log10( xhn(i) ) )
            trn(ndata) = sngl( log10( fftrn(i) ) )
	 end if    
 11   continue
      mwt = 0
      call fit( xin, trn, ndata, sig, mwt, acons, power_t, 
     &          siga, sigpi_t, chi2, q)
 15   continue
*
*                                           ! reflected spectrum
      indarf = 99
      ndata = 0
      do 21 i = 1, nx
         if( xhn(i) .ge. x_min .and. xhn(i) .lt. x_max ) then
            if( ffref(i) .le. 0.d0 ) then
               indarf = 0
               go to 25
            end if
            ndata = ndata + 1
            xin(ndata) = sngl( log10( xhn(i) ) )
            ref(ndata) = sngl( log10( ffref(i) ) )
	 end if    
 21   continue
      mwt = 0
      call fit( xin, ref, ndata, sig, mwt, acons, power_r, 
     &         siga, sigpi_r, chi2, q)     
 25   continue
*
      write(6,91)  e_min, e_max
      if( indatr .gt. 0 ) then
         write(6,92)  power_t, sigpi_t
       else
         write(6,103)
      end if
      if( indarf .gt. 0 ) then
         write(6,93)  power_r, sigpi_r
       else
         write(6,105)
      end if
 91   format(1h ,/1x,
     &   'power-law index in the energy range:',1pe10.3,1x,'-',
     &   e10.3,1x,'keV')
 92   format(1h ,'transmitted photons: index =',1pe11.4/
     &        1x,'       probable uncertaity =',e11.4)
 93   format(1h ,'reflected photons:   index =',1pe11.4/
     &        1x,'       probable uncertaity =',e11.4)
 103  format(1h ,'transmitted spectrum ... N/A')
 105  format(1h ,'reflected   spectrum ... N/A')
      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, mod
      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
*
*
      SUBROUTINE fit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
      INTEGER mwt,ndata
      REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
CU    USES gammq
      INTEGER i
      REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,gammq
      sx=0.
      sy=0.
      st2=0.
      b=0.
      if(mwt.ne.0) then
        ss=0.
        do 11 i=1,ndata
          wt=1./(sig(i)**2)
          ss=ss+wt
          sx=sx+x(i)*wt
          sy=sy+y(i)*wt
11      continue
      else
        do 12 i=1,ndata
          sx=sx+x(i)
          sy=sy+y(i)
12      continue
        ss=float(ndata)
      endif
      sxoss=sx/ss
      if(mwt.ne.0) then
        do 13 i=1,ndata
          t=(x(i)-sxoss)/sig(i)
          st2=st2+t*t
          b=b+t*y(i)/sig(i)
13      continue
      else
        do 14 i=1,ndata
          t=x(i)-sxoss
          st2=st2+t*t
          b=b+t*y(i)
14      continue
      endif
      b=b/st2
      a=(sy-sx*b)/ss
      siga=sqrt((1.+sx*sx/(ss*st2))/ss)
      sigb=sqrt(1./st2)
      chi2=0.
      if(mwt.eq.0) then
        do 15 i=1,ndata
          chi2=chi2+(y(i)-a-b*x(i))**2
15      continue
        q=1.
        sigdat=sqrt(chi2/(ndata-2))
        siga=siga*sigdat
        sigb=sigb*sigdat
      else
        do 16 i=1,ndata
          chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
16      continue
        q=gammq(0.5*(ndata-2),0.5*chi2)
      endif
      return
      END
      FUNCTION gammq(a,x)
      REAL a,gammq,x
CU    USES gcf,gser
      REAL gammcf,gamser,gln
      if(x.lt.0..or.a.le.0.)pause 'bad arguments in gammq'
      if(x.lt.a+1.)then
        call gser(gamser,a,x,gln)
        gammq=1.-gamser
      else
        call gcf(gammcf,a,x,gln)
        gammq=gammcf
      endif
      return
      END
      SUBROUTINE gcf(gammcf,a,x,gln)
      INTEGER ITMAX
      REAL a,gammcf,gln,x,EPS,FPMIN
      PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
CU    USES gammln
      INTEGER i
      REAL an,b,c,d,del,h,gammln
      gln=gammln(a)
      b=x+1.-a
      c=1./FPMIN
      d=1./b
      h=d
      do 11 i=1,ITMAX
        an=-i*(i-a)
        b=b+2.
        d=an*d+b
        if(abs(d).lt.FPMIN)d=FPMIN
        c=b+an/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        del=d*c
        h=h*del
        if(abs(del-1.).lt.EPS)goto 1
11    continue
      pause 'a too large, ITMAX too small in gcf'
1     gammcf=exp(-x+a*log(x)-gln)*h
      return
      END
      SUBROUTINE gser(gamser,a,x,gln)
      INTEGER ITMAX
      REAL a,gamser,gln,x,EPS
      PARAMETER (ITMAX=100,EPS=3.e-7)
CU    USES gammln
      INTEGER n
      REAL ap,del,sum,gammln
      gln=gammln(a)
      if(x.le.0.)then
        if(x.lt.0.)pause 'x < 0 in gser'
        gamser=0.
        return
      endif
      ap=a
      sum=1./a
      del=sum
      do 11 n=1,ITMAX
        ap=ap+1.
        del=del*x/ap
        sum=sum+del
        if(abs(del).lt.abs(sum)*EPS)goto 1
11    continue
      pause 'a too large, ITMAX too small in gser'
1     gamser=sum*exp(-x+a*log(x)-gln)
      return
      END
      FUNCTION gammln(xx)
      REAL gammln,xx
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     *-.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      END