User:Thierry Dugnolle/Python/High definition 3D line drawer

A Chua attractor


Remark

edit

I couldn't run and debug this program because my computer is hacked. The image above is an old one (2018) made with an old version (C#) of this program.

main.py

edit
# These files makes a drawer who draws images of 3D lines.
# To see the depth of space, what is in front of us shall hide what is behind.
# Ordinary visualizations of 3D lines don't do this.
# The trick is to surround a line with a thin halo, so that it hides the lines behind.

from LineDrawer import aLineDrawer
from Vector3D import aNew3Dvector
from ODE import aNewChuaSolution


print("Line3D_Drawer")

# The drawer draws a solution of a Chua system of ordinary differential equations (ODE).
TheDrawer = aLineDrawer()
    # The image:
TheWidth_inNumberOfPixels = 500
TheHeight_inNumberOfPixels = 500
    # The line of vision:
Theta = 90.0 # the angle (in degrees) of the line of vision relative to the z axis.
Phi = 120.0 # the angle (in degrees) between the projection of the line of vision on the xy plane and the x axis:
    # The position of the optical center is determined by the position of the object seen, the line of vision
    # and the distance between the object and the optical center:
TheObjectCenter = aNew3Dvector(0.0, 0.0, -0.15) # position of the center of the object seen.
TheDistance = 20.0 # ŧ distance between the object center and the optical center
    # The Zoom = 1.0 means that the image is neither reduced nor enlarged.
TheZoom = 2.3
    # Psi is the angle (in degrees) of rotation of the image around the line of vision
Psi = 0.0

# The line is a solution to a Chua system of ordinary differential equations (ODE)
TheNumberOfPoints = 2070 # number of points which determine the line
a = 15.568 # a and b are parameters in the Chua system
b = 28.0
TheInitialPosition = aNew3Dvector(-1.1752418506281723, 0.18422010169530834, 1.175655044221462)
dt = 5E-5 # the time between two steps of computation
TheNumberOfSteps = 100 # the number of steps between two points on the line
TheLine = aNewChuaSolution(TheNumberOfPoints, a, b, TheInitialPosition, dt, TheNumberOfSteps)

# The paintbrush (read the commentary at the top of LineDrawer.py)
TheHaloHalfWidth = 0.007
TheFringeWidth = 0.007
ThePowerOfSquaredCos = 1
ThePaintbrushColor = (255,215, 0) # (red, green, blue) : here gold
TheDrawer.takesAnewPaintbrush(TheHaloHalfWidth, TheFringeWidth, ThePowerOfSquaredCos)
TheDrawer.choosesAnewPaintbrushColor(ThePaintbrushColor)

# The canvas
TheCanvasStyle= "color" # "color", "black on white" or "white on black"
TheCanvasColor = (0, 0, 0) # (red, green, blue) : here black
TheDrawer.choosesAnewCanvasColor(TheCanvasColor)

# The animation
for i in range(360):
    TheDrawer.takesAnewCanvas(TheWidth_inNumberOfPixels, TheHeight_inNumberOfPixels,
                              TheCanvasStyle)
    Phi=120+i
    TheDrawer.choosesApointOfView(TheObjectCenter, TheDistance, Phi, Theta, Psi, TheZoom)
    TheDrawer.drawsAline(TheLine)
    TheDrawer.givesHisDrawing("Chua"+str(100+i)+".bmp")

print("Good bye")

LineDrawer.py

edit
# The line is black surrounded by a white halo (or the converse), so that when a segment is in front of others,
# we see it in the foreground (the line can also be in any color on any background color).
# The halo is itself surrounded by a fringe, so that the border of the halo is without crenelations.
# Inside the halo, the line is drawn more or less thinly. A pair power of cos(pi*(distance to the line)/haloHalfWidth))
# determines this thinness. The higher the power of cos^2, the thinner the line relative to its halo.
# The fringe width and the maximum distance of the halo to the (center of) the line are parameters. Therefore we have
# three parameters: haloHalfWidth (the maximum distance of the halo to the center of the line), fringeWidth, and
# powerOfSquaredCos.
# The grey shade of a pixel is a real number between zero and one. Zero means white and one, black, if the line is in
# black on a white background.
# The depth matrix associates with each pixel an ordered list of segments of the line which project on it. The first
# segment of the list is the nearest to the optical center.
# For the nearest segment of each pixel, the drawer determines if the center of the pixel is inside the halo, or inside
# a fringe of the halo. If it is inside the halo, the drawer considers only those segments which are near the first one.
# Those which are further than haloHalfWidth are ignored. Then it calculates a weighted average of the grey shades
# (calculated with a power of squared cos) of these segments. The weight of each segment decreases linearly with its
# distance to the nearest segment. The weighted average is the grey of the halo at the center of the pixel.
# If the center of the pixel is inside a fringe of the halo, the drawer considers the next segment on the list of the
# pixel. If there is none, the pixel is white. If there is one, the drawer calculates its shade as if it was a nearest
# segment. This shade is then considered as the shade of the outer border of the fringe. The fringe of its inner border
# (the outer border of the halo) is a weighted average of all the shades behind the first segment and very near.
# The shade of the pixel in the fringe is a weighted average of the two previous shades.

