def flatWorm():
    curveObject = rs.GetObject("pick a backbone curve", 4, True, False)
    samples = rs.GetInteger("# of crosssections", 100, 5)
    bend_radius = rs.GetReal("bend radius", 0.5, 0.001) # r1
    perp_radius = rs.GetReal("ribbon plan radius", 2.0, 0.001) #r2
    crvDom = rs.CurveDomain(curveObject) # domain != length // why do we use domains? == "The domain is the set of all possible input values to the function that defines the curve or surface."

    crossSections = [] # empty array to store sections
    t_step = (crvDom[1] - crvDom[0]) / samples # this is starting to be a pattern!
    t = crvDom[0] # start pt for loop at the start of the line

    for t in rs.frange(crvDom[0], crvDom[1], t_step): # loop thru entire domain w/ floats
        crvCurvature = rs.CurveCurvature(curveObject, t) # evaluate curve with a circle - gives 3d normals & gives radius info
        crossSecPlane = None
        if not crvCurvature:
            crvPoint = rs.EvaluateCurve(curveObject, t)
            crvTan = rs.CurveTangent(curveObject, t) # get tangent vector
            crvPerp = (0,0,1)
            crvNorm = rs.VectorCrossProduct(crvTan, crvPerp)# = product of 2 vectors
            crossSecPlane = rs.PlaneFromFrame(crvPoint, crvPerp, crvNorm)
        else:
            crvPoint = crvCurvature[0]
            crvTan = crvCurvature[1]
            crvPerp = rs.VectorUnitize(crvCurvature[4])
            crvNorm = rs.VectorCrossProduct(crvTan, crvPerp)# look up
            crossSecPlane = rs.PlaneFromFrame(crvPoint, crvPerp, crvNorm)
        if crossSecPlane:
            csec = rs.AddEllipse(crossSecPlane, bend_radius, perp_radius) # draw ellipse at tan/normal to point along curve with radii
            crossSections.append(csec) # add ellipses to an array
        t += t_step # step through domain

    rs.AddLoftSrf(crossSections) # loft list of curves
    rs.DeleteObjects(crossSections) # delete original list of curves as cleanup

  def create_cross_sections(self):
    crvdomain = rs.CurveDomain(self.curve_object)
    self.cross_sections = []
    self.cross_section_planes = []

    t_step = (crvdomain[1]-crvdomain[0])/self.SAMPLES
    pi_step_size = math.pi/self.SAMPLES
    pi_step = 0

    prev_normal = None
    prev_perp = None

    for t in rs.frange(crvdomain[0], crvdomain[1], t_step):
      crvcurvature = rs.CurveCurvature(self.curve_object, t)
      crosssectionplane = None

      if not crvcurvature:
        crosssectionplane = self.cross_section_plane_no_curvature(t, prev_normal, prev_perp)
      else:
        crosssectionplane = self.cross_section_plane_curvature(crvcurvature, prev_normal, prev_perp)

      if crosssectionplane:
        prev_perp = crosssectionplane.XAxis
        prev_normal = crosssectionplane.YAxis
        pi_scalar = self.create_scalar(pi_step)
        radii = self.ellipse_radii(pi_scalar)
        csec = rs.AddEllipse(crosssectionplane, radii[0], radii[1])
        self.cross_sections.append(csec)
        self.cross_section_planes.append(crosssectionplane)
      pi_step += pi_step_size

def ptsOnSrf ():
    surfaceId = rs.GetObject("pick surface", 8, True, True)
    uVal = rs.GetInteger("pick u/row count",4, 1, 20)
    vVal = rs.GetInteger("pick v/col count",4, 1, 20)
    uDomain = rs.SurfaceDomain(surfaceId, 0)
    vDomain = rs.SurfaceDomain(surfaceId, 1)
    uStep = (uDomain[1] - uDomain[0])/ uVal
    vStep = (vDomain[1] - vDomain[0])/ vVal
    count = 0
    
    for i in rs.frange(uDomain[0], uDomain[1], uStep):
        for j in rs.frange(vDomain[0], vDomain[1], vStep):
            point = rs.EvaluateSurface(surfaceId, i, j)
            newPt = rs.AddPoint(point)
            count += 1
            rs.AddText(count,newPt, .25)

