ltpth/HDECAY
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
4884213286e1edbdf3b2bb9fe30095b29e336dff
HDECAY/zzz/elw.f
Basil Bruhn 4884213286 initial 
Signed-off-by: Basil Bruhn <basil.bruhn@psi.ch>
2025-11-17 10:04:29 +01:00

2225 lines
104 KiB
Fortran
Raw Blame History

*******************************************************************************
* *
* EWgint *
* *
* Interpolating the EW grid in gg --> H *
* *
* November 2009, G. Passarino, C. Sturm, S, Uccirati *
* *
* *
* S.~Actis, G.~Passarino, C.~Sturm and S.~Uccirati, *
* Nucl.\ Phys.\ B {\bf 811} (2009) 182 *
* [arXiv:0809.3667 [hep-ph]]. *
* *
* S.~Actis, G.~Passarino, C.~Sturm and S.~Uccirati, *
* Phys.\ Lett.\ B {\bf 670} (2008) 12 *
* [arXiv:0809.1301 [hep-ph]] *
* *
* S.~Actis, G.~Passarino, C.~Sturm and S.~Uccirati, *
* Phys.\ Lett.\ B {\bf 669} (2008) 62 *
* [arXiv:0809.1302 [hep-ph]]. *
* *
* *
* Email: giampiero@to.infn.it *
* *
*******************************************************************************
*
double precision function glgl_elw(mtop,mh)
*
IMPLICIT NONE
*
INTEGER i,top,gdim
REAL*8 u,ui,value
REAL*8 vp,vpp,vppp
REAL*8 mtop,mh
REAL*8 x,x0,x1,x2,y0,y1,y2,a0,a1,a2
REAL*8 value0,value1,value2
REAL*8 bl(154),cl(154),dl(154)
REAL*8 bc(151),cc(151),dc(151)
REAL*8 bu(152),cu(152),du(152)
*
* u value of M_H at which the spline is to be evaluated
* top= -1,0,1 lower, central, upper value for m_top
*
u = mh
top = -1
gdim= 154
CALL FMMsplineSingle(bl,cl,dl,top,gdim)
CALL Seval3Single(u,bl,cl,dl,top,gdim,value0,vp,vpp,vppp)
top = 0
gdim= 151
CALL FMMsplineSingle(bc,cc,dc,top,gdim)
CALL Seval3Single(u,bc,cc,dc,top,gdim,value1,vp,vpp,vppp)
top = 1
gdim= 152
CALL FMMsplineSingle(bu,cu,du,top,gdim)
CALL Seval3Single(u,bu,cu,du,top,gdim,value2,vp,vpp,vppp)
x = mtop
x0 = 171.06d0
x1 = 172.64d0
x2 = 174.22d0
y0 = value0
y1 = value1
y2 = value2
A0=(X-X1)*(X-X2)/(X0-X1)/(X0-X2)
A1=(X-X0)*(X-X2)/(X1-X0)/(X1-X2)
A2=(X-X0)*(X-X1)/(X2-X0)/(X2-X1)
glgl_elw=(A0*Y0+A1*Y1+A2*Y2)/100.d0
RETURN
END
*
*-----------------------------------------------------------------------
*
c CONTAINS
*
SUBROUTINE FMMsplineSingle(b,c,d,top,gdim)
*
* ---------------------------------------------------------------------------
*
* PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline
* so that the interpolated value is given by
* s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3
* when x(k) <= x <= x(k+1)
* The end conditions match the third derivatives of the interpolated curve to
* the third derivatives of the unique polynomials thru the first four and
* last four points.
* Use Seval or Seval3 to evaluate the spline.
*
INTEGER k,n,i,top,gdim,l
*
REAL*8 xl(154),yl(154)
REAL*8 xc(151),yc(151)
REAL*8 xu(152),yu(152)
REAL*8 x(154),y(154)
*
REAL*8 b(gdim)
* linear coeff
*
REAL*8 c(gdim)
* quadratic coeff.
*
REAL*8 d(gdim)
* cubic coeff.
*
REAL*8 t
REAL*8 ZERO, TWO, THREE
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
*
* The grid
*
DATA (xl(i),i=1,154)/100.0d0,110.0d0,120.0d0,130.0d0,
# 140.0d0,145.0d0,150.0d0,151.0d0,
# 152.0d0,153.0d0,154.0d0,155.0d0,
# 156.0d0,157.0d0,158.0d0,159.0d0,
# 160.0d0,161.0d0,162.0d0,163.0d0,
# 164.0d0,165.0d0,166.0d0,
# 167.0d0,168.0d0,169.0d0,170.0d0,
# 171.0d0,172.0d0,173.0d0,174.0d0,
# 175.0d0,176.0d0,177.0d0,178.0d0,
# 179.0d0,180.0d0,181.0d0,182.0d0,
# 183.0d0,184.0d0,185.0d0,186.0d0,
# 187.0d0,188.0d0,189.0d0,190.0d0,
# 195.0d0,197.0d0,200.0d0,210.0d0,
# 220.0d0,230.0d0,240.0d0,250.0d0,
# 260.0d0,270.0d0,280.0d0,290.0d0,
# 300.0d0,310.0d0,320.0d0,330.0d0,
# 335.0d0,340.0d0,341.0d0,342.0d0,
# 343.0d0,344.0d0,345.0d0,345.3d0,
# 345.5d0,345.8d0,346.0d0,347.0d0,
# 348.0d0,349.0d0,349.5d0,349.75d0,
# 350.0d0,351.0d0,352.0d0,353.0d0,
# 354.0d0,355.0d0,356.0d0,357.0d0,
# 358.0d0,359.0d0,360.0d0,370.0d0,
# 380.0d0,390.0d0,400.0d0,410.0d0,
# 420.0d0,430.0d0,440.0d0,450.0d0,
# 460.0d0,470.0d0,480.0d0,490.0d0,
# 500.0d0,510.0d0,520.0d0,530.0d0,
# 540.0d0,550.0d0,560.0d0,570.0d0,
# 580.0d0,590.0d0,600.0d0,610.0d0,
# 620.0d0,630.0d0,640.0d0,650.0d0,
# 660.0d0,670.0d0,680.0d0,690.0d0,
# 700.0d0,710.0d0,720.0d0,730.0d0,
# 740.0d0,750.0d0,760.0d0,770.0d0,
# 780.0d0,790.0d0,800.0d0,810.0d0,
# 820.0d0,830.0d0,840.0d0,850.0d0,
# 860.0d0,870.0d0,880.0d0,890.0d0,
# 900.0d0,910.0d0,920.0d0,930.0d0,
# 940.0d0,950.0d0,960.0d0,970.0d0,
# 980.0d0,990.0d0,1000.0d0/
*
DATA (yl(i),i=1,154)/4.179467154d0,4.542088751d0,
# 4.938014960d0,5.335142457d0,
# 5.707514034d0,5.868498547d0,5.978310949d0,5.986920781d0,
# 5.987460634d0,5.977179301d0,5.951667414d0,5.904562092d0,
# 5.826081147d0,5.700679828d0,5.506218392d0,5.220567269d0,
# 4.861752146d0,4.527674442d0,4.222683291d0,3.880374785d0,
# 3.529354583d0,3.204116062d0,2.915522583d0,
# 2.662502577d0,2.440241875d0,2.243618198d0,2.067597988d0,
# 1.907471856d0,1.759641285d0,1.620735593d0,1.487565904d0,
# 1.355817302d0,1.222813267d0,1.081630216d0,0.926313642d0,
# 0.747981132d0,0.538696787d0,0.306541839d0,0.085775708d0,
# -0.092817756d0,-0.259991178d0,-0.443225249d0,-0.632350007d0,
# -0.811227830d0,-0.975263093d0,-1.126478218d0,-1.262721067d0,
# -1.757119775d0,-1.895339688d0,-2.057961996d0,-2.370870160d0,
# -2.488446977d0,-2.490894901d0,-2.433191749d0,-2.359142379d0,
# -2.260728908d0,-2.151155918d0,-2.036327859d0,-1.933794725d0,
# -1.862582453d0,-1.802070746d0,-1.812386184d0,-1.938240382d0,
# -2.108006605d0,-2.488306576d0,-2.655000325d0,-3.016952793d0,
# -3.680237527d0,-3.870995189d0,-3.959387796d0,-3.943637421d0,
# -3.982342165d0,-3.967290244d0,-4.002598626d0,-4.008430095d0,
# -4.010824003d0,-4.021419122d0,-4.018822218d0,-4.040141549d0,
# -4.024343762d0,-4.015666188d0,-3.965077563d0,-3.968563552d0,
# -3.963667776d0,-3.942755435d0,-3.932719674d0,-3.874191897d0,
# -3.876081328d0,-3.833375576d0,-3.796241844d0,-3.406242887d0,
# -3.088048772d0,-2.738617943d0,-2.446609831d0,-2.131731318d0,
# -1.878917579d0,-1.589864293d0,-1.358374992d0,-1.103676110d0,
# -0.884448436d0,-0.679633948d0,-0.472673942d0,-0.281736254d0,
# -0.108375079d0,0.093529154d0,0.257607461d0,0.430152846d0,
# 0.578396867d0,0.763320409d0,0.934522954d0,1.075378845d0,
# 1.249326065d0,1.426785290d0,1.535796936d0,1.715981764d0,
# 1.871449481d0,2.020126942d0,2.154595978d0,2.322995662d0,
# 2.456079279d0,2.608032219d0,2.749902433d0,2.905996823d0,
# 3.064612486d0,3.205141310d0,3.344347828d0,3.506879165d0,
# 3.660117972d0,3.783439131d0,3.947284161d0,4.097179155d0,
# 4.248755764d0,4.406763712d0,4.563193913d0,4.705927695d0,
# 4.853626937d0,5.015246327d0,5.167579841d0,5.341215028d0,
# 5.512112004d0,5.650306940d0,5.805087521d0,5.980800439d0,
# 6.125262600d0,6.272280377d0,6.477594572d0,6.625428154d0,
# 6.810107710d0,6.954753946d0,7.131768272d0,7.304668734d0,
# 7.458544663d0,7.653458644d0,7.821154640d0/
*
DATA (xc(i),i=1,151)/100.0d0,110.0d0,120.0d0,130.0d0,140.0d0,
# 145.0d0,150.0d0,151.0d0,152.0d0,153.0d0,154.0d0,155.0d0,
# 156.0d0,157.0d0,158.0d0,159.0d0,160.0d0,161.0d0,162.0d0,
# 163.0d0,164.0d0,165.0d0,166.0d0,167.0d0,168.0d0,169.0d0,
# 170.0d0,171.0d0,172.0d0,173.0d0,174.0d0,175.0d0,176.0d0,
# 177.0d0,178.0d0,179.0d0,180.0d0,181.0d0,182.0d0,183.0d0,
# 184.0d0,185.0d0,186.0d0,187.0d0,188.0d0,
# 189.0d0,190.0d0,195.0d0,197.0d0,200.0d0,
# 210.0d0,220.0d0,230.0d0,240.0d0,250.0d0,
# 260.0d0,270.0d0,280.0d0,290.0d0,300.0d0,
# 320.0d0,330.0d0,335.0d0,340.0d0,341.0d0,
# 342.0d0,343.0d0,344.0d0,345.0d0,345.3d0,
# 345.5d0,345.8d0,346.0d0,347.0d0,348.0d0,
# 349.0d0,350.0d0,351.0d0,352.0d0,353.0d0,
# 354.0d0,355.0d0,356.0d0,357.0d0,358.0d0,
# 359.0d0,360.0d0,370.0d0,380.0d0,390.0d0,
# 400.0d0,410.0d0,420.0d0,430.0d0,440.0d0,
# 450.0d0,460.0d0,470.0d0,480.0d0,490.0d0,
# 500.0d0,510.0d0,520.0d0,530.0d0,540.0d0,
# 550.0d0,560.0d0,570.0d0,580.0d0,590.0d0,
# 600.0d0,610.0d0,620.0d0,630.0d0,640.0d0,
# 650.0d0,660.0d0,670.0d0,680.0d0,690.0d0,
# 700.0d0,710.0d0,720.0d0,730.0d0,740.0d0,
# 750.0d0,760.0d0,770.0d0,780.0d0,790.0d0,
# 800.0d0,810.0d0,820.0d0,830.0d0,840.0d0,
# 850.0d0,860.0d0,870.0d0,880.0d0,890.0d0,
# 900.0d0,910.0d0,920.0d0,930.0d0,940.0d0,
# 950.0d0,960.0d0,970.0d0,980.0d0,990.0d0, 1000.0d0/
*
DATA (yc(i),i=1,151)/4.183581334d0,4.545978356d0,
# 4.941666452d0,5.338524432d0,
# 5.710552332d0,5.871305639d0,5.980798334d0,5.989325219d0,
# 5.989771990d0,5.979384655d0,5.953749466d0,5.906497163d0,
# 5.827835443d0,5.702202128d0,5.507427195d0,5.221327470d0,
# 4.861849085d0,4.526867019d0,4.220897196d0,3.877727401d0,
# 3.525995801d0,3.200170282d0,2.911100907d0,2.657698720d0,
# 2.435102574d0,2.238185021d0,2.061849416d0,1.901483110d0,
# 1.753411931d0,1.614301638d0,1.480965572d0,1.348997661d0,
# 1.215803604d0,1.074412250d0,0.919019083d0,0.740332966d0,
# 0.530749433d0,0.298400521d0,0.077023976d0,-0.102260900d0,
# -0.270062463d0,-0.453870620d0,-0.643469731d0,-0.822931558d0,
# -0.987443985d0,-1.138691526d0,-1.275416891d0,-1.770985289d0,
# -1.910067996d0,-2.073780511d0,-2.387106650d0,-2.506786934d0,
# -2.507793900d0,-2.451529121d0,-2.379171774d0,-2.280179838d0,
# -2.173130895d0,-2.052273974d0,-1.953341364d0,-1.870193084d0,
# -1.795179493d0,-1.874006942d0,-1.991907667d0,-2.212751797d0,
# -2.281770963d0,-2.366744787d0,-2.479318144d0,-2.637428316d0,
# -2.936192030d0,-3.392136424d0,-3.518873910d0,-3.682965171d0,
# -3.718363034d0,-3.872185129d0,-3.990484562d0,-4.003435184d0,
# -4.047831469d0,-4.056871742d0,-4.096278566d0,-4.096147175d0,
# -4.069660494d0,-4.092229117d0,-4.050384852d0,-4.043310923d0,
# -4.021791606d0,-4.007601020d0,-3.922141199d0,-3.576132878d0,
# -3.231059452d0,-2.868110119d0,-2.588448565d0,-2.286823317d0,
# -1.981203740d0,-1.697922238d0,-1.465949078d0,-1.221405908d0,
# -0.997896019d0,-0.787814864d0,-0.603630999d0,-0.385428569d0,
# -0.191489219d0,0.012545468d0,0.180344393d0,0.350524042d0,
# 0.535367220d0,0.686993292d0,0.894786143d0,1.012194233d0,
# 1.192054257d0,1.364029746d0,1.515391057d0,1.655859873d0,
# 1.814367096d0,1.961412767d0,2.115811847d0,2.260602653d0,
# 2.429881738d0,2.579117247d0,2.705063788d0,2.852512294d0,
# 3.013412959d0,3.190250179d0,3.334985916d0,3.458011579d0,
# 3.596519369d0,3.772420987d0,3.927207591d0,4.083443432d0,
# 4.220565084d0,4.367619511d0,4.540226035d0,4.671279216d0,
# 4.822978517d0,5.004941301d0,5.133890543d0,5.317216628d0,
# 5.516475363d0,5.649275363d0,5.771911413d0,5.988071980d0,
# 6.100250110d0,6.265613527d0,6.437260016d0,6.622288364d0,
# 6.781564747d0,6.953709972d0,7.135417771d0,7.291456432d0,
# 7.475317109d0,7.623304078d0,7.829501710d0/
*
DATA (xu(i),i=1,152)/100.0d0,110.0d0,120.0d0,130.0d0,
# 140.0d0,145.0d0,150.0d0,151.0d0,
# 152.0d0,153.0d0,154.0d0,155.0d0,
# 156.0d0,157.0d0,158.0d0,159.0d0,
# 160.0d0,161.0d0,162.0d0,163.0d0,
# 164.0d0,165.0d0,166.0d0,167.0d0,
# 168.0d0,169.0d0,170.0d0,171.0d0,
# 172.0d0,173.0d0,174.0d0,175.0d0,
# 176.0d0,177.0d0,178.0d0,179.0d0,
# 180.0d0,181.0d0,182.0d0,183.0d0,
# 184.0d0,185.0d0,186.0d0,187.0d0,
# 188.0d0,189.0d0,190.0d0,195.0d0,
# 197.0d0,200.0d0,210.0d0,220.0d0,
# 230.0d0,240.0d0,250.0d0,260.0d0,
# 270.0d0,280.0d0,290.0d0,300.0d0,
# 310.0d0,320.0d0,330.0d0,335.0d0,
# 340.0d0,341.0d0,342.0d0,343.0d0,
# 344.0d0,345.0d0,345.3d0,345.5d0,
# 345.8d0,346.0d0,347.0d0,348.0d0,
# 349.0d0,350.0d0,351.0d0,352.0d0,
# 353.0d0,354.0d0,355.0d0,356.0d0,
# 357.0d0,358.0d0,359.0d0,360.0d0,
# 370.0d0,380.0d0,390.0d0,400.0d0,
# 410.0d0,420.0d0,430.0d0,440.0d0,
# 450.0d0,460.0d0,470.0d0,480.0d0,
# 490.0d0,500.0d0,510.0d0,520.0d0,
# 530.0d0,540.0d0,550.0d0,560.0d0,
# 570.0d0,580.0d0,590.0d0,600.0d0,
# 610.0d0,620.0d0,630.0d0,640.0d0,
# 650.0d0,660.0d0,670.0d0,680.0d0,
# 690.0d0,700.0d0,710.0d0,720.0d0,
# 730.0d0,740.0d0,750.0d0,760.0d0,
# 770.0d0,780.0d0,790.0d0,800.0d0,
# 810.0d0,820.0d0,830.0d0,840.0d0,
# 850.0d0,860.0d0,870.0d0,880.0d0,
# 890.0d0,900.0d0,910.0d0,920.0d0,
# 930.0d0,940.0d0,950.0d0,960.0d0,
# 970.0d0,980.0d0,990.0d0,1000.0d0/
*
DATA (yu(i),i=1,152)/4.187720723d0,4.549885500d0,
# 4.945324580d0,5.341899188d0,
# 5.713567160d0,5.874080782d0,5.983246138d0,5.991688858d0,
# 5.992041524d0,5.981547445d0,5.955788600d0,5.908389532d0,
# 5.829548098d0,5.703685183d0,5.508601342d0,5.222061322d0,
# 4.861934021d0,4.526068905d0,4.219144123d0,3.875133708d0,
# 3.522706643d0,3.196306140d0,2.906767625d0,2.652991136d0,
# 2.430066986d0,2.232861152d0,2.056219593d0,1.895596836d0,
# 1.747309568d0,1.607964744d0,1.474493488d0,1.342295851d0,
# 1.208919594d0,1.067309617d0,0.911812113d0,0.732879892d0,
# 0.522870581d0,0.290366470d0,0.068493581d0,-0.111577414d0,
# -0.279962552d0,-0.464332982d0,-0.654447961d0,-0.834336934d0,
# -0.999445210d0,-1.150871949d0,-1.287684579d0,-1.784688615d0,
# -1.924511090d0,-2.089081933d0,-2.402982024d0,-2.524611659d0,
# -2.525387652d0,-2.470417238d0,-2.399507117d0,-2.298953258d0,
# -2.192459462d0,-2.069700408d0,-1.973189437d0,-1.880904758d0,
# -1.807288585d0,-1.785086586d0,-1.826974607d0,-1.915990447d0,
# -2.062701208d0,-2.104793465d0,-2.152577512d0,-2.208540347d0,
# -2.276262386d0,-2.361739355d0,-2.391880716d0,-2.413272954d0,
# -2.447493257d0,-2.471820336d0,-2.622996766d0,-2.885017381d0,
# -3.696189261d0,-3.893536695d0,-4.011725974d0,-4.074946087d0,
# -4.137470131d0,-4.089349430d0,-4.136099918d0,-4.123934918d0,
# -4.136446777d0,-4.132255956d0,-4.118820717d0,-4.088009318d0,
# -3.788430916d0,-3.406565873d0,-3.050120268d0,-2.724817324d0,
# -2.410462214d0,-2.120474142d0,-1.845522666d0,-1.583137592d0,
# -1.337585290d0,-1.094144564d0,-0.897029424d0,-0.671551439d0,
# -0.478508120d0,-0.269184704d0,-0.072526771d0,0.115732816d0,
# 0.274371119d0,0.457987609d0,0.619147501d0,0.809852863d0,
# 0.968131131d0,1.121473977d0,1.290603896d0,1.462164252d0,
# 1.583836738d0,1.785815006d0,1.912473989d0,2.068564737d0,
# 2.235729083d0,2.398009567d0,2.538161792d0,2.681468958d0,
# 2.803535096d0,2.966400538d0,3.143956306d0,3.310124390d0,
# 3.429746138d0,3.580977910d0,3.727340256d0,3.890426351d0,
# 4.045983092d0,4.203320822d0,4.348209516d0,4.497115304d0,
# 4.654911988d0,4.823398533d0,4.983830013d0,5.127442851d0,
# 5.265056876d0,5.448819268d0,5.642707192d0,5.785441771d0,
# 5.936113103d0,6.108286154d0,6.258324908d0,6.415575321d0,
# 6.591607190d0,6.761124005d0,6.951304838d0,7.117437285d0,
# 7.283053264d0,7.453490877d0,7.637739897d0,7.799659928d0/
*
IF(top.eq.-1) THEN
n= 154
DO l=1,154
x(l)= xl(l)
y(l)= yl(l)
ENDDO
ELSEIF(top.eq.0) THEN
n= 151
DO l=1,151
x(l)= xc(l)
y(l)= yc(l)
ENDDO
ELSEIF(top.eq.1) THEN
n= 152
DO l=1,152
x(l)= xu(l)
y(l)= yu(l)
ENDDO
ENDIF
*.....Set up tridiagonal system.........................................
* b=diagonal, d=offdiagonal, c=right-hand side
*
d(1)= x(2)-x(1)
c(2)= (y(2)-y(1))/d(1)
DO k= 2,n-1
d(k)= x(k+1)-x(k)
b(k)= TWO*(d(k-1)+d(k))
c(k+1)= (y(k+1)-y(k))/d(k)
c(k)= c(k+1)-c(k)
END DO
*
*.....End conditions. third derivatives at x(1) and x(n) obtained
* from divided differences.......................................
*
b(1)= -d(1)
b(n)= -d(n-1)
c(1)= ZERO
c(n)= ZERO
IF (n.gt.3) THEN
c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
END IF
*
DO k=2,n ! forward elimination
t= d(k-1)/b(k-1)
b(k)= b(k)-t*d(k-1)
c(k)= c(k)-t*c(k-1)
END DO
*
c(n)= c(n)/b(n)
*
* back substitution ( makes c the sigma of text)
*
DO k=n-1,1,-1
c(k)= (c(k)-d(k)*c(k+1))/b(k)
END DO
*
*.....Compute polynomial coefficients...................................
*
b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
DO k=1,n-1
b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
d(k)= (c(k+1)-c(k))/d(k)
c(k)= THREE*c(k)
END DO
c(n)= THREE*c(n)
d(n)= d(n-1)
*
RETURN
*
END
*
*------------------------------------------------------------------------
*
SUBROUTINE Seval3Single(u,b,c,d,top,gdim,f,fp,fpp,fppp)
*
* ---------------------------------------------------------------------------
*
* PURPOSE - Evaluate the cubic spline function
* Seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
* where x(i) <= u < x(i+1)
*
* NOTES- if u<x(1), i=1 is used;if u>x(n), i=n is used
*
REAL*8 u
* abscissa at which the spline is to be evaluated
*
INTEGER j,k,n,l,top,gdim
*
REAL*8 xl(154),yl(154)
REAL*8 xc(151),yc(151)
REAL*8 xu(152),yu(152)
REAL*8 x(154),y(154)
REAL*8 b(gdim),c(gdim),d(gdim)
* linear,quadratic,cubic coeff
*
REAL*8 f,fp,fpp,fppp
* function, 1st,2nd,3rd deriv
*
INTEGER i
REAL*8 dx
REAL*8 ZERO, TWO, THREE
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
SAVE i
DATA i/1/
*
* The grid
*
DATA (xl(l),l= 1,154)/100.0d0,110.0d0,120.0d0,130.0d0,
# 140.0d0,145.0d0,150.0d0,151.0d0,
# 152.0d0,153.0d0,154.0d0,155.0d0,
# 156.0d0,157.0d0,158.0d0,159.0d0,
# 160.0d0,161.0d0,162.0d0,163.0d0,
# 164.0d0,165.0d0,166.0d0,
# 167.0d0,168.0d0,169.0d0,170.0d0,
# 171.0d0,172.0d0,173.0d0,174.0d0,
# 175.0d0,176.0d0,177.0d0,178.0d0,
# 179.0d0,180.0d0,181.0d0,182.0d0,
# 183.0d0,184.0d0,185.0d0,186.0d0,
# 187.0d0,188.0d0,189.0d0,190.0d0,
# 195.0d0,197.0d0,200.0d0,210.0d0,
# 220.0d0,230.0d0,240.0d0,250.0d0,
# 260.0d0,270.0d0,280.0d0,290.0d0,
# 300.0d0,310.0d0,320.0d0,330.0d0,
# 335.0d0,340.0d0,341.0d0,342.0d0,
# 343.0d0,344.0d0,345.0d0,345.3d0,
# 345.5d0,345.8d0,346.0d0,347.0d0,
# 348.0d0,349.0d0,349.5d0,349.75d0,
# 350.0d0,351.0d0,352.0d0,353.0d0,
# 354.0d0,355.0d0,356.0d0,357.0d0,
# 358.0d0,359.0d0,360.0d0,370.0d0,
# 380.0d0,390.0d0,400.0d0,410.0d0,
# 420.0d0,430.0d0,440.0d0,450.0d0,
# 460.0d0,470.0d0,480.0d0,490.0d0,
# 500.0d0,510.0d0,520.0d0,530.0d0,
# 540.0d0,550.0d0,560.0d0,570.0d0,
# 580.0d0,590.0d0,600.0d0,610.0d0,
# 620.0d0,630.0d0,640.0d0,650.0d0,
# 660.0d0,670.0d0,680.0d0,690.0d0,
# 700.0d0,710.0d0,720.0d0,730.0d0,
# 740.0d0,750.0d0,760.0d0,770.0d0,
# 780.0d0,790.0d0,800.0d0,810.0d0,
# 820.0d0,830.0d0,840.0d0,850.0d0,
# 860.0d0,870.0d0,880.0d0,890.0d0,
# 900.0d0,910.0d0,920.0d0,930.0d0,
# 940.0d0,950.0d0,960.0d0,970.0d0,
# 980.0d0,990.0d0,1000.0d0/
*
DATA (yl(l),l=1,154)/4.179467154d0,4.542088751d0,
# 4.938014960d0,5.335142457d0,
# 5.707514034d0,5.868498547d0,5.978310949d0,5.986920781d0,
# 5.987460634d0,5.977179301d0,5.951667414d0,5.904562092d0,
# 5.826081147d0,5.700679828d0,5.506218392d0,5.220567269d0,
# 4.861752146d0,4.527674442d0,4.222683291d0,3.880374785d0,
# 3.529354583d0,3.204116062d0,2.915522583d0,
# 2.662502577d0,2.440241875d0,2.243618198d0,2.067597988d0,
# 1.907471856d0,1.759641285d0,1.620735593d0,1.487565904d0,
# 1.355817302d0,1.222813267d0,1.081630216d0,0.926313642d0,
# 0.747981132d0,0.538696787d0,0.306541839d0,0.085775708d0,
# -0.092817756d0,-0.259991178d0,-0.443225249d0,-0.632350007d0,
# -0.811227830d0,-0.975263093d0,-1.126478218d0,-1.262721067d0,
# -1.757119775d0,-1.895339688d0,-2.057961996d0,-2.370870160d0,
# -2.488446977d0,-2.490894901d0,-2.433191749d0,-2.359142379d0,
# -2.260728908d0,-2.151155918d0,-2.036327859d0,-1.933794725d0,
# -1.862582453d0,-1.802070746d0,-1.812386184d0,-1.938240382d0,
# -2.108006605d0,-2.488306576d0,-2.655000325d0,-3.016952793d0,
# -3.680237527d0,-3.870995189d0,-3.959387796d0,-3.943637421d0,
# -3.982342165d0,-3.967290244d0,-4.002598626d0,-4.008430095d0,
# -4.010824003d0,-4.021419122d0,-4.018822218d0,-4.040141549d0,
# -4.024343762d0,-4.015666188d0,-3.965077563d0,-3.968563552d0,
# -3.963667776d0,-3.942755435d0,-3.932719674d0,-3.874191897d0,
# -3.876081328d0,-3.833375576d0,-3.796241844d0,-3.406242887d0,
# -3.088048772d0,-2.738617943d0,-2.446609831d0,-2.131731318d0,
# -1.878917579d0,-1.589864293d0,-1.358374992d0,-1.103676110d0,
# -0.884448436d0,-0.679633948d0,-0.472673942d0,-0.281736254d0,
# -0.108375079d0,0.093529154d0,0.257607461d0,0.430152846d0,
# 0.578396867d0,0.763320409d0,0.934522954d0,1.075378845d0,
# 1.249326065d0,1.426785290d0,1.535796936d0,1.715981764d0,
# 1.871449481d0,2.020126942d0,2.154595978d0,2.322995662d0,
# 2.456079279d0,2.608032219d0,2.749902433d0,2.905996823d0,
# 3.064612486d0,3.205141310d0,3.344347828d0,3.506879165d0,
# 3.660117972d0,3.783439131d0,3.947284161d0,4.097179155d0,
# 4.248755764d0,4.406763712d0,4.563193913d0,4.705927695d0,
# 4.853626937d0,5.015246327d0,5.167579841d0,5.341215028d0,
# 5.512112004d0,5.650306940d0,5.805087521d0,5.980800439d0,
# 6.125262600d0,6.272280377d0,6.477594572d0,6.625428154d0,
# 6.810107710d0,6.954753946d0,7.131768272d0,7.304668734d0,
# 7.458544663d0,7.653458644d0,7.821154640d0/
*
DATA (xc(l),l= 1,151)/100.0d0,110.0d0,120.0d0,130.0d0,140.0d0,
# 145.0d0,150.0d0,151.0d0,152.0d0,153.0d0,154.0d0,155.0d0,
# 156.0d0,157.0d0,158.0d0,159.0d0,160.0d0,161.0d0,162.0d0,
# 163.0d0,164.0d0,165.0d0,166.0d0,167.0d0,168.0d0,169.0d0,
# 170.0d0,171.0d0,172.0d0,173.0d0,174.0d0,175.0d0,176.0d0,
# 177.0d0,178.0d0,179.0d0,180.0d0,181.0d0,182.0d0,183.0d0,
# 184.0d0,185.0d0,186.0d0,187.0d0,188.0d0,
# 189.0d0,190.0d0,195.0d0,197.0d0,200.0d0,
# 210.0d0,220.0d0,230.0d0,240.0d0,250.0d0,
# 260.0d0,270.0d0,280.0d0,290.0d0,300.0d0,
# 320.0d0,330.0d0,335.0d0,340.0d0,341.0d0,
# 342.0d0,343.0d0,344.0d0,345.0d0,345.3d0,
# 345.5d0,345.8d0,346.0d0,347.0d0,348.0d0,
# 349.0d0,350.0d0,351.0d0,352.0d0,353.0d0,
# 354.0d0,355.0d0,356.0d0,357.0d0,358.0d0,
# 359.0d0,360.0d0,370.0d0,380.0d0,390.0d0,
# 400.0d0,410.0d0,420.0d0,430.0d0,440.0d0,
# 450.0d0,460.0d0,470.0d0,480.0d0,490.0d0,
# 500.0d0,510.0d0,520.0d0,530.0d0,540.0d0,
# 550.0d0,560.0d0,570.0d0,580.0d0,590.0d0,
# 600.0d0,610.0d0,620.0d0,630.0d0,640.0d0,
# 650.0d0,660.0d0,670.0d0,680.0d0,690.0d0,
# 700.0d0,710.0d0,720.0d0,730.0d0,740.0d0,
# 750.0d0,760.0d0,770.0d0,780.0d0,790.0d0,
# 800.0d0,810.0d0,820.0d0,830.0d0,840.0d0,
# 850.0d0,860.0d0,870.0d0,880.0d0,890.0d0,
# 900.0d0,910.0d0,920.0d0,930.0d0,940.0d0,
# 950.0d0,960.0d0,970.0d0,980.0d0,990.0d0, 1000.0d0/
*
DATA (yc(l),l= 1,151)/4.183581334d0,4.545978356d0,
# 4.941666452d0,5.338524432d0,
# 5.710552332d0,5.871305639d0,5.980798334d0,5.989325219d0,
# 5.989771990d0,5.979384655d0,5.953749466d0,5.906497163d0,
# 5.827835443d0,5.702202128d0,5.507427195d0,5.221327470d0,
# 4.861849085d0,4.526867019d0,4.220897196d0,3.877727401d0,
# 3.525995801d0,3.200170282d0,2.911100907d0,2.657698720d0,
# 2.435102574d0,2.238185021d0,2.061849416d0,1.901483110d0,
# 1.753411931d0,1.614301638d0,1.480965572d0,1.348997661d0,
# 1.215803604d0,1.074412250d0,0.919019083d0,0.740332966d0,
# 0.530749433d0,0.298400521d0,0.077023976d0,-0.102260900d0,
# -0.270062463d0,-0.453870620d0,-0.643469731d0,-0.822931558d0,
# -0.987443985d0,-1.138691526d0,-1.275416891d0,-1.770985289d0,
# -1.910067996d0,-2.073780511d0,-2.387106650d0,-2.506786934d0,
# -2.507793900d0,-2.451529121d0,-2.379171774d0,-2.280179838d0,
# -2.173130895d0,-2.052273974d0,-1.953341364d0,-1.870193084d0,
# -1.795179493d0,-1.874006942d0,-1.991907667d0,-2.212751797d0,
# -2.281770963d0,-2.366744787d0,-2.479318144d0,-2.637428316d0,
# -2.936192030d0,-3.392136424d0,-3.518873910d0,-3.682965171d0,
# -3.718363034d0,-3.872185129d0,-3.990484562d0,-4.003435184d0,
# -4.047831469d0,-4.056871742d0,-4.096278566d0,-4.096147175d0,
# -4.069660494d0,-4.092229117d0,-4.050384852d0,-4.043310923d0,
# -4.021791606d0,-4.007601020d0,-3.922141199d0,-3.576132878d0,
# -3.231059452d0,-2.868110119d0,-2.588448565d0,-2.286823317d0,
# -1.981203740d0,-1.697922238d0,-1.465949078d0,-1.221405908d0,
# -0.997896019d0,-0.787814864d0,-0.603630999d0,-0.385428569d0,
# -0.191489219d0,0.012545468d0,0.180344393d0,0.350524042d0,
# 0.535367220d0,0.686993292d0,0.894786143d0,1.012194233d0,
# 1.192054257d0,1.364029746d0,1.515391057d0,1.655859873d0,
# 1.814367096d0,1.961412767d0,2.115811847d0,2.260602653d0,
# 2.429881738d0,2.579117247d0,2.705063788d0,2.852512294d0,
# 3.013412959d0,3.190250179d0,3.334985916d0,3.458011579d0,
# 3.596519369d0,3.772420987d0,3.927207591d0,4.083443432d0,
# 4.220565084d0,4.367619511d0,4.540226035d0,4.671279216d0,
# 4.822978517d0,5.004941301d0,5.133890543d0,5.317216628d0,
# 5.516475363d0,5.649275363d0,5.771911413d0,5.988071980d0,
# 6.100250110d0,6.265613527d0,6.437260016d0,6.622288364d0,
# 6.781564747d0,6.953709972d0,7.135417771d0,7.291456432d0,
# 7.475317109d0,7.623304078d0,7.829501710d0/
*
DATA (xu(l),l=1,152)/100.0d0,110.0d0,120.0d0,130.0d0,
# 140.0d0,145.0d0,150.0d0,151.0d0,
# 152.0d0,153.0d0,154.0d0,155.0d0,
# 156.0d0,157.0d0,158.0d0,159.0d0,
# 160.0d0,161.0d0,162.0d0,163.0d0,
# 164.0d0,165.0d0,166.0d0,167.0d0,
# 168.0d0,169.0d0,170.0d0,171.0d0,
# 172.0d0,173.0d0,174.0d0,175.0d0,
# 176.0d0,177.0d0,178.0d0,179.0d0,
# 180.0d0,181.0d0,182.0d0,183.0d0,
# 184.0d0,185.0d0,186.0d0,187.0d0,
# 188.0d0,189.0d0,190.0d0,195.0d0,
# 197.0d0,200.0d0,210.0d0,220.0d0,
# 230.0d0,240.0d0,250.0d0,260.0d0,
# 270.0d0,280.0d0,290.0d0,300.0d0,
# 310.0d0,320.0d0,330.0d0,335.0d0,
# 340.0d0,341.0d0,342.0d0,343.0d0,
# 344.0d0,345.0d0,345.3d0,345.5d0,
# 345.8d0,346.0d0,347.0d0,348.0d0,
# 349.0d0,350.0d0,351.0d0,352.0d0,
# 353.0d0,354.0d0,355.0d0,356.0d0,
# 357.0d0,358.0d0,359.0d0,360.0d0,
# 370.0d0,380.0d0,390.0d0,400.0d0,
# 410.0d0,420.0d0,430.0d0,440.0d0,
# 450.0d0,460.0d0,470.0d0,480.0d0,
# 490.0d0,500.0d0,510.0d0,520.0d0,
# 530.0d0,540.0d0,550.0d0,560.0d0,
# 570.0d0,580.0d0,590.0d0,600.0d0,
# 610.0d0,620.0d0,630.0d0,640.0d0,
# 650.0d0,660.0d0,670.0d0,680.0d0,
# 690.0d0,700.0d0,710.0d0,720.0d0,
# 730.0d0,740.0d0,750.0d0,760.0d0,
# 770.0d0,780.0d0,790.0d0,800.0d0,
# 810.0d0,820.0d0,830.0d0,840.0d0,
# 850.0d0,860.0d0,870.0d0,880.0d0,
# 890.0d0,900.0d0,910.0d0,920.0d0,
# 930.0d0,940.0d0,950.0d0,960.0d0,
# 970.0d0,980.0d0,990.0d0,1000.0d0/
*
DATA (yu(l),l=1,152)/4.187720723d0,4.549885500d0,
# 4.945324580d0,5.341899188d0,
# 5.713567160d0,5.874080782d0,5.983246138d0,5.991688858d0,
# 5.992041524d0,5.981547445d0,5.955788600d0,5.908389532d0,
# 5.829548098d0,5.703685183d0,5.508601342d0,5.222061322d0,
# 4.861934021d0,4.526068905d0,4.219144123d0,3.875133708d0,
# 3.522706643d0,3.196306140d0,2.906767625d0,2.652991136d0,
# 2.430066986d0,2.232861152d0,2.056219593d0,1.895596836d0,
# 1.747309568d0,1.607964744d0,1.474493488d0,1.342295851d0,
# 1.208919594d0,1.067309617d0,0.911812113d0,0.732879892d0,
# 0.522870581d0,0.290366470d0,0.068493581d0,-0.111577414d0,
# -0.279962552d0,-0.464332982d0,-0.654447961d0,-0.834336934d0,
# -0.999445210d0,-1.150871949d0,-1.287684579d0,-1.784688615d0,
# -1.924511090d0,-2.089081933d0,-2.402982024d0,-2.524611659d0,
# -2.525387652d0,-2.470417238d0,-2.399507117d0,-2.298953258d0,
# -2.192459462d0,-2.069700408d0,-1.973189437d0,-1.880904758d0,
# -1.807288585d0,-1.785086586d0,-1.826974607d0,-1.915990447d0,
# -2.062701208d0,-2.104793465d0,-2.152577512d0,-2.208540347d0,
# -2.276262386d0,-2.361739355d0,-2.391880716d0,-2.413272954d0,
# -2.447493257d0,-2.471820336d0,-2.622996766d0,-2.885017381d0,
# -3.696189261d0,-3.893536695d0,-4.011725974d0,-4.074946087d0,
# -4.137470131d0,-4.089349430d0,-4.136099918d0,-4.123934918d0,
# -4.136446777d0,-4.132255956d0,-4.118820717d0,-4.088009318d0,
# -3.788430916d0,-3.406565873d0,-3.050120268d0,-2.724817324d0,
# -2.410462214d0,-2.120474142d0,-1.845522666d0,-1.583137592d0,
# -1.337585290d0,-1.094144564d0,-0.897029424d0,-0.671551439d0,
# -0.478508120d0,-0.269184704d0,-0.072526771d0,0.115732816d0,
# 0.274371119d0,0.457987609d0,0.619147501d0,0.809852863d0,
# 0.968131131d0,1.121473977d0,1.290603896d0,1.462164252d0,
# 1.583836738d0,1.785815006d0,1.912473989d0,2.068564737d0,
# 2.235729083d0,2.398009567d0,2.538161792d0,2.681468958d0,
# 2.803535096d0,2.966400538d0,3.143956306d0,3.310124390d0,
# 3.429746138d0,3.580977910d0,3.727340256d0,3.890426351d0,
# 4.045983092d0,4.203320822d0,4.348209516d0,4.497115304d0,
# 4.654911988d0,4.823398533d0,4.983830013d0,5.127442851d0,
# 5.265056876d0,5.448819268d0,5.642707192d0,5.785441771d0,
# 5.936113103d0,6.108286154d0,6.258324908d0,6.415575321d0,
# 6.591607190d0,6.761124005d0,6.951304838d0,7.117437285d0,
# 7.283053264d0,7.453490877d0,7.637739897d0,7.799659928d0/
*
IF(top.eq.-1) THEN
n= 154
DO l=1,154
x(l)= xl(l)
y(l)= yl(l)
ENDDO
ELSEIF(top.eq.0) THEN
n= 151
DO l=1,151
x(l)= xc(l)
y(l)= yc(l)
ENDDO
ELSEIF(top.eq.1) THEN
n= 152
DO l=1,152
x(l)= xu(l)
y(l)= yu(l)
ENDDO
ENDIF
*
*.....First check if u is in the same interval found on the
* last call to Seval.............................................
*
IF ( (i.lt.1) .OR. (i.ge.n) ) i=1
IF ( (u.lt.x(i)) .OR. (u.ge.x(i+1)) ) THEN
i=1
*
* binary search
*
j= n+1
2468 k= (i+j)/2
IF (u.lt.x(k)) THEN
j= k
ELSE
i= k
ENDIF
IF (j.le.i+1) GOTO 2469
GOTO 2468
2469 CONTINUE
ENDIF
*
dx= u-x(i)
*
* evaluate the spline
*
f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i))
fpp= TWO*c(i) + dx*SIX*d(i)
fppp= SIX*d(i)
*
RETURN
*
END
*****************************************************************
C---electroweak corrections to H --> gamma gamma
*****************************************************************
double precision function gaga_elw(mtop,mh)
*
IMPLICIT NONE
*
INTEGER i,top,gdim
REAL*8 u,ui,value
REAL*8 vp,vpp,vppp
REAL*8 mtop,mh
REAL*8 x,x0,x1,x2,y0,y1,y2,a0,a1,a2
REAL*8 value0,value1,value2
REAL*8 bl(226),cl(226),dl(226)
REAL*8 bc(223),cc(223),dc(223)
REAL*8 bu(225),cu(225),du(225)
*
* u value of M_H at which the spline is to be evaluated
* top= -1,0,1 lower, central, upper value for m_top
*
u = mh
top = -1
gdim= 226
CALL FMMsplineSingle0(bl,cl,dl,top,gdim)
CALL Seval3Single0(u,bl,cl,dl,top,gdim,value0,vp,vpp,vppp)
top = 0
gdim= 223
CALL FMMsplineSingle0(bc,cc,dc,top,gdim)
CALL Seval3Single0(u,bc,cc,dc,top,gdim,value1,vp,vpp,vppp)
top = 1
gdim= 225
CALL FMMsplineSingle0(bu,cu,du,top,gdim)
CALL Seval3Single0(u,bu,cu,du,top,gdim,value2,vp,vpp,vppp)
x = mtop
x0 = 171.06d0
x1 = 172.64d0
x2 = 174.22d0
y0 = value0
y1 = value1
y2 = value2
A0=(X-X1)*(X-X2)/(X0-X1)/(X0-X2)
A1=(X-X0)*(X-X2)/(X1-X0)/(X1-X2)
A2=(X-X0)*(X-X1)/(X2-X0)/(X2-X1)
gaga_elw=(A0*Y0+A1*Y1+A2*Y2)/100.d0
c if(mh.gt.170.d0) gaga_elw = 0
RETURN
END
*
*-----------------------------------------------------------------------
*
c CONTAINS
*
SUBROUTINE FMMsplineSingle0(b,c,d,top,gdim)
*
* ---------------------------------------------------------------------------
*
* PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline
* so that the interpolated value is given by
* s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3
* when x(k) <= x <= x(k+1)
* The end conditions match the third derivatives of the interpolated curve to
* the third derivatives of the unique polynomials thru the first four and
* last four points.
* Use Seval or Seval3 to evaluate the spline.
*
INTEGER k,n,i,top,gdim,l
*
REAL*8 xl(226),yl(226)
REAL*8 xc(223),yc(223)
REAL*8 xu(225),yu(225)
REAL*8 x(226),y(226)
*
REAL*8 b(gdim)
* linear coeff
*
REAL*8 c(gdim)
* quadratic coeff.
*
REAL*8 d(gdim)
* cubic coeff.
*
REAL*8 t
REAL*8 ZERO, TWO, THREE
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
*
* The grid
*
DATA (xl(i),i=1,226)/100.0000d0,105.0000d0,110.0000d0,
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,179.0000d0,
# 180.0000d0,181.0000d0,182.0000d0,183.0000d0,184.0000d0,
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
# 340.5000d0,340.8000d0,341.0000d0,341.2000d0,341.4000d0,
# 341.6000d0,341.7000d0,341.8000d0,341.8500d0,341.9000d0,
# 341.9500d0,342.0000d0,342.0200d0,342.0400d0,342.0600d0,
# 342.0800d0,342.1000d0,342.1100d0,342.1120d0,342.1140d0,
# 342.1160d0,342.1180d0,342.1190d0,342.1196d0,342.1198d0,
# 342.1199d0,342.1204d0,342.1220d0,342.2000d0,342.3000d0,
# 342.4000d0,342.5000d0,342.6000d0,342.7000d0,342.8000d0,
# 342.9000d0,343.0000d0,343.2000d0,343.4000d0,343.6000d0,
# 343.8000d0,344.0000d0,345.0000d0,346.0000d0,347.0000d0,
# 348.0000d0,349.0000d0,350.0000d0,351.0000d0,352.0000d0,
# 353.0000d0,354.0000d0,355.0000d0,356.0000d0,357.0000d0,
# 358.0000d0,359.0000d0,360.0000d0,370.0000d0,380.0000d0,
# 390.0000d0,400.0000d0,410.0000d0,420.0000d0,430.0000d0,
# 440.0000d0,450.0000d0,460.0000d0,470.0000d0,480.0000d0,
# 490.0000d0,500.0000d0,510.0000d0,520.0000d0,530.0000d0,
# 540.0000d0,550.0000d0,560.0000d0,565.0000d0,570.0000d0,
# 575.0000d0,580.0000d0,585.0000d0,590.0000d0,595.0000d0,
# 600.0000d0,605.0000d0,610.0000d0,615.0000d0,620.0000d0,
# 625.0000d0,630.0000d0,635.0000d0,640.0000d0,641.0000d0,
# 642.0000d0,643.0000d0,644.0000d0,645.0000d0,646.0000d0,
# 647.0000d0,648.0000d0,649.0000d0,650.0000d0,651.0000d0,
# 652.0000d0,653.0000d0,654.0000d0,655.0000d0,656.0000d0,
# 657.0000d0,658.0000d0,659.0000d0,660.0000d0,670.0000d0,
# 680.0000d0,690.0000d0,700.0000d0,710.0000d0,720.0000d0,
# 730.0000d0,740.0000d0,750.0000d0,760.0000d0,770.0000d0,
# 780.0000d0,790.0000d0,800.0000d0,810.0000d0,820.0000d0,
# 830.0000d0,840.0000d0,850.0000d0,860.0000d0,870.0000d0,
# 880.0000d0,890.0000d0,900.0000d0,910.0000d0,920.0000d0,
# 930.0000d0,940.0000d0,950.0000d0,960.0000d0,970.0000d0,
# 980.0000d0,990.0000d0,1000.0000d0/
*
DATA (yl(i),i=1,226)/-2.927515523d0,-2.713875277d0,-2.477263056d0,
# -2.215702562d0, -1.926872830d0, -1.607794154d0, -1.254642316d0,
# -0.862759386d0, -0.424192821d0, 0.071393692d0, 0.642180822d0,
# 1.311566449d0, 1.458302972d0, 1.610089467d0, 1.771001228d0,
# 1.958842368d0, 1.980856139d0, 2.003769992d0, 2.027683661d0,
# 2.052701483d0, 2.078930176d0, 2.106476183d0, 2.135441928d0,
# 2.165921345d0, 2.197994437d0, 2.231721476d0, 2.323417121d0,
# 2.425109533d0, 2.534671123d0, 2.648366651d0, 2.869451789d0,
# 3.056555683d0, 3.197055532d0, 3.294656543d0, 3.399442397d0,
# 3.438265245d0, 3.450036365d0, 3.449871888d0, 3.442815310d0,
# 3.428262691d0, 3.411797242d0, 3.363961271d0, 3.295179282d0,
# 3.287861332d0, 3.325219720d0, 3.470009593d0, 3.737173473d0,
# 4.021846379d0, 4.241836705d0, 4.394931599d0, 4.501421364d0,
# 4.577551451d0, 4.633363121d0, 4.675411596d0, 4.707771996d0,
# 4.732593619d0, 4.751501682d0, 4.766245197d0, 4.778023555d0,
# 4.787199047d0, 4.793849566d0, 4.798265576d0, 4.801110206d0,
# 4.802711433d0, 4.767182955d0, 4.677220707d0, 4.558671105d0,
# 4.423406368d0, 4.276574444d0, 4.119375511d0, 3.974083777d0,
# 3.819733346d0, 3.671217773d0, 3.517370222d0, 3.344189966d0,
# 3.116802834d0, 2.754066836d0, 1.809839704d0, 1.691727991d0,
# 1.607871864d0, 1.544493194d0, 1.473234377d0, 1.390992309d0,
# 1.292466985d0, 1.234537350d0, 1.167169430d0, 1.128863846d0,
# 1.085963030d0, 1.037416400d0, 0.979461824d0, 0.952439442d0,
# 0.922628624d0, 0.888386832d0, 0.846809865d0, 0.791742524d0,
# 0.752347274d0, 0.742175952d0, 0.730422928d0, 0.716838887d0,
# 0.697956938d0, 0.684261787d0, 0.672163545d0, 0.665924128d0,
# 0.661453423d0, 0.652439216d0, 0.656998002d0, 0.677275168d0,
# 0.699459467d0, 0.721833458d0, 0.738243489d0, 0.751895487d0,
# 0.769698326d0, 0.779779847d0, 0.799512664d0, 0.816893187d0,
# 0.842794338d0, 0.874519045d0, 0.916496501d0, 0.950264630d0,
# 0.968728477d0, 1.105576625d0, 1.278360270d0, 1.391736231d0,
# 1.542276157d0, 1.675869257d0, 1.795595314d0, 1.923303286d0,
# 2.025352084d0, 2.157383377d0, 2.260417650d0, 2.361640519d0,
# 2.479888170d0, 2.572642417d0, 2.688809606d0, 2.796999182d0,
# 2.879392617d0, 3.848596335d0, 4.732837134d0, 5.525766484d0,
# 6.243989288d0, 6.975830059d0, 7.607337162d0, 8.167782071d0,
# 8.790374186d0, 9.422028844d0, 9.922844338d0, 10.426856521d0,
# 10.933372206d0, 11.329434560d0, 11.735173309d0, 11.795158880d0,
# 11.720854043d0, 11.455441492d0, 10.516061928d0, 9.017310079d0,
# 7.000279246d0, 5.560369962d0, 3.715234200d0, 1.373399680d0,
# -1.698051464d0, -5.394210410d0, -9.477907268d0,-13.677692505d0,
# -18.256012300d0,-24.029072069d0,-29.992304980d0,-35.639341296d0,
# -41.517060283d0,-46.430415213d0,-50.471411203d0,-53.247623713d0,
# -55.184140125d0,-55.523807703d0,-55.749322940d0,-55.785869793d0,
# -56.086101947d0,-56.289055861d0,-56.369321707d0,-56.555818993d0,
# -56.632365789d0,-56.616832544d0,-56.709629986d0,-56.539655434d0,
# -56.366567734d0,-56.272247546d0,-56.170833793d0,-55.984866548d0,
# -55.776968735d0,-55.354292076d0,-55.136594171d0,-54.856130174d0,
# -54.390381982d0,-50.585126291d0,-46.161103219d0,-41.187366158d0,
# -36.554683053d0,-32.185245783d0,-28.097939638d0,-24.608441494d0,
# -21.643195994d0,-18.984106588d0,-16.998672281d0,-15.273879570d0,
# -13.373334219d0,-11.294832079d0, -9.584938331d0, -8.286389135d0,
# -6.838755740d0, -5.409090746d0, -4.588339018d0, -3.705250749d0,
# -2.476432078d0, -1.367536281d0, -0.414988031d0, 0.510347996d0,
# 1.376925204d0, 2.276754999d0, 3.203065379d0, 3.629617497d0,
# 3.821185286d0, 4.122178708d0, 4.804270233d0, 5.929961061d0,
# 7.072311664d0, 7.458878393d0, 7.823110180d0/
*
DATA (xc(i),i=1,223)/100.0000d0,105.0000d0,110.0000d0,
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,178.5000d0,
# 179.0000d0,179.5000d0,180.0000d0,180.5000d0,181.0000d0,
# 181.5000d0,182.0000d0,182.5000d0,183.0000d0,184.0000d0,
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
# 341.0000d0,342.0000d0,343.0000d0,344.0000d0,344.2000d0,
# 344.4000d0,344.6000d0,344.8000d0,344.9000d0,345.0000d0,
# 345.0500d0,345.1000d0,345.1500d0,345.2000d0,345.2200d0,
# 345.2400d0,345.2600d0,345.2700d0,345.2720d0,345.2740d0,
# 345.2760d0,345.2780d0,345.2790d0,345.2796d0,345.2798d0,
# 345.2799d0,345.2804d0,345.2820d0,345.2900d0,345.3500d0,
# 345.4500d0,345.5200d0,345.6000d0,345.7000d0,345.8000d0,
# 346.0000d0,347.0000d0,348.0000d0,349.0000d0,350.0000d0,
# 351.0000d0,352.0000d0,353.0000d0,354.0000d0,355.0000d0,
# 356.0000d0,357.0000d0,358.0000d0,359.0000d0,360.0000d0,
# 370.0000d0,380.0000d0,390.0000d0,400.0000d0,410.0000d0,
# 420.0000d0,430.0000d0,440.0000d0,450.0000d0,460.0000d0,
# 470.0000d0,480.0000d0,490.0000d0,500.0000d0,510.0000d0,
# 520.0000d0,530.0000d0,540.0000d0,550.0000d0,560.0000d0,
# 565.0000d0,570.0000d0,575.0000d0,580.0000d0,585.0000d0,
# 590.0000d0,595.0000d0,600.0000d0,605.0000d0,610.0000d0,
# 615.0000d0,620.0000d0,625.0000d0,630.0000d0,635.0000d0,
# 640.0000d0,641.0000d0,642.0000d0,643.0000d0,644.0000d0,
# 645.0000d0,646.0000d0,647.0000d0,648.0000d0,649.0000d0,
# 650.0000d0,651.0000d0,652.0000d0,653.0000d0,654.0000d0,
# 655.0000d0,656.0000d0,657.0000d0,658.0000d0,659.0000d0,
# 660.0000d0,670.0000d0,680.0000d0,690.0000d0,700.0000d0,
# 710.0000d0,720.0000d0,730.0000d0,740.0000d0,750.0000d0,
# 760.0000d0,770.0000d0,780.0000d0,790.0000d0,800.0000d0,
# 810.0000d0,820.0000d0,830.0000d0,840.0000d0,850.0000d0,
# 860.0000d0,870.0000d0,880.0000d0,890.0000d0,900.0000d0,
# 910.0000d0,920.0000d0,930.0000d0,940.0000d0,950.0000d0,
# 960.0000d0,970.0000d0,980.0000d0,990.0000d0,1000.0000d0/
*
DATA (yc(i),i=1,223)/-2.969458664d0,-2.755568498d0,-2.518689112d0,
# -2.256840599d0, -1.967703151d0, -1.648299088d0, -1.294794285d0,
# -0.902525232d0, -0.463541061d0, 0.032517848d0, 0.603841958d0,
# 1.273887882d0, 1.420789352d0, 1.572747980d0, 1.733843353d0,
# 1.921879203d0, 1.943912705d0, 1.966846336d0, 1.990779880d0,
# 2.015817728d0, 2.042066621d0, 2.069633010d0, 2.098619322d0,
# 2.129119489d0, 2.161213532d0, 2.194961744d0, 2.286711412d0,
# 2.388459445d0, 2.498078035d0, 2.611830784d0, 2.833024715d0,
# 3.020227220d0, 3.160820150d0, 3.258502948d0, 3.363381747d0,
# 3.402305539d0, 3.414208425d0, 3.414049050d0, 3.407053730d0,
# 3.392629803d0, 3.376195597d0, 3.328499111d0, 3.268651902d0,
# 3.259860594d0, 3.253486002d0, 3.252546948d0, 3.262154622d0,
# 3.289924633d0, 3.345085115d0, 3.434732330d0, 3.557729513d0,
# 3.701897111d0, 3.986569984d0, 4.206557658d0, 4.359682402d0,
# 4.466176116d0, 4.542300703d0, 4.598184529d0, 4.640273713d0,
# 4.672604523d0, 4.697438255d0, 4.716361927d0, 4.731103154d0,
# 4.742925539d0, 4.752092918d0, 4.758670785d0, 4.763130070d0,
# 4.766074843d0, 4.767748196d0, 4.732269356d0, 4.642413616d0,
# 4.524181528d0, 4.389510854d0, 4.242822112d0, 4.086217524d0,
# 3.941680893d0, 3.789691168d0, 3.644241595d0, 3.494351920d0,
# 3.331867639d0, 3.122929519d0, 2.816374576d0, 2.186038507d0,
# 2.068538445d0, 1.926511260d0, 1.746359564d0, 1.494945494d0,
# 1.428923128d0, 1.354291504d0, 1.267892764d0, 1.163289817d0,
# 1.100587592d0, 1.027248191d0, 0.984553881d0, 0.936444199d0,
# 0.879742362d0, 0.808975906d0, 0.774145495d0, 0.731724792d0,
# 0.675557165d0, 0.635406715d0, 0.625004410d0, 0.613030568d0,
# 0.599164524d0, 0.579800865d0, 0.565857087d0, 0.553503488d0,
# 0.547057831d0, 0.542411799d0, 0.533743538d0, 0.538671942d0,
# 0.542814415d0, 0.562187598d0, 0.575771207d0, 0.586370628d0,
# 0.605546923d0, 0.622661850d0, 0.639021282d0, 0.683852301d0,
# 0.846885709d0, 0.995547199d0, 1.139428677d0, 1.270752033d0,
# 1.403862461d0, 1.540060220d0, 1.676615197d0, 1.791154110d0,
# 1.909977528d0, 2.023870438d0, 2.138393119d0, 2.233199625d0,
# 2.358706648d0, 2.446271802d0, 3.459751114d0, 4.411864477d0,
# 5.220622691d0, 5.988967245d0, 6.706065217d0, 7.359912016d0,
# 7.993290467d0, 8.675262256d0, 9.327550412d0, 9.917129519d0,
# 10.373037136d0, 11.009210530d0, 11.516244354d0, 11.985967329d0,
# 12.198319618d0, 12.357011718d0, 12.311558072d0, 11.754465713d0,
# 10.617092695d0, 9.043395280d0, 8.002055503d0, 6.601203443d0,
# 4.587232698d0, 1.906399387d0, -1.281177539d0, -5.052281649d0,
# -9.213421228d0,-13.868726558d0,-19.489194521d0,-25.988103189d0,
# -32.734525726d0,-39.402198052d0,-45.903352541d0,-51.189352305d0,
# -55.322071720d0,-58.186905136d0,-58.653812374d0,-59.048160142d0,
# -59.279596622d0,-59.992120796d0,-60.271154555d0,-60.549383981d0,
# -60.729563975d0,-60.954540486d0,-61.084961955d0,-61.164811836d0,
# -61.186911575d0,-61.109847454d0,-61.108201824d0,-60.992397754d0,
# -60.768938455d0,-60.645603901d0,-60.224362756d0,-59.860911894d0,
# -59.548720101d0,-59.301307840d0,-55.034929471d0,-49.745997496d0,
# -44.181440809d0,-39.061582838d0,-34.158470096d0,-29.603188102d0,
# -25.799341483d0,-22.634533844d0,-19.804942358d0,-17.654326276d0,
# -15.913630982d0,-13.830699918d0,-11.633319645d0,-10.006868842d0,
# -8.532937676d0, -7.132846889d0, -5.627566626d0, -4.758444504d0,
# -3.887401469d0, -2.729580751d0, -1.512994295d0, -0.550651909d0,
# 0.355022046d0, 1.248689020d0, 2.123963438d0, 3.014801413d0,
# 3.484406765d0, 3.696913629d0, 4.008426964d0, 4.651255618d0,
# 5.840647264d0, 6.935940541d0, 7.364646872d0, 7.670905366d0/
*
DATA (xu(i),i=1,225)/100.0000d0,105.0000d0,110.0000d0,
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,178.5000d0,
# 179.0000d0,179.5000d0,180.0000d0,180.5000d0,181.0000d0,
# 181.5000d0,182.0000d0,182.5000d0,183.0000d0,184.0000d0,
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
# 341.0000d0,342.0000d0,343.0000d0,344.0000d0,345.0000d0,
# 345.5000d0,345.8000d0,346.0000d0,347.0000d0,347.2000d0,
# 347.4000d0,347.6000d0,347.8000d0,348.0000d0,348.1000d0,
# 348.2000d0,348.2500d0,348.3000d0,348.3200d0,348.3400d0,
# 348.3600d0,348.3800d0,348.4000d0,348.4100d0,348.4200d0,
# 348.4300d0,348.4320d0,348.4340d0,348.4360d0,348.4380d0,
# 348.4390d0,348.4396d0,348.4398d0,348.4399d0,348.4404d0,
# 348.4420d0,348.4500d0,348.5000d0,348.7000d0,349.0000d0,
# 349.4000d0,350.0000d0,351.0000d0,352.0000d0,353.0000d0,
# 354.0000d0,355.0000d0,356.0000d0,357.0000d0,358.0000d0,
# 359.0000d0,360.0000d0,370.0000d0,380.0000d0,390.0000d0,
# 400.0000d0,410.0000d0,420.0000d0,430.0000d0,440.0000d0,
# 450.0000d0,460.0000d0,470.0000d0,480.0000d0,490.0000d0,
# 500.0000d0,510.0000d0,520.0000d0,530.0000d0,540.0000d0,
# 550.0000d0,560.0000d0,565.0000d0,570.0000d0,575.0000d0,
# 580.0000d0,585.0000d0,590.0000d0,595.0000d0,600.0000d0,
# 605.0000d0,610.0000d0,615.0000d0,620.0000d0,625.0000d0,
# 630.0000d0,635.0000d0,640.0000d0,641.0000d0,642.0000d0,
# 643.0000d0,644.0000d0,645.0000d0,646.0000d0,647.0000d0,
# 648.0000d0,649.0000d0,650.0000d0,651.0000d0,652.0000d0,
# 653.0000d0,654.0000d0,655.0000d0,656.0000d0,657.0000d0,
# 658.0000d0,659.0000d0,660.0000d0,670.0000d0,680.0000d0,
# 690.0000d0,700.0000d0,710.0000d0,720.0000d0,730.0000d0,
# 740.0000d0,750.0000d0,760.0000d0,770.0000d0,780.0000d0,
# 790.0000d0,800.0000d0,810.0000d0,820.0000d0,830.0000d0,
# 840.0000d0,850.0000d0,860.0000d0,870.0000d0,880.0000d0,
# 890.0000d0,900.0000d0,910.0000d0,920.0000d0,930.0000d0,
# 940.0000d0,950.0000d0,960.0000d0,970.0000d0,980.0000d0,
# 990.0000d0,1000.0000d0/
*
DATA (yu(i),i=1,225)/-3.011799533d0,-2.797658814d0,-2.560511342d0,
# -2.298374579d0, -2.008927458d0, -1.689194746d0, -1.335334782d0,
# -0.942678003d0, -0.503267763d0, -0.006741983d0, 0.565124142d0,
# 1.235840116d0, 1.382905273d0, 1.535044791d0, 1.696328604d0,
# 1.884567714d0, 1.906622347d0, 1.929577239d0, 1.953532153d0,
# 1.978591436d0, 2.004861792d0, 2.032449634d0, 2.061457355d0,
# 2.091978868d0, 2.124094156d0, 2.157863499d0, 2.249665540d0,
# 2.351465763d0, 2.461136573d0, 2.574941494d0, 2.796240093d0,
# 2.983545301d0, 3.124228166d0, 3.221982342d0, 3.326972382d0,
# 3.365993805d0, 3.378023193d0, 3.377887062d0, 3.370931083d0,
# 3.356634680d0, 3.340265619d0, 3.292682936d0, 3.232961811d0,
# 3.224194020d0, 3.217828260d0, 3.216889093d0, 3.226496714d0,
# 3.254271720d0, 3.309440914d0, 3.399097855d0, 3.522103315d0,
# 3.666276367d0, 3.950949463d0, 4.170926862d0, 4.324070476d0,
# 4.430588837d0, 4.506706461d0, 4.562622634d0, 4.604788444d0,
# 4.637112788d0, 4.661928180d0, 4.680880990d0, 4.695625131d0,
# 4.707459275d0, 4.716670603d0, 4.723178694d0, 4.727608289d0,
# 4.730650318d0, 4.732442748d0, 4.697030978d0, 4.607168713d0,
# 4.489397175d0, 4.355888392d0, 4.208987612d0, 4.052707904d0,
# 3.908019896d0, 3.759968739d0, 3.615348691d0, 3.469932507d0,
# 3.315332978d0, 3.122286809d0, 2.855004727d0, 2.377722123d0,
# 2.300655190d0, 2.213908305d0, 2.115631830d0, 2.000162289d0,
# 1.860840425d0, 1.779119559d0, 1.724912832d0, 1.685992661d0,
# 1.446048759d0, 1.384730617d0, 1.316264190d0, 1.238304389d0,
# 1.147571101d0, 1.036738128d0, 0.969143795d0, 0.888640045d0,
# 0.841031364d0, 0.785526212d0, 0.759880178d0, 0.731970096d0,
# 0.701136254d0, 0.665654946d0, 0.622492994d0, 0.596679361d0,
# 0.565315972d0, 0.524528227d0, 0.513928834d0, 0.501826683d0,
# 0.487735928d0, 0.468019598d0, 0.454013349d0, 0.441470198d0,
# 0.434889761d0, 0.430148561d0, 0.421464835d0, 0.429200773d0,
# 0.429255984d0, 0.437861500d0, 0.480919141d0, 0.544499260d0,
# 0.606035270d0, 0.710025546d0, 0.864539937d0, 1.009890014d0,
# 1.156589498d0, 1.268264088d0, 1.408793136d0, 1.528293889d0,
# 1.642389803d0, 1.771929787d0, 1.875426734d0, 1.979445634d0,
# 3.092059170d0, 4.047318673d0, 4.937427406d0, 5.684883866d0,
# 6.473364254d0, 7.141842873d0, 7.823212681d0, 8.484967861d0,
# 9.214907909d0, 9.784981877d0, 10.383431544d0, 11.054527450d0,
# 11.588357911d0, 12.165403348d0, 12.507103232d0, 12.848157153d0,
# 12.975158246d0, 12.770167330d0, 11.991429721d0, 11.123821384d0,
# 10.241165300d0, 9.157568585d0, 7.611314332d0, 5.499567320d0,
# 2.603705830d0, -0.579201122d0, -4.348429690d0, -8.823492981d0,
# -14.416229606d0,-20.848895478d0,-28.011934618d0,-36.077657061d0,
# -43.726529167d0,-50.682419374d0,-56.593917330d0,-60.848713817d0,
# -61.464423821d0,-62.283661600d0,-62.930084228d0,-63.549395305d0,
# -64.133867221d0,-64.465495513d0,-64.909336911d0,-65.479678788d0,
# -65.703953492d0,-66.113782243d0,-66.222507344d0,-66.278786408d0,
# -66.393237110d0,-66.366867383d0,-66.202199439d0,-66.029041364d0,
# -65.774112798d0,-65.485847916d0,-65.247943137d0,-64.801055924d0,
# -60.101647018d0,-54.084245796d0,-47.571709079d0,-41.661629810d0,
# -36.168693788d0,-31.212363345d0,-27.000704767d0,-23.605327875d0,
# -20.580041179d0,-18.274370612d0,-16.434309906d0,-14.327218251d0,
# -12.021979554d0,-10.331796926d0, -8.883031936d0, -7.348435651d0,
# -5.780362152d0, -4.953756827d0, -4.070387325d0, -2.840184650d0,
# -1.678475825d0, -0.679391004d0, 0.223829200d0, 1.086204496d0,
# 2.026633865d0, 2.923942775d0, 3.306978920d0, 3.579966278d0,
# 3.863101723d0, 4.490999475d0, 5.727923150d0, 6.829991246d0,
# 7.212458092d0, 7.551889308d0/
*
IF(top.eq.-1) THEN
n= 226
DO l=1,226
x(l)= xl(l)
y(l)= yl(l)
c write(35,*)l,x(l),y(l)
ENDDO
ELSEIF(top.eq.0) THEN
n= 223
DO l=1,223
x(l)= xc(l)
y(l)= yc(l)
c write(36,*)l,x(l),y(l)
ENDDO
ELSEIF(top.eq.1) THEN
n= 225
DO l=1,225
x(l)= xu(l)
y(l)= yu(l)
c write(37,*)l,x(l),y(l)
ENDDO
ENDIF
*.....Set up tridiagonal system.........................................
* b=diagonal, d=offdiagonal, c=right-hand side
*
d(1)= x(2)-x(1)
c(2)= (y(2)-y(1))/d(1)
DO k= 2,n-1
d(k)= x(k+1)-x(k)
b(k)= TWO*(d(k-1)+d(k))
c(k+1)= (y(k+1)-y(k))/d(k)
c(k)= c(k+1)-c(k)
END DO
*
*.....End conditions. third derivatives at x(1) and x(n) obtained
* from divided differences.......................................
*
b(1)= -d(1)
b(n)= -d(n-1)
c(1)= ZERO
c(n)= ZERO
IF (n.gt.3) THEN
c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
END IF
*
DO k=2,n ! forward elimination
t= d(k-1)/b(k-1)
b(k)= b(k)-t*d(k-1)
c(k)= c(k)-t*c(k-1)
END DO
*
c(n)= c(n)/b(n)
*
* back substitution ( makes c the sigma of text)
*
DO k=n-1,1,-1
c(k)= (c(k)-d(k)*c(k+1))/b(k)
END DO
*
*.....Compute polynomial coefficients...................................
*
b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
DO k=1,n-1
b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
d(k)= (c(k+1)-c(k))/d(k)
c(k)= THREE*c(k)
END DO
c(n)= THREE*c(n)
d(n)= d(n-1)
*
RETURN
*
END
*
*------------------------------------------------------------------------
*
SUBROUTINE Seval3Single0(u,b,c,d,top,gdim,f,fp,fpp,fppp)
*
* ---------------------------------------------------------------------------
*
* PURPOSE - Evaluate the cubic spline function
* Seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
* where x(i) <= u < x(i+1)
*
* NOTES- if u<x(1), i=1 is used;if u>x(n), i=n is used
*
REAL*8 u
* abscissa at which the spline is to be evaluated
*
INTEGER j,k,n,l,top,gdim
*
REAL*8 xl(226),yl(226)
REAL*8 xc(223),yc(223)
REAL*8 xu(225),yu(225)
REAL*8 x(226),y(226)
REAL*8 b(gdim),c(gdim),d(gdim)
* linear,quadratic,cubic coeff
*
REAL*8 f,fp,fpp,fppp
* function, 1st,2nd,3rd deriv
*
INTEGER i
REAL*8 dx
REAL*8 ZERO, TWO, THREE
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
SAVE i
DATA i/1/
*
* The grid
*
DATA (xl(i),i=1,226)/100.0000d0,105.0000d0,110.0000d0,
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,179.0000d0,
# 180.0000d0,181.0000d0,182.0000d0,183.0000d0,184.0000d0,
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
# 340.5000d0,340.8000d0,341.0000d0,341.2000d0,341.4000d0,
# 341.6000d0,341.7000d0,341.8000d0,341.8500d0,341.9000d0,
# 341.9500d0,342.0000d0,342.0200d0,342.0400d0,342.0600d0,
# 342.0800d0,342.1000d0,342.1100d0,342.1120d0,342.1140d0,
# 342.1160d0,342.1180d0,342.1190d0,342.1196d0,342.1198d0,
# 342.1199d0,342.1204d0,342.1220d0,342.2000d0,342.3000d0,
# 342.4000d0,342.5000d0,342.6000d0,342.7000d0,342.8000d0,
# 342.9000d0,343.0000d0,343.2000d0,343.4000d0,343.6000d0,
# 343.8000d0,344.0000d0,345.0000d0,346.0000d0,347.0000d0,
# 348.0000d0,349.0000d0,350.0000d0,351.0000d0,352.0000d0,
# 353.0000d0,354.0000d0,355.0000d0,356.0000d0,357.0000d0,
# 358.0000d0,359.0000d0,360.0000d0,370.0000d0,380.0000d0,
# 390.0000d0,400.0000d0,410.0000d0,420.0000d0,430.0000d0,
# 440.0000d0,450.0000d0,460.0000d0,470.0000d0,480.0000d0,
# 490.0000d0,500.0000d0,510.0000d0,520.0000d0,530.0000d0,
# 540.0000d0,550.0000d0,560.0000d0,565.0000d0,570.0000d0,
# 575.0000d0,580.0000d0,585.0000d0,590.0000d0,595.0000d0,
# 600.0000d0,605.0000d0,610.0000d0,615.0000d0,620.0000d0,
# 625.0000d0,630.0000d0,635.0000d0,640.0000d0,641.0000d0,
# 642.0000d0,643.0000d0,644.0000d0,645.0000d0,646.0000d0,
# 647.0000d0,648.0000d0,649.0000d0,650.0000d0,651.0000d0,
# 652.0000d0,653.0000d0,654.0000d0,655.0000d0,656.0000d0,
# 657.0000d0,658.0000d0,659.0000d0,660.0000d0,670.0000d0,
# 680.0000d0,690.0000d0,700.0000d0,710.0000d0,720.0000d0,
# 730.0000d0,740.0000d0,750.0000d0,760.0000d0,770.0000d0,
# 780.0000d0,790.0000d0,800.0000d0,810.0000d0,820.0000d0,
# 830.0000d0,840.0000d0,850.0000d0,860.0000d0,870.0000d0,
# 880.0000d0,890.0000d0,900.0000d0,910.0000d0,920.0000d0,
# 930.0000d0,940.0000d0,950.0000d0,960.0000d0,970.0000d0,
# 980.0000d0,990.0000d0,1000.0000d0/
*
DATA (yl(i),i=1,226)/-2.927515523d0,-2.713875277d0,-2.477263056d0,
# -2.215702562d0, -1.926872830d0, -1.607794154d0, -1.254642316d0,
# -0.862759386d0, -0.424192821d0, 0.071393692d0, 0.642180822d0,
# 1.311566449d0, 1.458302972d0, 1.610089467d0, 1.771001228d0,
# 1.958842368d0, 1.980856139d0, 2.003769992d0, 2.027683661d0,
# 2.052701483d0, 2.078930176d0, 2.106476183d0, 2.135441928d0,
# 2.165921345d0, 2.197994437d0, 2.231721476d0, 2.323417121d0,
# 2.425109533d0, 2.534671123d0, 2.648366651d0, 2.869451789d0,
# 3.056555683d0, 3.197055532d0, 3.294656543d0, 3.399442397d0,
# 3.438265245d0, 3.450036365d0, 3.449871888d0, 3.442815310d0,
# 3.428262691d0, 3.411797242d0, 3.363961271d0, 3.295179282d0,
# 3.287861332d0, 3.325219720d0, 3.470009593d0, 3.737173473d0,
# 4.021846379d0, 4.241836705d0, 4.394931599d0, 4.501421364d0,
# 4.577551451d0, 4.633363121d0, 4.675411596d0, 4.707771996d0,
# 4.732593619d0, 4.751501682d0, 4.766245197d0, 4.778023555d0,
# 4.787199047d0, 4.793849566d0, 4.798265576d0, 4.801110206d0,
# 4.802711433d0, 4.767182955d0, 4.677220707d0, 4.558671105d0,
# 4.423406368d0, 4.276574444d0, 4.119375511d0, 3.974083777d0,
# 3.819733346d0, 3.671217773d0, 3.517370222d0, 3.344189966d0,
# 3.116802834d0, 2.754066836d0, 1.809839704d0, 1.691727991d0,
# 1.607871864d0, 1.544493194d0, 1.473234377d0, 1.390992309d0,
# 1.292466985d0, 1.234537350d0, 1.167169430d0, 1.128863846d0,
# 1.085963030d0, 1.037416400d0, 0.979461824d0, 0.952439442d0,
# 0.922628624d0, 0.888386832d0, 0.846809865d0, 0.791742524d0,
# 0.752347274d0, 0.742175952d0, 0.730422928d0, 0.716838887d0,
# 0.697956938d0, 0.684261787d0, 0.672163545d0, 0.665924128d0,
# 0.661453423d0, 0.652439216d0, 0.656998002d0, 0.677275168d0,
# 0.699459467d0, 0.721833458d0, 0.738243489d0, 0.751895487d0,
# 0.769698326d0, 0.779779847d0, 0.799512664d0, 0.816893187d0,
# 0.842794338d0, 0.874519045d0, 0.916496501d0, 0.950264630d0,
# 0.968728477d0, 1.105576625d0, 1.278360270d0, 1.391736231d0,
# 1.542276157d0, 1.675869257d0, 1.795595314d0, 1.923303286d0,
# 2.025352084d0, 2.157383377d0, 2.260417650d0, 2.361640519d0,
# 2.479888170d0, 2.572642417d0, 2.688809606d0, 2.796999182d0,
# 2.879392617d0, 3.848596335d0, 4.732837134d0, 5.525766484d0,
# 6.243989288d0, 6.975830059d0, 7.607337162d0, 8.167782071d0,
# 8.790374186d0, 9.422028844d0, 9.922844338d0, 10.426856521d0,
# 10.933372206d0, 11.329434560d0, 11.735173309d0, 11.795158880d0,
# 11.720854043d0, 11.455441492d0, 10.516061928d0, 9.017310079d0,
# 7.000279246d0, 5.560369962d0, 3.715234200d0, 1.373399680d0,
# -1.698051464d0, -5.394210410d0, -9.477907268d0,-13.677692505d0,
# -18.256012300d0,-24.029072069d0,-29.992304980d0,-35.639341296d0,
# -41.517060283d0,-46.430415213d0,-50.471411203d0,-53.247623713d0,
# -55.184140125d0,-55.523807703d0,-55.749322940d0,-55.785869793d0,
# -56.086101947d0,-56.289055861d0,-56.369321707d0,-56.555818993d0,
# -56.632365789d0,-56.616832544d0,-56.709629986d0,-56.539655434d0,
# -56.366567734d0,-56.272247546d0,-56.170833793d0,-55.984866548d0,
# -55.776968735d0,-55.354292076d0,-55.136594171d0,-54.856130174d0,
# -54.390381982d0,-50.585126291d0,-46.161103219d0,-41.187366158d0,
# -36.554683053d0,-32.185245783d0,-28.097939638d0,-24.608441494d0,
# -21.643195994d0,-18.984106588d0,-16.998672281d0,-15.273879570d0,
# -13.373334219d0,-11.294832079d0, -9.584938331d0, -8.286389135d0,
# -6.838755740d0, -5.409090746d0, -4.588339018d0, -3.705250749d0,
# -2.476432078d0, -1.367536281d0, -0.414988031d0, 0.510347996d0,
# 1.376925204d0, 2.276754999d0, 3.203065379d0, 3.629617497d0,
# 3.821185286d0, 4.122178708d0, 4.804270233d0, 5.929961061d0,
# 7.072311664d0, 7.458878393d0, 7.823110180d0/
*
DATA (xc(i),i=1,223)/100.0000d0,105.0000d0,110.0000d0,
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,178.5000d0,
# 179.0000d0,179.5000d0,180.0000d0,180.5000d0,181.0000d0,
# 181.5000d0,182.0000d0,182.5000d0,183.0000d0,184.0000d0,
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
# 341.0000d0,342.0000d0,343.0000d0,344.0000d0,344.2000d0,
# 344.4000d0,344.6000d0,344.8000d0,344.9000d0,345.0000d0,
# 345.0500d0,345.1000d0,345.1500d0,345.2000d0,345.2200d0,
# 345.2400d0,345.2600d0,345.2700d0,345.2720d0,345.2740d0,
# 345.2760d0,345.2780d0,345.2790d0,345.2796d0,345.2798d0,
# 345.2799d0,345.2804d0,345.2820d0,345.2900d0,345.3500d0,
# 345.4500d0,345.5200d0,345.6000d0,345.7000d0,345.8000d0,
# 346.0000d0,347.0000d0,348.0000d0,349.0000d0,350.0000d0,
# 351.0000d0,352.0000d0,353.0000d0,354.0000d0,355.0000d0,
# 356.0000d0,357.0000d0,358.0000d0,359.0000d0,360.0000d0,
# 370.0000d0,380.0000d0,390.0000d0,400.0000d0,410.0000d0,
# 420.0000d0,430.0000d0,440.0000d0,450.0000d0,460.0000d0,
# 470.0000d0,480.0000d0,490.0000d0,500.0000d0,510.0000d0,
# 520.0000d0,530.0000d0,540.0000d0,550.0000d0,560.0000d0,
# 565.0000d0,570.0000d0,575.0000d0,580.0000d0,585.0000d0,
# 590.0000d0,595.0000d0,600.0000d0,605.0000d0,610.0000d0,
# 615.0000d0,620.0000d0,625.0000d0,630.0000d0,635.0000d0,
# 640.0000d0,641.0000d0,642.0000d0,643.0000d0,644.0000d0,
# 645.0000d0,646.0000d0,647.0000d0,648.0000d0,649.0000d0,
# 650.0000d0,651.0000d0,652.0000d0,653.0000d0,654.0000d0,
# 655.0000d0,656.0000d0,657.0000d0,658.0000d0,659.0000d0,
# 660.0000d0,670.0000d0,680.0000d0,690.0000d0,700.0000d0,
# 710.0000d0,720.0000d0,730.0000d0,740.0000d0,750.0000d0,
# 760.0000d0,770.0000d0,780.0000d0,790.0000d0,800.0000d0,
# 810.0000d0,820.0000d0,830.0000d0,840.0000d0,850.0000d0,
# 860.0000d0,870.0000d0,880.0000d0,890.0000d0,900.0000d0,
# 910.0000d0,920.0000d0,930.0000d0,940.0000d0,950.0000d0,
# 960.0000d0,970.0000d0,980.0000d0,990.0000d0,1000.0000d0/
*
DATA (yc(i),i=1,223)/-2.969458664d0,-2.755568498d0,-2.518689112d0,
# -2.256840599d0, -1.967703151d0, -1.648299088d0, -1.294794285d0,
# -0.902525232d0, -0.463541061d0, 0.032517848d0, 0.603841958d0,
# 1.273887882d0, 1.420789352d0, 1.572747980d0, 1.733843353d0,
# 1.921879203d0, 1.943912705d0, 1.966846336d0, 1.990779880d0,
# 2.015817728d0, 2.042066621d0, 2.069633010d0, 2.098619322d0,
# 2.129119489d0, 2.161213532d0, 2.194961744d0, 2.286711412d0,
# 2.388459445d0, 2.498078035d0, 2.611830784d0, 2.833024715d0,
# 3.020227220d0, 3.160820150d0, 3.258502948d0, 3.363381747d0,
# 3.402305539d0, 3.414208425d0, 3.414049050d0, 3.407053730d0,
# 3.392629803d0, 3.376195597d0, 3.328499111d0, 3.268651902d0,
# 3.259860594d0, 3.253486002d0, 3.252546948d0, 3.262154622d0,
# 3.289924633d0, 3.345085115d0, 3.434732330d0, 3.557729513d0,
# 3.701897111d0, 3.986569984d0, 4.206557658d0, 4.359682402d0,
# 4.466176116d0, 4.542300703d0, 4.598184529d0, 4.640273713d0,
# 4.672604523d0, 4.697438255d0, 4.716361927d0, 4.731103154d0,
# 4.742925539d0, 4.752092918d0, 4.758670785d0, 4.763130070d0,
# 4.766074843d0, 4.767748196d0, 4.732269356d0, 4.642413616d0,
# 4.524181528d0, 4.389510854d0, 4.242822112d0, 4.086217524d0,
# 3.941680893d0, 3.789691168d0, 3.644241595d0, 3.494351920d0,
# 3.331867639d0, 3.122929519d0, 2.816374576d0, 2.186038507d0,
# 2.068538445d0, 1.926511260d0, 1.746359564d0, 1.494945494d0,
# 1.428923128d0, 1.354291504d0, 1.267892764d0, 1.163289817d0,
# 1.100587592d0, 1.027248191d0, 0.984553881d0, 0.936444199d0,
# 0.879742362d0, 0.808975906d0, 0.774145495d0, 0.731724792d0,
# 0.675557165d0, 0.635406715d0, 0.625004410d0, 0.613030568d0,
# 0.599164524d0, 0.579800865d0, 0.565857087d0, 0.553503488d0,
# 0.547057831d0, 0.542411799d0, 0.533743538d0, 0.538671942d0,
# 0.542814415d0, 0.562187598d0, 0.575771207d0, 0.586370628d0,
# 0.605546923d0, 0.622661850d0, 0.639021282d0, 0.683852301d0,
# 0.846885709d0, 0.995547199d0, 1.139428677d0, 1.270752033d0,
# 1.403862461d0, 1.540060220d0, 1.676615197d0, 1.791154110d0,
# 1.909977528d0, 2.023870438d0, 2.138393119d0, 2.233199625d0,
# 2.358706648d0, 2.446271802d0, 3.459751114d0, 4.411864477d0,
# 5.220622691d0, 5.988967245d0, 6.706065217d0, 7.359912016d0,
# 7.993290467d0, 8.675262256d0, 9.327550412d0, 9.917129519d0,
# 10.373037136d0, 11.009210530d0, 11.516244354d0, 11.985967329d0,
# 12.198319618d0, 12.357011718d0, 12.311558072d0, 11.754465713d0,
# 10.617092695d0, 9.043395280d0, 8.002055503d0, 6.601203443d0,
# 4.587232698d0, 1.906399387d0, -1.281177539d0, -5.052281649d0,
# -9.213421228d0,-13.868726558d0,-19.489194521d0,-25.988103189d0,
# -32.734525726d0,-39.402198052d0,-45.903352541d0,-51.189352305d0,
# -55.322071720d0,-58.186905136d0,-58.653812374d0,-59.048160142d0,
# -59.279596622d0,-59.992120796d0,-60.271154555d0,-60.549383981d0,
# -60.729563975d0,-60.954540486d0,-61.084961955d0,-61.164811836d0,
# -61.186911575d0,-61.109847454d0,-61.108201824d0,-60.992397754d0,
# -60.768938455d0,-60.645603901d0,-60.224362756d0,-59.860911894d0,
# -59.548720101d0,-59.301307840d0,-55.034929471d0,-49.745997496d0,
# -44.181440809d0,-39.061582838d0,-34.158470096d0,-29.603188102d0,
# -25.799341483d0,-22.634533844d0,-19.804942358d0,-17.654326276d0,
# -15.913630982d0,-13.830699918d0,-11.633319645d0,-10.006868842d0,
# -8.532937676d0, -7.132846889d0, -5.627566626d0, -4.758444504d0,
# -3.887401469d0, -2.729580751d0, -1.512994295d0, -0.550651909d0,
# 0.355022046d0, 1.248689020d0, 2.123963438d0, 3.014801413d0,
# 3.484406765d0, 3.696913629d0, 4.008426964d0, 4.651255618d0,
# 5.840647264d0, 6.935940541d0, 7.364646872d0, 7.670905366d0/
*
DATA (xu(i),i=1,225)/100.0000d0,105.0000d0,110.0000d0,
# 115.0000d0,120.0000d0,125.0000d0,130.0000d0,135.0000d0,
# 140.0000d0,145.0000d0,150.0000d0,155.0000d0,156.0000d0,
# 157.0000d0,158.0000d0,159.0000d0,159.1000d0,159.2000d0,
# 159.3000d0,159.4000d0,159.5000d0,159.6000d0,159.7000d0,
# 159.8000d0,159.9000d0,160.0000d0,160.2500d0,160.5000d0,
# 160.7500d0,161.0000d0,161.5000d0,162.0000d0,162.5000d0,
# 163.0000d0,164.0000d0,165.0000d0,166.0000d0,167.0000d0,
# 168.0000d0,169.0000d0,170.0000d0,175.0000d0,178.5000d0,
# 179.0000d0,179.5000d0,180.0000d0,180.5000d0,181.0000d0,
# 181.5000d0,182.0000d0,182.5000d0,183.0000d0,184.0000d0,
# 185.0000d0,186.0000d0,187.0000d0,188.0000d0,189.0000d0,
# 190.0000d0,191.0000d0,192.0000d0,193.0000d0,194.0000d0,
# 195.0000d0,196.0000d0,197.0000d0,198.0000d0,199.0000d0,
# 200.0000d0,210.0000d0,220.0000d0,230.0000d0,240.0000d0,
# 250.0000d0,260.0000d0,270.0000d0,280.0000d0,290.0000d0,
# 300.0000d0,310.0000d0,320.0000d0,330.0000d0,340.0000d0,
# 341.0000d0,342.0000d0,343.0000d0,344.0000d0,345.0000d0,
# 345.5000d0,345.8000d0,346.0000d0,347.0000d0,347.2000d0,
# 347.4000d0,347.6000d0,347.8000d0,348.0000d0,348.1000d0,
# 348.2000d0,348.2500d0,348.3000d0,348.3200d0,348.3400d0,
# 348.3600d0,348.3800d0,348.4000d0,348.4100d0,348.4200d0,
# 348.4300d0,348.4320d0,348.4340d0,348.4360d0,348.4380d0,
# 348.4390d0,348.4396d0,348.4398d0,348.4399d0,348.4404d0,
# 348.4420d0,348.4500d0,348.5000d0,348.7000d0,349.0000d0,
# 349.4000d0,350.0000d0,351.0000d0,352.0000d0,353.0000d0,
# 354.0000d0,355.0000d0,356.0000d0,357.0000d0,358.0000d0,
# 359.0000d0,360.0000d0,370.0000d0,380.0000d0,390.0000d0,
# 400.0000d0,410.0000d0,420.0000d0,430.0000d0,440.0000d0,
# 450.0000d0,460.0000d0,470.0000d0,480.0000d0,490.0000d0,
# 500.0000d0,510.0000d0,520.0000d0,530.0000d0,540.0000d0,
# 550.0000d0,560.0000d0,565.0000d0,570.0000d0,575.0000d0,
# 580.0000d0,585.0000d0,590.0000d0,595.0000d0,600.0000d0,
# 605.0000d0,610.0000d0,615.0000d0,620.0000d0,625.0000d0,
# 630.0000d0,635.0000d0,640.0000d0,641.0000d0,642.0000d0,
# 643.0000d0,644.0000d0,645.0000d0,646.0000d0,647.0000d0,
# 648.0000d0,649.0000d0,650.0000d0,651.0000d0,652.0000d0,
# 653.0000d0,654.0000d0,655.0000d0,656.0000d0,657.0000d0,
# 658.0000d0,659.0000d0,660.0000d0,670.0000d0,680.0000d0,
# 690.0000d0,700.0000d0,710.0000d0,720.0000d0,730.0000d0,
# 740.0000d0,750.0000d0,760.0000d0,770.0000d0,780.0000d0,
# 790.0000d0,800.0000d0,810.0000d0,820.0000d0,830.0000d0,
# 840.0000d0,850.0000d0,860.0000d0,870.0000d0,880.0000d0,
# 890.0000d0,900.0000d0,910.0000d0,920.0000d0,930.0000d0,
# 940.0000d0,950.0000d0,960.0000d0,970.0000d0,980.0000d0,
# 990.0000d0,1000.0000d0/
*
DATA (yu(i),i=1,225)/-3.011799533d0,-2.797658814d0,-2.560511342d0,
# -2.298374579d0, -2.008927458d0, -1.689194746d0, -1.335334782d0,
# -0.942678003d0, -0.503267763d0, -0.006741983d0, 0.565124142d0,
# 1.235840116d0, 1.382905273d0, 1.535044791d0, 1.696328604d0,
# 1.884567714d0, 1.906622347d0, 1.929577239d0, 1.953532153d0,
# 1.978591436d0, 2.004861792d0, 2.032449634d0, 2.061457355d0,
# 2.091978868d0, 2.124094156d0, 2.157863499d0, 2.249665540d0,
# 2.351465763d0, 2.461136573d0, 2.574941494d0, 2.796240093d0,
# 2.983545301d0, 3.124228166d0, 3.221982342d0, 3.326972382d0,
# 3.365993805d0, 3.378023193d0, 3.377887062d0, 3.370931083d0,
# 3.356634680d0, 3.340265619d0, 3.292682936d0, 3.232961811d0,
# 3.224194020d0, 3.217828260d0, 3.216889093d0, 3.226496714d0,
# 3.254271720d0, 3.309440914d0, 3.399097855d0, 3.522103315d0,
# 3.666276367d0, 3.950949463d0, 4.170926862d0, 4.324070476d0,
# 4.430588837d0, 4.506706461d0, 4.562622634d0, 4.604788444d0,
# 4.637112788d0, 4.661928180d0, 4.680880990d0, 4.695625131d0,
# 4.707459275d0, 4.716670603d0, 4.723178694d0, 4.727608289d0,
# 4.730650318d0, 4.732442748d0, 4.697030978d0, 4.607168713d0,
# 4.489397175d0, 4.355888392d0, 4.208987612d0, 4.052707904d0,
# 3.908019896d0, 3.759968739d0, 3.615348691d0, 3.469932507d0,
# 3.315332978d0, 3.122286809d0, 2.855004727d0, 2.377722123d0,
# 2.300655190d0, 2.213908305d0, 2.115631830d0, 2.000162289d0,
# 1.860840425d0, 1.779119559d0, 1.724912832d0, 1.685992661d0,
# 1.446048759d0, 1.384730617d0, 1.316264190d0, 1.238304389d0,
# 1.147571101d0, 1.036738128d0, 0.969143795d0, 0.888640045d0,
# 0.841031364d0, 0.785526212d0, 0.759880178d0, 0.731970096d0,
# 0.701136254d0, 0.665654946d0, 0.622492994d0, 0.596679361d0,
# 0.565315972d0, 0.524528227d0, 0.513928834d0, 0.501826683d0,
# 0.487735928d0, 0.468019598d0, 0.454013349d0, 0.441470198d0,
# 0.434889761d0, 0.430148561d0, 0.421464835d0, 0.429200773d0,
# 0.429255984d0, 0.437861500d0, 0.480919141d0, 0.544499260d0,
# 0.606035270d0, 0.710025546d0, 0.864539937d0, 1.009890014d0,
# 1.156589498d0, 1.268264088d0, 1.408793136d0, 1.528293889d0,
# 1.642389803d0, 1.771929787d0, 1.875426734d0, 1.979445634d0,
# 3.092059170d0, 4.047318673d0, 4.937427406d0, 5.684883866d0,
# 6.473364254d0, 7.141842873d0, 7.823212681d0, 8.484967861d0,
# 9.214907909d0, 9.784981877d0, 10.383431544d0, 11.054527450d0,
# 11.588357911d0, 12.165403348d0, 12.507103232d0, 12.848157153d0,
# 12.975158246d0, 12.770167330d0, 11.991429721d0, 11.123821384d0,
# 10.241165300d0, 9.157568585d0, 7.611314332d0, 5.499567320d0,
# 2.603705830d0, -0.579201122d0, -4.348429690d0, -8.823492981d0,
# -14.416229606d0,-20.848895478d0,-28.011934618d0,-36.077657061d0,
# -43.726529167d0,-50.682419374d0,-56.593917330d0,-60.848713817d0,
# -61.464423821d0,-62.283661600d0,-62.930084228d0,-63.549395305d0,
# -64.133867221d0,-64.465495513d0,-64.909336911d0,-65.479678788d0,
# -65.703953492d0,-66.113782243d0,-66.222507344d0,-66.278786408d0,
# -66.393237110d0,-66.366867383d0,-66.202199439d0,-66.029041364d0,
# -65.774112798d0,-65.485847916d0,-65.247943137d0,-64.801055924d0,
# -60.101647018d0,-54.084245796d0,-47.571709079d0,-41.661629810d0,
# -36.168693788d0,-31.212363345d0,-27.000704767d0,-23.605327875d0,
# -20.580041179d0,-18.274370612d0,-16.434309906d0,-14.327218251d0,
# -12.021979554d0,-10.331796926d0, -8.883031936d0, -7.348435651d0,
# -5.780362152d0, -4.953756827d0, -4.070387325d0, -2.840184650d0,
# -1.678475825d0, -0.679391004d0, 0.223829200d0, 1.086204496d0,
# 2.026633865d0, 2.923942775d0, 3.306978920d0, 3.579966278d0,
# 3.863101723d0, 4.490999475d0, 5.727923150d0, 6.829991246d0,
# 7.212458092d0, 7.551889308d0/
*
IF(top.eq.-1) THEN
n= 226
DO l=1,226
x(l)= xl(l)
y(l)= yl(l)
c write(35,*)l,x(l),y(l)
ENDDO
ELSEIF(top.eq.0) THEN
n= 223
DO l=1,223
x(l)= xc(l)
y(l)= yc(l)
c write(36,*)l,x(l),y(l)
ENDDO
ELSEIF(top.eq.1) THEN
n= 225
DO l=1,225
x(l)= xu(l)
y(l)= yu(l)
c write(37,*)l,x(l),y(l)
ENDDO
ENDIF
*
*.....First check if u is in the same interval found on the
* last call to Seval.............................................
*
IF ( (i.lt.1) .OR. (i.ge.n) ) i=1
IF ( (u.lt.x(i)) .OR. (u.ge.x(i+1)) ) THEN
i=1
*
* binary search
*
j= n+1
2468 k= (i+j)/2
IF (u.lt.x(k)) THEN
j= k
ELSE
i= k
ENDIF
IF (j.le.i+1) GOTO 2469
GOTO 2468
2469 CONTINUE
ENDIF
*
dx= u-x(i)
*
* evaluate the spline
*
f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i))
fpp= TWO*c(i) + dx*SIX*d(i)
fppp= SIX*d(i)
*
RETURN
*
END
*****************************************************************
C---electroweak corrections to H --> gg in SM4
C m_t' = m_b' + 50*(1+1/5*log(M_H/(115 GeV)) GeV
*****************************************************************
double precision function glgl_elw4(im,mh)
*
IMPLICIT NONE
*
INTEGER i,top,gdim,im
REAL*8 u,ui,value
REAL*8 vp,vpp,vppp
REAL*8 mtop,mh
REAL*8 x,x0,x1,x2,y0,y1,y2,a0,a1,a2
REAL*8 value0,value1,value2
REAL*8 bl(94),cl(94),dl(94)
REAL*8 bc(120),cc(120),dc(120)
*
* u value of M_H at which the spline is to be evaluated
* top= -1,0,1 lower, central, upper value for m_top
*
u = mh
if(im.eq.1)then
gdim= 94
CALL F4MMsplineSingle1(bl,cl,dl,gdim)
CALL S4eval3Single1(u,bl,cl,dl,gdim,value,vp,vpp,vppp)
else
gdim= 120
CALL F4MMsplineSingle2(bc,cc,dc,gdim)
CALL S4eval3Single2(u,bc,cc,dc,gdim,value,vp,vpp,vppp)
endif
glgl_elw4=value/100.d0
RETURN
END
*
*-----------------------------------------------------------------------
*
c CONTAINS
*
SUBROUTINE F4MMsplineSingle1(b,c,d,gdim)
*
* ---------------------------------------------------------------------------
*
* PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline
* so that the interpolated value is given by
* s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3
* when x(k) <= x <= x(k+1)
* The end conditions match the third derivatives of the interpolated curve to
* the third derivatives of the unique polynomials thru the first four and
* last four points.
* Use Seval or Seval3 to evaluate the spline.
*
INTEGER k,n,i,top,gdim,l
*
REAL*8 xc(94),yc(94)
REAL*8 x(94),y(94)
*
REAL*8 b(gdim)
* linear coeff
*
REAL*8 c(gdim)
* quadratic coeff.
*
REAL*8 d(gdim)
* cubic coeff.
*
REAL*8 t
REAL*8 ZERO, TWO, THREE
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
*
* The grid
*
DATA (xc(i),i=1,94)/
.100.000000000, 110.000000000, 120.000000000, 130.000000000,
.140.000000000, 145.000000000, 150.000000000, 151.000000000,
.152.000000000, 153.000000000, 154.000000000, 155.000000000,
.156.000000000, 157.000000000, 158.000000000, 159.000000000,
.160.000000000, 161.000000000, 162.000000000, 163.000000000,
.164.000000000, 165.000000000, 166.000000000, 167.000000000,
.168.000000000, 169.000000000, 170.000000000, 171.000000000,
.172.000000000, 173.000000000, 174.000000000, 175.000000000,
.176.000000000, 177.000000000, 178.000000000, 179.000000000,
.180.000000000, 181.000000000, 182.000000000, 183.000000000,
.184.000000000, 185.000000000, 186.000000000, 187.000000000,
.188.000000000, 189.000000000, 190.000000000, 195.000000000,
.200.000000000, 210.000000000, 220.000000000, 230.000000000,
.240.000000000, 250.000000000, 260.000000000, 270.000000000,
.280.000000000, 290.000000000, 300.000000000, 310.000000000,
.325.000000000, 330.000000000, 335.000000000, 337.000000000,
.339.000000000, 341.000000000, 342.000000000, 343.000000000,
.344.000000000, 345.000000000, 346.000000000, 347.000000000,
.348.000000000, 349.000000000, 350.000000000, 355.000000000,
.360.000000000, 365.000000000, 370.000000000, 380.000000000,
.390.000000000, 395.000000000, 400.000000000, 420.000000000,
.430.000000000, 440.000000000, 450.000000000, 460.000000000,
.470.000000000, 480.000000000, 490.000000000, 500.000000000,
.550.000000000, 600.000000000/
*
DATA (yc(i),i=1,94)/
. 7.081784858, 7.005040906, 6.909638450, 6.767663371,
. 6.553481059, 6.394170663, 6.161453100, 6.100169493,
. 6.031159987, 5.952322286, 5.860793898, 5.752611211,
. 5.622148850, 5.462299560, 5.267015468, 5.045635485,
. 4.867448493, 4.858996855, 4.957159679, 4.991788742,
. 4.949709854, 4.869086008, 4.773624254, 4.672977010,
. 4.574106306, 4.473957550, 4.375216682, 4.277835718,
. 4.180106440, 4.082355068, 3.981303187, 3.876580238,
. 3.765612293, 3.650290738, 3.517994735, 3.374532294,
. 3.224190879, 3.081402574, 3.003176187, 3.023510652,
. 3.058333994, 3.061527486, 3.032524574, 2.980103496,
. 2.916221523, 2.850606479, 2.785863031, 2.465569133,
. 2.204423840, 1.758167155, 1.356606770, 1.009560117,
. 0.701093096, 0.393933499, 0.079112288, -0.206013527,
.-0.518586518, -0.811157249, -1.129144131, -1.439334674,
.-1.957798896, -2.212619484, -2.441510248, -2.573014363,
.-2.728219541, -2.899530309, -2.993431692, -3.158243996,
.-3.351666868, -3.920928883, -4.098697186, -4.132701122,
.-4.144213562, -4.088523452, -4.049143398, -3.913706462,
.-3.749072205, -3.669404489, -3.618383929, -3.575928277,
.-3.688589851, -3.744277795, -3.841702830, -4.446704386,
.-4.761423299, -5.287974189, -5.744301071, -6.211157195,
.-6.778672157, -7.445798763, -8.061705151, -8.714706599,
.-12.512297879, -16.999567032/
*
n= 94
DO l=1,n
x(l)= xc(l)
y(l)= yc(l)
ENDDO
*.....Set up tridiagonal system.........................................
* b=diagonal, d=offdiagonal, c=right-hand side
*
d(1)= x(2)-x(1)
c(2)= (y(2)-y(1))/d(1)
DO k= 2,n-1
d(k)= x(k+1)-x(k)
b(k)= TWO*(d(k-1)+d(k))
c(k+1)= (y(k+1)-y(k))/d(k)
c(k)= c(k+1)-c(k)
END DO
*
*.....End conditions. third derivatives at x(1) and x(n) obtained
* from divided differences.......................................
*
b(1)= -d(1)
b(n)= -d(n-1)
c(1)= ZERO
c(n)= ZERO
IF (n.gt.3) THEN
c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
END IF
*
DO k=2,n ! forward elimination
t= d(k-1)/b(k-1)
b(k)= b(k)-t*d(k-1)
c(k)= c(k)-t*c(k-1)
END DO
*
c(n)= c(n)/b(n)
*
* back substitution ( makes c the sigma of text)
*
DO k=n-1,1,-1
c(k)= (c(k)-d(k)*c(k+1))/b(k)
END DO
*
*.....Compute polynomial coefficients...................................
*
b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
DO k=1,n-1
b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
d(k)= (c(k+1)-c(k))/d(k)
c(k)= THREE*c(k)
END DO
c(n)= THREE*c(n)
d(n)= d(n-1)
*
RETURN
*
END
*
*------------------------------------------------------------------------
*
SUBROUTINE S4eval3Single1(u,b,c,d,gdim,f,fp,fpp,fppp)
*
* ---------------------------------------------------------------------------
*
* PURPOSE - Evaluate the cubic spline function
* Seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
* where x(i) <= u < x(i+1)
*
* NOTES- if u<x(1), i=1 is used;if u>x(n), i=n is used
*
REAL*8 u
* abscissa at which the spline is to be evaluated
*
INTEGER j,k,n,l,top,gdim
*
REAL*8 xc(94),yc(94)
REAL*8 x(94),y(94)
REAL*8 b(gdim),c(gdim),d(gdim)
* linear,quadratic,cubic coeff
*
REAL*8 f,fp,fpp,fppp
* function, 1st,2nd,3rd deriv
*
INTEGER i
REAL*8 dx
REAL*8 ZERO, TWO, THREE
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
SAVE i
DATA i/1/
*
* The grid
*
DATA (xc(l),l=1,94)/
.100.000000000, 110.000000000, 120.000000000, 130.000000000,
.140.000000000, 145.000000000, 150.000000000, 151.000000000,
.152.000000000, 153.000000000, 154.000000000, 155.000000000,
.156.000000000, 157.000000000, 158.000000000, 159.000000000,
.160.000000000, 161.000000000, 162.000000000, 163.000000000,
.164.000000000, 165.000000000, 166.000000000, 167.000000000,
.168.000000000, 169.000000000, 170.000000000, 171.000000000,
.172.000000000, 173.000000000, 174.000000000, 175.000000000,
.176.000000000, 177.000000000, 178.000000000, 179.000000000,
.180.000000000, 181.000000000, 182.000000000, 183.000000000,
.184.000000000, 185.000000000, 186.000000000, 187.000000000,
.188.000000000, 189.000000000, 190.000000000, 195.000000000,
.200.000000000, 210.000000000, 220.000000000, 230.000000000,
.240.000000000, 250.000000000, 260.000000000, 270.000000000,
.280.000000000, 290.000000000, 300.000000000, 310.000000000,
.325.000000000, 330.000000000, 335.000000000, 337.000000000,
.339.000000000, 341.000000000, 342.000000000, 343.000000000,
.344.000000000, 345.000000000, 346.000000000, 347.000000000,
.348.000000000, 349.000000000, 350.000000000, 355.000000000,
.360.000000000, 365.000000000, 370.000000000, 380.000000000,
.390.000000000, 395.000000000, 400.000000000, 420.000000000,
.430.000000000, 440.000000000, 450.000000000, 460.000000000,
.470.000000000, 480.000000000, 490.000000000, 500.000000000,
.550.000000000, 600.000000000/
*
DATA (yc(l),l=1,94)/
. 7.081784858, 7.005040906, 6.909638450, 6.767663371,
. 6.553481059, 6.394170663, 6.161453100, 6.100169493,
. 6.031159987, 5.952322286, 5.860793898, 5.752611211,
. 5.622148850, 5.462299560, 5.267015468, 5.045635485,
. 4.867448493, 4.858996855, 4.957159679, 4.991788742,
. 4.949709854, 4.869086008, 4.773624254, 4.672977010,
. 4.574106306, 4.473957550, 4.375216682, 4.277835718,
. 4.180106440, 4.082355068, 3.981303187, 3.876580238,
. 3.765612293, 3.650290738, 3.517994735, 3.374532294,
. 3.224190879, 3.081402574, 3.003176187, 3.023510652,
. 3.058333994, 3.061527486, 3.032524574, 2.980103496,
. 2.916221523, 2.850606479, 2.785863031, 2.465569133,
. 2.204423840, 1.758167155, 1.356606770, 1.009560117,
. 0.701093096, 0.393933499, 0.079112288, -0.206013527,
.-0.518586518, -0.811157249, -1.129144131, -1.439334674,
.-1.957798896, -2.212619484, -2.441510248, -2.573014363,
.-2.728219541, -2.899530309, -2.993431692, -3.158243996,
.-3.351666868, -3.920928883, -4.098697186, -4.132701122,
.-4.144213562, -4.088523452, -4.049143398, -3.913706462,
.-3.749072205, -3.669404489, -3.618383929, -3.575928277,
.-3.688589851, -3.744277795, -3.841702830, -4.446704386,
.-4.761423299, -5.287974189, -5.744301071, -6.211157195,
.-6.778672157, -7.445798763, -8.061705151, -8.714706599,
.-12.512297879, -16.999567032/
*
n= 94
DO l=1,n
x(l)= xc(l)
y(l)= yc(l)
ENDDO
*
*.....First check if u is in the same interval found on the
* last call to Seval.............................................
*
IF ( (i.lt.1) .OR. (i.ge.n) ) i=1
IF ( (u.lt.x(i)) .OR. (u.ge.x(i+1)) ) THEN
i=1
*
* binary search
*
j= n+1
2468 k= (i+j)/2
IF (u.lt.x(k)) THEN
j= k
ELSE
i= k
ENDIF
IF (j.le.i+1) GOTO 2469
GOTO 2468
2469 CONTINUE
ENDIF
*
dx= u-x(i)
*
* evaluate the spline
*
f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i))
fpp= TWO*c(i) + dx*SIX*d(i)
fppp= SIX*d(i)
*
RETURN
*
END
SUBROUTINE F4MMsplineSingle2(b,c,d,gdim)
*
* ---------------------------------------------------------------------------
*
* PURPOSE - Compute the coefficients b,c,d for a cubic interpolating spline
* so that the interpolated value is given by
* s(x) = y(k) + b(k)*(x-x(k)) + c(k)*(x-x(k))**2 + d(k)*(x-x(k))**3
* when x(k) <= x <= x(k+1)
* The end conditions match the third derivatives of the interpolated curve to
* the third derivatives of the unique polynomials thru the first four and
* last four points.
* Use Seval or Seval3 to evaluate the spline.
*
INTEGER k,n,i,top,gdim,l
*
REAL*8 xc(120),yc(120)
REAL*8 x(120),y(120)
*
REAL*8 b(gdim)
* linear coeff
*
REAL*8 c(gdim)
* quadratic coeff.
*
REAL*8 d(gdim)
* cubic coeff.
*
REAL*8 t
REAL*8 ZERO, TWO, THREE
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
*
* The grid
*
DATA (xc(i),i=1,120)/
.100.000000000, 110.000000000, 120.000000000, 130.000000000,
.140.000000000, 145.000000000, 150.000000000, 151.000000000,
.152.000000000, 153.000000000, 154.000000000, 155.000000000,
.156.000000000, 157.000000000, 158.000000000, 159.000000000,
.160.000000000, 161.000000000, 162.000000000, 163.000000000,
.164.000000000, 165.000000000, 166.000000000, 167.000000000,
.168.000000000, 169.000000000, 170.000000000, 171.000000000,
.172.000000000, 173.000000000, 174.000000000, 175.000000000,
.176.000000000, 177.000000000, 178.000000000, 179.000000000,
.180.000000000, 181.000000000, 182.000000000, 183.000000000,
.184.000000000, 185.000000000, 186.000000000, 187.000000000,
.188.000000000, 189.000000000, 190.000000000, 195.000000000,
.200.000000000, 210.000000000, 220.000000000, 230.000000000,
.240.000000000, 250.000000000, 260.000000000, 270.000000000,
.280.000000000, 290.000000000, 300.000000000, 310.000000000,
.320.000000000, 325.000000000, 330.000000000, 335.000000000,
.336.000000000, 337.000000000, 338.000000000, 339.000000000,
.340.000000000, 341.000000000, 342.000000000, 343.000000000,
.344.000000000, 345.000000000, 346.000000000, 347.000000000,
.348.000000000, 349.000000000, 350.000000000, 351.000000000,
.353.000000000, 355.000000000, 360.000000000, 365.000000000,
.370.000000000, 375.000000000, 380.000000000, 385.000000000,
.390.000000000, 400.000000000, 410.000000000, 420.000000000,
.430.000000000, 440.000000000, 450.000000000, 460.000000000,
.470.000000000, 480.000000000, 490.000000000, 500.000000000,
.550.000000000, 600.000000000, 650.000000000, 700.000000000,
.750.000000000, 800.000000000, 850.000000000, 900.000000000,
.950.000000000, 1000.000000000, 1100.000000000, 1150.000000000,
.1200.000000000, 1250.000000000, 1300.000000000, 1350.000000000,
.1400.000000000, 1600.000000000, 1800.000000000, 2000.000000000/
*
DATA (yc(i),i=1,120)/
.12.283000000, 12.212000000, 12.118000000, 11.982000000,
.11.780000000, 11.620000000, 11.394000000, 11.333000000,
.11.265000000, 11.191000000, 11.099000000, 10.994000000,
.10.862000000, 10.705000000, 10.510000000, 10.290000000,
.10.112000000, 10.104000000, 10.202000000, 10.238000000,
.10.197000000, 10.115000000, 10.019000000, 9.919000000,
. 9.822000000, 9.726000000, 9.632000000, 9.531000000,
. 9.432000000, 9.334000000, 9.234000000, 9.134000000,
. 9.027000000, 8.907000000, 8.786000000, 8.639000000,
. 8.485000000, 8.351000000, 8.278000000, 8.286000000,
. 8.320000000, 8.336000000, 8.300000000, 8.246000000,
. 8.176000000, 8.107000000, 8.057000000, 7.757000000,
. 7.465000000, 7.044000000, 6.720000000, 6.399000000,
. 6.065000000, 5.815000000, 5.600000000, 5.273000000,
. 5.038000000, 4.795000000, 4.541000000, 4.349000000,
. 4.106000000, 3.978000000, 3.847000000, 3.688000000,
. 3.663000000, 3.589000000, 3.539000000, 3.463000000,
. 3.417000000, 3.317000000, 3.223000000, 3.169000000,
. 3.048000000, 2.449000000, 2.361000000, 2.327000000,
. 2.356000000, 2.470000000, 2.570000000, 2.629000000,
. 2.721000000, 2.834000000, 3.049000000, 3.300000000,
. 3.405000000, 3.492000000, 3.589000000, 3.521000000,
. 3.463000000, 3.260000000, 3.020000000, 2.680000000,
. 2.211000000, 1.710000000, 1.139000000, 0.444000000,
.-0.223000000, -1.060000000, -1.740000000, -2.542000000,
.-7.281000000, -12.708000000, -18.593000000, -24.904000000,
.-31.227000000, -37.790000000, -44.322000000, -50.768000000,
.-57.516000000, -64.187000000, -80.842000000, -92.865000000,
.-131.639000000, -117.769000000, -109.903000000, -122.847000000,
.-101.490000000, -45.186000000, -8.716000000, 17.884000000/
*
n= 120
DO l=1,n
x(l)= xc(l)
y(l)= yc(l)
ENDDO
*.....Set up tridiagonal system.........................................
* b=diagonal, d=offdiagonal, c=right-hand side
*
d(1)= x(2)-x(1)
c(2)= (y(2)-y(1))/d(1)
DO k= 2,n-1
d(k)= x(k+1)-x(k)
b(k)= TWO*(d(k-1)+d(k))
c(k+1)= (y(k+1)-y(k))/d(k)
c(k)= c(k+1)-c(k)
END DO
*
*.....End conditions. third derivatives at x(1) and x(n) obtained
* from divided differences.......................................
*
b(1)= -d(1)
b(n)= -d(n-1)
c(1)= ZERO
c(n)= ZERO
IF (n.gt.3) THEN
c(1)= c(3)/(x(4)-x(2))-c(2)/(x(3)-x(1))
c(n)= c(n-1)/(x(n)-x(n-2))-c(n-2)/(x(n-1)-x(n-3))
c(1)= c(1)*d(1)*d(1)/(x(4)-x(1))
c(n)= -c(n)*d(n-1)*d(n-1)/(x(n)-x(n-3))
END IF
*
DO k=2,n ! forward elimination
t= d(k-1)/b(k-1)
b(k)= b(k)-t*d(k-1)
c(k)= c(k)-t*c(k-1)
END DO
*
c(n)= c(n)/b(n)
*
* back substitution ( makes c the sigma of text)
*
DO k=n-1,1,-1
c(k)= (c(k)-d(k)*c(k+1))/b(k)
END DO
*
*.....Compute polynomial coefficients...................................
*
b(n)= (y(n)-y(n-1))/d(n-1)+d(n-1)*(c(n-1)+c(n)+c(n))
DO k=1,n-1
b(k)= (y(k+1)-y(k))/d(k)-d(k)*(c(k+1)+c(k)+c(k))
d(k)= (c(k+1)-c(k))/d(k)
c(k)= THREE*c(k)
END DO
c(n)= THREE*c(n)
d(n)= d(n-1)
*
RETURN
*
END
*
*------------------------------------------------------------------------
*
SUBROUTINE S4eval3Single2(u,b,c,d,gdim,f,fp,fpp,fppp)
*
* ---------------------------------------------------------------------------
*
* PURPOSE - Evaluate the cubic spline function
* Seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3
* where x(i) <= u < x(i+1)
*
* NOTES- if u<x(1), i=1 is used;if u>x(n), i=n is used
*
REAL*8 u
* abscissa at which the spline is to be evaluated
*
INTEGER j,k,n,l,top,gdim
*
REAL*8 xc(120),yc(120)
REAL*8 x(120),y(120)
REAL*8 b(gdim),c(gdim),d(gdim)
* linear,quadratic,cubic coeff
*
REAL*8 f,fp,fpp,fppp
* function, 1st,2nd,3rd deriv
*
INTEGER i
REAL*8 dx
REAL*8 ZERO, TWO, THREE
PARAMETER(ZERO=0.0, TWO=2.0, THREE=3.0)
SAVE i
DATA i/1/
*
* The grid
*
DATA (xc(l),l=1,120)/
.100.000000000, 110.000000000, 120.000000000, 130.000000000,
.140.000000000, 145.000000000, 150.000000000, 151.000000000,
.152.000000000, 153.000000000, 154.000000000, 155.000000000,
.156.000000000, 157.000000000, 158.000000000, 159.000000000,
.160.000000000, 161.000000000, 162.000000000, 163.000000000,
.164.000000000, 165.000000000, 166.000000000, 167.000000000,
.168.000000000, 169.000000000, 170.000000000, 171.000000000,
.172.000000000, 173.000000000, 174.000000000, 175.000000000,
.176.000000000, 177.000000000, 178.000000000, 179.000000000,
.180.000000000, 181.000000000, 182.000000000, 183.000000000,
.184.000000000, 185.000000000, 186.000000000, 187.000000000,
.188.000000000, 189.000000000, 190.000000000, 195.000000000,
.200.000000000, 210.000000000, 220.000000000, 230.000000000,
.240.000000000, 250.000000000, 260.000000000, 270.000000000,
.280.000000000, 290.000000000, 300.000000000, 310.000000000,
.320.000000000, 325.000000000, 330.000000000, 335.000000000,
.336.000000000, 337.000000000, 338.000000000, 339.000000000,
.340.000000000, 341.000000000, 342.000000000, 343.000000000,
.344.000000000, 345.000000000, 346.000000000, 347.000000000,
.348.000000000, 349.000000000, 350.000000000, 351.000000000,
.353.000000000, 355.000000000, 360.000000000, 365.000000000,
.370.000000000, 375.000000000, 380.000000000, 385.000000000,
.390.000000000, 400.000000000, 410.000000000, 420.000000000,
.430.000000000, 440.000000000, 450.000000000, 460.000000000,
.470.000000000, 480.000000000, 490.000000000, 500.000000000,
.550.000000000, 600.000000000, 650.000000000, 700.000000000,
.750.000000000, 800.000000000, 850.000000000, 900.000000000,
.950.000000000, 1000.000000000, 1100.000000000, 1150.000000000,
.1200.000000000, 1250.000000000, 1300.000000000, 1350.000000000,
.1400.000000000, 1600.000000000, 1800.000000000, 2000.000000000/
*
DATA (yc(l),l=1,120)/
.12.283000000, 12.212000000, 12.118000000, 11.982000000,
.11.780000000, 11.620000000, 11.394000000, 11.333000000,
.11.265000000, 11.191000000, 11.099000000, 10.994000000,
.10.862000000, 10.705000000, 10.510000000, 10.290000000,
.10.112000000, 10.104000000, 10.202000000, 10.238000000,
.10.197000000, 10.115000000, 10.019000000, 9.919000000,
. 9.822000000, 9.726000000, 9.632000000, 9.531000000,
. 9.432000000, 9.334000000, 9.234000000, 9.134000000,
. 9.027000000, 8.907000000, 8.786000000, 8.639000000,
. 8.485000000, 8.351000000, 8.278000000, 8.286000000,
. 8.320000000, 8.336000000, 8.300000000, 8.246000000,
. 8.176000000, 8.107000000, 8.057000000, 7.757000000,
. 7.465000000, 7.044000000, 6.720000000, 6.399000000,
. 6.065000000, 5.815000000, 5.600000000, 5.273000000,
. 5.038000000, 4.795000000, 4.541000000, 4.349000000,
. 4.106000000, 3.978000000, 3.847000000, 3.688000000,
. 3.663000000, 3.589000000, 3.539000000, 3.463000000,
. 3.417000000, 3.317000000, 3.223000000, 3.169000000,
. 3.048000000, 2.449000000, 2.361000000, 2.327000000,
. 2.356000000, 2.470000000, 2.570000000, 2.629000000,
. 2.721000000, 2.834000000, 3.049000000, 3.300000000,
. 3.405000000, 3.492000000, 3.589000000, 3.521000000,
. 3.463000000, 3.260000000, 3.020000000, 2.680000000,
. 2.211000000, 1.710000000, 1.139000000, 0.444000000,
.-0.223000000, -1.060000000, -1.740000000, -2.542000000,
.-7.281000000, -12.708000000, -18.593000000, -24.904000000,
.-31.227000000, -37.790000000, -44.322000000, -50.768000000,
.-57.516000000, -64.187000000, -80.842000000, -92.865000000,
.-131.639000000, -117.769000000, -109.903000000, -122.847000000,
.-101.490000000, -45.186000000, -8.716000000, 17.884000000/
*
n= 120
DO l=1,n
x(l)= xc(l)
y(l)= yc(l)
ENDDO
*
*.....First check if u is in the same interval found on the
* last call to Seval.............................................
*
IF ( (i.lt.1) .OR. (i.ge.n) ) i=1
IF ( (u.lt.x(i)) .OR. (u.ge.x(i+1)) ) THEN
i=1
*
* binary search
*
j= n+1
2468 k= (i+j)/2
IF (u.lt.x(k)) THEN
j= k
ELSE
i= k
ENDIF
IF (j.le.i+1) GOTO 2469
GOTO 2468
2469 CONTINUE
ENDIF
*
dx= u-x(i)
*
* evaluate the spline
*
f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i)))
fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i))
fpp= TWO*c(i) + dx*SIX*d(i)
fppp= SIX*d(i)
*
RETURN
*
END
Reference in New Issue View Git Blame Copy Permalink