LHAPDF 5.2.2 source code.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / wrapa02m.f
1       subroutine A02Mevolve(xb,Q,fset)
2       implicit none
3       include 'parmsetup.inc'
4       integer npdf,npar,kschem,i,k,n,m,kx,nxbb
5       integer nxb,nq,np,nvar,pdfmem,imem
6       parameter(nxb=99,nq=20,np=9,nvar=17)
7       integer nexp(0:np)
8       data nexp / 0, 3, 4, 5, 5, 5, 5, 5, 5, 5 /
9       character*16 name(nmxset)
10       integer nmem(nmxset),ndef(nmxset),mem
11       common/NAME/name,nmem,ndef,mem
12 c      integer iset,iimem
13 c      common/SET/iset,iimem
14       integer nset
15       real*8 f(0:nvar,nxb,nq+1,0:np),xx(nxb)
16       real*8 fsp(nxb),bs(nxb),cs(nxb),ds(nxb)
17       real*8 bsp(0:nvar,nxb,nq+1,0:np),csp(0:nvar,nxb,nq+1,0:np)
18      +      ,dsp(0:nvar,nxb,nq+1,0:np)
19       real*8 x1,xd,del,dels,delx,delx1,xlog1
20       real*8 pdfs(0:np),fset(-6:6),alfas
21       real*8 x,Q,xb,q2,xmin,xmax,qsq,qsqmin,qsqmax,b,ss
22       real*8 aa,f0,fp,fm
23       data xmin,xmax,qsqmin,qsqmax/1d-7,1d0,0.8d0,2d8/
24
25       save f,npdf,npar,pdfmem,dels,delx,x1,delx1,xlog1,nxbb,xx
26      +,fsp,bsp,csp,dsp
27 c      save 
28       
29       q2=Q*Q
30       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
31       if(xb.lt.xmin.or.xb.gt.xmax)       print 98,x
32   99  format('  WARNING:  Q^2 VALUE IS OUT OF RANGE   ')
33   98  format('  WARNING:   X  VALUE IS OUT OF RANGE   ')
34
35       x=max(xb,xmin)
36       x=min(xb,xmax)
37       qsq=max(q2,qsqmin)
38       qsq=min(q2,qsqmax)
39
40       if (x.gt.x1) then
41         xd=(1d0-x1)**2-(1d0-x)**2
42         n=int(xd/delx1)+nxbb
43       else
44         xd=dlog(x)-xlog1
45         n=nxbb+int(xd/DELX)-1
46       end if
47       aa=x-xx(n)
48
49       ss=dlog(dlog(qsq/0.04d0))-dlog(dlog(qsqmin/0.04d0))
50       m=int(ss/dels)+1
51       b=ss/dels-dble(m)+1.d0
52  
53       k=pdfmem
54
55       do i=0,npdf
56         f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i) 
57      +              + aa**3*dsp(k,n,m,i)
58         fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i)
59      +                + aa**3*dsp(k,n,m+1,i)
60         if (m.ge.2) then 
61           fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i)
62      +                   +aa**3*dsp(k,n,m-1,i)
63           pdfs(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
64         else 
65           pdfs(i)= f0*(1d0-b) + fp*b
66          end if
67         pdfs(i) = pdfs(i)*(1d0-x)**dble(nexp(i))
68       end do
69
70       fset(-6)=pdfs(9)
71       fset(-5)=pdfs(8)
72       fset(-4)=pdfs(7)
73       fset(-3)=pdfs(5)
74 c--reversed mrs 7/7/04 due
75       fset(-1)=pdfs(6)
76       fset(-2)=pdfs(4)
77       fset(0)=pdfs(3)
78 c--reversed mrs 7/7/04 due
79       fset(2)=pdfs(1)+pdfs(4)
80       fset(1)=pdfs(2)+pdfs(6)
81       fset(3)=pdfs(5)
82       fset(4)=pdfs(7)
83       fset(5)=pdfs(8)
84       fset(6)=pdfs(9)
85       return
86 *
87       entry A02Malfa(alfas,Q)
88       q2=Q*Q
89       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
90       qsq=max(q2,qsqmin)
91       qsq=min(q2,qsqmax)
92
93       ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
94       m=int(ss/dels)+1
95       b=ss/dels-m+1
96
97       k=pdfmem
98       alfas=(1d0-b)*f(k,1,m,0)+b*f(k,1,m+1,0)
99       return
100 *
101       entry A02Mread(nset)
102 c following fix because members are 0-nvar
103 c      nmem = nvar + 1
104
105       nmem(nset) = nvar
106       ndef(nset) = 0
107
108
109       read(1,*) npdf
110 c      print *,'number of members',npdf
111       npar=nvar
112       do k=0,npar
113          do n=1,nxb-1
114             do m=1,nq
115                read(1,100) (f(k,n,m,i),i=0,npdf)
116 c              print 100,(f(k,n,m,i),i=0,npdf)
117             enddo
118          enddo
119       enddo
120   100 format (13f11.5)
121
122       nxbb=nxb/2
123       x1=0.3d0
124       xlog1=dlog(x1)
125       delx=(dlog(x1)-dlog(xmin))/dble(nxbb-1)
126       DELX1=(1.d0-x1)**2/dble(nxbb+1)
127
128 *...X GRID
129       do kx=1,nxbb
130         xx(kx)=dexp(dlog(xmin)+delx*dble(kx-1))
131       end do
132       do kx=nxbb+1,nxb-1
133         xx(kx)=1.d0-dsqrt(dabs((1.d0-x1)**2-delx1*dble(kx-nxbb)))
134       end do
135       xx(nxb)=1d0
136
137       do k=0,npar
138       do i=0,npdf
139         do m=1,nq
140           if (i.ne.0) then 
141             f(k,nxb,m,i)=0d0
142           else 
143             f(k,nxb,m,i)=f(k,nxb-1,m,i)
144           end if
145           do n=1,nxb-1
146             f(k,n,m,i)=f(k,n,m,i)/(1d0-xx(n))**nexp(i)
147           end do
148           do n=1,nxb
149             fsp(n)=f(k,n,m,i)
150           end do
151           call a02mspline (nxb,xx,fsp,bs,cs,ds)
152           do n=1,nxb
153             bsp(k,n,m,i)=bs(n)
154             csp(k,n,m,i)=cs(n)
155             dsp(k,n,m,i)=ds(n)
156           end do
157         end do
158       end do
159       end do
160
161       return
162 *
163       entry A02Minit
164
165       dels=(dlog(dlog(qsqmax/0.04d0))-
166      +      dlog(dlog(qsqmin/0.04d0)))/dble(nq-1)
167
168
169       return
170 *
171       entry A02Mpdf(imem)
172       pdfmem=imem
173       if ((pdfmem.lt.0).or.(pdfmem.gt.npar)) then
174          write(*,*) 'A02M PDF set:'
175          write(*,*) 'PDF member out of range:'
176          write(*,*) 'member = ',pdfmem,'    member range = (0,',npar,')'
177          stop
178       endif
179       return
180       end
181
182 * ---------------------------------------------------------------------
183       SUBROUTINE A02MSPLINE(N,X,Y,B,C,D)
184 * ---------------------------------------------------------------------
185 * CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
186 * INTERPOLATION SUBROUTINES ARE TAKEN FROM
187 * G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
188 * COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
189 *
190       IMPLICIT REAL*8(A-H,O-Z)
191
192       DIMENSION X(N), Y(N), B(N), C(N), D(N)
193 *
194       NM1=N-1
195       IF(N.LT.2) RETURN
196       IF(N.LT.3) GO TO 250
197       D(1)=X(2)-X(1)
198       C(2)=(Y(2)-Y(1))/D(1)
199       DO 210 K=2,NM1
200          D(K)=X(K+1)-X(K)
201          B(K)=2.0D0*(D(K-1)+D(K))
202          C(K+1)=(Y(K+1)-Y(K))/D(K)
203          C(K)=C(K+1)-C(K)
204   210 CONTINUE
205       B(1)=-D(1)
206       B(N)=-D(N-1)
207       C(1)=0.0D0
208       C(N)=0.0D0
209       IF(N.EQ.3) GO TO 215
210       C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
211       C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
212       C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1))
213       C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3))
214  215  CONTINUE
215       DO 220 K=2,N
216          T=D(K-1)/B(K-1)
217          B(K)=B(K)-T*D(K-1)
218          C(K)=C(K)-T*C(K-1)
219  220  CONTINUE
220       C(N)=C(N)/B(N)
221       DO 230 IB=1,NM1
222          K=N-IB
223          C(K)=(C(K)-D(K)*C(K+1))/B(K)
224  230  CONTINUE
225       B(N)=(Y(N)-Y(NM1))/D(NM1)
226      1     +D(NM1)*(C(NM1)+2.0D0*C(N))
227       DO 240 K=1,NM1
228          B(K)=(Y(K+1)-Y(K))/D(K)
229      1        -D(K)*(C(K+1)+2.0D0*C(K))
230          D(K)=(C(K+1)-C(K))/D(K)
231          C(K)=3.0D0*C(K)
232  240  CONTINUE
233       C(N)=3.0D0*C(N)
234       D(N)=D(N-1)
235       RETURN
236  250  CONTINUE
237       B(1)=(Y(2)-Y(1))/(X(2)-X(1))
238       C(1)=0.0D0
239       D(1)=0.0D0
240       B(2)=B(1)
241       C(2)=0.0D0
242       D(2)=0.0D0
243       RETURN
244       END