* -------------------------------------------- * 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.