/*
Include file: curve.inc
Spline and curve interpolation for many variants.

Based upon the works of Steve H. Noskowicz !

http://web.archive.org/web/20151002232205/http://home.comcast.net/~k9dci/site/?/page/Piecewise_Polynomial_Interpolation/
Get the book PDF and the appendix.

Author : Ingo Janssen
Date   : 2019-06-16

Version: 0.2 Alpha, 2021-05-24
         complete re-write from earlier attempt.
         cubic Lagrange still not right.

Run the file without including to render the tests at the end of the file.
These also show how to use the macro's. 

Curve Interpolation macros:
  CIlinear(VectorArray)
  CIsimple3(VectorArray)
  CIbSpline2(VectorArray)
  CIbezier2(VectorArray)
  CIgolden(VectorArray)
  CIbezier3(VectorArray)
  CI4p3(VectorArray)
  CIhermite(VectorArray)
  CIbSpline3(VectorArray)
  CIcatmul(VectorArray)
  CIsimple3ESC(VectorArray, EndSlope)
  CIsimple3MSC(VectorArray, MidSlope)
  CIsimple3IESC(VectorArray, EndSlope, EndSlope)
  CIbezier2EVC(VectorArray, Velocity)
  CIlagrange2(VectorArray, T2)
  CIlagrange3(VectorArray, T1, T2) !!broken!!
  CIbeta(VectorArray, Tension, Bias)
  CItcb(VectorArray, Tension, Continuity, Bias)
  CIbezier3EVC(VectorArray, Velocity)
  CIpalmer(VectorArray, Tension, Bias)

#declare BzCurve = CIbezier3(array[n]{....});
#declare PointOnBz = CIget(BzCurve, 0.5);

*/

#version 3.7;
#ifndef(CURVE_INC_TEMP)
  #declare CURVE_INC_TEMP = version;
  #ifdef(View_POV_Include_Stack)
    #debug "including curve.inc\n"
  #end


#macro arrFill(Len, Val)
  /*:Create and fill an array of length Len with value Val*/
  #local Arr = array[Len];
  #for(i,0,Len-1)
    #local Arr[i] = Val;
  #end
  Arr
#end 


//enum CICurve items
#declare _curve_ctrl = 0;
#declare _curve_point = 1;
#declare _curve_curmult = 2;
#declare _curve_def = 3;
//enum curve ctrl items
#declare _curve_type = 0;
#declare _curve_len = 1;
#declare _curve_segments = 2;
#declare _curve_curseg = 3;


//enumerate primary array  curve type definitions
#declare _ci_linear = 0;
#declare _ci_simple3 = 1;
#declare _ci_b_spline2 = 2;
#declare _ci_bezier2 = 3;
#declare _ci_golden = 4;
#declare _ci_bezier3 = 5;
#declare _ci_4p3 = 6;
#declare _ci_hermite= 7;
#declare _ci_b_spline3 = 8;
#declare _ci_catmul = 9;
#declare _ci_simple3_esc = 10;
#declare _ci_simple3_msc = 11;
#declare _ci_simple3_iesc = 12;
#declare _ci_bezier2_evc = 13;
#declare _ci_lagrange2 = 14;
#declare _ci_lagrange3 = 15;
#declare _ci_beta = 16;
#declare _ci_tcb = 17;
#declare _ci_bezier3_evc = 18;
#declare _ci_palmer = 19;

//enumerate secondary curve type definitions
#declare _ci_ctrl = 0;
#declare _ci_var = 1;
#declare _ci_matrix = 2; 
//enumerate variables in _ci_ctrl
#declare _ci_dimseg = 0;
#declare _ci_step = 1;
#declare _ci_end = 2;
//enumerate variables in _ci_var
#declare _ci_div = 0;


#declare _CIdef = array;
#declare _CIdef[_ci_linear] = array[3]{
  array[3]{2, 1, 1},  //dimseg, step, end
  array[1]{1},        //div
  array[2]{
    array[2]{-1, 1},
    array[1]{ 1}
  }
};
#macro CIlinear(CIarr)
  _CIinit(CIarr, _ci_linear, _CIdef[_ci_linear])
#end


#declare _CIdef[_ci_simple3] = array[3]{
  array[3]{2, 1, 1},  //dimseg, step, end
  array[1]{1},        //div              
  array[4]{
    array[2]{ 2,-2},
    array[2]{-3, 3},
    array[2]{ 0, 0},
    array[1]{ 1}
  }
};
#macro CIsimple3(CIarr)
  _CIinit(CIarr, _ci_simple3, _CIdef[_ci_simple3])