def GenerateCrossSection(sectionPoints, radius, rotation, smooth, curve, samples, p, q):
    points = []
    avoidRoundoff = 0.01
    for angle in rs.frange(0.0, 360.0+avoidRoundoff, 360.0/sectionPoints):
        points.append(rs.Polar((0.0, 0.0, 0.0), angle, radius))
    
    rotXform = rs.XformRotation2(rotation, (0, 0, 1), (0, 0, 0))
    points = rs.PointArrayTransform(points, rotXform)
    
    t = curve.Domain[0]
    crossSections = []
    curveCurvature = curve.CurvatureAt(t)
    crossSectionPlane = None
    if not curveCurvature:
        crvPoint = curve.PointAt(t)
        crvTangent = curve.TangentAt(t)
        crvPerp = (0,0,1)
        crvNormal = Rhino.Geometry.Vector3d.CrossProduct(crvTangent, crvPerp)
        crossSectionPlane = Rhino.Geometry.Plane(crvPoint, crvPerp, crvNormal)
    else:
        crvPoint = curve.PointAt(t)
        crvTangent = curve.TangentAt(t)
        crvPerp = curve.CurvatureAt(t)
        crvPerp.Unitize
        crvNormal = Rhino.Geometry.Vector3d.CrossProduct(crvTangent, crvPerp)
        crossSectionPlane = Rhino.Geometry.Plane(crvPoint, crvPerp, crvNormal)
    if crossSectionPlane:
        xform = rs.XformChangeBasis(crossSectionPlane, rs.WorldXYPlane())
        sectionVerts = rs.PointArrayTransform(points, xform)
        if (smooth): # Degree 3 curve to smooth it
            sectionCurve = Rhino.Geometry.Curve.CreateControlPointCurve(sectionVerts, 3)
        else: # Degree 1 curve (polyline)
            sectionCurve = Rhino.Geometry.Curve.CreateControlPointCurve(sectionVerts, 1)
        crossSection = rs.coercecurve(sectionCurve)
    return crossSection

def getHOYs(hours, days, months, timeStep, lb_preparation, method = 0):
    
    if method == 1: stDay, endDay = days
        
    numberOfDaysEachMonth = lb_preparation.numOfDaysEachMonth
    
    if timeStep != 1: hours = rs.frange(hours[0], hours[-1] + 1 - 1/timeStep, 1/timeStep)
    
    HOYS = []
    
    for monthCount, m in enumerate(months):
        # just a single day
        if method == 1 and len(months) == 1 and stDay - endDay == 0:
            days = [stDay]
        elif method == 1:
            #based on analysis period
            if monthCount == 0:
                days = range(stDay, numberOfDaysEachMonth[monthCount] + 1)
            elif monthCount == len(months) - 1: days = range(1, lb_preparation.checkDay(endDay, m) + 1)
            else: days = range(1, numberOfDaysEachMonth[monthCount] + 1)
        
        for d in days:
            for h in hours:
                h = lb_preparation.checkHour(float(h))
                m  = lb_preparation.checkMonth(int(m))
                d = lb_preparation.checkDay(int(d), m)
                HOY = lb_preparation.date2Hour(m, d, h)
                if HOY not in HOYS: HOYS.append(HOY)
    
    return HOYS

def flange(typeofFlange):

    #variables
    points = [] #collection
    nP = 100
    nSin = 10


    #conditional assignment

    if (typeofFlange == "floppy"):
        t = 3
        r = 15
        amp = 10
        
    if (typeofFlange == "extrafloppy"):
        t = 3
        r = 30
        amp = 20

#iteration to calculate points
    for i in rs.frange(0, nP, 1):
        x = r * math.sin(i*360/(nP-1))
        y = r * math.cos(i*360/(nP-1))
        z = amp*math.sin(i*nSin*360/(nP-1))
#Feature
        point = rs.AddPoint([x,y,z])
        points.append(point)

