Option Explicit
'Script written by Ted Ngai Sept 2006 - www.tedngai.net
'Script based on code written by Paul Bourke http://local.wasp.uwa.edu.au/~pbourke/surfaces_curves/supershape3d/ 
'Equations of Superformula are developed by Johan Gielis
'This work is licensed under a  Creative Commons Attribution-Share Alike 3.0 United States License.
'http://creativecommons.org/licenses/by-sa/3.0/us/
'Rhino.Command "_SelAll Delete"
Call Main()
Sub Main() 
	Dim Pi : Pi = 4 * Atn(1)
	Dim arrPoints
	Dim u, v, maxV, maxU : maxV = 40 : maxU = 30
	Dim r, vCyc, uCyc : vCyc = 1 : uCyc = 1
	Dim x,y,z,xt,yt,zt,xa,ya,za
	Dim R1, R2, m, n
	Dim c(1)
	Call Rhino.EnableRedraw(False)
	For v = 0 To maxV-1
		For u =-maxU/2*Pi*uCyc To maxU/2*Pi*uCyc Step Pi/32
			If Abs(v) Mod 2 = 0 Then
				xt=u*2
				yt=Cos(4*u)/2+v
				zt=Cos(4*u)/2+v
			ElseIf Abs(v) Mod 2 = 1 Then
				xt = u*2
				yt=Cos(4*u+Pi)/2+v
				zt=Cos(4*u+Pi)/2+v
			End If
			'Spherical mapping
			'x=SuperRad1(xt/maxU)*cos(xt/maxU)*SuperRad2(yt*vCyc*pi/maxV)*cos(yt*vCyc*pi/maxV)
			'y=SuperRad1(xt/maxU)*sin(xt/maxU)*SuperRad2(yt*vCyc*pi/maxV)*cos(yt*vCyc*pi/maxV)
			'z=SuperRad2(yt*vCyc*pi/maxV)*sin(yt*vCyc*pi/maxV)

			'Toroidal mapping
			R1=18
			R2=8
			x=Cos(xt/maxU)*(R1*SuperRad1(xt/maxU) + R2*SuperRad2(yt*vCyc*Pi/maxV)*Cos(yt*vCyc*Pi/maxV))
			y=Sin(xt/maxU)*(R1*SuperRad1(xt/maxU) + R2*SuperRad2(yt*vCyc*Pi/maxV)*Cos(yt*vCyc*Pi/maxV))
			z=R2*SuperRad2(yt*vCyc*Pi/maxV)*Sin(yt*vCyc*Pi/maxV)*2

			If IsArray(arrPoints) Then 
				ReDim Preserve arrPoints( UBound(arrPoints) + 1 ) 
			Else
				ReDim arrPoints(0) 
			End If 
			arrPoints(UBound(arrPoints))= Array(x,y,z) 
		Next
			Rhino.AddCurve(arrPoints)
			arrPoints = vbNull
	Next
	Call Rhino.EnableRedraw(True)
End Sub

Function SuperRad1(theta)
	Dim n1, n2, n3, m, a, b, t1, t2, Rad1
	'Shpae 3
	n1=0.15
	n2=1.9
	n3=1.6
	m=4
	a=1
	b=1
  
	t1 = (Abs(Cos(m * theta / 4) / a))^n2
	t2 = (Abs(Sin(m * theta / 4) / b))^n3
	Rad1= 1/(t1 + t2)^(1/n1)
	SuperRad1 = Rad1  
End Function

Function SuperRad2(theta)

	Dim n1, n2, n3, m, a, b, t1, t2, Rad2
  
	'Shape 3
	n1=0.2
	n2=1.9
	n3=1.9
	m=5
	a=1
	b=1
	
	t1 = (Abs(Cos(m * theta / 4) / a))^n2
	t2 = (Abs(Sin(m * theta / 4) / b))^n3
	Rad2= 1/(t1+ t2)^(1/n1)
	SuperRad2 = Rad2
End Function