Version update.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapa02.f
CommitLineData
4e9e3152 1 subroutine A02evolve(xb,Q,fset)
2 implicit none
3 include 'parmsetup.inc'
4 integer npdf,npar,kschem,i,k,n,m,nxbb
5 integer nxb,nq,np,nvar,pdfmem,imem
6 parameter(nxb=99,nq=20,np=9,nvar=15)
7 character*16 name(nmxset)
8 integer nmem(nmxset),ndef(nmxset),mem
9 common/NAME/name,nmem,ndef,mem
10 integer nset
11 real*4 f(0:nvar,nxb,nq+1,0:np)
12 real*8 x1,xd,del,dels,delx,delx1,xlog1
13 real*8 pdfs(np),fset(-6:6),alfas
14 real*8 x,Q,xb,q2,xmin,xmax,qsq,qsqmin,qsqmax,a,b,ss
15 data xmin,xmax,qsqmin,qsqmax/1d-7,1d0,0.8d0,2d8/
16
17 save f,npdf,npar,pdfmem,dels,delx,x1,delx1,xlog1,nxbb
18
19 q2=Q*Q
20 if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
21 if(xb.lt.xmin.or.xb.gt.xmax) print 98,x
22 99 format(' WARNING: Q^2 VALUE IS OUT OF RANGE ')
23 98 format(' WARNING: X VALUE IS OUT OF RANGE ')
24
25 x=max(xb,xmin)
26 x=min(x,xmax)
27 qsq=max(q2,qsqmin)
28 qsq=min(qsq,qsqmax)
29
30 if (x.gt.x1) then
31 xd=(1-x1)**2-(1-x)**2
32 n=int(xd/delx1)+nxbb
33 a=xd/delx1-n+nxbb
34 else
35 xd=log(x)-xlog1
36 n=nxbb+int(xd/DELX)-1
37 a=xd/delx-n+nxbb
38 end if
39
40 ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
41 m=int(ss/dels)+1
42 b=ss/dels-m+1
43
44 k=pdfmem
45
46 do i=1,npdf
47 if (n.gt.1.and.m.gt.1.and.n.ne.49) then
48 pdfs(i)= f(k,n,m,i)*(1+a*b-a**2-b**2)
49 + + f(k,n+1,m+1,i)*a*b
50 + + f(k,n+1,m,i)*a*(a-2*b+1)/2.
51 + + f(k,n,m+1,i)*b*(b-2*a+1)/2.
52 + + f(k,n-1,m,i)*a*(a-1)/2.
53 + + f(k,n,m-1,i)*b*(b-1)/2.
54 else
55 pdfs(i)= (1-a)*(1-b)*f(k,n,m,i)
56 . + (1-a)*b*f(k,n,m+1,i)
57 . + a*(1-b)*f(k,n+1,m,i)
58 . + a*b*f(k,n+1,m+1,i)
59 end if
60 end do
61
62
63 fset(-6)=pdfs(9)
64 fset(-5)=pdfs(8)
65 fset(-4)=pdfs(7)
66 fset(-3)=pdfs(5)
67c--reversed mrs 7/7/04 due
68 fset(-1)=pdfs(6)
69 fset(-2)=pdfs(4)
70 fset(0)=pdfs(3)
71c--reversed mrs 7/7/04 due
72 fset(2)=pdfs(1)+pdfs(4)
73 fset(1)=pdfs(2)+pdfs(6)
74 fset(3)=pdfs(5)
75 fset(4)=pdfs(7)
76 fset(5)=pdfs(8)
77 fset(6)=pdfs(9)
78 return
79*
80 entry A02alfa(alfas,Q)
81 q2=Q*Q
82 if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
83 qsq=max(q2,qsqmin)
84 qsq=min(qsq,qsqmax)
85
86 ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
87 m=int(ss/dels)+1
88 b=ss/dels-m+1
89
90 k=pdfmem
91 alfas=(1d0-b)*f(k,1,m,0)+b*f(k,1,m+1,0)
92 return
93*
94 entry A02read(nset)
95c following fix because members are 0-nvar
96c nmem = nvar + 1
97 nmem(nset) = nvar
98 ndef(nset) = 0
99 read(1,*) npdf
100 npar=nvar
101 do k=0,npar
102 do n=1,nxb-1
103 do m=1,nq
104 read(1,100) (f(k,n,m,i),i=0,npdf)
105 enddo
106 enddo
107 enddo
108 100 format (13f11.5)
109 return
110*
111 entry A02init
112
113 dels=(log(log(qsqmax/0.04))-log(log(qsqmin/0.04)))/(nq-1)
114
115 nxbb=nxb/2
116 x1=0.3
117 xlog1=log(x1)
118 delx=(log(x1)-log(xmin))/(nxbb-1)
119 DELX1=(1-x1)**2/(nxbb+1)
120
121 do i=1,npdf
122 do m=1,nq
123 do k=1,npar
124 f(k,nxb,m,i)=0d0
125 end do
126 end do
127 do m=1,nq
128 do k=1,npar
129 f(k,nxb,m,0)=f(k,nxb-1,m,0)
130 end do
131 end do
132 end do
133 return
134*
135 entry A02pdf(imem)
136 pdfmem=imem
137 if ((pdfmem.lt.0).or.(pdfmem.gt.npar)) then
138 write(*,*) 'A02 PDF set:'
139 write(*,*) 'PDF member out of range:'
140 write(*,*) 'member = ',pdfmem,' member range = (0,',npar,')'
141 stop
142 endif
143 return
144 end