| 0caf84a5 |
1 | ! -*- F90 -*- |
| 2 | |
| 3 | |
| 4 | subroutine CTEQ6evolve(x,Q,pdf) |
| 5 | implicit real*8(a-h,o-z) |
| 6 | include 'parmsetup.inc' |
| 7 | character*16 name(nmxset) |
| 8 | integer nmem(nmxset),ndef(nmxset),mmem |
| 9 | common/NAME/name,nmem,ndef,mmem |
| 10 | real*8 pdf(-6:6) |
| 11 | integer nset |
| 12 | Character Line*80 |
| 13 | PARAMETER (MXX = 96, MXQ = 20, MXF = 5, nhess = 40) |
| 14 | PARAMETER (MXPQX = (MXF + 3) * MXQ * MXX) |
| 15 | Common & |
| 16 | & / CtqPar1nhess / & |
| 17 | & Al(nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset), & |
| 18 | & UPD(0:nhess,MXPQX,nmxset) & |
| 19 | & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset) & |
| 20 | & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset) & |
| 21 | & / QCDtable / Alambda(nmxset), Nfl(nmxset), Iorder(nmxset) & |
| 22 | & / Masstbl / Amass(6,nmxset) |
| 23 | common/masses_LHA/cMass(nmxset),bMass(nmxset),tMass(nmxset) |
| 24 | data pi / 3.141592653589793d0 / |
| 25 | save |
| 26 | ! |
| 27 | call getnset(iset) |
| 28 | call getnmem(iset,imem) |
| 29 | ! |
| 30 | U = X * CtLhCtq6Pdf(imem,1,X,Q) |
| 31 | D = X * CtLhCtq6Pdf(imem,2,X,Q) |
| 32 | USEA = X * CtLhCtq6Pdf(imem,-1,X,Q) |
| 33 | DSEA = X * CtLhCtq6Pdf(imem,-2,X,Q) |
| 34 | STR = X * CtLhCtq6Pdf(imem,3,X,Q) |
| 35 | CHM = X * CtLhCtq6Pdf(imem,4,X,Q) |
| 36 | BOT = X * CtLhCtq6Pdf(imem,5,X,Q) |
| 37 | GLU = X * CtLhCtq6Pdf(imem,0,X,Q) |
| 38 | ! |
| 39 | pdf(0) = glu |
| 40 | pdf(1) = d |
| 41 | pdf(-1) = dsea |
| 42 | pdf(2) = u |
| 43 | pdf(-2) = usea |
| 44 | pdf(3) = str |
| 45 | pdf(-3) = str |
| 46 | pdf(4) = chm |
| 47 | pdf(-4) = chm |
| 48 | pdf(5) = bot |
| 49 | pdf(-5) = bot |
| 50 | pdf(6) = 0.0d0 |
| 51 | pdf(-6) = 0.0d0 |
| 52 | |
| 53 | return |
| 54 | ! |
| 55 | entry CTEQ6read(nset) |
| 56 | |
| 57 | call CtLhbldat1 |
| 58 | call CtLhbldat2 |
| 59 | |
| 60 | call LHct6set |
| 61 | !*** nmem+1=number of member |
| 62 | read(1,*)nmem(nset),ndef(nset) |
| 63 | if(nmem(nset) .gt. nhess) then |
| 64 | print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess |
| 65 | stop |
| 66 | endif |
| 67 | Read (1, '(A)') Line |
| 68 | Read (1, '(A)') Line |
| 69 | Read (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6) |
| 70 | Iorder(nset) = Nint(Dr) |
| 71 | Nfl(nset) = Nint(Fl) |
| 72 | Alambda(nset) = Al(nset) |
| 73 | |
| 74 | cMass(nset) = Amass(4,nset) |
| 75 | bMass(nset) = Amass(5,nset) |
| 76 | tMass(nset) = Amass(6,nset) |
| 77 | |
| 78 | Read (1, '(A)') Line |
| 79 | Read (1, *) NX(nset), NT(nset), NfMx(nset) |
| 80 | |
| 81 | Read (1, '(A)') Line |
| 82 | Read (1, *) & |
| 83 | & QINI(nset), QMAX(nset), (TV(I,nset), I =0, NT(nset)) |
| 84 | |
| 85 | Read (1, '(A)') Line |
| 86 | Read (1, *) XMIN(nset), (XV(I,nset), I =0, NX(nset)) |
| 87 | |
| 88 | Do 11 Iq = 0, NT(nset) |
| 89 | TV(Iq,nset) = Log(Log (TV(Iq,nset) /Al(nset))) |
| 90 | 11 Continue |
| 91 | ! |
| 92 | ! Since quark = anti-quark for nfl>2 at this stage, |
| 93 | ! we Read out only the non-redundent data points |
| 94 | ! No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence) |
| 95 | |
| 96 | Nblk = (NX(nset)+1) * (NT(nset)+1) |
| 97 | Npts = Nblk * (NfMx(nset)+3) |
| 98 | |
| 99 | !*** new version: allows nm |
| 100 | do ihess = 0,nmem(nset) |
| 101 | |
| 102 | Read (1, '(A)') Line |
| 103 | Read (1, '(A)') Line |
| 104 | Read (1, *, IOSTAT=IRET) (UPD(ihess,I,nset), I=1,Npts) |
| 105 | |
| 106 | enddo |
| 107 | return |
| 108 | ! |
| 109 | entry CTEQ6alfa(alfas,Qalfa) |
| 110 | alfas = pi*CtLhALPI(Qalfa) |
| 111 | return |
| 112 | ! |
| 113 | entry CTEQ6init(Eorder,Q2fit) |
| 114 | |
| 115 | return |
| 116 | ! |
| 117 | entry CTEQ6pdf(mem) |
| 118 | ! imem = mem |
| 119 | call getnset(iset) |
| 120 | call setnmem(iset,mem) |
| 121 | return |
| 122 | ! |
| 123 | END |
| 124 | |
| 125 | subroutine LHct6set |
| 126 | Implicit Double Precision (A-H,O-Z) |
| 127 | common /ctq6co/ xlast, qlast, nxsave |
| 128 | nxsave = -1000 |
| 129 | xlast = -2. |
| 130 | qlast = -2. |
| 131 | return |
| 132 | END |
| 133 | |
| 134 | !======================================================================= |
| 135 | Function CtLhPartonX6 (imem,IPRTN, XX, QQ) |
| 136 | ! Given the parton distribution function in the array U in |
| 137 | ! COMMON / PEVLDT / , this routine interpolates to find |
| 138 | ! the parton distribution at an arbitray point in x and q. |
| 139 | ! |
| 140 | Implicit Double Precision (A-H,O-Z) |
| 141 | |
| 142 | include 'parmsetup.inc' |
| 143 | Parameter (MXX = 96, MXQ = 20, MXF = 5, nhess = 40) |
| 144 | Parameter (MXQX= MXQ * MXX, MXPQX = MXQX * (MXF+3)) |
| 145 | |
| 146 | Common & |
| 147 | & / CtqPar1nhess / & |
| 148 | & Al(nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset), & |
| 149 | & UPD(0:nhess,MXPQX,nmxset) & |
| 150 | & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset) & |
| 151 | & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset) |
| 152 | |
| 153 | common /ctq6co/ xlast, qlast, nxsave |
| 154 | parameter(nqvec = 4) |
| 155 | |
| 156 | Dimension fvec(4), fij(4) |
| 157 | Dimension xvpow(0:mxx) |
| 158 | Data OneP / 1.00001 / |
| 159 | !**** choice of interpolation variable |
| 160 | Data xpow / 0.3d0 / |
| 161 | Save xvpow |
| 162 | Data ixprint,iqprint/0,0/ |
| 163 | save ixprint,iqprint |
| 164 | |
| 165 | save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3 |
| 166 | save const1 , const2, const3, const4, const5, const6 |
| 167 | save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet |
| 168 | |
| 169 | |
| 170 | call getnset(iset) |
| 171 | call getnmem(iset,imem) |
| 172 | |
| 173 | ! store the powers used for interpolation on first call... |
| 174 | if(nx(iset).ne. nxsave) then |
| 175 | nxsave = nx(iset) |
| 176 | xvpow(0) = 0.D0 |
| 177 | do i = 1, nx(iset) |
| 178 | xvpow(i) = xv(i,iset)**xpow |
| 179 | enddo |
| 180 | endif |
| 181 | |
| 182 | X = XX |
| 183 | Q = QQ |
| 184 | |
| 185 | if((x.lt.xmin(iset)).or.(x.gt.1.d0)) then |
| 186 | ixprint=ixprint+1 |
| 187 | if(ixprint.lt.11) print 98,x |
| 188 | if(ixprint.eq.10) print *, & |
| 189 | & 'more warning messages like the last suppressed.' |
| 190 | endif |
| 191 | 98 format(' WARNING: X=',e12.5,' OUT OF RANGE') |
| 192 | if((q.lt.qini(iset)).or.(q.gt.qmax(iset))) then |
| 193 | iqprint=iqprint+1 |
| 194 | if(iqprint.lt.11) print 99,q |
| 195 | if(iqprint.eq.10) print *, & |
| 196 | & 'more warning messages like the last suppressed.' |
| 197 | endif |
| 198 | 99 format(' WARNING: Q=',e12.5,' OUT OF RANGE') |
| 199 | |
| 200 | ! skip the initialization in x if same as in the previous call. |
| 201 | if(x .eq. xlast) goto 100 |
| 202 | xlast = x |
| 203 | |
| 204 | ! ------------- find lower end of interval containing x, i.e., |
| 205 | ! get jx such that xv(jx) .le. x .le. xv(jx+1)... |
| 206 | JLx = -1 |
| 207 | JU = Nx(iset)+1 |
| 208 | 11 If (JU-JLx .GT. 1) Then |
| 209 | JM = (JU+JLx) / 2 |
| 210 | If (X .Ge. XV(JM,iset)) Then |
| 211 | JLx = JM |
| 212 | Else |
| 213 | JU = JM |
| 214 | Endif |
| 215 | Goto 11 |
| 216 | Endif |
| 217 | ! Ix 0 1 2 Jx JLx Nx-2 Nx |
| 218 | ! |---|---|---|...|---|-x-|---|...|---|---| |
| 219 | ! x 0 Xmin x 1 |
| 220 | ! |
| 221 | If (JLx .LE. -1) Then |
| 222 | Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX6 x=',x |
| 223 | Stop |
| 224 | ElseIf (JLx .Eq. 0) Then |
| 225 | Jx = 0 |
| 226 | Elseif (JLx .LE. Nx(iset)-2) Then |
| 227 | |
| 228 | ! For interior points, keep x in the middle, as shown abo |
| 229 | Jx = JLx - 1 |
| 230 | Elseif (JLx.Eq.Nx(iset)-1 .or. x.LT.OneP) Then |
| 231 | |
| 232 | ! We tolerate a slight over-shoot of one (OneP=1.00001) |
| 233 | ! perhaps due to roundoff or whatever, but not more than th |
| 234 | ! Keep at least 4 points >= Jx |
| 235 | Jx = JLx - 2 |
| 236 | Else |
| 237 | Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX6 x=',x |
| 238 | Stop |
| 239 | Endif |
| 240 | ! ---------- Note: JLx uniquely identifies the x-bin; Jx does n |
| 241 | |
| 242 | ! This is the variable to be interpolated in |
| 243 | ss = x**xpow |
| 244 | |
| 245 | If (JLx.Ge.2 .and. JLx.Le.Nx(iset)-2) Then |
| 246 | |
| 247 | ! initiation work for "interior bins": store the lattice points in s |
| 248 | svec1 = xvpow(jx) |
| 249 | svec2 = xvpow(jx+1) |
| 250 | svec3 = xvpow(jx+2) |
| 251 | svec4 = xvpow(jx+3) |
| 252 | |
| 253 | s12 = svec1 - svec2 |
| 254 | s13 = svec1 - svec3 |
| 255 | s23 = svec2 - svec3 |
| 256 | s24 = svec2 - svec4 |
| 257 | s34 = svec3 - svec4 |
| 258 | |
| 259 | sy2 = ss - svec2 |
| 260 | sy3 = ss - svec3 |
| 261 | |
| 262 | ! constants needed for interpolating in s at fixed t lattice points... |
| 263 | const1 = s13/s23 |
| 264 | const2 = s12/s23 |
| 265 | const3 = s34/s23 |
| 266 | const4 = s24/s23 |
| 267 | s1213 = s12 + s13 |
| 268 | s2434 = s24 + s34 |
| 269 | sdet = s12*s34 - s1213*s2434 |
| 270 | tmp = sy2*sy3/sdet |
| 271 | const5 = (s34*sy2-s2434*sy3)*tmp/s12 |
| 272 | const6 = (s1213*sy2-s12*sy3)*tmp/s34 |
| 273 | |
| 274 | EndIf |
| 275 | |
| 276 | 100 continue |
| 277 | |
| 278 | ! skip the initialization in q if same as in the previous call. |
| 279 | if(q .eq. qlast) goto 110 |
| 280 | qlast = q |
| 281 | |
| 282 | tt = log(log(Q/Al(iset))) |
| 283 | |
| 284 | ! --------------Now find lower end of interval containing Q, i.e |
| 285 | ! get jq such that qv(jq) .le. q .le. qv(jq+1). |
| 286 | JLq = -1 |
| 287 | JU = NT(iset)+1 |
| 288 | 12 If (JU-JLq .GT. 1) Then |
| 289 | JM = (JU+JLq) / 2 |
| 290 | If (tt .GE. TV(JM,iset)) Then |
| 291 | JLq = JM |
| 292 | Else |
| 293 | JU = JM |
| 294 | Endif |
| 295 | Goto 12 |
| 296 | Endif |
| 297 | |
| 298 | If (JLq .LE. 0) Then |
| 299 | Jq = 0 |
| 300 | Elseif (JLq .LE. Nt(iset)-2) Then |
| 301 | ! keep q in the middle, as shown above |
| 302 | Jq = JLq - 1 |
| 303 | Else |
| 304 | ! JLq .GE. Nt-1 case: Keep at least 4 points >= |
| 305 | Jq = Nt(iset) - 3 |
| 306 | |
| 307 | Endif |
| 308 | ! This is the interpolation variable i |
| 309 | |
| 310 | If (JLq.GE.1 .and. JLq.LE.Nt(iset)-2) Then |
| 311 | ! store the lattice points in t.. |
| 312 | tvec1 = Tv(jq,iset) |
| 313 | tvec2 = Tv(jq+1,iset) |
| 314 | tvec3 = Tv(jq+2,iset) |
| 315 | tvec4 = Tv(jq+3,iset) |
| 316 | |
| 317 | t12 = tvec1 - tvec2 |
| 318 | t13 = tvec1 - tvec3 |
| 319 | t23 = tvec2 - tvec3 |
| 320 | t24 = tvec2 - tvec4 |
| 321 | t34 = tvec3 - tvec4 |
| 322 | |
| 323 | ty2 = tt - tvec2 |
| 324 | ty3 = tt - tvec3 |
| 325 | |
| 326 | tmp1 = t12 + t13 |
| 327 | tmp2 = t24 + t34 |
| 328 | |
| 329 | tdet = t12*t34 - tmp1*tmp2 |
| 330 | |
| 331 | EndIf |
| 332 | |
| 333 | 110 continue |
| 334 | |
| 335 | ! get the pdf function values at the lattice points... |
| 336 | ! In this code, we store 8 flavors: u,ubar,d,dbar,s=sbar,c=cbar,b=bbar,g |
| 337 | ! hence Iprtn=3,4,5 (s,c,b) are obtained from -3,-4,-5 (sbar,cbar,bbar) |
| 338 | |
| 339 | If (Iprtn .GE. 3) Then |
| 340 | Ip = - Iprtn |
| 341 | Else |
| 342 | Ip = Iprtn |
| 343 | EndIf |
| 344 | jtmp = ((Ip + NfMx(iset))*(NT(iset)+1)+(jq-1))*(NX(iset)+1)+jx+1 |
| 345 | |
| 346 | Do it = 1, nqvec |
| 347 | |
| 348 | J1 = jtmp + it*(NX(iset)+1) |
| 349 | |
| 350 | If (Jx .Eq. 0) Then |
| 351 | ! For the first 4 x points, interpolate x^2*f(x |
| 352 | ! This applies to the two lowest bins JLx = 0, |
| 353 | ! We cannot put the JLx.eq.1 bin into the "interior" section |
| 354 | ! (as we do for q), since Upd(J1) is undefined |
| 355 | fij(1) = 0 |
| 356 | fij(2) = Upd(imem,J1+1,iset) * XV(1,iset)**2 |
| 357 | fij(3) = Upd(imem,J1+2,iset) * XV(2,iset)**2 |
| 358 | fij(4) = Upd(imem,J1+3,iset) * XV(3,iset)**2 |
| 359 | ! |
| 360 | ! Use CtLhPolint which allows x to be anywhere w.r.t. th |
| 361 | |
| 362 | Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx) |
| 363 | |
| 364 | If (x .GT. 0D0) Fvec(it) = Fx / x**2 |
| 365 | ! Pdf is undefined for x.eq |
| 366 | ElseIf (JLx .Eq. Nx(iset)-1) Then |
| 367 | ! This is the highest x b |
| 368 | |
| 369 | !** fix allow 4 consecutive elements with imem... mrw 19.9.2005 |
| 370 | fij(1) = Upd(imem,j1,iset) |
| 371 | fij(2) = Upd(imem,j1+1,iset) |
| 372 | fij(3) = Upd(imem,j1+2,iset) |
| 373 | fij(4) = Upd(imem,j1+3,iset) |
| 374 | Call CtLhPolint4 (XVpow(Nx(iset)-3), Fij(1), 4, ss, Fx, Dfx) |
| 375 | |
| 376 | Fvec(it) = Fx |
| 377 | |
| 378 | Else |
| 379 | |
| 380 | ! for all interior points, use Jon's in-line function |
| 381 | ! This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2) |
| 382 | ! (This is cubic spline interpolation, as used by cteq; it was |
| 383 | ! changed to polint in previous Durham releases (jcp).) |
| 384 | sf2 = Upd(imem,J1+1,iset) |
| 385 | sf3 = Upd(imem,J1+2,iset) |
| 386 | |
| 387 | Fvec(it) = (const5*(Upd(imem,J1,iset) & |
| 388 | & - sf2*const1 + sf3*const2) & |
| 389 | & + const6*(Upd(imem,J1+3,iset) & |
| 390 | & + sf2*const3 - sf3*const4) & |
| 391 | & + sf2*sy3 - sf3*sy2) / s23 |
| 392 | |
| 393 | Endif |
| 394 | |
| 395 | enddo |
| 396 | ! We now have the four values Fvec(1:4 |
| 397 | ! interpolate in t... |
| 398 | |
| 399 | If (JLq .LE. 0) Then |
| 400 | ! 1st Q-bin, as well as extrapolation to lower Q |
| 401 | Call CtLhPolint4(TV(0,iset), Fvec(1), 4, tt, ff, Dfq) |
| 402 | |
| 403 | ElseIf (JLq .GE. Nt(iset)-1) Then |
| 404 | ! Last Q-bin, as well as extrapolation to higher |
| 405 | Call CtLhPolint4(TV(Nt(iset)-3,iset), Fvec(1), 4, tt, ff, Dfq) |
| 406 | Else |
| 407 | ! Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2) |
| 408 | ! which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for |
| 409 | ! the full range QV(0:Nt) (in contrast to XV) |
| 410 | tf2 = fvec(2) |
| 411 | tf3 = fvec(3) |
| 412 | |
| 413 | g1 = ( tf2*t13 - tf3*t12) / t23 |
| 414 | g4 = (-tf2*t34 + tf3*t24) / t23 |
| 415 | |
| 416 | h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12 & |
| 417 | & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34) |
| 418 | |
| 419 | ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23 |
| 420 | EndIf |
| 421 | |
| 422 | CtLhPartonX6 = ff |
| 423 | |
| 424 | Return |
| 425 | ! ******************** |
| 426 | END |
| 427 | !======================================================================= |
| 428 | Function CtLhCtq6Pdf (imem,Iparton, X, Q) |
| 429 | Implicit Double Precision (A-H,O-Z) |
| 430 | include 'parmsetup.inc' |
| 431 | Logical Warn |
| 432 | Common & |
| 433 | & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset) & |
| 434 | & / QCDtable / Alambda(nmxset), Nfl(nmxset), Iorder(nmxset) |
| 435 | |
| 436 | Data Warn /.true./ |
| 437 | save Warn |
| 438 | |
| 439 | call getnset(iset) |
| 440 | |
| 441 | If (X .lt. 0D0 .or. X .gt. 1D0) Then |
| 442 | Print *, 'X out of range in CtLhCtq6Pdf: ', X |
| 443 | Stop |
| 444 | Endif |
| 445 | If (Q .lt. Alambda(iset)) Then |
| 446 | Print *, 'Q out of range in CtLhCtq6Pdf: ', Q |
| 447 | Stop |
| 448 | Endif |
| 449 | |
| 450 | ! added to force pdf = 0.0 at x=1.0 exactly - mrw |
| 451 | if(x .eq. 1.0d0) then |
| 452 | CtLhCtq6Pdf = 0.0d0 |
| 453 | return |
| 454 | endif |
| 455 | ! |
| 456 | If ((Iparton .lt. -NfMx(iset) .or. Iparton .gt. NfMx(iset))) Then |
| 457 | If (Warn) Then |
| 458 | ! put a warning for calling extra flavor. |
| 459 | Warn = .false. |
| 460 | Print *, 'Warning: Iparton out of range in CtLhCtq6Pdf: ' & |
| 461 | & , Iparton |
| 462 | Endif |
| 463 | CtLhCtq6Pdf = 0D0 |
| 464 | Return |
| 465 | Endif |
| 466 | |
| 467 | CtLhCtq6Pdf = CtLhPartonX6 (imem,Iparton, X, Q) |
| 468 | if(CtLhCtq6Pdf.lt.0.D0) CtLhCtq6Pdf = 0.D0 |
| 469 | |
| 470 | Return |
| 471 | |
| 472 | ! ******************** |
| 473 | END |
| 474 | SUBROUTINE CtLhPOLINT4 (XA,YA,N,X,Y,DY) |
| 475 | ! fast version of polint, valid only for N=4 |
| 476 | ! Have explicitly unrolled the loops. |
| 477 | |
| 478 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) |
| 479 | |
| 480 | PARAMETER (NMAX=4) |
| 481 | DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) |
| 482 | |
| 483 | if(n .ne. 4) then |
| 484 | print *,'fatal CtLhPolint4 call',n |
| 485 | stop |
| 486 | endif |
| 487 | |
| 488 | NS=1 |
| 489 | DIF=ABS(X-XA(1)) |
| 490 | |
| 491 | DIFT=ABS(X-XA(1)) |
| 492 | IF (DIFT.LT.DIF) THEN |
| 493 | NS=1 |
| 494 | DIF=DIFT |
| 495 | ENDIF |
| 496 | C(1)=YA(1) |
| 497 | D(1)=YA(1) |
| 498 | |
| 499 | DIFT=ABS(X-XA(2)) |
| 500 | IF (DIFT.LT.DIF) THEN |
| 501 | NS=2 |
| 502 | DIF=DIFT |
| 503 | ENDIF |
| 504 | C(2)=YA(2) |
| 505 | D(2)=YA(2) |
| 506 | |
| 507 | DIFT=ABS(X-XA(3)) |
| 508 | IF (DIFT.LT.DIF) THEN |
| 509 | NS=3 |
| 510 | DIF=DIFT |
| 511 | ENDIF |
| 512 | C(3)=YA(3) |
| 513 | D(3)=YA(3) |
| 514 | |
| 515 | DIFT=ABS(X-XA(4)) |
| 516 | IF (DIFT.LT.DIF) THEN |
| 517 | NS=4 |
| 518 | DIF=DIFT |
| 519 | ENDIF |
| 520 | C(4)=YA(4) |
| 521 | D(4)=YA(4) |
| 522 | |
| 523 | |
| 524 | Y=YA(NS) |
| 525 | NS=NS-1 |
| 526 | |
| 527 | |
| 528 | HO=XA(1)-X |
| 529 | HP=XA(2)-X |
| 530 | W=C(2)-D(1) |
| 531 | DEN=W/(HO-HP) |
| 532 | D(1)=HP*DEN |
| 533 | C(1)=HO*DEN |
| 534 | |
| 535 | |
| 536 | HO=XA(2)-X |
| 537 | HP=XA(3)-X |
| 538 | W=C(3)-D(2) |
| 539 | DEN=W/(HO-HP) |
| 540 | D(2)=HP*DEN |
| 541 | C(2)=HO*DEN |
| 542 | |
| 543 | |
| 544 | HO=XA(3)-X |
| 545 | HP=XA(4)-X |
| 546 | W=C(4)-D(3) |
| 547 | DEN=W/(HO-HP) |
| 548 | D(3)=HP*DEN |
| 549 | C(3)=HO*DEN |
| 550 | |
| 551 | IF (2*NS.LT.3)THEN |
| 552 | DY=C(NS+1) |
| 553 | ELSE |
| 554 | DY=D(NS) |
| 555 | NS=NS-1 |
| 556 | ENDIF |
| 557 | Y=Y+DY |
| 558 | |
| 559 | HO=XA(1)-X |
| 560 | HP=XA(3)-X |
| 561 | W=C(2)-D(1) |
| 562 | DEN=W/(HO-HP) |
| 563 | D(1)=HP*DEN |
| 564 | C(1)=HO*DEN |
| 565 | |
| 566 | HO=XA(2)-X |
| 567 | HP=XA(4)-X |
| 568 | W=C(3)-D(2) |
| 569 | DEN=W/(HO-HP) |
| 570 | D(2)=HP*DEN |
| 571 | C(2)=HO*DEN |
| 572 | |
| 573 | IF (2*NS.LT.2)THEN |
| 574 | DY=C(NS+1) |
| 575 | ELSE |
| 576 | DY=D(NS) |
| 577 | NS=NS-1 |
| 578 | ENDIF |
| 579 | Y=Y+DY |
| 580 | |
| 581 | HO=XA(1)-X |
| 582 | HP=XA(4)-X |
| 583 | W=C(2)-D(1) |
| 584 | DEN=W/(HO-HP) |
| 585 | D(1)=HP*DEN |
| 586 | C(1)=HO*DEN |
| 587 | |
| 588 | IF (2*NS.LT.4-3)THEN |
| 589 | DY=C(NS+1) |
| 590 | ELSE |
| 591 | DY=D(NS) |
| 592 | NS=NS-1 |
| 593 | ENDIF |
| 594 | Y=Y+DY |
| 595 | |
| 596 | RETURN |
| 597 | END |
| 598 | SUBROUTINE CTLHPOLINT3 (XA,YA,N,X,Y,DY) |
| 599 | ! fast version of polint, valid only for N=3 |
| 600 | ! Have explicitly unrolled the loops. |
| 601 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) |
| 602 | PARAMETER (NMAX=3) |
| 603 | DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) |
| 604 | if(n .ne. 3) then |
| 605 | print *,'fatal CtLhPolint3 call',n |
| 606 | stop |
| 607 | endif |
| 608 | NS=1 |
| 609 | DIF=ABS(X-XA(1)) |
| 610 | DIFT=ABS(X-XA(1)) |
| 611 | IF (DIFT.LT.DIF) THEN |
| 612 | NS=1 |
| 613 | DIF=DIFT |
| 614 | ENDIF |
| 615 | C(1)=YA(1) |
| 616 | D(1)=YA(1) |
| 617 | DIFT=ABS(X-XA(2)) |
| 618 | IF (DIFT.LT.DIF) THEN |
| 619 | NS=2 |
| 620 | DIF=DIFT |
| 621 | ENDIF |
| 622 | C(2)=YA(2) |
| 623 | D(2)=YA(2) |
| 624 | DIFT=ABS(X-XA(3)) |
| 625 | IF (DIFT.LT.DIF) THEN |
| 626 | NS=3 |
| 627 | DIF=DIFT |
| 628 | ENDIF |
| 629 | C(3)=YA(3) |
| 630 | D(3)=YA(3) |
| 631 | Y=YA(NS) |
| 632 | NS=NS-1 |
| 633 | HO=XA(1)-X |
| 634 | HP=XA(2)-X |
| 635 | W=C(2)-D(1) |
| 636 | DEN=W/(HO-HP) |
| 637 | D(1)=HP*DEN |
| 638 | C(1)=HO*DEN |
| 639 | HO=XA(2)-X |
| 640 | HP=XA(3)-X |
| 641 | W=C(3)-D(2) |
| 642 | DEN=W/(HO-HP) |
| 643 | D(2)=HP*DEN |
| 644 | C(2)=HO*DEN |
| 645 | IF (2*NS.LT.2)THEN |
| 646 | DY=C(NS+1) |
| 647 | ELSE |
| 648 | DY=D(NS) |
| 649 | NS=NS-1 |
| 650 | ENDIF |
| 651 | Y=Y+DY |
| 652 | HO=XA(1)-X |
| 653 | HP=XA(3)-X |
| 654 | W=C(2)-D(1) |
| 655 | DEN=W/(HO-HP) |
| 656 | D(1)=HP*DEN |
| 657 | C(1)=HO*DEN |
| 658 | IF (2*NS.LT.1)THEN |
| 659 | DY=C(NS+1) |
| 660 | ELSE |
| 661 | DY=D(NS) |
| 662 | NS=NS-1 |
| 663 | ENDIF |
| 664 | Y=Y+DY |
| 665 | RETURN |
| 666 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame