]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TUHKMgen/PYQUEN/progs_fortran.f
Corrected warnings, coding conventions, bug fixes (Ionut)
[u/mrichter/AliRoot.git] / TUHKMgen / PYQUEN / progs_fortran.f
CommitLineData
b1c2e580 1c*******************************************
2
3 SUBROUTINE myini
4
5 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
6 double precision nbcol,npart,npart0
7 real p,v,plu
8 external ludata,pydata
9 double precision numpar,npar0,nbco0
7b7936e9 10c common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
b1c2e580 11 common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
12 common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
13 common /hyflow/ ytfl,ylfl,Tf,fpart
14 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
15 common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
16 common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
17 common /pypars/ mstp(200),parp(200),msti(200),pari(200)
18 common /pyqpar/ T0,tau0,nf,ienglu,ianglu
19 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
20 common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
21
22 COMMON/PYDATR/MRPY(6),RRPY(100)
23 common/SERVICE/iseed_fromC,iPythDecay,parPYTH(100)
24
25
26 common /hypart/ppart(10,50000),bmin,bmax,njp
7b7936e9 27c save /lujets/,/hyjets/,/hyflow/,/hyjpar/,/hyfpar/,/pyqpar/,
28c > /pysubs/,/pypars/,/pydat1/
b1c2e580 29
30* ----------- INITIALIZATION OF RUNDOM GENERATOR
31 MRPY(1)=iseed_fromC
32 write(*,*)'--initialization of high-pt part-- : '
33* write(*,*)'evnt generator phase ',iseed_fromC
34
35* initialize HYINIT with the input parameters
36 write(*,*)' AW= ', AW,' energy= ', ene, ' ptmin= ',ptmin
37 write(*,*)' bminh= ',bminh,' bmaxh= ', bmaxh
38 write(*,*)' ifb= ', ifb, ' bfix= ', bfix
39 write(*,*)' ishad= ',ishad, ' nhsel= ',nhsel
40 write(*,*)' ienglu= ',ienglu, ' ianglu= ',ianglu
41 write(*,*)' T0= ',T0,' tau0= ',tau0,' nf= ', nf
42
43* set input PYTHIA parameters from SERVICE common:
44 CKIN(3)=ptmin ! minimum pt in initial hard scattering, GeV
45 MSEL=0 ! QCD-dijet production
46 MSTU(21)=1 ! avoid stopping run
47 PARU(14)=1. ! tolerance parameter to adjust fragmentation
48 MSTP(81)=1 ! pp multiple scattering off
49 MSTJ(21)=0 ! hadron decays off
50 MSTP(2)=1 !which order running alphaS
51 MSTP(33)=0 !unclusion of k factor in cross section
52 MSTP(51)=7 !structure function chosen - CTEQ5M pdf
53 MSTP(82)=4 !Defines the multi-parton model
54 PARP(67)=1 !amount of initial-state radiation
55 PARP(82)=2. !pt cutoff for multiparton interactions
56 MSTJ(11)=3 !Choice of the fragmentation function
57
58* flags to select dijet production channels (active if msel=0)
59 MSUB(11)=1 ! q+q->q+q
60 MSUB(12)=1 ! q+qbar->q+qbar
61 MSUB(13)=1 ! q+qbar->g+g
62 MSUB(28)=1 ! q+g->q+g
63 MSUB(53)=1 ! g+g->q+qbar
64 MSUB(68)=1 ! g+g->g+g
65
66* flags to select prompt photon production channels (active if msel=0)
67 MSUB(14)=1 ! q+qbar->g+gamma
68 MSUB(18)=1 ! q+qbar->gamma+gamma
69 MSUB(29)=1 ! q+g->q+gamma
70 MSUB(114)=1 ! g+g->gamma+gamma
71 MSUB(115)=1 ! g+g->g+gamma
72
73
74 write(*,*)'--initialization of PYTHIA parameters-- : '
75
76 write(*,*)'CKIN(3) MSEL MSTU(21) PARU(14)'
77 write(*,*)CKIN(3), MSEL, MSTU(21), PARU(14)
78
79 write(*,*)'MSTP(81) MSTJ(21) MSTP(2) MSTP(33)'
80 write(*,*)MSTP(81), MSTJ(21), MSTP(2), MSTP(33)
81
82 write(*,*)'MSTP(51) MSTP(82) PARP(67) PARP(82) MSTJ(11)'
83 write(*,*)MSTP(51),MSTP(82),PARP(67),PARP(82),MSTJ(11)
84
85 write(*,*)'MSUB(11) MSUB(12) MSUB(13) MSUB(28) MSUB(53) MSUB(68)'
86 write(*,*)MSUB(11),MSUB(12),MSUB(13),MSUB(28),MSUB(53),MSUB(68)
87
88 write(*,*)'MSUB(14) MSUB(18) MSUB(29) MSUB(114) MSUB(115)'
89 write(*,*)MSUB(14),MSUB(18),MSUB(29),MSUB(114),MSUB(115)
90
91 call hyinit ! (ene,AW,ifb,bmin,bmax,bfix1,nh)
92
93 end
94****************end myini***************************************************************
95
96
97
98********************************* HYINIT ***************************
99 SUBROUTINE HYINIT
100
101*ml-- (energy,A,ifb1,bmin,bmax,bfix1,nh1)
102* PYTHIA inizialization, calculation of total and hard cross sections and
103* # of participants and binary sub-collisions at "reference point" (Pb,b=0),
104* tabulation of nuclear thickness function and nuclear overlap function,
105* test of input parametes
106 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
107 double precision numpar,npar0,nbco0
108 external numpar,rhoaa,hfunc3
109 common /pyint7/ sigt(0:6,0:6,0:5)
110 common /pypars/ mstp(200),parp(200),msti(200),pari(200)
111 common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
112 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
113 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
114 common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
115 common /hyflow/ ytfl,ylfl,Tf,fpart
116 common /hygeom/ BC
117 common /hythic/ BAB(110),TAB(110),TAAB(110)
118 common /hynup1/ bp,x
119
120*
121
7b7936e9 122c-- save /pyint7/,/pypars/,/pysubs/,/hyjpar/,/hyipar/,/hypyin/,
123c-- > /hyflow/,/hygeom/,/hythic/
b1c2e580 124
125* start HYDJET initialization
126 init=1
127 ipr=1
128
129* set beam paramters
130c-ml ene=energy ! c.m.s. energy per nucleon
131c-ml AW=A ! atomic weight
132
133
134 RA=1.15d0*AW**0.333333d0 ! nuclear radius
135 rzta=1.d0/(1.98d0+0.015d0*A**0.666667d0) ! fraction of protons in nucleus
136 rnta=1.d0-rzta ! fraction of neutrons in nucleus
137c-ml ifb=ifb1 ! centrality flag
138c-ml bfix=bfix1 ! fixed impact parameter
139c-ml nh=nh1 ! mean soft mult. in central PbPb
140
141c-ml ptmin=ckin(3) ! minimum pt of hard scattering
142
143c* printout of HYDJET versioning
144c write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145c >%%%%%%%%%%%%%%%%%%%%%%%%'
146c write(6,*) '%% This is HYDJET version 1.5
147c > %%'
148c write(6,*) '%% Corresponding author: Igor Lokhtin (Igor.Lo
149c >khtin@cern.ch) %%'
150c write(6,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151c >%%%%%%%%%%%%%%%%%%%%%%%%'
152
153* Pythia inizialization for pp collisions
154 call pyinit('cms','p','p',ene)
155
156* calculation of hard scattering cross section
157 if(nhsel.ne.0) then
158
159 mstp(111)=0
160* no printout of Pythia initialization information hereinafter
161 mstp(122)=0
162
163
164* initialize HYINIT with the input parameters
165 write(*,*)' in hyinit AW= ', AW,' energy= ', ene, ' ptmin= ',ptmin
166 write(*,*)' bminh= ',bminh,' bmaxh= ', bmaxh
167 write(*,*)' ifb= ', ifb, ' bfix= ', bfix
168
169
170* Pythia test pp event run
171 do i=1,1000
172 call pyevnt
173 end do
174* hard scattering pp cross section
175 sjpp=pari(1)
176
177* Pythia inizialization for pn collisions
178 call pyinit('cms','p','n',ene)
179* Pythia test pn event run
180 do i=1,1000
181 call pyevnt
182 end do
183* hard scattering pn cross section
184 sjpn=pari(1)
185
186* Pythia inizialization for nn collisions
187 call pyinit('cms','n','n',ene)
188* Pythia test nn event run
189 do i=1,1000
190 call pyevnt
191 end do
192* hard scattering nn cross section
193 sjnn=pari(1)
194
195 sigjet=rzta*rzta*sjpp+rnta*rnta*sjnn+2.d0*rzta*rnta*sjpn
196
197 end if
198
199* total inelastic cross section
200 if(sigin.lt.10.d0.or.sigin.gt.200.d0)
201 > sigin=sigt(0,0,0)-sigt(0,0,1)
202
203* # of nucelons-participants and NN sub-collisions at "reference point" (Pb,b=0)
204 Apb=207.d0
205 Rpb=1.15d0*Apb**0.333333d0
206 EPS=0.005d0
207 Z2=4.d0*Rpb
208 Z1=-1.d0*Z2
209 H=0.01d0*(Z2-Z1)
210 do ib=1,110
211 BC=3.d0*Rpb*(ib-1)/109.d0
212 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS)
213 BAB(ib)=BC
214 TAB(ib)=Apb*RES
215 end do
216 Z1=0.d0
217 Z2=6.28318d0
218 H=0.01d0*(Z2-Z1)
219 bp=0.d0
220 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,TAAPB0,AIH,AIABS)
221
222 npar0=numpar(0.d0) ! Npart(Pb,b=0)
223 nbco0=0.1d0*sigin*TAAPB0 ! Nsub(Pb,b=0)
224
225 init=0
226
227* creation of arrays for tabulation of beam/target nuclear thickness function
228 Z2=4.d0*RA
229 Z1=-1.d0*Z2
230 H=0.01d0*(Z2-Z1)
231 do ib=1,110
232 BC=3.d0*RA*(ib-1)/109.d0
233 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,rhoaa,Z,RES,AIH,AIABS)
234 BAB(ib)=BC
235 TAB(ib)=AW*RES
236 end do
237
238* creation of arrays for tabulation of nuclear overlap function
239 Z1=0.d0
240 Z2=6.28318d0
241 H=0.01d0*(Z2-Z1)
242 do ib=1,110
243 bp=BAB(ib)
244 CALL SIMPA(Z1,Z2,H,EPS,1.d-8,HFUNC3,X,RES,AIH,AIABS)
245 TAAB(ib)=RES
246 end do
247
248* test of centrality selection
249 if(ifb.eq.0) then
250 if(bfix.lt.0.d0) then
251 write(6,*) 'Impact parameter less than zero!'
252 bfix=0.d0
253 end if
254 if (bfix.gt.3.d0) then
255 write(6,*) 'Impact parameter larger than three nuclear radius!'
256 bfix=3.d0
257 end if
258 else
259 if(bmin.lt.0.d0) then
260 write(6,*) 'Impact parameter less than zero!'
261 bmin=0.d0
262 end if
263 if(bmax.gt.3.d0) then
264 write(6,*) 'Impact parameter larger than three nuclear radius!'
265 bmax=3.d0
266 end if
267c-ml bminh=bmin
268c-ml bmaxh=bmax
269 end if
270
271* test of flow parameter selection
272 if (Tf.lt.0.08d0.or.Tf.gt.0.2d0) Tf=0.1d0 ! freeze-out temperature
273 if (ylfl.lt.0.01d0.or.ylfl.gt.7.d0) ylfl=4.d0 ! longitudinal flow rapidity
274 if (ytfl.lt.0.01d0.or.ytfl.gt.3.d0) ytfl=1.5d0 ! transverse flow rapidity
275 if (fpart.le.0.d0.or.fpart.gt.1.d0) fpart=1.d0 ! fraction of soft multiplicity
276 ! proport. to # of participants
277* test of 'nhsel' selection
278 if(nhsel.ne.1.and.nhsel.ne.2.and.nhsel.ne.3.and.nhsel.ne.4)
279 > nhsel=0
280
281 return
282 end
283****************************** END HYINIT ***************************
284 SUBROUTINE MYDELTA
285 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
286 common/SERVICEEV/delta,KC,ipdg
287 COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
288 integer KC, PYCOMP
289 real delta
290* allows to transmit delta -limits for BW from PYTHIA to c++ part
291 KC=PYCOMP(ipdg)
292 delta=PMAS(int(KC),3)
293c write(*,*)" ipdg ", ipdg, " KC ",KC," delta ",delta
294 return
295 end
296
297
298********************************* HYEVNT ****************************
299 SUBROUTINE HYEVNT
300* generation of single HYDJET event (soft+hard parts) at given parameters
301 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
302 double precision numpar,npar0,nbco0,npart,nbcol,npart0
303 real p,v, delta
304 external hsin,hftaa,numpar,hyhard,hipsear,pyr,pymass,PYCOMP
305 external ludata
7b7936e9 306c common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
b1c2e580 307 common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
308 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
309 common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
310 common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
311 common /hyflow/ ytfl,ylfl,Tf,fpart
312 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
313c-ml
314 common /hypart/ppart(10,50000),bmin,bmax,njp
315 COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
316
7b7936e9 317c-- save/lujets/,/hyjets/,/hyipar/,/hyfpar/,/hyflow/,/hyjpar/,/hypyin/
b1c2e580 318
3fa37a65 319* reset hyjets arrays before event generation
b1c2e580 320
321 n=0
322 nhj=0
7b7936e9 323
b1c2e580 324 do ncl=1,150000
325 do j=1,5
7b7936e9 326c p(ncl,j)=0.
b1c2e580 327 phj(ncl,j)=0.d0
7b7936e9 328c v(ncl,j)=0.
b1c2e580 329 vhj(ncl,j)=0.d0
7b7936e9 330c k(ncl,j)=0
b1c2e580 331 khj(ncl,j)=0
332 enddo
333 end do
334*
335 pi=3.14159d0
336
337* generate impact parameter of A-A collision
338 if(ifb.eq.0) then
339 b1=bfix*RA
340 bgen=bfix
341 else
342 call hipsear(fmax1,xmin1)
343 fmax=fmax1
344 xmin=xmin1
345 3 bb1=xmin*pyr(0)+bminh*RA
346 ff1=fmax*pyr(0)
347 fb=hsin(bb1)
348 if(ff1.gt.fb) goto 3
349 b1=bb1
350 bgen=bb1/RA
351 end if
352
7b7936e9 353 write(*,*)'in hyevnt bminh bgen sigin',RA,bminh,bgen,sigin
b1c2e580 354
355* calculate # of nucelons-participants and binary NN sub-collisions
356 npart=numpar(b1) ! Npart(b)
357 npart0=numpar(0.d0) ! Npart(b)
358
359c write(*,*)'in hyevnt npart',npart
360 nbcol=0.1d0*sigin*hftaa(b1) ! Nsub(b)
361c write(*,*)'in hyevnt nbcol',nbcol
362* generate hard parton-parton scatterings (Q>ptmin) 'njet' times
363 njet=0
364 if(nhsel.ne.0) then
365 pjet=sigjet/sigin
366c write(*,*)'in hyevnt pjet',pjet
367 do i=1,int(nbcol)
368 if(pyr(0).lt.pjet) njet=njet+1
369 end do
370c write(*,*)'before hyhard'
371 call hyhard
372 end if
373
374 npyt=nhj
375
376c write(*,*)'in hyevnt pjet njet=', pjet,njet
377c-ml
378* fill array 'ppart'
379
380 do ih=1,nhj
381 if(i.le.50000)then
382 ppart(1,ih)= khj(ih,1) ! pdg ?
383 ppart(2,ih)= khj(ih,2) ! pdg -mother
384 ppart(3,ih)=phj(ih,1) ! px
385 ppart(4,ih)=phj(ih,2) ! py
386 ppart(5,ih)=phj(ih,3) ! pz
387 ppart(6,ih)=phj(ih,4) ! E
388 ppart(7,ih)= vhj(ih,1)
389 ppart(8,ih)= vhj(ih,2)
390 ppart(9,ih)= vhj(ih,3)
391 ppart(10,ih)= vhj(ih,4)
392
393 endif
394 end do
395
396 njp=nhj ! fill number of jet particles for InitialStateBjorken
397c--
398
399* fill array 'lujets'
400
7b7936e9 401c do ih=1,nhj
402c do jh=1,5
403c p(ih,jh)=phj(ih,jh)
404c k(ih,jh)=khj(ih,jh)
405c v(ih,jh)=vhj(ih,jh)
406c end do
407c end do
b1c2e580 408
409c-ml
3fa37a65 410c call luedit(2) ! remove unstable particles and partons
b1c2e580 411
412
413
414 return
415 end
416****************************** END HYEVNT **************************
417********************************* HYHARD ***************************
418 SUBROUTINE HYHARD
419* generate 'njet' number of hard parton-parton scatterings with PYTHIA
420 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
421 double precision npar0,nbco0,npart,nbcol,npart0
422 double precision gaus
423
424 INTEGER PYK,PYCHGE,PYCOMP
425 CHARACTER beam*2,targ*2
426 external pydata
427 external pyp,pyr,pyk,pyquen,shad1,gaus
428
429 common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
430 common /hyjets/ nhj,nhp,khj(150000,5),phj(150000,5),vhj(150000,5)
431 COMMON /PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
432 COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
433 COMMON /PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
434 COMMON /PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
435 common /pypars/ mstp(200),parp(200),msti(200),pari(200)
436 common /parimp/ b1, psib1, r0, rb1, rb2, noquen
437 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
438 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
439 common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
440 common /hypyin/ ene,rzta,rnta,bfix,ifb,nh
7b7936e9 441c-- save /pyjets/,/pypars/,/pydat1/,/pydat2/,/pydat3/,/pysubs/,
442c-- + /hyjets/,/parimp/,/hyjpar/,/hyipar/,/hyfpar/,/hypyin/
b1c2e580 443
444* generate 'njet' PYTHIA events and fill arrays for partons and hadrons
445
446c write(*,*)'in hyhard 0'
447
448 nshad=0
449 noquen=0
450 if(nhsel.eq.1.or.nhsel.eq.3) noquen=1
451 if(njet.ge.1) then
452 mdcy(pycomp(111),1)=0 ! no pi0 decay
453 mdcy(pycomp(310),1)=0 ! no K_S0 decay
454
455 ifbp=0 ! fix impact parameter
456 bfixp=RA*bgen
457 Ap=AW ! atomic weight
458
459 do ihard=1,njet
460 mstp(111)=0
461
462* generate type of NN sub-collision (pp, pn or nn)
463 rand1=pyr(0)
464 if(rand1.lt.rzta) then
465 beam='p'
466 else
467 beam='n'
468 end if
469 rand2=pyr(0)
470 if(rand2.lt.rzta) then
471 targ='p'
472 else
473 targ='n'
474 end if
475 call pyinit('cms',beam,targ,ene)
476c mstj(41)=0 ! vacuum showering off
477
478
479c write(*,*)'in hyhard 1'
480
481 call pyevnt ! generate hard scattering
482
483c write(*,*)'in hyhard 2'
484
485
486* PYQUEN: quenched jets if noquen=0 or non-quenched jets if noquen=1 and ishad=1
487 if(ishad.eq.1.or.nhsel.eq.2.or.nhsel.eq.4)
488 > call pyquen(Ap,ifbp,bfixp)
489
490c-ml coordinate info if we need of it
491
492c write(*,*)'in hyhard pari(21)',pari(21)
493
494 Q=pari(21)
495 x=r0*cos(psib1) !fm
496 y=r0*sin(psib1) !fm
497c write(*,*)'in hyhard ',x,y
498
499 if(Q.gt.0)then
500 tau=1./Q ! 1/GeV
501 tau=tau*0.197 !fm/c
502 rm=0.0d0
503 sig=1.0d0
504
505c write(*,*)'in hyhard tau',tau
506
507 etaLj=gaus(rm,sig) !fm
508
509c write(*,*)'in hyhard after gauss',tau
510 z=tau*sinh(etaLj)
511 t=tau*cosh(etaLj)
512c write(*,*)'in hyhard z t',z,t
513 endif
514
515c write(*,*)'x y z t', x,y,z,t,etaLj,Q
516
517* treatment of "nuclear shadowing" (for Pb, Au, Pd or Ca beams only)
518
519
520 if(ishad.eq.1) then
521 kfh1=abs(k(3,2))
522 kfh2=abs(k(4,2))
523 xh1=pari(33)
524 xh2=pari(34)
525 Q2=pari(22)
526 shad=shad1(kfh1,xh1,Q2,rb1)*shad1(kfh2,xh2,Q2,rb2)
527 if(pyr(0).gt.shad) then
528 nshad=nshad+1
529 goto 53
530 end if
531 end if
532
533
534c write(*,*)'in hyhard 4'
535
536 call pyexec ! hadronization done
537
538c write(*,*)'in hyhard 5'
539
540c-ml
541 call pyedit(2) ! remove partons & leave hadrons
542
543c write(*,*)'in hyhard 6'
544
545
546* fill array of final particles
547 nu=nhj+n
548 if(nu.gt.150000-np) then
549 write(6,*) 'Warning, multiplicity too large! Cut hard part.'
550 goto 52
551 end if
552 nhj=nu
553
554
555 do i=nhj-n+1,nhj
556 ip=i+n-nhj
557c write(*,*)'in hyhard 7'
558
559 do j=1,5
560 phj(i,j)=p(ip,j)
561 end do
562
563 vhj(i,1)=x
564 vhj(i,2)=y
565 vhj(i,3)=z
566 vhj(i,4)=t
567
568 do j=1,5
569 khj(i,j)=k(ip,j)
570 end do
571 do j=3,5
572 kk=khj(i,j)
573 if(kk.gt.0) then
574 khj(i,j)=kk+nhj-n
575 end if
576 end do
577 end do
578
579
580 53 continue
581
582c write(*,*)'in hyhard 8'
583
584 end do
585 52 njet=ihard-1
586
587c write(*,*)'in hyhard 9', njet,nshad
588
589 end if
590 njet=njet-nshad
591
592
593c write(*,*)'in hyhard OUT'
594
595
596 return
597 end
598****************************** END HYHARD **************************
599********************************* HIPSEAR ***************************
600 SUBROUTINE HIPSEAR (fmax,xmin)
601* find maximum and 'sufficient minimum' of differential inelasic AA cross
602* section as a function of impact paramater (xm, fm are outputs)
603 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
604 double precision npar0,nbco0
605 external hsin
606 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
7b7936e9 607c-- save /hyipar/
b1c2e580 608 xmin=(bmaxh-bminh)*RA
609
610c write(*,*)'bmaxh bminh',bmaxh,bminh
611
612 fmax=0.d0
613 do j=1,1000
614 x=bminh*RA+xmin*(j-1)/999.d0
615 f=hsin(x)
616 if(f.gt.fmax) then
617 fmax=f
618 endif
619 end do
620 return
621 end
622****************************** END HIPSEAR **************************
623
624* differential inelastic AA cross section
625 double precision function hsin(x)
626 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
627 external hftaa
628 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
7b7936e9 629c-- save /hyjpar/
b1c2e580 630 br=x
631 hsin=br*(1.d0-dexp(-0.1d0*hftaa(br)*sigin))
632 return
633 end
634
635* number of nucleons-participants at impact parameter b
636 double precision function numpar(c)
637 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
638 external HFUNC1
639 common /hynup1/ bp,x
640 EPS=0.005d0
641 A=0.d0
642 B=6.28318d0
643 H=0.01d0*(B-A)
644 bp=c
645 CALL SIMPA(A,B,H,EPS,1.d-8,HFUNC1,X,RES,AIH,AIABS)
646 numpar=2.d0*RES
647 return
648 end
649*
650 double precision function HFUNC1(x)
651 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
652 double precision npar0,nbco0
653 external HFUNC2
654 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
655 common /hynup1/ bp,xx
7b7936e9 656c-- save /hyipar/
b1c2e580 657 if(init.eq.1) then
658 Rl=Rpb
659 else
660 Rl=RA
661 end if
662 xx=x
663 EPS=0.005d0
664 A=0.d0
665 B=3.d0*Rl
666 H=0.01d0*(B-A)
667 CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC2,Y,RES,AIH,AIABS)
668 HFUNC1=RES
669 return
670 end
671*
672 double precision function HFUNC2(y)
673 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
674 double precision npar0,nbco0
675 external hythik
676 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
677 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
678 common /hynup1/ bp,x
7b7936e9 679c-- save /hyipar/,/hyjpar/
b1c2e580 680 r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x)))
681 r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x)))
682 s=1.d0-dexp(-0.1d0*sigin*hythik(r2))
683 HFUNC2=y*hythik(r1)*s
684 return
685 end
686
687* nuclear overlap function at impact parameter b
688 double precision function hftaa(c)
689 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
690 common /hythic/ BAB(110),TAB(110),TAAB(110)
7b7936e9 691c-- save /hythic/
b1c2e580 692 call parinv(c,BAB,TAAB,110,RES)
693 hftaa=RES
694 return
695 end
696*
697 double precision function HFUNC3(x)
698 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
699 double precision npar0,nbco0
700 external HFUNC4
701 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
702 common /hynup1/ bp,xx
7b7936e9 703c-- save /hyipar/
b1c2e580 704 if(init.eq.1) then
705 Rl=Rpb
706 else
707 Rl=RA
708 end if
709 xx=x
710 EPS=0.005d0
711 A=0.d0
712 B=3.d0*Rl
713 H=0.01d0*(B-A)
714 CALL SIMPB(A,B,H,EPS,1.d-8,HFUNC4,Y,RES,AIH,AIABS)
715 HFUNC3=RES
716 return
717 end
718*
719 double precision function HFUNC4(y)
720 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
721 double precision npar0,nbco0
722 external hythik
723 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
724 common /hyjpar/ ptmin,sigin,sigjet,nhsel,ishad,njet
725 common /hynup1/ bp,x
7b7936e9 726c-- save /hyipar/,/hyjpar/
b1c2e580 727 r1=dsqrt(abs(y*y+bp*bp/4.d0+y*bp*dcos(x)))
728 r2=dsqrt(abs(y*y+bp*bp/4.d0-y*bp*dcos(x)))
729 HFUNC4=y*hythik(r1)*hythik(r2)
730 return
731 end
732
733* nuclear thickness function
734 double precision function hythik(r)
735 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
736 common /hythic/ BAB(110),TAB(110),TAAB(110)
7b7936e9 737c-- save /hythic/
b1c2e580 738 call parinv(r,BAB,TAB,110,RES)
739 hythik=RES
740 return
741 end
742
743* Wood-Saxon nucleon distrubution
744 double precision function rhoaa(z)
745 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
746 double precision npar0,nbco0
747 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
748 common /hygeom/ BC
7b7936e9 749c-- save /hyipar/,/hygeom/
b1c2e580 750 if(init.eq.1) then
751 Rl=Rpb
752 else
753 Rl=RA
754 end if
755 pi=3.14159d0
756 df=0.54d0
757 r=sqrt(bc*bc+z*z)
758 rho0=3.d0/(4.d0*pi*Rl**3)/(1.d0+(pi*df/Rl)**2)
759 rhoaa=rho0/(1.d0+dexp((r-Rl)/df))
760 return
761 end
762
763* function to calculate nuclear shadowing factor (for Pb, Au, Pd or Ca beams)
764 double precision function shad1(kfh,xbjh,Q2h,r)
765 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
766 double precision npar0,nbco0,nbcol,npart,npart0
767 external ggshad
768 common /hyipar/ bminh,bmaxh,AW,RA,npar0,nbco0,Apb,Rpb,np,init,ipr
769 common /hyfpar/ bgen,nbcol,npart,npart0,npyt,nhyd
770 common /hyshad/ bbmin,bbmax,inuc
7b7936e9 771c-- save /hyipar/,/hyfpar/,/hyshad/
b1c2e580 772 dimension res(2)
773 kf=kfh
774 xbj=xbjh
775 Q2=Q2h
776 bb=r
777 inuc=0
778 if(AW.gt.205.d0.and.AW.lt.209.d0) inuc=4 ! Pb-206, 207 or 208
779 if(AW.gt.196.d0.and.AW.lt.198.d0) inuc=3 ! Au-197
780 if(AW.gt.109.d0.and.AW.lt.111.d0) inuc=2 ! Pd-110
781 if(AW.gt.39.d0.and.AW.lt.41.d0) inuc=1 ! Ca-40
782 if(inuc.eq.0.and.ipr.eq.1) then
783 write(6,*)
784 > 'Warning! Shadowing is not foreseen for atomic weigth A ='
785 > ,AW
786 write(6,*)'******************************************************
787 >************************'
788 ipr=0
789 end if
790 if(inuc.ne.0) then
791 xbj=max(5.d-5,xbj)
792 xbj=min(0.95d0,xbj)
793 Q2=max(4.d0,Q2)
794 Q2=min(520.d0,Q2)
795 bb=max(0.d0,bb)
796 call ggshad(inuc,xbj,Q2,bb,res,ta)
797 if(kf.eq.21) then
798 shad1=res(1)
799 elseif(kf.eq.1.or.kf.eq.2.or.kf.eq.3) then
800 shad1=res(2)
801 else
802 shad1=1.d0
803 end if
804 else
805 shad1=1.d0
806 end if
807 return
808 end
809
810******************************************************************************
811* The part of the code which follows below, includes nuclear shadowing model *
812* and has been written by Konrad Tywoniuk (Oslo University, Norway) *
813******************************************************************************
814c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$
815c$$$
816c$$$
817c$$$ Shadowing from Glauber-Gribov theory for
818c$$$ gluons and light quarks (u,d,s and their antiquarks).
819c$$$ All Pomeron tree diagrams have been summed (Schwimmer model).
820c$$$ More details about the model in
821c$$$ K. Tywoniuk, I.C. Arsene, L. Bravina, A. Kaidalov and E.
822c$$$ Zabrodin, Phys. Lett. B 657 (2007) 170
823c$$$
824c$$$ We use FIT B from the H1 parameterization published in
825c$$$ A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 715
826c$$$ A. Aktas et al. (H1 Collaboration), Eur. Phys. J C 48 (2006) 749
827c$$$
828c$$$ Main routine is GGSHAD which provides the user with the
829c$$$ ratio of nuclear and nucleon parton distribution function
830c$$$ normalized by the atomic number for a given
831c$$$ INUCL: nucleus (1=Ca,2=Pd,3=Au,4=Pb)
832c$$$ X: Bjorken x
833c$$$ Q2: q**2, scale squared
834c$$$ B: impact parameter
835c$$$
836c$$$ Then RES(1): gluon shadowing
837c$$$ RES(2): sea quark shadowing
838c$$$ TA: nuclear profile function
839c$$$
840c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$$$c$
841 SUBROUTINE GGSHAD(INUCL,X,Q2,B,RES,TAF)
842 IMPLICIT NONE
843 DOUBLE PRECISION TA(31,4),IMPAR(31),ANUCL(4)
844 DOUBLE PRECISION XB(36),Q2V(13)
845 DOUBLE PRECISION XMAX,XMAXX,Q2MIN,Q2MAX,BMAX
846 DOUBLE PRECISION G(36,13,4),LQ(36,13,4)
847 DOUBLE PRECISION TATMP(31),SHAD(2)
848 DOUBLE PRECISION C(100),D(100),E(100)
849 DOUBLE PRECISION X,Q2,B
850 DOUBLE PRECISION RES(2)
851 DOUBLE PRECISION TAF
852 DOUBLE PRECISION SEVAL
853 INTEGER INUCL,IREAD,TMAX
854 INTEGER I,J,K,IK
855
856 PARAMETER(TMAX=31)
857 PARAMETER(XMAX=0.1)
858 PARAMETER(XMAXX=0.95)
859 PARAMETER(Q2MIN=4.)
860 PARAMETER(Q2MAX=520.)
861 PARAMETER(BMAX=30.)
862
863 COMMON/GG07/XB,Q2V,G,LQ
864
865 DATA ANUCL/40.,110.,197.,206./
866 DATA IMPAR
867 > /0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,6.5,7.,7.5,8.,
868 > 8.5,9.,9.5,10.,10.5,11.,11.5,12.,12.5,13.,13.5,14.,14.5,15./
869 DATA TA
870 > /0.0291662,0.0288636,0.0279354,0.0263044,0.0238341,0.0203565,
871 > 0.0158333,0.0107343,0.00617758,0.00307859,0.00139738,
872 > 0.000603839,0.000254843,0.000106329,4.40965e-05,1.82206e-05,
873 > 7.50949e-06,3.08881e-06,1.26837e-06,5.20082e-07,2.12978e-07,
874 > 8.71154e-08,3.55956e-08,1.45304e-08,5.92622e-09,2.41504e-09,
875 > 9.83429e-10,4.00184e-10,1.62741e-10,6.61407e-11,2.68657e-11,
876 > 0.0153487,0.0152775,0.0150614,0.0146926,0.0141561,0.0134268,
877 > 0.0124651,0.0112115,0.00959258,0.00756907,0.00527514,
878 > 0.00313482,0.00159852,0.000732479,0.000316628,0.000133116,
879 > 5.52508e-05,2.27904e-05,9.3692e-06,3.84349e-06,1.57423e-06,
880 > 6.43962e-07,2.63132e-07,1.07414e-07,4.38089e-08,1.78529e-08,
881 > 7.2699e-09,2.95832e-09,1.20304e-09,4.8894e-10,1.98603e-10,
882 > 0.0105299,0.010498,0.0104017,0.010239,0.010006,0.00969677,
883 > 0.00930205,0.00880781,0.00819264,0.00742455,0.00646058,
884 > 0.00526252,0.00385835,0.00243988,0.00131252,0.000621079,
885 > 0.000272322,0.000114993,4.77312e-05,1.96569e-05,8.06377e-06,
886 > 3.30066e-06,1.34904e-06,5.50751e-07,2.24633e-07,9.15434e-08,
887 > 3.72777e-08,1.51694e-08,6.16885e-09,2.50714e-09,1.01838e-09,
888 > 0.0102179,0.0101879,0.0100975,0.00994467,0.00972606,
889 > 0.00943628,0.0090671,0.00860605,0.00803416,0.00732289,
890 > 0.00643266,0.00532296,0.00400025,0.00261264,0.00144993,
891 > 0.00070098,0.000310867,0.000131949,5.48887e-05,2.26247e-05,
892 > 9.28457e-06,3.80091e-06,1.55358e-06,6.34273e-07,2.58701e-07,
893 > 1.05427e-07,4.29316e-08,1.74701e-08,7.10447e-09,2.88739e-09,
894 > 1.17283e-09/
895 DATA Q2V
896 > /4.000, 6.000, 9.000, 13.500, 20.250, 30.375, 45.562,
897 > 68.344, 102.516, 153.773, 230.660, 345.990, 518.985/
898 DATA XB
899 > /0.99999998E-05, 0.13000000E-04, 0.16899999E-04,
900 > 0.21970000E-04, 0.28561000E-04, 0.37129299E-04,
901 > 0.48268099E-04, 0.62748499E-04, 0.81573096E-04,
902 > 0.10604500E-03, 0.13785800E-03, 0.17921600E-03,
903 > 0.23298099E-03, 0.30287501E-03, 0.39373801E-03,
904 > 0.51185902E-03, 0.66541700E-03, 0.86504198E-03,
905 > 0.11245500E-02, 0.14619200E-02, 0.19005000E-02,
906 > 0.24706500E-02, 0.32118401E-02, 0.41753901E-02,
907 > 0.54280101E-02, 0.70564100E-02, 0.91733299E-02,
908 > 0.11925300E-01, 0.15502900E-01, 0.20153800E-01,
909 > 0.26200000E-01, 0.34059901E-01, 0.44277899E-01,
910 > 0.57561301E-01, 0.74829698E-01, 0.97278602E-01/
911 DATA G
912 > /2.542140,2.457280,2.372200,2.283920,2.195320,2.105540,
913 > 2.014040,1.921480,1.827600,1.733930,1.639990,1.546170,
914 > 1.453220,1.361040,1.268750,1.178500,1.088730,1.000740,
915 > 0.914606,0.830842,0.749315,0.671419,0.596883,0.526173,
916 > 0.460067,0.398489,0.341895,0.289998,0.242628,0.199405,
917 > 0.159625,0.122521,0.087252,0.053699,0.022804,0.000414,
918 > 1.825600,1.768440,1.712550,1.654660,1.595830,1.533980,
919 > 1.472860,1.409910,1.346060,1.282030,1.217520,1.153050,
920 > 1.089030,1.025930,0.962709,0.900094,0.837787,0.776733,
921 > 0.716072,0.656570,0.598295,0.541607,0.486661,0.433846,
922 > 0.383421,0.335585,0.290685,0.248542,0.209305,0.172614,
923 > 0.138166,0.105552,0.074246,0.044592,0.017935,0.000291,
924 > 1.445600,1.405920,1.364150,1.322050,1.276900,1.232020,
925 > 1.184670,1.136620,1.086950,1.037370,0.987225,0.937488,
926 > 0.888244,0.839061,0.790370,0.741798,0.693426,0.645511,
927 > 0.598298,0.551516,0.505500,0.460400,0.416326,0.373488,
928 > 0.332281,0.292614,0.254859,0.218944,0.184976,0.152759,
929 > 0.122077,0.092747,0.064477,0.037873,0.014527,0.000216,
930 > 1.209700,1.179300,1.146560,1.112700,1.077790,1.041400,
931 > 1.003150,0.964598,0.923825,0.882423,0.841304,0.800381,
932 > 0.759486,0.718947,0.678587,0.638579,0.598671,0.558978,
933 > 0.519524,0.480617,0.442116,0.404231,0.367057,0.330599,
934 > 0.295291,0.261100,0.228220,0.196605,0.166367,0.137331,
935 > 0.109496,0.082677,0.056835,0.032696,0.012013,0.000165,
936 > 1.049390,1.024810,0.998839,0.971705,0.942633,0.912223,
937 > 0.880504,0.847236,0.812716,0.777477,0.741670,0.706578,
938 > 0.671021,0.636183,0.601384,0.566923,0.532196,0.497956,
939 > 0.463884,0.430069,0.396612,0.363536,0.330915,0.298900,
940 > 0.267627,0.237155,0.207697,0.179169,0.151566,0.124961,
941 > 0.099266,0.074436,0.050544,0.028462,0.010013,0.000127,
942 > 0.943086,0.922388,0.900700,0.877793,0.853134,0.826808,
943 > 0.798937,0.769284,0.738709,0.706829,0.674818,0.643041,
944 > 0.611250,0.580168,0.548982,0.517908,0.486720,0.456021,
945 > 0.425312,0.394899,0.364616,0.334682,0.305173,0.276042,
946 > 0.247465,0.219584,0.192474,0.166033,0.140391,0.115481,
947 > 0.091347,0.068001,0.045636,0.025180,0.008510,0.000101,
948 > 0.858923,0.841883,0.823827,0.803952,0.782124,0.759447,
949 > 0.734659,0.708141,0.680281,0.651690,0.622372,0.593451,
950 > 0.564610,0.536045,0.507521,0.479248,0.450845,0.422731,
951 > 0.394687,0.366774,0.339095,0.311580,0.284350,0.257527,
952 > 0.231101,0.205195,0.179929,0.155178,0.131077,0.107555,
953 > 0.084712,0.062629,0.041557,0.022497,0.007326,0.000082,
954 > 0.790086,0.776082,0.760408,0.743282,0.724441,0.703812,
955 > 0.681673,0.657893,0.632354,0.605918,0.578976,0.552364,
956 > 0.525820,0.499665,0.473306,0.447125,0.421027,0.394926,
957 > 0.369017,0.343250,0.317596,0.292061,0.266805,0.241770,
958 > 0.217121,0.192830,0.169123,0.145809,0.122984,0.100685,
959 > 0.078967,0.057983,0.038065,0.020236,0.006364,0.000067,
960 > 0.732957,0.720881,0.707633,0.692613,0.676200,0.657914,
961 > 0.637626,0.615846,0.592373,0.567813,0.542877,0.518036,
962 > 0.493386,0.468899,0.444596,0.420291,0.395908,0.371641,
963 > 0.347423,0.323358,0.299396,0.275511,0.251830,0.228313,
964 > 0.205146,0.182274,0.159813,0.137693,0.115983,0.094703,
965 > 0.073974,0.053979,0.035068,0.018327,0.005580,0.000056,
966 > 0.684126,0.674389,0.662432,0.649473,0.634783,0.618458,
967 > 0.600173,0.580108,0.558464,0.535421,0.512051,0.488760,
968 > 0.465717,0.442837,0.419997,0.397309,0.374393,0.351598,
969 > 0.328943,0.306305,0.283685,0.261210,0.238875,0.216696,
970 > 0.194737,0.173035,0.151692,0.130607,0.109864,0.089483,
971 > 0.069628,0.050489,0.032492,0.016714,0.004940,0.000047,
972 > 0.642684,0.634012,0.624311,0.612687,0.599751,0.584885,
973 > 0.568165,0.549716,0.529352,0.507582,0.485673,0.463844,
974 > 0.442031,0.420480,0.398892,0.377446,0.355834,0.334383,
975 > 0.312903,0.291520,0.270128,0.248790,0.227632,0.206518,
976 > 0.185624,0.164939,0.144527,0.124348,0.104415,0.084843,
977 > 0.065754,0.047403,0.030217,0.015305,0.004396,0.000040,
978 > 0.606784,0.599501,0.591216,0.581053,0.569262,0.555644,
979 > 0.540298,0.522979,0.503888,0.483449,0.462725,0.441941,
980 > 0.421308,0.400965,0.380515,0.360166,0.339667,0.319306,
981 > 0.298896,0.278531,0.258224,0.237932,0.217694,0.197574,
982 > 0.177567,0.157731,0.138168,0.118766,0.099587,0.080708,
983 > 0.062310,0.044656,0.028210,0.014081,0.003935,0.000035,
984 > 0.575250,0.569055,0.561834,0.552756,0.542358,0.530173,
985 > 0.515880,0.499697,0.481592,0.462197,0.442485,0.422845,
986 > 0.403137,0.383704,0.364332,0.344943,0.325467,0.306020,
987 > 0.286496,0.267106,0.247644,0.228251,0.208881,0.189557,
988 > 0.170377,0.151316,0.132475,0.113748,0.095226,0.076992,
989 > 0.059214,0.042195,0.026424,0.013004,0.003541,0.000030,
990 > 2.510860,2.424890,2.338800,2.251430,2.161030,2.070080,
991 > 1.976500,1.884380,1.789520,1.695430,1.600900,1.507110,
992 > 1.413980,1.321430,1.229430,1.138950,1.049380,0.961914,
993 > 0.875863,0.792403,0.711728,0.634220,0.560088,0.490193,
994 > 0.424574,0.363647,0.307681,0.256400,0.209654,0.167105,
995 > 0.128273,0.092911,0.060892,0.033178,0.011804,0.000172,
996 > 1.797890,1.740890,1.683730,1.625330,1.564220,1.502940,
997 > 1.439440,1.375560,1.311540,1.247060,1.182260,1.117820,
998 > 1.053940,0.990598,0.927193,0.864982,0.802868,0.741596,
999 > 0.681232,0.621931,0.563990,0.507817,0.452982,0.400714,
1000 > 0.350809,0.303514,0.259139,0.217505,0.178849,0.142987,
1001 > 0.109689,0.079022,0.051121,0.027212,0.009208,0.000121,
1002 > 1.421860,1.380980,1.338550,1.294150,1.249120,1.201850,
1003 > 1.153940,1.105280,1.054910,1.004680,0.954677,0.905086,
1004 > 0.855141,0.806467,0.757283,0.708817,0.660443,0.613233,
1005 > 0.565935,0.519307,0.473770,0.428741,0.385002,0.342676,
1006 > 0.301797,0.262591,0.225442,0.189974,0.156683,0.125290,
1007 > 0.095921,0.068630,0.043899,0.022874,0.007411,0.000090,
1008 > 1.187690,1.155530,1.122280,1.087030,1.051020,1.013270,
1009 > 0.973852,0.933953,0.892905,0.851713,0.810455,0.769192,
1010 > 0.728379,0.687798,0.647423,0.607517,0.567475,0.528048,
1011 > 0.488850,0.450119,0.411897,0.374353,0.337466,0.301396,
1012 > 0.266522,0.232713,0.200387,0.169339,0.139725,0.111661,
1013 > 0.085284,0.060612,0.038324,0.019572,0.006096,0.000069,
1014 > 1.028240,1.002610,0.975714,0.947100,0.916605,0.885287,
1015 > 0.852214,0.817988,0.782800,0.747390,0.711939,0.676325,
1016 > 0.641102,0.606200,0.571466,0.536994,0.502546,0.468438,
1017 > 0.434547,0.400992,0.367663,0.334910,0.302578,0.270890,
1018 > 0.240102,0.210041,0.181153,0.153175,0.126406,0.100826,
1019 > 0.076605,0.054074,0.033774,0.016896,0.005055,0.000053,
1020 > 0.923239,0.901357,0.878307,0.853928,0.827994,0.800448,
1021 > 0.770833,0.740864,0.709208,0.677436,0.645505,0.613815,
1022 > 0.582211,0.551121,0.519863,0.488847,0.457919,0.427291,
1023 > 0.396728,0.366593,0.336576,0.306991,0.277689,0.248974,
1024 > 0.220914,0.193485,0.166892,0.141077,0.116297,0.092513,
1025 > 0.069960,0.049024,0.030239,0.014840,0.004277,0.000042,
1026 > 0.840109,0.821572,0.802221,0.780737,0.758007,0.733570,
1027 > 0.707391,0.680060,0.651741,0.622598,0.593525,0.564790,
1028 > 0.535923,0.507607,0.479180,0.450770,0.422703,0.394808,
1029 > 0.366922,0.339201,0.311718,0.284555,0.257727,0.231240,
1030 > 0.205291,0.179821,0.155134,0.131145,0.107877,0.085551,
1031 > 0.064423,0.044815,0.027334,0.013173,0.003668,0.000034,
1032 > 0.771898,0.756403,0.739454,0.720917,0.700462,0.678667,
1033 > 0.655360,0.630295,0.604073,0.577412,0.550894,0.524237,
1034 > 0.497613,0.471552,0.445325,0.419311,0.393354,0.367582,
1035 > 0.341816,0.316292,0.290960,0.265746,0.240859,0.216190,
1036 > 0.191974,0.168255,0.145100,0.122534,0.100715,0.079595,
1037 > 0.059666,0.041209,0.024869,0.011778,0.003176,0.000028,
1038 > 0.715432,0.701697,0.687335,0.671043,0.652774,0.633196,
1039 > 0.611618,0.588703,0.564555,0.539941,0.515046,0.490273,
1040 > 0.465666,0.441417,0.417189,0.393083,0.368842,0.344918,
1041 > 0.320853,0.296990,0.273346,0.249853,0.226554,0.203403,
1042 > 0.180727,0.158303,0.136484,0.115160,0.094448,0.074509,
1043 > 0.055551,0.038122,0.022769,0.010610,0.002775,0.000023,
1044 > 0.666960,0.655627,0.643169,0.628497,0.612317,0.594375,
1045 > 0.574670,0.553372,0.531024,0.507805,0.484742,0.461630,
1046 > 0.438595,0.415762,0.393171,0.370479,0.347803,0.325380,
1047 > 0.302957,0.280541,0.258199,0.236074,0.214141,0.192375,
1048 > 0.170880,0.149710,0.129011,0.108739,0.089084,0.069999,
1049 > 0.052020,0.035461,0.020980,0.009629,0.002450,0.000020,
1050 > 0.626089,0.616427,0.605000,0.592191,0.577470,0.561221,
1051 > 0.543098,0.523316,0.502350,0.480560,0.458749,0.436948,
1052 > 0.415413,0.393933,0.372648,0.351193,0.329853,0.308661,
1053 > 0.287323,0.266214,0.245090,0.224191,0.203371,0.182736,
1054 > 0.162324,0.142202,0.122478,0.103078,0.084237,0.066064,
1055 > 0.048902,0.033114,0.019409,0.008778,0.002175,0.000017,
1056 > 0.590664,0.582220,0.572479,0.560768,0.547507,0.532490,
1057 > 0.515657,0.497052,0.477369,0.456893,0.436199,0.415686,
1058 > 0.395192,0.374819,0.354542,0.334339,0.314120,0.293919,
1059 > 0.273897,0.253708,0.233747,0.213772,0.193977,0.174220,
1060 > 0.154737,0.135524,0.116625,0.098102,0.080014,0.062548,
1061 > 0.046083,0.031041,0.018032,0.008042,0.001942,0.000014,
1062 > 0.559453,0.552473,0.543647,0.533124,0.521070,0.507229,
1063 > 0.491452,0.474143,0.455506,0.435945,0.416256,0.396977,
1064 > 0.377326,0.358052,0.338792,0.319598,0.300262,0.281122,
1065 > 0.261883,0.242706,0.223579,0.204566,0.185577,0.166710,
1066 > 0.148088,0.129570,0.111434,0.093593,0.076187,0.059415,
1067 > 0.043600,0.029214,0.016816,0.007398,0.001744,0.000012,
1068 > 2.494430,2.407750,2.320680,2.230570,2.139410,2.045820,
1069 > 1.954300,1.859160,1.764370,1.669990,1.574490,1.481120,
1070 > 1.386830,1.294720,1.202200,1.111670,1.021960,0.934479,
1071 > 0.848525,0.765125,0.684539,0.607178,0.533502,0.463759,
1072 > 0.398700,0.338188,0.282593,0.231728,0.185473,0.143600,
1073 > 0.105860,0.072262,0.043334,0.020543,0.005838,0.000062,
1074 > 1.781830,1.724140,1.666660,1.605560,1.544230,1.481280,
1075 > 1.418440,1.353130,1.288210,1.223120,1.158330,1.093950,
1076 > 1.029490,0.965787,0.902413,0.840091,0.777784,0.716620,
1077 > 0.656180,0.597156,0.539338,0.483075,0.428748,0.376710,
1078 > 0.327154,0.280125,0.236105,0.194961,0.156820,0.121629,
1079 > 0.089495,0.060685,0.035892,0.016616,0.004508,0.000044,
1080 > 1.407480,1.365110,1.320730,1.276090,1.228930,1.181910,
1081 > 1.132930,1.082570,1.032770,0.982143,0.932084,0.882066,
1082 > 0.832465,0.783250,0.734272,0.685821,0.637505,0.589736,
1083 > 0.542652,0.496248,0.450545,0.405906,0.362451,0.320187,
1084 > 0.279650,0.240757,0.203901,0.169026,0.136216,0.105650,
1085 > 0.077490,0.052176,0.030469,0.013808,0.003599,0.000032,
1086 > 1.173640,1.140320,1.105880,1.069540,1.032320,0.993897,
1087 > 0.953913,0.913021,0.871617,0.829958,0.788517,0.747238,
1088 > 0.706258,0.666094,0.625534,0.585422,0.545512,0.506014,
1089 > 0.466972,0.428374,0.390208,0.352773,0.316057,0.280282,
1090 > 0.245648,0.212195,0.180160,0.149604,0.120645,0.093445,
1091 > 0.068267,0.045631,0.026330,0.011698,0.002940,0.000025,
1092 > 1.014900,0.987942,0.959684,0.930270,0.899044,0.866678,
1093 > 0.832642,0.797620,0.762091,0.726232,0.690987,0.655077,
1094 > 0.620091,0.585191,0.550543,0.515986,0.481432,0.447342,
1095 > 0.413584,0.380034,0.346933,0.314295,0.282265,0.250834,
1096 > 0.220232,0.190535,0.162003,0.134530,0.108400,0.083726,
1097 > 0.060878,0.040604,0.022978,0.010002,0.002422,0.000019,
1098 > 0.909560,0.887017,0.862965,0.837572,0.810211,0.781875,
1099 > 0.751746,0.720666,0.688703,0.656776,0.624847,0.593122,
1100 > 0.561398,0.530391,0.499138,0.468280,0.437229,0.406831,
1101 > 0.376391,0.346350,0.316458,0.286984,0.257979,0.229524,
1102 > 0.201720,0.174606,0.148473,0.123248,0.099111,0.076308,
1103 > 0.055184,0.036277,0.020402,0.008711,0.002038,0.000015,
1104 > 0.827323,0.808240,0.787194,0.764785,0.741102,0.715710,
1105 > 0.688724,0.660472,0.631563,0.602307,0.573344,0.544275,
1106 > 0.515643,0.487203,0.458901,0.430660,0.402505,0.374606,
1107 > 0.346855,0.319408,0.292040,0.265141,0.238553,0.212326,
1108 > 0.186666,0.161641,0.137398,0.113940,0.091458,0.070180,
1109 > 0.050481,0.032923,0.018295,0.007672,0.001739,0.000012,
1110 > 0.759389,0.742705,0.724897,0.705170,0.684041,0.661178,
1111 > 0.636646,0.610897,0.584460,0.557556,0.530715,0.504308,
1112 > 0.477779,0.451618,0.425597,0.399640,0.373702,0.348041,
1113 > 0.322356,0.297020,0.271779,0.246872,0.222152,0.197859,
1114 > 0.173964,0.150642,0.127976,0.105994,0.084903,0.064934,
1115 > 0.046463,0.030071,0.016527,0.006814,0.001499,0.000010,
1116 > 0.703019,0.688814,0.673059,0.655722,0.636535,0.615667,
1117 > 0.593245,0.569605,0.545086,0.520194,0.495411,0.470863,
1118 > 0.446201,0.421921,0.397793,0.373652,0.349528,0.325620,
1119 > 0.301845,0.278188,0.254649,0.231349,0.208318,0.185538,
1120 > 0.163153,0.141239,0.119901,0.099172,0.079281,0.060422,
1121 > 0.043027,0.027647,0.015032,0.006099,0.001305,0.000008,
1122 > 0.655549,0.643053,0.628966,0.613794,0.596440,0.577313,
1123 > 0.556691,0.534653,0.511805,0.488542,0.465448,0.442353,
1124 > 0.419418,0.396807,0.374153,0.351526,0.328993,0.306561,
1125 > 0.284259,0.262079,0.240000,0.218086,0.196361,0.174927,
1126 > 0.153790,0.133102,0.112916,0.093252,0.074414,0.056530,
1127 > 0.040069,0.025572,0.013769,0.005504,0.001148,0.000007,
1128 > 0.614886,0.603957,0.591709,0.577655,0.561834,0.544400,
1129 > 0.525193,0.504891,0.483435,0.461525,0.439729,0.418166,
1130 > 0.396478,0.375064,0.353725,0.332565,0.311308,0.290188,
1131 > 0.269125,0.248171,0.227323,0.206600,0.186070,0.165736,
1132 > 0.145686,0.125993,0.106807,0.088112,0.070133,0.053112,
1133 > 0.037464,0.023751,0.012665,0.004991,0.001015,0.000006,
1134 > 0.579569,0.570010,0.559085,0.546465,0.532125,0.516026,
1135 > 0.498089,0.478853,0.458681,0.438173,0.417404,0.396958,
1136 > 0.376499,0.356301,0.336115,0.316014,0.295930,0.275886,
1137 > 0.255871,0.236010,0.216223,0.196524,0.176991,0.157626,
1138 > 0.138515,0.119742,0.101404,0.083545,0.066339,0.050078,
1139 > 0.035167,0.022152,0.011703,0.004550,0.000904,0.000005,
1140 > 0.548753,0.540580,0.530520,0.519276,0.506004,0.491061,
1141 > 0.474436,0.456132,0.436996,0.417488,0.397755,0.378421,
1142 > 0.359085,0.339851,0.320700,0.301518,0.282357,0.263296,
1143 > 0.244273,0.225311,0.206421,0.187624,0.168959,0.150424,
1144 > 0.132145,0.114172,0.096557,0.079446,0.062967,0.047378,
1145 > 0.033126,0.020736,0.010987,0.004166,0.000809,0.000004,
1146 > 2.486670,2.401390,2.313150,2.223100,2.132210,2.040620,
1147 > 1.946240,1.852960,1.758660,1.663790,1.569300,1.475400,
1148 > 1.381630,1.288980,1.197380,1.106900,1.017420,0.930068,
1149 > 0.844430,0.760955,0.680549,0.603582,0.530075,0.460570,
1150 > 0.395641,0.335338,0.279919,0.229212,0.183101,0.141385,
1151 > 0.103826,0.070754,0.041898,0.019580,0.005430,0.000056,
1152 > 1.776150,1.719900,1.661320,1.600440,1.539110,1.477010,
1153 > 1.412790,1.348550,1.283220,1.218900,1.153830,1.088770,
1154 > 1.025080,0.961812,0.898582,0.836002,0.773940,0.712963,
1155 > 0.652650,0.593610,0.536099,0.480001,0.425916,0.373940,
1156 > 0.324454,0.277574,0.233716,0.192654,0.154624,0.119622,
1157 > 0.087676,0.059112,0.034659,0.015816,0.004189,0.000039,
1158 > 1.403340,1.361000,1.316930,1.271610,1.225380,1.177160,
1159 > 1.128240,1.078580,1.028640,0.977947,0.928063,0.878521,
1160 > 0.828550,0.779706,0.730569,0.682199,0.633948,0.586522,
1161 > 0.539401,0.493171,0.447639,0.403145,0.359673,0.317729,
1162 > 0.277222,0.238454,0.201716,0.166913,0.134251,0.103840,
1163 > 0.075864,0.050772,0.029388,0.013129,0.003341,0.000029,
1164 > 1.169950,1.136380,1.102100,1.065820,1.028190,0.989956,
1165 > 0.950189,0.908974,0.867893,0.826368,0.784940,0.743815,
1166 > 0.702872,0.662550,0.622187,0.582284,0.542431,0.503183,
1167 > 0.464112,0.425583,0.387521,0.350185,0.313576,0.277943,
1168 > 0.243398,0.210036,0.178125,0.147636,0.118821,0.091774,
1169 > 0.066794,0.044372,0.025371,0.011112,0.002727,0.000022,
1170 > 1.011370,0.984724,0.956435,0.926814,0.895443,0.862964,
1171 > 0.829274,0.794479,0.758875,0.722999,0.687335,0.652160,
1172 > 0.616882,0.582106,0.547373,0.512881,0.478596,0.444627,
1173 > 0.410813,0.377446,0.344502,0.311951,0.279953,0.248594,
1174 > 0.218124,0.188529,0.160091,0.132732,0.106696,0.082174,
1175 > 0.059499,0.039204,0.022122,0.009492,0.002245,0.000017,
1176 > 0.906689,0.883943,0.860162,0.834079,0.807067,0.778589,
1177 > 0.748446,0.717460,0.685716,0.653595,0.621726,0.590095,
1178 > 0.558631,0.527321,0.496368,0.465511,0.434680,0.404249,
1179 > 0.373821,0.343908,0.314123,0.284707,0.255819,0.227414,
1180 > 0.199673,0.172714,0.146637,0.121510,0.097519,0.074849,
1181 > 0.053902,0.035215,0.019623,0.008259,0.001888,0.000014,
1182 > 0.824084,0.805081,0.784198,0.761980,0.737877,0.712401,
1183 > 0.685185,0.657419,0.628573,0.599438,0.570456,0.541554,
1184 > 0.512854,0.484478,0.456174,0.428019,0.399930,0.372215,
1185 > 0.344523,0.317092,0.289826,0.262998,0.236396,0.210320,
1186 > 0.184715,0.159801,0.135667,0.112306,0.089938,0.068791,
1187 > 0.049274,0.031936,0.017585,0.007270,0.001610,0.000011,
1188 > 0.756740,0.739909,0.722048,0.702184,0.681164,0.658293,
1189 > 0.633526,0.608140,0.581469,0.554572,0.527904,0.501517,
1190 > 0.475145,0.448970,0.422970,0.397121,0.371248,0.345564,
1191 > 0.320059,0.294788,0.269623,0.244755,0.220148,0.195923,
1192 > 0.172134,0.148871,0.126312,0.104417,0.083450,0.063624,
1193 > 0.045337,0.029156,0.015873,0.006450,0.001387,0.000009,
1194 > 0.700680,0.686067,0.670562,0.652740,0.633579,0.612848,
1195 > 0.590476,0.566706,0.542207,0.517556,0.492696,0.468105,
1196 > 0.443562,0.419350,0.395237,0.371204,0.347185,0.323339,
1197 > 0.299543,0.276011,0.252614,0.229355,0.206336,0.183677,
1198 > 0.161386,0.139527,0.118283,0.097686,0.077903,0.059183,
1199 > 0.041956,0.026789,0.014429,0.005771,0.001207,0.000008,
1200 > 0.652991,0.640812,0.626563,0.610980,0.593505,0.574586,
1201 > 0.553939,0.531965,0.509196,0.485938,0.462796,0.439793,
1202 > 0.416876,0.394269,0.371716,0.349126,0.326720,0.304355,
1203 > 0.282067,0.259966,0.237971,0.216140,0.194490,0.173138,
1204 > 0.152089,0.131483,0.111379,0.091839,0.073081,0.055347,
1205 > 0.039055,0.024764,0.013208,0.005204,0.001061,0.000006,
1206 > 0.612493,0.601527,0.589260,0.574987,0.559293,0.541793,
1207 > 0.522714,0.502071,0.480580,0.459021,0.437202,0.415483,
1208 > 0.394052,0.372703,0.351511,0.330271,0.309080,0.288012,
1209 > 0.267062,0.246156,0.225353,0.204703,0.184235,0.163954,
1210 > 0.144013,0.124421,0.105317,0.086709,0.068861,0.051973,
1211 > 0.036498,0.022991,0.012143,0.004716,0.000938,0.000005,
1212 > 0.577224,0.567718,0.556768,0.543992,0.529614,0.513264,
1213 > 0.495593,0.476332,0.456118,0.435643,0.414958,0.394475,
1214 > 0.374159,0.353948,0.333917,0.313860,0.293761,0.273842,
1215 > 0.253865,0.234066,0.214285,0.194711,0.175202,0.155932,
1216 > 0.136903,0.118198,0.099962,0.082192,0.065107,0.048988,
1217 > 0.034247,0.021430,0.011215,0.004297,0.000835,0.000005,
1218 > 0.546580,0.538342,0.528410,0.516932,0.503593,0.488594,
1219 > 0.471844,0.453700,0.434588,0.415081,0.395496,0.376054,
1220 > 0.356702,0.337597,0.318488,0.299383,0.280225,0.261271,
1221 > 0.242274,0.223404,0.204576,0.185839,0.167243,0.148787,
1222 > 0.130586,0.112676,0.095184,0.078140,0.061781,0.046329,
1223 > 0.032244,0.020050,0.010400,0.003932,0.000747,0.000004/
1224 DATA LQ
1225 > /0.514104,0.515832,0.517091,0.517536,0.517832,0.517221,
1226 > 0.516574,0.515224,0.513287,0.510858,0.508000,0.504917,
1227 > 0.501246,0.497403,0.492983,0.488152,0.482358,0.476122,
1228 > 0.468649,0.460339,0.450645,0.439812,0.427166,0.412836,
1229 > 0.396439,0.377753,0.356388,0.331566,0.302870,0.269780,
1230 > 0.231439,0.187668,0.138736,0.086850,0.036211,0.000604,
1231 > 0.579741,0.577476,0.575013,0.571386,0.567464,0.562828,
1232 > 0.557414,0.551410,0.544332,0.536891,0.528953,0.521049,
1233 > 0.512612,0.503967,0.494945,0.485619,0.475531,0.464958,
1234 > 0.453985,0.442176,0.429817,0.416392,0.402026,0.386267,
1235 > 0.369251,0.350413,0.329631,0.306106,0.279316,0.248539,
1236 > 0.213030,0.172492,0.127054,0.078865,0.032177,0.000513,
1237 > 0.617205,0.612889,0.607448,0.601280,0.593987,0.585725,
1238 > 0.577239,0.567645,0.557436,0.546468,0.535124,0.523690,
1239 > 0.511913,0.500192,0.488271,0.476132,0.463482,0.450653,
1240 > 0.437460,0.423920,0.409767,0.395178,0.379804,0.363620,
1241 > 0.346492,0.327958,0.307882,0.285606,0.260283,0.231441,
1242 > 0.198228,0.160224,0.117620,0.072457,0.029003,0.000445,
1243 > 0.635918,0.629589,0.622092,0.614001,0.604751,0.594717,
1244 > 0.583683,0.572041,0.559188,0.546115,0.532523,0.518899,
1245 > 0.505201,0.491782,0.477675,0.463752,0.449738,0.435623,
1246 > 0.421062,0.406484,0.391587,0.376279,0.360631,0.344349,
1247 > 0.327299,0.309280,0.289970,0.268596,0.244665,0.217365,
1248 > 0.185991,0.150104,0.109826,0.067166,0.026423,0.000392,
1249 > 0.642211,0.634584,0.625818,0.616705,0.605951,0.594718,
1250 > 0.582351,0.569037,0.554986,0.540276,0.525303,0.510521,
1251 > 0.495596,0.480637,0.465827,0.451058,0.435872,0.421066,
1252 > 0.405880,0.390796,0.375397,0.359903,0.344185,0.327987,
1253 > 0.311213,0.293614,0.274949,0.254454,0.231575,0.205516,
1254 > 0.175652,0.141477,0.103138,0.062618,0.024219,0.000348,
1255 > 0.639180,0.631092,0.621893,0.612137,0.600935,0.589006,
1256 > 0.576122,0.561754,0.547048,0.531226,0.515512,0.499878,
1257 > 0.484089,0.468662,0.453258,0.437735,0.422292,0.407026,
1258 > 0.391544,0.376258,0.360735,0.345208,0.329470,0.313579,
1259 > 0.297158,0.280066,0.261996,0.242200,0.220305,0.195335,
1260 > 0.166765,0.134023,0.097367,0.058717,0.022356,0.000312,
1261 > 0.631156,0.622788,0.613580,0.603527,0.592277,0.579779,
1262 > 0.566484,0.551813,0.536423,0.520508,0.504251,0.488021,
1263 > 0.472093,0.456329,0.440551,0.424908,0.409169,0.393745,
1264 > 0.378229,0.362803,0.347473,0.331992,0.316442,0.300771,
1265 > 0.284762,0.268174,0.250682,0.231645,0.210500,0.186498,
1266 > 0.158976,0.127569,0.092372,0.055348,0.020770,0.000282,
1267 > 0.619698,0.611519,0.602336,0.591943,0.580814,0.568514,
1268 > 0.554982,0.540509,0.524970,0.508623,0.492412,0.476046,
1269 > 0.459935,0.444036,0.428166,0.412462,0.396815,0.381305,
1270 > 0.365892,0.350546,0.335395,0.320121,0.304877,0.289437,
1271 > 0.273841,0.257707,0.240774,0.222305,0.201894,0.178726,
1272 > 0.152209,0.121898,0.087984,0.052399,0.019398,0.000257,
1273 > 0.606343,0.598499,0.589365,0.579821,0.568434,0.556327,
1274 > 0.542872,0.528454,0.512900,0.496651,0.480390,0.464149,
1275 > 0.447972,0.432310,0.416261,0.400748,0.385116,0.369727,
1276 > 0.354542,0.339434,0.324280,0.309328,0.294433,0.279362,
1277 > 0.264086,0.248363,0.231853,0.214025,0.194261,0.171833,
1278 > 0.146197,0.116876,0.084113,0.049805,0.018206,0.000236,
1279 > 0.591896,0.584256,0.575797,0.566248,0.555494,0.543797,
1280 > 0.530664,0.516396,0.501105,0.484951,0.468541,0.452368,
1281 > 0.436358,0.420764,0.405075,0.389617,0.374238,0.359020,
1282 > 0.343972,0.329043,0.314268,0.299607,0.284884,0.270198,
1283 > 0.255337,0.239949,0.223985,0.206644,0.187466,0.165701,
1284 > 0.140843,0.112412,0.080664,0.047512,0.017165,0.000218,
1285 > 0.576697,0.569585,0.561998,0.552972,0.542790,0.531174,
1286 > 0.518335,0.504287,0.489160,0.473285,0.457139,0.441230,
1287 > 0.425520,0.409903,0.394534,0.379177,0.364141,0.349046,
1288 > 0.334214,0.319602,0.305086,0.290690,0.276241,0.261890,
1289 > 0.247381,0.232411,0.216817,0.199969,0.181340,0.160182,
1290 > 0.135997,0.108375,0.077551,0.045446,0.016236,0.000202,
1291 > 0.561779,0.555491,0.548244,0.539719,0.529917,0.518731,
1292 > 0.506334,0.492791,0.477798,0.462049,0.446208,0.430534,
1293 > 0.414938,0.399619,0.384500,0.369572,0.354611,0.339912,
1294 > 0.325296,0.310898,0.296649,0.282465,0.268439,0.254375,
1295 > 0.240158,0.225578,0.210376,0.193910,0.175777,0.155185,
1296 > 0.131618,0.104727,0.074745,0.043586,0.015406,0.000188,
1297 > 0.547406,0.541522,0.534772,0.526637,0.517320,0.506665,
1298 > 0.494929,0.481541,0.466958,0.451581,0.435977,0.420382,
1299 > 0.405163,0.390039,0.375242,0.360562,0.345699,0.331352,
1300 > 0.316934,0.302875,0.288897,0.275075,0.261211,0.247494,
1301 > 0.233546,0.219335,0.204477,0.188465,0.170694,0.150613,
1302 > 0.127623,0.101395,0.072188,0.041898,0.014659,0.000176,
1303 > 0.510949,0.512461,0.513125,0.513775,0.513627,0.512953,
1304 > 0.511676,0.510153,0.508026,0.505586,0.502488,0.499514,
1305 > 0.495918,0.491777,0.487221,0.482318,0.476322,0.469865,
1306 > 0.462254,0.453539,0.443579,0.431957,0.418405,0.402714,
1307 > 0.384513,0.363211,0.338441,0.309525,0.276023,0.237418,
1308 > 0.194095,0.146991,0.098692,0.053921,0.018675,0.000251,
1309 > 0.573840,0.571539,0.568660,0.564636,0.560125,0.555429,
1310 > 0.549250,0.542654,0.535827,0.528339,0.520207,0.512007,
1311 > 0.503505,0.494854,0.485721,0.476289,0.466253,0.455649,
1312 > 0.444534,0.432471,0.419685,0.405709,0.390526,0.373656,
1313 > 0.355038,0.334208,0.310578,0.283443,0.252376,0.217005,
1314 > 0.177144,0.133976,0.089700,0.048645,0.016531,0.000213,
1315 > 0.610274,0.604837,0.598760,0.591967,0.584049,0.575825,
1316 > 0.566634,0.556589,0.545957,0.534847,0.523533,0.511998,
1317 > 0.500328,0.488702,0.476646,0.464532,0.451851,0.438900,
1318 > 0.425698,0.412022,0.397498,0.382650,0.366763,0.349574,
1319 > 0.330976,0.310680,0.288132,0.262614,0.233680,0.200673,
1320 > 0.163728,0.123628,0.082475,0.044455,0.014854,0.000185,
1321 > 0.627192,0.620322,0.611968,0.603362,0.593426,0.582683,
1322 > 0.571079,0.558798,0.545824,0.532347,0.519069,0.505318,
1323 > 0.491597,0.478201,0.464023,0.450314,0.436183,0.422097,
1324 > 0.407576,0.392963,0.377856,0.362363,0.346096,0.329071,
1325 > 0.310942,0.291304,0.269699,0.245618,0.218210,0.187279,
1326 > 0.152695,0.115069,0.076563,0.041003,0.013496,0.000163,
1327 > 0.631625,0.623836,0.614536,0.604467,0.593362,0.581189,
1328 > 0.568126,0.554473,0.540005,0.525342,0.510320,0.495556,
1329 > 0.480543,0.465856,0.450901,0.436047,0.421171,0.406166,
1330 > 0.391154,0.375993,0.360616,0.344868,0.328695,0.311826,
1331 > 0.294125,0.275057,0.254485,0.231491,0.205428,0.176164,
1332 > 0.143362,0.107930,0.071534,0.038049,0.012338,0.000145,
1333 > 0.628674,0.619289,0.609887,0.598892,0.587603,0.574461,
1334 > 0.560634,0.545957,0.530648,0.515019,0.499229,0.483499,
1335 > 0.467857,0.452587,0.437187,0.421817,0.406489,0.391076,
1336 > 0.375987,0.360550,0.345003,0.329399,0.313389,0.296812,
1337 > 0.279488,0.261191,0.241276,0.219209,0.194576,0.166551,
1338 > 0.135508,0.101696,0.067217,0.035540,0.011363,0.000130,
1339 > 0.619655,0.610842,0.600741,0.589506,0.577444,0.564098,
1340 > 0.550163,0.534861,0.519180,0.503427,0.486955,0.471011,
1341 > 0.455040,0.439305,0.423732,0.408132,0.392629,0.377169,
1342 > 0.361831,0.346498,0.331073,0.315724,0.299859,0.283638,
1343 > 0.266810,0.249063,0.229985,0.208716,0.185019,0.158309,
1344 > 0.128580,0.096355,0.063475,0.033376,0.010535,0.000117,
1345 > 0.607474,0.598863,0.589070,0.578071,0.565959,0.552353,
1346 > 0.538238,0.523141,0.507028,0.491211,0.474478,0.458296,
1347 > 0.442260,0.426619,0.410731,0.395250,0.379629,0.364356,
1348 > 0.348971,0.333841,0.318608,0.303265,0.287813,0.272025,
1349 > 0.255702,0.238406,0.219910,0.199549,0.176755,0.151092,
1350 > 0.122529,0.091666,0.060228,0.031496,0.009822,0.000107,
1351 > 0.593760,0.585285,0.576038,0.564716,0.553007,0.539720,
1352 > 0.525715,0.510544,0.494716,0.478059,0.461888,0.445882,
1353 > 0.429808,0.413922,0.398545,0.382937,0.367557,0.352300,
1354 > 0.337238,0.322197,0.307168,0.292182,0.277095,0.261674,
1355 > 0.245761,0.228939,0.211113,0.191442,0.169413,0.144749,
1356 > 0.117225,0.087581,0.057361,0.029844,0.009202,0.000098,
1357 > 0.579324,0.570883,0.561883,0.551297,0.539642,0.526840,
1358 > 0.512914,0.497903,0.481998,0.465838,0.449788,0.433693,
1359 > 0.417703,0.402286,0.386792,0.371437,0.356215,0.341283,
1360 > 0.326313,0.311637,0.296947,0.282140,0.267321,0.252270,
1361 > 0.236761,0.220580,0.203236,0.184230,0.162954,0.139137,
1362 > 0.112557,0.083990,0.054824,0.028391,0.008663,0.000090,
1363 > 0.564315,0.556562,0.547782,0.537905,0.526572,0.514011,
1364 > 0.500033,0.485378,0.470013,0.453915,0.438022,0.422153,
1365 > 0.406530,0.391112,0.375843,0.360792,0.345939,0.331089,
1366 > 0.316383,0.301913,0.287442,0.273075,0.258564,0.243883,
1367 > 0.228827,0.213013,0.196181,0.177712,0.157167,0.134036,
1368 > 0.108376,0.080670,0.052555,0.027085,0.008183,0.000084,
1369 > 0.549474,0.542540,0.534001,0.524222,0.513538,0.501429,
1370 > 0.487998,0.473665,0.458530,0.442760,0.426794,0.411187,
1371 > 0.395839,0.380807,0.365652,0.350917,0.336154,0.321565,
1372 > 0.307328,0.293072,0.278886,0.264801,0.250668,0.236339,
1373 > 0.221517,0.206147,0.189835,0.171902,0.151860,0.129450,
1374 > 0.104545,0.077707,0.050505,0.025913,0.007753,0.000078,
1375 > 0.535066,0.528539,0.520546,0.511350,0.500893,0.489256,
1376 > 0.476407,0.462092,0.447135,0.431648,0.416302,0.400932,
1377 > 0.385885,0.370871,0.356128,0.341562,0.327354,0.312933,
1378 > 0.298890,0.284864,0.270948,0.257177,0.243326,0.229321,
1379 > 0.214967,0.199986,0.183994,0.166574,0.147078,0.125313,
1380 > 0.101089,0.075065,0.048640,0.024855,0.007368,0.000073,
1381 > 0.509409,0.510722,0.511674,0.511738,0.511286,0.510766,
1382 > 0.509198,0.507564,0.505152,0.502704,0.499650,0.496452,
1383 > 0.492756,0.488747,0.483986,0.478817,0.472835,0.466164,
1384 > 0.458412,0.449362,0.438892,0.426796,0.412430,0.395391,
1385 > 0.375497,0.352110,0.324651,0.292275,0.254919,0.212461,
1386 > 0.165886,0.117521,0.071431,0.033526,0.009193,0.000091,
1387 > 0.570888,0.568118,0.564921,0.561074,0.555970,0.550478,
1388 > 0.544438,0.537826,0.530394,0.523005,0.514685,0.506342,
1389 > 0.497778,0.489102,0.479852,0.470348,0.460202,0.449558,
1390 > 0.438121,0.426046,0.412891,0.398444,0.382493,0.364765,
1391 > 0.344689,0.321925,0.295886,0.265915,0.231610,0.192817,
1392 > 0.150406,0.106352,0.064426,0.030028,0.008099,0.000077,
1393 > 0.605452,0.599834,0.593296,0.586627,0.578039,0.569200,
1394 > 0.559754,0.549624,0.538672,0.527535,0.516054,0.504396,
1395 > 0.492657,0.481104,0.468974,0.456682,0.443949,0.431100,
1396 > 0.417683,0.403750,0.389161,0.373852,0.357151,0.339250,
1397 > 0.319633,0.297684,0.273111,0.245042,0.213167,0.177272,
1398 > 0.138115,0.097498,0.058870,0.027271,0.007249,0.000067,
1399 > 0.621534,0.614157,0.605835,0.596272,0.585946,0.574335,
1400 > 0.562822,0.550233,0.536841,0.523598,0.509920,0.496293,
1401 > 0.482476,0.468781,0.455150,0.441275,0.427053,0.412946,
1402 > 0.398426,0.383702,0.368436,0.352438,0.335827,0.317985,
1403 > 0.298878,0.277758,0.254411,0.228040,0.198146,0.164659,
1404 > 0.128112,0.090261,0.054331,0.025017,0.006563,0.000059,
1405 > 0.625646,0.616575,0.607120,0.596395,0.584780,0.572447,
1406 > 0.558731,0.544391,0.529948,0.515041,0.499832,0.485344,
1407 > 0.470145,0.455555,0.440783,0.425851,0.410922,0.396061,
1408 > 0.381038,0.365825,0.350234,0.334195,0.317648,0.300202,
1409 > 0.281570,0.261329,0.238985,0.213898,0.185675,0.154079,
1410 > 0.119726,0.084159,0.050488,0.023108,0.005981,0.000052,
1411 > 0.621816,0.611958,0.601503,0.590515,0.577713,0.564667,
1412 > 0.550311,0.535181,0.519529,0.503933,0.488235,0.472641,
1413 > 0.456837,0.441572,0.426210,0.410935,0.395680,0.380387,
1414 > 0.365098,0.349710,0.334114,0.318243,0.301798,0.284744,
1415 > 0.266642,0.247149,0.225778,0.201848,0.175065,0.145114,
1416 > 0.112575,0.078976,0.047163,0.021485,0.005492,0.000047,
1417 > 0.612456,0.602615,0.592112,0.580629,0.567547,0.554000,
1418 > 0.539266,0.523505,0.507494,0.491317,0.475170,0.459335,
1419 > 0.443356,0.427588,0.412089,0.396613,0.381018,0.365836,
1420 > 0.350418,0.335168,0.319712,0.303995,0.287942,0.271269,
1421 > 0.253737,0.234846,0.214374,0.191514,0.165911,0.137420,
1422 > 0.106436,0.074522,0.044401,0.020096,0.005079,0.000042,
1423 > 0.599906,0.590658,0.579871,0.568198,0.555552,0.541752,
1424 > 0.526573,0.511088,0.494690,0.478321,0.462067,0.446064,
1425 > 0.430078,0.414174,0.398700,0.383180,0.367731,0.352392,
1426 > 0.337125,0.322012,0.306729,0.291406,0.275582,0.259369,
1427 > 0.242384,0.224158,0.204420,0.182500,0.157984,0.130669,
1428 > 0.101087,0.070603,0.041960,0.018893,0.004724,0.000039,
1429 > 0.586385,0.576832,0.566175,0.554843,0.542371,0.528570,
1430 > 0.513481,0.497995,0.481584,0.465454,0.449166,0.432955,
1431 > 0.417165,0.401569,0.386012,0.370517,0.355274,0.340119,
1432 > 0.324992,0.310107,0.295066,0.280058,0.264679,0.248876,
1433 > 0.232336,0.214712,0.195681,0.174542,0.150960,0.124776,
1434 > 0.096383,0.067222,0.039821,0.017842,0.004416,0.000035,
1435 > 0.571562,0.562230,0.553483,0.541106,0.528820,0.515154,
1436 > 0.500724,0.484959,0.468895,0.452796,0.436478,0.420645,
1437 > 0.404877,0.389321,0.373971,0.358762,0.343687,0.328774,
1438 > 0.313968,0.299216,0.284619,0.269828,0.254780,0.239398,
1439 > 0.223345,0.206321,0.187886,0.167537,0.144752,0.119523,
1440 > 0.092229,0.064212,0.037929,0.016920,0.004149,0.000033,
1441 > 0.556430,0.547716,0.538118,0.527507,0.515330,0.502167,
1442 > 0.487523,0.472378,0.456637,0.440513,0.424536,0.408900,
1443 > 0.393308,0.378059,0.362874,0.347899,0.333001,0.318366,
1444 > 0.303866,0.289409,0.275001,0.260570,0.245889,0.230899,
1445 > 0.215327,0.198784,0.180938,0.161191,0.139240,0.114867,
1446 > 0.088515,0.061527,0.036242,0.016094,0.003912,0.000030,
1447 > 0.541088,0.533114,0.524472,0.513843,0.502197,0.489243,
1448 > 0.475392,0.460194,0.444773,0.428968,0.413337,0.397800,
1449 > 0.382488,0.367387,0.352568,0.337799,0.323189,0.308771,
1450 > 0.294545,0.280383,0.266247,0.252107,0.237863,0.223231,
1451 > 0.208050,0.191998,0.174706,0.155516,0.134237,0.110630,
1452 > 0.085174,0.059101,0.034727,0.015356,0.003700,0.000028,
1453 > 0.526746,0.519306,0.510794,0.500680,0.489449,0.477026,
1454 > 0.463177,0.448622,0.433350,0.417867,0.402498,0.387171,
1455 > 0.372235,0.357466,0.342866,0.328463,0.314082,0.299986,
1456 > 0.285991,0.272106,0.258305,0.244422,0.230511,0.216246,
1457 > 0.201466,0.185822,0.169011,0.150402,0.129717,0.106819,
1458 > 0.082146,0.056897,0.033348,0.014690,0.003511,0.000026,
1459 > 0.508294,0.509537,0.510232,0.510364,0.510028,0.509066,
1460 > 0.507574,0.506096,0.503764,0.501269,0.498261,0.494940,
1461 > 0.491357,0.487139,0.482585,0.477518,0.471487,0.464685,
1462 > 0.456954,0.448138,0.437547,0.425359,0.410953,0.393893,
1463 > 0.373877,0.350282,0.322653,0.290145,0.252556,0.209889,
1464 > 0.163221,0.114885,0.069163,0.031966,0.008546,0.000081,
1465 > 0.569358,0.566803,0.562983,0.559001,0.554723,0.548955,
1466 > 0.542742,0.536150,0.528626,0.520981,0.513074,0.504898,
1467 > 0.496364,0.487653,0.478428,0.468815,0.458593,0.447996,
1468 > 0.436563,0.424615,0.411428,0.396895,0.381014,0.363131,
1469 > 0.342953,0.320173,0.294051,0.263852,0.229272,0.190367,
1470 > 0.147875,0.103902,0.062322,0.028615,0.007526,0.000069,
1471 > 0.603670,0.598438,0.591834,0.584727,0.576222,0.567680,
1472 > 0.558196,0.547591,0.536787,0.525654,0.514148,0.502844,
1473 > 0.491102,0.479283,0.467390,0.455046,0.442422,0.429358,
1474 > 0.416062,0.402212,0.387705,0.372222,0.355668,0.337678,
1475 > 0.318005,0.295922,0.271145,0.242997,0.210936,0.174968,
1476 > 0.135724,0.095201,0.057278,0.025971,0.006733,0.000060,
1477 > 0.619938,0.612660,0.603237,0.594399,0.584220,0.572945,
1478 > 0.560906,0.547974,0.535000,0.521640,0.508074,0.494503,
1479 > 0.480788,0.467052,0.453309,0.439430,0.425326,0.411154,
1480 > 0.396834,0.381957,0.366731,0.350903,0.334155,0.316435,
1481 > 0.297197,0.275976,0.252539,0.225982,0.196020,0.162415,
1482 > 0.125838,0.088092,0.052512,0.023816,0.006094,0.000053,
1483 > 0.623752,0.615162,0.604853,0.594199,0.582617,0.570319,
1484 > 0.556672,0.542512,0.527993,0.513108,0.498185,0.483234,
1485 > 0.468378,0.453629,0.438899,0.424139,0.409158,0.394331,
1486 > 0.379329,0.364110,0.348592,0.332662,0.316045,0.298521,
1487 > 0.279860,0.259548,0.237136,0.211942,0.183600,0.151947,
1488 > 0.117534,0.082104,0.048767,0.021982,0.005551,0.000047,
1489 > 0.619727,0.610281,0.600069,0.588100,0.575830,0.562394,
1490 > 0.548038,0.533038,0.517667,0.501730,0.486231,0.470564,
1491 > 0.454906,0.439606,0.424254,0.409150,0.393803,0.378587,
1492 > 0.363406,0.347957,0.332413,0.316594,0.300201,0.283135,
1493 > 0.264987,0.245411,0.223977,0.199987,0.173072,0.143020,
1494 > 0.110486,0.077018,0.045586,0.020431,0.005096,0.000042,
1495 > 0.610495,0.600588,0.589924,0.578181,0.565270,0.551795,
1496 > 0.537196,0.521624,0.505437,0.489150,0.472983,0.457321,
1497 > 0.441422,0.425796,0.410213,0.394710,0.379352,0.363971,
1498 > 0.348802,0.333422,0.317927,0.302300,0.286286,0.269619,
1499 > 0.252021,0.233168,0.212605,0.189631,0.163957,0.135368,
1500 > 0.104416,0.072632,0.042854,0.019101,0.004711,0.000038,
1501 > 0.597856,0.588522,0.577800,0.566289,0.553520,0.539436,
1502 > 0.524435,0.508812,0.492453,0.476231,0.459950,0.443865,
1503 > 0.427942,0.412263,0.396763,0.381179,0.365851,0.350734,
1504 > 0.335459,0.320164,0.305048,0.289712,0.274016,0.257736,
1505 > 0.240735,0.222511,0.202688,0.180671,0.156046,0.128701,
1506 > 0.099140,0.068812,0.040482,0.017948,0.004381,0.000035,
1507 > 0.584084,0.574752,0.564446,0.552705,0.540225,0.526223,
1508 > 0.511434,0.495863,0.479694,0.463424,0.447226,0.431058,
1509 > 0.415167,0.399560,0.384026,0.368648,0.353514,0.338375,
1510 > 0.323406,0.308441,0.293400,0.278400,0.263020,0.247257,
1511 > 0.230713,0.213051,0.193943,0.172756,0.149112,0.122832,
1512 > 0.094502,0.065475,0.038404,0.016946,0.004094,0.000032,
1513 > 0.569517,0.560385,0.550236,0.539113,0.526421,0.513222,
1514 > 0.498378,0.482790,0.466696,0.450744,0.434441,0.418603,
1515 > 0.402901,0.387328,0.372129,0.356950,0.341855,0.326972,
1516 > 0.312231,0.297539,0.282901,0.268148,0.253147,0.237801,
1517 > 0.221735,0.204645,0.186195,0.165769,0.142932,0.117672,
1518 > 0.090399,0.062530,0.036573,0.016065,0.003846,0.000029,
1519 > 0.554361,0.545724,0.536138,0.525127,0.513348,0.500011,
1520 > 0.485538,0.470383,0.454424,0.438520,0.422677,0.406888,
1521 > 0.391373,0.376154,0.360890,0.345988,0.331296,0.316666,
1522 > 0.302083,0.287727,0.273311,0.258904,0.244293,0.229265,
1523 > 0.213676,0.197142,0.179286,0.159481,0.137449,0.113024,
1524 > 0.086731,0.059884,0.034929,0.015279,0.003625,0.000027,
1525 > 0.539566,0.531375,0.522021,0.511690,0.499982,0.487119,
1526 > 0.474212,0.458150,0.442576,0.426822,0.411177,0.395600,
1527 > 0.380520,0.365481,0.350710,0.336033,0.321423,0.307074,
1528 > 0.292885,0.278746,0.264660,0.250555,0.236263,0.221649,
1529 > 0.206449,0.190350,0.173011,0.154255,0.132811,0.108865,
1530 > 0.083444,0.057506,0.033461,0.014573,0.003428,0.000025,
1531 > 0.524865,0.517319,0.508721,0.498692,0.487348,0.474698,
1532 > 0.461245,0.446504,0.431299,0.416072,0.400507,0.385278,
1533 > 0.370252,0.355605,0.341042,0.326661,0.312354,0.298301,
1534 > 0.284289,0.270473,0.256711,0.242903,0.228985,0.214684,
1535 > 0.199885,0.184233,0.167367,0.148740,0.128003,0.105093,
1536 > 0.080442,0.055358,0.032128,0.013936,0.003252,0.000024/
1537
1538C IF (X.GE.XMAXX.OR.Q2.LT.Q2MIN.OR.Q2.GT.Q2MAX) THEN
1539C WRITE(*,*) 'X or Q2 out of range!!'
1540C WRITE(*,*) 'Valid range is:'
1541C WRITE(*,*) '1E-05 < X < 0.95 and'
1542C WRITE(*,*) '4 < Q2 < 520 GeV**2'
1543C RETURN
1544C END IF
1545
1546 IF (X.GE.XMAX.OR.B.GT.BMAX) THEN
1547 RES(1) = 1.0
1548 RES(2) = 1.0
1549 RETURN
1550 END IF
1551 IF(Q2.LT.Q2MIN) THEN
1552 Q2 = Q2MIN
1553 ELSE IF(Q2.GT.Q2MAX) THEN
1554 Q2 = Q2MAX
1555 ENDIF
1556
1557 CALL GGINTER(INUCL,X,Q2,SHAD)
1558
1559 DO I=1,31
1560 TATMP(I) = TA(I,INUCL)
1561 ENDDO
1562 CALL SPLINE(TMAX,IMPAR,TATMP,C,D,E)
1563 TAF = SEVAL(TMAX,B,IMPAR,TATMP,C,D,E)
1564
1565 DO IK=1,2
1566 RES(IK) = 1.0/(1.0 + (ANUCL(INUCL)-1.0)*SHAD(IK)*TAF)
1567 ENDDO
1568
1569 RETURN
1570 END
1571
1572 SUBROUTINE GGINTER(INUCL,X,Q2,SHAD)
1573 IMPLICIT DOUBLE PRECISION(A-H,L-Z)
1574 IMPLICIT INTEGER(I-K)
1575
1576 DIMENSION SHAD(2)
1577 DIMENSION G(36,13,4),LQ(36,13,4)
1578 DIMENSION XB(36),Q2V(13)
1579 DIMENSION GQTMP(13),LQQTMP(13)
1580 DIMENSION GXTMP(36),LQXTMP(36)
1581 DIMENSION C(100),D(100),E(100)
1582
1583 PARAMETER(INMAX=13,IMMAX=36)
1584
1585 COMMON/GG07/ XB,Q2V,G,LQ
1586
1587 DO K=1,36
1588 DO J=1,13
1589 GQTMP(J) = G(K,J,INUCL)
1590 LQQTMP(J) = LQ(K,J,INUCL)
1591 ENDDO
1592 CALL SPLINE(INMAX,Q2V,GQTMP,C,D,E)
1593 GXTMP(K) = SEVAL(INMAX,Q2,Q2V,GQTMP,C,D,E)
1594 CALL SPLINE(INMAX,Q2V,LQQTMP,C,D,E)
1595 LQXTMP(K) = SEVAL(INMAX,Q2,Q2V,LQQTMP,C,D,E)
1596 ENDDO
1597
1598 CALL SPLINE(IMMAX,XB,GXTMP,C,D,E)
1599 SHAD(1) = SEVAL(IMMAX,X,XB,GXTMP,C,D,E)
1600 CALL SPLINE(IMMAX,XB,LQXTMP,C,D,E)
1601 SHAD(2) = SEVAL(IMMAX,X,XB,LQXTMP,C,D,E)
1602
1603 RETURN
1604 END
1605
1606C ---------------------------------------------------------------------
1607 SUBROUTINE SPLINE(N,X,Y,B,C,D)
1608C ---------------------------------------------------------------------
1609c***************************************************************************
1610C CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
1611C INTERPOLATION SUBROUTINES ARE TAKEN FROM
1612C G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
1613C COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
1614 IMPLICIT double precision (A-H,O-Z)
1615 DIMENSION X(100),Y(100),B(100),C(100),D(100)
1616 NM1=N-1
1617 IF(N.LT.2) RETURN
1618 IF(N.LT.3) GO TO 250
1619 D(1)=X(2)-X(1)
1620 C(2)=(Y(2)-Y(1))/D(1)
1621 DO 210 I=2,NM1
1622 D(I)=X(I+1)-X(I)
1623 B(I)=2.0D0*(D(I-1)+D(I))
1624 C(I+1)=(Y(I+1)-Y(I))/D(I)
1625 C(I)=C(I+1)-C(I)
1626 210 CONTINUE
1627 B(1)=-D(1)
1628 B(N)=-D(N-1)
1629 C(1)=0.0D0
1630 C(N)=0.0D0
1631 IF(N.EQ.3) GO TO 215
1632 C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
1633 C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
1634 C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1))
1635 C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3))
1636 215 CONTINUE
1637 DO 220 I=2,N
1638 T=D(I-1)/B(I-1)
1639 B(I)=B(I)-T*D(I-1)
1640 C(I)=C(I)-T*C(I-1)
1641 220 CONTINUE
1642 C(N)=C(N)/B(N)
1643 DO 230 IB=1,NM1
1644 I=N-IB
1645 C(I)=(C(I)-D(I)*C(I+1))/B(I)
1646 230 CONTINUE
1647 B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.0D0*C(N))
1648 DO 240 I=1,NM1
1649 B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.0D0*C(I))
1650 D(I)=(C(I+1)-C(I))/D(I)
1651 C(I)=3.0D0*C(I)
1652 240 CONTINUE
1653 C(N)=3.0D0*C(N)
1654 D(N)=D(N-1)
1655 RETURN
1656 250 CONTINUE
1657 B(1)=(Y(2)-Y(1))/(X(2)-X(1))
1658 C(1)=0.0D0
1659 D(1)=0.0D0
1660 B(2)=B(1)
1661 C(2)=0.0D0
1662 D(2)=0.0D0
1663 RETURN
1664 END
1665c
1666c***************************************************************************
1667C ---------------------------------------------------------------------
1668 double precision FUNCTION SEVAL(N,XX,X,Y,B,C,D)
1669C ---------------------------------------------------------------------
1670c***************************************************************************
1671C CALCULATE THE DISTRIBUTION AT XX BY CUBIC SPLINE INTERPOLATION.
1672 implicit double precision(A-H,O-Z)
1673 DIMENSION X(100),Y(100),B(100),C(100),D(100)
1674 DATA I/1/
1675 IF(I.GE.N) I=1
1676 IF(XX.LT.X(I)) GO TO 310
1677 IF(XX.LE.X(I+1)) GO TO 330
1678 310 CONTINUE
1679 I=1
1680 J=N+1
1681 320 CONTINUE
1682 K=(I+J)/2
1683 IF(XX.LT.X(K)) J=K
1684 IF(XX.GE.X(K)) I=K
1685 IF(J.GT.I+1) GO TO 320
1686 330 CONTINUE
1687 DX=XX-X(I)
1688 SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))
1689 RETURN
1690 END
1691c
1692c***************************************************************************
1693
1694
1695* function to generate gauss distribution
1696 double precision function gaus(x0,sig)
1697 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1698 41 u1=pyr(0)
1699 u2=pyr(0)
1700 v1=2.d0*u1-1.d0
1701 v2=2.d0*u2-1.d0
1702 s=v1**2+v2**2
1703 if(s.gt.1) go to 41
1704 gaus=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
1705 return
1706 end