1 SUBROUTINE v_spline(k_bc1,k_bcn,n,x,f,wk)
110 INTEGER k_bc1, k_bcn,
123 integer :: i_bc1,i_bcn
125 integer :: iord1,iord2
127 real :: h1,h2,h3,dels
128 real :: f1,f2,f3,aa,bb
139 if((i_bc1.lt.-1).or.(i_bc1.gt.7)) i_bc1=0
140 if((i_bcn.lt.0).or.(i_bcn.gt.7)) i_bcn=0
142 if(i_bc1.eq.-1) i_bcn=-1
149 if(i_bc1.eq.0) i3knots = i3knots + 1
150 if(i_bcn.eq.0) i3knots = i3knots + 1
151 if(i_bc1.eq.-1) i3perio = 1
159 if((i_bc1.eq.0).or.(i_bc1.gt.5)) i_bc1=5
160 if((i_bcn.eq.0).or.(i_bcn.gt.5)) i_bcn=5
162 if((i_bc1.eq.1).or.(i_bc1.eq.3).or.(i_bc1.eq.5))
then
168 if((i_bcn.eq.1).or.(i_bcn.eq.3).or.(i_bcn.eq.5))
then
185 ELSEIF(i_bc1.eq.2)
THEN
187 ELSEIF(i_bc1.eq.5)
THEN
188 a1=(f(1,2)-f(1,1))/(x(2)-x(1))
189 ELSEIF(i_bc1.eq.6)
THEN
190 b1=2.0*((f(1,3)-f(1,2))/(x(3)-x(2))
191 & -(f(1,2)-f(1,1))/(x(2)-x(1)))/(x(3)-x(1))
195 ELSEIF(i_bcn.eq.2)
THEN
197 ELSEIF(i_bcn.eq.5)
THEN
198 an=(f(1,n)-f(1,n-1))/(x(n)-x(n-1))
199 ELSEIF(i_bcn.eq.6)
THEN
200 bn=2.0*((f(1,n)-f(1,n-1))/(x(n)-x(n-1))
201 & -(f(1,n-1)-f(1,n-2))/(x(n-1)-x(n-2)))/(x(n)-x(n-2))
208 if((i_bc1.eq.5).and.(i_bcn.eq.5))
then
210 f(2,1)=(f(1,2)-f(1,1))/(x(2)-x(1))
216 else if((iord1.eq.1).and.(iord2.eq.1))
then
234 f(3,1) = (3*(fh-f0)/(h*h) - (2*a1 + an)/h)*2
235 f(4,1) = (-2*(fh-f0)/(h*h*h) + (a1 + an)/(h*h))*6
238 f(3,2) = f(4,1)*h + f(3,1)
240 else if((iord1.eq.1).and.(iord2.eq.2))
then
258 f(3,1) = (-bn/4 + 3*(fh-f0)/(2*h*h) - 3*a1/(2*h))*2
259 f(4,1) = (bn/(4*h) - (fh-f0)/(2*h*h*h) + a1/(2*h*h))*6
262 f(2,2) = f(4,1)*h*h/2 + f(3,1)*h + a1
263 else if((iord1.eq.2).and.(iord2.eq.1))
then
281 f(2,1) = 3*(fh-f0)/(2*h) - b1*h/4 - an/2
282 f(4,1) = (an/(2*h*h) - (fh-f0)/(2*h*h*h) - b1/(4*h))*6
285 f(3,2) = f(4,1)*h + f(3,1)
286 else if((iord1.eq.2).and.(iord2.eq.2))
then
304 f(2,1) = (fh-f0)/h -b1*h/3 -bn*h/6
308 f(2,2) = f(4,1)*h*h/2 + b1*h + f(2,1)
311 ELSE IF(i3perio.eq.1)
then
317 dels=(f(1,3)-f(1,2))/h2 - (f(1,2)-f(1,1))/h1
319 f(2,1)= (f(1,2)-f(1,1))/h1 + (h1*dels)/h
321 f(4,1)= 12*dels/(h1*h)
323 f(2,2)= (f(1,3)-f(1,2))/h2 - (h2*dels)/h
325 f(4,2)= -12*dels/(h2*h)
332 ELSE IF(i3knots.eq.2)
then
345 aa = (f2*h1 + f1*h2)/(h1*h2*h)
346 bb = (f2*h1*h1 - f1*h2*h2)/(h1*h2*h)
355 ELSE IF(i3knots.eq.1)
then
357 if((i_bc1.eq.1).or.(i_bc1.eq.3).or.(i_bc1.eq.5))
then
370 aa=a1/(h2*h3) + f3/(h3*h3*(h3-h2)) - f2/(h2*h2*(h3-h2))
371 bb=-a1*(h3*h3-h2*h2)/(h2*h3*(h3-h2))
372 > + f2*h3/(h2*h2*(h3-h2)) - f3*h2/(h3*h3*(h3-h2))
378 f(2,2)=3*aa*h2*h2 + 2*bb*h2 + a1
379 f(3,2)=6*aa*h2 + 2*bb
382 f(2,3)=3*aa*h3*h3 + 2*bb*h3 + a1
383 f(3,3)=6*aa*h3 + 2*bb
386 else if((i_bc1.eq.2).or.(i_bc1.eq.4).or.(i_bc1.eq.6))
then
399 aa= -(b1/2)*(h3-h2)/(h3*h3-h2*h2)
400 > -f2/(h2*(h3*h3-h2*h2)) + f3/(h3*(h3*h3-h2*h2))
401 bb= -(b1/2)*h2*h3*(h3-h2)/(h3*h3-h2*h2)
402 > +f2*h3*h3/(h2*(h3*h3-h2*h2))
403 > -f3*h2*h2/(h3*(h3*h3-h2*h2))
409 f(2,2)=3*aa*h2*h2 + b1*h2 + bb
413 f(2,3)=3*aa*h3*h3 + b1*h3 + bb
417 else if((i_bcn.eq.1).or.(i_bcn.eq.3).or.(i_bcn.eq.5))
then
430 aa=an/(h2*h3) + f3/(h3*h3*(h3-h2)) - f2/(h2*h2*(h3-h2))
431 bb=-an*(h3*h3-h2*h2)/(h2*h3*(h3-h2))
432 > + f2*h3/(h2*h2*(h3-h2)) - f3*h2/(h3*h3*(h3-h2))
438 f(2,2)=3*aa*h2*h2 + 2*bb*h2 + an
439 f(3,2)=6*aa*h2 + 2*bb
442 f(2,1)=3*aa*h3*h3 + 2*bb*h3 + an
443 f(3,1)=6*aa*h3 + 2*bb
446 else if((i_bcn.eq.2).or.(i_bcn.eq.4).or.(i_bcn.eq.6))
then
459 aa= -(bn/2)*(h3-h2)/(h3*h3-h2*h2)
460 > -f2/(h2*(h3*h3-h2*h2)) + f3/(h3*(h3*h3-h2*h2))
461 bb= -(bn/2)*h2*h3*(h3-h2)/(h3*h3-h2*h2)
462 > +f2*h3*h3/(h2*(h3*h3-h2*h2))
463 > -f3*h2*h2/(h3*(h3*h3-h2*h2))
469 f(2,2)=3*aa*h2*h2 + bn*h2 + bb
473 f(2,1)=3*aa*h3*h3 + bn*h3 + bb
485 f(3,2)=(f(1,2)-f(1,1))/f(4,1)
488 f(2,i)=2.0*(f(4,i-1)+f(4,i))
489 f(3,i+1)=(f(1,i+1)-f(1,i))/f(4,i)
490 f(3,i)=f(3,i+1)-f(3,i)
501 f(2,1)=2.0*(f(4,1)+f(4,n-1))
502 f(3,1)=(f(1,2)-f(1,1))/f(4,1)-(f(1,n)-f(1,n-1))/f(4,n-1)
509 ELSEIF(i_bc1.eq.1.or.i_bc1.eq.3.or.i_bc1.eq.5)
THEN
511 f(3,1)=(f(1,2)-f(1,1))/f(4,1)-a1
512 ELSEIF(i_bc1.eq.2.or.i_bc1.eq.4.or.i_bc1.eq.6)
THEN
516 ELSEIF(i_bc1.eq.7)
THEN
518 f(3,1)=f(3,3)/(x(4)-x(2))-f(3,2)/(x(3)-x(1))
519 f(3,1)=f(3,1)*f(4,1)**2/(x(4)-x(1))
522 f(2,2)=f(4,1)+2.0*f(4,2)
523 f(3,2)=f(3,2)*f(4,2)/(f(4,1)+f(4,2))
526 IF(i_bcn.eq.1.or.i_bcn.eq.3.or.i_bcn.eq.5)
THEN
528 f(3,n)=-(f(1,n)-f(1,n-1))/f(4,n-1)+an
529 ELSEIF(i_bcn.eq.2.or.i_bcn.eq.4.or.i_bcn.eq.6)
THEN
531 f(3,n)=f(4,n-1)*bn/3.0
534 ELSEIF(i_bcn.eq.7)
THEN
536 f(3,n)=f(3,n-1)/(x(n)-x(n-2))-f(3,n-2)/(x(n-1)-x(n-3))
537 f(3,n)=-f(3,n)*f(4,n-1)**2/(x(n)-x(n-3))
538 ELSEIF(i_bc1.ne.-1)
THEN
540 f(2,n-1)=2.0*f(4,n-2)+f(4,n-1)
541 f(3,n-1)=f(3,n-1)*f(4,n-2)/(f(4,n-1)+f(4,n-2))
549 f(2,i)=f(2,i)-t*f(4,i-1)
550 f(3,i)=f(3,i)-t*f(3,i-1)
551 wk(i)=wk(i)-t*wk(i-1)
554 f(2,n-1)=f(2,n-1)-q*wk(i-1)
555 f(3,n-1)=f(3,n-1)-q*f(3,i-1)
558 wk(n-1)=wk(n-1)+f(4,n-2)
562 f(2,n-1)=f(2,n-1)-t*wk(n-2)
563 f(3,n-1)=f(3,n-1)-t*f(3,n-2)
565 f(3,n-1)=f(3,n-1)/f(2,n-1)
566 f(3,n-2)=(f(3,n-2)-wk(n-2)*f(3,n-1))/f(2,n-2)
569 f(3,i)=(f(3,i)-f(4,i)*f(3,i+1)-wk(i)*f(3,n-1))/f(2,i)
577 IF((i.eq.n-1).and.(imax.eq.n-1))
THEN
578 t=(f(4,i-1)-f(4,i))/f(2,i-1)
588 IF((i.eq.imin+1).and.(imin.eq.2))
THEN
589 f(2,i)=f(2,i)-t*(f(4,i-1)-f(4,i-2))
591 f(2,i)=f(2,i)-t*f(4,i-1)
593 f(3,i)=f(3,i)-t*f(3,i-1)
596 f(3,imax)=f(3,imax)/f(2,imax)
599 IF((i.eq.2).and.(imin.eq.2))
THEN
600 f(3,i)=(f(3,i)-(f(4,i)-f(4,i-1))*f(3,i+1))/f(2,i)
602 f(3,i)=(f(3,i)-f(4,i)*f(3,i+1))/f(2,i)
609 IF(i_bc1.le.0.or.i_bc1.gt.7)
THEN
610 f(3,1)=(f(3,2)*(f(4,1)+f(4,2))-f(3,3)*f(4,1))/f(4,2)
613 IF(i_bcn.le.0.or.i_bcn.gt.7)
THEN
614 f(3,n)=f(3,n-1)+(f(3,n-1)-f(3,n-2))*f(4,n-1)/f(4,n-2)
621 > (f(1,i+1)-f(1,i))/f(4,i)-f(4,i)*(f(3,i+1)+2.0*f(3,i))
622 f(4,i)=(f(3,i+1)-f(3,i))/f(4,i)
632 f(2,n)=f(2,n-1)+hn*(f(3,n-1)+0.5*hn*f(4,n-1))
633 f(3,n)=f(3,n-1)+hn*f(4,n-1)
635 IF(i_bcn.eq.1.or.i_bcn.eq.3.or.i_bcn.eq.5)
THEN
637 ELSE IF(i_bcn.eq.2.or.i_bcn.eq.4.or.i_bcn.eq.6)
THEN