DataMuseum.dk

top - metrics - download

    Length: 2304 (0x900)
    Types: TextFile
    Names: »tpolfit1«

└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
    └─⟦b2ec5d50f⟧ 
        └─⟦09b4e9619⟧ »thcømat« 
            └─⟦this⟧ 

TextFile

; 28 04 80 Ole J. Heilmann
polfit1=algol
external
real procedure polfit1( X, R, n, C, S, r);
value n, r; integer n, r; array X, R, C, S;

if n<=0 or r<0  then begin
   for r:=r step -1 until 0 do begin
      C(r):=0; S(r):=maxreal end;
   polfit1:=-1
end else
begin integer j, k , l;
real f, fj, fx, f1, rf, a, b, c, Q, x, y, xn, xk;
array F, F1(1:n), A, B, D(0:r);
comment initialisering;
if r>=n then begin
   for j:=r step -1 until n do begin
      C(j):=0; S(j):=maxreal end;
   r:=n-1 end;
fx:= rf:= 0;
for j:= 1 step 1 until n do
begin
      F(j):=1; F1(j):= 0; fx:=fx+X(j);  rf:= rf + R(j);
end j;
f:=n; f1:=1;
comment beregn polynomiet;
for k:= 0 step 1 until r do
begin a:=A(k):=fx/f;
      b:=B(k):=f/f1;
      c:=C(k):=rf/f;
      S(k):=1/f;
      f1:= f; f:= fx:= rf:= 0;
      for j:= 1 step 1 until n do
      begin fj:= (X(j) - a) * F(j) - b * F1(j);
            R(j):= R(j) - F(j) * c;
            F1(j):= F(j); F(j):= fj;
            rf:= rf + R(j) * fj;
            f:= f + fj * fj; fx:= fx + X(j) * fj * fj
      end j
end beregn;
comment omform til sædvanligt polynomium og beregn usikkerheder;
if n-r>1 then begin
   Q:=0;
   for j:=1 step 1 until n do Q:=R(j)**2+Q;
   Q:=polfit1:=Q/(n-r-1) end
else begin
   Q:=1; polfit1:=0 end;
for j:=0 step 1 until r do D(j):=0;
for k:=0 step 1 until r do begin
   c:=D(k); D(k):=a:=1; b:=0; fj:=C(k); rf:=S(k);
   for j:=k+1 step 1 until r do begin
      f:=-A(j-1)*a-B(j-1)*b+c;
      c:=D(j); b:=a; D(j):=a:=f;
      fj:=C(j)*a+fj; rf:=S(j)*a*a+rf end;
   C(k):=fj; S(k):=sqrt(rf*Q) end
end polfit1
; end
▶EOF◀