#end


#declare _CIdef[_ci_b_spline2] = array[3]{
  array[3]{3, 1, 1},  //dimseg, step, end
  array[1]{2},         //div              
  array[3]{
    array[3]{ 1,-2, 1},
    array[2]{-2, 2},
    array[2]{ 1, 1}
  }
};
#macro CIbSpline2(CIarr)
  _CIinit(CIarr, _ci_b_spline2, _CIdef[_ci_b_spline2])
#end


#declare _CIdef[_ci_bezier2] = array[3]{
  array[3]{3, 2, 2},  //dimseg, step, end
  array[1]{1},        //div              
  array[3]{
    array[3]{ 1,-2, 1},
    array[2]{-2, 2},
    array[1]{ 1}
  }
};
#macro CIbezier2(CIarr)
  _CIinit(CIarr, _ci_bezier2, _CIdef[_ci_bezier2])
#end


#declare _CIdef[_ci_golden] = array[3]{
  array[3]{3, 2, 2},  //dimseg, step, end
  array[1]{1},         //div              
  array[3]{
    array[3]{ 2,-4, 2},
    array[3]{-3, 4,-1},
    array[3]{ 1, 0, 0}
  }
};
#macro CIgolden(CIarr)
  _CIinit(CIarr, _ci_golden, _CIdef[_ci_golden])
#end


#declare _CIdef[_ci_bezier3] = array[3]{
  array[3]{4, 3, 3},  //dimseg, step, end
  array[1]{1},        //div              
  array[4]{
    array[4]{-1, 3,-3, 1},
    array[3]{ 3,-6, 3},
    array[2]{-3, 3,},
    array[1]{ 1}
  }
};
#macro CIbezier3(CIarr)
  _CIinit(CIarr, _ci_bezier3, _CIdef[_ci_bezier3])
#end


#declare _CIdef[_ci_4p3] = array[3]{
  array[3]{4, 3, 3},  //dimseg, step, end
  array[1]{2},         //div              
  array[4]{
    array[4]{ -9, 27,-27, 9},
    array[4]{ 18,-45, 36,-9},
    array[4]{-11, 18, -9, 2},
    array[4]{  2,  0,  0, 0}
  }
};
#macro CI4p3(CIarr)
  _CIinit(CIarr, _ci_4p3, _CIdef[_ci_4p3])
#end


#declare _CIdef[_ci_hermite] = array[3]{
  array[3]{4, 3, 3},  //dimseg, step, end
  array[1]{1},         //div              
  array[4]{
    array[4]{ 2, 1, 1,-2},
    array[4]{-3,-2,-1, 3},
    array[4]{ 0, 1, 0, 0},
    array[4]{ 1, 0, 0, 0}
  }
};
#macro CIhermite(CIarr)
  _CIinit(CIarr, _ci_hermite, _CIdef[_ci_hermite])
#end


#declare _CIdef[_ci_b_spline3] = array[3]{
  array[3]{4, 1, 2},  //dimseg, step, end
  array[1]{6},         //div              
  array[4]{
    array[4]{-1, 3,-3, 1},
    array[3]{ 3,-6, 3},
    array[3]{-3, 0, 3},
    array[3]{ 1, 4, 1}
  }
};
#macro CIbSpline3(CIarr)
  _CIinit(CIarr, _ci_b_spline3, _CIdef[_ci_b_spline3])
#end


#declare _CIdef[_ci_catmul] = array[3]{
  array[3]{4, 1, 2},  //dimseg, step, end
  array[1]{2},         //div              
  array[4]{
    array[4]{-1, 3,-3, 1},
    array[4]{ 2,-5, 4,-1},
    array[4]{-1, 0, 1, 0},
    array[4]{ 0, 2, 0, 0}
  }
};
#macro CIcatmul(CIarr)
  _CIinit(CIarr, _ci_catmul, _CIdef[_ci_catmul])
#end