def ptsOnSrf():
    surfaceId = rs.GetObject("pick surface", 8, True, True)
    uVal = rs.GetInteger("pick u/row count", 4, 1, 20)
    vVal = rs.GetInteger("pick v/col count", 4, 1, 20)
    uDomain = rs.SurfaceDomain(surfaceId, 0)
    vDomain = rs.SurfaceDomain(surfaceId, 1)
    uStep = (uDomain[1] - uDomain[0]) / uVal
    vStep = (vDomain[1] - vDomain[0]) / vVal
    count = 0
    allPts = []

    for i in rs.frange(uDomain[0], uDomain[1], uStep):
        for j in rs.frange(vDomain[0], vDomain[1], vStep):
            point = rs.EvaluateSurface(surfaceId, i, j)
            newPt = rs.AddPoint(point)
            allPts.append(newPt)

    rs.AddInterpCrvOnSrf(surfaceId, allPts)
    rs.DeleteObjects(allPts)

def getHOYs(hours, days, months, timeStep, lb_preparation, method = 0):
    
    if method == 1: stDay, endDay = days
        
    numberOfDaysEachMonth = lb_preparation.numOfDaysEachMonth
    
    numOfHours = timeStep * len(hours) 
    
    if timeStep != 1:
        step = 1/timeStep
        hours = rs.frange(hours[0], hours[-1] + 1, step)
        
        # make sure hours are generated correctly
        if len(hours) > numOfHours:
            hours = hours[:numOfHours]
        elif len(hours) < numOfHours:
            newHour = hours[-1] + step
    
    HOYS = []
    
    for monthCount, m in enumerate(months):
        # just a single day
        if method == 1 and len(months) == 1 and stDay - endDay == 0:
            days = [stDay]
        # few days in a single month
        
        elif method == 1 and len(months) == 1:
            days = range(stDay, endDay + 1)
        
        elif method == 1:
            #based on analysis period
            if monthCount == 0:
                # first month
                days = range(stDay, numberOfDaysEachMonth[m-1] + 1)
            elif monthCount == len(months) - 1:
                # last month
                days = range(1, lb_preparation.checkDay(endDay, m) + 1)
            else:
                #rest of the months
                days = range(1, numberOfDaysEachMonth[m-1] + 1)
        
        for d in days:
            for h in hours:
                h = lb_preparation.checkHour(float(h))
                m  = lb_preparation.checkMonth(int(m))
                d = lb_preparation.checkDay(int(d), m)
                HOY = lb_preparation.date2Hour(m, d, h)
                if HOY not in HOYS: HOYS.append(HOY)
    
    return HOYS

 def getValueBasedOnOrientation(valueList):
     angles = []
     if valueList == None or len(valueList) == 0: value = None
     if len(valueList) == 1:
         value = valueList[0]
     elif len(valueList) > 1:
         initAngles = rs.frange(0, 360, 360/len(valueList))
         for an in initAngles: angles.append(an-(360/(2*len(valueList))))
         angles.append(360)
         for angleCount in range(len(angles)-1):
             if angles[angleCount] <= (getAngle2North(normalVector))%360 <= angles[angleCount +1]:
                 targetValue = valueList[angleCount%len(valueList)]
         value = targetValue
     return value

def orientaionStr(divisionAngle, totalAngle):
    divisionAngle = float(divisionAngle)
    totalAngle = float(totalAngle)
    # check the numbers
    if divisionAngle> totalAngle:
        print "Division angle cannot be bigger than the total angle!"
        w = gh.GH_RuntimeMessageLevel.Error
        ghenv.Component.AddRuntimeMessage(w, "Division angle cannot be bigger than the total angle!")
        return 0
    elif totalAngle%divisionAngle!=0:
        print "Total angle should be divisible by the division angle!"
        w = gh.GH_RuntimeMessageLevel.Error
        ghenv.Component.AddRuntimeMessage(w, "Total angle should be divisible by the division angle!")
        return 0
    else:
        return [0] + rs.frange(0, totalAngle, divisionAngle)

  def points_from_ellipse(self, index):
    ellipse = self.cross_sections[index]
    even = index % 2 == 0

    ellipse_domain = rs.CurveDomain(ellipse)
    ellipse_step = (ellipse_domain[1] - ellipse_domain[0])/self.POINTS_ON_CROSS_SECTION
    points = []

    j = 0
    for i in rs.frange(ellipse_domain[0], ellipse_domain[1] - ellipse_step, ellipse_step):
      if even:
        if j % 2 == 0:
          points.append(rs.EvaluateCurve(ellipse, i))
      else:
        if (j + 1) % 2 == 0:
          points.append(rs.EvaluateCurve(ellipse, i))
      j += 1
    return points

