From z96020@kgupyr.kwansei.ac.jp Thu Jun 6 02:12:51 1996 Received: from mailgate.kwansei.ac.jp by spacsun.rice.edu (CAA03939); Thu, 6 Jun 1996 02:11:46 -0500 (CDT) From: z96020@kgupyr.kwansei.ac.jp Received: from kgupyr.kwansei.ac.jp ([192.218.161.23]) by mailgate.kwansei.ac.jp (8.6.12+2.4W/3.3W9-95090519) with ESMTP id QAA27052 for ; Thu, 6 Jun 1996 16:05:47 +0900 Received: (from z96020@localhost) by kgupyr.kwansei.ac.jp (8.6.12+2.4W/3.4W-95092211) id QAA29449 for liang@spacsun.rice.edu; Thu, 6 Jun 1996 16:09:21 +0900 Date: Thu, 6 Jun 1996 16:09:21 +0900 Message-Id: <199606060709.QAA29449@kgupyr.kwansei.ac.jp> To: liang@spacsun.rice.edu Subject: disk code Status: RO * * -------------------------------------------- * 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 * = 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