LHAPDF veraion 5.9.1
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf-5.9.1 / src / wrapabkm09.f
1       subroutine ABKM09evolve(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,nvar2,pdfmem,imem
6       parameter(nxb=99,nq=20,np=9,nvar=30,nvar2=25)
7       character*16 name(nmxset)
8       integer nmem(nmxset),ndef(nmxset),mem
9       common/NAME/name,nmem,ndef,mem
10       double precision gridx(nmxgridx),gridq(nmxgridq)
11       integer ngridx,ngridq,jx,jq
12 !      integer iset,iimem
13 !      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       real*8 pdfs0(0:np),pdfsp(0:np),pdfsm(0:np),q2plus,q2minus,q2zero
24       real*8 q2mat,delq2,f2p,f2pp,aq2,bq2,cq2
25       data xmin,xmax,qsqmin,qsqmax/1d-7,1d0,0.8d0,2d8/
26       data q2mat /1d0/
27
28       save f,npdf,npar,pdfmem,dels,delx,x1,delx1,xlog1,nxbb,xx &
29      &,fsp,bsp,csp,dsp
30 !      save 
31       
32       q2=Q*Q
33       if(q2.gt.qsqmax) print 99,q2
34       if(xb.lt.xmin.or.xb.gt.xmax)    print 98,xb
35
36   99  format('  WARNING:  Q^2 VALUE IS OUT OF RANGE   ',g12.3)
37   98  format('  WARNING:   X  VALUE IS OUT OF RANGE   ',g12.3)
38
39       x=max(xb,xmin)
40       x=min(x,xmax)
41       qsq=min(q2,qsqmax)
42
43       if (x.gt.x1) then
44         xd=(1d0-x1)**2-(1d0-x)**2
45         n=int(xd/delx1)+nxbb
46       else
47         xd=dlog(x)-xlog1
48         n=nxbb+int(xd/DELX)-1
49       end if
50       aa=x-xx(n)
51
52       k=pdfmem
53
54       if (qsq.ge.q2mat) then 
55         ss=dlog(dlog(qsq/0.04d0))-dlog(dlog(qsqmin/0.04d0))
56         m=int(ss/dels)+1
57         b=ss/dels-dble(m)+1.d0
58  
59         do i=1,npdf
60           f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i) & 
61      &                  + aa**3*dsp(k,n,m,i)
62           fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i) &
63      &                    + aa**3*dsp(k,n,m+1,i)
64           if (m.ge.2) then 
65             fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i) &
66      &                      +aa**3*dsp(k,n,m-1,i)
67             pdfs(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
68           else 
69             pdfs(i)= f0*(1d0-b) + fp*b
70            end if
71         end do
72       else 
73         delq2=q2mat*0.1
74
75         do i=1,npdf
76
77           q2plus=q2mat+delq2
78           ss=dlog(dlog(q2plus/0.04d0))-dlog(dlog(qsqmin/0.04d0))
79           m=int(ss/dels)+1
80           b=ss/dels-dble(m)+1.d0
81
82           f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i)  &
83      &                  + aa**3*dsp(k,n,m,i)
84           fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i) &
85      &                    + aa**3*dsp(k,n,m+1,i)
86           if (m.ge.2) then 
87             fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i) &
88      &                      +aa**3*dsp(k,n,m-1,i)
89             pdfsp(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
90           else 
91             pdfsp(i)= f0*(1d0-b) + fp*b
92            end if
93
94           q2minus=q2mat-delq2
95           ss=dlog(dlog(q2minus/0.04d0))-dlog(dlog(qsqmin/0.04d0))
96           m=int(ss/dels)+1
97           b=ss/dels-dble(m)+1.d0
98
99           f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i) & 
100      &                  + aa**3*dsp(k,n,m,i)
101           fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i) &
102      &                    + aa**3*dsp(k,n,m+1,i)
103           if (m.ge.2) then 
104             fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i) &
105      &                      +aa**3*dsp(k,n,m-1,i)
106             pdfsm(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
107           else 
108             pdfsm(i)= f0*(1d0-b) + fp*b
109            end if
110
111           q2zero=q2mat
112           ss=dlog(dlog(q2zero/0.04d0))-dlog(dlog(qsqmin/0.04d0))
113           m=int(ss/dels)+1
114           b=ss/dels-dble(m)+1.d0
115
116           f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i) & 
117      &                  + aa**3*dsp(k,n,m,i)
118           fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i) &
119      &                    + aa**3*dsp(k,n,m+1,i)
120           if (m.ge.2) then 
121             fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i) &
122      &                      +aa**3*dsp(k,n,m-1,i)
123             pdfs0(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
124           else 
125             pdfs0(i)= f0*(1d0-b) + fp*b
126            end if
127
128           f2p=(pdfsp(i)-pdfsm(i))/2./delq2
129           f2pp=(pdfsp(i)+pdfsm(i)-2*pdfs0(i))/delq2**2
130           cq2=(f2pp*q2mat**2-2*f2p*q2mat+2*pdfs0(i))/2/q2mat**3
131           bq2=f2pp/2.-3*cq2*q2mat
132           aq2=f2p-2*bq2*q2mat-3*cq2*q2mat**2
133           pdfs(i)=aq2*q2 + bq2*q2**2 + cq2*q2**3
134         end do
135       end if
136
137       fset(-6)=pdfs(9)
138       fset(-5)=pdfs(8)
139       fset(-4)=pdfs(7)
140       fset(-3)=pdfs(5)
141 !--reversed mrs 7/7/04 due
142       fset(-1)=pdfs(6)
143       fset(-2)=pdfs(4)
144       fset(0)=pdfs(3)
145 !--reversed mrs 7/7/04 due
146       fset(2)=pdfs(1)+pdfs(4)
147       fset(1)=pdfs(2)+pdfs(6)
148       fset(3)=pdfs(5)
149       fset(4)=pdfs(7)
150       fset(5)=pdfs(8)
151       fset(6)=pdfs(9)
152       return
153 !
154       entry ABKM09alfa(alfas,Q)
155       q2=Q*Q
156       if(q2.gt.qsqmax) print 99,q2
157       qsq=max(q2,q2mat)
158       qsq=min(qsq,qsqmax)
159
160       ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
161       m=int(ss/dels)+1
162       b=ss/dels-m+1
163
164       k=pdfmem
165
166       if (m.ge.2) then 
167         fm=f(k,1,m-1,0)
168         fp=f(k,1,m+1,0)
169         f0=f(k,1,m,0)
170         alfas=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
171       else 
172         alfas= f0*(1d0-b) + fp*b
173       end if
174
175       return
176 !
177       entry ABKM09getgrid(nset,ngridx,ngridq,gridx,gridq)
178      
179       ngridx=nxb
180       do jx=1,nxb
181           gridx(jx)=xx(jx)
182       enddo
183
184       ngridq=nq
185       do jq=1,ngridq
186           gridq(jq)= 0.04*exp(exp( &
187      & log(log(qsqmin/0.04))+(float(jq-1)/19)*( log(log(qsqmax/0.04))-log(log(qsqmin/0.04)) ) & 
188      & ))
189       enddo 
190        
191       return
192       
193       entry ABKM09read(nset)
194 ! following fix because members are 0-nvar
195 !      nmem = nvar + 1
196
197       read(1,fmt="(i2,i2)") npdf,npar
198       if (npar.eq.0) npar=nvar2
199
200       nmem(nset) = npar
201       ndef(nset) = 0
202
203       !print *,'number of members',npdf
204       do k=0,npar
205          do n=1,nxb-1
206             do m=1,nq
207                read(1,*) (f(k,n,m,i),i=0,npdf)
208 !              print 100,(f(k,n,m,i),i=0,npdf)
209             enddo
210          enddo
211       enddo
212   100 format (13f11.5)
213
214       nxbb=nxb/2
215       x1=0.3d0
216       xlog1=dlog(x1)
217       delx=(dlog(x1)-dlog(xmin))/dble(nxbb-1)
218       DELX1=(1.d0-x1)**2/dble(nxbb+1)
219
220 !...X GRID
221       do kx=1,nxbb
222         xx(kx)=dexp(dlog(xmin)+delx*dble(kx-1))
223       end do
224       do kx=nxbb+1,nxb-1
225         xx(kx)=1.d0-dsqrt(dabs((1.d0-x1)**2-delx1*dble(kx-nxbb)))
226       end do
227       xx(nxb)=1d0
228
229       do k=0,npar
230       do i=0,npdf
231         do m=1,nq
232           if (i.ne.0) then 
233             f(k,nxb,m,i)=0d0
234           else 
235             f(k,nxb,m,i)=f(k,nxb-1,m,i)
236           end if
237           do n=1,nxb-1
238           end do
239           do n=1,nxb
240             fsp(n)=f(k,n,m,i)
241           end do
242           call ABKM09spline (nxb,xx,fsp,bs,cs,ds)
243           do n=1,nxb
244             bsp(k,n,m,i)=bs(n)
245             csp(k,n,m,i)=cs(n)
246             dsp(k,n,m,i)=ds(n)
247           end do
248         end do
249       end do
250       end do
251
252       return
253 !
254       entry ABKM09init
255
256       dels=(dlog(dlog(qsqmax/0.04d0))- &
257      &      dlog(dlog(qsqmin/0.04d0)))/dble(nq-1)
258
259
260       return
261 !
262       entry ABKM09pdf(imem)
263       pdfmem=imem
264       if ((pdfmem.lt.0).or.(pdfmem.gt.npar)) then
265          write(*,*) 'ABKM09 PDF set:'
266          write(*,*) 'PDF member out of range:'
267          write(*,*) 'member = ',pdfmem,'    member range = (0,',npar,')'
268          stop
269       endif
270       return
271       end
272
273 ! ---------------------------------------------------------------------
274       SUBROUTINE ABKM09SPLINE(N,X,Y,B,C,D)
275 ! ---------------------------------------------------------------------
276 ! CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
277 ! INTERPOLATION SUBROUTINES ARE TAKEN FROM
278 ! G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
279 ! COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
280 !
281       IMPLICIT REAL*8(A-H,O-Z)
282
283       DIMENSION X(N), Y(N), B(N), C(N), D(N)
284 !
285       NM1=N-1
286       IF(N.LT.2) RETURN
287       IF(N.LT.3) GO TO 250
288       D(1)=X(2)-X(1)
289       C(2)=(Y(2)-Y(1))/D(1)
290       DO 210 K=2,NM1
291          D(K)=X(K+1)-X(K)
292          B(K)=2.0D0*(D(K-1)+D(K))
293          C(K+1)=(Y(K+1)-Y(K))/D(K)
294          C(K)=C(K+1)-C(K)
295   210 CONTINUE
296       B(1)=-D(1)
297       B(N)=-D(N-1)
298       C(1)=0.0D0
299       C(N)=0.0D0
300       IF(N.EQ.3) GO TO 215
301       C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
302       C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
303       C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1))
304       C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3))
305  215  CONTINUE
306       DO 220 K=2,N
307          T=D(K-1)/B(K-1)
308          B(K)=B(K)-T*D(K-1)
309          C(K)=C(K)-T*C(K-1)
310  220  CONTINUE
311       C(N)=C(N)/B(N)
312       DO 230 IB=1,NM1
313          K=N-IB
314          C(K)=(C(K)-D(K)*C(K+1))/B(K)
315  230  CONTINUE
316       B(N)=(Y(N)-Y(NM1))/D(NM1) &
317      &     +D(NM1)*(C(NM1)+2.0D0*C(N))
318       DO 240 K=1,NM1
319          B(K)=(Y(K+1)-Y(K))/D(K) &
320      &        -D(K)*(C(K+1)+2.0D0*C(K))
321          D(K)=(C(K+1)-C(K))/D(K)
322          C(K)=3.0D0*C(K)
323  240  CONTINUE
324       C(N)=3.0D0*C(N)
325       D(N)=D(N-1)
326       RETURN
327  250  CONTINUE
328       B(1)=(Y(2)-Y(1))/(X(2)-X(1))
329       C(1)=0.0D0
330       D(1)=0.0D0
331       B(2)=B(1)
332       C(2)=0.0D0
333       D(2)=0.0D0
334       RETURN
335       END