from PIL import Image
import math
from Vector2D import a2Dvector
from Vector3D import a3Dvector
from DepthMatrix import aDepthMatrix, aNewDepthMatrix, Order

class aCanvas:
    width_inNumberOfPixels = 1000
    height_inNumberOfPixels = 500
    realWidth = 1.0
    realHeight = 0.5
    style = "color" # "color", "black on white" or "white on black"
    color = (255, 255, 255)  # white in RGB

    def convertedToPixelCoordinates(self,aPoint):
        i = math.floor((aPoint.x + 0.5 * self.realWidth) * self.width_inNumberOfPixels / self.realWidth)
        j = math.floor((aPoint.y + 0.5 * self.realHeight) * self.height_inNumberOfPixels / self.realHeight)
        j = self.height_inNumberOfPixels - j - 1
        return [i,j]
 
    def convertedToRealCoordinates(self,i,j):
        aPoint = a2Dvector()
        aPoint.x = (i+0.5)*self.realWidth/self.width_inNumberOfPixels-0.5*self.realWidth
        aPoint.y = 0.5*self.realHeight-((j+0.5)*self.realHeight/self.height_inNumberOfPixels)
        return aPoint
        
class aPaintbrush:
    haloHalfWidth = 0.003 # 3 pixels if the width in pixels of the image is 1000, and its real width is 1.
    fringeWidth = 0.003 # the same
    powerOfSquaredCos = 1 # cos^2 determine the width of the line relative to the width of its halo
    color = (0, 0, 0) # Red, Green, Blue

class aLineDrawer:
    canvas = aCanvas()
    paintbrush = aPaintbrush()
    drawing = Image.new("L", (canvas.width_inNumberOfPixels, canvas.height_inNumberOfPixels), 255)
    pixel = drawing.load()
    # The line of vision:
    # phi is the angle between the projection of line of vision on the xy plane and the x axis.
    # theta is the angle of the line of vision relative to the z axis.
    phi = math.pi / 2  # the line of vision is parallel to the y axis
    theta = math.pi / 2  # the line of vision is horizontal
    # The line of vision is directed towards the object center
    objectCenter = a3Dvector()  # (0, 0, 0)
    # psi is the angle of rotation of the image around the line of vision.
    psi = 0.0
    # The position of the optical center is determined by the object center, the line of vision
    # and the distance between the object center and the optical center:
    opticalCenter = a3Dvector()
    distance = 10.0 # the distance between the object center and the optical center.
    # zoom = 1 means that the image is neither reduced nor enlarged.
    zoom = 1.0
    # depthMatrix is the memory of the perception of the depth of the object seen (its
    # distance relative to the observer).
    # Read the commentary in DepthMatrix.py for more explanations.
    depthMatrix = aDepthMatrix()
    # The following parameters are used to calculate the shade of the color painted by the drawer paintbrush:
    diffRed = paintbrush.color[0] - canvas.color[0]
    diffGreen = paintbrush.color[1] - canvas.color[1]
    diffBlue = paintbrush.color[2] - canvas.color[2]

    def takesAnewCanvas(self, TheWidth_inNumberOfPixels, TheHeight_inNumberOfPixels, style):
        self.canvas.width_inNumberOfPixels = TheWidth_inNumberOfPixels
        self.canvas.height_inNumberOfPixels = TheHeight_inNumberOfPixels
        self.canvas.realWidth = 1.0
        self.canvas.realHeight = self.canvas.realWidth * self.canvas.height_inNumberOfPixels / self.canvas.width_inNumberOfPixels
        self.canvas.style = style
        if style == "black on white":
            self.drawing = Image.new("L", (self.canvas.width_inNumberOfPixels, self.canvas.height_inNumberOfPixels), 255)
        if style == "white on black":
            self.drawing = Image.new("L", (self.canvas.width_inNumberOfPixels, self.canvas.height_inNumberOfPixels), 0)
        if style == "color":
            self.drawing = Image.new("RGB", (self.canvas.width_inNumberOfPixels,
                                        self.canvas.height_inNumberOfPixels), self.canvas.color)
        self.pixel = self.drawing.load()

    def choosesAnewCanvasColor(self, TheCanvasColor):
        self.canvas.color = TheCanvasColor
        diffRed = self.paintbrush.color[0] - self.canvas.color[0]
        diffGreen = self.paintbrush.color[1] - self.canvas.color[1]
        diffBlue = self.paintbrush.color[2] - self.canvas.color[2]

    def takesAnewPaintbrush(self, haloHalfWidth, fringeWidth, powerOfSquaredCos):
        self.paintbrush.halohalfWidth = haloHalfWidth
        self.paintbrush.fringeWidth = fringeWidth
        self.paintbrush.powerOfSquaredCos = powerOfSquaredCos

    def choosesAnewPaintbrushColor(self, ThePaintbrushColor):
        self.paintbrush.color = ThePaintbrushColor
        diffRed = self.paintbrush.color[0] - self.canvas.color[0]
        diffGreen = self.paintbrush.color[1] - self.canvas.color[1]
        diffBlue = self.paintbrush.color[2] - self.canvas.color[2]

    def choosesApointOfView(self, TheObjectCenter, TheDistance, Phi, Theta, Psi, zoom):
        # Phi, Theta, and Psi are in degrees not in radians.
        self.objectCenter = TheObjectCenter
        self.distance=TheDistance
        self.phi = Phi*math.pi/180 # phi is Phi in radians
        self.theta = Theta*math.pi/180
        self.psi = Psi*math.pi/180
        self.opticalCenter.x = 0.0
        self.opticalCenter.y = -self.distance
        self.opticalCenter.z = 0.0
        self.opticalCenter = self.opticalCenter.XaxisRotated(math.pi/2 - self.theta)
        self.opticalCenter = self.opticalCenter.ZaxisRotated(self.phi - math.pi/2)
        self.opticalCenter = self.opticalCenter.plus(self.objectCenter)
        self.zoom = zoom
        self.depthMatrix = aNewDepthMatrix(self.canvas.width_inNumberOfPixels, self.canvas.height_inNumberOfPixels)

    def givesHisDrawing(self,TheFileName):
        self.drawing.save(TheFileName)

    def pointImageOf(self,aPoint):
        p = a2Dvector()
        w = a3Dvector()
        w = aPoint.minus(self.opticalCenter)
        w = w.ZaxisRotated(math.pi/2-self.phi)
        w = w.XaxisRotated(self.theta-math.pi/2)
        w = w.YaxisRotated(self.psi)
        p = w.centeredProjectedOnXZ()
        p = p.times(self.zoom)
        return p

    # perceivesTheDepthOfAsegment fills the depth matrix for a segment.
    def perceivesTheDepthOfAsegment(self,TheFirstPoint, TheSecondPoint, TheFirstPointNumber):
        TheTotalHalfWidth = self.paintbrush.haloHalfWidth + self.paintbrush.fringeWidth
        imPt1 = self.pointImageOf(TheFirstPoint)
        imPt2 = self.pointImageOf(TheSecondPoint)
        # The drawing window is determined by its top left and bottom right corners.
        TopLeft = a2Dvector()
        BottomRight = a2Dvector()
        if imPt1.x < imPt2.x:
            TopLeft.x = imPt1.x
            BottomRight.x = imPt2.x
        else:
            TopLeft.x = imPt2.x
            BottomRight.x = imPt1.x
        if imPt1.y < imPt2.y:
            TopLeft.y = imPt1.y
            BottomRight.y = imPt2.y
        else:
            TopLeft.y = imPt2.y
            BottomRight.y = imPt1.y
        P1 = self.canvas.convertedToPixelCoordinates(TopLeft)
        P2 = self.canvas.convertedToPixelCoordinates(BottomRight)
        TheTotalHalfWidth_inPixels = math.floor(
            TheTotalHalfWidth * self.canvas.width_inNumberOfPixels / self.canvas.realWidth)
        # The window is enlarged :
        left = P1[0] - TheTotalHalfWidth_inPixels
        top = P2[1] - TheTotalHalfWidth_inPixels
        right = P2[0] + TheTotalHalfWidth_inPixels
        bottom = P1[1] + TheTotalHalfWidth_inPixels
        for i in range(left, right + 1):
            if i >= 0 and i < self.canvas.width_inNumberOfPixels:
                for j in range(top, bottom + 1):
                    if j >= 0 and j < self.canvas.height_inNumberOfPixels:
                        # The center of a pixel in real coordinates:
                        ThePixelCenter = self.canvas.convertedToRealCoordinates(i,j)
                        # The distance between the pixel center and the first point:
                        d1 = ThePixelCenter.minus(imPt1).norm()
                        # The distance between the pixel center and the second point:
                        d2 = ThePixelCenter.minus(imPt2).norm()
                        proj = ThePixelCenter.orthogonallyProjectedOnTheLine(imPt1, imPt2)
                        if (proj.x>=TopLeft.x and proj.x<=BottomRight.x and proj.y>=TopLeft.y and proj.y<=BottomRight.y
                                and d2>=TheTotalHalfWidth or d1<=TheTotalHalfWidth):
                            # dist is the distance between the center of the pixel and the line:
                            dist = ThePixelCenter.minus(proj).norm()
                            if dist < TheTotalHalfWidth:
                                self.depthMatrix.numberOfSegments[i][j] += 1
                                self.depthMatrix.segment[i][j] += [TheFirstPointNumber]

    def perceivesTheDepthOfAline(self, TheLine):
        # The depth matrix is filled but not ordered:
        for i in range(TheLine.numberOfPoints - 2):
            self.perceivesTheDepthOfAsegment(TheLine.point[i], TheLine.point[i+1], i)
        # The depth matrix is ordered:
        for i in range(self.canvas.width_inNumberOfPixels):
            for j in range(self.canvas.height_inNumberOfPixels):
                numberOfSegments = self.depthMatrix.numberOfSegments[i][j]
                depthList = [0.0 for i in range(numberOfSegments)]
                pointNumberList = [0 for i in range(numberOfSegments)]
                for k in range(numberOfSegments):
                    pointNumber = self.depthMatrix.segment[i][j][k]
                    pointNumberList[k] =pointNumber
                    TheFirstPoint = TheLine.point[pointNumber]
                    TheSecondPoint = TheLine.point[self.depthMatrix.segment[i][j][k] +1]
                    depthList[k] = TheFirstPoint.plus(TheSecondPoint.minus(TheFirstPoint).times(0.5)).minus(
                        self.opticalCenter).norm()
                Order(depthList, self.depthMatrix.segment[i][j], numberOfSegments)

    def drawsAline(self, TheLine):
        self.perceivesTheDepthOfAline(TheLine)
        for i in range(self.canvas.width_inNumberOfPixels):
            for j in range(self.canvas.height_inNumberOfPixels):
                shade = 0.0
                if self.depthMatrix.numberOfSegments[i][j] > 0:
                    ThePixelCenter = self.canvas.convertedToRealCoordinates(i, j)
                    shade = self.shade(i, j, ThePixelCenter, 0, TheLine)
                if self.canvas.style == "black on white":
                    intShade = 255 - math.floor(shade*255)
                    self.pixel[i, j] = intShade
                if self.canvas.style == "white on black":
                    intShade = math.floor(shade*255)
                    self.pixel[i, j] = intShade
                if self.canvas.style == "color":
                    intRedShade = math.floor(self.canvas.color[0] + shade*self.diffRed)
                    intGreenShade = math.floor(self.canvas.color[1] + shade*self.diffGreen)
                    intBlueShade = math.floor(self.canvas.color[2] + shade*self.diffBlue)
                    self.pixel[i, j] = [intRedShade, intGreenShade, intBlueShade]

    # The shade (a real number between 0 and 1) of the segment (i,j, depthRank) in the depth
    # matrix, calculated as if it was a nearest segment:
    def shade(self, i, j, ThePixelCenter, depthRank, TheLine):
        segmentNumber = self.depthMatrix.segment[i][j][depthRank]
        TheFirstPoint = TheLine.point[segmentNumber]
        TheSecondPoint = TheLine.point[segmentNumber+1]
        imPt1 = self.pointImageOf(TheFirstPoint)
        imPt2 = self.pointImageOf(TheSecondPoint)
        proj = ThePixelCenter.orthogonallyProjectedOnTheLine(imPt1, imPt2)
        # dist is the distance between the center of the pixel and the line:
        dist = ThePixelCenter.minus(proj).norm()
        depth = TheFirstPoint.plus(TheSecondPoint.minus(TheFirstPoint).times(0.5)).minus(self.opticalCenter).norm()
        if dist<self.paintbrush.haloHalfWidth:
            return self.shadeHalo(i, j, ThePixelCenter, depthRank, depth, dist, TheLine)
        else:
            return self.shadeFringe(i, j, ThePixelCenter, depthRank, depth, dist, TheLine)

    def shadeHalo(self, i, j, ThePixelCenter, depthRank, depth, dist, TheLine):
        shade = pow(math.cos(dist*math.pi/2/self.paintbrush.haloHalfWidth), 2*self.paintbrush.powerOfSquaredCos)
        shadeNumberOfShades= self.shadeNear(i,j, ThePixelCenter, depthRank, depth, TheLine)
        return shadeNumberOfShades[0]/(shadeNumberOfShades[1]+1)

    def shadeFringe(self, i, j, ThePixelCenter, depthRank, depth, dist, TheLine):
        if depthRank+1 < self.depthMatrix.numberOfSegments[i][j]:
            weight = (1.0 - ( (self.paintbrush.haloHalfWidth+self.paintbrush.fringeWidth)-dist)/self.paintbrush.fringeWidth)
            shade = weight*self.shade(i,j, ThePixelCenter, depthRank+1, TheLine)
            shade += (1.0-weight)*self.shadeNear(i,j, ThePixelCenter, depthRank, depth, TheLine)[0]
            return shade
        else:
            return 0.0

    # A weighted average of the shades of the segments very near and behind the segment (i,j,depthRank) in the
    # depth matrix. depth is the the depth of this segment:
    def shadeNear(self, i,j, ThePixelCenter, depthRank, depth, TheLine):
        out = 0
        shade = 0.0
        numberOfShades = 0
        while depthRank+1 < self.depthMatrix.numberOfSegments[i][j] and out < 1:
            depthRank +=1
            TheFirstPoint = TheLine.point[self.depthMatrix.segment[i][j][depthRank]]
            TheSecondPoint = TheLine.point[self.depthMatrix.segment[i][j][depthRank] + 1]
            depthNext = TheFirstPoint.plus(TheSecondPoint.minus(TheFirstPoint).times(0.5)).minus(
                self.opticalCenter).norm()
            if depthNext < depth + self.paintbrush.haloHalfWidth:
                numberOfShades += 1
                shade += (1.0-(depthNext-depth)/self.paintbrush.haloHalfWidth)*self.shade(i, j, ThePixelCenter, depthRank, TheLine)
            else:
                out = 1
        if numberOfShades>0:
            return [shade, numberOfShades]
        else:
            return [0.0, 0]

DepthMatrix.py

edit
# For the drawer to see depth, he needs to associate to each pixel the ordered (according to proximity)
# list of segments which are projected on it. The depth matrix associates to each pixel this ordered list of segments.
import math

class aDepthMatrix:
	numberOfSegments = []
	segment = []

def aNewDepthMatrix(TheWidth, TheHeight):
	TheMatrix = aDepthMatrix()
	TheMatrix.numberOfSegments = [[] for i in range(TheWidth)]
	TheMatrix.segment = [[] for i in range(TheWidth)]
	TheMatrix.dist = [[] for i in range(TheWidth)]
	for i in range(TheWidth):
		TheMatrix.numberOfSegments[i] = [0 for j in range(TheHeight)]
		TheMatrix.segment[i] = [[] for j in range(TheHeight)]
	return TheMatrix

# The following method orders the first list in its natural order and the second one
# according to the order of the first one. n is the number of items in both list
def Order(list1, list2, n):
	if n == 2:
		if list1[1] < list1[0]:
			storage = list1[0]
			list1[0] = list1[1]
			list1[1] = storage
	if n > 2:
		n1 = n // 2
		n2 = n - n1
		l11 = [0.0 for i in range(n1)]
		l21 = [0 for i in range(n1)]
		for i in range(n1):
			l11[i] = list1[i]
			l21[i] = list2[i]
		l12 = [0.0 for i in range(n2)]
		l22 = [0 for i in range(n2)]
		for i in range(n1):
			l12[i] = list1[n1 + i]
			l22[i] = list2[n1 + i]
		Order(l11, l21, n1)
		Order(l12, l22, n2)
		i1 = 0
		i2 = 0
		i3 = 0
		while i3 < n:
			if i1 == n1:
				list1[i3] = l12[i2]
				list2[i3] = l22[i2]
				i2 += 1
			else:
				if i2 == n2:
					list1[i3] = l11[i1]
					list2[i3] = l21[i1]
					i1 += 1
				else:
					if l11[i1] < l12[i2]:
						list1[i3] = l11[i1]
						list2[i3] = l21[i1]
						i1 += 1
					else:
						list1[i3] = l12[i2]
						list2[i3] = l22[i2]
						i2 += 1
			i3 += 1

ODE.py

edit
# ODE means Ordinary Differential Equation

from Vector3D import a3Dvector
from Line import aNewLine

# Runge-Kutta method
# A flow is determined by a system of ordinary differential equations.
# For example flow(x ,y ,z) = (dx/dt, dy/dt, dz/dt).
# The Leibniz method calculates an approximation of x(t+dt) with:
# x(t+dt) = x(t) + f(x)*dt where f = dx/dt is the flow.
# The Runge-Kutta method is a much better approximation:
# x(t+dt) = x(t) + (k1 + 2*k2 +2*k3 + k4)*dt/6 where
# k1 = f(x(t))
# k2 = f(x(t) + dt*k1/2)
# k3 = f(x(t) + dt*k2/2)
# k4 = f(x(t) + dt*k3)
def TheNextPosition(TheInitialPosition, TheFlow, dt):
    k1 = TheFlow(TheInitialPosition)
    k2 = TheFlow(TheInitialPosition.plus(k1.times(dt/2)))
    k3 = TheFlow(TheInitialPosition.plus(k2.times(dt/2)))
    k4 = TheFlow(TheInitialPosition.plus(k3.times(dt)))
    k5 = k2.times(2).plus(k3.times(2))
    k6 = k1.plus(k5).plus(k4).times(dt/6)
    return TheInitialPosition.plus(k6)

