spline.f

Go to the documentation of this file.

      subroutine spline(s,x,y,n,in,t,il,iu,vl,vu,e,u)                   
c     x,y array of ind. & depen. var.,s the argument to be interpolated 
c     t the interpolated value,n the dimension of (x,y)                 
c     in=1 determines spline func,in=2 interpolates spline fnc          
c     il,iu=1 parabolic runout conditions at lower & upper bound.       
c     il,iu=2 1st derivative  (vl,vu) at lower or upper bound resp.     
c     il,iu=3 2nd derivative  (vl,vu) at lower or upper bound resp.     
      implicit real*8 (a-h,o-z)
      dimension x(*), y(*), e(*),u(*)                                   
      go to (10,40),in                                                  
   10 continue                                                          
      n1=n-1                                                            
      b1=x(2)-x(1)                                                      
      c1=(y(2)-y(1))/b1                                                 
      go to (12,14,16),il                                               
   12 e(1)=1.0                                                          
      u(1)=0.0                                                          
      go to 18                                                          
   14 e(1)=-.5                                                          
      u(1)=(c1-vl)/2/b1                                                 
      go to 18                                                          
   16 e(1)=0                                                            
      u(1)=vl/12                                                        
   18 continue                                                          
      do 20 j=2,n1                                                      
      b2=x(j+1)-x(j)                                                    
      c2=(y(j+1)-y(j))/b2                                               
      b=x(j+1)-x(j-1)                                                   
      d=(c2-c1)/b                                                       
      c=b1/b                                                            
      b1=b2                                                             
      c1=c2                                                             
      p=c*e(j-1)+2.0                                                    
      e(j)=(c-1.0)/p                                                    
   20 u(j)=(d-c*u(j-1))/p                                               
      go to (22,24,26),iu                                               
   22 continue                                                          
      e(n)=u(n1)/(1.0-e(n1))                                            
      go to 28                                                          
   24 c2=vu-(y(n)-y(n1))/(x(n)-x(n1))                                   
      e(n)=(c2/(x(n)-x(n1))-u(n1))/(2.0d0+e(n1))                        
      go to 28                                                          
   26 e(n)=vu/12.0                                                      
   28 continue                                                          
      do 30 kk=1,n1                                                     
      k=n-kk                                                            
      e(k)=e(k)*e(k+1)+u(k)                                             
c     to obtain the derivatives at the knots remove c from the          
c     following comments cards.then the array u will contain            
c     the derivatives of the spline function at the knots               
      b2=x(k+1)-x(k)                                                    
      u(k)=(y(k+1)-y(k))/b2-b2*(e(k)*2+e(k+1))                          
   30 continue                                                          
      b2=x(n)-x(n1)                                                     
      u(n)=(e(n1)+2.0d0*e(n))*b2+(y(n)-y(n1))/b2                        
   40 if (x(1).gt.x(n)) go to 50                                        
      idir=0                                                            
      mlb=0                                                             
      mub=n                                                             
      go to 60                                                          
   50 idir=1                                                            
      mlb=n                                                             
      mub=0                                                             
   60 if (s.ge.x(mub+idir)) go to 100                                   
      if (s.le.x(mlb+1-idir)) go to 110                                 
      ml=mlb                                                            
      mu=mub                                                            
      go to 80                                                          
   70 if (iabs(mu-ml).le.1) go to 120                                   
   80 mav=(ml+mu)/2                                                     
      if (s.lt.x(mav)) go to 90                                         
      ml=mav                                                            
      go to 70                                                          
   90 mu=mav                                                            
      go to 70                                                          
  100 mu=mub+2*idir                                                     
      go to 130                                                         
  110 mu=mlb+2*(1-idir)                                                 
      go to 130                                                         
  120 mu=mu+idir                                                        
130    continue                                                         
      t=(e(mu-1)*((x(mu)-s)**3)+e(mu)*((s-x(mu-1))**3)+                 
     1   (y(mu-1)-e(mu-1)*((x(mu)-x(mu-1))**2))*(x(mu)-s)+              
     2   (y(mu)-e(mu)*((x(mu)-x(mu-1))**2))*(s-x(mu-1)))/               
     3   (x(mu)-x(mu-1))                                                
       return                                                           
       end                                                              
c***********************************************************************