]>
Commit | Line | Data |
---|---|---|
4e9e3152 | 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(x,xmax) | |
37 | qsq=max(q2,qsqmin) | |
38 | qsq=min(qsq,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(qsq,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,*) (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 |