Page MenuHomeHEPForge
Files F11221930

sample.f
No OneTemporary
Actions

sample.f
View Options

double precision function sigma (s)
implicit none
double precision s
sigma = 1d0 / s
end
double precision function d12 (t1, t2)
implicit none
double precision t1, t2, x1, x2, sigma, circee
double precision EPS, PWR
parameter (EPS = 1d-6, PWR = 5d0)
double precision KIREPS
parameter (KIREPS = 1D-6)
x1 = 1d0 - t1**PWR
x2 = 1d0 - t2**PWR
d12 = PWR*PWR * (t1*t2)**(PWR-1d0)
$ * sigma (x1*x2) * circee (x1, x2)
end
double precision function d1 (t1)
implicit none
double precision t1, x1, sigma, circee
double precision EPS, PWR
parameter (EPS = 1d-6, PWR = 5d0)
double precision KIREPS
parameter (KIREPS = 1D-6)
x1 = 1d0 - t1**PWR
d1 = PWR * t1**(PWR-1d0) * sigma (x1) * circee (x1, 1d0)
end
double precision function d2 (t2)
implicit none
double precision t2, x2, sigma, circee
double precision EPS, PWR
parameter (EPS = 1d-6, PWR = 5d0)
double precision KIREPS
parameter (KIREPS = 1D-6)
x2 = 1d0 - t2**PWR
d2 = PWR * t2**(PWR-1d0) * sigma (x2) * circee (1d0, x2)
end
double precision function d12a (x1, x2)
implicit none
double precision x1, x2, sigma, kirkee
d12a = sigma (x1*x2) * kirkee (x1, x2)
end
program sample
implicit none
integer SBAND, TESLA, XBAND
parameter (SBAND = 1, TESLA = 2, XBAND = 3)
integer JLCNLC
parameter (JLCNLC = 3)
integer SBNDEE, TESLEE, XBNDEE
parameter (SBNDEE = 4, TESLEE = 5, XBNDEE = 6)
integer NACC
parameter (NACC = 6)
double precision EPS, PWR
parameter (EPS = 1d-6, PWR = 5d0)
double precision KIREPS
parameter (KIREPS = 1D-6)
double precision s
double precision gauss1, gauss2, circee, sigma, d1, d2, d12, d12a
external d1, d2, d12, d12a
double precision w, s2, x1, x2
external random
integer NEVENT, n
parameter (NEVENT = 10000)
integer acc, ver, i
double precision roots(6)
data roots / 90D0, 170D0, 350D0, 500D0, 800D0, 1000D0 /
do 10 acc = TESLA, JLCNLC
do 11 ver = 7, 7
do 12 i = 1, 6
call circes (0d0, 0d0, roots(i), acc, ver, 20000501, 1)
s = sigma (1d0) * circee (1d0, 1d0)
$ + gauss1 (d1, 0d0, 1d0, EPS)
$ + gauss1 (d2, 0d0, 1d0, EPS)
$ + gauss2 (d12, 0d0, 1d0, 0d0, 1d0, EPS)
write (*, 1000) 'delta(sigma) (Gauss) =', (s-1d0)*100d0
1000 format (1X, A22, 1X, F6.2, '%')
s = gauss2 (d12a, 0d0, 1d0-KIREPS, 0d0, 1d0-KIREPS, EPS)
$ + gauss2 (d12a, 0d0, 1d0-KIREPS, 1d0-KIREPS, 1d0, EPS)
$ + gauss2 (d12a, 1d0-KIREPS, 1d0, 0d0, 1d0-KIREPS, EPS)
$ + gauss2 (d12a, 1d0-KIREPS, 1d0, 1d0-KIREPS, 1d0, EPS)
write (*, 1000) 'delta(sigma) (Gauss) =', (s-1d0)*100d0
s = 0d0
s2 = 0d0
do 100 n = 1, NEVENT
call gircee (x1, x2, random)
w = sigma (x1*x2)
s = s + w
s2 = s2 + w*w
100 continue
s = s / dble(NEVENT)
s2 = s2 / dble(NEVENT)
write (*, 1000) 'delta(sigma) (MC) =', (s-1d0)*100d0
write (*, 1000) ' +/-',
$ sqrt((s2-s*s)/dble(NEVENT))*100d0
14 continue
12 continue
13 continue
11 continue
10 continue
end
subroutine random (r)
implicit none
double precision r
integer m, a, c
parameter (M = 259200, A = 7141, C = 54773)
integer n
save n
data n /0/
n = mod(n*a+c,m)
r = dble (n) / dble (m)
end
double precision function gauss1 (f, a, b, eps)
implicit none
double precision f, a, b, eps
external f
double precision Z1, HF, CST
parameter (Z1 = 1, HF = Z1/2, CST = 5*Z1/1000)
integer i
double precision h, const, aa, bb, c1, c2, s8, s16, u
double precision w(12), x(12)
data x( 1) /9.6028985649753623d-1/, w( 1) /1.0122853629037626d-1/
data x( 2) /7.9666647741362674d-1/, w( 2) /2.2238103445337447d-1/
data x( 3) /5.2553240991632899d-1/, w( 3) /3.1370664587788729d-1/
data x( 4) /1.8343464249564980d-1/, w( 4) /3.6268378337836198d-1/
data x( 5) /9.8940093499164993d-1/, w( 5) /2.7152459411754095d-2/
data x( 6) /9.4457502307323258d-1/, w( 6) /6.2253523938647893d-2/
data x( 7) /8.6563120238783174d-1/, w( 7) /9.5158511682492785d-2/
data x( 8) /7.5540440835500303d-1/, w( 8) /1.2462897125553387d-1/
data x( 9) /6.1787624440264375d-1/, w( 9) /1.4959598881657673d-1/
data x(10) /4.5801677765722739d-1/, w(10) /1.6915651939500254d-1/
data x(11) /2.8160355077925891d-1/, w(11) /1.8260341504492359d-1/
data x(12) /9.5012509837637440d-2/, w(12) /1.8945061045506850d-1/
h = 0
if (b .eq. a) go to 99
const = CST / dabs(b-a)
bb = a
1 continue
aa = bb
bb = b
2 continue
c1 = HF*(bb+aa)
c2 = HF*(bb-aa)
s8 = 0
do 3 i = 1, 4
u = c2*x(i)
s8 = s8 + w(i) * (f (c1+u) + f (c1-u))
3 continue
s16 = 0
do 4 i = 5, 12
u = c2*x(i)
s16 = s16 + w(i) * (f (c1+u) + f (c1-u))
4 continue
s16 = c2*s16
if (dabs(s16-c2*s8) .le. eps*(1+dabs(s16))) then
h = h + s16
if (bb .ne. b) go to 1
else
bb = c1
if (1 + const*dabs(c2) .ne. 1) go to 2
h = 0
print *, 'gauss: too high accuracy required'
go to 99
end if
99 continue
gauss1 = h
end
double precision function gaussx (f, y, a, b, eps)
implicit none
double precision y
double precision f, a, b, eps
external f
double precision Z1, HF, CST
parameter (Z1 = 1, HF = Z1/2, CST = 5*Z1/1000)
integer i
double precision h, const, aa, bb, c1, c2, s8, s16, u
double precision w(12), x(12)
data x( 1) /9.6028985649753623d-1/, w( 1) /1.0122853629037626d-1/
data x( 2) /7.9666647741362674d-1/, w( 2) /2.2238103445337447d-1/
data x( 3) /5.2553240991632899d-1/, w( 3) /3.1370664587788729d-1/
data x( 4) /1.8343464249564980d-1/, w( 4) /3.6268378337836198d-1/
data x( 5) /9.8940093499164993d-1/, w( 5) /2.7152459411754095d-2/
data x( 6) /9.4457502307323258d-1/, w( 6) /6.2253523938647893d-2/
data x( 7) /8.6563120238783174d-1/, w( 7) /9.5158511682492785d-2/
data x( 8) /7.5540440835500303d-1/, w( 8) /1.2462897125553387d-1/
data x( 9) /6.1787624440264375d-1/, w( 9) /1.4959598881657673d-1/
data x(10) /4.5801677765722739d-1/, w(10) /1.6915651939500254d-1/
data x(11) /2.8160355077925891d-1/, w(11) /1.8260341504492359d-1/
data x(12) /9.5012509837637440d-2/, w(12) /1.8945061045506850d-1/
h = 0
if (b .eq. a) go to 99
const = CST / dabs(b-a)
bb = a
1 continue
aa = bb
bb = b
2 continue
c1 = HF*(bb+aa)
c2 = HF*(bb-aa)
s8 = 0
do 3 i = 1, 4
u = c2*x(i)
s8 = s8 + w(i) * (f (y, c1+u) + f (y, c1-u))
3 continue
s16 = 0
do 4 i = 5, 12
u = c2*x(i)
s16 = s16 + w(i) * (f (y, c1+u) + f (y, c1-u))
4 continue
s16 = c2*s16
if (dabs(s16-c2*s8) .le. eps*(1+dabs(s16))) then
h = h + s16
if (bb .ne. b) go to 1
else
bb = c1
if (1 + const*dabs(c2) .ne. 1) go to 2
h = 0
print *, 'gauss: too high accuracy required'
go to 99
end if
99 continue
gaussx = h
end
double precision function gauss2 (f, a, b, a1, b1, eps)
implicit none
double precision a1, b1, gaussx
double precision f, a, b, eps
external f
double precision Z1, HF, CST
parameter (Z1 = 1, HF = Z1/2, CST = 5*Z1/1000)
integer i
double precision h, const, aa, bb, c1, c2, s8, s16, u
double precision w(12), x(12)
data x( 1) /9.6028985649753623d-1/, w( 1) /1.0122853629037626d-1/
data x( 2) /7.9666647741362674d-1/, w( 2) /2.2238103445337447d-1/
data x( 3) /5.2553240991632899d-1/, w( 3) /3.1370664587788729d-1/
data x( 4) /1.8343464249564980d-1/, w( 4) /3.6268378337836198d-1/
data x( 5) /9.8940093499164993d-1/, w( 5) /2.7152459411754095d-2/
data x( 6) /9.4457502307323258d-1/, w( 6) /6.2253523938647893d-2/
data x( 7) /8.6563120238783174d-1/, w( 7) /9.5158511682492785d-2/
data x( 8) /7.5540440835500303d-1/, w( 8) /1.2462897125553387d-1/
data x( 9) /6.1787624440264375d-1/, w( 9) /1.4959598881657673d-1/
data x(10) /4.5801677765722739d-1/, w(10) /1.6915651939500254d-1/
data x(11) /2.8160355077925891d-1/, w(11) /1.8260341504492359d-1/
data x(12) /9.5012509837637440d-2/, w(12) /1.8945061045506850d-1/
h = 0
if (b .eq. a) go to 99
const = CST / dabs(b-a)
bb = a
1 continue
aa = bb
bb = b
2 continue
c1 = HF*(bb+aa)
c2 = HF*(bb-aa)
s8 = 0
do 3 i = 1, 4
u = c2*x(i)
s8 = s8 + w(i) * (gaussx (f, c1+u, a1, b1, eps)
$ + gaussx (f, c1-u, a1, b1, eps))
3 continue
s16 = 0
do 4 i = 5, 12
u = c2*x(i)
s16 = s16 + w(i) * (gaussx (f, c1+u, a1, b1, eps)
$ + gaussx (f, c1-u, a1, b1, eps))
4 continue
s16 = c2*s16
if (dabs(s16-c2*s8) .le. eps*(1+dabs(s16))) then
h = h + s16
if (bb .ne. b) go to 1
else
bb = c1
if (1 + const*dabs(c2) .ne. 1) go to 2
h = 0
print *, 'gauss: too high accuracy required'
go to 99
end if
99 continue
gauss2 = h
end

File Metadata

Mime Type
text/plain
Expires
Wed, May 14, 11:08 AM (21 h, 4 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111354
Default Alt Text
sample.f (9 KB)

Event Timeline