#declare _CIdef[_ci_simple3_esc] = array[3]{
  array[3]{2, 1, 1},  //dimseg, step, end
  array[1]{1},        //div              
  array[4]{
    array[2]{
      function(ESC){2-(2*ESC)}, 
      function(ESC){(2*ESC)-2}
    },
    array[2]{
      function(ESC){(3*ESC)-3}, 
      function(ESC){3-(3*ESC)}
    },
    array[2]{
      function(ESC){-ESC},
      function(ESC){ESC}
    },
    array[1]{
      function(ESC){1}
    }
  }
};
#macro CIsimple3ESC(CIarr, ESC)
  #local Type = _ci_simple3_esc;  
  #local Def = _CIoneVar(Type, ESC)
  _CIinit(CIarr, Type, Def)
#end


#declare _CIdef[_ci_simple3_msc] = array[3]{
  array[3]{2, 1, 1},  //dimseg, step, end
  array[1]{1},        //div              
  array[4]{
    array[2]{
      function(MSC){(4*MSC)-4}, 
      function(MSC){4-(4*MSC)}
    },
    array[2]{
      function(MSC){6-(6*MSC)}, 
      function(MSC){(6*MSC)-6}
    },
    array[2]{
      function(MSC){2*MSC-3},   
      function(MSC){3-(2*MSC)}
    },
    array[1]{
      function(MSC){1}
    }
  }
};
#macro CIsimple3MSC(CIarr, MSC)
  #local Type = _ci_simple3_msc;
  #local Def = _CIoneVar(Type, MSC)
  _CIinit(CIarr, Type, Def)
#end


#declare _CIdef[_ci_simple3_iesc] = array[3]{
  array[3]{2, 1, 1},  //dimseg, step, end
  array[1]{1},        //div              
  array[4]{
    array[2]{
      function(ESC0, ESC1){2 - ESC0 - ESC1},
      function(ESC0, ESC1){ESC0 + ESC1 - 2}
    },
    array[2]{
      function(ESC0, ESC1){(2*ESC0) + ESC1 - 3}, 
      function(ESC0, ESC1){3-(2 * ESC0) - ESC1}
    },
    array[2]{
      function(ESC0, ESC1){-ESC0},               
      function(ESC0, ESC1){ESC0}
    },
    array[1]{
      function(ESC0, ESC1){1}
    }
  }
};
#macro CIsimple3IESC(CIarr, ESC0, ESC0)
  #local Type = _ci_simple3_iesc;
  #local Def = _CItwoVar(Type, ESC0, ESC0)
  _CIinit(CIarr, Type, Def)
#end


#declare _CIdef[_ci_bezier2_evc] = array[3]{
  array[3]{3, 2, 2},  //dimseg, step, end
  array[1]{1},        //div              
  array[4]{
    array[3]{
      function(V){-V},          
      function(V){0},    
      function(V){V}
    },
    array[3]{
      function(V){1 + (2 * V)}, 
      function(V){-2-V}, 
      function(V){1-V}
    },
    array[2]{
      function(V){-2-V},        
      function(V){2+V}
    },
    array[1]{
      function(V){1}
    }
  }
};
#macro CIbezier2EVC(CIarr, V)
  #local Type = _ci_bezier2_evc;
  #local Def = _CIoneVar(Type, V)
  _CIinit(CIarr, Type, Def)
#end


#declare _CIdef[_ci_lagrange2] = array[3]{
  array[3]{3, 2, 2},   //dimseg, step, end
  array[1]{1},         //div              
  array[3]{
    array[3]{
      function(N1,N2,N3,T2){N1},    
      function(N1,N2,N3,T2){N2},  
      function(N1,N2,N3,T2){N3}
    },
    array[3]{
      function(N1,N2,N3,T2){-N1-1}, 
      function(N1,N2,N3,T2){-N2}, 
      function(N1,N2,N3,T2){-N3*T2}
    },
    array[2]{
      function(N1,N2,N3,T2){1},      
      function(N1,N2,N3,T2){0}
    }
  }
};
#macro CIlagrange2(CIarr, T2)
  /* 0<T2<1 but not really, negative values work, values bigger than 1 work. Experiment.*/
  #local Type = _ci_lagrange2;
  #local N1 = 1/T2;
  #local N2 = 1/(T2*(T2-1));
  #local N3 = 1/(1-T2);
  #local Def = _CIfourVar(Type,N1,N2,N3,T2);
  _CIinit(CIarr, Type, Def)
