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