//-------------------------------------------------------------//
// File created by Michael Andrews (M.C.Andrews@reading.ac.uk) //
//                13/8/99 Release version 1.0                  //
//-------------------------------------------------------------//

//=========================================
// macro mAlign - from Gilles Trans' Maketree.pov
//-----------------------------------------
// returns a matrix operation that aligns an object or texture along P2-P1
// the object is translated to P1
// translate to P2-P1 if you want the object to be on P2
#macro mAlign(P1,P2)
	#local yV1=vnormalize(P2-P1);
	#local xV1=vnormalize(vcross(yV1,z));
	#local zV1=vcross(xV1,yV1);
	matrix <xV1.x,xV1.y,xV1.z,yV1.x,yV1.y,yV1.z,zV1.x,zV1.y,zV1.z,P1.x,P1.y,P1.z>
#end   
//-----------------------------------------
// end of mAlign macro
//=========================================

// spline interpolation arrays
// defined for 1 to 4 values
//
#declare INTERP = 0;
#declare mA = array[4]
#declare mA[0] = array[1][1] {{1}}
#declare mA[1] = array[2][2] {{1,0},{-1,1}}
#declare mA[2] = array[3][3] {{1,0,0},{-3,4,-1},{2,-4,2}}
#declare mA[3] = array[4][4] {{1,0,0,0},{-5.5,9,-4.5,1},{9,-22.5,18,-4.5},{-4.5,13.5,-13.5,4.5}}

#declare BEZIER = 1;
#declare bA = array[4]
#declare bA[0] = array[1][1] {{1}}
#declare bA[1] = array[2][2] {{1,0},{-1,1}}
#declare bA[2] = array[3][3] {{1,0,0},{-2,2,0},{1,-2,1}}
#declare bA[3] = array[4][4] {{1,0,0,0},{-3,3,0,0},{3,-6,3,0},{-1,3,-3,1}}

// M_interp(vP, fT, iType)
// 
// macro which performs the actual matrix mult
//
// vP = array to interpolate, vector or float, size 1 to 4
// fT = interpolation point, float [0,1]
// iType = interpolation type: 0 = INTERP = interpolating, 1 = BEZIER = bezier
//
#macro M_interp ( vP, fT, iType )
	#local dS = dimension_size(vP, 1);
	#local viT = array[dS]
	#local viT[0] = 1;
	#local viC = 1; #while(viC < dS)
		#local viT[viC] = viT[viC-1]*fT;
	#local viC = viC + 1; #end
	#local viTemp = vP[0]*0;
	#switch (iType)
		#case (INTERP)
			#local vmR = 0; #while(vmR < dS)
				#local vmC = 0; #while(vmC < dS)
					#local viTemp = viTemp + viT[vmR] * mA[dS-1][vmR][vmC] * vP[vmC];
				#local vmC = vmC + 1; #end
			#local vmR = vmR + 1; #end
		#break
		#case (BEZIER)
			#local vmR = 0; #while(vmR < dS)
				#local vmC = 0; #while(vmC < dS)
					#local viTemp = viTemp + viT[vmR] * bA[dS-1][vmR][vmC] * vP[vmC];
				#local vmC = vmC + 1; #end
			#local vmR = vmR + 1; #end
		#break
	#else
		#error concat("Unknown spline type: ", str(iType,0,0), ".\n")
	#end
	viTemp
#end

// M_BSphStack ( vA, vIType, rA, rIType, sA, sIType, S )
// 
// macro which builds spline stack of blob elipsoid components
//
// vA = position array to interpolate, vector or float, size 1 to 4
// vIType = interpolation type for positions
// rA = component scaling array to interpolate, vector or float, size 1 to 4
// rIType = interpolation type for scaling
// sA = component strength array to interpolate, vector or float, size 1 to 4
// sIType = interpolation type for strengths
// S = component separation, 1 = components touching, 0 = no separation of centres, best to use 0.1 to 0.3
//
#macro M_BSphStack ( vA, vIType, rA, rIType, sA, sIType, S )
	#local F = 2/S+2;
	#local DIV = (1/5-2/3+1)*F-1;
	#local T = 0.0; #while ( T <= 1.0 )
		#local VP = M_interp ( vA, T, vIType );
		#local RP = M_interp ( rA, T, rIType );
		#local SP = M_interp ( sA, T, sIType );
		#local VLP = M_interp ( vA, T+0.01, vIType );
		#local VLM = M_interp ( vA, T-0.01, vIType );
		#local VG = (VLP-VLM)/0.02;
		#if (abs(SP/DIV) > 0.01 & vlength(RP) > 0.01)
			sphere { VP/RP, 1, SP/DIV scale RP }
		#end
	#local T = T + S * (0.01+vlength(RP)) / vlength(VG); #end
