ipf.pas 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949
  1. {
  2. This file is part of the Numlib package.
  3. Copyright (c) 1986-2000 by
  4. Kees van Ginneken, Wil Kortsmit and Loek van Reij of the
  5. Computational centre of the Eindhoven University of Technology
  6. FPC port Code by Marco van de Voort ([email protected])
  7. documentation by Michael van Canneyt ([email protected])
  8. Interpolate and (curve) fitting.
  9. Slegpb in this unit patched parameters slightly. Units IPF and sle
  10. were not in the same revision in this numlib copy (which was a
  11. copy of the work directory of the author) .
  12. Contains two undocumented functions. If you recognize the algoritm,
  13. mail us.
  14. See the file COPYING.FPC, included in this distribution,
  15. for details about the copyright.
  16. This program is distributed in the hope that it will be useful,
  17. but WITHOUT ANY WARRANTY; without even the implied warranty of
  18. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  19. **********************************************************************}
  20. {
  21. }
  22. unit ipf;
  23. {$I direct.inc}
  24. interface
  25. uses typ, mdt, dsl, sle, spe;
  26. { Determine natural cubic spline "s" for data set (x,y), output to (a,d2a)
  27. term=1 success,
  28. =2 failure calculating "s"
  29. =3 wrong input (e.g. x,y is not sorted increasing on x)}
  30. procedure ipffsn(n: ArbInt; var x, y, a, d2a: ArbFloat; var term: ArbInt);
  31. {calculate d2s from x,y, which can be used to calculate s}
  32. procedure ipfisn(n: ArbInt; var x, y, d2s: ArbFloat; var term: ArbInt);
  33. {Calculate function value for dataset (x,y), with n.c. spline d2s for
  34. x value t. Return (corrected) y value.
  35. s calculated from x,y, with e.g. ipfisn}
  36. function ipfspn(n: ArbInt; var x, y, d2s: ArbFloat; t: ArbFloat;
  37. var term: ArbInt): ArbFloat;
  38. {Calculate minimum and maximum values for the n.c. spline d2s.
  39. Does NOT take source points into account.}
  40. procedure ipfsmm(n: ArbInt; var x, y, d2s, minv, maxv: ArbFloat;
  41. var term: ArbInt);
  42. {Calculate n-degree polynomal b for dataset (x,y) with m elements
  43. using the least squares method.}
  44. procedure ipfpol(m, n: ArbInt; var x, y, b: ArbFloat; var term: ArbInt);
  45. {**** undocumented ****}
  46. function spline( n : ArbInt;
  47. x : complex;
  48. var ac : complex;
  49. var gammar: ArbFloat;
  50. u1 : ArbFloat;
  51. pf : complex): ArbFloat;
  52. {**** undocumented ****}
  53. procedure splineparameters
  54. ( n : ArbInt;
  55. var ac, alfadc : complex;
  56. var lambda,
  57. gammar, u1,
  58. kwsom, energie : ArbFloat;
  59. var pf : complex);
  60. implementation
  61. procedure ipffsn(n: ArbInt; var x, y, a, d2a: ArbFloat; var term: ArbInt);
  62. var i, j, sr, n1s, ns1, ns2: ArbInt;
  63. s, lam, lam0, lam1, lambda, ey, ca, p, q, r: ArbFloat;
  64. px, py, pd, pa, pd2a,
  65. h, z, diagb, dinv, qty, qtdinvq, c, t, tl: ^arfloat1;
  66. ub: boolean;
  67. procedure solve; {n, py, qty, h, qtdinvq, dinv, lam, t, pa, pd2a, term}
  68. var i: ArbInt;
  69. p, q, r, ca: ArbFloat;
  70. f, c: ^arfloat1;
  71. begin
  72. getmem(f, 3*ns1); getmem(c, ns1);
  73. for i:=1 to n-1 do
  74. begin
  75. f^[3*i]:=qtdinvq^[3*i]+lam*t^[2*i];
  76. if i > 1
  77. then
  78. f^[3*i-1]:=qtdinvq^[3*i-1]+lam*t^[2*i-1];
  79. if i > 2
  80. then
  81. f^[3*i-2]:=qtdinvq^[3*i-2];
  82. if lam=0
  83. then
  84. c^[i]:=qty^[i]
  85. else
  86. c^[i]:=lam*qty^[i]
  87. end;
  88. slegpb(n-1, 2,{ 3,} f^[1], c^[1], pd2a^[1], ca, term);
  89. if term=2
  90. then
  91. begin
  92. freemem(f, 3*ns1); freemem(c, ns1);
  93. exit
  94. end;
  95. p:=1/h^[1];
  96. if lam=0
  97. then
  98. r:=1
  99. else
  100. r:=1/lam;
  101. q:=1/h^[2]; pa^[1]:=py^[1]-r*dinv^[1]*p*pd2a^[1];
  102. pa^[2]:=py^[2]-r*dinv^[2]*(pd2a^[2]*q-(p+q)*pd2a^[1]); p:=q;
  103. for i:=3 to n-1 do
  104. begin
  105. q:=1/h^[i];
  106. pa^[i]:=py^[i]-r*dinv^[i]*
  107. (p*pd2a^[i-2]-(p+q)*pd2a^[i-1]+q*pd2a^[i]);
  108. p:=q
  109. end;
  110. q:=1/h^[n];
  111. pa^[n]:=py^[n]-r*dinv^[n]*(p*pd2a^[n-2]-(p+q)*pd2a^[n-1]);
  112. pa^[n+1]:=py^[n+1]-r*dinv^[n+1]*q*pd2a^[n-1];
  113. if lam=0
  114. then
  115. for i:=1 to n-1 do
  116. pd2a^[i]:=0;
  117. freemem(f, 3*ns1); freemem(c, ns1);
  118. end; {solve}
  119. function e(var c, h: ArbFloat; n:ArbInt): ArbFloat;
  120. var i:ArbInt;
  121. s:ArbFloat;
  122. pc, ph: ^arfloat1;
  123. begin
  124. ph:=@h; pc:=@c;
  125. s:=ph^[1]*pc^[1]*pc^[1];
  126. for i:=1 to n-2 do
  127. s:=s+(pc^[i]*(pc^[i]+pc^[i+1])+pc^[i+1]*pc^[i+1])*ph^[i+1];
  128. e:=(s+pc^[n-1]*pc^[n-1]*ph^[n])/3
  129. end; {e}
  130. function cr(lambda: ArbFloat): ArbFloat;
  131. var s, crs: ArbFloat;
  132. i: ArbInt;
  133. begin
  134. cr:=0; lam:=lambda;
  135. solve; { n, py, qty, h, qtdinvq, dinv, lam, t, pa, pd2a, term }
  136. if term=2
  137. then
  138. exit;
  139. crs:=ey;
  140. if lam <> 0
  141. then
  142. begin
  143. crs:=crs+e(pd2a^[1], h^[1], n);
  144. s:=0;
  145. for i:=1 to n-1 do
  146. s:=s+pd2a^[i]*qty^[i];
  147. crs:=crs-2*s
  148. end;
  149. s:=0;
  150. for i:=1 to n+1 do
  151. s:=s+sqr(pa^[i]-py^[i])*diagb^[i];
  152. cr:=crs-s
  153. end; {cr}
  154. procedure roof1r(a, b, ae, re: ArbFloat; var x: ArbFloat);
  155. var fa, fb, c, fc, m, tol, w1, w2 : ArbFloat;
  156. k : ArbInt;
  157. stop : boolean;
  158. begin
  159. fa:=cr(a);
  160. if term=2
  161. then
  162. exit;
  163. fb:=cr(b);
  164. if term=2
  165. then
  166. exit;
  167. if abs(fb)>abs(fa)
  168. then
  169. begin
  170. c:=b; fc:=fb; x:=a; b:=a; fb:=fa; a:=c; fa:=fc
  171. end
  172. else
  173. begin
  174. c:=a; fc:=fa; x:=b
  175. end;
  176. k:=0;
  177. tol:=ae+re*spemax(abs(a), abs(b));
  178. w1:=abs(b-a); stop:=false;
  179. while (abs(b-a)>tol) and (fb<>0) and (not stop) do
  180. begin
  181. m:=(a+b)/2;
  182. if (k>=2) or (fb=fc)
  183. then
  184. x:=m
  185. else
  186. begin
  187. x:=(b*fc-c*fb)/(fc-fb);
  188. if abs(b-x)<tol
  189. then
  190. x:=b-tol*spesgn(b-a);
  191. if spesgn(x-m)=spesgn(x-b)
  192. then
  193. x:=m
  194. end;
  195. c:=b; fc:=fb; b:=x; fb:=cr(x);
  196. if term=2
  197. then
  198. exit;
  199. if spesgn(fa)*spesgn(fb)>0
  200. then
  201. begin
  202. a:=c; fa:=fc; k:=0
  203. end
  204. else
  205. k:=k+1;
  206. if abs(fb)>=abs(fa)
  207. then
  208. begin
  209. c:=b; fc:=fb; x:=a; b:=a; fb:=fa; a:=c; fa:=fc; k:=0
  210. end;
  211. tol:=ae+re*spemax(abs(a), abs(b));
  212. w2:=abs(b-a);
  213. if w2>=w1
  214. then
  215. stop:=true;
  216. w1:=w2
  217. end
  218. end; {roof1r}
  219. procedure NoodGreep;
  220. var I, j: ArbInt;
  221. begin
  222. i:=1;
  223. while i <= n do
  224. begin
  225. if (pd^[i] <= 0) or (px^[i+1] <= px^[i])
  226. then
  227. begin
  228. term:=3;
  229. exit
  230. end;
  231. i:=i+1
  232. end;
  233. if pd^[n+1] <= 0
  234. then
  235. begin
  236. term:=3;
  237. exit
  238. end;
  239. for i:=1 to n+1 do
  240. dinv^[i]:=1/pd^[i];
  241. for i:=1 to n do
  242. h^[i]:=px^[i+1]-px^[i];
  243. t^[2]:=(h^[1]+h^[2])/3;
  244. for i:=2 to n-1 do
  245. begin
  246. t^[2*i]:=(h^[i]+h^[i+1])/3; t^[2*i-1]:=h^[i]/6
  247. end;
  248. move(t^[1], tl^[1], ns2);
  249. mdtgpb(n-1, 1, 2, tl^[1], ca, term);
  250. if term=2
  251. then
  252. exit;
  253. z^[1]:=1/(h^[1]*tl^[2]);
  254. for j:=2 to n-1 do
  255. z^[j]:=-(tl^[2*j-1]*z^[j-1])/tl^[2*j];
  256. s:=0;
  257. for j:=1 to n-1 do
  258. s:=s+sqr(z^[j]);
  259. diagb^[1]:=s;
  260. z^[1]:=(-1/h^[1]-1/h^[2])/tl^[2];
  261. if n>2
  262. then
  263. z^[2]:=(1/h^[2]-tl^[3]*z^[1])/tl^[4];
  264. for j:=3 to n-1 do
  265. z^[j]:=-tl^[2*j-1]*z^[j-1]/tl^[2*j];
  266. s:=0;
  267. for j:=1 to n-1 do
  268. s:=s+sqr(z^[j]);
  269. diagb^[2]:=s;
  270. for i:=2 to n-2 do
  271. begin
  272. z^[i-1]:=1/(h^[i]*tl^[2*(i-1)]);
  273. z^[i]:=(-1/h^[i]-1/h^[i+1]-tl^[2*i-1]*z^[i-1])/tl^[2*i];
  274. z^[i+1]:=(1/h^[i+1]-tl^[2*i+1]*z^[i])/tl^[2*(i+1)];
  275. for j:=i+2 to n-1 do
  276. z^[j]:=-tl^[2*j-1]*z^[j-1]/tl^[2*j];
  277. s:=0;
  278. for j:=i-1 to n-1 do
  279. s:=s+sqr(z^[j]);
  280. diagb^[i+1]:=s
  281. end;
  282. z^[n-2]:=1/(h^[n-1]*tl^[2*(n-2)]);
  283. z^[n-1]:=(-1/h^[n-1]-1/h^[n]-tl^[2*n-3]*z^[n-2])/tl^[2*(n-1)];
  284. s:=0;
  285. for j:=n-2 to n-1 do
  286. s:=s+sqr(z^[j]);
  287. diagb^[n]:=s;
  288. diagb^[n+1]:=1/sqr(h^[n]*tl^[2*(n-1)]);
  289. p:=1/h^[1];
  290. for i:=2 to n do
  291. begin
  292. q:=1/h^[i]; qty^[i-1]:=py^[i+1]*q-py^[i]*(p+q)+py^[i-1]*p;
  293. p:=q
  294. end;
  295. p:=1/h^[1]; q:=1/h^[2]; r:=1/h^[3];
  296. qtdinvq^[3]:=dinv^[1]*p*p+dinv^[2]*(p+q)*(p+q)+dinv^[3]*q*q;
  297. if n>2
  298. then
  299. begin
  300. qtdinvq^[6]:=dinv^[2]*q*q+dinv^[3]*(q+r)*(q+r)+dinv^[4]*r*r;
  301. qtdinvq^[5]:=-(dinv^[2]*(p+q)+dinv^[3]*(q+r))*q;
  302. p:=q; q:=r;
  303. for i:=3 to n-1 do
  304. begin
  305. r:=1/h^[i+1];
  306. qtdinvq^[3*i]:=dinv^[i]*q*q+dinv^[i+1]*(q+r)*(q+r)+dinv^[i+2]*r*r;
  307. qtdinvq^[3*i-1]:=-(dinv^[i]*(p+q)+dinv^[i+1]*(q+r))*q;
  308. qtdinvq^[3*i-2]:=dinv^[i]*p*q;
  309. p:=q; q:=r
  310. end
  311. end;
  312. dslgpb(n-1, 1, 2, tl^[1], qty^[1], c^[1], term);
  313. if term=2
  314. then
  315. exit;
  316. ey:=e(c^[1], h^[1], n);
  317. lam0:=0;
  318. s:=cr(lam0);
  319. if term=2
  320. then
  321. exit;
  322. if s >= 0
  323. then
  324. begin
  325. lambda:=0; term:=4
  326. end
  327. else
  328. begin
  329. lam1:=1e-8; ub:=false;
  330. while (not ub) and (lam1<=1.1e8) do
  331. begin
  332. s:=cr(lam1);
  333. if term=2
  334. then
  335. exit;
  336. if s >= 0
  337. then
  338. ub:=true
  339. else
  340. begin
  341. lam0:=lam1; lam1:=10*lam1
  342. end
  343. end;
  344. if not ub
  345. then
  346. begin
  347. term:=4; lambda:=lam0
  348. end
  349. else
  350. roof1r(lam0, lam1, 0, 1e-6, lambda);
  351. if term=2
  352. then
  353. exit
  354. end;
  355. end;
  356. begin
  357. term:=1;
  358. if n < 2
  359. then
  360. begin
  361. term:=3; exit
  362. end;
  363. sr:=sizeof(ArbFloat);
  364. n1s:=(n+1)*sr;
  365. ns2:=2*(n-1)*sr;
  366. ns1:=(n-1)*sr;
  367. getmem(dinv, n1s);
  368. getmem(h, n*sr);
  369. getmem(t, ns2);
  370. getmem(tl, ns2);
  371. getmem(z, ns1);
  372. getmem(diagb, n1s);
  373. getmem(qtdinvq, 3*ns1);
  374. getmem(c, ns1);
  375. getmem(qty, ns1);
  376. getmem(pd, n1s);
  377. { pd:=@d; }
  378. px:=@x;
  379. py:=@y;
  380. pa:=@a;
  381. pd2a:=@d2a;
  382. { de gewichten van de punten worden op 1 gezet}
  383. for i:=1 to n+1 do
  384. pd^[i]:=1;
  385. NoodGreep;
  386. freemem(dinv, n1s);
  387. freemem(h, n*sr);
  388. freemem(t, ns2);
  389. freemem(tl, ns2);
  390. freemem(z, ns1);
  391. freemem(diagb, n1s);
  392. freemem(qtdinvq, 3*ns1);
  393. freemem(c, ns1);
  394. freemem(qty, ns1);
  395. freemem(pd, n1s);
  396. end; {ipffsn}
  397. procedure ortpol(m, n: ArbInt; var x, alfa, beta: ArbFloat);
  398. // this function used to use mark/release.
  399. var
  400. i, j, ms : ArbInt;
  401. xppn1, ppn1, ppn, p, alfaj, betaj : ArbFloat;
  402. px, pal, pbe, pn, pn1 : ^arfloat1;
  403. begin
  404. px:=@x; pal:=@alfa; pbe:=@beta; ms:=m*sizeof(ArbFloat);
  405. getmem(pn, ms); getmem(pn1, ms);
  406. xppn1:=0; ppn1:=m;
  407. for i:=1 to m do
  408. begin
  409. pn^[i]:=0; pn1^[i]:=1; xppn1:=xppn1+px^[i]
  410. end;
  411. pal^[1]:=xppn1/ppn1; pbe^[1]:=0;
  412. for j:=2 to n do
  413. begin
  414. alfaj:=pal^[j-1]; betaj:=pbe^[j-1];
  415. ppn:=ppn1; ppn1:=0; xppn1:=0;
  416. for i:=1 to m do
  417. begin
  418. p:=(px^[i]-alfaj)*pn1^[i]-betaj*pn^[i];
  419. pn^[i]:=pn1^[i]; pn1^[i]:=p; p:=p*p;
  420. ppn1:=ppn1+p; xppn1:=xppn1+px^[i]*p
  421. end; {i}
  422. pal^[j]:=xppn1/ppn1; pbe^[j]:=ppn1/ppn
  423. end; {j}
  424. freemem(pn); freemem(pn1);
  425. end; {ortpol}
  426. procedure ortcoe(m, n: ArbInt; var x, y, alfa, beta, a: ArbFloat);
  427. // this function used to use mark/release.
  428. var i, j, mr : ArbInt;
  429. fpn, ppn, p, alphaj, betaj : ArbFloat;
  430. px, py, pal, pbe, pa, pn, pn1 : ^arfloat1;
  431. begin
  432. mr:=m*sizeof(ArbFloat);
  433. px:=@x; py:=@y; pal:=@alfa; pbe:=@beta; pa:=@a;
  434. getmem(pn, mr); getmem(pn1, mr);
  435. fpn:=0;
  436. for i:=1 to m do
  437. begin
  438. pn^[i]:=0; pn1^[i]:=1; fpn:=fpn+py^[i]
  439. end; {i}
  440. pa^[1]:=fpn/m;
  441. for j:=1 to n do
  442. begin
  443. fpn:=0; ppn:=0; alphaj:=pal^[j]; betaj:=pbe^[j];
  444. for i:=1 to m do
  445. begin
  446. p:=(px^[i]-alphaj)*pn1^[i]-betaj*pn^[i];
  447. pn^[i]:=pn1^[i]; pn1^[i]:=p;
  448. fpn:=fpn+py^[i]*p; ppn:=ppn+p*p
  449. end; {i}
  450. pa^[j+1]:=fpn/ppn
  451. end; {j}
  452. freemem(pn); freemem(pn1);
  453. end; {ortcoe}
  454. procedure polcoe(n:ArbInt; var alfa, beta, a, b: ArbFloat);
  455. var k, j : ArbInt;
  456. pal, pbe : ^arfloat1;
  457. pa, pb : ^arfloat0;
  458. begin
  459. pal:=@alfa; pbe:=@beta; pa:=@a; pb:=@b;
  460. move(pa^[0], pb^[0], (n+1)*sizeof(ArbFloat));
  461. for j:=0 to n-1 do
  462. for k:=n-j-1 downto 0 do
  463. begin
  464. pb^[k+j]:=pb^[k+j]-pal^[k+1]*pb^[k+j+1];
  465. if k+j<>n-1
  466. then
  467. pb^[k+j]:=pb^[k+j]-pbe^[k+2]*pb^[k+j+2]
  468. end
  469. end; {polcoe}
  470. procedure ipfpol(m, n: ArbInt; var x, y, b: ArbFloat; var term: ArbInt);
  471. var i, ns: ArbInt;
  472. fsum: ArbFloat;
  473. px, py, alfa, beta: ^arfloat1;
  474. pb, a: ^arfloat0;
  475. begin
  476. if (n<0) or (m<1)
  477. then
  478. begin
  479. term:=3; exit
  480. end;
  481. term:=1;
  482. if n = 0
  483. then
  484. begin
  485. py:=@y; pb:=@b;
  486. fsum:=0;
  487. for i:=1 to m do
  488. fsum:=fsum+py^[i];
  489. pb^[0]:=fsum/m
  490. end
  491. else
  492. begin
  493. if n>m-1
  494. then
  495. begin
  496. pb:=@b;
  497. fillchar(pb^[m], (n-m+1)*sizeof(ArbFloat), 0);
  498. n:=m-1
  499. end;
  500. ns:=n*sizeof(ArbFloat);
  501. getmem(alfa, ns); getmem(beta, ns);
  502. getmem(a, (n+1)*sizeof(ArbFloat));
  503. ortpol(m, n, x, alfa^[1], beta^[1]);
  504. ortcoe(m, n, x, y, alfa^[1], beta^[1], a^[0]);
  505. polcoe(n, alfa^[1], beta^[1], a^[0], b);
  506. freemem(alfa, ns); freemem(beta, ns);
  507. freemem(a, (n+1)*sizeof(ArbFloat));
  508. end
  509. end; {ipfpol}
  510. procedure ipfisn(n: ArbInt; var x, y, d2s: ArbFloat; var term: ArbInt);
  511. var
  512. s, i : ArbInt;
  513. p, q, ca : ArbFloat;
  514. px, py, h, b, t : ^arfloat0;
  515. pd2s : ^arfloat1;
  516. begin
  517. px:=@x; py:=@y; pd2s:=@d2s;
  518. term:=1;
  519. if n < 2
  520. then
  521. begin
  522. term:=3; exit
  523. end; {n<2}
  524. s:=sizeof(ArbFloat);
  525. getmem(h, n*s);
  526. getmem(b, (n-1)*s);
  527. getmem(t, 2*(n-1)*s);
  528. for i:=0 to n-1 do
  529. h^[i]:=px^[i+1]-px^[i];
  530. q:=1/6; p:=2*q;
  531. t^[1]:=p*(h^[0]+h^[1]);
  532. for i:=2 to n-1 do
  533. begin
  534. t^[2*i-1]:=p*(h^[i-1]+h^[i]); t^[2*i-2]:=q*h^[i-1]
  535. end; {i}
  536. p:=1/h^[0];
  537. for i:=2 to n do
  538. begin
  539. q:=1/h^[i-1]; b^[i-2]:=py^[i]*q-py^[i-1]*(p+q)+py^[i-2]*p; p:=q
  540. end;
  541. slegpb(n-1, 1, {2,} t^[1], b^[0], pd2s^[1], ca, term);
  542. freemem(h, n*s);
  543. freemem(b, (n-1)*s);
  544. freemem(t, 2*(n-1)*s);
  545. end; {ipfisn}
  546. function ipfspn(n: ArbInt; var x, y, d2s: ArbFloat; t:ArbFloat;
  547. var term: ArbInt): ArbFloat;
  548. var
  549. px, py : ^arfloat0;
  550. pd2s : ^arfloat1;
  551. i, j, m : ArbInt;
  552. d, s3, h, dy : ArbFloat;
  553. begin
  554. i:=1; term:=1;
  555. if n<2
  556. then
  557. begin
  558. term:=3; exit
  559. end; {n<2}
  560. px:=@x; py:=@y; pd2s:=@d2s;
  561. if t <= px^[0]
  562. then
  563. begin
  564. h:=px^[1]-px^[0];
  565. dy:=(py^[1]-py^[0])/h-h*pd2s^[1]/6;
  566. ipfspn:=py^[0]+(t-px^[0])*dy
  567. end { t <= x[0] }
  568. else
  569. if t >= px^[n]
  570. then
  571. begin
  572. h:=px^[n]-px^[n-1];
  573. dy:=(py^[n]-py^[n-1])/h+h*pd2s^[n-1]/6;
  574. ipfspn:=py^[n]+(t-px^[n])*dy
  575. end { t >= x[n] }
  576. else
  577. begin
  578. i:=0; j:=n;
  579. while j <> i+1 do
  580. begin
  581. m:=(i+j) div 2;
  582. if t>=px^[m]
  583. then
  584. i:=m
  585. else
  586. j:=m
  587. end; {j}
  588. h:=px^[i+1]-px^[i];
  589. d:=t-px^[i];
  590. if i=0
  591. then
  592. begin
  593. s3:=pd2s^[1]/h;
  594. dy:=(py^[1]-py^[0])/h-h*pd2s^[1]/6;
  595. ipfspn:=py^[0]+d*(dy+d*d*s3/6)
  596. end
  597. else
  598. if i=n-1
  599. then
  600. begin
  601. s3:=-pd2s^[n-1]/h;
  602. dy:=(py^[n]-py^[n-1])/h-h*pd2s^[n-1]/3;
  603. ipfspn:=py^[n-1]+d*(dy+d*(pd2s^[n-1]/2+d*s3/6))
  604. end
  605. else
  606. begin
  607. s3:=(pd2s^[i+1]-pd2s^[i])/h;
  608. dy:=(py^[i+1]-py^[i])/h-h*(2*pd2s^[i]+pd2s^[i+1])/6;
  609. ipfspn:=py^[i]+d*(dy+d*(pd2s^[i]/2+d*s3/6))
  610. end
  611. end { x[0] < t < x[n] }
  612. end; {ipfspn}
  613. procedure ipfsmm(
  614. n: ArbInt; var x, y, d2s, minv, maxv: ArbFloat; var term: ArbInt);
  615. var
  616. i: ArbInt;
  617. h: ArbFloat;
  618. px, py: ^arfloat0;
  619. pd2s: ^arfloat1;
  620. procedure UpdateMinMax(v: ArbFloat);
  621. begin
  622. if (0 >= v) or (v >= h) then exit;
  623. v := ipfspn(n, x, y, d2s, px^[i]+v, term);
  624. if v < minv then
  625. minv := v;
  626. if v > maxv then
  627. maxv := v;
  628. end;
  629. procedure MinMaxOnSegment;
  630. var
  631. a, b, c: ArbFloat;
  632. d: ArbFloat;
  633. begin
  634. h:=px^[i+1]-px^[i];
  635. if i=0
  636. then
  637. begin
  638. a:=pd2s^[1]/h/2;
  639. b:=0;
  640. c:=(py^[1]-py^[0])/h-h*pd2s^[1]/6;
  641. end
  642. else
  643. if i=n-1
  644. then
  645. begin
  646. a:=-pd2s^[n-1]/h/2;
  647. b:=pd2s^[n-1];
  648. c:=(py^[n]-py^[n-1])/h-h*pd2s^[n-1]/3;
  649. end
  650. else
  651. begin
  652. a:=(pd2s^[i+1]-pd2s^[i])/h/2;
  653. b:=pd2s^[i];
  654. c:=(py^[i+1]-py^[i])/h-h*(2*pd2s^[i]+pd2s^[i+1])/6;
  655. end;
  656. if a=0 then exit;
  657. d := b*b-4*a*c;
  658. if d<0 then exit;
  659. d:=Sqrt(d);
  660. UpdateMinMax((-b+d)/(2*a));
  661. UpdateMinMax((-b-d)/(2*a));
  662. end;
  663. begin
  664. term:=1;
  665. if n<2 then begin
  666. term:=3;
  667. exit;
  668. end;
  669. px:=@x; py:=@y; pd2s:=@d2s;
  670. for i:=0 to n-1 do
  671. MinMaxOnSegment;
  672. end;
  673. function p(x, a, z:complex): ArbFloat;
  674. begin
  675. x.sub(a);
  676. p:=x.Inp(z)
  677. end;
  678. function e(x, y: complex): ArbFloat;
  679. const c1: ArbFloat = 0.01989436788646;
  680. var s: ArbFloat;
  681. begin x.sub(y);
  682. s := x.norm;
  683. if s=0 then e:=0 else e:=c1*s*ln(s)
  684. end;
  685. function spline( n : ArbInt;
  686. x : complex;
  687. var ac : complex;
  688. var gammar: ArbFloat;
  689. u1 : ArbFloat;
  690. pf : complex): ArbFloat;
  691. var i : ArbInt;
  692. s : ArbFloat;
  693. a : arcomp0 absolute ac;
  694. gamma : arfloat0 absolute gammar;
  695. begin
  696. s := u1 + p(x, a[n-2], pf);
  697. for i:=0 to n do s := s + gamma[i]*e(x,a[i]);
  698. spline := s
  699. end;
  700. procedure splineparameters
  701. ( n : ArbInt;
  702. var ac, alfadc : complex;
  703. var lambda,
  704. gammar, u1,
  705. kwsom, energie : ArbFloat;
  706. var pf : complex);
  707. procedure SwapC(var v, w: complex);
  708. var x: complex;
  709. begin
  710. x := v; v := w; w := x
  711. end;
  712. procedure pxpy(a, b, c: complex; var p:complex);
  713. var det: ArbFloat;
  714. begin
  715. b.sub(a); c.sub(a); det := b.xreal*c.imag-b.imag*c.xreal;
  716. b.sub(c); p.Init(b.imag/det, -b.xreal/det)
  717. end;
  718. procedure pfxpfy(a, b, c: complex; f: vector; var pf: complex);
  719. begin
  720. b.sub(a); c.sub(a);
  721. f.j := f.j-f.i; f.k := f.k-f.i;
  722. pf.init(f.j*c.imag - f.k*b.imag, -f.j*c.xreal + f.k*b.xreal);
  723. pf.scale(1/(b.xreal*c.imag - b.imag*c.xreal))
  724. end;
  725. function InpV(n: ArbInt; var v1, v2: ArbFloat): ArbFloat;
  726. var i: ArbInt;
  727. a1: arfloat0 absolute v1;
  728. a2: arfloat0 absolute v2;
  729. s : ArbFloat;
  730. begin
  731. s := 0;
  732. for i:=0 to n-1 do s := s + a1[i]*a2[i];
  733. InpV := s
  734. end;
  735. PROCEDURE SPDSOL( N : INTEGER;
  736. VAR AP : pointer;
  737. VAR B : ArbFloat);
  738. VAR I, J, K : INTEGER;
  739. H : ArbFloat;
  740. a : ^ar2dr absolute ap;
  741. bx : arfloat0 absolute b;
  742. BEGIN
  743. for k:=0 to n do
  744. BEGIN
  745. h := sqrt(a^[k]^[k]-InpV(k, a^[k]^[0], a^[k]^[0]));
  746. a^[k]^[k] := h;
  747. FOR I:=K+1 TO N do a^[i]^[k] := (a^[i]^[k] - InpV(k, a^[k]^[0], a^[i]^[0]))/h;
  748. BX[K] := (bx[k] - InpV(k, a^[k]^[0], bx[0]))/h
  749. END;
  750. FOR I:=N DOWNTO 0 do
  751. BEGIN
  752. H := BX[I];
  753. FOR J:=I+1 TO N DO H := H - A^[J]^[I]*BX[J];
  754. BX[I] := H/A^[I]^[I]
  755. END
  756. END;
  757. var i, j, i1 : ArbInt;
  758. x, h,
  759. absdet,
  760. absdetmax,
  761. s, s1, ca: ArbFloat;
  762. alfa, dv, hulp,
  763. u, v, w : vector;
  764. e22 : array[0..2] of vector;
  765. e21, b : ^arvect0;
  766. k, c : ^ar2dr;
  767. gamma : arfloat0 absolute gammar;
  768. an2, an1, an, z,
  769. vz, wz : complex;
  770. a : arcomp0 absolute ac;
  771. alfad : arcomp0 absolute alfadc;
  772. begin
  773. i1:=0;
  774. x:=a[0].xreal;
  775. for i:=1 to n do
  776. begin
  777. h:=a[i].xreal;
  778. if h<x then begin i1:=i; x:=h end
  779. end;
  780. SwapC(a[n-2], a[i1]);
  781. SwapC(alfad[n-2], alfad[i1]);
  782. x:=a[0].xreal;
  783. i1 := 0;
  784. for i:=1 to n do
  785. begin
  786. h:=a[i].xreal;
  787. if h>x then begin i1:=i; x:=h end
  788. end;
  789. SwapC(a[n-1], a[i1]);
  790. SwapC(alfad[n-1], alfad[i1]);
  791. vz := a[n-2]; vz.sub(a[n-1]);
  792. absdetmax := -1;
  793. for i:=0 to n do
  794. begin
  795. wz := a[i]; wz.sub(a[n-2]);
  796. absdet := abs(wz.imag*vz.xreal-wz.xreal*vz.imag);
  797. if absdet>absdetmax then begin i1:=i; absdetmax:=absdet end
  798. end;
  799. SwapC(a[n], a[i1]);
  800. SwapC(alfad[n], alfad[i1]);
  801. an2 := a[n-2]; an1 := a[n-1]; an := a[n];
  802. alfa.i := alfad[n-2].xreal; dv.i := alfad[n-2].imag;
  803. alfa.j := alfad[n-1].xreal; dv.j := alfad[n-1].imag;
  804. alfa.k := alfad[n ].xreal; dv.k := alfad[n ].imag;
  805. n := n - 3;
  806. GetMem(k, (n+1)*SizeOf(pointer));
  807. for j:=0 to n do GetMem(k^[j], (j+1)*SizeOf(ArbFloat));
  808. GetMem(e21, (n+1)*SizeOf(vector));
  809. GetMem(b, (n+1)*SizeOf(vector));
  810. pxpy(an2,an1,an,z); for i:=0 to n do b^[i].i:=1+p(a[i],an2,z);
  811. pxpy(an1,an,an2,z); for i:=0 to n do b^[i].j:=1+p(a[i],an1,z);
  812. pxpy(an,an2,an1,z); for i:=0 to n do b^[i].k:=1+p(a[i],an,z);
  813. e22[0].init(0,e(an1,an2),e(an,an2));
  814. e22[1].init(e(an1,an2),0,e(an,an1));
  815. e22[2].init(e(an,an2),e(an,an1),0);
  816. for j:=0 to n do e21^[j].init(e(an2,a[j]),e(an1,a[j]),e(an,a[j]));
  817. GetMem(c, (n+1)*SizeOf(pointer));
  818. for j:=0 to n do GetMem(c^[j], (j+1)*SizeOf(ArbFloat));
  819. for i:=0 to n do
  820. for j:=0 to i do
  821. begin
  822. if j=i then s:=0 else s:=e(a[i],a[j]);
  823. hulp.init(b^[j].Inprod(e22[0]), b^[j].Inprod(e22[1]), b^[j].Inprod(e22[2]));
  824. hulp.sub(e21^[j]);
  825. k^[i]^[j] := s+b^[i].InProd(hulp)-b^[j].Inprod(e21^[i]);
  826. if j=i then s:=1/alfad[i].imag else s:=0;
  827. hulp.init(b^[j].i/dv.i, b^[j].j/dv.j, b^[j].k/dv.k);
  828. c^[i]^[j] := k^[i]^[j] + (s + b^[i].Inprod(hulp))/lambda
  829. end;
  830. for i:=0 to n do gamma[i]:=alfad[i].xreal - b^[i].Inprod(alfa);
  831. SpdSol(n, pointer(c), gamma[0]);
  832. for j:=n downto 0 do FreeMem(c^[j], (j+1)*SizeOf(ArbFloat));
  833. FreeMem(c, (n+1)*SizeOf(pointer));
  834. s:=0; for j:=0 to n do s:=s+b^[j].i*gamma[j]; w.i:=s; gamma[n+1] := -s;
  835. s:=0; for j:=0 to n do s:=s+b^[j].j*gamma[j]; w.j:=s; gamma[n+2] := -s;
  836. s:=0; for j:=0 to n do s:=s+b^[j].k*gamma[j]; w.k:=s; gamma[n+3] := -s;
  837. FreeMem(b, (n+1)*SizeOf(vector));
  838. u.init(w.i/dv.i, w.j/dv.j, w.k/dv.k);
  839. u.scale(1/lambda);
  840. u.add(alfa);
  841. s:=0; for j:=0 to n do s:=s+e21^[j].i*gamma[j]; v.i := e22[0].inprod(w)-s;
  842. s:=0; for j:=0 to n do s:=s+e21^[j].j*gamma[j]; v.j := e22[1].inprod(w)-s;
  843. s:=0; for j:=0 to n do s:=s+e21^[j].k*gamma[j]; v.k := e22[2].inprod(w)-s;
  844. FreeMem(e21, (n+1)*SizeOf(vector));
  845. u.add(v);
  846. pfxpfy(an2, an1, an, u, pf); u1:=u.i;
  847. kwsom := 0; for j:=0 to n do kwsom:=kwsom+sqr(gamma[j])/alfad[j].imag;
  848. kwsom := kwsom+sqr(w.i)/dv.i+sqr(w.j)/dv.j+sqr(w.k)/dv.k;
  849. kwsom := kwsom/sqr(lambda);
  850. s:=0;
  851. for i:=0 to n do
  852. begin s1:=0;
  853. for j:=0 to i do s1:=s1+k^[i]^[j]*gamma[j];
  854. for j:=i+1 to n do s1:=s1+k^[j]^[i]*gamma[j];
  855. s := gamma[i]*s1+s
  856. end;
  857. for j:=n downto 0 do FreeMem(k^[j], (j+1)*SizeOf(ArbFloat));
  858. FreeMem(k, (n+1)*SizeOf(pointer));
  859. energie := s
  860. end {splineparameters};
  861. end.