LHAPDF veraion 5.9.1
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf-5.9.1 / src / wraphera.f
1       subroutine HERAevolve(x,Q,pdf)
2       implicit double precision (a-h,o-z)
3       include 'parmsetup.inc'
4       character*64 gridname
5       character*16 name(nmxset)
6       integer nmem(nmxset),ndef(nmxset),mem
7       common/NAME/name,nmem,ndef,mem
8       integer nset,iset,isetlast
9       data isetlast/-1/
10       integer Eorder
11       real*8 mc,mc2,mb,mb2,mt,mt2
12       real*8 f(-6:6),pdf(-6:6)
13       integer qnerr
14       parameter(nstartp=7)
15       DIMENSION QSP(NSTARTP)
16       DATA QSP/10.,20.,30.,40.,50.,80.,100./
17       real*8 xval(45)
18       logical HEAVY,VFN
19       real*8 pwgt(20)
20       save
21 !
22       call getnset(iset)
23 !      print *,'iset=',iset,' now calling get_pdfqcd'
24       if(iset.ne.isetlast) then
25         call get_pdfqcd(iset)
26         isetlast = iset
27       endif
28 !
29       Q2=Q*Q
30       UCENT=QPDFXQ('UPVAL',X,Q2,IFAIL)
31       DCENT=QPDFXQ('DNVAL',X,Q2,IFAIL)
32       GCENT=QPDFXQ('GLUON',X,Q2,IFAIL)
33       UBCEN=QPDFXQ('UB',X,Q2,IFAIL)
34       DBCEN=QPDFXQ('DB',X,Q2,IFAIL)
35       STCEN=QPDFXQ('SB',X,Q2,IFAIL)
36       IF(Q2.GE.Q2C)THEN
37         CHCEN=QPDFXQ('CB',X,Q2,IFAIL)
38       ELSE
39         CHCEN=0.0
40       ENDIF
41       IF(Q2.GE.Q2B)THEN
42         BTCEN=QPDFXQ('BB',X,Q2,IFAIL)
43       ELSE
44         BTCEN=0.0
45       ENDIF
46       pdf(1)=dcent+dbcen
47       pdf(2)=ucent+ubcen
48       pdf(3)=stcen
49       pdf(4)=chcen
50       pdf(5)=btcen
51       pdf(6)=0d0
52       pdf(0)=gcent
53       pdf(-1)=dbcen
54       pdf(-2)=ubcen
55       pdf(-3)=stcen
56       pdf(-4)=chcen
57       pdf(-5)=btcen
58       pdf(-6)=0d0
59 !
60       return
61 !
62 !===========================================================================
63       entry HERAalfa(alfas,Q)
64       Q2=Q*Q
65       nf=6
66       if (Q2.lt.mt2) nf=5
67       if (Q2.lt.mb2) nf=4
68       if (Q2.lt.mc2) nf=3
69       alfas=QALFAS(Q2,Qlam,nf,iflag)
70 !      print *,q2,alfas,qlam,nf,iflag
71       return
72 !
73 !===========================================================================
74       entry HERAread(nset)
75          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
76 !        print *,gridname,nx,xmin,xmax,nq,qmin,qmax,' from HERAread'
77       return
78 !
79 !===========================================================================
80       entry HERAinit(nset,Eorder,Q2fit)
81 !      print *,name(nset), ' from HERAinit'
82       if(name(nset).eq.'QCDNUM_HERA_TR') then
83          HEAVY = .FALSE.
84          VFN = .TRUE.
85       else if(name(nset).eq.'QCDNUM_HERA_FF') then
86          HEAVY = .TRUE.
87          VFN = .FALSE.
88       else if(name(nset).eq.'QCDNUM_HERA_ZM') then
89          HEAVY = .FALSE.
90          VFN = .FALSE.
91       else
92         print *,'name/scheme not recognized'
93         stop 1
94       endif
95 !--try 3 way logic ffn/zm-vfn/rt-vfn
96       IF(HEAVY)THEN
97       IVFN=1
98       ELSE
99       IVFN=0
100       ENDIF
101       IF(VFN)THEN
102       IVFN=IVFN+2
103       ELSE
104       IVFN=IVFN
105       ENDIF
106 ! IVFN=0 IS ZM-VFN, 1 IS FFN,2 IS RT-VFN, 3 IS NOT ALLOWED  
107
108       IF(IVFN.EQ.3)THEN
109       WRITE(*,*)'IVFN=3 SO STOP',IVFN
110       STOP
111       ENDIF    
112
113       
114       
115 !--qcdnum initialisation
116       CALL QNINIT
117 !--se thresholds
118       Q0=Q2fit
119       ZM=91.187D0
120       ZM2=ZM*ZM
121       ALPHAS=QNALFA(ZM2)
122
123       call getQmassM(nset,4,mc)
124       mc2=mc*mc
125       call getQmassM(nset,5,mb)
126       mb2=mb*mb
127       call getQmassM(nset,6,mt)
128       mt2=mt*mt
129     
130 !      Q2C=1.8225
131       Q2C=mc2
132 !      Q2B=18.49
133       Q2B=mb2 
134 !      print  *,q2c,q2b    
135       q2t=mt2
136
137       IF (Q0.LT.Q2C) THEN
138         NACT=3
139       ELSE
140         NACT=4
141       ENDIF
142 !--this merely defines nact where we startevolution
143 !--namely at q0
144       IF (HEAVY) NACT=3
145
146       CALL QNRSET('MCSTF',SQRT(Q2C))
147       CALL QNRSET('MBSTF',SQRT(Q2B))
148       CALL QNRSET('MCALF',SQRT(Q2C))
149       CALL QNRSET('MBALF',SQRT(Q2B))
150       CALL QNRSET('MTALF',SQRT(Q2T))
151     
152       IF (HEAVY) THEN
153         CALL QTHRES(1D10,2D10)
154 !        CALL QTHRES(1D6,2D6)
155       ELSE
156         CALL QTHRES(Q2C,Q2B)
157       ENDIF
158  
159       DO I=1,NSTARTP
160         CALL GRQINP(QSP(I),1)
161       ENDDO
162       CALL GRQINP(Q0,1)
163       CALL GRQINP(Q2C,1)
164       CALL GRQINP(Q2B,1)
165 !      qcdnum grid not my grid
166
167 !      CALL GRXLIM(120,97D-8)
168 !      print *,'nx = ',nx
169       CALL GRXLIM(nx,xmin)
170
171 !      CALL GRQLIM(61,29D-2,200D3) 
172 !      print *,'nq = ',nq
173       CALL GRQLIM(nq,qmin,qmax) 
174
175 !--   Get final grid definitions and grid indices of Q0, Q2C and Q2B
176
177       CALL GRGIVE(NXGRI,XMI,XMA,NQGRI,QMI,QMA)
178 !      WRITE(*,*)'NX,XL,XH,NQ,QL,QH',NXGRI,XMI,XMA,NQGRI,QMI,QMA
179       IQ0 = IQFROMQ(Q0)
180       IQC = IQFROMQ(Q2C)
181       IQB = IQFROMQ(Q2B)
182 !--   Allow for heavy weights
183
184       IF (HEAVY) THEN
185         CALL QNLSET('WTF2C',.TRUE.)
186         CALL QNLSET('WTF2B',.TRUE.)
187         CALL QNLSET('CLOWQ',.FALSE.)
188         CALL QNLSET('WTFLC',.TRUE.)
189         CALL QNLSET('WTFLB',.TRUE.)
190       ENDIF
191
192 !--   Compute weights and dump, or read in
193 !
194 !      IF (READIN) THEN 
195 !        OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
196 !     .                                  STATUS='UNKNOWN')
197 !        CALL QNREAD(24,ISTOP,IERR)
198 !      ELSE
199 !        CALL QNFILW(0,0)
200 !        IF (HEAVY) THEN
201 !          OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
202 !     .                                    STATUS='UNKNOWN')
203 !          CALL QNDUMP(24)
204 !        ENDIF
205 !      ENDIF
206
207
208       if (index(gridname,'none').eq.1) then
209          call qnfilw(0,0)
210       else
211          qnerr=-1
212          open(unit=2,status='old',file=gridname, &
213      &        form='unformatted',err=1)
214          call QNREAD(2,1,qnerr)
215  1       close(2)
216          if (qnerr.ne.0) then
217             write(*,*) 'Grid file problem: ',gridname
218             if (qnerr.lt.0) then 
219                write(*,*) 'Grid file does not exist'
220                write(*,*) 'Calculating and creating grid file'
221                call qnfilw(0,0)
222                open(unit=2,status='unknown',file=gridname, &
223      &              form='unformatted')
224                call QNDUMP(2)
225                close(2)
226             else
227                write(*,*) 'Existing grid file is inconsistent'
228                if (qnerr.eq.1) &
229      &              write(*,*) 'Defined grid different'
230                if (qnerr.eq.2) &
231      &              write(*,*) 'Heavy quark weight table different'
232                if (qnerr.eq.3) &
233      &              write(*,*) 'Charm mass different'
234                if (qnerr.eq.4) &
235      &              write(*,*) 'Bottom mass different'
236                stop
237             endif
238          endif
239       endif
240
241 !--   Apply cuts to grid
242 !--taking away the s cut at 600d0
243       CALL GRCUTS(-1D0,-1D0,-1D0,-1D0)
244
245
246
247
248 !--   Choose renormalisation and factorisation scales
249
250       CALL QNRSET('AAAR2',1D0)  ! renormalisation
251       CALL QNRSET('BBBR2',0D0)
252       CALL QNRSET('AAM2L',1D0)  ! factorisation (light)
253       CALL QNRSET('BBM2L',0D0)
254       CALL QNRSET('AAM2H',1D0)  ! factorisation (heavy)
255       CALL QNRSET('BBM2H',0D0)
256  
257 !       ZM=91.187D0
258       imem=0
259 !      print *,imem
260 ! -- only need call to listPDF here to get alphas
261       call listPDF(nset,imem,xval)
262 !      print *,xval
263         AS=XVAL(13)
264 !        AS=0.118d0
265       CALL QNRSET('ALFQ0',ZM*ZM)
266       CALL QNRSET('ALFAS',AS)
267
268 !      ZM2=ZM*ZM
269       ALPHAS=QNALFA(ZM2)
270       WRITE(*,*)'ALPHAS AT Mz2',ALPHAS
271
272 !--   Book non-singlet distributions
273
274       CALL QNBOOK(2,'UPLUS')
275       CALL QNBOOK(3,'DPLUS')
276       CALL QNBOOK(4,'SPLUS')
277       CALL QNBOOK(5,'CPLUS')
278       CALL QNBOOK(6,'BPLUS')
279       CALL QNBOOK(7,'UPVAL')
280       CALL QNBOOK(8,'DNVAL')
281  
282 !--   Book linear combinations for proton for f = 3,4,5 flavours
283
284 !--define some quark pdfs
285          CALL dVZERO(PWGT,20)        
286         PWGT(2) = 0.5
287
288         PWGT(7) = -0.5
289
290         PWGT(1) = 0.5/3.
291         CALL QNLINC(17,'UB',3,PWGT)
292         PWGT(1) = 0.5/4.
293         CALL QNLINC(17,'UB',4,PWGT)
294         PWGT(1) = 0.5/5.
295         CALL QNLINC(17,'UB',5,PWGT) 
296         CALL dVZERO(PWGT,20) 
297
298         PWGT(4) = 0.5
299         PWGT(1) = 0.5/3.
300         CALL QNLINC(18,'SB',3,PWGT)
301         PWGT(1) = 0.5/4.
302         CALL QNLINC(18,'SB',4,PWGT)
303         PWGT(1) = 0.5/5.
304         CALL QNLINC(18,'SB',5,PWGT)
305          CALL dVZERO(PWGT,20)        
306         CALL QNLINC(19,'CB',3,PWGT)
307         PWGT(5) = 0.5
308         PWGT(1) = 0.5/4.
309
310         CALL QNLINC(19,'CB',4,PWGT)
311         PWGT(1) = 0.5/5.
312         CALL QNLINC(19,'CB',5,PWGT) 
313         CALL dVZERO(PWGT,20) 
314         PWGT(3) = 0.5
315
316         PWGT(8) = -0.5
317
318         PWGT(1) = 0.5/3.
319         CALL QNLINC(20,'DB',3,PWGT)
320         PWGT(1) = 0.5/4.
321         CALL QNLINC(20,'DB',4,PWGT)
322         PWGT(1) = 0.5/5.
323         CALL QNLINC(20,'DB',5,PWGT)
324          CALL dVZERO(PWGT,20)        
325         CALL QNLINC(21,'BB',3,PWGT)
326        CALL QNLINC(21,'BB',4,PWGT)  
327         PWGT(6) = 0.5
328   
329         PWGT(1) = 0.5/5.
330         CALL QNLINC(21,'BB',5,PWGT) 
331 !---
332
333       return
334 !
335 !==========================================================================
336       entry HERApdf(nset)
337 !      print *,'calling HERApdf ',nset
338 !      ZM = 91.187d0
339 !      zm2 = zm*zm
340
341 !      ALPHAS=QNALFA(ZM2)
342       
343       
344 !      imem=mem
345       call getnmem(nset,imem)    
346       call listPDF(nset,imem,xval)
347 !      print *,xval
348 !      print *,imem,xval
349
350         FS=XVAL(14)
351         FC=XVAL(15)
352 !-- iremeber to change fs and fc with choices of Q20, and fc with =
353 !-- choices of mc.
354 !
355 !-- the free parameters are XVAL(1),2,3,4,6,9,10,11,12,13,14 which are :
356
357 !- 1,2,3,4 the parameters of the u_valence quark
358 !--ignore the names UA,Ub in the code, in the write-up these are called
359 !--x^B(1-x)^C (1+D x+E x^2) where 1=B,2=C,3=E,4=D.. you can see =
360 !--this from reading the code, but i thought it c--best to spell it out, =
361 !--normalisation comes from number sum-rule and is done in the code
362
363
364       UA=XVAL(1)
365       UB=XVAL(2)
366       UE=XVAL(3)
367       UC=XVAL(4)
368 !... 6 is the only free parameter for the d-valence: 6=C, where x^B =
369 !(1-x)^C, the B parameter is set equal to c--that of the u-valence, where =
370 !-- DA=UA means that B=1 again. normalisation comes from number sum-rule =
371 !c-- and is done in the code
372       DA=UA
373       DB=XVAL(5)
374       DE=0.0
375       DC=0.0
376 !------sea is in two bits Ubar and Dbar: parameter XVAL(9) gives an =
377 !overall sea normalisation, which is split up
378 !-- into Ubar (UBN)and Dbar(DBN) using the fractions fc and fs in the =
379 !-- code below
380       SN=XVAL(6)
381 !-- Ubar and Dbar share a common low -x slope parameter XVAL(10), so =
382 !-- both Ubar and Dbar take the form
383 !--  X^B(1-x)C, B=10 is the same for both C=11 for Ubar and C=12 =
384 !--for Dbar
385       SA=XVAL(7)
386       SB=XVAL(8)
387       SC=XVAL(9)
388       SE=0.
389 !--- make sea more complex
390       UBA=SA
391       DBA=SA
392       UBB=SB
393       DBB=SC
394       UBC=0.
395       DBC=0.
396       UBE=0.
397       DBE=0.
398 !      UBN=SN/4.
399 !      DBN=SN/4.
400        UBN=SN/2.*(1.-FS)/(2.-FS-FC)
401        DBN=SN/2.*(1.-FC)/(2.-FS-FC)
402 !-- now for a simple gluon param. x^B(1-x)^C, where B=13 and C=14, =
403 !-- and norm comes from the momentum sum-rule c--below.
404       GA=XVAL(10)
405       GB=XVAL(11)
406       GC=XVAL(12)
407       GE=0.0d0
408       
409       AS=XVAL(13)
410 !-- alphas is 0.1176 and parameter 17 is not free in the fit
411        CALL QNRSET('ALFAS',AS)
412
413 !--   Input quark distributions at Q2 = Q02 GeV2
414        UN=2./AREA(UA-1,UB,UE,UC)
415        DN=1./AREA(DA-1,DB,DE,DC)
416         GN=(1-UN*AREA(UA,UB,UE,UC)- &
417      &        DN*AREA(DA,DB,DE,DC)-2.*UBN*AREA(UBA,UBB,UBE,UBC) &
418      &       -2.*DBN*AREA(DBA,DBB,DBE,DBC))/AREA(GA,GB,GE,GC)
419
420       nxgri = nx
421       DO IX = 1,NXGRI
422         xX = XFROMIX(IX)
423         UPVAL=UN*FF_LHA2(xX,UA,UB,UE,UC)
424         DNVAL=DN*FF_LHA2(xX,DA,DB,DE,DC)
425         UBAR=UBN*FF_LHA2(xX,UBA,UBB,UBE,UBC)
426         DBAR=DBN*FF_LHA2(xX,DBA,DBB,DBE,DBC)
427         SEA=2.*(UBAR+DBAR)
428
429         SINGL=UPVAL+DNVAL+SEA  
430 !
431         GLUON=GN*FF_LHA2(xX,GA,GB,GE,GC)
432
433         UBQ=(1.-FC)*UBAR
434         CBQ=FC*UBAR
435         DBQ=(1.-FS)*DBAR
436         SBQ=FS*DBAR
437         USEA=2*UBQ
438         DSEA=2*DBQ
439         SSEA=2*SBQ
440         CSEA=2*CBQ
441 !
442         UPLUS=UPVAL+USEA-1./NACT*SINGL
443         DPLUS=DNVAL+DSEA-1./NACT*SINGL
444         SPLUS=SSEA-1./NACT*SINGL
445         CPLUS=CSEA-1./NACT*SINGL
446 !-- the fraction of charm is simply imposed- not dynamically generated
447         CALL QNPSET('GLUON',IX,IQ0,GLUON)
448         CALL QNPSET('SINGL',IX,IQ0,SINGL)
449         CALL QNPSET('UPLUS',IX,IQ0,UPLUS)
450         CALL QNPSET('DPLUS',IX,IQ0,DPLUS)
451         CALL QNPSET('SPLUS',IX,IQ0,SPLUS)
452         CALL QNPSET('CPLUS',IX,IQ0,CPLUS)
453         CALL QNPSET('UPVAL',IX,IQ0,UPVAL)
454         CALL QNPSET('DNVAL',IX,IQ0,DNVAL)
455       ENDDO
456 !--THINGS ARE FINE FOR HEAVY SO DO IT
457 !       print *, 'HEAVY is: ',heavy, 'NACT is: ',nact
458        IF (HEAVY) THEN
459
460 !--     Evolve over whole grid
461
462         CALL EVOLSG(IQ0,1,NQGRI)
463         CALL EVPLUS('UPLUS',IQ0,1,NQGRI)
464         CALL EVPLUS('DPLUS',IQ0,1,NQGRI)
465         CALL EVPLUS('SPLUS',IQ0,1,NQGRI)
466         CALL EVOLNM('UPVAL',IQ0,1,NQGRI)
467         CALL EVOLNM('DNVAL',IQ0,1,NQGRI)
468
469       ELSE
470
471 !--     Evolve gluon, singlet and valence over whole grid
472         CALL EVOLSG(IQ0,1,NQGRI)
473         CALL EVOLNM('UPVAL',IQ0,1,NQGRI)
474         CALL EVOLNM('DNVAL',IQ0,1,NQGRI)
475 !--     Be more careful with plus distributions
476         IF (NACT.EQ.3) THEN
477 !--THINGS ARE ALSO FINE IF 1Q0 IS BELOW 1QC THEN CLEARLY CSEA=0. IS OK
478 !--SITUATION CD BE 1<Q0<Q2C<Q2B ETC
479 !--GO DOWN QO TO 1 UP Q0 TO  Q2C
480           CALL EVPLUS('UPLUS',IQ0,1,IQC)
481           CALL EVPLUS('DPLUS',IQ0,1,IQC)
482           CALL EVPLUS('SPLUS',IQ0,1,IQC)
483 !--DEAL WITH CHARM THRESH
484           FACTOR=1./3.-1./4.
485           CALL QADDSI('UPLUS',IQC,FACTOR)
486           CALL QADDSI('DPLUS',IQC,FACTOR)
487           CALL QADDSI('SPLUS',IQC,FACTOR)
488           CALL QNPNUL('CPLUS')
489           FACTOR=-1/4.
490           CALL QADDSI('CPLUS',IQC,FACTOR)
491 !--GO UP Q2C TO Q2B
492           CALL EVPLUS('UPLUS',IQC,IQC,IQB)
493           CALL EVPLUS('DPLUS',IQC,IQC,IQB)
494           CALL EVPLUS('SPLUS',IQC,IQC,IQB)
495           CALL EVPLUS('CPLUS',IQC,IQC,IQB)
496
497
498         ELSEIF(NACT.EQ.4)THEN
499
500 !-NOW DO MIDDLE FOR UPLUS,DPLUS,SPLUS,CPLUS
501           CALL EVPLUS('UPLUS',IQ0,IQC,IQB)
502           CALL EVPLUS('DPLUS',IQ0,IQC,IQB)
503           CALL EVPLUS('SPLUS',IQ0,IQC,IQB)
504           CALL EVPLUS('CPLUS',IQ0,IQC,IQB)
505 !--THEN GO DOWN IQC TO 1
506              IF(IQC.GT.1)THEN
507              FACTOR=1./4.-1./3.
508              CALL QADDSI('UPLUS',IQC,FACTOR)
509              CALL QADDSI('DPLUS',IQC,FACTOR)
510              CALL QADDSI('SPLUS',IQC,FACTOR)
511              CALL EVPLUS('UPLUS',IQC,1,IQC)
512              CALL EVPLUS('DPLUS',IQC,1,IQC)
513              CALL EVPLUS('SPLUS',IQC,1,IQC)
514              ENDIF
515
516         ENDIF
517 !--THEN DEAL WITH B THRESHOLD FOR ALL4
518         FACTOR=1./4.-1./5.
519         CALL QADDSI('UPLUS',IQB,FACTOR)
520         CALL QADDSI('DPLUS',IQB,FACTOR)
521         CALL QADDSI('SPLUS',IQB,FACTOR)
522         CALL QADDSI('CPLUS',IQB,FACTOR)
523 !--THEN GO UP
524         IF(IQB.LT.NQGRI)THEN
525        CALL EVPLUS('UPLUS',IQB,IQB,NQGRI)
526         CALL EVPLUS('DPLUS',IQB,IQB,NQGRI)
527         CALL EVPLUS('SPLUS',IQB,IQB,NQGRI)
528         CALL EVPLUS('CPLUS',IQB,IQB,NQGRI)
529         ENDIF
530 !C--THEN DEAL WITH B TURNING ON AT IQB AND GO UP
531         CALL QNPNUL('BPLUS')
532         FACTOR=-1./5.
533        CALL QADDSI('BPLUS',IQB,FACTOR)
534         CALL EVPLUS('BPLUS',IQB,IQB,NQGRI)
535       ENDIF 
536       AMSR= GN*AREA(GA,GB,GE,GC)+UN*AREA(UA,UB,UE,UC) &
537      &      +  DN*AREA(DA,DB,DE,DC)+2.*UBN*AREA(UBA,UBB,UBE,UBC) &
538      &      +  2.*DBN*AREA(DBA,DBB,DBE,DBC)   
539
540       call getnset(iset)
541       call save_pdfqcd(iset)
542
543       RETURN
544   
545       END
546
547 !------------------------------------------------------------------------=
548       DOUBLE PRECISION FUNCTION AREA(A,B,E,C)
549 !------------------------------------------------------------------------=
550
551       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
552       IF (A.LE.-0.99.OR.B.LE.-0.99) THEN
553         AREA=1D6
554         RETURN
555       ENDIF
556
557       AR1=(DGAMMA_LHA(A+1)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2)
558 !      AR2=E*(DGAMMA(A+1.5)*DGAMMA(B+1))/DGAMMA(A+B+2.5)
559       AR2=E*(DGAMMA_LHA(A+3)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+4)
560       AR3=C*(DGAMMA_LHA(A+2)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+3)
561       AREA=AR1+AR2+AR3
562       IF (AREA.LE.1D-6) AREA=1D-6
563
564       RETURN
565       END
566
567 !------------------------------------------------------------------------=
568       DOUBLE PRECISION FUNCTION FF_LHA2(X,A,B,E,C)
569 !------------------------------------------------------------------------=
570
571       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
572 !      FF=X**A*(1D0-X)**B*(1+E*SQRT(X)+C*X)
573       FF_LHA2=X**A*(1D0-X)**B*(1+C*X + E*X*X)
574       RETURN
575       END
576
577 !------------------------------------------------------------------------ 
578       
579
580
581