#end
		
// M_BSphAlign ( vA, vIType, rA, rIType, sA, sIType, S )
// 
// macro which builds spline stack of blob elipsoid components,
//   aligns y scaling along tangent of spline
//
// vA = position array to interpolate, vector or float, size 1 to 4
// vIType = interpolation type for positions
// rA = component scaling array to interpolate, vector or float, size 1 to 4
// rIType = interpolation type for scaling
// sA = component strength array to interpolate, vector or float, size 1 to 4
// sIType = interpolation type for strengths
// S = component separation, 1 = components touching, 0 = no separation of centres, best to use 0.1 to 0.3
//
#macro M_BSphAlign ( vA, vIType, rA, rIType, sA, sIType, S )
	#local F = 2/S+2;
	#local DIV = (1/5-2/3+1)*F-1;
	#local T = 0.0; #while ( T <= 1.0 )
		#local VP = M_interp ( vA, T, vIType );
		#local RP = M_interp ( rA, T, rIType );
		#local SP = M_interp ( sA, T, sIType );
		#local VLP = M_interp ( vA, T+0.01, vIType );
		#local VLM = M_interp ( vA, T-0.01, vIType );
		#local VG = (VLP-VLM)/0.02;
		#if (abs(SP/DIV) > 0.01 & vlength(RP) > 0.01)
			sphere {
				0, 1, SP/DIV
				scale RP
				mAlign(VP, VP+VG)
			}
		#end
	#local T = T + S * (0.01+abs(RP.y)) / vlength(VG); #end
#end
		
// M_BCylStack ( vAa, vAb, vIType, rA, rIType, sA, sIType, S )
// 
// macro which builds stack of blob cylinder components connecting two splines
//
// vAa = first position array to interpolate, vector or float, size 1 to 4
// vAb = secondposition array to interpolate, vector or float, size 1 to 4
// vIType = interpolation type for positions
// rA = component scaling array to interpolate, vector or float, size 1 to 4
// rIType = interpolation type for scaling
// sA = component strength array to interpolate, vector or float, size 1 to 4
// sIType = interpolation type for strengths
// S = component separation, 1 = components touching, 0 = no separation of centres, best to use 0.1 to 0.3
//
#macro M_BCylStack ( vAa, vAb, vIType, rA, rIType, sA, sIType, S )
	#local F = 2/S+2;
	#local DIV = (1/5-2/3+1)*F-1;
	#local T = 0.0; #while ( T <= 1.0 )
		#local VPA = M_interp ( vAa, T, vIType );
		#local VPB = M_interp ( vAb, T, vIType );
		#local RP = M_interp ( rA, T, rIType );
		#local SP = M_interp ( sA, T, sIType );
		#local VLPA = M_interp ( vAa, T+0.01, vIType );
		#local VLMA = M_interp ( vAa, T-0.01, vIType );
		#local VLPB = M_interp ( vAb, T+0.01, vIType );
		#local VLMB = M_interp ( vAb, T-0.01, vIType );
		#local VGA = (VLPA-VLMA)/0.02;
		#local VGB = (VLPB-VLMB)/0.02;
		#local VG = max(vlength(VGA),vlength(VGB));
		#if (abs(SP/DIV) > 0.01 & vlength(RP) > 0.01 & vlength(VPB-VPA) > 0.01)
			cylinder { VPA/RP, VPB/RP, 1, SP/DIV scale RP }
		#end
	#local T = T + S * (0.01+vlength(RP)) / VG; #end
#end