#end

 
///!!!!!!!!!NOT RIGHT!!!!!!!!!!!! 
#declare _CIdef[_ci_lagrange3] = array[3]{
  array[3]{4, 3, 3},   //dimseg, step, end
  array[1]{1},         //div              
  array[4]{
    array[4]{
      function(N0,N1,N2,N3,T1,T2){-N0},
      function(N0,N1,N2,N3,T1,T2){ N1},
      function(N0,N1,N2,N3,T1,T2){ N2},
      function(N0,N1,N2,N3,T1,T2){ N3}
    },
    array[4]{
      function(N0,N1,N2,N3,T1,T2){ N0*(T1+T2+1)},
      function(N0,N1,N2,N3,T1,T2){-N1*(T2+1)},
      function(N0,N1,N2,N3,T1,T2){-N2*(T1+1)},
      function(N0,N1,N2,N3,T1,T2){-N3*(T1+T2)}
    },
    array[4]{
      function(N0,N1,N2,N3,T1,T2){-N0*(T1+T2+(T1*T2))},
      function(N0,N1,N2,N3,T1,T2){ N1*T2},
      function(N0,N1,N2,N3,T1,T2){ N2*T1},
      function(N0,N1,N2,N3,T1,T2){ N3*T1*T2}      
    },
    array[1]{function(N0,N1,N2,N3,T1,T2){1}}
  }
};
#macro CIlagrange3(CIarr, T1, T2)
  #local Type = _ci_lagrange3;
  #local N0 = 1/(T1 * T2);
  #local N1 = 1/(T1 * (T1 - T2) * (T1 - 1));
  #local N2 = 1/(T2 * (T2 - T1) * (T2 - 1));
  #local N3 = 1/(1 - T1) * (1 - T2);
  #local Def = _CIsixVar(Type,N0,N1,N2,N3,T1,T2);
  _CIinit(CIarr, Type, Def)
#end
///!!!!!!!!!NOT RIGHT!!!!!!!!!!!! 


#declare _CIdef[_ci_beta] = array[3]{
  array[3]{4, 1, 2},  //dimseg, step, end
  array[1]{function(T,B){(T + (2 * pow(B, 3)) + (4*pow(B,2)) + (4*B) + 2)}}, //div              
  array[4]{
    array[4]{
      function(T,B){-2 * pow(B, 3)},
      function(T,B){ 2 * (T + pow(B, 3) + pow(B, 2) + B)},
      function(T,B){-2 * (T + pow(B, 2) + B + 1)},
      function(T,B){ 2}
    },
    array[3]{
      function(T,B){ 6 * pow(B, 3)},
      function(T,B){-3 * (T + (2 * pow(B, 3) + (2 * pow(B, 2))))},
      function(T,B){ 3 * (T + (2 * pow(B, 2)))},
    },
    array[3]{
      function(T,B){-6 * pow(B, 3)},
      function(T,B){ 6 * (pow(B, 3) - B)},
      function(T,B){ 6 * B},
    },
    array[3]{
      function(T,B){2 * pow(B, 3)},
      function(T,B){T + (4 * (pow(B, 2) + B))},
      function(T,B){2},      
    }
  }
};
#macro CIbeta(CIarr, T, B)
  #local Type = _ci_beta;
  #local Def  = _CItwoVarDiv(Type, T, B);
  _CIinit(CIarr, Type, Def)
#end


#declare _CIdef[_ci_tcb] = array[3]{
  array[3]{4, 1, 2},  //dimseg, step, end
  array[1]{2},        //div              
  array[4]{
    array[4]{
      function(A,B,C,D){-A},
      function(A,B,C,D){ 4 + A - B - C},
      function(A,B,C,D){-4 + B + C - D},
      function(A,B,C,D){D}
    },
    array[4]{
      function(A,B,C,D){ 2 * A},
      function(A,B,C,D){-6 - (2 * A) + (2 * B) + C},
      function(A,B,C,D){ 6 - (2 * B) - C + D},
      function(A,B,C,D){-D}    
    },
    array[3]{
      function(A,B,C,D){-A},
      function(A,B,C,D){A - B},
      function(A,B,C,D){B}    
    },
    array[2]{
      function(A,B,C,D){0},
      function(A,B,C,D){2},
    }
  }
};
#macro CItcb(CIarr, Te, Co, Bi)
  #local Type = _ci_tcb;
  #local A = (1 - Te) * (1 + Co) * (1 + Bi);
  #local B = (1 - Te) * (1 - Co) * (1 - Bi);
  #local C = (1 - Te) * (1 - Co) * (1 + Bi);
  #local D = (1 - Te) * (1 + Co) * (1 - Bi);  
  #local Def  = _CIfourVar(Type, A, B, C, D);
  _CIinit(CIarr, Type, Def)
