Following documentation and some math calculations I have replicated shape of
sor as spline/function. For given height (value along y axis) it returns
radius. Like for other splines heights below first valid point return first
valid radius and above last valid point return last valid radius. Call below
macro with the same set of points as in sor object. NOTE: Macro/function not
validate points! If you want array validation just make sor{} before this
macro and parser will tell you if points are incorrect. I have used the same
symbols as in 6.5.1.12 chapter so reading/investigating source should be easy.
Of course 3.5 compatible. Any comments ?

#macro sor_spline(Array)
  #local N=dimension_size(Array,1);
  #local a=array[N][2];
  #local n=N;
  #while (n)
    #local n=n-1;
    #local a[n][0]=Array[n].u;
    #local a[n][1]=Array[n].v;
  #end
  function(h){
    select(
      h-a[1][1],
      a[1][0],
      #local n=2;
      #while (n<N-1)
        #local b0=a[n-1][0]^2;
        #local b1=a[n][0]^2;
        #local b2=2*a[n-1][0]*(a[n][0]-a[n-2][0])/
                  (a[n][1]-a[n-2][1]);
        #local b3=2*a[n][0]*(a[n+1][0]-a[n-1][0])/
                  (a[n+1][1]-a[n-1][1]);
        #local M00=a[n-1][1]^3;
        #local M01=a[n-1][1]^2;
        #local M02=a[n-1][1];
        #local M10=a[n][1]^3;
        #local M11=a[n][1]^2;
        #local M12=a[n][1];
        #local M20=3*a[n-1][1]^2;
        #local M21=2*a[n-1][1];
        #local M30=3*a[n][1]^2;
        #local M31=2*a[n][1];
        #local M2131=M21-M31;
        #local M1101=M11-M01;
        #local M1202=M12-M02;
        #local M1101_1=M1101-M31*M1202;
        #local M1101_2=M21*M1202-M1101;
        #local A=((b0-b1)*M2131+b2*M1101_1+b3*M1101_2)/
                 ((M00-M10)*M2131+M20*M1101_1+M30*M1101_2);
        #local B=(b2-b3+A*(M30-M20))/M2131;
        #local C=b3-M30*A-M31*B;
        #local D=b1-M10*A-M11*B-M12*C;
        select(
          h-a[n][1],
          (A*h^3+B*h^2+C*h+D)^.5,
        #local n=n+1;
      #end
      a[N-2][0]
      #local n=2;#while(n<N-1))#local n=n+1;#end
    )
  }
#end

ABX

Post a reply to this message