# An ODE solution for a 3D flow is a 3D line.
# TheNumberOfPoints is the number of points of the line.
# TheNumberOfSteps is the number of dt between two points.
def aNewODEsolution(TheNumberOfPoints, TheFlow, TheInitialPosition, dt, TheNumberOfSteps):
    TheSolution = aNewLine(TheNumberOfPoints)
    TheSolution.point[0] = TheInitialPosition
    TheCurrentPosition = TheInitialPosition
    for i in range(TheNumberOfPoints-1):
        for j in range(TheNumberOfSteps):
            TheCurrentPosition = TheNextPosition(TheCurrentPosition, TheFlow, dt)
        TheSolution.point[i+1] = TheCurrentPosition
    return TheSolution

#The Chua system (first version)
# dx/dt = a(y - x^3/16 + x/6)
# dy/dt = x - y + z
# dz/dt = -by
# For example b=28 and a =15.4
# Chua(pos, a, b) is the speed vector v of a moving point at position pos in a Chua system with a and b as parameters:
def Chua(pos, a, b):
    v = a3Dvector()
    v.x = a*(pos.y - pos.x*pos.x*pos.x/16 + pos.x/6)
    v.y = pos.x - pos.y + pos.z
    v.z = -b*pos.y
    return v

def aNewChuaSolution(TheNumberOfPoints, a, b, TheInitialPosition, dt, TheNumberOfSteps):
    def TheChuaFlow(pos):
        return Chua(pos, a, b)
    return aNewODEsolution(TheNumberOfPoints, TheChuaFlow, TheInitialPosition, dt, TheNumberOfSteps)

References:

Differential equations, dynamical systems and an introduction to chaos, Hirsch, Morris William., Smale, Stephen, Devaney, Robert Luke (2004)

Chapitres supplémentaires de la théorie des équations différentielles ordinaires, V. Arnold (1980)

Line.py

edit
from Vector3D import a3Dvector
import math

class aLine:
    numberOfPoints = 0  # the number of successive points on the line
    point = []

def aNewLine(TheNumberOfPoints):
    line = aLine()
    line.numberOfPoints=TheNumberOfPoints
    line.point = [a3Dvector() for i in range(TheNumberOfPoints)]
    return line