#end


#declare _CIdef[_ci_bezier3_evc] = array[3]{
  array[3]{4, 3, 3},  //dimseg, step, end
  array[1]{1},         //div              
  array[4]{
    array[4]{
      function(V){-1 - V},   
      function(V){ 3 + V},   
      function(V){-3 - V},   
      function(V){ 1 + V}   
    },
    array[4]{ 
      function(V){ 3 + (2 * V)},   
      function(V){-6 - (2 * V)},   
      function(V){ 3 + V},   
      function(V){-V}   
    },
    array[2]{ 
      function(V){-3 - V},   
      function(V){ 3 + V}   
    },
    array[1]{ 
      function(V){1}   
    }
  }
};
#macro CIbezier3EVC(CIarr, V)
  #local Type = _ci_bezier3_evc;
  #local Def  = _CIoneVar(Type, V);
  _CIinit(CIarr, Type, Def)
#end


#declare _CIdef[_ci_palmer] = array[3]{
  array[3]{3, 2, 2},  //dimseg, step, end
  array[1]{1},         //div              
  array[4]{
    array[3]{
      function(N0,N1,N2,T){-2 * T},
      function(N0,N1,N2,T){ 0},
      function(N0,N1,N2,T){ 2 * T},
    },
    array[3]{
      function(N0,N1,N2,T){ N0 + (3 * T)},
      function(N0,N1,N2,T){-N1},
      function(N0,N1,N2,T){ N2 - (3 * T)},
    },
    array[3]{
      function(N0,N1,N2,T){-1 - N0 - T},
      function(N0,N1,N2,T){ N1},
      function(N0,N1,N2,T){ T - (N2 / N0)},
    },
    array[1]{
      function(N0,N1,N2,T){1},
    }
  }
};
#macro CIpalmer(CIarr, T, B)
  #local Type = _ci_palmer;
  #local N0 = 1 / B;
  #local N1 = 1 / (B * (1 - B));
  #local N2 = 1 / (1 - B);
  #local Def = _CIfourVar(Type, N0, N1, N2, T);
  _CIinit(CIarr, Type, Def)
#end


/*
//template
#declare _CIdef[xxxx] = array[3]{
  array[3]{0, 0, 0},  //dimseg, step, end
  array[1]{1},         //div              
  array[4]{
    array[4]{ 0, 0, 0, 0},
    array[4]{ 0, 0, 0, 0},
    array[4]{ 0, 0, 0, 0},
    array[4]{ 0, 0, 0, 0}
  }
};
#macro CIxxxx(CIarr)
  _CIinit(CIarr, xxxx, _CIdef[xxx])
#end
*/


#macro _CIdefDef(Type)
  #local Def = array[3];
  #local Def[_ci_ctrl] = _CIdef[Type][_ci_ctrl];
  #local Def[_ci_var] = _CIdef[Type][_ci_var];
  #local Def[_ci_matrix] = array[dimension_size(_CIdef[Type][_ci_matrix], 1)];  
  Def
#end


#macro _CIdefDefDiv(Type, Var0, Var1)
  #local Def = array[3];
  #local Def[_ci_ctrl] = _CIdef[Type][_ci_ctrl];
  #local Def[_ci_var] = array[1];
  #local Def[_ci_var][_ci_div] = _CIdef[Type][_ci_var][_ci_div](Var0, Var1);
  #local Def[_ci_matrix] = array[dimension_size(_CIdef[Type][_ci_matrix], 1)];  
  Def
#end


#macro _CIoneVar(Type, Var)
  #local Def = _CIdefDef(Type);
  #for(i, 0, dimension_size(_CIdef[Type][_ci_matrix], 1) - 1)
    #local Def[_ci_matrix][i] = array[dimension_size(_CIdef[Type][_ci_matrix][i], 1)];
    #for(j, 0, dimension_size(Def[_ci_matrix][i], 1) - 1)
      #local Def[_ci_matrix][i][j] = _CIdef[Type][_ci_matrix][i][j](Var);
    #end
  #end
  Def
#end


