DataMuseum.dk

top - metrics - download

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

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

TextFile

message spln3inv

spln3inv=algol message.no
cubic spline inverse function
external
boolean procedure spln3inv(x,y,X,Y,M,n);
value y,n; integer n; real x,y; array X,Y,M;
begin
  real ym,xm,yp,xp,h,a0,a1,a2,a3;
  own real xh,yh,x0,y0;
  own integer i;
  spln3inv:=false;
  if n=0 then i:=0
  else 
  begin
    if n<0 then 
    begin
      n:=-n; i:=0 
    end; 
    if x<X(1) or x>X(n) then alarm(<:<10>***spln3inv illegal x = :>,
    string exactlay(x,i,x),x);
    xp:=1;
    if i>=n then i:=0;
    if x<>x0 or y<>y0 or i=0 then xh:=x
    else
    if x=xh then xp:=-1;
    for i:=i while xp>0 do 
    begin
      for xp:=X(i+1)-xh while xp<=0 and i<n-1 do i:=i+1;
      for xm:=xh-X(i) while xm<0 and i>1 do 
      begin
        xp:=-xm; i:=i-1
      end;
      h:=xm+xp;
      a0:=Y(i);
      a2:=M(i)/2;
      a1:=(Y(i+1)-a0)/h-(4*a2+M(i+1))*h/6;
      a3:=(M(i+1)-M(i))/h/6;
      ym:=(if x<>xh then yh else
      if x=x0 and i>0 then y0
      else ((a3*xm+a2)*xm+a1)*xm+a0);
      if a3=0 then xh:=(if a2<>0 then -a1/a2/2 else h)
      else 
      begin
        xh:=(a2*a2/a3/3-a1)/a3/3;
        if xh>0 then 
        begin
          xh:=-a2/3/a3-sqrt(xh);
          if xh<=xm then xh:=-xh-2*a2/3/a3 
        end
        else
        if xh=0 then xh:=-a2/3/a3
      end;
      if xh<=xm or xh>=h then 
      begin
        xh:=h; yh:=Y(i+1) 
      end 
      else
      yh:=((a3*xh+a2)*xh+a1)*xh+a0;

      if y=ym or (ym<y and y<yh) or (ym>y and y>yh)
      or (y=yh and xp=0) then 
      begin
        if xp=0 then x0:=xh
        else
        if y=ym then x0:=xm
        else 
        begin
          if y<ym then
          begin
            xp:=xm; xm:=xh; yp:=ym; ym:=yh;
          end
          else 
          begin
            xp:=xh; yp:=yh 
          end;
          for x0:=xm+(xp-xm)/(yp-ym)*(y-ym) while x0<>h do
          begin
            h:=x0;
            y0:=((a3*x0+a2)*x0+a1)*x0+a0;
            if y0<y then 
            begin
              xm:=x0; ym:=y0 
            end
            else
            if y<y0 then 
            begin
              xp:=x0; yp:=y0 
            end
          end
        end;
        x:=x0:=x0+X(i); y0:=y;
        xp:=-1; spln3inv:=true 
      end;
      xh:=xh+X(i);
    end for i; 
  end i<>0; 
end spln3inv
; end
▶EOF◀