******************************************************************************* * * * 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 * * * ******************************************************************************* * PROGRAM exa_ewgrid * IMPLICIT NONE * INTEGER i,top,gdim REAL*8 u,ui,value REAL*8, dimension(154) :: bl,cl,dl REAL*8, dimension(151) :: bc,cc,dc REAL*8, dimension(152) :: bu,cu,du * * u value of M_H at which the spline is to be evaluated * top= -1,0,1 lower, central, upper value for m_top * top= 1 * IF(top.eq.-1) THEN gdim= 154 ELSEIF(top.eq.0) THEN gdim= 151 ELSEIF(top.eq.1) THEN gdim= 152 ENDIF * IF(top.eq.-1) THEN CALL FMMsplineSingle(bl,cl,dl,top,gdim) ELSEIF(top.eq.0) THEN CALL FMMsplineSingle(bc,cc,dc,top,gdim) ELSEIF(top.eq.1) THEN CALL FMMsplineSingle(bu,cu,du,top,gdim) ENDIF * ui= 307.5d0 DO i=1,11 u= ui+(i-1)/10.d0 IF(top.eq.-1) THEN CALL Seval3Single(u,bl,cl,dl,top,gdim,value) ELSEIF(top.eq.0) THEN CALL Seval3Single(u,bc,cc,dc,top,gdim,value) ELSEIF(top.eq.1) THEN CALL Seval3Single(u,bu,cu,du,top,gdim,value) ENDIF print 1,u,value ENDDO * 1 format(1x,2e20.5) * STOP * *----------------------------------------------------------------------- * 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, dimension(154) :: xl,yl REAL*8, dimension(151) :: xc,yc REAL*8, dimension(152) :: xu,yu REAL*8, dimension(154) :: x,y * REAL*8, DIMENSION(gdim) :: b * linear coeff * REAL*8, DIMENSION(gdim) :: c * quadratic coeff. * REAL*8, DIMENSION(gdim) :: d * cubic coeff. * REAL*8 :: t REAL*8,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 x= xl y= yl ELSEIF(top.eq.0) THEN n= 151 FORALL(l=1:151) x(l)= xc(l) y(l)= yc(l) ENDFORALL ELSEIF(top.eq.1) THEN n= 152 FORALL(l=1:152) x(l)= xu(l) y(l)= yu(l) ENDFORALL 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 > 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 FMMsplineSingle * *------------------------------------------------------------------------ * 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 ux(n), i=n is used * REAL*8,INTENT(IN) :: u * abscissa at which the spline is to be evaluated * INTEGER j,k,n,l,top,gdim * REAL*8, dimension(154) :: xl,yl REAL*8, dimension(151) :: xc,yc REAL*8, dimension(152) :: xu,yu REAL*8, dimension(154) :: x,y REAL*8, DIMENSION(gdim) :: b,c,d * linear,quadratic,cubic coeff * REAL*8,INTENT(OUT),OPTIONAL:: f,fp,fpp,fppp * function, 1st,2nd,3rd deriv * INTEGER, SAVE :: i=1 REAL*8 :: dx REAL*8,PARAMETER:: TWO=2.0, THREE=3.0, SIX=6.0 * * 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 x= xl y= yl ELSEIF(top.eq.0) THEN n= 151 FORALL(l=1:151) x(l)= xc(l) y(l)= yc(l) ENDFORALL ELSEIF(top.eq.1) THEN n= 152 FORALL(l=1:151) x(l)= xu(l) y(l)= yu(l) ENDFORALL ENDIF * *.....First check if u is in the same interval found on the * last call to Seval............................................. * IF ( (i<1) .OR. (i >= n) ) i=1 IF ( (u < x(i)) .OR. (u >= x(i+1)) ) THEN i=1 * * binary search * j= n+1 DO k= (i+j)/2 IF (u < x(k)) THEN j= k ELSE i= k ENDIF IF (j <= i+1) EXIT ENDDO ENDIF * dx= u-x(i) * * evaluate the spline * IF (Present(f)) f= y(i)+dx*(b(i)+dx*(c(i)+dx*d(i))) IF (Present(fp)) fp= b(i)+dx*(TWO*c(i) + dx*THREE*d(i)) IF (Present(fpp)) fpp= TWO*c(i) + dx*SIX*d(i) IF (Present(fppp)) fppp= SIX*d(i) * RETURN * END Subroutine Seval3Single * END PROGRAM exa_ewgrid