#macro _CItwoVar(Type, Var0, Var1)
  #local Def = _CIdefDef(Type);
  #for(i, 0, dimension_size(_CIdef[Type][_ci_matrix], 1) - 1)
    #local Def[_ci_matrix][i] = array[dimension_size(_CIdef[Type][_ci_matrix][i], 1)];
    #for(j, 0, dimension_size(Def[_ci_matrix][i], 1) - 1)
      #local Def[_ci_matrix][i][j] = _CIdef[Type][_ci_matrix][i][j](Var0, Var1);
    #end
  #end
  Def
#end


#macro _CItwoVarDiv(Type, Var0, Var1)
  #local Def = _CIdefDefDiv(Type, Var0, Var1);
  #for(i, 0, dimension_size(_CIdef[Type][_ci_matrix], 1) - 1)
    #local Def[_ci_matrix][i] = array[dimension_size(_CIdef[Type][_ci_matrix][i], 1)];
    #for(j, 0, dimension_size(Def[_ci_matrix][i], 1) - 1)
      #local Def[_ci_matrix][i][j] = _CIdef[Type][_ci_matrix][i][j](Var0, Var1);
    #end
  #end
  Def
#end


#macro _CIfourVar(Type, Var0, Var1, Var2, Var3)
  #local Def = _CIdefDef(Type);
  #for(i, 0, dimension_size(_CIdef[Type][_ci_matrix], 1) - 1)
    #local Def[_ci_matrix][i] = array[dimension_size(_CIdef[Type][_ci_matrix][i], 1)];
    #for(j, 0, dimension_size(Def[_ci_matrix][i], 1) - 1)
      #local Def[_ci_matrix][i][j] = _CIdef[Type][_ci_matrix][i][j](Var0, Var1, Var2, Var3);
    #end
  #end
  Def
#end
 

#macro _CIsixVar(Type, Var0, Var1, Var2, Var3, Var4, Var5)
  #local Def = _CIdefDef(Type);
  #for(i, 0, dimension_size(_CIdef[Type][_ci_matrix], 1) - 1)
    #local Def[_ci_matrix][i] = array[dimension_size(_CIdef[Type][_ci_matrix][i], 1)];
    #for(j, 0, dimension_size(Def[_ci_matrix][i], 1) - 1)
      #local Def[_ci_matrix][i][j] = _CIdef[Type][_ci_matrix][i][j](Var0, Var1, Var2, Var3, Var4, Var5);
    #end
  #end
  Def
#end


#macro CIget(CICurve, T)
  /*: Gets the point on the curve at T
  
  parameter:
    CICurve (array, vector): an array of vectors and other parameters as created via one 
                              of the spline initialisation macros. 
    T (float, [0,1]): position at which to interpolate the curve.
  
  return: 
    (vector) returns a point on the curve at T.
  
  use example:
    #declare TheCurve = CIlinear(PointArray);
    #declare P = CIget(TheCurve, 0.1);    
  */
  
  /* T goes from 0 to 1 for the whole curve. From current T find the proper curve segment.
     Tseg also goes from 0 to 1. Once the proper segment is found, get the proper Tseg.
  */
  #local Type = CICurve[_curve_ctrl][_curve_type];
  #local Len = CICurve[_curve_ctrl][_curve_len];
  #local Segments = CICurve[_curve_ctrl][_curve_segments];
  #local CImult = CICurve[_curve_curmult];
  #if(T = 1)
    #local Segment = Segments-1;
    #local Tseg = 1;
  #else
    #local Segment = floor(T * Segments);
    #local Tseg =  (T * Segments) - Segment;    
  #end
  #if(Segment != CICurve[_curve_ctrl][_curve_curseg])
    #local CICurve[_curve_ctrl][_curve_curseg] = Segment;  
    #local CurrPos = Segment * CICurve[_curve_def][_ci_ctrl][_ci_step];
    #local CImatrix = CICurve[_curve_def][_ci_matrix];
    #local Div = CICurve[_curve_def][_ci_var][_ci_div];  
    #local MaxPos = Len - CICurve[_curve_def][_ci_ctrl][_ci_dimseg] + CICurve[_curve_def][_ci_ctrl][_ci_step];
    #if(CurrPos < MaxPos)
      #for(i, 0, dimension_size(CImatrix, 1) - 1)
        #local CImult[i] = <0, 0, 0>;
        #for(j, 0, dimension_size(CImatrix[i], 1) - 1)
          #local jSegment = j + CurrPos;
          #local CImult[i] = CImult[i] + (CImatrix[i][j] / Div) * CICurve[_curve_point][jSegment];    
        #end
      #end
    #end
    #local CICurve[_curve_curmult] = CImult; 
  #end
  #local F = CImult[0];
  #for(i, 1, dimension_size(CImult, 1) - 1)
    #local F = (F * Tseg) + CImult[i];
  #end
  F
