subroutine spl1d1(n,x,f,w,iop,ij,a,b,c) c where n= number of points in the interpolation c x= origin of table of independent variable c f= origin of table of dependent variable c w= an array of dimension n which contains the calculated c second derivatives upon return c iop= an array of dimension 2 which contains combinations of c the integers 1 thru 5 used to specify the boundary c conditions c ij= spacing in the f and w tables c a,b,c= arrays of dimension n used for temporary storage c c-------------------------------------------------------------------------- c include 'ucom.h' c dimension iop(2),x(n),f(n*ij),w(n*ij),a(n),b(n),c(n) c-------------------------------------------------------------------------- c* k=n-1 a(2)=-(x(2)-x(1))/6. b(2)=(x(3)-x(1))/3. w(ij+1)=(f(2*ij+1)-f(ij+1))/(x(3)-x(2))-(f(ij+1)-f(1)) 1/(x(2)-x(1)) if (n-3)3,4,3 3 do 10 i=3,k m=(i-1)*ij+1 j1=m+ij j2=m-ij con=(x(i+1)-x(i-1))/3. don=(x(i)-x(i-1))/6. b(i)=con-(don**2)/b(i-1) e=(f(j1)-f(m))/(x(i+1)-x(i))-(f(m)-f(j2))/(x(i)-x(i-1)) w(m)=e-(don*w(j2))/b(i-1) 10 a(i)=-(don*a(i-1))/b(i-1) 4 k1=(n-2)*ij+1 c(n-1)=-((x(n)-x(n-1))/6.)/b(n-1) w(k1)=w(k1)/b(n-1) a(n-1)=a(n-1)/b(n-1) k2=k-1 if (n-3)7,8,7 7 do 20 i=2,k2 j=n-i con=(x(j+1)-x(j))/6. a(j)=(a(j)-con*a(j+1))/b(j) c(j)=-(con*c(j+1))/b(j) k3=(j-1)*ij+1 m=k3+ij 20 w(k3)=(w(k3)-con*w(m))/b(j) 8 k4=(n-1)*ij+1 if (iop(1)-5) 201,200,201 201 c1=w(1) if (iop(2)-5) 203,202,203 203 c2=w(k4) go to 205 200 if (n-4)300,302,302 302 a1=x(1)-x(2) a2=x(1)-x(3) a3=x(1)-x(4) a4=x(2)-x(3) a5=x(2)-x(4) a6=x(3)-x(4) c w(1)=f(1)*(1./a1+1./a2+1./a3)-a2*a3*f(ij+1)/(a1*a4*a5)+ c 1 a1*a3*f(2*ij+1)/(a2*a4*a6 )-a1*a2*f(3*ij+1)/(a3*a5*a6) w(1)=f(1)*(1./a1+1./a2+1./a3)-a2*a3*f(ij+1)/(a1*a4*a5) w(1)=w(1)+a1*a3*f(2*ij+1)/(a2*a4*a6 )-a1*a2*f(3*ij+1)/(a3*a5*a6) go to 201 202 if (n-4)300,303,303 303 b1=x(n)-x(n-3) b2=x(n)-x(n-2) b3=x(n)-x(n-1) b4=x(n-1)-x(n-3) b5=x(n-1)-x(n-2) b6=x(n-2)-x(n-3) l1=k4-ij l2=l1-ij l3=l2-ij c w(k4)=-b2*b3*f(l3)/(b6*b4*b1)+b1*b3*f(l2)/(b6*b5*b2) c 1 -b1*b2*f(l1)/(b4*b5*b3)+f(k4)*(1./b1+1./b2+1./b3) w(k4)=-b2*b3*f(l3)/(b6*b4*b1)+b1*b3*f(l2)/(b6*b5*b2) w(k4)=w(k4)-b1*b2*f(l1)/(b4*b5*b3)+f(k4)*(1./b1+1./b2+1./b3) go to 203 205 do 50 i=1,k m=(i-1)*ij+1 60 mk=iop(1) go to (62,64,66,68,66),mk 62 if (i-1)71,63,71 63 a(1)=-1. c(1)=0. go to500 71 bob=0. go to 500 64 if (i-1)73,76,73 76 a(1)=-1. c(1)=0. w(1)=0. go to 500 73 if (i-2)81,81,82 81 bob=-c1 go to 500 82 bob=0. go to 500 66 if (i-1)83,84,83 84 a(1)=-(x(2)-x(1))/3. c(1)=0. w(1)=-c1+(f(ij+1)-f(1))/(x(2)-x(1)) go to 500 83 if (i-2)85,85,86 85 bob=(x(2)-x(1))/6. go to 500 86 bob=0. go to 500 68 if (i-1)87,88,87 88 a(1)=-1. c(1)=1. w(1)=0. go to 500 87 bob=0. 500 ml=iop(2) go to (120,130,140,150,140),ml 120 if (i-1)121,122,121 122 a(n)=0. c(n)=-1. go to 70 121 bill=0. go to 70 130 if (i-1)131,132,131 132 a(n)=0. c(n)=-1. w(k4)=0. go to 70 131 if (i-k)134,133,134 133 bill=-c2 go to 70 134 bill=0. go to 70 140 if (i-1)141,142,141 142 a(n)=0. c(n)=(x(n-1)-x(n))/3. w(k4)=c2-(f(k4)-f(k1))/(x(n)-x(n-1)) go to 70 141 if (i-k)143,144,143 144 bill=(x(n)-x(n-1))/6. go to 70 143 bill=0. go to 70 150 if (i-1)151,152,151 152 a(n)=0. c(n)=(x(n-1)+x(1)-x(n)-x(2))/3. w(k4)=(f(ij+1)-f(1))/(x(2)-x(1))-(f(k4)-f(k1))/(x(n)-x(n-1)) go to 70 151 if (i-2)153,154,153 154 bill=(x(2)-x(1))/6. go to 70 153 if (i-k)155,156,155 156 bill=(x(n)-x(n-1))/6. go to 70 155 bill=0. 70 if (i-1)80,50,80 80 w(1)=w(1)-bob*w(m) w(k4)=w(k4)-bill*w(m) a(1)=a(1)-bob*a(i) a(n)=a(n)-bill*a(i) c(1)=c(1)-bob*c(i) c(n)=c(n)-bill*c(i) 50 continue go to 100 100 con=a(1)*c(n)-c(1)*a(n) d1=-w(1) d2=-w(k4) w(1)=(d1*c(n)-c(1)*d2)/con w(k4)=(a(1)*d2-d1*a(n))/con do 110 i=2,k m=(i-1)*ij+1 110 w(m)=w(m)+a(i)*w(1)+c(i)*w(k4) go to 305 300 write(ndiag,*)' hsp1d1, n.lt.4 results incorrect ' call appendparm stop ' Abnormal stop spl1d1 ' 305 return end c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++