DataMuseum.dk

top - metrics - download

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

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

TextFile

message spln3int

spln3int=algol message.no
cubic spline integral
external
real procedure spln3int(x,X,Y,M,n);
value n,x; integer n; real x; array X,Y,M;
begin
  own integer i;
  real xm,xp,h;
  own real int;
  if n=0 then 
  begin
    i:=0; int:=0 
  end
  else
  begin
    if x<X(1) or x>X(n) then alarm(<:<10>***spln3int illegal x = :>,
    string exactlay(x,i,x),x);
    if i>=n then 
    begin
      i:=0; int:=0 
    end;
    for xp:=X(i+1)-x while xp<=0 and i<n-1 do
    begin
      h:=(xp+x-X(i))/2;
      int:=int-h*(h*h*(M(i+1)+M(i))/3-Y(i+1)-Y(i));
      i:=i+1
    end;
    for xm:=x-X(i) while xm<0 do
    begin
      i:=i-1;
      h:=(x-xm-X(i))/2;
      int:=(if i=1 then 0
      else int+h*(h*h*(M(i+1)+M(i))/3-Y(i+1)-Y(i)));
      xp:=-xm
    end;
    spln3int:=int+xm*(Y(i)+xm*((Y(i+1)-Y(i))/2/(xm+xp)
    -(M(i+1)+M(i)+M(i))/12*(xm+xp)+xm*(M(i)/6
    +xm*(M(i+1)-M(i))/24/(xm+xp))))
  end
end spln3int
; end
▶EOF◀