-!****************************************************************************\r
-!\r
-! EPS09.f\r
-!\r
-! An interface for the scale dependent nuclear modifications\r
-! R_f^A(x,Q) = f_A(x,Q)/f_p(x,Q) \r
-! where f_A is the distribution of the parton flavour f for a PROTON in a\r
-! nucleus A, and f_p is the corresponding parton distribution in the \r
-! free proton.\r
-! \r
-! When using this interface, please refer to:\r
-! \r
-! K.J. Eskola, H. Paukkunen and C.A. Salgado,\r
-! "EPS09 - a New Generation of NLO and LO Nuclear Parton Distribution Functions,"\r
-! Published as JHEP04(2009) 065.\r
-! Eprint: arXiv:0902.4154 [hep-ph].\r
-!\r
-! Questions & comments to:\r
-! hannu.paukkunen@phys.jyu.fi\r
-! kari.eskola@phys.jyu.fi\r
-! carlos.salgado@usc.es\r
-! \r
-! ***************************************************************************\r
-! Instructions:\r
-!\r
-! For given input values of\r
-!\r
-! order: 1=LO, 2=NLO ; integer\r
-! pset : 1...31 ; integer\r
-! 1 = central fit\r
-! 2,3 = error sets S{+1}, S{-1}\r
-! 4,5 = error sets S{+2}, S{-2}\r
-! ... ...\r
-! 30,31 = error sets {S+15}, {S-15}\r
-! A : atomic number ; integer\r
-! x : Bjorken-x ; double precision\r
-! Q : scale in GeV ; double precision\r
-!\r
-! the command \r
-!\r
-! Call EPS09(order, pset, A, x, Q, ruv, rdv, ru, rd, rs, rc, rb, rg)\r
-!\r
-! returns the bound proton nuclear corrections R_f^A(x,Q)\r
-! (in double precision) for\r
-! \r
-! ruv = up valence\r
-! rdv = down valence\r
-! ru = up sea\r
-! rd = down sea\r
-! rs = strange\r
-! rc = charm\r
-! rb = bottom\r
-! rg = gluons\r
-!\r
-! The nuclear corrections for bound neutrons can be obtained\r
-! by the isospin symmetry, e.g. the total up quark distribution\r
-! per nucleon in a nucleus A with Z protons is\r
-!\r
-! u_A(x,Q) = Z/A * [ruv*uV_p(x,Q) + ru*uSea_p(x,Q)] +\r
-! (A-Z)/A * [rdv*dV_p(x,Q) + rd*dSea_p(x,Q)]\r
-!\r
-! Note that the parametrization should only be applied at the\r
-! kinematical domain\r
-!\r
-! 1e-6 <= x <= 1\r
-! 1.3 <= Q <= 1000 GeV.\r
-!\r
-! No warning message is displayed if these limits are\r
-! exceeded, and outside these boundaries the modifications\r
-! are frozen to the boundary values, i.e\r
-!\r
-! for Q > 1000, the modifications at Q=1000 are returned,\r
-! for Q < 1.3, the modifications at Q=1.3 are returned,\r
-! for x > 1, the modifications at x=1 are returned\r
-! for x < 1e-6, the modifications at x=1e-6 are returned,\r
-!\r
-! The data used by the program for required order\r
-! and atomic number A, are stored in separate files\r
-!\r
-! LO : EPS09LOR_A\r
-! NLO: EPS09NLOR_A\r
-!\r
-! which must be located in the working directory.\r
-!\r
-! The error bands for absolute cross-sections and for\r
-! their nuclear ratios should be computed as explained\r
-! in Secs. 2.5 and 4 of arXiv:0902.4154 [hep-ph]. For\r
-! the absolute cross sections, both the errors in the\r
-! free-proton PDFs f_p(x,Q) and the errors in\r
-! the modifications R_f^A(x,Q) should be accounted for.\r
-! For the nuclear ratios, it is sufficient to account only\r
-! for the errors in the modifications R_f^A(x,Q).\r
-!\r
-! *********************************************************\r
-! *********************************************************\r
-\r
-\r
-\r
- Subroutine EPS09(order, pset, AAA, xxx, QQQ,ruv, rdv, ru, rd, rs, rc, rb, rg)\r
-\r
- Implicit none\r
- Double precision :: ruv, rdv, ru, rd, rs, rc, rb, rg, QQQ, xxx\r
- Double precision :: LSTEP, x, Q, Q2, allvalues(1:31,1:8,0:50,0:50)\r
- Double precision :: x_i=0.000001, arg(4), fu(4), res, fg(3)\r
- Double precision :: result(9), dummy\r
- Double precision :: realQ, Q2min=1.69, Q2max=1000000.0, Qsteps=50.0\r
- Double precision :: n_x, zero=0.0\r
-\r
- Character (Len=50) filenimi\r
-\r
- Integer :: xlinsteps=25, xlogsteps=25, startline, lineno\r
- Integer :: k, p, t, Qpoint, xpoint, pset, iovar\r
- Integer :: setnumber,j, A, openchannel, order, AAA\r
- Integer :: psetlast = -10, Alast = -10, orderlast = -10\r
- \r
- character*512 dirpath,setpath\r
-\r
- save Alast\r
- save psetlast\r
- save orderlast\r
- save allvalues\r
-\r
-! *********************************************\r
-! Stop if the set specifications are wrong ones\r
-! *********************************************\r
-\r
- If (order .NE. 1 .and. order .NE. 2) then\r
- Write(*,*) 'Wrong order!'\r
- Write(*,*) 'LO : order = 1'\r
- Write(*,*) 'NLO: order = 2'\r
- Stop\r
- End If\r
-\r
- If (pset .LT. 1 .or. pset .GT. 31) then\r
- Write(*,*) 'Wrong set!'\r
- Write(*,*) 'Central set: pset = 1'\r
- Write(*,*) 'Error sets : pset = 2...31'\r
- Stop\r
- End If\r
-\r
-! ********************************\r
-! Make sure not to change any\r
-! specifications given by the user\r
-! ********************************\r
-\r
- A = AAA\r
- x = xxx\r
- Q = QQQ\r
- Q2 = Q*Q \r
-\r
-! *******************************\r
-! Freeze x if it's < 10E-6 or > 1\r
-! *******************************\r
-\r
- If (x .LT. x_i) Then\r
- x = x_i\r
- End If\r
- If (x .GT. 1) Then\r
- x = 1.0\r
- End If\r
-\r
-! ************************************\r
-! Freeze Q^2 if it's < 1.69 or > 10E+6\r
-! ************************************\r
-\r
- If (Q2 .LT. Q2min) Then\r
- Q2 = Q2min\r
- End If\r
- If (Q2 .GT. Q2max) Then\r
- Q2 = Q2max\r
- End If\r
-\r
-! If the set specifications have been changed, read the tables again\r
-\r
- If (A .NE. Alast .or. order .NE. orderlast) Then\r
-\r
-! Write(*,*) 'Set changed!'\r
-\r
-! Read the table1\r
-\r
- If (order .EQ. 1) then\r
-\r
- If (A < 10) Then\r
- Write(filenimi,'("EPS09LOR_", I1)'), A\r
- Else If (A < 100) Then\r
- Write(filenimi,'("EPS09LOR_", I2)'), A\r
- Else If (A < 1000) Then\r
- Write(filenimi,'("EPS09LOR_", I3)'), A\r
- End If\r
-\r
- Else\r
-\r
- If (A < 10) Then\r
- Write(filenimi,'("EPS09NLOR_", I1)'), A\r
- Else If (A < 100) Then\r
- Write(filenimi,'("EPS09NLOR_", I2)'), A\r
- Else If (A < 1000) Then\r
- Write(filenimi,'("EPS09NLOR_", I3)'), A\r
- End If\r
-\r
- End If\r
-\r
- call getdirpath(dirpath) \r
- setpath=dirpath(:LEN_TRIM(dirpath))//"/"//filenimi(:LEN_TRIM(filenimi))//".LHgrid"\r
- print *,setpath(:LEN_TRIM(setpath))\r
- OPEN (11, file = setpath(:LEN_TRIM(setpath)), status='OLD', IOSTAT=iovar)\r
-\r
- If (iovar .NE. 0) Then\r
- Write(*,*) 'Missing file: ',filenimi\r
- stop\r
- End If\r
-\r
- Do setnumber = 1, 31\r
-\r
- Do k = 0,50\r
-\r
- Read(11,*) dummy\r
-\r
- Do t = 0,50\r
-\r
- Read(11,*) (allvalues(setnumber,p,k,t), p=1,8)\r
-\r
- End Do\r
- End Do\r
-\r
- End Do\r
-\r
- Close(11)\r
-\r
- psetlast = pset\r
- Alast = A\r
- orderlast = order\r
-\r
- End If\r
-\r
-! Find out the position in the loglog Q^2-grid\r
-\r
- realQ = Qsteps * (log(log(Q2)/log(Q2min)))/(log(log(Q2max)/log(Q2min)))\r
- Qpoint = Aint(realQ)\r
-\r
- If (Qpoint .LE. 0) Then\r
- Qpoint = 1\r
- End If\r
- If (Qpoint .GE. Anint(Qsteps)-1) Then\r
- Qpoint = Anint(Qsteps)-1\r
- End If\r
-\r
- LSTEP = (1.0/(xlogsteps)) * LOG(0.1/x_i)\r
-\r
-! *********************\r
-! Interpolate the grids \r
-! *********************\r
-\r
- Do t=1,8\r
-\r
-! Find the position in the x-grid\r
-\r
- If (x .LE. 0.1) then\r
- n_x = ((1.0/LSTEP) * Log(x/x_i))\r
- xpoint = Aint(n_x)\r
- Else\r
- n_x = ((x-0.1)*xlinsteps/(1.0-0.1) + xlogsteps)\r
- xpoint = Aint(n_x)\r
- End If\r
-\r
- If (xpoint .LE. 0) Then\r
- xpoint = 1\r
- End If\r
-\r
- If (t .EQ. 1 .or. t .EQ. 2) Then\r
- If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then\r
- xpoint = (xlinsteps+xlogsteps)-4\r
- End If\r
- End If\r
-\r
- If (t .EQ. 3 .or. t .EQ. 4 .or. t .EQ. 5 .or. t .EQ. 6 .or. t .EQ. 7) Then\r
- If (xpoint .GE. (xlinsteps+xlogsteps)-7) Then\r
- xpoint = (xlinsteps+xlogsteps)-7\r
- End If\r
- End If\r
- \r
- If (t .EQ. 8) Then\r
- If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then\r
- xpoint = (xlinsteps+xlogsteps)-4\r
- End If\r
- End If\r
-\r
- Do k = 1, 4\r
- If (xpoint-2+k .LT. xlogsteps) Then\r
- arg(k) = (x_i) * exp(LSTEP * (xpoint-2+k))\r
- Else\r
- arg(k) = 0.1 + (xpoint-2+k-xlogsteps) * (1-0.1)/xlinsteps\r
- End If\r
- End Do\r
- \r
- Do j=1,3\r
-\r
- fu(1) = allvalues(pset,t,Qpoint-2+j,xpoint-1)\r
- fu(2) = allvalues(pset,t,Qpoint-2+j,xpoint)\r
- fu(3) = allvalues(pset,t,Qpoint-2+j,xpoint+1)\r
- fu(4) = allvalues(pset,t,Qpoint-2+j,xpoint+2)\r
- Call luovi(fu,arg,4,x,res)\r
- fg(j) = res\r
-\r
-! *****************************************\r
-! *****************************************\r
-\r
- End Do\r
-\r
- arg(1) = Qpoint-1\r
- arg(2) = Qpoint\r
- arg(3) = Qpoint+1\r
-\r
- Call luovi(fg,arg,3,realQ,res)\r
- \r
- result(t) = res\r
-\r
- End Do\r
-\r
- ruv = max(result(1),zero)\r
- rdv = max(result(2),zero)\r
- ru = max(result(3),zero)\r
- rd = max(result(4),zero)\r
- rs = max(result(5),zero)\r
- rc = max(result(6),zero)\r
- rb = max(result(7),zero)\r
- rg = max(result(8),zero)\r
-\r
-200 Continue\r
-\r
- End Subroutine EPS09\r
-\r
-! ********************************\r
-! Modified version of Cern Library\r
-! interpolation routine E100\r
-! ********************************\r
-\r
- SUBROUTINE luovi(F,ARG,MMM,Z,SUM)\r
-\r
- Implicit none\r
- INTEGER :: MMM\r
- Double precision :: F(MMM), ARG(MMM), COF(MMM), SUM, Z\r
- INTEGER :: M, MM, I, J, JNDEX, INDEX\r
-\r
- MM = MIN(MMM, 20)\r
- M = MM - 1\r
- DO 1780 I= 1, MM\r
- COF(I) = F(I)\r
- 1780 Continue\r
- DO 1800 I= 1, M\r
- DO 1790 J= I, M\r
- JNDEX = MM - J\r
- INDEX = JNDEX + I\r
- COF(INDEX) = (COF(INDEX)-COF(INDEX-1))/(ARG(INDEX)-ARG(JNDEX))\r
- 1790 CONTINUE\r
- 1800 CONTINUE\r
- SUM = COF(MM)\r
- DO 1810 I= 1, M\r
- INDEX = MM - I\r
- SUM = (Z-ARG(INDEX))*SUM + COF(INDEX)\r
- 1810 CONTINUE\r
-\r
- End SUBROUTINE luovi\r
-\r
+!****************************************************************************
+!
+! EPS09.f
+!
+! An interface for the scale dependent nuclear modifications
+! R_f^A(x,Q) = f_A(x,Q)/f_p(x,Q)
+! where f_A is the distribution of the parton flavour f for a PROTON in a
+! nucleus A, and f_p is the corresponding parton distribution in the
+! free proton.
+!
+! When using this interface, please refer to:
+!
+! K.J. Eskola, H. Paukkunen and C.A. Salgado,
+! "EPS09 - a New Generation of NLO and LO Nuclear Parton Distribution Functions,"
+! Published as JHEP04(2009) 065.
+! Eprint: arXiv:0902.4154 [hep-ph].
+!
+! Questions & comments to:
+! hannu.paukkunen@phys.jyu.fi
+! kari.eskola@phys.jyu.fi
+! carlos.salgado@usc.es
+!
+! ***************************************************************************
+! Instructions:
+!
+! For given input values of
+!
+! order: 1=LO, 2=NLO ; integer
+! pset : 1...31 ; integer
+! 1 = central fit
+! 2,3 = error sets S{+1}, S{-1}
+! 4,5 = error sets S{+2}, S{-2}
+! ... ...
+! 30,31 = error sets {S+15}, {S-15}
+! A : atomic number ; integer
+! x : Bjorken-x ; double precision
+! Q : scale in GeV ; double precision
+!
+! the command
+!
+! Call EPS09(order, pset, A, x, Q, ruv, rdv, ru, rd, rs, rc, rb, rg)
+!
+! returns the bound proton nuclear corrections R_f^A(x,Q)
+! (in double precision) for
+!
+! ruv = up valence
+! rdv = down valence
+! ru = up sea
+! rd = down sea
+! rs = strange
+! rc = charm
+! rb = bottom
+! rg = gluons
+!
+! The nuclear corrections for bound neutrons can be obtained
+! by the isospin symmetry, e.g. the total up quark distribution
+! per nucleon in a nucleus A with Z protons is
+!
+! u_A(x,Q) = Z/A * [ruv*uV_p(x,Q) + ru*uSea_p(x,Q)] +
+! (A-Z)/A * [rdv*dV_p(x,Q) + rd*dSea_p(x,Q)]
+!
+! Note that the parametrization should only be applied at the
+! kinematical domain
+!
+! 1e-6 <= x <= 1
+! 1.3 <= Q <= 1000 GeV.
+!
+! No warning message is displayed if these limits are
+! exceeded, and outside these boundaries the modifications
+! are frozen to the boundary values, i.e
+!
+! for Q > 1000, the modifications at Q=1000 are returned,
+! for Q < 1.3, the modifications at Q=1.3 are returned,
+! for x > 1, the modifications at x=1 are returned
+! for x < 1e-6, the modifications at x=1e-6 are returned,
+!
+! The data used by the program for required order
+! and atomic number A, are stored in separate files
+!
+! LO : EPS09LOR_A
+! NLO: EPS09NLOR_A
+!
+! which must be located in the working directory.
+!
+! The error bands for absolute cross-sections and for
+! their nuclear ratios should be computed as explained
+! in Secs. 2.5 and 4 of arXiv:0902.4154 [hep-ph]. For
+! the absolute cross sections, both the errors in the
+! free-proton PDFs f_p(x,Q) and the errors in
+! the modifications R_f^A(x,Q) should be accounted for.
+! For the nuclear ratios, it is sufficient to account only
+! for the errors in the modifications R_f^A(x,Q).
+!
+! *********************************************************
+! *********************************************************
+
+
+
+ Subroutine EPS09(order, pset, AAA, xxx, QQQ,ruv, rdv, ru, rd, rs, rc, rb, rg)
+
+ Implicit none
+ Double precision :: ruv, rdv, ru, rd, rs, rc, rb, rg, QQQ, xxx
+ Double precision :: LSTEP, x, Q, Q2, allvalues(1:31,1:8,0:50,0:50)
+ Double precision :: x_i=0.000001, arg(4), fu(4), res, fg(3)
+ Double precision :: result(9), dummy
+ Double precision :: realQ, Q2min=1.69, Q2max=1000000.0, Qsteps=50.0
+ Double precision :: n_x, zero=0.0
+
+ Character (Len=50) filenimi
+
+ Integer :: xlinsteps=25, xlogsteps=25, startline, lineno
+ Integer :: k, p, t, Qpoint, xpoint, pset, iovar
+ Integer :: setnumber,j, A, openchannel, order, AAA
+ Integer :: psetlast = -10, Alast = -10, orderlast = -10
+
+ character*512 dirpath,setpath
+
+ save Alast
+ save psetlast
+ save orderlast
+ save allvalues
+
+! *********************************************
+! Stop if the set specifications are wrong ones
+! *********************************************
+
+ If (order .NE. 1 .and. order .NE. 2) then
+ Write(*,*) 'Wrong order!'
+ Write(*,*) 'LO : order = 1'
+ Write(*,*) 'NLO: order = 2'
+ Stop
+ End If
+
+ If (pset .LT. 1 .or. pset .GT. 31) then
+ Write(*,*) 'Wrong set!'
+ Write(*,*) 'Central set: pset = 1'
+ Write(*,*) 'Error sets : pset = 2...31'
+ Stop
+ End If
+
+! ********************************
+! Make sure not to change any
+! specifications given by the user
+! ********************************
+
+ A = AAA
+ x = xxx
+ Q = QQQ
+ Q2 = Q*Q
+
+! *******************************
+! Freeze x if it's < 10E-6 or > 1
+! *******************************
+
+ If (x .LT. x_i) Then
+ x = x_i
+ End If
+ If (x .GT. 1) Then
+ x = 1.0
+ End If
+
+! ************************************
+! Freeze Q^2 if it's < 1.69 or > 10E+6
+! ************************************
+
+ If (Q2 .LT. Q2min) Then
+ Q2 = Q2min
+ End If
+ If (Q2 .GT. Q2max) Then
+ Q2 = Q2max
+ End If
+
+! If the set specifications have been changed, read the tables again
+
+ If (A .NE. Alast .or. order .NE. orderlast) Then
+
+! Write(*,*) 'Set changed!'
+
+! Read the table1
+
+ If (order .EQ. 1) then
+
+ If (A < 10) Then
+ Write(filenimi,'("EPS09LOR_", I1)'), A
+ Else If (A < 100) Then
+ Write(filenimi,'("EPS09LOR_", I2)'), A
+ Else If (A < 1000) Then
+ Write(filenimi,'("EPS09LOR_", I3)'), A
+ End If
+
+ Else
+
+ If (A < 10) Then
+ Write(filenimi,'("EPS09NLOR_", I1)'), A
+ Else If (A < 100) Then
+ Write(filenimi,'("EPS09NLOR_", I2)'), A
+ Else If (A < 1000) Then
+ Write(filenimi,'("EPS09NLOR_", I3)'), A
+ End If
+
+ End If
+
+ call getdirpath(dirpath)
+ setpath=dirpath(:LEN_TRIM(dirpath))//"/"//filenimi(:LEN_TRIM(filenimi))//".LHgrid"
+ print *,setpath(:LEN_TRIM(setpath))
+ OPEN (11, file = setpath(:LEN_TRIM(setpath)), status='OLD', IOSTAT=iovar)
+
+ If (iovar .NE. 0) Then
+ Write(*,*) 'Missing file: ',filenimi
+ stop
+ End If
+
+ Do setnumber = 1, 31
+
+ Do k = 0,50
+
+ Read(11,*) dummy
+
+ Do t = 0,50
+
+ Read(11,*) (allvalues(setnumber,p,k,t), p=1,8)
+
+ End Do
+ End Do
+
+ End Do
+
+ Close(11)
+
+ psetlast = pset
+ Alast = A
+ orderlast = order
+
+ End If
+
+! Find out the position in the loglog Q^2-grid
+
+ realQ = Qsteps * (log(log(Q2)/log(Q2min)))/(log(log(Q2max)/log(Q2min)))
+ Qpoint = Aint(realQ)
+
+ If (Qpoint .LE. 0) Then
+ Qpoint = 1
+ End If
+ If (Qpoint .GE. Anint(Qsteps)-1) Then
+ Qpoint = Anint(Qsteps)-1
+ End If
+
+ LSTEP = (1.0/(xlogsteps)) * LOG(0.1/x_i)
+
+! *********************
+! Interpolate the grids
+! *********************
+
+ Do t=1,8
+
+! Find the position in the x-grid
+
+ If (x .LE. 0.1) then
+ n_x = ((1.0/LSTEP) * Log(x/x_i))
+ xpoint = Aint(n_x)
+ Else
+ n_x = ((x-0.1)*xlinsteps/(1.0-0.1) + xlogsteps)
+ xpoint = Aint(n_x)
+ End If
+
+ If (xpoint .LE. 0) Then
+ xpoint = 1
+ End If
+
+ If (t .EQ. 1 .or. t .EQ. 2) Then
+ If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then
+ xpoint = (xlinsteps+xlogsteps)-4
+ End If
+ End If
+
+ If (t .EQ. 3 .or. t .EQ. 4 .or. t .EQ. 5 .or. t .EQ. 6 .or. t .EQ. 7) Then
+ If (xpoint .GE. (xlinsteps+xlogsteps)-7) Then
+ xpoint = (xlinsteps+xlogsteps)-7
+ End If
+ End If
+
+ If (t .EQ. 8) Then
+ If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then
+ xpoint = (xlinsteps+xlogsteps)-4
+ End If
+ End If
+
+ Do k = 1, 4
+ If (xpoint-2+k .LT. xlogsteps) Then
+ arg(k) = (x_i) * exp(LSTEP * (xpoint-2+k))
+ Else
+ arg(k) = 0.1 + (xpoint-2+k-xlogsteps) * (1-0.1)/xlinsteps
+ End If
+ End Do
+
+ Do j=1,3
+
+ fu(1) = allvalues(pset,t,Qpoint-2+j,xpoint-1)
+ fu(2) = allvalues(pset,t,Qpoint-2+j,xpoint)
+ fu(3) = allvalues(pset,t,Qpoint-2+j,xpoint+1)
+ fu(4) = allvalues(pset,t,Qpoint-2+j,xpoint+2)
+ Call luovi(fu,arg,4,x,res)
+ fg(j) = res
+
+! *****************************************
+! *****************************************
+
+ End Do
+
+ arg(1) = Qpoint-1
+ arg(2) = Qpoint
+ arg(3) = Qpoint+1
+
+ Call luovi(fg,arg,3,realQ,res)
+
+ result(t) = res
+
+ End Do
+
+ ruv = max(result(1),zero)
+ rdv = max(result(2),zero)
+ ru = max(result(3),zero)
+ rd = max(result(4),zero)
+ rs = max(result(5),zero)
+ rc = max(result(6),zero)
+ rb = max(result(7),zero)
+ rg = max(result(8),zero)
+
+200 Continue
+
+ End Subroutine EPS09
+
+! ********************************
+! Modified version of Cern Library
+! interpolation routine E100
+! ********************************
+
+ SUBROUTINE luovi(F,ARG,MMM,Z,SUM)
+
+ Implicit none
+ INTEGER :: MMM
+ Double precision :: F(MMM), ARG(MMM), COF(MMM), SUM, Z
+ INTEGER :: M, MM, I, J, JNDEX, INDEX
+
+ MM = MIN(MMM, 20)
+ M = MM - 1
+ DO 1780 I= 1, MM
+ COF(I) = F(I)
+ 1780 Continue
+ DO 1800 I= 1, M
+ DO 1790 J= I, M
+ JNDEX = MM - J
+ INDEX = JNDEX + I
+ COF(INDEX) = (COF(INDEX)-COF(INDEX-1))/(ARG(INDEX)-ARG(JNDEX))
+ 1790 CONTINUE
+ 1800 CONTINUE
+ SUM = COF(MM)
+ DO 1810 I= 1, M
+ INDEX = MM - I
+ SUM = (Z-ARG(INDEX))*SUM + COF(INDEX)
+ 1810 CONTINUE
+
+ End SUBROUTINE luovi
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff