#......................
import Blender
import math
from math import *

class XYZ:
   x=0.0
   y=0.0
   z=0.0
   
p=[XYZ,XYZ,XYZ,XYZ]
for n in range(4): print p[n].x,p[n].y,p[n].z

N=40
M=18
PI=3.141592653589793238462643
TWOPI=6.283185307179586476925287

a = 0.9
b = 1.0

# limites pour le parametre t 
tstart = 0
if (a <= b):
   tstop = PI
else:
   tstop = 0.5 * asin(b*b/(a*a))


#fait pivoter un point autour de l'axe x 
def XRotate(p,theta):
   q=XYZ 
   q.x =  p.x
   q.y =  p.y * cos(theta) + p.z * sin(theta)
   q.z = -p.y * sin(theta) + p.z * cos(theta)
   return q

"""
   Evalue la courbe de Cassini à t pour a et b
   t est suppose etre définit dans un ecart correct
   pour a <= b, 0 < t < TWOPI
   pour a >  b, -to < t < to ou to = 0.5 * asin(b^2/a^2)
"""
def Eval(t,a,b,sign):
   p=XYZ
   c1 = b*b*b*b - a*a*a*a*sin(2*t)*sin(2*t)
   if (c1 <= 0): #ne devrait pas survenir dans une utilisation correcte
      c2 = a * a * cos(2*t)
   else:
      c2 = a * a * cos(2*t) + sign * sqrt(c1)
   p.x = cos(t) * sqrt(c2)
   p.y = sin(t) * sqrt(c2)
   p.z = 0.0
   return p

me=Blender.NMesh.GetRaw()

for j in range(M):# Rotation autour de x 
      theta1 = TWOPI*j/M
      theta2 = TWOPI*((j+1)%M)/M
      for i in range(N): # Points sur le plan x-y  
            thesign = 1
            if (a <= b):
                  t1 = tstart + (tstop - tstart) * i / N
                  t2 = tstart + (tstop - tstart) * (i+1) / N
            else:
               if (i < N/2):
                  t1 = tstart + 2 * (tstop - tstart) * i / N
                  t2 = tstart + 2 * (tstop - tstart) * (i+1) / N
               else:
                  t1 = tstart + 2 * (tstop - tstart) * (i-N/2) /N
                  t2 = tstart + 2 * (tstop - tstart) * (i+1-N/2) / N
                  thesign = -1
                 #/* Calcul les 4 points pour construire une facette */
            p[0] = Eval(t1,a,b,thesign)
            p[1] = Eval(t2,a,b,thesign)
            p[2] = p[1]
            p[3] = p[0]
            p[0] = XRotate(p[0],theta1)
            p[1] = XRotate(p[1],theta1)
            p[2] = XRotate(p[2],theta2)
            p[3] = XRotate(p[3],theta2)
            for n in range(4):
                v=Blender.NMesh.Vert(p[n].x,p[n].y,p[n].z)
                me.verts.append(v)
            n=len(me.verts)-4    
            f=Blender.NMesh.Face()  
            f.v.append(me.verts[n])  
            f.v.append(me.verts[n+1])
            f.v.append(me.verts[n+2])
            f.v.append(me.verts[n+3])
            #on l'adoucit
            f.smooth=1 
            me.faces.append(f)
            
Blender.NMesh.PutRaw(me,'cassini')