Option Explicit
'Script written by Luis Gil

'This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. 
'To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/3.0

'www.legildesign.com

Call Main()
Sub Main()
	
	Dim dblEdge
	dblEdge = Rhino.GetReal("Set edge length for genesis triangle.", 100)
	If IsNull(dblEdge) Then Exit Sub
	
	Dim intIterations
	intIterations = Rhino.GetInteger("How Many iterations?", 5)
	If IsNull(intIterations) Then Exit Sub
	
	Dim dblAngle
	dblAngle = Rhino.GetReal("Set incident angle for adjacent faces.", 30)
	If IsNull(dblAngle) Then Exit Sub
	
	Dim dblH
	dblH = dblEdge * sin(Rhino.Pi/3)
	
	Dim arrPts(2)
	arrPts(0) = Array(0,0,0)
	arrPts(1) = Array(dblEdge, 0,0)
	arrPts(2) = Array(dblEdge/2, dblH, 0)
	
	Dim arrStartSrf(0)
	arrStartSrf(0) = Rhino.AddSrfPt(arrPts)
	
	Dim ID : ID = 0
	
	Rhino.EnableRedraw False
	Call SplitAndGrow(arrStartSrf, dblAngle, intIterations, ID)
	Rhino.EnableRedraw True

End Sub

Sub SplitAndGrow(arrSrf, dblAngle, intLimit, intID)
	
	If intID > intLimit Then Exit Sub
	
	Rhino.Print "On iteration " & intID & " of " & intLimit
	
	Dim arrSurfaces()
	Dim index : index = -1
	
	Dim s
	For Each s In arrSrf
		'the 4 new surfaces. 0 is the central; 1,2, & 3 are peripheral with 1 being
		'the lower left and progressing counter-clockwise
		ReDim arrNewSrfs(3)
		
		Dim arrSrfPts
		arrSrfPts = Rhino.SurfacePoints(s)
		
		'find the midpoints each edge to construct the four smaller triangles
		Dim arrMidPt1, arrMidPt2, arrMidPt3
		arrMidPt1 = Rhino.PointDivide(Rhino.PointAdd(arrSrfPts(0), arrSrfPts(1)), 2)
		arrMidPt2 = Rhino.PointDivide(Rhino.PointAdd(arrSrfPts(1), arrSrfPts(2)), 2)
		arrMidPt3 = Rhino.PointDivide(Rhino.PointAdd(arrSrfPts(0), arrSrfPts(2)), 2)
		
		arrNewSrfs(0) = Rhino.AddSrfPt(Array(arrMidPt1, arrMidPt3, arrMidPt2))
		arrNewSrfs(1) = Rhino.AddSrfPt(Array(arrSrfPts(0), arrMidPt3, arrMidPt1))
		arrNewSrfs(2) = Rhino.AddSrfPt(Array(arrSrfPts(1), arrMidPt1, arrMidPt2))
		arrNewSrfs(3) = Rhino.AddSrfPt(Array(arrSrfPts(2), arrMidPt2, arrMidPt3))
		
		
		Dim dblAngleRAD : dblAngleRAD = dblAngle/180 * Rhino.Pi
		'set up the rotation axes for the 3 peripheral surfaces
		Dim arrRotVec1, arrRotVec2, arrRotVec3
		arrRotVec1 = Rhino.VectorCreate(arrMidPt3, arrMidPt1)
		arrRotVec2 = Rhino.VectorCreate(arrMidPt1, arrMidPt2)
		arrRotVec3 = Rhino.VectorCreate(arrMidPt2, arrMidPt3)
		
		arrNewSrfs(1) = Rhino.RotateObject(arrNewSrfs(1),arrMidPt1,dblAngle, arrRotVec1)
		arrNewSrfs(2) = Rhino.RotateObject(arrNewSrfs(2),arrMidPt2,dblAngle, arrRotVec2)
		arrNewSrfs(3) = Rhino.RotateObject(arrNewSrfs(3),arrMidPt3,dblAngle, arrRotVec3)
		
		Dim arrCenPt 'centroid
		arrCenPt = Rhino.SurfaceAreaCentroid(s)
		
		Dim arrCenParam
		arrCenParam = Rhino.SurfaceClosestPoint(s, arrCenPt(0))
		
		Dim arrNormal
		arrNormal = Rhino.SurfaceNormal(s,arrCenParam)
		
		'the offset height is set off of the desired angle of contact between faces
		'first the distance must be found from a point to a perp bisector for the smaller faces
		'then this number is used to calculate the offset
		
		Dim dblNH : dblNH = Rhino.Distance(arrSrfPts(0), arrMidPt1) * sin(dblAngleRAD)
		Dim dblOffset
		dblOffset = dblNH * sin(dblAngleRAD)
				
		arrNormal = Rhino.VectorScale(arrNormal, dblOffset)
		
		Rhino.MoveObjects arrNewSrfs, arrNormal
		
		index = index + 4 'adding 4 new surfaces to the total array at a time
		
		ReDim Preserve arrSurfaces(index)
		arrSurfaces(UBound(arrSurfaces) - 3) = arrNewSrfs(0)
		arrSurfaces(UBound(arrSurfaces) - 2) = arrNewSrfs(1)
		arrSurfaces(UBound(arrSurfaces) - 1) = arrNewSrfs(2)
		arrSurfaces(UBound(arrSurfaces)) = arrNewSrfs(3)

		
		Rhino.DeleteObject s
		
	Next
	
	Dim intNewID : intNewID = intID + 1
	
	Call SplitAndGrow(arrSurfaces,dblAngle, intLimit, intNewID)


End Sub