def distortArray(fuseList, rec):
    bFuse = rs.BoundingBox(fuseList[0])
    lenFuse = bFuse[1][0] - bFuse[0][0]
    
    bRec = rs.BoundingBox(rec)
    lenRec = bRec[1][0] - bRec[0][0]
    
    xRange = lenRec - lenFuse
    step = xRange/ (len(fuseList))
    
    distortion = rs.frange(0, xRange, step)
   
    for i in range(0, len(fuseList)):
        x = choice(distortion)
        distortion.remove(x)
        
        xform = [x,0,0]
        rs.MoveObject(fuseList[i], xform)
        
    return fuseList 

def TorusKnotVerts(p, q, res, zScale):
    # Determine if we need to loop less because of non coprime p and q
    upLimit = 2.0*math.pi
    P = int(p); Q = int(q)
    hcf = HighestCommonFactor(P, Q)
    if (hcf != 1):
        upLimit /= hcf
    
    verts = []
    roundTol = 0.01 # Fudge factor to make sure loop closes
    tStep = math.pi/(res/2.0)
    for t in rs.frange(0.0, upLimit+roundTol, tStep):
        r = math.cos(q*t)+2.0
        x = r*math.cos(p*t)
        y = r*math.sin(p*t)
        z = -math.sin(q*t)*2.0
        z *= zScale
        pt = Rhino.Geometry.Point3d(x, y, z)
        verts.append(pt)
    return verts

	def __generate_form_levels(self, spine_curve):
		crvdomain = rs.CurveDomain(spine_curve)
		printedState = ""
		crosssection_planes = []
		crosssection_plane_nums = []
		crosssections = []
		t_step = (crvdomain[1] - crvdomain[0]) / (self.object_properties["number_of_lofts"]-1)
		t = crvdomain[0]
		for t in rs.frange(crvdomain[0], crvdomain[1], t_step):
			if(self.emotion_properties["vertical_AR"][str(int(t+1))] != None):		
				crvcurvature = rs.CurveCurvature(spine_curve, t)
				crosssectionplane = None
				if not crvcurvature:
					crvPoint = rs.EvaluateCurve(spine_curve, t)
					crvTangent = rs.CurveTangent(spine_curve, t)
					crvPerp = (0,0,1)
					crvNormal = rs.VectorCrossProduct(crvTangent, crvPerp)
					printedState = printedState + str(crvNormal)
					crosssectionplane = rs.PlaneFromNormal(crvPoint, crvTangent)
					if(t==0):
						crosssectionplane = rs.PlaneFromNormal([0,0,0], [0,0,1])
				else:
					crvPoint = crvcurvature[0]
					crvTangent = crvcurvature[1]
					crvPerp = rs.VectorUnitize(crvcurvature[4])
					crvNormal = rs.VectorCrossProduct(crvTangent, crvPerp)
					printedState = printedState + str(crvNormal)
					crosssectionplane = rs.PlaneFromNormal(crvPoint, crvTangent, crvNormal)
				if crosssectionplane:
					crosssection_planes.append(crosssectionplane)
					crosssection_plane_nums.append(str(int(t+1)))
		if len(crosssection_plane_nums) > 0:
			last_element = crosssection_plane_nums.pop(len(crosssection_plane_nums)-1)
			crosssection_plane_nums.insert(0, last_element)
		for index in xrange(len(crosssection_plane_nums)):
			crosssections.append(self.__generate_individual_levels(crosssection_planes[index], crosssection_plane_nums[index]))
		if not crosssections: return
		crosssections.append(crosssections.pop(0))
		rs.AddLoftSrf(crosssections,closed=False,loft_type=int(round(self.emotion_properties["vertical_wrapping"])))
				
		return crosssection_planes[0]

