LHAPDF veraion 5.9.1
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf-5.9.1 / src / wrapabm11-lite.f
1       subroutine ABM11evolve(xb,Q,fset)
2       implicit none
3       include 'parmsetup.inc'
4       integer npdf,npar,kschem,i,k,kk,n,m,kx,nxbb,nxb1,nxb2
5       integer nxb,nq,np,nvar,pdfmem,imem,nvar2
6       parameter(nxb=99,nq=20,np=9,nvar=0,nvar2=25)
7       character*16 name(nmxset)
8       character*80 line 
9       character*512 setpath 
10       integer nmem(nmxset),ndef(nmxset),mem
11       common/NAME/name,nmem,ndef,mem
12       double precision gridx(nmxgridx),gridq(nmxgridq)
13       integer ngridx,ngridq,jx,jq
14 !      integer iset,iimem
15 !      common/SET/iset,iimem
16       integer nset,iset
17       real*8 f(0:nvar,nxb,nq+1,0:np),xx(nxb)
18       real*8 fsp(nxb),bs(nxb),cs(nxb),ds(nxb)
19       real*8 bsp(0:nvar,nxb,nq+1,0:np),csp(0:nvar,nxb,nq+1,0:np) &
20      &      ,dsp(0:nvar,nxb,nq+1,0:np)
21       real*8 x1,xd,del,dels,delx,delx1,xlog1
22       real*8 pdfs(0:np),fset(-6:6),alfas
23       real*8 x,Q,xb,q2,xmin,xmax,xmax1,qsq,qsqmin,qsqmax,b,ss
24       real*8 aa,f0,fp,fm
25       real*8 pdfs0(0:np),pdfsp(0:np),pdfsm(0:np),q2plus,q2minus,q2zero
26       real*8 q2mat,delq2,f2p,f2pp,aq2,bq2,cq2
27       integer ios
28       data xmin,xmax,qsqmin,qsqmax/1d-7,1d0,0.8d0,2d8/
29       data q2mat /1d0/
30
31       save f,npdf,npar,pdfmem,dels,delx,x1,delx1,xlog1,nxbb,xx &
32      &,fsp,bsp,csp,dsp,nxb1,nxb2,xmax1
33 !      save 
34       
35       q2=Q*Q
36       if(q2.gt.qsqmax) print 99,q2
37       if(xb.lt.xmin.or.xb.gt.xmax)    print 98,xb
38
39   99  format('  WARNING:  Q^2 VALUE IS OUT OF RANGE   ',g12.3)
40   98  format('  WARNING:   X  VALUE IS OUT OF RANGE   ',g12.3)
41
42       x=max(xb,xmin)
43       x=min(x,xmax)
44       qsq=min(q2,qsqmax)
45
46 !      if (x.gt.x1) then
47 !        xd=(1d0-x1)**2-(1d0-x)**2
48 !        n=int(xd/delx1)+nxbb
49 !      else
50 !        xd=dlog(x)-xlog1
51 !        n=nxbb+int(xd/DELX)-1
52 !      end if
53
54       do n=1,nxb-1
55         if (x.lt.xx(n+1)) goto 300
56         end do
57  300  aa=x-xx(n)
58
59 !      print *,x,xx(n),xx(n+1),aa
60       !k=pdfmem
61       k=0
62
63       if (qsq.ge.q2mat) then 
64         ss=dlog(dlog(qsq/0.04d0))-dlog(dlog(qsqmin/0.04d0))
65         m=int(ss/dels)+1
66         b=ss/dels-dble(m)+1.d0
67  
68         do i=1,npdf
69           f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i) & 
70      &                  + aa**3*dsp(k,n,m,i)
71           fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i) &
72      &                    + aa**3*dsp(k,n,m+1,i)
73           if (m.ge.2) then 
74             fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i) &
75      &                      +aa**3*dsp(k,n,m-1,i)
76             pdfs(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
77           else 
78             pdfs(i)= f0*(1d0-b) + fp*b
79            end if
80         end do
81       else 
82         delq2=q2mat*0.1
83
84         do i=1,npdf
85
86           q2plus=q2mat+delq2
87           ss=dlog(dlog(q2plus/0.04d0))-dlog(dlog(qsqmin/0.04d0))
88           m=int(ss/dels)+1
89           b=ss/dels-dble(m)+1.d0
90
91           f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i)  &
92      &                  + aa**3*dsp(k,n,m,i)
93           fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i) &
94      &                    + aa**3*dsp(k,n,m+1,i)
95           if (m.ge.2) then 
96             fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i) &
97      &                      +aa**3*dsp(k,n,m-1,i)
98             pdfsp(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
99           else 
100             pdfsp(i)= f0*(1d0-b) + fp*b
101            end if
102
103           q2minus=q2mat-delq2
104           ss=dlog(dlog(q2minus/0.04d0))-dlog(dlog(qsqmin/0.04d0))
105           m=int(ss/dels)+1
106           b=ss/dels-dble(m)+1.d0
107
108           f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i) & 
109      &                  + aa**3*dsp(k,n,m,i)
110           fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i) &
111      &                    + aa**3*dsp(k,n,m+1,i)
112           if (m.ge.2) then 
113             fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i) &
114      &                      +aa**3*dsp(k,n,m-1,i)
115             pdfsm(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
116           else 
117             pdfsm(i)= f0*(1d0-b) + fp*b
118           end if
119
120           q2zero=q2mat
121           ss=dlog(dlog(q2zero/0.04d0))-dlog(dlog(qsqmin/0.04d0))
122           m=int(ss/dels)+1
123           b=ss/dels-dble(m)+1.d0
124
125           f0=f(k,n,m,i) + aa*bsp(k,n,m,i) + aa**2*csp(k,n,m,i) & 
126      &                  + aa**3*dsp(k,n,m,i)
127           fp=f(k,n,m+1,i) + aa*bsp(k,n,m+1,i) + aa**2*csp(k,n,m+1,i) &
128      &                    + aa**3*dsp(k,n,m+1,i)
129           if (m.ge.2) then 
130             fm=f(k,n,m-1,i) + aa*bsp(k,n,m-1,i) + aa**2*csp(k,n,m-1,i) &
131      &                      +aa**3*dsp(k,n,m-1,i)
132             pdfs0(i)=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
133           else 
134             pdfs0(i)= f0*(1d0-b) + fp*b
135            end if
136
137           f2p=(pdfsp(i)-pdfsm(i))/2./delq2
138           f2pp=(pdfsp(i)+pdfsm(i)-2*pdfs0(i))/delq2**2
139           cq2=(f2pp*q2mat**2-2*f2p*q2mat+2*pdfs0(i))/2/q2mat**3
140           bq2=f2pp/2.-3*cq2*q2mat
141           aq2=f2p-2*bq2*q2mat-3*cq2*q2mat**2
142           pdfs(i)=aq2*q2 + bq2*q2**2 + cq2*q2**3
143 !          print *,i,xb,q2,aq2,bq2,cq2
144         end do
145       end if
146
147 !      fset(-6)=pdfs(9)
148 !      fset(-5)=pdfs(8)
149 !      fset(-4)=pdfs(7)
150       fset(-3)=pdfs(5)
151 !--reversed mrs 7/7/04 due
152       fset(-1)=pdfs(6)
153       fset(-2)=pdfs(4)
154       fset(0)=pdfs(3)
155 !--reversed mrs 7/7/04 due
156       fset(2)=pdfs(1)+pdfs(4)
157       fset(1)=pdfs(2)+pdfs(6)
158       fset(3)=pdfs(5)
159 !      fset(4)=pdfs(7)
160 !      fset(5)=pdfs(8)
161 !      fset(6)=pdfs(9)
162
163       do kk=7,npdf
164         fset(kk-3)=pdfs(kk)
165         fset(-kk+3)=pdfs(kk)
166       end do
167       do kk=npdf+1,9
168         fset(kk-3)=0.
169         fset(-kk+3)=0.
170       end do
171
172       return
173 !
174       entry ABM11alfa(alfas,Q)
175       q2=Q*Q
176       if(q2.gt.qsqmax) print 99,q2
177       qsq=max(q2,q2mat)
178       qsq=min(qsq,qsqmax)
179
180       ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
181       m=int(ss/dels)+1
182       b=ss/dels-m+1
183
184       k=pdfmem
185
186       fp=f(k,1,m+1,0)
187       f0=f(k,1,m,0)
188       if (m.ge.2) then 
189         fm=f(k,1,m-1,0)
190         alfas=fm*b*(b-1d0)/2d0 + f0*(1d0-b**2) + fp*b*(b+1d0)/2d0
191       else 
192         alfas= f0*(1d0-b) + fp*b
193       end if
194
195 !      print *,q2,q2mat,m,alfas
196
197
198       return
199 !
200       entry ABM11getgrid(nset,ngridx,ngridq,gridx,gridq)
201      
202       ngridx=nxb
203       do jx=1,nxb
204           gridx(jx)=xx(jx)
205       enddo
206
207       ngridq=nq
208       do jq=1,ngridq
209           gridq(jq)= 0.04*exp(exp( &
210      & log(log(qsqmin/0.04))+(float(jq-1)/19)*( log(log(qsqmax/0.04))-log(log(qsqmin/0.04)) ) & 
211      & ))
212       enddo 
213        
214       return
215 !
216       entry ABM11read(nset)
217 ! following fix because members are 0-nvar
218 !      nmem = nvar + 1
219
220       nmem(nset) = nvar
221       ndef(nset) = 0
222
223 ! - dummy read in to get to End: (stream 1 is still open)               
224       read(1,fmt="(i2,i2)") npdf,npar
225 !      print *,'number of members',npdf,npar
226       if (npar.eq.0) npar=nvar2
227
228       do k=0,npar
229          do n=1,nxb-1
230             do m=1,nq
231                read(1,*) (f(k,n,m,i),i=0,npdf)
232 !              print 100,k,n,m,i,(f(k,n,m,i),i=0,npdf)
233             enddo
234          enddo
235       enddo
236   100 format (4i4,13f11.5)
237
238       return
239 !
240       entry ABM11init
241
242       dels=(dlog(dlog(qsqmax/0.04d0))- &
243      &      dlog(dlog(qsqmin/0.04d0)))/dble(nq-1)
244
245
246       return
247 !
248       entry ABM11pdf(imem)
249       pdfmem=imem
250       if ((pdfmem.lt.0).or.(pdfmem.gt.npar)) then
251          write(*,*) 'ABM11 PDF set:'
252          write(*,*) 'PDF member out of range:'
253          write(*,*) 'member = ',pdfmem,'    member range = (0,',npar,')'
254          stop
255       endif
256
257       call getnset(iset)
258       nmem(iset) = npar 
259       ndef(iset) = 0 
260       
261 ! have to reopen stream 1                                               
262       call getsetpath(setpath) 
263       open(1,file=setpath(1:len_trim(setpath)),action='READ') 
264       line = '' 
265       do while (line(2:11).ne.'Evolution:') 
266          read(1,'(a)') line 
267       enddo 
268
269       read(1,'(a)'),line 
270       read(1,'(a)'),line 
271
272       read(1,fmt="(i2,i2)") npdf,npar
273       if (npar.eq.0) npar=nvar2
274
275       do k=0,imem-1
276          do n=1,nxb-1
277             do m=1,nq
278                read(1,*) (f(0,n,m,i),i=0,npdf)
279             enddo
280          enddo
281       enddo
282       
283       do n=1,nxb-1 
284          do m=1,nq 
285             read(1,*) (f(0,n,m,i),i=0,npdf) 
286          enddo 
287       enddo 
288
289 !...X GRID
290
291       nxb1=30
292       nxb2=nxb-nxb1
293       x1=0.3d0
294       xmax1=0.99d0
295       delx=(log(x1)-log(xmin))/(nxb1-1)
296       DELX1=(log(1-xmax1)-log(1-x1))/(nxb2-1)
297
298       do kx=1,nxb1
299         xx(kx)=exp(log(xmin)+delx*(kx-1))
300       end do
301
302       do kx=nxb-1,nxb1,-1
303         xx(kx)=1-exp(log(1-xmax1)+delx1*(kx+1-nxb))
304       end do
305
306 !      nxbb=nxb/2
307 !      x1=0.3d0
308 !      xlog1=dlog(x1)
309 !      delx=(dlog(x1)-dlog(xmin))/dble(nxbb-1)
310 !      DELX1=(1.d0-x1)**2/dble(nxbb+1)
311 !      do kx=1,nxbb
312 !        xx(kx)=dexp(dlog(xmin)+delx*dble(kx-1))
313 !      end do
314 !      do kx=nxbb+1,nxb-1
315 !        xx(kx)=1.d0-dsqrt(dabs((1.d0-x1)**2-delx1*dble(kx-nxbb)))
316 !      end do
317
318       xx(nxb)=1d0
319
320       do i=0,npdf
321         do m=1,nq
322           if (i.ne.0) then 
323             f(0,nxb,m,i)=0d0
324           else 
325             f(0,nxb,m,i)=f(0,nxb-1,m,i)
326           end if
327           do n=1,nxb-1
328           end do
329           do n=1,nxb
330             fsp(n)=f(0,n,m,i)
331           end do
332           call ABM11spline (nxb,xx,fsp,bs,cs,ds)
333           do n=1,nxb
334             bsp(0,n,m,i)=bs(n)
335             csp(0,n,m,i)=cs(n)
336             dsp(0,n,m,i)=ds(n)
337           end do
338         end do
339       end do
340       close(1)
341       return
342       end
343
344 ! ---------------------------------------------------------------------
345       SUBROUTINE ABM11SPLINE(N,X,Y,B,C,D)
346 ! ---------------------------------------------------------------------
347 ! CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
348 ! INTERPOLATION SUBROUTINES ARE TAKEN FROM
349 ! G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
350 ! COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
351 !
352       IMPLICIT REAL*8(A-H,O-Z)
353
354       DIMENSION X(N), Y(N), B(N), C(N), D(N)
355 !
356       NM1=N-1
357       IF(N.LT.2) RETURN
358       IF(N.LT.3) GO TO 250
359       D(1)=X(2)-X(1)
360       C(2)=(Y(2)-Y(1))/D(1)
361       DO 210 K=2,NM1
362          D(K)=X(K+1)-X(K)
363          B(K)=2.0D0*(D(K-1)+D(K))
364          C(K+1)=(Y(K+1)-Y(K))/D(K)
365          C(K)=C(K+1)-C(K)
366   210 CONTINUE
367       B(1)=-D(1)
368       B(N)=-D(N-1)
369       C(1)=0.0D0
370       C(N)=0.0D0
371       IF(N.EQ.3) GO TO 215
372       C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
373       C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
374       C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1))
375       C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3))
376  215  CONTINUE
377       DO 220 K=2,N
378          T=D(K-1)/B(K-1)
379          B(K)=B(K)-T*D(K-1)
380          C(K)=C(K)-T*C(K-1)
381  220  CONTINUE
382       C(N)=C(N)/B(N)
383       DO 230 IB=1,NM1
384          K=N-IB
385          C(K)=(C(K)-D(K)*C(K+1))/B(K)
386  230  CONTINUE
387       B(N)=(Y(N)-Y(NM1))/D(NM1) &
388      &     +D(NM1)*(C(NM1)+2.0D0*C(N))
389       DO 240 K=1,NM1
390          B(K)=(Y(K+1)-Y(K))/D(K) &
391      &        -D(K)*(C(K+1)+2.0D0*C(K))
392          D(K)=(C(K+1)-C(K))/D(K)
393          C(K)=3.0D0*C(K)
394  240  CONTINUE
395       C(N)=3.0D0*C(N)
396       D(N)=D(N-1)
397       RETURN
398  250  CONTINUE
399       B(1)=(Y(2)-Y(1))/(X(2)-X(1))
400       C(1)=0.0D0
401       D(1)=0.0D0
402       B(2)=B(1)
403       C(2)=0.0D0
404       D(2)=0.0D0
405       RETURN
406       END