#end


#macro _CIinit(CIarr, Type, Def)
  #local Len = dimension_size(CIarr, 1);
  #local Segments = (Len-1) / _CIdef[Type][_ci_ctrl][_ci_step];
  #if(Len < _CIdef[Type][_ci_ctrl][_ci_dimseg] | Segments - int(Segments) != 0)
    #error "Not the right amount of points in curve array"
  #end
  #local CICurve = array[4];
  #local CICurve[_curve_ctrl] = array{Type, Len, Segments, -1};
  #local CICurve[_curve_point] = CIarr; 
  #local CICurve[_curve_curmult] = arrFill(dimension_size(_CIdef[Type][_ci_matrix],1), <0,0,0>);
  #local CICurve[_curve_def] = Def;
  CICurve
#end 


#if(input_file_name="curve.inc")
  //tests
  global_settings{ assumed_gamma 1.0 }
  #default{ finish{ ambient 0.1 diffuse 0.9 emission 0}} 

  #declare Lin = CIlinear(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Simp3 = CIsimple3(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Bspl2 = CIbSpline2(
    array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Bez2 = CIbezier2(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Gold = CIgolden(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Bez3 = CIbezier3(
    array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Four3 = CI4p3(
    array[7]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>,<2,0,0>,<2,1,0>,<3,1,0>}
  );
  #declare Herm = CIhermite(
    array[4]{<0, 0, 0>, <0, 2, 0>, <1, 2, 0>, <1, 0, 0>}
  );
  #declare Bspl3 = CIbSpline3(
    array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Catmul = CIcatmul(
    array[6]{<-1, 0, 0>,<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>,<2,0,0>}
  );
  #declare ESC = CIsimple3ESC(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    3 //End point Slope control  parameter
  );
  #declare MSC = CIsimple3MSC(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    0.5 //Mid point Slope control  parameter
  );          
  #declare IESC = CIsimple3IESC(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    -1, 5  //Independent End point Slope control parameters
  );
  #declare Bez2EVC = CIbezier2EVC(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    -3 //End point Velocity Control parameter
  );
  #declare Lag2 = CIlagrange2(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    0.3 
  );
  #declare Lag3 = CIlagrange3( //!!!!Not Right!!!!!!!!!??
    array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>}
    1/3, 2/3
  );
  #declare Beta = CIbeta(
    array[6]{<0,0,0>, <0,1,0>, <1,1,0>, <1,0,0>, <2,0,0>, <2,1,0>}
    0, 1  //Tension, Bias
  )
  #declare TCB = CItcb(
    array[6]{<0,0,0>, <0,1,0>, <1,1,0>, <1,0,0>, <2,0,0>, <2,1,0>}
    0, 0, 0 //Tension, Continuity, Bias
  )
  #declare Bez3EVC = CIbezier3EVC(
    array[7]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>, <1,-1,0>, <2,-1,0>, <2,0,0>},
    2 //End point Velocity Control parameter
  );
  #declare Palmer = CIpalmer(
    array[5]{<0, 0, 0>, <0, 1, 0>, <1, 0.5, 0>, <2,1,0>,<2,0,0>},
    0, 0.5 //Tension, Bias 
  );

  
  // Set curve for test render //
  #declare Test =  M_B_Simp3;
  
  //get the points on the curve
  #for(T, 0, 1, 0.01)
    sphere{CIget(Test, T), 0.02 pigment{rgb 1}}
  #end
  
  //control points
  sphere{Test[1][0] 0.03 pigment{rgb <1,0,0>}}
  #for(i, 1, dimension_size(Test[1],1) - 1)
    sphere{Test[1][i] 0.03 pigment{rgb <1,0,0>}}
    cylinder{Test[1][i-1], Test[1][i], 0.01 pigment{rgb <1,0,0>}}
  #end
  
  camera {
    perspective angle 95 
    location  <0.5, 0.5,-2>
    right     x*image_width/image_height
    look_at   <0.5, 0.5, 0.0>
  }
  
  light_source{<0, 0,-3000> color rgb 1}  

#end

#version CURVE_INC_TEMP;
#end