def contour(crvOffset):
    # get geometry
    surfaceId = rs.GetObject("pick surface to contour", 0, True, True)
    startPt = rs.GetPoint("base point of contours")
    endPt = rs.GetPoint("end point of contours")
    count = 0
    reference = []
    target = []

    # make contours & store in newCrvs
    newCrvs = rs.AddSrfContourCrvs(
        surfaceId, (startPt, endPt), crvOffset
    )  # output is a list of GUIDs. can't access raw points

    # divide the target surface
    printBed = rs.GetObject("pick surface for layout", 8, True, True)
    intCount = len(newCrvs)
    uDomain = rs.SurfaceDomain(printBed, 0)
    vDomain = rs.SurfaceDomain(printBed, 1)
    uStep = (uDomain[1] - uDomain[0]) / intCount

    for u in rs.frange(uDomain[0], uDomain[1], uStep):
        layout = rs.SurfaceFrame(printBed, [u, 1])
        target1 = layout[0]  # set target to point inside of object - note this is a single point
        target2 = rs.PointAdd(target1, (0, 10, 0))  # this could be more parametric
        target.extend([target1, target2])
    # print target

    # add text, reference and orient!
    # for orient, we need a list 3 points to reference and 3 points to target!
    # maybe reference should be curve origin crvPl[0] and midpoint? or else polyline verticies -- need to convert curve to polyline first?
    for crv in newCrvs:
        count += 1  # for label
        crvPl = rs.CurvePlane(crv)  # plane for text
        rs.AddText(count, crvPl, 0.25)  # should you label on the ground?
        # crvPl = rs.PointAdd(crvPl[0], (0,0,-1)) # if you wanted to offset text
        ref1 = rs.CurveMidPoint(crv)
        ref2 = crvPl[0]
        reference.extend([ref1, ref2])

def analyzeGlz(glzSrf, distBetween, numOfShds, horOrVertical, lb_visualization, normalVector):
    # find the bounding box
    bbox = glzSrf.GetBoundingBox(True)
    if horOrVertical == None:
        horOrVertical = True
    if numOfShds == None and distBetween == None:
        numOfShds = 1
    
    if numOfShds == 0 or distBetween == 0:
        sortedPlanes = []
    
    elif horOrVertical == True:
        # Horizontal
        #Define a bounding box for use in calculating the number of shades to generate
        minZPt = bbox.Corner(False, True, True)
        minZPt = rc.Geometry.Point3d(minZPt.X, minZPt.Y, minZPt.Z)
        maxZPt = bbox.Corner(False, True, False)
        maxZPt = rc.Geometry.Point3d(maxZPt.X, maxZPt.Y, maxZPt.Z - sc.doc.ModelAbsoluteTolerance)
        centerPt = bbox.Center 
        #glazing hieghts
        glzHeight = minZPt.DistanceTo(maxZPt)
        
        # find number of shadings
        try:
            numOfShd = int(numOfShds)
            shadingHeight = glzHeight/numOfShd
            shadingRemainder = shadingHeight
        except:
            shadingHeight = distBetween
            shadingRemainder = (((glzHeight/distBetween) - math.floor(glzHeight/distBetween))*distBetween)
            if shadingRemainder == 0:
                shadingRemainder = shadingHeight
        
        # find shading base planes
        planeOrigins = []
        planes = []
        X, Y, z = minZPt.X, minZPt.Y, minZPt.Z
        zHeights = rs.frange(minZPt.Z + shadingRemainder, maxZPt.Z + 0.5*sc.doc.ModelAbsoluteTolerance, shadingHeight)
        try:
            for Z in zHeights:
                planes.append(rc.Geometry.Plane(rc.Geometry.Point3d(X, Y, Z), rc.Geometry.Vector3d.ZAxis))
        except:
            # single shading
            planes.append(rc.Geometry.Plane(rc.Geometry.Point3d(maxZPt), rc.Geometry.Vector3d.ZAxis))
        # sort the planes
        sortedPlanes = sorted(planes, key=lambda a: a.Origin.Z)
    
    elif horOrVertical == False:
        # Vertical
        # Define a vector to be used to generate the planes
        planeVec = rc.Geometry.Vector3d(normalVector.X, normalVector.Y, 0)
        planeVec.Rotate(1.570796, rc.Geometry.Vector3d.ZAxis)
        
        
        #Define a bounding box for use in calculating the number of shades to generate
        minXYPt = bbox.Corner(True, True, True)
        minXYPt = rc.Geometry.Point3d(minXYPt.X, minXYPt.Y, minXYPt.Z)
        maxXYPt = bbox.Corner(False, False, True)
        maxXYPt = rc.Geometry.Point3d(maxXYPt.X, maxXYPt.Y, maxXYPt.Z)
        centerPt = bbox.Center
        
        #Test to be sure that the values are parallel to the correct vector.
        testVec = rc.Geometry.Vector3d.Subtract(rc.Geometry.Vector3d(minXYPt.X, minXYPt.Y, minXYPt.Z), rc.Geometry.Vector3d(maxXYPt.X, maxXYPt.Y, maxXYPt.Z))
        if testVec.IsParallelTo(planeVec) == 0:
            minXYPt = bbox.Corner(False, True, True)
            minXYPt = rc.Geometry.Point3d(minXYPt.X, minXYPt.Y, minXYPt.Z)
            maxXYPt = bbox.Corner(True, False, True)
            maxXYPt = rc.Geometry.Point3d(maxXYPt.X, maxXYPt.Y, maxXYPt.Z)
        
        #Adjust the points to ensure the creation of the correct number of shades starting from the northernmost side of the window.
        tolVec = rc.Geometry.Vector3d.Subtract(rc.Geometry.Vector3d(minXYPt.X, minXYPt.Y, minXYPt.Z), rc.Geometry.Vector3d(maxXYPt.X, maxXYPt.Y, maxXYPt.Z))
        tolVec.Unitize()
        tolVec = rc.Geometry.Vector3d.Multiply(sc.doc.ModelAbsoluteTolerance*2, tolVec)
        
        if tolVec.X > 0 and  tolVec.Y > 0:
            tolVec = rc.Geometry.Vector3d.Multiply(1, tolVec)
            norOrient = False
        if tolVec.X < 0 and  tolVec.Y > 0:
            tolVec = rc.Geometry.Vector3d.Multiply(1, tolVec)
            norOrient = False
        if tolVec.X < 0 and  tolVec.Y < 0:
            tolVec = rc.Geometry.Vector3d.Multiply(-1, tolVec)
            norOrient = True
        else:
            tolVec = rc.Geometry.Vector3d.Multiply(-1, tolVec)
            norOrient = True
        
        maxXYPt = rc.Geometry.Point3d.Subtract(maxXYPt, tolVec)
        minXYPt = rc.Geometry.Point3d.Subtract(minXYPt, tolVec)
        
        #glazing distance
        glzHeight = minXYPt.DistanceTo(maxXYPt)
        
        # find number of shadings
        try:
            numOfShd = int(numOfShds)
            shadingHeight = glzHeight/numOfShd
            shadingRemainder = shadingHeight
        except:
            shadingHeight = distBetween
