diff options
Diffstat (limited to 'packages/numlib/src/spl.pas')
| -rw-r--r-- | packages/numlib/src/spl.pas | 1114 |
1 files changed, 1114 insertions, 0 deletions
diff --git a/packages/numlib/src/spl.pas b/packages/numlib/src/spl.pas new file mode 100644 index 0000000000..68c76c399b --- /dev/null +++ b/packages/numlib/src/spl.pas @@ -0,0 +1,1114 @@ +{ + This file is part of the Numlib package. + Copyright (c) 1986-2000 by + Kees van Ginneken, Wil Kortsmit and Loek van Reij of the + Computational centre of the Eindhoven University of Technology + + FPC port Code by Marco van de Voort (marco@freepascal.org) + documentation by Michael van Canneyt (Michael@freepascal.org) + + Undocumented unit. B- and other Splines. Not imported by the other units + afaik. + + See the file COPYING.FPC, included in this distribution, + for details about the copyright. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + + **********************************************************************} +unit spl; +{$I direct.inc} + +interface + +uses typ, sle; + +function spl1bspv(q: ArbInt; var kmin1, c1: ArbFloat; x: ArbFloat; var term: ArbInt): ArbFloat; +function spl2bspv(qx, qy: ArbInt; var kxmin1, kymin1, c11: ArbFloat; x, y: ArbFloat; var term: ArbInt): ArbFloat; +procedure spl1bspf(M, Q: ArbInt; var XYW1: ArbFloat; + var Kmin1, C1, residu: ArbFloat; + var term: ArbInt); +procedure spl2bspf(M, Qx, Qy: ArbInt; var XYZW1: ArbFloat; + var Kxmin1, Kymin1, C11, residu: ArbFloat; + var term: ArbInt); +procedure spl1nati(n: ArbInt; var xyc1: ArbFloat; var term: ArbInt); +procedure spl1naki(n: ArbInt; var xyc1: ArbFloat; var term: ArbInt); +procedure spl1cmpi(n: ArbInt; var xyc1: ArbFloat; dy1, dyn: ArbFloat; + var term: ArbInt); +procedure spl1peri(n: ArbInt; var xyc1: ArbFloat; var term: ArbInt); +function spl1pprv(n: ArbInt; var xyc1: ArbFloat; t: ArbFloat; var term: ArbInt): ArbFloat; + +procedure spl1nalf(n: ArbInt; var xyw1: ArbFloat; lambda:ArbFloat; + var xac1, residu: ArbFloat; var term: ArbInt); +function spl2natv(n: ArbInt; var xyg0: ArbFloat; u, v: ArbFloat): ArbFloat; +procedure spl2nalf(n: ArbInt; var xyzw1: ArbFloat; lambda:ArbFloat; + var xyg0, residu: ArbFloat; var term: ArbInt); + { term = 1: succes, + term = 2: set linear equations is not "PD" + term = 4: Approx. number of points? On a line. + term = 3: wrong input n<3 or a weight turned out to be <=0 } +implementation +{$goto on} + +type + Krec = record K1, K2, K3, K4, K5, K6 : ArbFloat end; + +function spl1bspv(q: ArbInt; var kmin1, c1: ArbFloat; x: ArbFloat; var term: ArbInt): ArbFloat; +var c : arfloat1 absolute c1; + k : arfloat_1 absolute kmin1; + D1, D2, D3, + E2, E3, E4, E5: ArbFloat; + pk : ^Krec; + l, r, m : ArbInt; +begin + spl1bspv := NaN; + term := 3; { q >=4 ! } + if q<4 then exit; { at least 1 interval } + if (x<k[2]) or (x>k[q-1]) then exit; { x inside the interval } + term := 1; { Let's hope the params are good :-)} + l := 2; r := q-1; + while l+1<r do { after this loop goes: } + begin { k[l]<=x<=k[l+1] with } + m := (l+r) div 2; { k[l] < k[l+1] } + if x>=k[m] then l := m else r := m + end; + pk := @k[l-2]; { the (de) Boor algoritm .. } + with pk^ do + begin + E2 := X - K2; E3 := X - K3; E4 := K4 - X; E5 := K5 - X; + D2 := C[l]; D3 := C[l+1]; + D1 := ((X-K1)*D2+E4*C[l-1])/(K4-K1); + D2 := (E2*D3+E5*D2)/(K5-K2); + D3 := (E3*C[l+2]+(K6-X)*D3)/(K6-K3); + D1 := (E2*D2+E4*D1)/(K4-K2); + D2 := (E3*D3+E5*D2)/(K5-K3); + spl1bspv := (E3*D2+E4*D1)/(K4-K3) + end; +end; + +function spl2bspv(qx, qy: ArbInt; var kxmin1, kymin1, c11: ArbFloat; x, y: ArbFloat; var term: ArbInt): ArbFloat; +var pd: ^arfloat1; + i, iy: ArbInt; + c: arfloat1 absolute c11; +begin + GetMem(pd, qx*SizeOf(ArbFloat)); + i := 0; + iy := 1; + repeat + i := i + 1; + pd^[i] := spl1bspv(qy, kymin1, c[iy], y, term); + Inc(iy, qy) + until (i=qx) or (term<>1); + if term=1 + then spl2bspv := spl1bspv(qx, kxmin1, pd^[1], x, term) + else spl2bspv := NaN; + FreeMem(pd, qx*SizeOf(ArbFloat)); +end; + +(* Bron: NAG LIBRARY SUBROUTINE E02BAF *) + +function Imin(x, y: ArbInt): ArbInt; +begin if x<y then Imin := x else Imin := y end; + +type ar4 = array[1..$ffe0 div (4*SizeOf(ArbFloat)),1..4] of ArbFloat; + ar3 = array[1..$ffe0 div (3*SizeOf(ArbFloat)),1..3] of ArbFloat; + r_3 = record x, y, w: ArbFloat end; + r3Ar= array[1..$ffe0 div SizeOf(r_3)] of r_3; + r_4 = record x, y, z, w: ArbFloat end; + r4Ar= array[1..$ffe0 div SizeOf(r_4)] of r_4; + r4 = array[1..4] of ArbFloat; + r2 = array[1..2] of ArbFloat; + + r4x = record xy: R2; alfa, d: ArbFloat end; + r4xAr= array[1..$ffe0 div SizeOf(r4x)] of r4x; + nsp2rec = array[0..$ff80 div (3*SizeOf(ArbFloat))] of + record xy: R2; gamma: ArbFloat end; + +procedure spl1bspf(M, Q: ArbInt; var XYW1: ArbFloat; + var Kmin1, C1, residu: ArbFloat; + var term: ArbInt); +var work1: ^arfloat1; + work2: ^ar4; + c : arfloat1 absolute c1; + k : arfloat_1 absolute kmin1; + xyw : r3Ar absolute XYW1; + r, j, jmax, l, lplus1, i, iplusj, jold, jrev, + jplusl, iu, lplusu : ArbInt; + s, k0, k4, sigma, + d, d4, d5, d6, d7, d8, d9, + e2, e3, e4, e5, + n1, n2, n3, + relemt, dprime, cosine, sine, + acol, arow, crow, ccol, ss : ArbFloat; + pk : ^Krec; + +label einde; +(* + DOUBLE PRECISION C(NCAP7), K(NCAP7), W(M), WORK1(M), + * WORK2(4,NCAP7), X(M), Y(M) + .. Local Scalars .. + DOUBLE PRECISION ACOL, AROW, CCOL, COSINE, CROW, D, D4, D5, D6, + * D7, D8, D9, DPRIME, E2, E3, E4, E5, K0, K1, K2, + * K3, K4, K5, K6, N1, N2, N3, RELEMT, S, SIGMA, + * SINE, WI, XI + INTEGER I, IERROR, IPLUSJ, IU, J, JOLD, JPLUSL, JREV, L, + * L4, LPLUS1, LPLUSU, NCAP, NCAP3, NCAPM1, R +*) +begin + term := 3; + if q<4 then exit; + if m<q then exit; +(* + CHECK THAT THE VALUES OF M AND NCAP7 ARE REASONABLE + IF (NCAP7.LT.8 .OR. M.LT.NCAP7-4) GO TO 420 + NCAP = NCAP7 - 7 + NCAPM1 = NCAP - 1 + NCAP3 = NCAP + 3 + + IN ORDER TO DEFINE THE FULL B-SPLINE BASIS, AUGMENT THE + PRESCRIBED INTERIOR KNOTS BY KNOTS OF MULTIPLICITY FOUR + AT EACH END OF THE DATA RANGE. + +*) + for j:=-1 to 2 do k[j] := xyw[1].x; + for j:=q-1 to q+2 do k[j] := xyw[m].x; + + if (k[3]<=xyw[1].x) or (k[q-2]>=xyw[m].x) then exit; +(* + CHECK THAT THE KNOTS ARE ORDERED AND ARE INTERIOR + TO THE DATA INTERVAL. +*) + j := 3; while (k[j]<=k[j+1]) and (j<q-2) do Inc(j); + if j<q-2 then exit; +(* + CHECK THAT THE WEIGHTS ARE STRICTLY POSITIVE. +*) + j := 1; + while (xyw[j].w>0) and (j<m) do Inc(j); + if xyw[j].w<=0 then exit; +(* + CHECK THAT THE DATA ABSCISSAE ARE ORDERED, THEN FORM THE + ARRAY WORK1 FROM THE ARRAY X. THE ARRAY WORK1 CONTAINS + THE + SET OF DISTINCT DATA ABSCISSAE. +*) + GetMem(Work1, m*SizeOf(ArbFloat)); + GetMem(Work2, q*4*SizeOf(ArbFloat)); + r := 1; work1^[1] := xyw[1].x; + j := 1; + while (j<m) do + begin + Inc(j); + if xyw[j].x>work1^[r] + then begin Inc(r); work1^[r] := xyw[j].x end + else if xyw[j].x<work1^[r] then goto einde; + end; + if r<q then goto einde; + +(* + CHECK THE FIRST S AND THE LAST S SCHOENBERG-WHITNEY + CONDITIONS ( S = MIN(NCAP - 1, 4) ). +*) + jmax := Imin(q-4,4); + j := 1; + while (j<=jmax) do + begin + if (work1^[j]>=k[j+2]) or (k[q-j-1]>=work1^[r-j+1]) then goto einde; + Inc(j) + end; +(* + CHECK ALL THE REMAINING SCHOENBERG-WHITNEY CONDITIONS. +*) + Dec(r, 4); i := 4; j := 5; + while j<=q-4 do + begin + K0 := K[j+2]; K4 := K[J-2]; + repeat Inc(i) until (Work1^[i]>k4); + if (I>R) or (WORK1^[I]>=K0) then goto einde; + Inc(j) + end; + +(* + INITIALISE A BAND TRIANGULAR SYSTEM (I.E. A + MATRIX AND A RIGHT HAND SIDE) TO ZERO. THE + PROCESSING OF EACH DATA POINT IN TURN RESULTS + IN AN UPDATING OF THIS SYSTEM. THE SUBSEQUENT + SOLUTION OF THE RESULTING BAND TRIANGULAR SYSTEM + YIELDS THE COEFFICIENTS OF THE B-SPLINES. +*) + FillChar(Work2^, q*4*SizeOf(ArbFloat), 0); + FillChar(c, q*SizeOf(ArbFloat), 0); + + SIGMA := 0; j := 0; jold := 0; + for i:=1 to m do + with xyw[i] do + begin +(* + FOR THE DATA POINT (X(I), Y(I)) DETERMINE AN INTERVAL + K(J + 3) .LE. X .LT. K(J + 4) CONTAINING X(I). (IN THE + CASE J + 4 .EQ. NCAP THE SECOND EQUALITY IS RELAXED TO + INCLUDE + EQUALITY). +*) + while (x>=k[j+2]) and (j<=q-4) do Inc(j); + if j<>jold then + begin + pk := @k[j-1]; + with pk^ do + begin + D4 := 1/(K4-K1); D5 := 1/(K5-K2); D6 := 1/(K6-K3); + D7 := 1/(K4-K2); D8 := 1/(K5-K3); D9 := 1/(K4-K3) + end; + JOLD := J; + end; +(* + COMPUTE AND STORE IN WORK1(L) (L = 1, 2, 3, 4) THE VALUES + OF + THE FOUR NORMALIZED CUBIC B-SPLINES WHICH ARE NON-ZERO AT + X=X(I). +*) with pk^ do + begin + E5 := k5 - X; + E4 := K4 - X; + E3 := X - K3; + E2 := X - K2; + N1 := W*D9; + N2 := E3*N1*D8; + N1 := E4*N1*D7; + N3 := E3*N2*D6; + N2 := (E2*N1+E5*N2)*D5; + N1 := E4*N1*D4; + WORK1^[4] := E3*N3; + WORK1^[3] := E2*N2 + (K6-X)*N3; + WORK1^[2] := (X-K1)*N1 + E5*N2; + WORK1^[1] := E4*N1; + CROW := Y*W; + end; +(* + ROTATE THIS ROW INTO THE BAND TRIANGULAR SYSTEM USING PLANE + ROTATIONS. +*) + for lplus1:=1 to 4 do + begin L := LPLUS1 - 1; + RELEMT := WORK1^[LPLUS1]; + if relemt<>0 then + begin JPLUSL := J + L; + D := WORK2^[JPLUSL,1]; + IF (ABS(RELEMT)>=D) + then DPRIME := ABS(RELEMT)*SQRT(1+sqr(D/RELEMT)) + else DPRIME := D*SQRT(1+sqr(RELEMT/D)); + WORK2^[JPLUSL,1] := DPRIME; + COSINE := D/DPRIME; SINE := RELEMT/DPRIME; + for iu :=2 to 4-l do + begin + LPLUSU := L + IU; + ACOL := WORK2^[JPLUSL,iu]; + AROW := WORK1^[LPLUSU]; + WORK2^[JPLUSL,iu] := COSINE*ACOL + SINE*AROW; + WORK1^[LPLUSU] := COSINE*AROW - SINE*ACOL + end; + + CCOL := C[JPLUSL]; + C[JPLUSL] := COSINE*CCOL + SINE*CROW; + CROW := COSINE*CROW - SINE*CCOL + end; + end; + SIGMA := SIGMA + sqr(CROW) + end; + + residu := SIGMA; +(* + SOLVE THE BAND TRIANGULAR SYSTEM FOR THE B-SPLINE + COEFFICIENTS. IF A DIAGONAL ELEMENT IS ZERO, AND HENCE + THE TRIANGULAR SYSTEM IS SINGULAR, THE IMPLICATION IS + THAT THE SCHOENBERG-WHITNEY CONDITIONS ARE ONLY JUST + SATISFIED. THUS IT IS APPROPRIATE TO EXIT IN THIS + CASE WITH THE SAME VALUE (IFAIL=5) OF THE ERROR + INDICATOR. +*) + term := 2; + L := -1; + for jrev:=1 to q do + begin + J := q - JREV + 1; D := WORK2^[J,1]; + if d=0 then goto einde; + IF l<3 then L := L + 1; + S := C[j]; + for i:=1 to l do + begin + IPLUSJ := I + J; + S := S - WORK2^[j,i+1]*C[IPLUSJ]; + end; + C[J] := S/D + end; + + term:=1; +einde: + FreeMem(Work2, q*4*SizeOf(ArbFloat)); + FreeMem(Work1, m*SizeOf(ArbFloat)) + +end; + +procedure spl2bspf(M, Qx, Qy: ArbInt; var XYZW1: ArbFloat; + var Kxmin1, Kymin1, C11, residu: ArbFloat; + var term: ArbInt); + +(* !!!!!!!! Test input !!!!!!!!!! *) + +(* +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c part 1: determination of the number of knots and their position. c +c **************************************************************** c +c given a set of knots we compute the least-squares spline sinf(x,y), c +c and the corresponding weighted sum of squared residuals fp=f(p=inf). c +c if iopt=-1 sinf(x,y) is the requested approximation. c +c if iopt=0 or iopt=1 we check whether we can accept the knots: c +c if fp <=s we will continue with the current set of knots. c +c if fp > s we will increase the number of knots and compute the c +c corresponding least-squares spline until finally fp<=s. c +c the initial choice of knots depends on the value of s and iopt. c +c if iopt=0 we first compute the least-squares polynomial of degree c +c 3 in x and 3 in y; nx=nminx=2*3+2 and ny=nminy=2*3+2. c +c fp0=f(0) denotes the corresponding weighted sum of squared c +c residuals c +c if iopt=1 we start with the knots found at the last call of the c +c routine, except for the case that s>=fp0; then we can compute c +c the least-squares polynomial directly. c +c eventually the independent variables x and y (and the corresponding c +c parameters) will be switched if this can reduce the bandwidth of the c +c system to be solved. c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc *) + +function Min(a, b:ArbInt): ArbInt; +begin if a<b then Min := a else Min := b end; + +procedure WisselR(var x, y: ArbFloat); +var h: ArbFloat; begin h := x; x := y; y := h end; + +procedure Wisseli(var x, y: ArbInt); +var h: ArbInt; begin h := x; x := y; y := h end; + +procedure fprota(var cos1, sin1, a, b: ArbFloat); +var store: ArbFloat; +begin + store := b; b := cos1*b+sin1*a; a := cos1*a-sin1*store +end; + +procedure fpgivs(var piv, ww, cos1, sin1: ArbFloat); +var store, dd: ArbFloat; +begin + store := abs(piv); + if store>=ww + then dd := store*sqrt(1+sqr(ww/piv)) + else dd := ww*sqrt(1+sqr(piv/ww)); + cos1 := ww/dd; sin1 := piv/dd; ww := dd +end; + +procedure fpback(var a11, z1: ArbFloat; n, k: ArbInt; var c1: ArbFloat); +(* + subroutine fpback calculates the solution of the system of + equations a*c = z with a a n x n upper triangular matrix + of bandwidth k. + ArbFloat a(.,k) +*) +var a: arfloat1 absolute a11; + z: arfloat1 absolute z1; + c: arfloat1 absolute c1; + i, l: ArbInt; + store : ArbFloat; +begin + for i:=n downto 1 do + begin + store := z[i]; + for l:=min(n+1-i,k)-1 downto 1 do store := store-c[i+l]*a[(i-1)*k+l+1]; + c[i] := store/a[(i-1)*k+1] + end; +end; + +procedure fpbspl(var kmin1: ArbFloat; x: ArbFloat; l: ArbInt; var h: r4); +(* + subroutine fpbspl evaluates the 4 non-zero b-splines of + degree 3 at t(l) <= x < t(l+1) using the stable recurrence + relation of de boor and cox. +*) +var k : arfloat_1 absolute kmin1; + f : ArbFloat; + hh: array[1..3] of ArbFloat; + i, j, li, lj : ArbInt; +begin + h[1] := 1; + for j:=1 to 3 do + begin + for i:=1 to j do hh[i] := h[i]; + h[1] := 0; + for i:=1 to j do + begin + li := l+i; lj := li-j; + f := hh[i]/(k[li]-k[lj]); + h[i] := h[i]+f*(k[li]-x); + h[i+1] := f*(x-k[lj]) + end; + end; +end; + +procedure fporde(m, qx, qy: ArbInt; var xyzw1, kxmin1, kymin1: ArbFloat; + var nummer1, index1: ArbInt); +var xi, yi : ArbFloat; + i, im, num, + k, l : ArbInt; + xyzw : r4Ar absolute xyzw1; + kx : arfloat_1 absolute kxmin1; + ky : arfloat_1 absolute kymin1; + nummer : arint1 absolute nummer1; + index : arint1 absolute index1; +begin + for i:=1 to (qx-3)*(qy-3) do index[i] := 0; + for im:=1 to m do + with xyzw[im] do + begin + l := 2; while (x>=kx[l+1]) and (l<qx-2) do Inc(l); + k := 2; while (y>=ky[k+1]) and (k<qy-2) do Inc(k); + num := (l-2)*(qy-3)+k-1; + nummer[im] := index[num]; index[num] := im + end; +end; + +label einde; + +var x0, x1, y0, y1, eps, cos1, sin1, dmax, sigma, + wi, zi, hxi, piv : ArbFloat; + i, j, l, l1, l2, lx, ly, nreg, ncof, jrot, + inpanel, i1, j1, + iband, num, irot : ArbInt; + xyzw : r4Ar absolute xyzw1; + kx, ky : ^arfloat_1; + a, f, h : ^arfloat1; + c : arfloat1 absolute c11; + nummer, index : ^arint1; + hx, hy : r4; + ichang, fullrank : boolean; +begin + + eps := 10*macheps; +(* find the position of the additional knots which are needed for the + b-spline representation of s(x,y) *) + iband := 1+min(3*qy+3,3*qx+3); + if qy>qx then + begin + ichang := true; + kx := @kymin1; ky := @kxmin1; + for i:=1 to m do with xyzw[i] do Wisselr(x, y); + WisselI(qx, qy) + end else + begin + ichang := false; + kx := @kxmin1; ky := @kymin1; + end; + with xyzw[1] do begin x0 := x; x1 := x; y0 := y; y1 := y end; + for i:=2 to m do with xyzw[i] do + begin if x<x0 then x0 := x; if x>x1 then x1 := x; + if y<y0 then y0 := y; if y>y1 then y1 := y + end; + for i:=-1 to 2 do kx^[i] := x0; + for i:=-1 to 2 do ky^[i] := y0; + for i:=qx-1 to qx+2 do kx^[i] := x1; + for i:=qy-1 to qy+2 do ky^[i] := y1; +(* arrange the data points according to the panel they belong to *) + nreg := (qx-3)*(qy-3); + ncof := qx*qy; + GetMem(nummer, m*SizeOf(ArbInt)); + GetMem(index, nreg*SizeOf(ArbInt)); + GetMem(h, iband*SizeOf(ArbFloat)); + GetMem(a, iband*ncof*SizeOf(ArbFloat)); + GetMem(f, ncof*SizeOf(ArbFloat)); + fporde(m, qx, qy, xyzw1, kx^[-1], ky^[-1], nummer^[1], index^[1]); + for i:=1 to ncof do f^[i] := 0; + for j:=1 to ncof*iband do a^[j] := 0; + residu := 0; +(* fetch the data points in the new order. main loop for the different panels *) + for num:=1 to nreg do + begin + lx := (num-1) div (qy-3); l1 := lx+2; + ly := (num-1) mod (qy-3); l2 := ly+2; + jrot := lx*qy+ly; + inpanel := index^[num]; + while inpanel<>0 do + with xyzw[inpanel] do + begin + wi := w; zi := z*wi; + fpbspl(kx^[-1], x, l1, hx); + fpbspl(ky^[-1], y, l2, hy); + for i:=1 to iband do h^[i] := 0; + i1 := 0; + for i:=1 to 4 do + begin + hxi := hx[i]; j1 := i1; + for j:=1 to 4 do begin Inc(j1); h^[j1] := hxi*hy[j]*wi end; + Inc(i1, qy) + end; + irot := jrot; + for i:=1 to iband do + begin + Inc(irot); piv := h^[i]; + if piv<>0 then + begin + fpgivs(piv, a^[(irot-1)*iband+1], cos1, sin1); + fprota(cos1, sin1, zi, f^[irot]); + for j:=i+1 to iband do + fprota(cos1, sin1, h^[j], a^[(irot-1)*iband+j-i+1]) + end; + end; + residu := residu+sqr(zi); + inpanel := nummer^[inpanel] + end; + end; + + dmax := 0; + i := 1; + while i<ncof*iband do + begin + if dmax<a^[i] then dmax:=a^[i]; Inc(i, iband) + end; + + sigma := eps*dmax; + i := 1; fullrank := true; + while fullrank and (i<ncof*iband) do + begin + fullrank := a^[i]>sigma; Inc(i, iband) + end; + + term := 2; if not fullrank then goto einde; + term := 1; + + fpback(a^[1], f^[1], ncof, iband, c11); + if ichang then + begin + l1 := 1; + for i:=1 to qx do + begin + l2 := i; + for j:=1 to qy do + begin + f^[l2] := c[l1]; Inc(l1); Inc(l2, qx) + end; + end; + for i:=1 to ncof do c[i] := f^[i] + end; + +einde: + if ichang then for i:=1 to m do with xyzw[i] do Wisselr(x, y); + FreeMem(f, ncof*SizeOf(ArbFloat)); + FreeMem(a, iband*ncof*SizeOf(ArbFloat)); + FreeMem(h, iband*SizeOf(ArbFloat)); + FreeMem(index, nreg*SizeOf(ArbInt)); + FreeMem(nummer, m*SizeOf(ArbInt)) +end; + + +procedure spl1nati(n: ArbInt; var xyc1: ArbFloat; var term: ArbInt); +var + xyc : r3Ar absolute XYC1; + l, b, d, u, c : ^arfloat1; + h2, h3, s2, s3: ArbFloat; + i, m : ArbInt; { afmeting van op te lossen stelsel } +begin + term:=3; + if n < 2 then exit; + for i:=2 to n do if xyc[i-1].x>=xyc[i].x then exit; + term:=1; + xyc[1].w := 0; xyc[n].w := 0; { c1=cn=0 } + m := n-2; + if m=0 then exit; + + getmem(u, n*SizeOf(ArbFloat)); + getmem(l, n*Sizeof(ArbFloat)); + getmem(d, n*SizeOf(ArbFloat)); + getmem(c, n*SizeOf(ArbFloat)); + getmem(b, n*SizeOf(ArbFloat)); + h3:=xyc[2].x-xyc[1].x; + s3:=(xyc[2].y-xyc[1].y)/h3; + + for i:=2 to n-1 do + begin + h2:=h3; h3:=xyc[i+1].x-xyc[i].x; + s2:=s3; s3:=(xyc[i+1].y-xyc[i].y)/h3; + l^[i]:=h2/6; + d^[i]:=(h2+h3)/3; + u^[i]:=h3/6; + b^[i]:=s3-s2 + end; + sledtr(m, l^[3], d^[2], u^[2], b^[2], c^[2], term); + for i:=2 to n-1 do xyc[i].w := c^[i]; + Freemem(b, n*SizeOf(ArbFloat)); + Freemem(c, n*SizeOf(ArbFloat)); + Freemem(d, n*SizeOf(ArbFloat)); + Freemem(l, n*Sizeof(ArbFloat)); + Freemem(u, n*SizeOf(ArbFloat)); +end; {spl1nati} + +procedure spl1naki(n: ArbInt; var xyc1: ArbFloat; var term: ArbInt); +var + xyc : r3Ar absolute XYC1; + l, b, d, u, c : ^arfloat1; + h2, h3, s2, s3: ArbFloat; + i, m : ArbInt; { Dimensions of set lin eqs to solve} +begin + term:=3; + if n < 4 then exit; + for i:=2 to n do if xyc[i-1].x>=xyc[i].x then exit; + term:=1; + m := n-2; + getmem(u, n*SizeOf(ArbFloat)); + getmem(l, n*Sizeof(ArbFloat)); + getmem(d, n*SizeOf(ArbFloat)); + getmem(c, n*SizeOf(ArbFloat)); + getmem(b, n*SizeOf(ArbFloat)); + h3:=xyc[2].x-xyc[1].x; + s3:=(xyc[2].y-xyc[1].y)/h3; + for i:=2 to n-1 do + begin + h2:=h3; h3:=xyc[i+1].x-xyc[i].x; + s2:=s3; s3:=(xyc[i+1].y-xyc[i].y)/h3; + l^[i]:=h2/6; + d^[i]:=(h2+h3)/3; + u^[i]:=h3/6; + b^[i]:=s3-s2 + end; + d^[n-1]:=d^[n-1]+h3/6*(1+h3/h2); l^[n-1]:=l^[n-1]-sqr(h3)/(6*h2); + h2:=xyc[2].x-xyc[1].x; h3:=xyc[3].x-xyc[2].x; + d^[2]:=d^[2]+h2/6*(1+h2/h3); u^[2]:=u^[2]-sqr(h2)/(6*h3); + + sledtr(m, l^[3], d^[2], u^[2], b^[2], c^[2], term); + for i:=2 to n-1 do xyc[i].w := c^[i]; + xyc[1].w := xyc[2].w + (h2/h3)*(xyc[2].w-xyc[3].w); + h2:=xyc[n-1].x-xyc[n-2].x; h3:=xyc[n].x-xyc[n-1].x; + xyc[n].w := xyc[n-1].w + (h3/h2)*(xyc[n-1].w-xyc[n-2].w); + Freemem(b, n*SizeOf(ArbFloat)); + Freemem(c, n*SizeOf(ArbFloat)); + Freemem(d, n*SizeOf(ArbFloat)); + Freemem(l, n*Sizeof(ArbFloat)); + Freemem(u, n*SizeOf(ArbFloat)); +end; {spl1naki} + +procedure spl1cmpi(n: ArbInt; var xyc1: ArbFloat; dy1, dyn: ArbFloat; + var term: ArbInt); +var + xyc : r3Ar absolute XYC1; + l, b, d, u, c : ^arfloat1; + h2, h3, s2, s3: ArbFloat; + i : ArbInt; { Dimensions of set lin eqs to solve} +begin + term:=3; + if n < 2 then exit; + for i:=2 to n do if xyc[i-1].x>=xyc[i].x then exit; + term:=1; + getmem(u, n*SizeOf(ArbFloat)); + getmem(l, n*Sizeof(ArbFloat)); + getmem(d, n*SizeOf(ArbFloat)); + getmem(c, n*SizeOf(ArbFloat)); + getmem(b, n*SizeOf(ArbFloat)); + h3:=xyc[2].x-xyc[1].x; + s3:=(xyc[2].y-xyc[1].y)/h3; + d^[1] := h3/3; u^[1] := h3/6; b^[1] := -dy1+s3; + for i:=2 to n-1 do + begin + h2:=h3; h3:=xyc[i+1].x-xyc[i].x; + s2:=s3; s3:=(xyc[i+1].y-xyc[i].y)/h3; + l^[i]:=h2/6; + d^[i]:=(h2+h3)/3; + u^[i]:=h3/6; + b^[i]:=s3-s2 + end; + d^[n] := h3/3; l^[n] := h3/6; b^[n] := dyn-s3; + + sledtr(n, l^[2], d^[1], u^[1], b^[1], c^[1], term); + for i:=1 to n do xyc[i].w := c^[i]; + Freemem(b, n*SizeOf(ArbFloat)); + Freemem(c, n*SizeOf(ArbFloat)); + Freemem(d, n*SizeOf(ArbFloat)); + Freemem(l, n*Sizeof(ArbFloat)); + Freemem(u, n*SizeOf(ArbFloat)); +end; {spl1cmpi} + +procedure spl1peri(n: ArbInt; var xyc1: ArbFloat; var term: ArbInt); +var + xyc : r3Ar absolute XYC1; + l, b, d, u, c, k : ^arfloat1; + k2, kn1, dy1, cn, + h2, h3, s2, s3: ArbFloat; + i, m : ArbInt; { Dimensions of set lin eqs to solve} +begin + term:=3; + if n < 2 then exit; + if xyc[1].y<>xyc[n].y then exit; + for i:=2 to n do if xyc[i-1].x>=xyc[i].x then exit; + term:=1; + m := n-2; + xyc[1].w := 0; xyc[n].w := 0; { c1=cn=0 } + if m=0 then exit; + if m=1 then + begin + h2:=xyc[2].x-xyc[1].x; + s2:=(xyc[2].y-xyc[1].y)/h2; + h3:=xyc[3].x-xyc[2].x; + s3:=(xyc[3].y-xyc[2].y)/h3; + xyc[2].w := 6*(s3-s2)/(h2+h3); + xyc[3].w := -xyc[2].w; + xyc[1].w := xyc[3].w; + exit + end; + + getmem(u, n*SizeOf(ArbFloat)); + getmem(l, n*Sizeof(ArbFloat)); + getmem(k, n*SizeOf(ArbFloat)); + getmem(d, n*SizeOf(ArbFloat)); + getmem(c, n*SizeOf(ArbFloat)); + getmem(b, n*SizeOf(ArbFloat)); + h3:=xyc[2].x-xyc[1].x; + s3:=(xyc[2].y-xyc[1].y)/h3; + k2 := h3/6; dy1 := s3; + for i:=2 to n-1 do + begin + h2:=h3; h3:=xyc[i+1].x-xyc[i].x; + s2:=s3; s3:=(xyc[i+1].y-xyc[i].y)/h3; + l^[i]:=h2/6; + d^[i]:=(h2+h3)/3; + u^[i]:=h3/6; + b^[i]:=s3-s2; + k^[i]:=0 + end; + kn1 := h3/6; k^[2] := k2; k^[n-1] := kn1; + sledtr(m, l^[3], d^[2], u^[2], k^[2], k^[2], term); + sledtr(m, l^[3], d^[2], u^[2], b^[2], c^[2], term); + cn := (dy1-s3-k2*c^[2]-kn1*c^[n-1])/(2*(k2+kn1)-k2*k^[2]-kn1*k^[n-1]); + for i:=2 to n-1 do xyc[i].w := c^[i] - cn*k^[i]; + xyc[1].w := cn; xyc[n].w := cn; + Freemem(b, n*SizeOf(ArbFloat)); + Freemem(c, n*SizeOf(ArbFloat)); + Freemem(d, n*SizeOf(ArbFloat)); + Freemem(l, n*Sizeof(ArbFloat)); + Freemem(k, n*SizeOf(ArbFloat)); + Freemem(u, n*SizeOf(ArbFloat)); +end; {spl1peri} + +function spl1pprv(n: ArbInt; var xyc1: ArbFloat; t: ArbFloat; var term: ArbInt): ArbFloat; +var + xyc : r3Ar absolute XYC1; + i, j, m : ArbInt; + d, d3, h, dy : ArbFloat; +begin { Assumption : x[i]<x[i+1] i=1..n-1 } + spl1pprv := NaN; + term:=3; if n<2 then exit; + if (t<xyc[1].x) or (t>xyc[n].x) then exit; + term:=1; + i:=1; j:=n; + while j <> i+1 do + begin + m:=(i+j) div 2; + if t>=xyc[m].x then i:=m else j:=m + end; { x[i]<= t <=x[i+1] } + h := xyc[i+1].x-xyc[i].x; + d := t-xyc[i].x; + d3 :=(xyc[i+1].w-xyc[i].w)/h; + dy :=(xyc[i+1].y-xyc[i].y)/h-h*(2*xyc[i].w+xyc[i+1].w)/6; + spl1pprv:= xyc[i].y+d*(dy+d*(xyc[i].w/2+d*d3/6)) + +end; {spl1pprv} + +procedure spl1nalf(n: ArbInt; var xyw1: ArbFloat; lambda:ArbFloat; + var xac1, residu: ArbFloat; var term: ArbInt); +var + xyw : r3Ar absolute xyw1; + xac : r3Ar absolute xac1; + i, j, ncd : ArbInt; + ca, crow : ArbFloat; + h, qty : ^arfloat1; + ch : ^arfloat0; + qtdq : ^arfloat1; +begin + term := 3; { testing input} + if n<2 then exit; + for i:=2 to n do if xyw[i-1].x>=xyw[i].x then exit; + for i:=1 to n do if xyw[i].w<=0 then exit; + if lambda<0 then exit; + term := 1; + Move(xyw, xac, n*SizeOf(r_3)); + if n=2 then begin xac[1].w := 0; xac[2].w := 0; exit end; + + Getmem(ch, (n+2)*SizeOf(ArbFloat)); FillChar(ch^, (n+2)*SizeOf(ArbFloat), 0); + Getmem(h, n*SizeOf(ArbFloat)); + Getmem(qty, n*SizeOf(ArbFloat)); + ncd := n-3; if ncd>2 then ncd := 2; + Getmem(qtdq, ((n-2)*(ncd+1)-(ncd*(ncd+1)) div 2)*SizeOf(ArbFloat)); + for i:=2 to n do h^[i] := 1/(xyw[i].x-xyw[i-1].x); h^[1] := 0; + for i:=1 to n-2 + do qty^[i] := (h^[i+1]*xyw[i].y - + (h^[i+1]+h^[i+2])*xyw[i+1].y + + h^[i+2]*xyw[i+2].y); + j := 1; i := 1; + qtdq^[j] := sqr(h^[i+1])/xyw[i].w + + sqr(h^[i+1]+h^[i+2])/xyw[i+1].w + + sqr(h^[i+2])/xyw[i+2].w + + lambda*(1/h^[i+1]+1/h^[i+2])/3; + Inc(j); + if ncd>0 then + begin i := 2; + qtdq^[j] := -h^[i+1]*(h^[i]+h^[i+1])/xyw[i].w + -h^[i+1]*(h^[i+1]+h^[i+2])/xyw[i+1].w + + lambda/h^[i+1]/6; + Inc(j); + qtdq^[j] := sqr(h^[i+1])/xyw[i].w + + sqr(h^[i+1]+h^[i+2])/xyw[i+1].w + + sqr(h^[i+2])/xyw[i+2].w + + lambda*(1/h^[i+1]+1/h^[i+2])/3; + Inc(j) + end; + for i:=3 to n-2 + do begin + qtdq^[j] := h^[i]*h^[i+1]/xyw[i].w; + Inc(j); + qtdq^[j] := -h^[i+1]*(h^[i]+h^[i+1])/xyw[i].w + -h^[i+1]*(h^[i+1]+h^[i+2])/xyw[i+1].w + + lambda/h^[i+1]/6; + Inc(j); + qtdq^[j] := sqr(h^[i+1])/xyw[i].w + + sqr(h^[i+1]+h^[i+2])/xyw[i+1].w + + sqr(h^[i+2])/xyw[i+2].w + + lambda*(1/h^[i+1]+1/h^[i+2])/3; + Inc(j) + end; + { Solving for c/lambda } + Slegpb(n-2, ncd, qtdq^[1], qty^[1], ch^[2], ca, term); + if term=1 then + begin + residu := 0; + for i:=1 to n do + begin + crow := (h^[i]*ch^[i-1] - (h^[i]+h^[i+1])*ch^[i]+h^[i+1]*ch^[i+1]) + /xyw[i].w; + xac[i].y := xyw[i].y - crow; + residu := residu + sqr(crow)*xyw[i].w + end; + xac[1].w := 0; + for i:=2 to n-1 do xac[i].w := lambda*ch^[i]; + xac[n].w := 0; + end; + Freemem(qtdq, ((n-2)*(ncd+1)-(ncd*(ncd+1)) div 2)*SizeOf(ArbFloat)); + Freemem(qty, n*SizeOf(ArbFloat)); + Freemem(h, n*SizeOf(ArbFloat)); + Freemem(ch, (n+2)*SizeOf(ArbFloat)); +end; + + +procedure spl2nalf(n: ArbInt; var xyzw1: ArbFloat; lambda:ArbFloat; + var xyg0, residu: ArbFloat; var term: ArbInt); +type R3 = array[1..3] of ArbFloat; + R33= array[1..3] of R3; + Rn3= array[1..$ffe0 div SizeOf(R3)] of R3; + +var b,e21t,ht :^Rn3; + pfac :par2dr1; + e22 :R33; + i,j,l,i1,i2,n3 :ArbInt; + s,s1,px,py,hr,ca, + x,absdet,x1,x2, + absdetmax :ArbFloat; + vr :R4x; + wr :R2; + w,u :R3; + a_alfa_d :R4xAr absolute xyzw1; + a_gamma :nsp2rec absolute xyg0; + gamma :^arfloat1; + + + function e(var x,y:R2):ArbFloat; + const c1:ArbFloat=1/(16*pi); + var s:ArbFloat; + begin s:=sqr(x[1]-y[1]) +sqr(x[2]-y[2]); + if s=0 then e:=0 else e:=c1*s*ln(s) + end {e}; + + procedure pfxpfy(var a,b,c:R2;var f:r3; var pfx,pfy:ArbFloat); + var det:ArbFloat; + begin det:=(b[1]-a[1])*(c[2]-a[2]) - (b[2]-a[2])*(c[1]-a[1]); + pfx:=((f[2]-f[1])*(c[2]-a[2]) - (f[3]-f[1])*(b[2]-a[2]))/det; + pfy:=(-(f[2]-f[1])*(c[1]-a[1]) + (f[3]-f[1])*(b[1]-a[1]))/det + end {pfxpfy}; + + procedure pxpy(var a,b,c:R2; var px,py:ArbFloat); + var det : ArbFloat; + begin det:=(b[1]-a[1])*(c[2]-a[2]) - (b[2]-a[2])*(c[1]-a[1]); + px:=(b[2]-c[2])/det; py:=(c[1]-b[1])/det + end {pxpy}; + + function p(var x,a:R2; var px,py:ArbFloat):ArbFloat; + begin p:=1 + (x[1]-a[1])*px +(x[2]-a[2])*py end {p}; + + procedure slegpdlown(n: ArbInt; var a1; var bx1: ArbFloat; + var term: ArbInt); + var i, j, k, kmin1 : ArbInt; + h, lkk : ArbFloat; + a : ar2dr1 absolute a1; + x : arfloat1 absolute bx1; + begin + k:=0; term := 2; + while (k<n) do + begin + kmin1:=k; k:=k+1; lkk:=a[k]^[k]; + for j:=1 to kmin1 do lkk:=lkk-sqr(a[k]^[j]); + if lkk<=0 then exit else + begin + a[k]^[k]:=sqrt(lkk); lkk:=a[k]^[k]; + for i:=k+1 to n do + begin + h:=a[i]^[k]; + for j:=1 to kmin1 do h:=h-a[k]^[j]*a[i]^[j]; + a[i]^[k]:=h/lkk + end; {i} + h:=x[k]; + for j:=1 to kmin1 do h:=h-a[k]^[j]*x[j]; + x[k]:=h/lkk + end {lkk > 0} + end; {k} + for i:=n downto 1 do + begin + h:=x[i]; + for j:=i+1 to n do h:=h-a[j]^[i]*x[j]; + x[i]:=h/a[i]^[i]; + end; {i} + term := 1 + end; + +begin + term := 3; if n<3 then exit; + n3 := n - 3; + i1:=1; x1:=a_alfa_d[1].xy[1]; i2:=1; x2:=x1; + for i:= 2 to n do + begin hr:=a_alfa_d[i].xy[1]; + if hr < x1 then begin i1:=i; x1:=hr end else + if hr > x2 then begin i2:=i; x2:=hr end; + end; + vr:=a_alfa_d[n-2]; a_alfa_d[n-2]:=a_alfa_d[i1]; a_alfa_d[i1]:=vr; + vr:=a_alfa_d[n-1]; a_alfa_d[n-1]:=a_alfa_d[i2]; a_alfa_d[i2]:=vr; + + for i:=1 to 2 do vr.xy[i]:=a_alfa_d[n-2].xy[i]-a_alfa_d[n-1].xy[i]; + absdetmax:=-1; i1:=0; + for i:=1 to n do + begin for j:=1 to 2 do wr[j]:=a_alfa_d[i].xy[j]-a_alfa_d[n-2].xy[j]; + if a_alfa_d[i].d<=0 then exit; + absdet:=abs(wr[1]*vr.xy[2]-wr[2]*vr.xy[1]); + if absdet > absdetmax then begin i1:=i; absdetmax:=absdet end; + end; + term := 4; + if absdetmax<=macheps*abs(x2-x1) then exit; + term := 1; + vr:=a_alfa_d[n]; a_alfa_d[n]:=a_alfa_d[i1]; a_alfa_d[i1]:=vr; + GetMem(e21t, n3*SizeOf(r3)); + GetMem(b, n3*SizeOf(r3)); + GetMem(gamma, n*SizeOf(ArbFloat)); + + pxpy(a_alfa_d[n-2].xy,a_alfa_d[n-1].xy,a_alfa_d[n].xy,px,py); + for i:=1 to n3 do b^[i][1]:=p(a_alfa_d[i].xy,a_alfa_d[n-2].xy,px,py); + pxpy(a_alfa_d[n-1].xy,a_alfa_d[n].xy,a_alfa_d[n-2].xy,px,py); + for i:=1 to n3 do b^[i][2]:=p(a_alfa_d[i].xy,a_alfa_d[n-1].xy,px,py); + pxpy(a_alfa_d[n].xy,a_alfa_d[n-2].xy,a_alfa_d[n-1].xy,px,py); + for i:=1 to n3 do b^[i][3]:=p(a_alfa_d[i].xy,a_alfa_d[n].xy,px,py); + e22[1,1]:=0; e22[2,2]:=0; e22[3,3]:=0; + e22[2,1]:=e(a_alfa_d[n-1].xy,a_alfa_d[n-2].xy); e22[1,2]:=e22[2,1]; + e22[3,1]:=e(a_alfa_d[n].xy,a_alfa_d[n-2].xy); e22[1,3]:=e22[3,1]; + e22[3,2]:=e(a_alfa_d[n].xy,a_alfa_d[n-1].xy); e22[2,3]:=e22[3,2]; + for i:=1 to 3 do + for j:=1 to n3 do e21t^[j,i]:=e(a_alfa_d[n3+i].xy,a_alfa_d[j].xy); + + GetMem(ht, n3*SizeOf(r3)); + for i:=1 to 3 do + for j:=1 to n3 do + begin s:=0; + for l:= 1 to 3 do s:=s+e22[i,l]*b^[j][l]; ht^[j][i]:=s + end; + AllocateL2dr(n3,pfac); + for i:= 1 to n3 do + for j:= 1 to i do + begin if j=i then s1:=0 else s1:=e(a_alfa_d[i].xy,a_alfa_d[j].xy); + for l:= 1 to 3 do s1:=s1+b^[i][l]*(ht^[j][l]-e21t^[j][l])-e21t^[i][l]*b^[j][l]; + if j=i then s:=1/a_alfa_d[i].d else s:=0; + for l:= 1 to 3 do s:=s+b^[i][l]*b^[j][l]/a_alfa_d[n3+l].d; + pfac^[i]^[j] := s1+s/lambda + end; + for i:= 1 to n3 do + gamma^[i]:=a_alfa_d[i].alfa-b^[i][1]*a_alfa_d[n-2].alfa-b^[i][2]*a_alfa_d[n-1].alfa-b^[i][3]*a_alfa_d[n].alfa; + slegpdlown(n3,pfac^[1],gamma^[1],term); + DeAllocateL2dr(n3,pfac); + FreeMem(ht, n3*SizeOf(r3)); + + if term=1 then + begin + for i:= 1 to 3 do + begin s:= 0; + for j:= 1 to n3 do + s:=s+b^[j][i]*gamma^[j]; w[i]:=s; + gamma^[n3+i]:=-w[i] + end;{w=btgamma} + for i:=1 to 3 do + begin s:=0; + for l:=1 to n3 do s:=s+e21t^[l][i]*gamma^[l]; + s1:=0; + for l:=1 to 3 do s1:=s1+e22[i,l]*w[l]; + u[i]:=a_alfa_d[n3+i].alfa+w[i]/(lambda*a_alfa_d[n3+i].d)+s1-s + end; + with a_gamma[0] do + pfxpfy(a_alfa_d[n-2].xy,a_alfa_d[n-1].xy,a_alfa_d[n].xy,u,xy[1],xy[2]); + residu:=0;for i:=1 to n3 do residu:=residu+sqr(gamma^[i])/a_alfa_d[i].d; + for i:= 1 to 3 do residu:=residu+sqr(w[i])/a_alfa_d[n3+i].d; + residu:=residu/sqr(lambda); + a_gamma[0].gamma := u[1]; + for i:=1 to n do + begin + a_gamma[i].xy := a_alfa_d[i].xy; + a_gamma[i].gamma := gamma^[i] + end; + end; + FreeMem(gamma, n*SizeOf(ArbFloat)); + FreeMem(b, n3*SizeOf(r3)); + FreeMem(e21t, n3*SizeOf(r3)) + end; + +function spl2natv(n: ArbInt; var xyg0: ArbFloat; u, v: ArbFloat): ArbFloat; + +const c1: ArbFloat=1/(16*pi); + + var i : ArbInt; + s : ArbFloat; + a_gamma : nsp2rec absolute xyg0; + z : R2; + + function e(var x,y:R2):ArbFloat; + var s:ArbFloat; + begin + s:=sqr(x[1]-y[1]) + sqr(x[2]-y[2]); + if s=0 then + e:= 0 + else + e:=s*ln(s) + end {e}; + + function pf(var x,a:R2;fa,pfx,pfy:ArbFloat):ArbFloat; + begin + pf:=fa + (x[1]-a[1])*pfx + (x[2]-a[2])*pfy + end {pf}; + + begin + s:=0; + z[1] := u; + z[2] := v; + for i:=1 to n do + s:=s+a_gamma[i].gamma*e(z, a_gamma[i].xy); + with a_gamma[0] do + spl2natv :=s*c1+pf(z,a_gamma[n-2].xy, gamma, xy[1], xy[2]) + end; + +begin + +end. |