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