DataMuseum.dk

top - metrics - download

    Length: 1536 (0x600)
    Types: TextFile
    Names: »tspln3«

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

TextFile

message spln3

spln3=algol message.no
cubic spline initialization
external
procedure spln3(type,X,Y,M,n);
value n; integer n;
array X,Y,M,type;
begin
  if abs(n)<2 then alarm(<:<10>***spln3 n<60>2:>);
  if n>1 then
  begin
    spln3int(.0,X,Y,M,0);
    spln3inv(.0,.0,X,Y,M,0)
  end
  else n:=-n;
  begin
    integer k;
    real dx,dx1,s,s1,p,q;
    array Q(1:n);
    for k:=2 step 1 until n do
    begin
      if X(k)<=X(k-1) then
      alarm(<:<10>***spln3 illegal X(i):>)
    end;
    dx:=X(2)-X(1);
    s:=(Y(2)-Y(1))/dx;
    case round type(0)+1 of
    begin
      begin
        Q(1):=q:=type(2);
        M(1):=type(1);
        M(n):=type(4);
        Q(n):=type(3)
      end;
      Q(1):=Q(n):=M(1):=M(n):=q:=0.0;
      begin
        Q(1):=q:=Q(n):=.5;
        M(1):=3/dx*(s-type(1));
        M(n):=3/(X(n)-X(n-1))*
        (type(4)-(Y(n)-Y(n-1))/(X(n)-X(n-1)))
      end;
      begin
        M(1):=M(n):=0;
        Q(1):=q:=Q(n):=-3/4
      end;
    end;
    for k:=2 step 1 until n-1 do
    begin
      dx1:=dx; dx:=X(k+1)-X(k);
      s1:=s; s:=(Y(k+1)-Y(k))/dx;
      p:=2*(dx+dx1)-dx1*q;
      M(k):=(6*(s-s1)-dx1*M(k-1))/p;
      Q(k):=q:=dx/p
    end k;
    M(n):=(M(n)-Q(n)*M(n-1))/(1-Q(n)*Q(n-1));
    for k:=n-1 step -1 until 1 do M(k):=M(k)-Q(k)*M(k+1);
  end 
end spln3
; end
▶EOF◀