| 3c5d1739 |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame