# https://blender.stackexchange.com/questions/82480/how-to-make-a-mobius-strip

import bpy
from math import *
from mathutils import *


def mobius_mesh(resolution, major_radius, minor_radius, thick):

    verts = []
    faces = []

    for i in range(resolution):

        theta = 2*pi * i/resolution
        phi = pi * i/resolution

        rot1 = Matrix.Rotation(phi, 3, [0,1,0])
        rot2 = Matrix.Rotation(theta, 3, [0,0,1])
        c1 = Vector([major_radius, 0, 0])
        v1 = rot2*(c1 + rot1 * Vector([-thick / 2, 0, minor_radius]) )
        v2 = rot2*(c1 + rot1 * Vector([thick / 2, 0, minor_radius]) )
        v3 = rot2*(c1 + rot1 * Vector([thick / 2, 0, -minor_radius]) )
        v4 = rot2*(c1 + rot1 * Vector([-thick / 2, 0, -minor_radius]) )

        i1 = len(verts)
        verts.extend([v1,v2,v3,v4])

        if i+1<resolution:
            ia = i1+4
            ib = i1+5
            ic = i1+6
            id = i1+7
        else:
            ia = 2
            ib = 3
            ic = 0
            id = 1

        # faces.append( [i1+j for j in range(4) ])
        faces.append( [i1,i1+1,ib,ia])
        faces.append( [i1+1,i1+2,ic,ib])
        faces.append( [i1+2,i1+3,id,ic])
        faces.append( [i1+3,i1,ia,id])

    #print (verts)
    #print (faces)
    mesh = bpy.data.meshes.new("mobius")
    mesh.from_pydata(verts, [], faces)
    mesh.validate(True)

    for p in mesh.polygons:
        p.use_smooth=True

    return mesh

def mission1(scn, resolution, major_radius, minor_radius, thick):

    mesh = mobius_mesh(resolution, major_radius, minor_radius, thick)

    obj = bpy.data.objects.new("mobius strip", mesh)
    scn.objects.link(obj)

    mod = obj.modifiers.new("edge split", 'EDGE_SPLIT')


mission1(bpy.context.scene, 36, 5, 1, 0.1)