class Lines:
    numberOfLines = 0
    line = []
    lineColor = []

Vector3D.py

edit
import math
from Vector2D import a2Dvector

class a3Dvector:
    x = 0.0
    y = 0.0
    z = 0.0
    def plus(self,v):
        sum = a3Dvector()
        sum.x = self.x + v.x
        sum.y = self.y + v.y
        sum.z = self.z + v.z
        return sum

    def minus(self,v):
        diff = a3Dvector()
        diff.x = self.x - v.x
        diff.y = self.y - v.y
        diff.z = self.z - v.z
        return diff

    def times(self,a):
        prod = a3Dvector()
        prod.x = a*self.x
        prod.y = a*self.y
        prod.z = a*self.z
        return prod

    def norm(self):
        return math.sqrt(self.x*self.x + self.y*self.y + self.z*self.z)

    def XaxisRotated(self,theta):
        v = a3Dvector()
        costh = math.cos(theta)
        sinth = math.sin(theta)
        v.x = self.x
        v.y = costh*self.y - sinth*self.z
        v.z = sinth*self.y + costh*self.z
        return v

    def YaxisRotated(self,theta):
        v = a3Dvector()
        costh = math.cos(theta)
        sinth = math.sin(theta)
        v.y = self.y
        v.x = costh*self.x + sinth*self.z
        v.z = -sinth*self.x + costh*self.z
        return v

    def ZaxisRotated(self,theta):
        v = a3Dvector()
        costh = math.cos(theta)
        sinth = math.sin(theta)
        v.z = self.z
        v.x = costh*self.x - sinth*self.y
        v.y = sinth*self.x + costh*self.y
        return v

    # Consider the projection of the point v = (x,y,z) on the plane XZ, y = -1 followed
    # by a half turn rotation around the y axis: aPoint.centeredProjectedOnXZ() is the
    # image of aPoint by such a projection-rotation on the plane XZ, y = -1.
    # The image is first projected upside down and then rotated, like in the eye and the brain.
    def centeredProjectedOnXZ(self):
        p = a2Dvector()
        if abs(self.y) < 0.000001:
            p.x = 1000000.0
            p.y = 1000000.0
        else:
            p.x = self.x/self.y
            p.y = self.z/self.y
        return p

def aNew3Dvector(x,y,z):
    v = a3Dvector()
    v.x = x
    v.y = y
    v.z = z
    return v

Vector2D.py

edit
import math

class a2Dvector:
    x = 0.0
    y = 0.0
    def plus(self, The2Dvector):
        TheVectorSum = a2Dvector()
        TheVectorSum.x = self.x + The2Dvector.x
        TheVectorSum.y = self.y + The2Dvector.y
        return TheVectorSum
    def times(self, TheScalar):
        TheVectorTimesTheScalar = a2Dvector()
        TheVectorTimesTheScalar.x = TheScalar*self.x
        TheVectorTimesTheScalar.y = TheScalar*self.y
        return TheVectorTimesTheScalar

    def minus(self, The2Dvector):
        TheVectorDifference = a2Dvector()
        TheVectorDifference.x = self.x - The2Dvector.x
        TheVectorDifference.y = self.y - The2Dvector.y
        return TheVectorDifference

    def norm(self):
        return math.sqrt(self.x*self.x + self.y*self.y)

    def orthogonallyProjectedOnTheLine(self, TheFirstPoint, TheSecondPoint):
        TheProjectedPoint = a2Dvector()
        dx = TheSecondPoint.x - TheFirstPoint.x
        dy = TheSecondPoint.y - TheFirstPoint.y
        dx2 = dx*dx
        dy2 = dy*dy
        d2 = dx2 + dy2
        m11 = 0.0
        m22 = 0.0
        m21 = 0.0
        m12 = 0.0
        if d2>1E-200:
            m11 = dx2/d2
            m22 = dy2/d2
            m12 = dx*dy/d2
            m21 = m12
        TheProjectedPoint.x = TheFirstPoint.x + m11*(self.x - TheFirstPoint.x) + m12*(self.y - TheFirstPoint.y)
        TheProjectedPoint.y = TheFirstPoint.y + m21*(self.x - TheFirstPoint.x) + m22*(self.y - TheFirstPoint.y)
        return TheProjectedPoint

def aNew2Dvector(x, y):
    The2Dvector = a2Dvector()
    The2Dvector.x = x
    The2Dvector.y = y
    return The2Dvector