* -------------------------------------------- * MONTE CARLO SIMULATION OF COMPTON SCATTERING * -------------------------------------------- * * @ GEOMETRY ... sphere * thickness of a cloud ... taues * * @ ELECTRONS ... isotropic, non-thermal * @ INPUT PHOTONS ... planckian * * + INPUT : * * alpha: power-law index of nonthermal electrons * 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 * in the center of a sphere * 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 * * cool ... compton cooling rate per unit volume (erg/s/cm^3) * This is calculated, when the physical radicu R(cm) is given * * * SYMBOLS : * xhn(j) ... photon energy in units of electron mass * ( dimension = nx0 ) * w(i) ... weight * -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- * -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- * 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, 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 ... Boltzmann 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, aff COMMON /statis/ snj(NX0), snjs(NX0), dsnj(NX0), esbj(NX0) COMMON /absvle/ aff(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, ntotin, nphoton DOUBLE PRECISION alpha, gam1, gam2, tradkv, taues, & epsct, tradk, trad, xlow, xhigh, dxl, c_nom, x0, & snjwa, x0mean, x1mean, acomp, elnorm DOUBLE PRECISION flinpt, cool, radius * * * ! Input parameters * c alpha = 2.d0 c gam1 = 30.d0 c gam2 = 1.d5 c tradkv = 511.d0*1.d-8 c taues = 1.d-5 c nphoton = 10000 c epsct = 1.d-12 c open(8,file='inp304') 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-10 xhigh = 1.d6 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 * ! initialize random number generator * the input value should be in the range: 0 < = ijkl < = 900 000 000 * ijkl = 199615 call rmarin(ijkl) * * -------------- BEGINNING OF SCATTERING --------------- * call sumofm call initfx do 51 i = 1, nphoton if( mod(i,1000) .eq. 0 ) then open(25,file='number') write(25,*) 'number =',i close(25) 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, snjwa ) * * ! mean photon energies call meanen( trad, x0mean, x1mean, acomp ) * * --------------------------------------------- * open(12,file='spect-sph.dat') 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) * 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-sph.dat') 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 ... Spherical 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/) 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( 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 initfx * ---------------------------------------- * Initialization of flux, etc. * ---------------------------------------- INTEGER NX0, MU PARAMETER( NX0 = 1026, MU = 10 ) 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 absflx( acomp, flinpt ) DOUBLE PRECISION acomp, flinpt * ------------------------------------------ * Absolute value of emission intensity aff * aff = a_c * b * a_c ... Compton amplification factor * b ... blackbody * ------------------------------------------ INTEGER NX0 PARAMETER( NX0 = 1026 ) DOUBLE PRECISION PI PARAMETER( PI = 3.1415927d0 ) INTEGER nx DOUBLE PRECISION xhn, aff, snj, snjs, dsnj, esbj COMMON /enrbin/ xhn(0:NX0), nx COMMON /absvle/ aff(NX0) COMMON /statis/ snj(NX0), snjs(NX0), dsnj(NX0), esbj(NX0) INTEGER j DOUBLE PRECISION dlx, wa, fbtrd * dlx = log( xhn(2)/xhn(1) ) wa = 0.d0 do 5 j = 1, nx wa = wa + snj(j)*xhn(j) 5 continue fbtrd = acomp * flinpt / dlx / wa do 11 j = 1, nx aff(j) = fbtrd*snj(j) 11 continue * return END * * SUBROUTINE phnxyz * ------------------------------------------------ * initial values of photon position and direction * + input at the bottom of a sphere * + outward 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 * + ( x, y, z ) * unit = radius of a sphere xph = 0.d0 yph = 0.d0 zph = 0.d0 * + propagation direction omg03 = 1.d0 omg02 = 0.d0 omg01 = 0.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 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 REAL gxi INTEGER j DOUBLE PRECISION tauinv, hph, rlmean, rovecs, rleng2, & freefr, fracli, rlamph, x0new, xob, pathli, & omgp1, omgp2, omgp3 * 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 if( gam2*hph .le. 0.01d0 ) then rlmean = tauinv else call meanfr( 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 - 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) REAL rvec(4) INTEGER i 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 DOUBLE PRECISION gmma, voc, rmu, etap, v01, v02, v03 COMMON /positn/ xph, yph, zph, omg01, omg02, omg03 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 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 = xp / 2.d0 / 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 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