DataMuseum.dk

top - metrics - download

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

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

TextFile

; 12 09 73 Anders Lindgård
; 28 04 80 Ole J. Heilmann
\f


message polfit1

polfit=algol message.no
external
real procedure polfit( 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;
  polfit:=-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:=polfit:=Q/(n-r-1) 
  end
  else 
  begin
    Q:=1; polfit:=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 polfit
; end
▶EOF◀