#.........这里部分代码省略.........

#3,8 Torus knot
import rhinoscriptsyntax as rs
import math


points = [] #declares the size of my point array

pi = math.pi

#rs.EnableRedraw(False)

for t in rs.frange(0,360,1):

    rads = t*pi/2 #defines the rotation angle which is equal to time t multiplied by angular velocity

#mathematical description of a 3,8 Torus knot
    x = (3+ 2*math.cos(rads))*math.cos(t)
    y = (3+2*math.cos(rads))*math.sin(t)
    z = 2*math.sin(rads)

    n = rs.AddPoint([x, y, z])
    points.append(n)

rs.AddInterpCurve(points)

#rs.EnableRedraw(True)

import math
import rhinoscriptsyntax as rs


dblA = -8.0
dblB = 8.0
dblStep = 0.25
for x in rs.frange(dblA, dblB, dblStep):
    y = 2*math.sin(x)
    rs.AddPoint([x, y, 0])

mpoints = rs.ObjectsByType(1, True)


mcurves = rs.AddCurve(mpoints)

rs.DeleteObjects(mpoints)

rs.CopyObject( mcurves, [0,1,0] )

#.........这里部分代码省略.........
                print 'Checking conditional statements...'
                # send all data and statement to a function and return back
                # True, False Pattern and condition statement
                titleStatement, patternList = checkConditionalStatement(annualHourlyData, conditionalStatement)
            
            if titleStatement != -1 and True not in patternList:
                warning = 'No hour meets the conditional statement.' 
                print warning
                ghenv.Component.AddRuntimeMessage(gh.GH_RuntimeMessageLevel.Warning, warning)
                return -1
            
            if titleStatement == -1:
                patternList = [[True]] * 8760
                titleStatement = False
            
           # check the scale
            try:
                if float(scale)!=0:
                    try:scale = 10*float(scale)/conversionFac
                    except: scale = 10/conversionFac
                else: scale = 10/conversionFac
            except: scale = 10/conversionFac
            
            cenPt = lb_preparation.getCenPt(centerPoint)
            
            if not numOfDirections or int(numOfDirections) < 4: numOfDirections = 16
            elif int(numOfDirections)> 360: numOfDirections = 360
            else:
                try: numOfDirections = int(numOfDirections)
                except: numOfDirections = 16
            
            # define angles
            segAngle = 360/numOfDirections
            roseAngles = rs.frange(0,360,segAngle);
            if round(roseAngles[-1]) == 360: roseAngles.remove(roseAngles[-1])
            
            movingVectors = []; sideVectors = []
            northAngle, northVector = lb_preparation.angle2north(north)
            for i, angle in enumerate(roseAngles):
                northVector1 = rc.Geometry.Vector3d(northVector)
                northVector2 = rc.Geometry.Vector3d(northVector)
                northVector1.Rotate(-float(math.radians(angle)), rc.Geometry.Vector3d.ZAxis)
                movingVectors.append(northVector1)
                northVector2.Rotate(-float(math.radians(angle + (segAngle/2))), rc.Geometry.Vector3d.ZAxis)
                sideVectors.append(northVector2)
            
            HOYS, months, days = lb_preparation.getHOYsBasedOnPeriod(analysisPeriod, 1)
            selectedWindDir = []
            for count in HOYS:
                selectedWindDir.append(windDir[count-1])
            #selectedWindDir = lb_preparation.selectHourlyData(windDir, analysisPeriod)[7:]
            # read analysis period
            stMonth, stDay, stHour, endMonth, endDay, endHour = lb_preparation.readRunPeriod(analysisPeriod, False)

            # find the study hours based on wind direction data
            startHour = lb_preparation.date2Hour(stMonth, stDay, stHour)
            endingHour =  lb_preparation.date2Hour(endMonth, endDay, endHour)
            if startHour <= endingHour: studyHours = range(startHour-1, endingHour)
            else: studyHours = range(startHour - 1, 8760) + range(0, endingHour)
            
            calmHour = [] # count hours with no wind
            separatedBasedOnAngle = []
            [separatedBasedOnAngle.append([]) for i in range(len(roseAngles))]
            #print len(studyHours)
            #print len(selectedWindDir)
            for hour, windDirection in enumerate(selectedWindDir):

