Option Explicit
'Script written by Will Sixel
'Script copyuppered by Intact Solutions, LLC
'Script version Monday, June 22, 2015 10:06:28 AM

'Function: Returns the maximum displacement found in the solid passed to it, for the given loading'
Function totalDisplacement(angle, SNS, solid, restraintFace, loadFace, resolution)
	totalDisplacement = -1
	Dim load
	
	'Define the loading based on angle'
	load = Array(Cos(angle * Rhino.PI / 180), Sin(angle * Rhino.PI / 180), 0)
	
	'Define material, restaint, load, and then find the solution'
	Call SNS.SetMaterial(solid, "Aluminum 6061")
	Call SNS.AddFaceRestraint(solid, "resEnd", restraintFace, "xyz")
	Call SNS.AddFaceVectorLoad(solid, "endLoad", loadFace, load, False)
	Call SNS.Solve(solid, resolution)
	
	'Query the solution to find maximum displacement'
	totalDisplacement = SNS.QuerySolutionMax(solid, "DSPT")
	
	'Reset SnS data for the next iteration'
	Call SNS.Reset(solid)
	
End Function


'Main loop'
Call Main()
Sub Main()
	'Create the angle iron geometry by extrusion of rectangles and boolean subtraction'
	Dim length, width, height, thickness, points, rec1, rec2, box1, box2, angleIron, extrusionLine
	length = 50
	width = 5
	height = 5
	thickness = 1
	
	'Make the line the extrude along'
	extrusionLine = Rhino.AddLine(Array(0, 0, 0), Array(0, 0, length))
	
	'Define the first, larger box'
	points = Array(Array(0, 0, 0), Array(width, 0, 0), Array(width, height, 0), Array(0, height, 0))
	rec1 = Rhino.AddSrfPt(points)
	box1 = Rhino.ExtrudeSurface(rec1, extrusionLine)
	
	'Define the second, smaller box to subtract with'
	points = Array(Array(thickness, thickness, 0), Array(width, thickness, 0), Array(width, height, 0), Array(thickness, height, 0))
	rec2 = Rhino.AddSrfPt(points)
	box2 = Rhino.ExtrudeSurface(rec2, extrusionLine)
	
	'Take the boolean differenece to create the L-shaped cross section'
	angleIron = Rhino.BooleanDifference(box1, box2)
	
	'Delete non-solid objects'
	Call Rhino.DeleteObjects(Array(rec1, rec2, extrusionLine))
	
	'Set up SnS'
	Dim restraintFace, loadFace, load, SNS, resolution
	Set SNS = Rhino.GetPluginObject("SnSScript")
	
	'Find the indices of restraint faces'
	restraintFace = SNS.GetNearestFace(angleIron(0), Array(thickness * 0.5, thickness * 0.5, Length * 1.01))
	loadFace = SNS.GetNearestFace(angleIron(0), Array(thickness * 0.5, thickness * 0.5, Length * -0.01))
	
	'Define load and maximum number of loops to run'
	Const loops = 20
	resolution = 2000
	
	'Declare variables to be used in loop'
	Dim upper, lower, phi, delta, tolerance, runs, tau
	Dim vUp, vLow, vX1, vX2, x1, x2, angle, vLast
	
	'Initialize variables:'
	runs = 0				'Iteration counter'
	delta = 100				'Convergence variable'
	tolerance = 0.0001		'Compared to delta between interior values'
	lower = 0				'Lower bound of domain:  0 degrees'
	upper = 180				'Upper bound of domain: 180 degrees'
	vX1 = 1					'Initialize'
	
	'Golden ratio'
	phi = (1 + Sqr(5)) / 2
	tau = 1 / phi
	
	'Boolean variables: used to note which boundary changes, '
	' and thus which boundary does not need to be re-computed'
	Dim upperChange, lowerChange
	upperChange = True
	lowerChange = True
	
	'Golden Search Loop: Run until tolerance is reached, or maximum iterations'
	While ( (delta > tolerance) And (runs < loops) )

		'Incrment runs on start of loop; print number to command line'
		runs = runs + 1
		Call Rhino.Print(runs)

		'Define the inner boundaries per golden ration algorithm'
		'x1 is located at a ratio of ~0.382 away from lower boundary'
		'x2 is located at ~0.618 of the upper boundary'
		x1 = lower + (1 - tau) * (upper - lower)
		x2 = upper - (1 - tau) * (upper - lower)
		
		'Check to see if the upper boundary changed; vX2 becomes old vX1, and need to recompute vX1
		If (upperChange) Then
			vX1 = totalDisplacement(x1, SNS, angleIron(0), restraintFace, loadFace, resolution)
		End If
		'Check to see if lower boundary has changed; then need to recompute vX2'
		If (lowerChange) Then
			vX2 = totalDisplacement(x2, SNS, angleIron(0), restraintFace, loadFace, resolution)
		End If
		
		'Delta: Calculate percent difference between the two inner values'
		delta = Abs((vX2 - vX1) / vX1)
		
		'Search algorithm:'
		If (vX1 > vX2) Then
			'Case 1:  X1 becomes new lower bound '
			'the value at x2 becomes the new value of x1.'
			lower = x1
			vX1 = vX2
			upperChange = False
			lowerChange = True
		ElseIf (vX1 < vX2) Then
			'Case 2:  X2 becomes new upper bound.'
			'Value at x2 becomes the new x1 value'
			upper = x2
			upperChange = True
			lowerChange = False
			vX2 = vX1

		End If
		 
		
	Wend
	
	x1 = (x1 + x2) / 2
	If (runs < loops) Then
		Call Rhino.MessageBox("Reached tolerance in " & runs & " iterations.")
	Else
		Call Rhino.MessageBox("Reached maximum number of runs. Aborting.")
	End If
		
	Call Rhino.MessageBox("Delta: " & delta & "Angle: " & x1) 
	
	
End Sub