DataMuseum.dk

top - metrics - download

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

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

TextFile

message spln3smth

spln3smth=algol message.no
cubic splinesmoothing
external
boolean procedure spln3smth(X,YS,SD,n);
value n; integer n;
array X,YS,SD;
begin
  integer i,j;
  real s,s0,sc,w,m,p,q,dx,dx1,a,b;
  array A(1:n,1:n),Q(1:n-1);
  integer array P(1:n);
  spln3smth:=false;
  if n<4 then goto END;
  for i:=2 step 1 until n do
  begin
    if X(i-1)>=X(i) then goto END
  end;
  p:=q:=dx:=dx1:=0;
  for i:=1 step 1 until n do
  begin
    p:=p+YS(i); q:=q+X(i)*YS(i);
    dx:=dx+X(i); dx1:=dx1+X(i)**2
  end;
  a:=(n*q-dx*p)/(n*dx1-dx*dx);
  b:=(dx1*p-dx*q)/(n*dx1-dx*dx);
  YS(1):=YS(1)-a*X(1)-b;
  s:=s0:=YS(1);
  for i:=2 step 1 until n do
  begin
    YS(i):=sc:=YS(i)-a*X(i)-b;
    if sc<s0 then s0:=sc
    else
    if sc>s then s:=sc
  end;
  if s=s0 then goto END;
  sc:=(s-s0)**2/(X(n)-X(1))**3;
  for i:=1 step 1 until n do
  begin
    dx:=X(2)-X(1);
    s:=(if i=1 then -1 else if i=2 then 1 else 0)/dx;
    q:=m:=Q(1):=A(i,1):=0;
    for j:=2 step 1 until n-1 do
    begin
      dx1:=dx; dx:=X(j+1)-X(j);
      s0:=s; s:=(if i=j then -1 else if i=j+1 then 1 else 0)/dx;
      p:=2*(dx+dx1)-dx1*q;
      A(i,j):=m:=(6*(s-s0)-dx1*m)/p;
      Q(j):=q:=dx/p
    end j;
    w:=SD(i)**2/sc;
    A(i,n):=s0:=m/dx*w;
    dx1:=X(n-1);
    for j:=n-2 step -1 until 1 do 
    begin
      q:=m; m:=A(i,j)-Q(j)*q;
      dx:=dx1; dx1:=X(j);
      s:=s0; s0:=(m-q)/(dx-dx1)*w;
      A(i,j+1):=s0-s
    end j;
    A(i,1):=-s0;
    A(i,i):=A(i,i)+1;
  end i;
  if decompose(A,P,1) then
  begin
    solve(A,P,1,YS);
    for i:=1 step 1 until n do YS(i):=YS(i)+a*X(i)+b;
    spln3smth:=true
  end;
END:
end spln3smth
; end
▶EOF◀