idSurface = rs.GetObject("srf to frame? ", 8, True, True)
intCount = rs.GetReal("# of itera. per direction")

#domains are rough dims of the surface - can check via _what 
uDomain = rs.SurfaceDomain(idSurface, 0)
vDomain = rs.SurfaceDomain(idSurface, 1)

# get size int from domain & / by intended count = increment between surfaces
uStep = (uDomain[1] - uDomain[0]) / intCount
vStep = (vDomain[1] - vDomain[0]) / intCount

#print uStep
#print vStep

rs.EnableRedraw(False)

# loop through size of surface by step increment in both u & v directions
for u in rs.frange(uDomain[0], uDomain[1], uStep):
    for v in rs.frange(vDomain[0], vDomain[1], vStep):
        pt = rs.EvaluateSurface(idSurface, u, v) # get points at each domain step
        if rs.Distance(pt, rs.BrepClosestPoint(idSurface, pt)[0]) < 0.1:
            srfFrame = rs.SurfaceFrame(idSurface, [u, v]) # create ref plane at each division point
            newPlane = rs.AddPlaneSurface(srfFrame, 1.0, 1.0) # create surface at each ref plane
            # add height guideline from plane's origin to "ceiling"
            centerPt = rs.SurfaceAreaCentroid(newPlane) # locate center of each new plane
            limit = 10 # set "ceiling"
            endPt = limit - centerPt[0][2] # access z coordinate with the tuple
            rs.AddLine(centerPt[0], rs.PointAdd(centerPt[0], (0,0,endPt)) )
            
rs.EnableRedraw(True) # refresh viewport at one time only