User:Thierry Dugnolle/Python/Mathematical painter

From Wikipedia, the free encyclopedia
A Chua attractor


First lesson in mathematical painting[edit]

The canvas and the palette[edit]

A mathematical painter draws and paints on electronic screens.

A canvas is a matrix of pixels. Each pixel is identified by its position (i, j) on the canvas, where i is the column number, j the row number. By convention (0, 0) is the left top pixel, and (Width-1, Height-1) the right bottom pixel, where Width is the number of pixels in width and Height the number of pixels in height.

A mathematical image is defined by assigning a brightness and a color to each pixel.

The brightness of a pixel is defined by an integer between 0 and 255. In black and white, 0 is black, 255 is white, 254 shades of gray separate them. 256 = 2 exponent 8. This is the number of numbers that can be encoded on one byte = 8 bits.

The color of a pixel is defined by a triplet of integers between 0 and 255: (red, green, blue).

(0, 0, 0) is black, (255, 255, 255) white, (255, 0, 0) brightest red, (0, 255, 0) brightest green, (0, 0, 255) brightest blue, (255, 255, 0) brightest yellow, (0, 255, 255) brightest cyan, and (255, 0, 255) brightest purple.

The color palette has 256*256*256 = 16,777,216 shades of color and brightness. A mathematical painter works with a palette of 16M colors. For computers, 1K = 210 = 1024 and 1M = 220 = 1,048,576. M is for Mega.

The palette can be represented by a cube. Black and White are two opposite vertices of the cube. Red, Green, Blue, Yellow, Cyan and Purple are the other six vertices.

The six faces of the cube of colors


Inside the cube of colors


The palette can also be represented by a double cone:

The two sides of the color double cone


Inside the color double cone


The axis of the color double cone is the line of grey shades (like the diagonal of the color cube) from white to black. Each point inside the double cone represents a color, defined by three coordinates, brightness, grey deviation and rainbow color, instead of red, green and blue in the color cube.

The canvas is a discrete field, because each of its pixels is identified by a pair of integers: (column number, row number). It is a field with discrete values, because each pixel is assigned a triplet of integers: (red, green, blue).

The paintbrush[edit]

A paintbrush draws continuous and differentiable (or rectifiable, this means that lines always have a tangent) lines except at a finite number of points: a line can be broken and have peaks.

A mathematical canvas is a discrete field with discrete values. A paintbrush must draw a line, therefore a continuous reality, in a discrete space.

If a line is one pixel wide, it is always with crenelations. To draw a high definition line without crenelation, we have to draw it a few pixels wide (at least 3) and to surround it by a halo:

The above triangle has been drawn as a succession of three corners, like the following one:

Any line can be drawn as a broken line, therefore as a succession of corners. To make rounded corners, a broken line must be divided at the middle of each of its segments.

An example of a broken line with 98 corners (and 100 points):

When a line comes from the back to the front of another line, there is a kind of jump from "line 1 is behind line 2" to "line 1 is in front of line 2". Animations are better with continuous transitions, like the following one:

Transparent and opaque paint[edit]

The paint can be more or less transparent, less or more opaque:

The vector image[edit]

A line is defined by a succession of a finite number of points. A point of the vector image is defined by its two coordinates in a plane (two coordinates is a vector). The vector image is the key to high definition computer graphics, because it is like a very high definition image: the position of each point in the vector image is defined with twice 15 digits or more.

A point of a line in space is defined by three coordinates. It is first projected on the plane of the camera or of the eye of the drawer, considered as a vector image, and then again projected on the canvas.

In a plane, the x axis is usually the left right axis and the y axis, the down up one. In the space, the x axis is still the left right axis, but the z axis is the down up one, and with a right handed frame, the y axis is the near far axis. The more y is positive, the farther.

The perception of depth[edit]

A 3D painter perceives the depth of space. To each point of a 3D line is assigned its depth, that is, its distance to the optical center. All the corners of all the lines which project on the canvas are identified in a depth matrix. In such a matrix, to each pixel of the canvas is assigned an array of corners which project on this pixel. The array is ordered according to the depth of the corners. The depth matrix is like an ordered memory of the depth of all the parts of all lines (or other objects) which project on the canvas.

Examples of main.py files[edit]

And God said, Let there be light: and there was light. (Genesis 1,4)

In the beginning was the Word (John 1,1)


A main.py file is like a script. It gives all the instructions to the painter to draw and paint an image or a succession of images.

High definition triangle

Continuous line crossing transition

A Chua line:


Cubic Galaxy:


A very little grain of salt:


Diffusion in one dimension:


Linear palettes[edit]

A linear palette assigns a color to each real number, usually between 0 and 1 (or -1 and 1). A linear palette is like a line in the cube of colors, the fundamental palette. The simplest line is the straight line between two colors. It is a gradation from the first color to the second one.

Linear palettes


Gradation from black to white

.

.

Gradation from gold to black

.

.

Heat colors

.

.

Depth colors

.

.

Psychedelic colors

.

.

More psychedelic colors

.

.

Periodic rainbow flag

.

.

Periodic bright rainbow

.

.

Periodic yellow green cyan

.

.


A periodic linear palette is like a loop in the cube of colors, a line without extremities.

Computer graphics is polycosmoscopy[edit]

Physics is monocosmic, mathematics is polycosmic.

Physics is a science of our universe, the Universe, the visible and tangible world in which we live. It is easy to prove that this universe is unique: two beings who act on each other are part of the same universe, by definition of a universe. A universe contains all beings that act on each other. But an observation always results from an action of an observed object on an observer. Everything that can be observed by a being from a universe is therefore also a being from this universe, which proves that our observable universe is unique. A being that is part of another universe cannot be observed from our universe.

Mathematics is a science of all theoretical beings. As soon as we can think of a being, it has a theoretical existence. That we can define it in words and reason about it is enough for it to exist. Theoretical beings are possibilities, their reality is only virtual. Their existence as a possibility is both more and less than existence in our universe, because theoretical existence is unaffected by the passage of time. Everything that exists today in our universe will eventually disappear, or almost everything, while a theoretically possible being is eternally possible.

Computers are universal machines. They can make all possible calculations on all beings that can be defined with a finite amount of information. Now all mathematical beings can be defined with a finite quantity of information, or be approached by an infinite series of steps which require a finite quantity of information. So computers can do all possible calculations on all mathematical beings. Computer graphics is polycosmoscopy, because it enables us to study all theoretically possible universes.

Reference: the cosmoscope is a theoretical invention of David Chalmers in Constructing The World (2012). Polycosmoscopy is a generalization.

See also[edit]

and

both in

The kernel[edit]

Painter.py[edit]

# 2023, October 11. Version 231011.2

# The class of painters is the fundamental class of all painters and drawers in this program.
# 2D and 3D lines drawers, 1D and 2D scalar fields painters and drawers are all
# members of this class.
# A painter has a canvas which is defined by a matrix width*height of pixels.
# The canvas is projected on a plane: the vector image.
# The image center is the position of the center of the canvas in the vector image.
# The length unit is the number of pixels of a unit of length in the vector image.
# A painter can take a new canvas, and gives his painting. All the work between
# requires more specific skills.
# A painter has also another canvas, where the color of each pixel is defined
# by a real number. This second canvas is like a matrix of grey shades, where there is an
# infinite number of shades.
# A painter is endowed with a few linear palettes. A linear palette converts a real number
# (between 0 and 1, or between -1 and 1) into a color defined by its three
# components: (red, green, blue)

import math
from PIL import Image
from time import process_time
from Vector2D import a2Dvector
from Palette import aPalette

class aCanvas:
    width = 1000 # number of pixels
    height = 500 # number of pixels
    imageCenter = a2Dvector()
    lengthUnit = 500 # number of pixels
    style = "black on white"  # "color", "black on white" or "white on black"
    color = (255, 255, 255)  # white in RGB
    drawing = Image.new("L", (width, height), 255)
    pixel = drawing.load()
    # pixelShade is a secondary canvas where the shade of a pixel is defined by a
    # real number between 0 and 1.
    pixelShade = []

    def convertedToPixelCoordinates(self, aPoint):
        i = math.floor((aPoint.x - self.imageCenter.x)*self.lengthUnit + 0.5*self.width-0.5)
        j = math.floor((aPoint.y - self.imageCenter.y)*self.lengthUnit + 0.5*self.height-0.5)
        j = self.height - j - 1
        return [i, j]

    def convertedToRealCoordinates(self, i, j):
        aPoint = a2Dvector()
        aPoint.x = (i + 0.5 - 0.5*self.width)/self.lengthUnit + self.imageCenter.x
        aPoint.y = (0.5*self.height - (j + 0.5))/self.lengthUnit + self.imageCenter.y
        return aPoint

class aPaintbrush:
    lineHalfWidth = 0.003
    color = (0, 0, 0) # (red, green, blue) are natural numbers between 0 and 255.
    # (0, 0, 0) is black. (255, 255, 255) is white. (255, 0, 0) is red.

class aPainter:
    canvas = aCanvas()
    palette = aPalette()
    paintbrush = aPaintbrush()

    def takesAnewCanvas(self, TheWidth, TheHeight, TheLengthUnit, style, color):
        self.canvas.width = TheWidth
        self.canvas.height = TheHeight
        self.canvas.lengthUnit = TheLengthUnit
        self.canvas.style = style
        if style == "black on white":
            self.canvas.color = (255, 255 ,255)
            self.paintbrush.color = (0, 0, 0)
            self.canvas.drawing = Image.new("L", (self.canvas.width, self.canvas.height),255)
        if style == "white on black":
            self.canvas.color = (0, 0, 0)
            self.paintbrush.color = (255, 255, 255)
            self.canvas.drawing = Image.new("L", (self.canvas.width, self.canvas.height), 0)
        if style == "color":
            self.canvas.color = color
            self.canvas.drawing = Image.new("RGB", (self.canvas.width, self.canvas.height), self.canvas.color)
        self.canvas.pixel = self.canvas.drawing.load()
        self.canvas.pixelShade = [[0.0 for j in range(self.canvas.height)] for i in range(self.canvas.width)]

    def drawsOnCanvas(self):
        print("Time:", process_time(), "The drawer draws on canvas, style", self.canvas.style)
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                shade = self.canvas.pixelShade[i][j]
                if shade < 1:
                    if shade < 0:
                        shade = 0
                else:
                    shade = 1
                if self.canvas.style == "black on white":
                    self.canvas.pixel[i, j] = 255 - math.floor(shade*255)
                if self.canvas.style == "white on black":
                    self.canvas.pixel[i, j] = math.floor(shade*255)
                if self.canvas.style == "color":
                    self.canvas.pixel[i, j] = self.palette.gradation(shade, self.canvas.color, self.paintbrush.color)

    def givesApainting(self,TheFileName):
        self.canvas.drawing.save(TheFileName)

Painter2D.py[edit]

# 2023, October 20. Version 231020

# A drawer always draws a broken line (a succession of a finite number of
# segments, defined by a succession of a finite number of points) corner
# after corner. A corner of a broken line is a succession of three points:
# the second point of the corner is a point i of the line, the first point
# is the middle between point i-1 and point i, the third point is the middle
# between point i and point i+1. The line is divided in corners and not in
# segments so that the drawing of rounded corners is easy.

import math
from Painter import aPainter
from Vector2D import The2Dvector, The2Dnorm

# A 2D line drawer has all the properties and methods of a painter.
class a2Dpainter(aPainter):

# cornerDrawingWindow returns the left top and right bottom corners (in pixel
# coordinates) of the drawing window of the line corner whose three points are
# linePoint1, linePoint2 and linePoint3 in the vector image.
    def cornerDrawingWindow(self, cornerPoint1, cornerPoint2, cornerPoint3, margin):
        leftTop = The2Dvector(min(cornerPoint1.x, cornerPoint2.x, cornerPoint3.x), max(cornerPoint1.y, cornerPoint2.y, cornerPoint3.y))
        rightBottom = The2Dvector(max(cornerPoint1.x, cornerPoint2.x, cornerPoint3.x), min(cornerPoint1.y, cornerPoint2.y, cornerPoint3.y))
        leftTopPixel = self.canvas.convertedToPixelCoordinates(leftTop)
        rightBottomPixel = self.canvas.convertedToPixelCoordinates(rightBottom)
        marginInPixels = math.ceil(margin*self.canvas.lengthUnit)
        left = leftTopPixel[0] - marginInPixels
        top = leftTopPixel[1] - marginInPixels
        right = rightBottomPixel[0] + marginInPixels + 1
        bottom = rightBottomPixel[1] + marginInPixels + 1
        if left < 0:
            left = 0
        if right > self.canvas.width - 1:
            right = self.canvas.width - 1
        if top < 0:
            top = 0
        if bottom > self.canvas.height - 1:
            bottom = self.canvas.height - 1
        return [left, top, right, bottom]

    def diskDrawingWindow(self, center, radius):
        centerPixel = self.canvas.convertedToPixelCoordinates(center)
        radiusInPixels = math.ceil(radius*self.canvas.lengthUnit)
        left = centerPixel[0] - radiusInPixels
        top = centerPixel[1] - radiusInPixels
        right = centerPixel[0] + radiusInPixels + 1
        bottom = centerPixel[1] + radiusInPixels + 1
        if left < 0:
            left = 0
        if right > self.canvas.width - 1:
            right = self.canvas.width - 1
        if top < 0:
            top = 0
        if bottom > self.canvas.height - 1:
            bottom = self.canvas.height - 1
        return [left, top, right, bottom]

    def distanceOfThePixelCenterFromAcorner(self, i, j, cornerPoint1, cornerPoint2, cornerPoint3):
        vF = cornerPoint2.minus(cornerPoint1)
        vS = cornerPoint3.minus(cornerPoint2)
        vF2 = vF.scalar(vF)
        vS2 = vS.scalar(vS)
        ThePixelCenter = self.canvas.convertedToRealCoordinates(i, j)  # The center of a pixel in real coordinates
        # Projection of the center of the pixel on the first segment of the corner:
        uF = ThePixelCenter.minus(cornerPoint1)
        pF = uF.scalar(vF)
        if vF2 > 1E-20:
            projF = cornerPoint1.plus(vF.times(pF/vF2))
        else:
            projF = cornerPoint2
        distF = The2Dnorm(ThePixelCenter.minus(projF))
        # Projection of the center of the pixel on the second segment of the corner:
        uS = ThePixelCenter.minus(cornerPoint2)
        pS = uS.scalar(vS)
        if vS2 > 1E-20:
            projS = cornerPoint2.plus(vS.times(pS/vS2))
        else:
            projS = cornerPoint2
        distS = The2Dnorm(ThePixelCenter.minus(projS))
        dist = 1E20
        if pF >= 0 and pF < vF2:
            dist = distF
        if pS >= 0 and pS < vS2:
            dist = distS
            if pF >= 0 and pF < vF2 and distF < distS:
                dist = distF
        if pF >= vF2 and pS < 0:
            dist = The2Dnorm(ThePixelCenter.minus(cornerPoint2))
        return dist

    def drawsA2Dcorner(self, cornerPoint1, cornerPoint2, cornerPoint3):
        fourIntegers = self.cornerDrawingWindow(cornerPoint1, cornerPoint2, cornerPoint3, self.paintbrush.lineHalfWidth)
        left = fourIntegers[0]
        top = fourIntegers[1]
        right = fourIntegers[2]
        bottom = fourIntegers[3]
        for i in range(left, right):
            for j in range(top, bottom):
                dist = self.distanceOfThePixelCenterFromAcorner(i, j, cornerPoint1, cornerPoint2, cornerPoint3)
                if dist < self.paintbrush.lineHalfWidth:
                    self.canvas.pixelShade[i][j] += pow(math.cos(dist*math.pi/2/self.paintbrush.lineHalfWidth), 2)

    def drawsA2Datom(self, center, radius):
        fourIntegers = self.diskDrawingWindow(center, radius)
        left = fourIntegers[0]
        top = fourIntegers[1]
        right = fourIntegers[2]
        bottom = fourIntegers[3]
        for i in range(left, right):
            for j in range(top, bottom):
                ThePixelCenter = self.canvas.convertedToRealCoordinates(i, j)
                dist = The2Dnorm(center.minus(ThePixelCenter))
                if dist < radius:
                    self.canvas.pixelShade[i][j] +=pow(math.cos(dist*math.pi/2/radius), 2)

    def drawsA2Dline(self, TheLine):
        for i in range(TheLine.numberOfPoints-2):
            v1 = TheLine.point[i+1].minus(TheLine.point[i]).times(0.5)
            v2 = TheLine.point[i+2].minus(TheLine.point[i+1]).times(0.5)
            cornerPoint1 = TheLine.point[i].plus(v1)
            cornerPoint2 = TheLine.point[i+1]
            cornerPoint3 = TheLine.point[i+1].plus(v2)
            self.drawsA2Dcorner(cornerPoint1, cornerPoint2, cornerPoint3)
        self.drawsOnCanvas()

Vector2D.py[edit]

# 2023, October 18. Version 231018

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 scalar(self, v):
        return self.x*v.x + self.y*v.y

def The2Dnorm(v):
    return math.sqrt(v.x*v.x + v.y*v.y)

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

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

def aNew2Dline(TheNumberOfPoints):
    line = a2Dline()
    line.numberOfPoints = TheNumberOfPoints
    line.point = [a2Dvector() for i in range(TheNumberOfPoints)]
    return line

Monocolor3Dpainter.py[edit]

# 2023, October 20. Version 231020

# The 3D brush:
# In 3D, a line can have a fringe, with which it hides and whitens the background (whitens
# if the canvas is white).
# Atoms or points are painted like dots. They have a width and they may have a fringe.
# A 3D brush is a 2D paintbrush, defined by a line half width and a color, completed
# by a line fringe relative width, a line depth, a line opacity, an atom fringe relative
# width, an atom opacity and an intensity. The line half width shall usually be chosen
# to a few pixels, for a thin line. The line depth is usually of
# the same order of magnitude than the line width.
# If the opacity = 0 the paint is transparent, and depth is not seen, if it
# is = 1/2 the paint is semi-transparent, if it is = 1 or more the paint is opaque.
# An image is first defined with an infinitely precise shade of color (of grey, if
# the line is black on a white canvas), a real number, usually the shade between 0 and 1.
# The intensity of the brush can be a useful parameter to lighten or darken an image. It
# multiplies the shade. If it is < 1, the intensity of the shade is reduced, if > 1 the
# intensity of the shade is augmented.

# 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.
# The line of vision is directed towards the object center.
# 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.

# The zoom = 1 means that the image is neither reduced nor enlarged, > 1 that it is enlarged,
# > 0 and < 1 that it is reduced.
# psi is the angle of rotation of the image around the line of vision.

# The depth matrix:
# The depth matrix is the memory of the perception of the depth of the object seen
# (its distance relative to the observer). A depth matrix associates to each pixel
# the ordered list, according to proximity, of line corners which are projected on
# it. A line corner is memorized in the depth matrix with the distance of its
# projected point to the pixel and its distance to the optical center (its depth).

from Vector2D import a2Dvector, The2Dnorm
from Vector3D import a3Dvector, The3Dnorm
from Painter2D import a2Dpainter
from DepthMatrix import aDepthMatrix, aNewDepthMatrix, Order2lists
import math
from time import process_time
import json

class aVision:
    phi = math.pi / 2  # The line of vision is parallel to the y axis.
    theta = math.pi / 2  # The line of vision is horizontal.
    objectCenter = a3Dvector()  # (0, 0, 0)
    opticalCenter = a3Dvector()
    distance = 10.0 # The distance between the object center and the optical center.
    zoom = 1.0
    psi = 0.0

class a3Dbrush:
    lineFringeRelativeWidth = 1
    lineDepth = 0.03
    lineOpacity = 2 # A real number > or = 0. opacity = 0 means transparency.
    atomFringeRelativeWidth = 0.1
    atomOpacity = 2
    atomStyle = "cloud" # or "bright"
    intensity = 1

# A 3D line drawer has all the properties and methods of 2D line drawer:
class aMonocolor3Dpainter(a2Dpainter):
    brush3D = a3Dbrush()
    vision = aVision
    depthMatrix = aDepthMatrix()

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

    # 'draws3Dlines' is a competence of the painter:
    def draws3Dlines(self, TheLines):
        print("Time:", process_time(), "The drawer fills the depth matrix with lines.")
        self.perceivesLines(TheLines)
        self.paintsAview()

    def drawsA3Dline(self, TheLine):
        print("Time:", process_time(), "The drawer fills the depth matrix with a line.")
        self.perceivesAline(TheLine)
        self.paintsAview()

    # 'perceives lines' fills the depth matrix with lines:
    def perceivesLines(self, TheLines):
        for i in range(TheLines.numberOfLines):
            if i % 50000 == 0:
                print("Time:", process_time(), "The drawer perceives line", i)
            self.perceivesAline(TheLines.line[i])

    # 'perceives a line' fills the depth matrix with a line:
    def perceivesAline(self, TheLine):
        for i in range(TheLine.numberOfPoints - 2):
            if i > 0 and i % 50000 == 0:
                print("Time:", process_time(), "The drawer perceives point", i)
            linePoint1 = TheLine.point[i]
            linePoint2 = TheLine.point[i + 1]
            linePoint3 = TheLine.point[i + 2]
            v1 = linePoint2.minus(linePoint1).times(0.5)
            v2 = linePoint3.minus(linePoint2).times(0.5)
            cornerPoint1 = linePoint1.plus(v1)
            cornerPoint2 = linePoint2
            cornerPoint3 = linePoint2.plus(v2)
            self.perceivesAcorner(cornerPoint1, cornerPoint2, cornerPoint3)

    # 'perceives a corner' fills the depth matrix with a corner:
    def perceivesAcorner(self, cornerPoint1, cornerPoint2, cornerPoint3):
        im1 = self.pointImageOf(cornerPoint1)
        im2 = self.pointImageOf(cornerPoint2)
        im3 = self.pointImageOf(cornerPoint3)
        if im1[1] == 1 and im2[1] == 1 and im3[1] == 1:
            cornerPoint1image = im1[0]
            cornerPoint2image = im2[0]
            cornerPoint3image = im3[0]
            fourIntegers = self.cornerDrawingWindow(cornerPoint1image, cornerPoint2image, cornerPoint3image,
                                        self.paintbrush.lineHalfWidth*(1+self.brush3D.lineFringeRelativeWidth))
            left = fourIntegers[0]
            top = fourIntegers[1]
            right = fourIntegers[2]
            bottom = fourIntegers[3]
            for i in range(left, right):
                for j in range(top, bottom):
                    vF = cornerPoint2image.minus(cornerPoint1image)
                    vS = cornerPoint3image.minus(cornerPoint2image)
                    vF2 = vF.scalar(vF)
                    vS2 = vS.scalar(vS)
                    ThePixelCenter = self.canvas.convertedToRealCoordinates(i,j)
                    # Projection of the center of the pixel on the first segment of the corner:
                    uF = ThePixelCenter.minus(cornerPoint1image)
                    pF = uF.scalar(vF)
                    if vF2 > 1E-20:
                        projF = cornerPoint1image.plus(vF.times(pF / vF2))
                    else:
                        projF = cornerPoint2image
                    distF = The2Dnorm(ThePixelCenter.minus(projF))
                    # Projection of the center of the pixel on the second segment of the corner:
                    uS = ThePixelCenter.minus(cornerPoint2image)
                    pS = uS.scalar(vS)
                    if vS2 > 1E-20:
                        projS = cornerPoint2image.plus(vS.times(pS / vS2))
                    else:
                        projS = cornerPoint2image
                    distS = The2Dnorm(ThePixelCenter.minus(projS))
                    dist = 1E20
                    cornerPoint4 = cornerPoint2
                    if pF >= 0 and pF < vF2:
                        dist = distF
                        if vF2 > 1E-20:
                            cornerPoint4 = cornerPoint1.plus(cornerPoint2.minus(cornerPoint1).times(pF / vF2))
                    if pS >= 0 and pS < vS2:
                        dist = distS
                        if vS2 > 1E-20:
                            cornerPoint4 = cornerPoint2.plus(cornerPoint3.minus(cornerPoint2).times(pS / vS2))
                        if pF >= 0 and pF < vF2 and distF < distS:
                            dist = distF
                            if vF2 > 1E-20:
                                cornerPoint4 = cornerPoint1.plus(cornerPoint2.minus(cornerPoint1).times(pF / vF2))
                    if pF >= vF2 and pS < 0:
                        dist = The2Dnorm(ThePixelCenter.minus(cornerPoint2image))
                    if dist < self.paintbrush.lineHalfWidth*(1+self.brush3D.lineFringeRelativeWidth):
                        self.depthMatrix.numberOfPoints[i][j] += 1
                        self.depthMatrix.dist[i][j] += [dist/self.paintbrush.lineHalfWidth]
                        self.depthMatrix.depth[i][j] += [The3Dnorm(cornerPoint4.minus(self.vision.opticalCenter))]

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

    # The atoms are of the same kind:
    def paintsAtoms(self, TheAtoms):
        print("Time:", process_time(), "The drawer fills the depth matrix with", TheAtoms.numberOfAtoms, "atoms.")
        self.perceivesAtoms(TheAtoms)
        self.paintsAview()

    def paintsBondedAtoms(self, TheBondedAtoms):
        print("Time:", process_time(), "The drawer fills the depth matrix with", TheBondedAtoms.atoms.numberOfAtoms,
              "atoms and", TheBondedAtoms.bonds.numberOfLines, "bonds.")
        self.perceivesAtoms(TheBondedAtoms.atoms)
        self.perceivesLines(TheBondedAtoms.bonds)
        self.paintsAview()

    # The atoms are of different kinds:
    def paintsAtomsOfDifferentKinds(self, TheAtoms):
        for i in range(TheAtoms.numberOfKinds):
            self.perceivesAtoms(TheAtoms.kind[i])
        self.paintsAview()

    def paintsBondedAtomsOfDifferentKinds(self, TheBondedAtoms):
        print("Time:", process_time(), "The drawer fills the depth matrix with atoms.")
        for i in range(TheBondedAtoms.atoms.numberOfKinds):
            self.perceivesAtoms(TheBondedAtoms.atoms.kind[i])
        print("Time:", process_time(), "The drawer fills the depth matrix with", TheBondedAtoms.bonds.numberOfLines, "bonds.")
        self.perceivesLines(TheBondedAtoms.bonds)
        self.paintsAview()


    def perceivesAtoms(self, TheAtoms):
        for i in range(TheAtoms.numberOfAtoms):
            if i % 500000 == 0:
                print("Time:", process_time(), "The drawer perceives atom", i)
            self.perceivesAnAtom(TheAtoms.position[i], TheAtoms.radius)

    def perceivesAnAtom(self, TheAtomCenter, TheAtomRadius):
        TheAtomCenterDepth = The3Dnorm(TheAtomCenter.minus(self.vision.opticalCenter))
        TheVisibleRadius = TheAtomRadius/TheAtomCenterDepth*self.vision.zoom
        if self.brush3D.atomStyle == "bright":
            TheVisibleRadius = max(TheVisibleRadius, self.paintbrush.lineHalfWidth*(self.brush3D.lineFringeRelativeWidth+1))
        im = self.pointImageOf(TheAtomCenter)
        if im[1] == 1:
            TheAtomImage = im[0]
            fourIntegers = self.diskDrawingWindow(TheAtomImage, TheVisibleRadius)
            left = fourIntegers[0]
            top = fourIntegers[1]
            right = fourIntegers[2]
            bottom = fourIntegers[3]
            for i in range(left, right):
                for j in range(top, bottom):
                    ThePixelCenter = self.canvas.convertedToRealCoordinates(i,j)
                    dist = The2Dnorm(TheAtomImage.minus(ThePixelCenter))
                    if dist < TheVisibleRadius:
                        self.depthMatrix.numberOfPoints[i][j] += 1
                        # The minus sign in front of 'dist' and the -1 distinguishes atomes from corners:
                        self.depthMatrix.dist[i][j] += [-dist/TheVisibleRadius - 1]
                        if self.brush3D.atomStyle == "bright":
                            depth = TheAtomCenterDepth - TheAtomRadius
                        else:
                            depth = TheAtomCenterDepth - TheAtomRadius*math.sqrt(pow(TheVisibleRadius, 2)
                                                                - dist*dist)/TheVisibleRadius
                        self.depthMatrix.depth[i][j] += [depth]

    def paintsAview(self):
        self.ordersTheDepthMatrix()
        self.paintsTheDepthMatrix()

    def ordersTheDepthMatrix(self):
        print("Time:", process_time(), "The drawer orders the depth matrix.")
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                numberOfPoints = self.depthMatrix.numberOfPoints[i][j]
                if numberOfPoints > 1:
                    orderDepth = Order2lists(self.depthMatrix.depth[i][j], self.depthMatrix.dist[i][j], numberOfPoints)
                    self.depthMatrix.depth[i][j] = orderDepth[0]
                    self.depthMatrix.dist[i][j] = orderDepth[1]

    def paintsTheDepthMatrix(self):
        print("Time:", process_time(), "The drawer paints the monocolor depth matrix.")
        max = 0.0
        for i in range(self.canvas.width):
            if i % 100 == 0:
                print("Time:", process_time(), "column", i)
            for j in range(self.canvas.height):
                if self.depthMatrix.numberOfPoints[i][j] > 0:
                    self.canvas.pixelShade[i][j] = self.shade(i, j)
                    if max < self.canvas.pixelShade[i][j]:
                        max = self.canvas.pixelShade[i][j]
        print("Time:", process_time(), "shade max =", max, "intensity =", self.brush3D.intensity)
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                self.canvas.pixelShade[i][j] *= self.brush3D.intensity
        self.drawsOnCanvas()

    def shade(self, i, j):
        numberOfPoints = self.depthMatrix.numberOfPoints[i][j]
        totalGiveColor = 0.0
        if numberOfPoints > 0:
            pointTakeColor = [1.0 for i in range(numberOfPoints)]
            totalTakeColor = 1.0
            totalTakeColorRank = -1
            depthRank1 = 0
            while depthRank1 < numberOfPoints:
                dist1 = self.depthMatrix.dist[i][j][depthRank1]
                if dist1 >= 0:
                    giveColor1 = self.giveColor(dist1)
                    pointTakeColor[depthRank1] = self.lineTakeColor(dist1)
                else:
                    giveColor1 = self.atomGiveColor(abs(dist1) - 1)
                    pointTakeColor[depthRank1] = self.atomTakeColor(abs(dist1) - 1)
                if depthRank1 > 0:
                    depth1 = self.depthMatrix.depth[i][j][depthRank1]
                    depthRank2 = depthRank1 - 1
                    depth2 = self.depthMatrix.depth[i][j][depthRank2]
                    while depthRank2 >= 0 and depth2 > depth1 - self.brush3D.lineDepth:
                        grade = (depth1 - depth2) / self.brush3D.lineDepth
                        giveColor1 *= 1 - grade * (1 - pointTakeColor[depthRank2])
                        depthRank2 -= 1
                        if depthRank2 >= 0:
                            depth2 = self.depthMatrix.depth[i][j][depthRank2]
                    totalTakeColorRankStorage = depthRank2
                    while depthRank2 >= 0 and depthRank2 > totalTakeColorRank:
                        totalTakeColor *= pointTakeColor[depthRank2]
                        depthRank2 -= 1
                    totalTakeColorRank = totalTakeColorRankStorage
                totalGiveColor += giveColor1 * totalTakeColor
                if totalTakeColor*(numberOfPoints - depthRank1)*self.brush3D.intensity < 1/255:
                    depthRank1 = numberOfPoints
                depthRank1 += 1
        return totalGiveColor

    # giveColor = 0 means no coloring, giveColor = 1 means full coloring,
    # because giveColor works additively.
    def giveColor(self, dist):
        giveColor = 0
        if dist < 1:
            giveColor = pow(math.cos(dist*math.pi/2), 2)
        return giveColor

    # takeColor = 1 means no change, takeColor = 0 means full color deleting,
    # because takeColor works multiplicatively.
    def atomGiveColor(self, dist):
        giveColor = 0
        if dist < 1:
            if self.brush3D.atomStyle == "bright":
                giveColor = 1
            else:
                giveColor = pow(math.cos(dist*math.pi/2), 2)
        return giveColor

    # takeColor = 1 means no change, takeColor = 0 means full color deleting,
    # because takeColor works multiplicatively.
    def lineTakeColor(self, dist):
        TakeColor = 1
        if dist < 1 + self.brush3D.lineFringeRelativeWidth:
            TakeColor = 1 - self.brush3D.lineOpacity*pow(math.cos(dist*math.pi/2/(1 + self.brush3D.lineFringeRelativeWidth)), 2)
            if TakeColor < 0:
                TakeColor = 0
        return TakeColor

    def atomTakeColor(self, dist):
        TakeColor = 1
        if dist < 1 + self.brush3D.atomFringeRelativeWidth:
            TakeColor = 1 - self.brush3D.atomOpacity*pow(math.cos(dist*math.pi/2/(1 + self.brush3D.atomFringeRelativeWidth)), 2)
            if TakeColor < 0:
                TakeColor = 0
        return TakeColor

    def paintsAtomsOfDifferentKindsSimple(self,TheAtoms):
        for k in range(TheAtoms.numberOfKinds):
            print("Time:", process_time(), "The drawer paints", TheAtoms.kind[k].numberOfAtoms, "atoms of kind", k)
            for i in range(TheAtoms.kind[k].numberOfAtoms):
                if i % 1000000 == 0:
                    print("Time:", process_time(), "The drawer paints atom", i)
                TheAtomCenterDepth = The3Dnorm(TheAtoms.kind[k].position[i].minus(self.vision.opticalCenter))
                TheVisibleRadius = TheAtoms.kind[k].radius / TheAtomCenterDepth * self.vision.zoom
                im = self.pointImageOf(TheAtoms.kind[k].position[i])
                if im[1] == 1:
                    TheAtomImage = im[0]
                    self.drawsA2Datom(TheAtomImage, TheVisibleRadius)
            max = 0
            for i in range(self.canvas.width):
                for j in range(self.canvas.height):
                    if max < self.canvas.pixelShade[i][j]:
                        max = self.canvas.pixelShade[i][j]
            print("Time:", process_time(), "shade max =", max, "intensity =", self.brush3D.intensity)
            for i in range(self.canvas.width):
                for j in range(self.canvas.height):
                    self.canvas.pixelShade[i][j] *= self.brush3D.intensity
        self.drawsOnCanvas()


    def paintsAtomsOfTheSameKindSimple(self, TheAtoms):
        print("Time:", process_time(), "The drawer paints", TheAtoms.numberOfAtoms, "atoms.")
        for i in range(TheAtoms.numberOfAtoms):
            if i % 1000000 == 0:
                print("Time:", process_time(), "The drawer paints atom", i)
            TheAtomCenterDepth = The3Dnorm(TheAtoms.position[i].minus(self.vision.opticalCenter))
            TheVisibleRadius = TheAtoms.radius/TheAtomCenterDepth*self.vision.zoom
            im = self.pointImageOf(TheAtoms.position[i])
            if im[1] == 1:
                TheAtomImage = im[0]
                self.drawsA2Datom(TheAtomImage, TheVisibleRadius)
        max = 0
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                if max < self.canvas.pixelShade[i][j]:
                    max = self.canvas.pixelShade[i][j]
        print("Time:", process_time(), "shade max =", max, "intensity =", self.brush3D.intensity)
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                self.canvas.pixelShade[i][j] *= self.brush3D.intensity
        self.drawsOnCanvas()

    def savesAperception(self, TheFileName):
        with open(TheFileName, 'w') as f:
            txtMat = []
            txtMat += [self.canvas.width]
            txtMat += [self.canvas.height]
            txtMat += [self.paintbrush.lineHalfWidth]
            txtMat += [self.brush3D.lineFringeRelativeWidth]
            txtMat += [self.brush3D.atomFringeRelativeWidth]
            for i in range(self.canvas.width):
                for j in range(self.canvas.height):
                    TheNumberOfPoints = self.depthMatrix.numberOfPoints[i][j]
                    txtMat += [TheNumberOfPoints]
                    for k in range(TheNumberOfPoints):
                        txtMat += [self.depthMatrix.dist[i][j][k]]
                        txtMat += [self.depthMatrix.depth[i][j][k]]
            json.dump(txtMat, f)

    def retrievesAmemory(self, TheFileName):
        with open(TheFileName, 'r') as f:
            txtMat = json.load(f)
            self.depthMatrix = aNewDepthMatrix(int(txtMat[0]), int(txtMat[1]))
            self.paintbrush.lineHalfWidth = float(txtMat[2])
            self.brush3D.lineFringeRelativeWidth = float(txtMat[3])
            self.brush3D.atomFringeRelativeWidth = float(txtMat[4])
            count = 5
            for i in range(self.depthMatrix.width):
                for j in range(self.depthMatrix.height):
                    TheNumberOfPoints = int(txtMat[count])
                    count += 1
                    self.depthMatrix.numberOfPoints[i][j] = TheNumberOfPoints
                    self.depthMatrix.dist[i][j] = [0.0 for i in range(TheNumberOfPoints)]
                    self.depthMatrix.depth[i][j] = [0.0 for i in range(TheNumberOfPoints)]
                    for k in range(TheNumberOfPoints):
                        self.depthMatrix.dist[i][j][k] = float(txtMat[count])
                        count += 1
                        self.depthMatrix.depth[i][j][k] = float(txtMat[count])
                        count += 1

Line3D.py[edit]

# 2023, October 16. Version 231016

import math
import json
from Vector3D import a3Dvector

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

    def plus(self, v):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].plus(v)
        return TheLine

    def minus(self, v):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].minus(v)
        return TheLine

    def times(self, a):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].times(a)
        return TheLine

    def XaxisRotated(self, theta):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].XaxisRotated(theta)
        return TheLine

    def YaxisRotated(self, theta):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].YaxisRotated(theta)
        return TheLine

    def ZaxisRotated(self, theta):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].ZaxisRotated(theta)
        return TheLine

    def Xdilated(self, a):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].Xdilated(a)
        return TheLine

    def Ydilated(self, a):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].Ydilated(a)
        return TheLine

    def Zdilated(self, a):
        TheLine = aNew3Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = self.point[i].Zdilated(a)
        return TheLine

    def isWrittenOnFile(self, TheFileName):
        with open(TheFileName, 'w') as f:
            txtLine = []
            txtLine += [self.numberOfPoints]
            for i in range(self.numberOfPoints):
                txtLine += [self.point[i].x, self.point[i].y, self.point[i].z]
            json.dump(txtLine, f)

def TheLineOnFile(TheFileName):
    with open(TheFileName, 'r') as f:
        txtLine = json.load(f)
        TheNumberOfPoints = int(txtLine[0])
        TheLine = aNew3Dline(TheNumberOfPoints)
        for i in range(TheLine.numberOfPoints):
            TheLine.point[i].x = float(txtLine[1 + 3*i])
            TheLine.point[i].y = float(txtLine[2 + 3*i])
            TheLine.point[i].z = float(txtLine[3 + 3*i])
    return TheLine

def aNew3Dline(TheNumberOfPoints):
    TheLine = a3Dline()
    TheLine.numberOfPoints = TheNumberOfPoints
    TheLine.point = [a3Dvector() for i in range(TheNumberOfPoints)]
    return TheLine

def TheStraight3Dline(TheFirstPoint, TheLastPoint, TheNumberOfPoints):
    TheLine = aNew3Dline(TheNumberOfPoints)
    TheVector = TheLastPoint.minus(TheFirstPoint)
    for i in range(TheNumberOfPoints):
        TheLine.point[i] = TheFirstPoint.plus(TheVector.times((i-0.5)/(TheNumberOfPoints - 2)))
    return TheLine

def TheDotted3Dline(TheLine, TheNumberOfDots, TheRelativeDotLength,
                             TheNumberOfPointsBeforeTheFirstDot, TheScale):
    TheDottedLine = new3Dlines(TheNumberOfDots)
    TheDotNumberOfPoints = math.floor(TheLine.numberOfPoints/TheScale/(TheNumberOfDots+1)*TheRelativeDotLength/(1+TheRelativeDotLength))
    TheVoidNumberOfPoints = math.floor(TheLine.numberOfPoints/TheScale/(TheNumberOfDots+1)/(1+TheRelativeDotLength))
    ThePeriod = TheDotNumberOfPoints + TheVoidNumberOfPoints
    for i in range(TheDottedLine.numberOfLines):
        TheDottedLine.line[i] = aNew3Dline(TheDotNumberOfPoints)
        for j in range(TheDotNumberOfPoints):
            pointNumber = TheNumberOfPointsBeforeTheFirstDot + (i*ThePeriod + j)*TheScale
            TheDottedLine.line[i].point[j] = TheLine.point[pointNumber]
    return TheDottedLine

class aColored3Dline(a3Dline):
    color = (0,0,0)

def aNewColored3Dline(TheNumberOfPoints, color):
    TheLine = aColored3Dline()
    TheLine.numberOfPoints = TheNumberOfPoints
    TheLine.point = [a3Dvector() for i in range(TheNumberOfPoints)]
    TheLine.color = color
    return TheLine

class lines3D:
    numberOfLines = 0
    line = []

    def plus(self, v):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].plus(v)
        return TheLines

    def minus(self, v):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].minus(v)
        return TheLines

    def times(self, a):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].times(a)
        return TheLines

    def XaxisRotated(self, theta):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].XaxisRotated(theta)
        return TheLines

    def YaxisRotated(self, theta):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].YaxisRotated(theta)
        return TheLines

    def ZaxisRotated(self, theta):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].ZaxisRotated(theta)
        return TheLines

    def Xdilated(self, theta):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].Xdilated(theta)
        return TheLines

    def Ydilated(self, theta):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].Ydilated(theta)
        return TheLines

    def Zdilated(self, theta):
        TheLines = new3Dlines(self.numberOfLines)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].Zdilated(theta)
        return TheLines

def new3Dlines(TheNumberOfLines):
    TheLines = lines3D()
    TheLines.numberOfLines = TheNumberOfLines
    TheLines.line = [a3Dline() for i in range(TheNumberOfLines)]
    return TheLines

class monocolored3Dlines(lines3D):
    color = (0, 0 ,0)

    def plus(self, v):
        TheLines = newMonocolored3Dlines(self.numberOfLines, self.color)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].plus(v)
        return TheLines

    def minus(self, v):
        TheLines = newMonocolored3Dlines(self.numberOfLines, self.color)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].minus(v)
        return TheLines

    def times(self, a):
        TheLines = newMonocolored3Dlines(self.numberOfLines, self.color)
        for i in range(self.numberOfLines):
            TheLines.line[i] = self.line[i].times(a)
        return TheLines

def newMonocolored3Dlines(TheNumberOfLines, TheColor):
    TheLines = monocolored3Dlines()
    TheLines.color = TheColor
    TheLines.numberOfLines = TheNumberOfLines
    TheLines.line = [a3Dline() for i in range(TheNumberOfLines)]
    return TheLines

class multicolored3Dlines:
    numberOfLines = 0
    line = []

def new3Dmulticolored3Dlines(TheNumberOfLines):
    TheLines = multicolored3Dlines()
    TheLines.numberOfLines = TheNumberOfLines
    TheLines.line = [aColored3Dline() for i in range(TheNumberOfLines)]
    return TheLines

Vector3D.py[edit]

# 2023, October 18. Version 231018

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 scalar(self, v):
        return self.x*v.x + self.y*v.y + self.z*v.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

    def Xdilated(self, a):
        v = a3Dvector()
        v.x = self.x*a
        v.y = self.y
        v.z = self.z
        return v

    def Ydilated(self, a):
        v = a3Dvector()
        v.x = self.x
        v.y = self.y*a
        v.z = self.z
        return v

    def Zdilated(self, a):
        v = a3Dvector()
        v.x = self.x
        v.y = self.y
        v.z = self.z*a
        return v

    # Consider the projection of the point v = (x,y,z) through
    # the center (0, 0, 0) on the plane XZ, y = -1 followed
    # by a half turn rotation around the y axis: a point centered
    # projected on XZ is the image of the point 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()
        p.x = 0.0
        p.y = 0.0
        proj = 0
        if self.y > 1e-10:
            proj = 1
            p.x = self.x/self.y
            p.y = self.z/self.y
        return [p, proj]

def The3Dnorm(v):
    return math.sqrt(v.x*v.x + v.y*v.y + v.z*v.z)

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

DepthMatrix.py[edit]

# 2023, October 11. Version 231011.2

# A depth matrix is the memory of the perception of the depth of the object seen
# (its distance relative to the observer). A depth matrix associates to each pixel
# the ordered list, according to proximity, of line corners which are projected on
# it. A line corner is memorized in the depth matrix with the distance of its
# projected point to the pixel and its distance to the optical center (its depth).

# A multicolor depth matrix associates with each pixel the ordered list, according to
# proximity, of color emitters (line corners, atoms, facets...) which are projected on
# # it. A color emitter is memorized in the depth matrix with the distance (which
# determine its brightness) of its projected point to the pixel, its color and its
# distance to the optical center (its depth).

import json
from PIL import ImageColor

class aDepthMatrix:
	width = 0
	height = 0
	numberOfPoints = []
	dist = [] # the distance to the pixel center.
	depth = [] # the distance to the optical center.

	def isWrittenOnFile(self, TheFileName):
		with open(TheFileName, 'w') as f:
			txtMat = []
			txtMat += [self.width]
			txtMat += [self.height]
			for i in range(self.width):
				for j in range(self.height):
					TheNumberOfPoints = self.numberOfPoints[i][j]
					txtMat += [TheNumberOfPoints]
					for k in range(TheNumberOfPoints):
						txtMat += [self.dist[i][j][k]]
						txtMat += [self.depth[i][j][k]]
			json.dump(txtMat, f)

def TheDepthMatrixOnFile(TheFileName):
	with open(TheFileName, 'r') as f:
		txtMat = json.load(f)
		TheWidth = int(txtMat[0])
		TheHeight = int(txtMat[1])
		TheMatrix = aNewDepthMatrix(TheWidth, TheHeight)
		count = 2
		for i in range(TheMatrix.width):
			for j in range(TheMatrix.height):
				TheNumberOfPoints = int(txtMat[count])
				count += 1
				TheMatrix.numberOfPoints[i][j] = TheNumberOfPoints
				TheMatrix.dist[i][j] = [0.0 for i in range(TheNumberOfPoints)]
				TheMatrix.depth[i][j] = [0.0 for i in range(TheNumberOfPoints)]
				for k in range(TheNumberOfPoints):
					TheMatrix.dist[i][j][k] = float(txtMat[count])
					count += 1
					TheMatrix.depth[i][j][k] = float(txtMat[count])
					count += 1
	return TheMatrix

def aNewDepthMatrix(TheWidth, TheHeight):
	TheMatrix = aDepthMatrix()
	TheMatrix.width = TheWidth
	TheMatrix.height = TheHeight
	TheMatrix.numberOfPoints = [[] for i in range(TheWidth)]
	TheMatrix.dist = [[] for i in range(TheWidth)]
	TheMatrix.depth = [[] for i in range(TheWidth)]
	for i in range(TheWidth):
		TheMatrix.numberOfPoints[i] = [0 for j in range(TheHeight)]
		TheMatrix.dist[i] = [[] for j in range(TheHeight)]
		TheMatrix.depth[i] = [[] for j in range(TheHeight)]
	return TheMatrix

class aMultiColorDepthMatrix:
	color = [] # The color of the perceived point.

# 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 the two lists.
def Order2lists(list1, list2, n):
	if n == 2:
		if list1[1] < list1[0]:
			storage1 = list1[0]
			list1[0] = list1[1]
			list1[1] = storage1
			storage2 = list2[0]
			list2[0] = list2[1]
			list2[1] = storage2
	if n > 2:
		n1 = n // 2
		n2 = n - n1
		l11 = [0.0 for i in range(n1)]
		l21 = [0.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.0 for i in range(n2)]
		for i in range(n2):
			l12[i] = list1[n1 + i]
			l22[i] = list2[n1 + i]
		order1 = Order2lists(l11, l21, n1)
		l11 = order1[0]
		l21 = order1[1]
		order2 = Order2lists(l12, l22, n2)
		l12 = order2[0]
		l22 = order2[1]
		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
	return [list1, list2]

# The following method orders the first list, of real numbers, in its
# natural order and the two next ones, of real numbers and of triplets of
# integers, according to the order of the first one. n is the number
# of items in the three lists.
def Order3lists(list1, list2, list3, n):
	if n == 2:
		if list1[1] < list1[0]:
			storage1 = list1[0]
			list1[0] = list1[1]
			list1[1] = storage1
			storage2 = list2[0]
			list2[0] = list2[1]
			list2[1] = storage2
			storage3 = list3[0]
			list3[0] = list3[1]
			list3[1] = storage3
	if n > 2:
		n1 = n // 2
		n2 = n - n1
		l11 = [0.0 for i in range(n1)]
		l21 = [0.0 for i in range(n1)]
		l31 = [(0, 0, 0) for i in range(n1)]
		for i in range(n1):
			l11[i] = list1[i]
			l21[i] = list2[i]
			l31[i] = list3[i]
		l12 = [0.0 for i in range(n2)]
		l22 = [0.0 for i in range(n2)]
		l32 = [(0, 0, 0) for i in range(n2)]
		for i in range(n2):
			l12[i] = list1[n1 + i]
			l22[i] = list2[n1 + i]
			l32[i] = list3[n1 + i]
		order1 = Order3lists(l11, l21, l31, n1)
		l11 = order1[0]
		l21 = order1[1]
		l31 = order1[2]
		order2 = Order3lists(l12, l22, l32, n2)
		l12 = order2[0]
		l22 = order2[1]
		l32 = order2[2]
		i1 = 0
		i2 = 0
		i3 = 0
		while i3 < n:
			if i1 == n1:
				list1[i3] = l12[i2]
				list2[i3] = l22[i2]
				list3[i3] = l32[i2]
				i2 += 1
			else:
				if i2 == n2:
					list1[i3] = l11[i1]
					list2[i3] = l21[i1]
					list3[i3] = l31[i1]
					i1 += 1
				else:
					if l11[i1] < l12[i2]:
						list1[i3] = l11[i1]
						list2[i3] = l21[i1]
						list3[i3] = l31[i1]
						i1 += 1
					else:
						list1[i3] = l12[i2]
						list2[i3] = l22[i2]
						list3[i3] = l32[i2]
						i2 += 1
			i3 += 1
	return [list1, list2, list3]


Atom.py[edit]

# 2023, October 18. Version 231018

from Vector3D import a3Dvector
from Line3D import monocolored3Dlines

class atoms:
    numberOfAtoms = 0
    radius = 1.0
    color = (0, 0, 0)
    position = []

    def minus(self, v):
        TheAtoms = newAtoms(self.numberOfAtoms, self.radius, self.color)
        for i in range(self.numberOfAtoms):
            TheAtoms.position[i] = self.position[i].minus(v)
        return TheAtoms

def newAtoms(TheNumberOfAtoms, TheRadius, TheColor):
    TheAtoms = atoms()
    TheAtoms.numberOfAtoms = TheNumberOfAtoms
    TheAtoms.radius = TheRadius
    TheAtoms.color = TheColor
    TheAtoms.position = [ a3Dvector() for i in range(TheNumberOfAtoms)]
    return TheAtoms


class atomsOfDifferentKinds:
    numberOfKinds = 0
    kind = []

    def minus(self, v):
        TheAtoms = newAtomsOfDifferentKinds(self.numberOfKinds)
        for i in range(self.numberOfKinds):
            TheAtoms.kind[i] = self.kind[i].minus(v)
        return TheAtoms

def newAtomsOfDifferentKinds(TheNumberOfKinds):
    TheAtoms = atomsOfDifferentKinds()
    TheAtoms.numberOfKinds = TheNumberOfKinds
    TheAtoms.kind = [atoms() for i in range(TheNumberOfKinds)]
    return TheAtoms

class bondedAtoms:
    atoms = atoms()
    bonds = monocolored3Dlines()

    def minus(self, v):
        TheBondedAtoms = bondedAtoms()
        TheBondedAtoms.atoms = self.atoms.minus(v)
        TheBondedAtoms.bonds = self.bonds.minus(v)
        return TheBondedAtoms

class bondedAtomsOfDifferentKinds:
    atoms = atomsOfDifferentKinds()
    bonds = monocolored3Dlines()

    def minus(self, v):
        TheBondedAtoms = bondedAtomsOfDifferentKinds()
        TheBondedAtoms.atoms = self.atoms.minus(v)
        TheBondedAtoms.bonds = self.bonds.minus(v)
        return TheBondedAtoms

Multicolor3Dpainter.py[edit]

# 2023, October 19. Version 231019

# The multicolor depth matrix:
# The multicolor depth matrix associates with each pixel the ordered list, according to
# proximity, of color emitters (line corners, atoms, facets...) which are projected on
# it. A color emitter is memorized in the depth matrix with the distance (which
# determine its brightness) of its projected point to the pixel, its color and its
# distance to the optical center (its depth).

# On a black screen, the superposition of red and green is yellow. On a white screen,
# the same superposition is black. On a black screen, colors are colored light and
# they are added: more colors, more light. On a white screen, colors are filters of
# light: more colors, more filters, less light.

from Vector2D import The2Dnorm
from Vector3D import The3Dnorm
from Monocolor3Dpainter import aMonocolor3Dpainter
from DepthMatrix import aMultiColorDepthMatrix, aNewDepthMatrix, Order3lists
import math
from time import process_time

class aMulticolorCanvas:
    redShade = []
    greenShade = []
    blueShade = []

# A multicolor 3D painter has all the properties and methods of a monocolor 3D painter:
class aMulticolor3Dpainter(aMonocolor3Dpainter):
    multicolorCanvas = aMulticolorCanvas()
    multicolorDepthMatrix = aMultiColorDepthMatrix()

    def takesAnewCanvasMulticolor(self, TheWidth, TheHeight, TheLengthUnit, style, color):
        self.takesAnewCanvas(TheWidth, TheHeight, TheLengthUnit, style, color)
        self.multicolorCanvas.redShade = [[0.0 for j in range(self.canvas.height)] for i in range(self.canvas.width)]
        self.multicolorCanvas.greenShade = [[0.0 for j in range(self.canvas.height)] for i in range(self.canvas.width)]
        self.multicolorCanvas.blueShade = [[0.0 for j in range(self.canvas.height)] for i in range(self.canvas.width)]
        # The depth matrix is refreshed:
        self.depthMatrix = aNewDepthMatrix(self.canvas.width, self.canvas.height)
        self.multicolorDepthMatrix.color = [ [] for i in range(self.canvas.width)]
        for i in range(self.canvas.width):
            self.multicolorDepthMatrix.color[i] = [ [] for i in range(self.canvas.height)]

    # 'draws3Dlines' is a competence of the painter:
    def drawsMulticolored3DlinesMulticolor(self, TheLines):
        print("Time:", process_time(), "The drawer fills the depth matrix.")
        self.perceivesMulticoloredLinesMulticolor(TheLines)
        self.paintsAviewMulticolor()

    def drawsMonocolored3DlinesMulticolor(self, TheLines):
        print("Time:", process_time(), "The drawer fills the depth matrix.")
        self.perceivesMonocoloredLinesMulticolor(TheLines)
        self.paintsAviewMulticolor()

    def drawsA3DlineMulticolor(self, TheLine):
        print("Time:", process_time(), "The drawer fills the depth matrix.")
        self.perceivesAlineMulticolor(TheLine, self.paintbrush.color)
        self.paintsAviewMulticolor()

    # The perception of lines fills the depth matrix with lines:
    def perceivesMulticoloredLinesMulticolor(self, TheLines):
        for i in range(TheLines.numberOfLines):
            if i % 50000 == 0:
                print("Time:", process_time(), "The drawer perceives line", i)
            self.perceivesAlineMulticolor(TheLines.line[i], TheLines.line[i].color)

    def perceivesMonocoloredLinesMulticolor(self, TheLines):
        for i in range(TheLines.numberOfLines):
            if i % 50000 == 0:
                print("Time:", process_time(), "The drawer perceives line", i)
            self.perceivesAlineMulticolor(TheLines.line[i], TheLines.color)

    # 'perceives a line' fills the depth matrix with a line:
    def perceivesAlineMulticolor(self, TheLine, color):
        for i in range(TheLine.numberOfPoints - 2):
            if i > 0 and i % 50000 == 0:
                print("Time:", process_time(), "The drawer perceives point", i)
            linePoint1 = TheLine.point[i]
            linePoint2 = TheLine.point[i + 1]
            linePoint3 = TheLine.point[i + 2]
            v1 = linePoint2.minus(linePoint1).times(0.5)
            v2 = linePoint3.minus(linePoint2).times(0.5)
            cornerPoint1 = linePoint1.plus(v1)
            cornerPoint2 = linePoint2
            cornerPoint3 = linePoint2.plus(v2)
            self.perceivesAcornerMulticolor(cornerPoint1, cornerPoint2, cornerPoint3, color)

    # 'perceives a corner' fills the depth matrix with a corner:
    def perceivesAcornerMulticolor(self, cornerPoint1, cornerPoint2, cornerPoint3, color):
        im1 = self.pointImageOf(cornerPoint1)
        im2 = self.pointImageOf(cornerPoint2)
        im3 = self.pointImageOf(cornerPoint3)
        if im1[1] == 1 and im2[1] == 1 and im3[1] == 1:
            cornerPoint1image = im1[0]
            cornerPoint2image = im2[0]
            cornerPoint3image = im3[0]
            fourIntegers = self.cornerDrawingWindow(cornerPoint1image, cornerPoint2image, cornerPoint3image,
                                        self.paintbrush.lineHalfWidth*(self.brush3D.lineFringeRelativeWidth +1))
            left = fourIntegers[0]
            top = fourIntegers[1]
            right = fourIntegers[2]
            bottom = fourIntegers[3]
            for i in range(left, right):
                for j in range(top, bottom):
                    vF = cornerPoint2image.minus(cornerPoint1image)
                    vS = cornerPoint3image.minus(cornerPoint2image)
                    vF2 = vF.scalar(vF)
                    vS2 = vS.scalar(vS)
                    ThePixelCenter = self.canvas.convertedToRealCoordinates(i, j)  # The center of a pixel in real coordinates
                    # Projection of the center of the pixel on the first segment of the corner:
                    uF = ThePixelCenter.minus(cornerPoint1image)
                    pF = uF.scalar(vF)
                    if vF2 > 1E-20:
                        projF = cornerPoint1image.plus(vF.times(pF / vF2))
                    else:
                        projF = cornerPoint2image
                    distF = The2Dnorm(ThePixelCenter.minus(projF))
                    # Projection of the center of the pixel on the second segment of the corner:
                    uS = ThePixelCenter.minus(cornerPoint2image)
                    pS = uS.scalar(vS)
                    if vS2 > 1E-20:
                        projS = cornerPoint2image.plus(vS.times(pS / vS2))
                    else:
                        projS = cornerPoint2image
                    distS = The2Dnorm(ThePixelCenter.minus(projS))
                    dist = 1E20
                    cornerPoint4 = cornerPoint2
                    if pF >= 0 and pF < vF2:
                        dist = distF
                        if vF2 > 1E-20:
                            cornerPoint4 = cornerPoint1.plus(cornerPoint2.minus(cornerPoint1).times(pF / vF2))
                    if pS >= 0 and pS < vS2:
                        dist = distS
                        if vS2 > 1E-20:
                            cornerPoint4 = cornerPoint2.plus(cornerPoint3.minus(cornerPoint2).times(pS / vS2))
                        if pF >= 0 and pF < vF2 and distF < distS:
                            dist = distF
                            if vF2 > 1E-20:
                                cornerPoint4 = cornerPoint1.plus(cornerPoint2.minus(cornerPoint1).times(pF / vF2))
                    if pF >= vF2 and pS < 0:
                        dist = The2Dnorm(ThePixelCenter.minus(cornerPoint2image))
                    if dist < self.paintbrush.lineHalfWidth*(self.brush3D.lineFringeRelativeWidth +1):
                        self.depthMatrix.numberOfPoints[i][j] += 1
                        self.depthMatrix.dist[i][j] += [dist/self.paintbrush.lineHalfWidth]
                        self.depthMatrix.depth[i][j] += [The3Dnorm(cornerPoint4.minus(self.opticalCenter))]
                        self.multicolorDepthMatrix.color[i][j] += [color]

    # The atoms are of the same kind:
    def paintsAtomsOfTheSameKindMulticolor(self, TheAtoms):
        print("Time:", process_time(), "The drawer fills the depth matrix with", TheAtoms.numberOfAtoms,
              "atoms.")
        self.perceivesAtomsMulticolor(TheAtoms)
        self.paintsAviewMulticolor()

    def paintsBondedAtomsOfTheSameKindMulicolor(self, TheBondedAtoms):
        print("Time:", process_time(), "The drawer fills the depth matrix with", TheBondedAtoms.atoms.numberOfAtoms,
              "atoms and", TheBondedAtoms.bonds.numberOfLines, "bonds.")
        self.perceivesAtomsMulticolor(TheBondedAtoms.atoms)
        self.perceivesMonocoloredLinesMulticolor(TheBondedAtoms.bonds)
        self.paintsAviewMulticolor()

    # The atoms are of different kinds:
    def paintsAtomsMulticolor(self, TheAtoms):
        for i in range(TheAtoms.numberOfKinds):
            self.perceivesAtomsMulticolor(TheAtoms.kind[i])
        self.paintsAviewMulticolor()

    def paintsBondedAtomsMulticolor(self, TheBondedAtoms):
        print("Time:", process_time(), "The drawer fills the depth matrix with atoms.")
        for i in range(TheBondedAtoms.atoms.numberOfKinds):
            self.perceivesAtomsMulticolor(TheBondedAtoms.atoms.kind[i])
        self.perceivesMonocoloredLinesMulticolor(TheBondedAtoms.bonds)
        self.paintsAviewMulticolor()

    def perceivesAtomsMulticolor(self, TheAtoms):
        for i in range(TheAtoms.numberOfAtoms):
            if i % 500000 == 0:
                print("Time:", process_time(), "The drawer perceives atom", i)
            self.perceivesAnAtomMulticolor(TheAtoms.position[i], TheAtoms.radius, TheAtoms.color)

    def perceivesAnAtomMulticolor(self, TheAtomCenter, TheAtomRadius, color):
        TheAtomCenterDepth = The3Dnorm(TheAtomCenter.minus(self.opticalCenter))
        TheVisibleRadius = TheAtomRadius/TheAtomCenterDepth*self.zoom
        if self.brush3D.atomStyle == "bright":
            TheVisibleRadius = max(TheVisibleRadius, self.paintbrush.lineHalfWidth*(self.brush3D.lineFringeRelativeWidth+1))
        im = self.pointImageOf(TheAtomCenter)
        if im[1] == 1:
            TheAtomImage = im[0]
            fourIntegers = self.diskDrawingWindow(TheAtomImage, TheVisibleRadius)
            left = fourIntegers[0]
            top = fourIntegers[1]
            right = fourIntegers[2]
            bottom = fourIntegers[3]
            for i in range(left, right):
                for j in range(top, bottom):
                    ThePixelCenter = self.canvas.convertedToRealCoordinates(i,j)
                    dist = The2Dnorm(TheAtomImage.minus(ThePixelCenter))
                    if dist < TheVisibleRadius:
                        self.depthMatrix.numberOfPoints[i][j] += 1
                        # The minus sign in front of 'dist' and the -1 distinguishes atomes from corners:
                        self.depthMatrix.dist[i][j] += [-dist/TheVisibleRadius - 1]
                        if self.brush3D.atomStyle == "bright":
                            depth = TheAtomCenterDepth - TheAtomRadius
                        else:
                            depth = TheAtomCenterDepth - TheAtomRadius*math.sqrt(pow(TheVisibleRadius, 2)
                                                                    - dist*dist)/TheVisibleRadius
                        self.depthMatrix.depth[i][j] += [depth]
                        self.multicolorDepthMatrix.color[i][j] += [color]

    def paintsAviewMulticolor(self):
        print("Time:", process_time(), "The drawer orders the depth matrix.")
        self.ordersTheDepthMatrixMulticolor()
        print("Time:", process_time(), "The drawer paints.")
        Max = 0.0
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                if self.depthMatrix.numberOfPoints[i][j] > 0:
                    self.multicolorCanvas.redShade[i][j] = self.shadeMulticolor(i, j)[0]
                    self.multicolorCanvas.greenShade[i][j] = self.shadeMulticolor(i, j)[1]
                    self.multicolorCanvas.blueShade[i][j] = self.shadeMulticolor(i, j)[2]
                    if Max < max(self.multicolorCanvas.redShade[i][j], self.multicolorCanvas.greenShade[i][j],
                                 self.multicolorCanvas.blueShade[i][j]):
                        Max = max(self.multicolorCanvas.redShade[i][j], self.multicolorCanvas.greenShade[i][j],
                                  self.multicolorCanvas.blueShade[i][j])
        print("Time:", process_time(), "shade max =", Max, "intensity =", self.brush3D.intensity)
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                self.multicolorCanvas.redShade[i][j] *= self.brush3D.intensity
                self.multicolorCanvas.greenShade[i][j] *= self.brush3D.intensity
                self.multicolorCanvas.blueShade[i][j] *= self.brush3D.intensity
        self.drawsOnCanvasMulticolor()

    def ordersTheDepthMatrixMulticolor(self):
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                numberOfPoints = self.depthMatrix.numberOfPoints[i][j]
                if numberOfPoints > 1:
                    orderDepth = Order3lists(self.depthMatrix.depth[i][j], self.depthMatrix.dist[i][j],
                                             self.multicolorDepthMatrix.color[i][j], numberOfPoints)
                    self.depthMatrix.depth[i][j] = orderDepth[0]
                    self.depthMatrix.dist[i][j] = orderDepth[1]
                    self.multicolorDepthMatrix.color[i][j] = orderDepth[2]

    def shadeMulticolor(self, i, j):
        numberOfPoints = self.depthMatrix.numberOfPoints[i][j]
        totalRed = 0.0
        totalGreen = 0.0
        totalBlue = 0.0
        if numberOfPoints > 0:
            pointTakeColor = [1.0 for i in range(numberOfPoints)]
            totalTakeColor = 1.0
            totalTakeColorRank = -1
            depthRank1 = 0
            while depthRank1 < numberOfPoints:
                dist1 = self.depthMatrix.dist[i][j][depthRank1]
                color1 = self.multicolorDepthMatrix.color[i][j][depthRank1]
                if dist1 >= 0:
                    giveColor1 = self.giveColor(dist1)
                    pointTakeColor[depthRank1] = self.lineTakeColor(dist1)
                else:
                    giveColor1 = self.atomGiveColor(abs(dist1) - 1)
                    pointTakeColor[depthRank1] = self.atomTakeColor(abs(dist1) - 1)
                if depthRank1 > 0:
                    depth1 = self.depthMatrix.depth[i][j][depthRank1]
                    depthRank2 = depthRank1 - 1
                    depth2 = self.depthMatrix.depth[i][j][depthRank2]
                    while depthRank2 >= 0 and depth2 > depth1 - self.brush3D.lineDepth:
                        grade = (depth1 - depth2) / self.brush3D.lineDepth
                        giveColor1 *= 1 - grade * (1 - pointTakeColor[depthRank2])
                        depthRank2 -= 1
                        if depthRank2 >= 0:
                            depth2 = self.depthMatrix.depth[i][j][depthRank2]
                    totalTakeColorRankStorage = depthRank2
                    while depthRank2 >= 0 and depthRank2 > totalTakeColorRank:
                        totalTakeColor *= pointTakeColor[depthRank2]
                        depthRank2 -= 1
                    totalTakeColorRank = totalTakeColorRankStorage
                giveColor1 *= totalTakeColor
                totalRed += giveColor1*(color1[0] - self.canvas.color[0])
                totalGreen += giveColor1*(color1[1] - self.canvas.color[1])
                totalBlue += giveColor1*(color1[2] - self.canvas.color[2])
                if totalTakeColor*(numberOfPoints - depthRank1)*self.brush3D.intensity < 1/255:
                    depthRank1 = numberOfPoints
                depthRank1 += 1
        return [totalRed, totalGreen, totalBlue]

    def drawsOnCanvasMulticolor(self):
        print("draws on canvas", self.canvas.style)
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                redShade = math.floor(self.canvas.color[0] + self.multicolorCanvas.redShade[i][j])
                if redShade < 255:
                    if redShade < 0:
                        redShade = 0
                else:
                    redShade = 255
                if self.canvas.style == "black on white":
                    self.canvas.pixel[i, j] = 255 - redShade
                else:
                    if self.canvas.style == "white on black":
                        self.canvas.pixel[i, j] = redShade
                    else:
                        if self.canvas.style == "color":
                            greenShade = math.floor(self.canvas.color[1] + self.multicolorCanvas.greenShade[i][j])
                            if greenShade < 255:
                                if greenShade < 0:
                                    greenShade = 0
                            else:
                                greenShade = 255
                            blueShade = math.floor(self.canvas.color[2] + self.multicolorCanvas.blueShade[i][j])
                            if blueShade < 255:
                                if blueShade < 0:
                                    blueShade = 0
                            else:
                                blueShade = 255
                            self.canvas.pixel[i, j] = (redShade, greenShade, blueShade)


Palette.py[edit]

# 2023, October 18. Version 231018

# A painter is endowed with a few palettes. A palette converts a real number
# (between 0 and 1, or between -1 and 1) into a color defined by its three
# components (red, green, blue) each an integer between 0 and 255.
# Black is defined by (0, 0, 0), white by (255, 255, 255).

import math

class aPalette:
    # Fundamental colors:
    white = (255, 255, 255)
    black = (0, 0, 0)
    red = (255, 0, 0)
    yellow = (255, 255, 0)
    green = (0, 255, 0)
    cyan = (0, 255, 255)
    blue = (0, 0, 255)
    purple = (255, 0, 255)

    gold = (255, 215, 0)
    quarterGrey = (192, 192, 192)
    middleGrey = (128, 128, 128)
    deepRed1 = (192, 0, 0)
    deepRed2 = (128, 0, 0)
    deepGreen1 = (0, 192, 0)
    deepGreen2 = (0, 128, 0 )
    deepBlue1 = (0, 0, 192)
    deepBlue2 = (0, 0, 128 )

    # Rainbow flag colors:
    redF = (229, 0, 0)
    orangeF = (255, 141, 0)
    yellowF = (255, 238, 0)
    greenF = (2, 129, 33)
    blueF = (0, 76, 255)
    purpleF = (119, 0, 136)

    # Periodic rainbow flag colors:
    redPF = (229, 0, 0)
    orangePF = (255, 160, 0)
    yellowPF = (255, 238, 0)
    greenPF = (2, 129, 33)
    bluePF = (0, 76, 255)
    purplePF = (220, 0, 220)

    # For the palettes: gradation, heatColors, depthColors, psychedelicColors,
    # morePsychedelicColors and WolpertFlagColors the shade is always supposed
    # to be a real number between 0 and 1.

    def gradation(self, shade, color1, color2):
        # Gradation from color1 (shade=0) to color2 (shade=1).
        if shade < 1:
            if shade<0:
                shade = 0
        else:
            shade = 1
        redShade = math.floor(color1[0] + shade * (color2[0] - color1[0]))
        greenShade = math.floor(color1[1] + shade * (color2[1] - color1[1]))
        blueShade = math.floor(color1[2] + shade * (color2[2] - color1[2]))
        return (redShade, greenShade, blueShade)

    def heatColors(self, shade):
        # black, red, yellow, white
        if shade > 2/3:
            return self.gradation(3*(shade-2/3), self.yellow, self.white)
        else:
            if shade > 1/3:
                return self.gradation(3*(shade - 1/3), self.red, self.yellow)
            else:
                if shade >= 0:
                    return self.gradation(3*shade, self.black, self.red)
                else:
                    return self.black

    def depthColors(self, shade):
        # black, blue, cyan, white
        if shade > 2/3:
            return self.gradation(3*(shade-2/3), self.cyan, self.white)
        else:
            if shade > 1/3:
                return self.gradation(3*(shade - 1/3), self.blue, self.cyan)
            else:
                if shade >= 0:
                    return self.gradation(3*shade, self.black, self.blue)
                else:
                    return self.black

    def psychedelicColors(self, shade):
        # black, blue, purple, red, yellow, white
        if shade > 4/5:
            return self.gradation(5*(shade- 4/5), self.yellow, self.white)
        else:
            if shade > 3/5:
                return self.gradation(5*(shade- 3/5), self.red, self.yellow)
            else:
                if shade > 2/5:
                    return self.gradation(5*(shade- 2/5), self.purple, self.red)
                else:
                    if shade > 1/5:
                        return self.gradation(5*(shade- 1/5), self.blue, self.purple)
                    else:
                        if shade >= 0:
                            return self.gradation(5*shade, self.black, self.blue)
                        else:
                            return self.black

    def rainbowFlagColors(self, shade):
        # red, orange, yellow, green, blue, purple
        if shade > 0.8:
            return self.gradation((shade - 0.8)*5, self.blueF, self.purpleF)
        else:
            if shade > 0.6:
                return self.gradation((shade - 0.6)*5, self.greenF, self.blueF)
            else:
                if shade > 0.4:
                   return self.gradation((shade - 0.4)*5, self.yellowF, self.greenF)
                else:
                    if shade > 0.2:
                       return self.gradation((shade - 0.2)*5, self.orangeF, self.yellowF)
                    else:
                        return self.gradation(shade*5, self.redF, self.orangeF)

    def morePsychedelicColors(self, shade):
        # black, green, blue, purple, red, yellow, white
        if shade > 6/7:
            return self.gradation(7*(shade- 6/7), self.yellow, self.white)
        else:
            if shade > 5/7:
                return self.gradation(7*(shade- 5/7), self.red, self.yellow)
            else:
                if shade > 4/7:
                    return self.gradation(7*(shade- 4/7), self.purple, self.red)
                else:
                    if shade > 3 / 7:
                        return self.gradation(7*(shade- 3/7), self.blue, self.purple)
                    else:
                        if shade > 2 / 7:
                            return self.gradation(7*(shade- 2/7), self.deepBlue1, self.blue)
                        else:
                            if shade > 1/7:
                                return self.gradation(7*(shade- 1/7), self.deepGreen1, self.deepBlue1)
                            else:
                                if shade >= 0:
                                    return self.gradation(7*shade, self.black, self.deepGreen1)
                                else:
                                    return self.black

    def WolpertFlagColors(self, shade):
        # black, red, grey, blue
        shadeInt = min(255, math.floor(shade*255))
        if shade > 2/3:
            return (0, 0, shadeInt) # blue
        else:
            if shade > 1/3:
                return (shadeInt, shadeInt, shadeInt) # grey
            else:
                if shade >= 0:
                    return (shadeInt, 0, 0) # red
                else:
                    return self.black

    def periodicThreeColors(self, shade, color1, color2, color3):
        # The shade is supposed to be a real number between 0 and 1.
        if shade < 1/3:
            return self.gradation(3*shade, color1, color2)
        else:
            if shade < 2/3:
                return self.gradation(3*(shade - 1/3), color2, color3)
            else:
                return self.gradation(3*(shade - 2/3), color3, color1)

    def periodicYellowGreenCyan(self, shade):
        # The color number and the brightness are supposed to be
        # real numbers between 0 and 1.
        return self.periodicThreeColors(shade, self.yellow, self.green, self.cyan)

    def periodicRainbowFlag(self, shade):
        # red, orange, yellow, green, blue, purple
        if shade > 5/6:
            return self.gradation((shade - 5/6)*6, self.bluePF, self.purplePF)
        else:
            if shade > 2/3:
                return self.gradation((shade - 2/3)*6, self.greenPF, self.bluePF)
            else:
                if shade > 0.5:
                   return self.gradation((shade - 0.5)*6, self.yellowPF, self.greenPF)
                else:
                    if shade > 1/3:
                       return self.gradation((shade - 1/3)*6, self.orangePF, self.yellowPF)
                    else:
                        if shade > 1/6:
                            return self.gradation((shade - 1/6)*6, self.redPF, self.orangePF)
                        else:
                            return self.gradation(shade*6, self.purplePF, self.redPF)

    def periodicBrightRainbow(self, shade):
        # red, orange, yellow, green, cyan, blue, purple
        if shade > 5/6:
            return self.gradation((shade - 5/6)*6, self.blue, self.purple)
        else:
            if shade > 4/6:
                return self.gradation((shade - 4/6)*6, self.cyan, self.blue)
            else:
                if shade > 3/6:
                   return self.gradation((shade - 3/6)*6, self.green, self.cyan)
                else:
                    if shade > 2/6:
                       return self.gradation((shade - 2/6)*6, self.yellow, self.green)
                    else:
                        if shade > 1/6:
                            return self.gradation((shade - 1/6)*6, self.red, self.yellow)
                        else:
                            return self.gradation(shade*6, self.purple, self.red)

    def silverBlackGold(self, shade):
        # The shade is supposed to be  a real number between
        # -1 (silver) and 1 (gold). 0 is black.
        return self.doubleGradation(shade, self.white, self.black, self.gold)

    def doubleGradation(self, shade, color1, color2, color3):
        # The shade is supposed to be  a real number between
        # -1 (color1) and 1 (color3). 0 is color2.
        if shade >= 0:
            shadeCut = min(1.0, shade)
            return self.gradation(shadeCut, color2, color3)
        if shade < 0:
            shadeCut = min(1.0, -shade)
            return self.gradation(shadeCut, color2, color1)

    def depthHeat(self, shade): # blue, cyan, white, yellow, red.
        # The shade is supposed to be  a real number between
        # -1 (blue) and 1 (red). 0 is white.
        if shade >= 0:
            if shade < 1/2:
                return self.gradation(2*shade, self.white, self.yellow)
            else:
                if shade < 1:
                    return self.gradation(2*(shade- 1/2), self.yellow, self.red)
                else:
                    return self.red
        if shade < 0:
            if -shade < 1/2:
                return self.gradation(-2*shade, self.white, self.cyan)
            else:
                if -shade < 1:
                    return self.gradation(2*(-shade-1/2), self.cyan, self.blue)
                else:
                    return self.blue

    def coneCoordinatesColor(self,greyShade, greyDeviation, periodicColor):
        grey = self.gradation(greyShade, (0, 0, 0), (255, 255, 255))
        brightColor = self.periodicBrightRainbow(periodicColor)
        if greyShade < 0.5 and greyShade >= 0:
            paleColor = self.gradation(2 * greyShade, (0, 0, 0), brightColor)
            return self.gradation(greyDeviation, grey, paleColor)
        else:
            if greyShade <= 1:
                paleColor = self.gradation(2 - 2 * greyShade, (255, 255, 255), brightColor)
                return self.gradation(greyDeviation, grey, paleColor)
            else:
                return (255, 255, 255)

The complements[edit]

SaveTheWorld.py[edit]

# 2023, October 19. Version 231019.3

# The drawer calculates, perceives and draws a world, here a line. Then it writes
# on file this world and its perception: the line and the depth matrix.
# The line and its perception can be retrieved with the file 'TheWorldRebirth.py'.

from Monocolor3Dpainter import aMonocolor3Dpainter
from Line3D import a3Dline
from Vector3D import The3Dvector
from time import process_time

print ("Time:", process_time(), "Mathematical painter")

# The drawer:
TheDrawer = aMonocolor3Dpainter()

# The canvas must be defined:
TheCanvasStyle =  "black on white" # "color", "black on white" or "white on black"
TheWidth =  500 # number of pixels
TheHeight = 200 # number of pixels
TheLengthUnit = 100 # number of pixels
TheCanvasColor = (255, 255, 255)# (red, green, blue)
TheDrawer.takesAnewCanvas(TheWidth, TheHeight, TheLengthUnit, TheCanvasStyle, TheCanvasColor)

# The line of vision must be defined:
Phi = 90 # in degrees. The line of vision is parallel to the y axis.
Theta = 90 # in degrees. The line of vision is horizontal.
TheObjectCenter = The3Dvector(0, 0, 0)
TheDistance = 4 # The distance between the object center and the optical center.
TheZoom = 2
Psi = 0 # in degrees.The angle of rotation of the image around the line of vision.

TheDrawer.choosesApointOfView(TheObjectCenter, TheDistance, Phi, Theta, Psi, TheZoom)

# The paintbrush must be defined :
TheDrawer.paintbrush.color = (0, 0, 0) # Black
TheDrawer.paintbrush.lineHalfWidth = 0.04 # (The line half width * the length unit) is the
# number of pixels in the half width of the line.
TheDrawer.brush3D.lineFringeRelativeWidth = 1
TheDrawer.brush3D.lineDepth = TheDrawer.paintbrush.lineHalfWidth*(1+TheDrawer.brush3D.lineFringeRelativeWidth)
TheDrawer.brush3D.lineOpacity = 1
TheDrawer.brush3D.intensity = 1


# The line must be defined:
TheLine = a3Dline()
# Here a definition of the line.

# The line is saved on file:
TheLine.isWrittenOnFile("lineName.txt")

print("Time:", process_time(), "The drawer draws the image.")
TheDrawer.drawsA3Dline(TheLine)

# The depth matrix is saved on file:
TheDrawer.savesAperception("depthMatrixName.txt")

TheDrawer.givesApainting("imageName.png")

print("Time:", process_time(), "Good bye")

TheWorldRebirth.py[edit]

# 2023, October 19. Version 231019.2

# The drawer paints the perception of the world written on the file
# 'depthMatrixName.txt' in the folder 'folderName', and gives the
# painting in the same folder with the name 'image.png'.
# A perception of the world can be written on file with the method
# 'saves a perception'. It can be read on file with the method 
# 'retrieves a memory'. Read the file 'SaveTheWorld.py' for an example.

from Monocolor3Dpainter import aMonocolor3Dpainter
from time import process_time

print ("Time:", process_time(), "Mathematical painter")

# The drawer:
TheDrawer = aMonocolor3Dpainter()

print("Time:", process_time(), "The drawer reads the depth matrix.")
folderName = "folderName"
depthMatrixName = folderName + "/depthMatrixName.txt"
TheDrawer.retrievesAmemory(depthMatrixName)

# The canvas:
TheCanvasStyle= "black on white" # "color", "black on white" or "white on black"
TheWidth = TheDrawer.depthMatrix.width
TheHeight =  TheDrawer.depthMatrix.height
TheCanvasColor = (255, 255, 255) # (red, green, blue)
TheDrawer.takesAnewCanvas(TheWidth, TheHeight, 1, TheCanvasStyle, TheCanvasColor)

# The paintbrush:
TheDrawer.paintbrush.color = (0, 0, 0) # Black
TheDrawer.brush3D.lineDepth = 0.05
TheDrawer.brush3D.lineOpacity = 1
TheDrawer.brush3D.intensity = 1

print("Time:", process_time(), "The drawer paints the depth matrix.")
TheDrawer.paintsTheDepthMatrix()
TheDrawer.givesApainting(folderName + "/image.png")

print("Time:", process_time(), "Good bye")

ODE.py[edit]

# 2023, October 25. Version 231025

# Ordinary Differential Equations

from Line3D import aNew3Dline, new3Dlines
import math
from time import process_time

# Runge-Kutta method:
# A system of ordinary differential equations is defined by a flow.
# If the space of states has three dimensions, the system of
# differential equations is:
# dX/dt = (dx/dt, dy/dt, dz/dt) = flow(x, y, z) = flow(X)
# The Leibniz method calculates an approximation of X(t+dt) with:
# X(t+dt) = X(t) + (dX/dt)dt = f(X)dt where f 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(TheFlow, TheInitialPosition, 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.
# The delay is the time between the initial position and the first point of the line.
def TheODEsolution(TheFlow, TheInitialPosition, The_dt, TheDelay, TheNumberOfPoints,
                   TheNumberOf_dt_betweenPoints):
    TheSolution = aNew3Dline(TheNumberOfPoints)
    TheCurrentPosition = TheInitialPosition
    # Delay before first point:
    TheDelayNumberOf_dt = math.floor(TheDelay/The_dt)
    for i in range(TheDelayNumberOf_dt):
        TheCurrentPosition = TheNextPosition(TheFlow, TheCurrentPosition, The_dt)
    # The line:
    TheSolution.point[0] = TheCurrentPosition
    for i in range(TheNumberOfPoints-1):
        if i % 50000 == 0:
            print("Time:", process_time(), "The drawer calculates point", i)
        for j in range(TheNumberOf_dt_betweenPoints):
            TheCurrentPosition = TheNextPosition(TheFlow, TheCurrentPosition, The_dt)
        TheSolution.point[i+1] = TheCurrentPosition
    return TheSolution

def TheDottedODEsolution(TheFlow, TheInitialPosition, TheDelay, The_dt, TheNumberOfDots,
                            TheDotNumberOfPoints,  TheNumberOfPointsBetweenDots,
                            TheNumberOf_dt_betweenPoints):
    TheSolution = new3Dlines(TheNumberOfDots)
    TheCurrentPosition = TheInitialPosition
    # Delay before first dot:
    TheDelayNumberOf_dt = math.floor(TheDelay/The_dt)
    for i in range(TheDelayNumberOf_dt):
        TheCurrentPosition = TheNextPosition(TheFlow, TheCurrentPosition, The_dt)
    # The line:
    for i in range(TheNumberOfDots):
        if i%100 == 0:
            print ("Time:", process_time(), "dot", i)
        TheSolution.line[i]= aNew3Dline(TheDotNumberOfPoints)
        for j in range(TheDotNumberOfPoints):
            TheSolution.line[i].point[j] = TheCurrentPosition
            for k in range(TheNumberOf_dt_betweenPoints):
                TheCurrentPosition = TheNextPosition(TheFlow, TheCurrentPosition, The_dt)
        for j in range(TheNumberOfPointsBetweenDots*TheNumberOf_dt_betweenPoints):
            TheCurrentPosition = TheNextPosition(TheFlow, TheCurrentPosition, The_dt)
    return TheSolution

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)

Lattice.py[edit]

# 2023, October 19. Version 231019

from Vector3D import a3Dvector, The3Dvector
from Atom import bondedAtoms, bondedAtomsOfDifferentKinds, newAtoms, newAtomsOfDifferentKinds
from Line3D import monocolored3Dlines, TheStraight3Dline
import math
from time import process_time

class a3Dcell(bondedAtomsOfDifferentKinds):
    numberOfAtomKinds = 1
    # X, Y and Z are the three vectors of the basis of the cell:
    X = a3Dvector()
    Y = a3Dvector()
    Z = a3Dvector()
    atomPresence = []
    numberOfBonds = 0
    bondColor = (0, 0, 0)
    bondFirstPoint = []
    bondSecndPoint = []
    bondPresence = []
    # presence = [left, right, near, far, down up]. left = 0 means that the bond or the atom
    # is not present when the cell is at a left extremity of the lattice, left = 1 means
    # that it is present. The same for right, near, far, down and up extremities.
    specialCase = "nothingToReport" # or "diamond"

    def times(self, scalar):
        hX = self.X.times(scalar)
        hY = self.Y.times(scalar)
        hZ = self.Z.times(scalar)
        TheHomotheticCell = aNew3Dcell(hX, hY, hZ, self.numberOfAtomKinds, self.numberOfBonds, self.bondColor)
        TheHomotheticCell.specialCase = self.specialCase
        for i in range(self.numberOfAtomKinds):
            TheHomotheticCell.atoms.kind[i] = newAtoms(self.atoms.kind[i].numberOfAtoms,
                                    self.atoms.kind[i].radius*scalar, self.atoms.kind[i].color)
            TheHomotheticCell.atomPresence[i] = [ [0, 0, 0, 0, 0, 0] for i in range(
                                                    self.atoms.kind[i].numberOfAtoms) ]
            for j in range(self.atoms.kind[i].numberOfAtoms):
                TheHomotheticCell.atoms.kind[i].position[j] = self.atoms.kind[i].position[j].times(scalar)
                TheHomotheticCell.atomPresence[i][j] = self.atomPresence[i][j]
        TheHomotheticCell.bondColor = self.bondColor
        for i in range(self.numberOfBonds):
            TheHomotheticCell.bondFirstPoint[i] = self.bondFirstPoint[i].times(scalar)
            TheHomotheticCell.bondSecndPoint[i] = self.bondSecndPoint[i].times(scalar)
            TheHomotheticCell.bondPresence[i] = self.bondPresence[i]
        return TheHomotheticCell

    def atomsLattice(self, width, depth, height):
        atomsOfTheLattice = newAtomsOfDifferentKinds(self.numberOfAtomKinds)
        TheLatticeCenter = self.X.times(width).plus(self.Y.times(depth)).plus(self.Z.times(height)).times(0.5)
        for kindNumber in range(self.numberOfAtomKinds):
            atomsOfTheLattice.kind[kindNumber].radius = self.atoms.kind[kindNumber].radius
            atomsOfTheLattice.kind[kindNumber].color = self.atoms.kind[kindNumber].color
            atomsOfTheLattice.kind[kindNumber].numberOfAtoms = 0
            for i in range(width + 1):
                if i % 20 == 0:
                    print("Time:", process_time(), "The drawer calculates the plane", i, "of the lattice.")
                for j in range(depth + 1):
                    for k in range(height + 1):
                        v = self.X.times(i).plus(self.Y.times(j)).plus(self.Z.times(k)).minus(TheLatticeCenter)
                        for l in range(self.atoms.kind[kindNumber].numberOfAtoms):
                            if ((i > 0 or self.atomPresence[kindNumber][l][0] == 1) and (
                                    i < width or self.atomPresence[kindNumber][l][1] == 1)
                                    and (j > 0 or self.atomPresence[kindNumber][l][2] == 1) and (
                                            j < depth or self.atomPresence[kindNumber][l][3] == 1)
                                    and (k > 0 or self.atomPresence[kindNumber][l][4] == 1) and (
                                            k < height or self.atomPresence[kindNumber][l][5] == 1)):
                                atomsOfTheLattice.kind[kindNumber].numberOfAtoms += 1
                                atomsOfTheLattice.kind[kindNumber].position = (atomsOfTheLattice.kind[kindNumber].position
                                                                    + [self.atoms.kind[kindNumber].position[l].plus(v)])
        return atomsOfTheLattice

    def bondedAtomsLattice(self, width, depth, height):
        TheLattice = bondedAtoms()
        TheLattice.atoms = newAtomsOfDifferentKinds(self.numberOfAtomKinds)
        TheLattice.bonds = monocolored3Dlines()
        TheLattice.bonds.color = self.bondColor
        TheLatticeCenter = self.X.times(width).plus(self.Y.times(depth)).plus(self.Z.times(height)).times(0.5)
        for kindNumber in range(self.numberOfAtomKinds):
            TheLattice.atoms.kind[kindNumber].radius = self.atoms.kind[kindNumber].radius
            TheLattice.atoms.kind[kindNumber].color = self.atoms.kind[kindNumber].color
            TheLattice.atoms.kind[kindNumber].numberOfAtoms = 0
        for i in range(width + 1):
            for j in range(depth + 1):
                for k in range(height + 1):
                    v = self.X.times(i).plus(self.Y.times(j)).plus(self.Z.times(k)).minus(TheLatticeCenter)
                    # Atoms:
                    for kindNumber in range(self.numberOfAtomKinds):
                        for l in range(self.atoms.kind[kindNumber].numberOfAtoms):
                            if ((i > 0 or self.atomPresence[kindNumber][l][0] == 1) and (
                                    i < width or self.atomPresence[kindNumber][l][1] == 1)
                                    and (j > 0 or self.atomPresence[kindNumber][l][2] == 1) and (
                                            j < depth or self.atomPresence[kindNumber][l][3] == 1)
                                    and (k > 0 or self.atomPresence[kindNumber][l][4] == 1) and (
                                            k < height or self.atomPresence[kindNumber][l][5] == 1)):
                                TheLattice.atoms.kind[kindNumber].numberOfAtoms += 1
                                TheLattice.atoms.kind[kindNumber].position = TheLattice.atoms.kind[kindNumber].position + [self.atoms.kind[kindNumber].position[l].plus(v)]
                    # Bonds:
                    for l in range(self.numberOfBonds):
                        if ((i > 0 or self.bondPresence[l][0] == 1) and (
                                i < width or self.bondPresence[l][1] == 1)
                                and (j > 0 or self.bondPresence[l][2] == 1) and (
                                        j < depth or self.bondPresence[l][3] == 1)
                                and (k > 0 or self.bondPresence[l][4] == 1) and (
                                        k < height or self.bondPresence[l][5] == 1)):
                            yes = True
                            if self.specialCase == "diamond":
                                yes = diamondBondPresence(i, j, k, l, width, depth, height)
                            if yes:
                                TheLattice.bonds.numberOfLines += 1
                                TheLattice.bonds.line += [TheStraight3Dline(self.bondFirstPoint[l].plus(
                                    v), self.bondSecndPoint[l].plus(v), 3)]
        return TheLattice

def aNew3Dcell(X, Y, Z, TheNumberOfAtomKinds, TheNumberOfBonds, TheBondColor):
    TheCell = a3Dcell()
    TheCell.X = X
    TheCell.Y = Y
    TheCell.Z = Z
    TheCell.numberOfAtomKinds = TheNumberOfAtomKinds
    TheCell.atoms = newAtomsOfDifferentKinds(TheNumberOfAtomKinds)
    TheCell.atomPresence = [[] for i in range(TheNumberOfAtomKinds)]
    TheCell.numberOfBonds = TheNumberOfBonds
    TheCell.bondColor = TheBondColor
    TheCell.bondFirstPoint = [a3Dvector() for i in range(TheNumberOfBonds)]
    TheCell.bondSecndPoint = [a3Dvector() for i in range(TheNumberOfBonds)]
    TheCell.bondPresence = [[0, 0, 0, 0, 0, 0] for i in range(TheNumberOfBonds)]
    return TheCell

def TheCubicCell(side, TheBondColor, TheAtomRadius, TheAtomColor):
    X = The3Dvector(1, 0, 0)
    Y = The3Dvector(0, 1, 0)
    Z = The3Dvector(0, 0, 1)
    TheCell = aNew3Dcell(X, Y, Z, 1, 3, TheBondColor)
    TheCell.atoms.kind[0] = newAtoms(1, TheAtomRadius, TheAtomColor)
    TheCell.atoms.kind[0].position[0] = The3Dvector(0, 0, 0)
    TheCell.atomPresence[0] = [ [0, 0, 0, 0, 0, 0] for i in range(1)]
    TheCell.atomPresence[0][0] =  [1, 1, 1, 1, 1, 1]
    TheCell.bondColor = TheBondColor
    TheCell.bondFirstPoint[0] = The3Dvector(0, 0, 0)
    TheCell.bondSecndPoint[0] = The3Dvector(1, 0, 0)
    TheCell.bondPresence[0] = [1, 0, 1, 1, 1, 1]
    TheCell.bondFirstPoint[1] = The3Dvector(0, 0, 0)
    TheCell.bondSecndPoint[1] = The3Dvector(0, 1, 0)
    TheCell.bondPresence[1] = [1, 1, 1, 0, 1, 1]
    TheCell.bondFirstPoint[2] = The3Dvector(0, 0, 0)
    TheCell.bondSecndPoint[2] = The3Dvector(0, 0, 1)
    TheCell.bondPresence[2] = [1, 1, 1, 1, 1, 0]
    return TheCell.times(side)

def TheBodyCenteredCubicCell(side, TheBondColor, TheAtomRadius, TheAtomColor):
    X = The3Dvector(2, 0, 0)
    Y = The3Dvector(0, 2, 0)
    Z = The3Dvector(0, 0, 2)
    TheCell = aNew3Dcell(X, Y, Z, 1, 8, TheBondColor)
    TheCell.atoms.kind[0] = newAtoms(2, TheAtomRadius, TheAtomColor)
    TheCell.atomPresence[0] = [ [0, 0, 0, 0, 0, 0] for i in range(2)]
    TheCell.atoms.kind[0].position[0] = The3Dvector(1, 1, 0)
    TheCell.atomPresence[0][0] =  [1, 0, 1, 0, 1, 1]
    TheCell.atoms.kind[0].position[1] = The3Dvector(0, 0, 1)
    TheCell.atomPresence[0][1] =  [1, 1, 1, 1, 1, 0]
    TheCell.bondColor = TheBondColor
    TheCell.bondFirstPoint[0] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[0] = The3Dvector(0, 0, 1)
    TheCell.bondPresence[0] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[1] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[1] = The3Dvector(0, 2, 1)
    TheCell.bondPresence[1] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[2] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[2] = The3Dvector(2, 2, 1)
    TheCell.bondPresence[2] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[3] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[3] = The3Dvector(2, 0, 1)
    TheCell.bondPresence[3] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[4] = The3Dvector(1, 1, 2)
    TheCell.bondSecndPoint[4] = The3Dvector(0, 0, 1)
    TheCell.bondPresence[4] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[5] = The3Dvector(1, 1, 2)
    TheCell.bondSecndPoint[5] = The3Dvector(0, 2, 1)
    TheCell.bondPresence[5] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[6] = The3Dvector(1, 1, 2)
    TheCell.bondSecndPoint[6] = The3Dvector(2, 2, 1)
    TheCell.bondPresence[6] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[7] = The3Dvector(1, 1, 2)
    TheCell.bondSecndPoint[7] = The3Dvector(2, 0, 1)
    TheCell.bondPresence[7] = [1, 0, 1, 0, 1, 0]
    return TheCell.times(side/2)

def TheFaceCenteredCubicCell(side, TheBondColor, TheAtomRadius, TheAtomColor):
    X = The3Dvector(2, 0, 0)
    Y = The3Dvector(0, 2, 0)
    Z = The3Dvector(0, 0, 2)
    TheCell = aNew3Dcell(X, Y, Z, 1, 12, TheBondColor)
    TheCell.atoms.kind[0] = newAtoms(4, TheAtomRadius, TheAtomColor)
    TheCell.atomPresence[0] = [ [0, 0, 0, 0, 0, 0] for i in range(4)]
    TheCell.atoms.kind[0].position[0] = The3Dvector(0, 0, 0)
    TheCell.atomPresence[0][0] =  [1, 1, 1, 1, 1, 1]
    TheCell.atoms.kind[0].position[1] = The3Dvector(1, 1, 0)
    TheCell.atomPresence[0][1] =  [1, 0, 1, 0, 1, 1]
    TheCell.atoms.kind[0].position[2] = The3Dvector(0, 1, 1)
    TheCell.atomPresence[0][2] =  [1, 1, 1, 0, 1, 0]
    TheCell.atoms.kind[0].position[3] = The3Dvector(1, 0, 1)
    TheCell.atomPresence[0][3] = [1, 0, 1, 1, 1, 0]
    TheCell.bondColor = TheBondColor
    TheCell.bondFirstPoint[0] = The3Dvector(0, 0, 0)
    TheCell.bondSecndPoint[0] = The3Dvector(1, 1, 0)
    TheCell.bondPresence[0] = [1, 0, 1, 0, 1, 1]
    TheCell.bondFirstPoint[1] = The3Dvector(0, 0, 0)
    TheCell.bondSecndPoint[1] = The3Dvector(0, 1, 1)
    TheCell.bondPresence[1] = [1, 1, 1, 0, 1, 0]
    TheCell.bondFirstPoint[2] = The3Dvector(0, 0, 0)
    TheCell.bondSecndPoint[2] = The3Dvector(1, 0, 1)
    TheCell.bondPresence[2] = [1, 0, 1, 1, 1, 0]
    TheCell.bondFirstPoint[3] = The3Dvector(0, 2, 0)
    TheCell.bondSecndPoint[3] = The3Dvector(1, 1, 0)
    TheCell.bondPresence[3] = [1, 0, 1, 0, 1, 1]
    TheCell.bondFirstPoint[4] = The3Dvector(0, 2, 0)
    TheCell.bondSecndPoint[4] = The3Dvector(0, 1, 1)
    TheCell.bondPresence[4] = [1, 1, 1, 0, 1, 0]
    TheCell.bondFirstPoint[5] = The3Dvector(2, 0, 0)
    TheCell.bondSecndPoint[5] = The3Dvector(1, 0, 1)
    TheCell.bondPresence[5] = [1, 0, 1, 1, 1, 0]
    TheCell.bondFirstPoint[6] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[6] = The3Dvector(2, 2, 0)
    TheCell.bondPresence[6] = [1, 0, 1, 0, 1, 1]
    TheCell.bondFirstPoint[7] = The3Dvector(0, 2, 2)
    TheCell.bondSecndPoint[7] = The3Dvector(0, 1, 1)
    TheCell.bondPresence[7] = [1, 1, 1, 0, 1, 0]
    TheCell.bondFirstPoint[8] = The3Dvector(2, 0, 2)
    TheCell.bondSecndPoint[8] = The3Dvector(1, 0, 1)
    TheCell.bondPresence[8] = [1, 0, 1, 1, 1, 0]
    TheCell.bondFirstPoint[9] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[9] = The3Dvector(2, 0, 0)
    TheCell.bondPresence[9] = [1, 0, 1, 0, 1, 1]
    TheCell.bondFirstPoint[10] = The3Dvector(0, 1, 1)
    TheCell.bondSecndPoint[10] = The3Dvector(0, 0, 2)
    TheCell.bondPresence[10] = [1, 1, 1, 0, 1, 0]
    TheCell.bondFirstPoint[11] = The3Dvector(1, 0, 1)
    TheCell.bondSecndPoint[11] = The3Dvector(0, 0, 2)
    TheCell.bondPresence[11] = [1, 0, 1, 1, 1, 0]
    return TheCell.times(side/2)

def TheDiamondCell(side, TheBondColor, TheAtomRadius, TheAtomColor):
    X = The3Dvector(4, 0, 0)
    Y = The3Dvector(0, 4, 0)
    Z = The3Dvector(0, 0, 4)
    TheCell = aNew3Dcell(X, Y, Z, 1, 16, TheBondColor)
    TheCell.specialCase = "diamond"
    TheCell.atoms.kind[0] = newAtoms(8, TheAtomRadius, TheAtomColor)
    TheCell.atomPresence[0] = [ [0, 0, 0, 0, 0, 0] for i in range(8)]
    TheCell.atoms.kind[0].position[0] = The3Dvector(0, 0, 0)
    TheCell.atomPresence[0][0] = [0, 0, 0, 0, 0, 0]
    TheCell.atoms.kind[0].position[1] = The3Dvector(2, 2, 0)
    TheCell.atomPresence[0][1] = [1, 0, 1, 0, 1, 1]
    TheCell.atoms.kind[0].position[2] = The3Dvector(0, 2, 2)
    TheCell.atomPresence[0][2] = [1, 1, 1, 0, 1, 0]
    TheCell.atoms.kind[0].position[3] = The3Dvector(2, 0, 2)
    TheCell.atomPresence[0][3] = [1, 0, 1, 1, 1, 0]
    TheCell.atoms.kind[0].position[4] = The3Dvector(3, 3, 3)
    TheCell.atomPresence[0][4] = [1, 0, 1, 0, 1, 0]
    TheCell.atoms.kind[0].position[5] = The3Dvector(3, 1, 1)
    TheCell.atomPresence[0][5] = [1, 0, 1, 0, 1, 0]
    TheCell.atoms.kind[0].position[6] = The3Dvector(1, 3, 1)
    TheCell.atomPresence[0][6] = [1, 0, 1, 0, 1, 0]
    TheCell.atoms.kind[0].position[7] = The3Dvector(1, 1, 3)
    TheCell.atomPresence[0][7] = [1, 0, 1, 0, 1, 0]
    TheCell.bondColor = TheBondColor
    TheCell.bondFirstPoint[0] = The3Dvector(1, 3, 1)
    TheCell.bondSecndPoint[0] = The3Dvector(0, 2, 2)
    TheCell.bondPresence[0] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[1] = The3Dvector(1, 1, 3)
    TheCell.bondSecndPoint[1] = The3Dvector(0, 2, 2)
    TheCell.bondPresence[1] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[2] = The3Dvector(3, 1, 1)
    TheCell.bondSecndPoint[2] = The3Dvector(2, 0, 2)
    TheCell.bondPresence[2] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[3] = The3Dvector(1, 1, 3)
    TheCell.bondSecndPoint[3] = The3Dvector(2, 0, 2)
    TheCell.bondPresence[3] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[4] = The3Dvector(3, 1, 1)
    TheCell.bondSecndPoint[4] = The3Dvector(2, 2, 0)
    TheCell.bondPresence[4] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[5] = The3Dvector(1, 3, 1)
    TheCell.bondSecndPoint[5] = The3Dvector(2, 2, 0)
    TheCell.bondPresence[5] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[6] = The3Dvector(3, 3, 3)
    TheCell.bondSecndPoint[6] = The3Dvector(4, 4, 4)
    TheCell.bondPresence[6] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[7] = The3Dvector(3, 3, 3)
    TheCell.bondSecndPoint[7] = The3Dvector(4, 2, 2)
    TheCell.bondPresence[7] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[8] = The3Dvector(3, 3, 3)
    TheCell.bondSecndPoint[8] = The3Dvector(2, 4, 2)
    TheCell.bondPresence[8] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[9] = The3Dvector(3, 3, 3)
    TheCell.bondSecndPoint[9] = The3Dvector(2, 2, 4)
    TheCell.bondPresence[9] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[10] = The3Dvector(3, 1, 1)
    TheCell.bondSecndPoint[10] = The3Dvector(4, 2, 2)
    TheCell.bondPresence[10] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[11] = The3Dvector(3, 1, 1)
    TheCell.bondSecndPoint[11] = The3Dvector(4, 0, 0)
    TheCell.bondPresence[11] = [1, 0, 0, 0, 0, 0]
    TheCell.bondFirstPoint[12] = The3Dvector(1, 3, 1)
    TheCell.bondSecndPoint[12] = The3Dvector(2, 4, 2)
    TheCell.bondPresence[12] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[13] = The3Dvector(1, 3, 1)
    TheCell.bondSecndPoint[13] = The3Dvector(0, 4, 0)
    TheCell.bondPresence[13] = [0, 0, 1, 0, 0, 0]
    TheCell.bondFirstPoint[14] = The3Dvector(1, 1, 3)
    TheCell.bondSecndPoint[14] = The3Dvector(2, 2, 4)
    TheCell.bondPresence[14] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[15] = The3Dvector(1, 1, 3)
    TheCell.bondSecndPoint[15] = The3Dvector(0, 0, 4)
    TheCell.bondPresence[15] = [0, 0, 0, 0, 1, 0]
    return TheCell.times(side/4)

def diamondBondPresence(i, j, k, l, width, depth, height):
    yes = True
    if (i == width - 1 or j == depth - 1 or k == height - 1) and l == 6:
        yes = False
    if i == width - 1 and l == 11:
        yes = False
    if k == height - 1 and l == 15:
        yes = False
    if j == depth - 1 and l == 13:
        yes = False
    return yes

def TheTwoAtomKindsBodyCenteredCubicCell(side, TheBondColor, TheAtomRadius1, TheAtomColor1,
                                         TheAtomRadius2, TheAtomColor2):
    X = The3Dvector(2, 0, 0)
    Y = The3Dvector(0, 2, 0)
    Z = The3Dvector(0, 0, 2)
    TheCell = aNew3Dcell(X, Y, Z, 2, 8, TheBondColor)
    TheCell.atoms.kind[0] = newAtoms(1, TheAtomRadius1, TheAtomColor1)
    TheCell.atomPresence[0] = [ [0, 0, 0, 0, 0, 0] for i in range(1)]
    TheCell.atoms.kind[0].position[0] = The3Dvector(1, 1, 0)
    TheCell.atomPresence[0][0] =  [1, 0, 1, 0, 1, 1]
    TheCell.atoms.kind[1] = newAtoms(1, TheAtomRadius2, TheAtomColor2)
    TheCell.atomPresence[1] = [ [0, 0, 0, 0, 0, 0] for i in range(1)]
    TheCell.atoms.kind[1].position[0] = The3Dvector(0, 0, 1)
    TheCell.atomPresence[1][0] =  [1, 1, 1, 1, 1, 0]
    TheCell.bondColor = TheBondColor
    TheCell.bondFirstPoint[0] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[0] = The3Dvector(0, 0, 1)
    TheCell.bondPresence[0] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[1] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[1] = The3Dvector(0, 2, 1)
    TheCell.bondPresence[1] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[2] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[2] = The3Dvector(2, 2, 1)
    TheCell.bondPresence[2] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[3] = The3Dvector(1, 1, 0)
    TheCell.bondSecndPoint[3] = The3Dvector(2, 0, 1)
    TheCell.bondPresence[3] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[4] = The3Dvector(1, 1, 2)
    TheCell.bondSecndPoint[4] = The3Dvector(0, 0, 1)
    TheCell.bondPresence[4] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[5] = The3Dvector(1, 1, 2)
    TheCell.bondSecndPoint[5] = The3Dvector(0, 2, 1)
    TheCell.bondPresence[5] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[6] = The3Dvector(1, 1, 2)
    TheCell.bondSecndPoint[6] = The3Dvector(2, 2, 1)
    TheCell.bondPresence[6] = [1, 0, 1, 0, 1, 0]
    TheCell.bondFirstPoint[7] = The3Dvector(1, 1, 2)
    TheCell.bondSecndPoint[7] = The3Dvector(2, 0, 1)
    TheCell.bondPresence[7] = [1, 0, 1, 0, 1, 0]
    return TheCell.times(side/2)

Scalar1DfieldPainter.py[edit]

## 2023, October 20. Version 231020

# A scalar 1D field drawer draws and paints a scalar 1D field.
# If the field is real, the drawing is a line which represents the field value.
# The painting represents the field value with brightness and or colors.
# If the field is complex, the drawing is a line which represents the squared
# modulus of the field. The painting represents the phase of the field with
# a color, and its squared modulus with brightness.

import math, cmath
from Vector2D import The2Dvector
from Painter2D import a2Dpainter


class aScalar1DfieldPainter(a2Dpainter): # This means that a scalar 1D field drawer has all
    # the properties and methods of a 2D line drawer.
    def drawsAreal1Dfield(self, TheField, TheFieldMaximumValue):
        self.drawsA2Dline(TheField.times(self.canvas.height/self.canvas.lengthUnit/TheFieldMaximumValue).line2D())

    def paintsAreal1Dfield(self, TheField, TheFieldMaximumValue, ThePalette):
        for i in range(self.canvas.width):
            x = math.floor(i/(self.canvas.width-1)*(TheField.numberOfPoints-1))
            shade = TheField.value[x]/TheFieldMaximumValue
            for j in range(self.canvas.height):
                self.canvas.pixel[i, j] = ThePalette(shade)

    # The Field is resized
    def paintsAndDrawsAreal1Dfield(self, TheField, TheFieldMaximumValue, ThePalette):
        self.takesAnewCanvas(self.canvas.width, self.canvas.height, self.canvas.lengthUnit, "color", self.canvas.color)
        self.drawsAreal1Dfield(TheField, TheFieldMaximumValue)
        self.paintsAreal1Dfield(TheField, TheFieldMaximumValue, ThePalette)
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                grey = 1 - self.canvas.pixelShade[i][j]
                red = math.floor(self.canvas.pixel[i, j][0]*grey)
                green = math.floor(self.canvas.pixel[i, j][1]*grey)
                blue = math.floor(self.canvas.pixel[i, j][2]*grey)
                self.canvas.pixel[i, j] = (red, green, blue)

    def drawsAsquaredModulusField(self, TheField, ThePixelNumberOfPoints):
        for i in range(math.floor(TheField.width/ThePixelNumberOfPoints)-2):
            TheFirstPoint = The2Dvector( (i+0.5)*ThePixelNumberOfPoints*TheField.dx,
                                         pow(abs(TheField.value[i*ThePixelNumberOfPoints]),2) )
            TheSecondPoint = The2Dvector( (i+1.5)*ThePixelNumberOfPoints*TheField.dx,
                                          pow(abs(TheField.value[(i+1)*ThePixelNumberOfPoints]),2) )
            TheThirdPoint = The2Dvector((i+2.5)*ThePixelNumberOfPoints*TheField.dx,
                                          TheField.value[(i+2)*ThePixelNumberOfPoints] )
            self.drawsA2Dcorner(TheFirstPoint, TheSecondPoint, TheThirdPoint)
        self.drawsOnCanvas()

    def paintsAcomplex1Dfield(self, TheField, ThePixelNumberOfPoints, maxField2, ThePalette):
        for i in range(self.canvas.width):
            psi = TheField.value[i*ThePixelNumberOfPoints]
            psiPhase = cmath.phase(psi)
            colorNumber = (psiPhase + math.pi)/(2*math.pi) # a real number between 0 and 1
            psiMod = abs(psi)
            brightness = psiMod*psiMod/maxField2 # a real number between 0 and 1
            color = ThePalette(colorNumber, brightness)
            for j in range(self.canvas.height):
                self.canvas.pixelShade[i][j] = color

    def paintsAndDrawsAcomplex1Dfield(self, TheField, ThePixelNumberOfPoints, maxField2, ThePalette):
        self.takesAnewCanvas(self.canvas.width, self.canvas.height, self.canvas.lengthUnit, "black on white",
                             self.canvas.color)
        self.drawsAsquaredModulusField(TheField, ThePixelNumberOfPoints)
        pixelMem = [ [] for i in range(self.canvas.width)]
        for i in range(self.canvas.width):
            pixelMem[i] = [0 for j in range(self.canvas.height)]
            for j in range(self.canvas.height):
                pixelMem[i][j] = self.canvas.pixelShade[i][j]
        self.takesAnewCanvas(self.canvas.width, self.canvas.height, self.canvas.lengthUnit, "color", self.canvas.color)
        self.paintsAcomplex1Dfield(TheField, ThePixelNumberOfPoints, maxField2, ThePalette)
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                grey = pixelMem[i][j]
                red = math.floor(self.canvas.pixelShade[i][j][0]*grey/255)
                green = math.floor(self.canvas.pixelShade[i][j][1]*grey/255)
                blue = math.floor(self.canvas.pixelShade[i][j][2]*grey/255)
                self.canvas.pixelShade[i][j] = (red, green, blue)

Scalar1Dfield.py[edit]

# 2023, October 18. Version 231018

from Vector2D import aNew2Dline, The2Dvector
import math


class aScalar1Dfield:
    numberOfPoints = 0  # number of points
    value = []  # the value (a real or complex number) of the field at a point.
    dx = 0.001  # the distance between two points in the real image.
    xStart = 0.0
    scalarType = "real"  # or "complex"
    limitCondition = "periodic"  # or "constant"

    def nextField(self, TheFlow, The_dt):  # This means that the field follows its flow.
        TheNextField = self.plus(TheFlow(self).times(The_dt))
        if self.limitCondition == "constant":
            TheNextField.value[0] = self.value[0]
            TheNextField.value[self.numberOfPoints - 1] = self.value[self.numberOfPoints - 1]
        return TheNextField

    # lc means limit condition. lc(i) is the value of the column number i according to
    # the chosen limit condition.
    def lc(self, i):
        j = i
        if self.limitCondition == "periodic":
            if i < 0:
                j += self.numberOfPoints
            else:
                if i > self.numberOfPoints - 1:
                    j -= self.numberOfPoints
        if self.limitCondition == "constant":
            if i < 0:
                j = 0
            else:
                if i > self.numberOfPoints - 1:
                    j = self.numberOfPoints - 1
        return j

    def ddx(self):  # ddx means d/dx, the first derivative (the gradient) of the field.
        TheDerivedField = aNew1DscalarField(self.numberOfPoints, self.dx, self.scalarType,
                                            self.limitCondition)
        for i in range(self.numberOfPoints):
            TheDerivedField.value[i] = (self.value[self.lc(i + 1)] -
                                        self.value[self.lc(i - 1)]) / 2 / self.dx
        return TheDerivedField

    # Only for a real 1D field:
    def line2D(self):
        TheLine = aNew2Dline(self.numberOfPoints)
        for i in range(self.numberOfPoints):
            TheLine.point[i] = The2Dvector(self.xStart + i*self.dx, self.value[i])
        return TheLine

    # Only for a complex 1D field:
    # def The3Dline(self):

    def plus(self, TheField):  # addition of two fields
        fieldSum = aNew1DscalarField(self.numberOfPoints, self.dx, self.scalarType, self.limitCondition)
        for i in range(self.numberOfPoints):
            fieldSum.value[i] = self.value[i] + TheField.value[i]
        return fieldSum

    def times(self, TheScalar):  # multiplication of a field by a scalar
        fieldProd = aNew1DscalarField(self.numberOfPoints, self.dx, self.scalarType, self.limitCondition)
        for i in range(self.numberOfPoints):
            fieldProd.value[i] = self.value[i] * TheScalar
        return fieldProd

    def verticallyShifted(self, shift):
        TheShiftedField = aNew1DscalarField(self.numberOfPoints, self.dx, self.scalarType,
                                            self.limitCondition)
        for i in range(self.numberOfPoints):
            TheShiftedField.value[i] = self.value[i] + shift
        return TheShiftedField

    def shifting1Dflow(self, c):  # df/dt = -c df/dx.
        # It is a flow for a progressive wave, with celerity c.
        return self.ddx().times(-c)

    def diffusion1Dflow(self, sigma):  # df/dt = sigma2 d2f/dx2. It is a flow
        # for a diffusive wave.
        return self.ddx().ddx().times(sigma * sigma)

    def free1DquantumParticleFlow(self, m):  # Schrödinger equation for a free 1D particle:
        # d psi / dt = (i hbar / 2m) (d2 psi / dx2)
        # Here hbar = 1
        return self.ddx().ddx().times(1j / (2 * m))


def aNew1DscalarField(TheNumberOfPoints, The_dx, TheScalarType, TheLimitCondition):
    TheField = aScalar1Dfield()
    TheField.numberOfPoints = TheNumberOfPoints
    TheField.dx = The_dx
    TheField.scalarType = TheScalarType
    TheField.limitCondition = TheLimitCondition
    if TheScalarType == "real":
        TheField.value = [0.0 for i in range(TheNumberOfPoints)]
    if TheScalarType == "complex":
        TheField.value = [0.0j for i in range(TheNumberOfPoints)]
    return TheField


def aNewSimpleSinus1DField(TheNumberOfPoints):
    TheField = aScalar1Dfield()
    TheField.width = TheNumberOfPoints
    TheField.dx = 2 * math.pi / (TheNumberOfPoints - 1)
    TheField.scalarType = "real"
    TheField.limitCondition = "periodic"
    TheField.value = [0.0 for i in range(TheNumberOfPoints)]
    for i in range(TheNumberOfPoints):
        TheField.value[i] = math.sin(i * TheField.dx)
    return TheField


def TheSinus1Dfield(TheNumberOfPoints, TheRealWidth, k, phi):
    TheField = aScalar1Dfield()
    TheField.numberOfPoints = TheNumberOfPoints
    TheField.dx = TheRealWidth / (TheNumberOfPoints - 1)
    TheField.scalarType = "real"
    TheField.limitCondition = "periodic"
    TheField.value = [0.0 for i in range(TheNumberOfPoints)]
    for i in range(TheNumberOfPoints):
        TheField.value[i] = math.sin(k * i * TheField.dx + phi)
    return TheField


def TheGaussian1Dfield(TheNumberOfPoints, TheRealWidth, sigma, mu, max):
    TheField = aScalar1Dfield()
    TheField.numberOfPoints = TheNumberOfPoints
    TheField.dx = TheRealWidth / (TheNumberOfPoints - 1)
    TheField.scalarType = "real"
    TheField.value = [0.0 for i in range(TheNumberOfPoints)]
    for i in range(TheNumberOfPoints):
        x = i * TheField.dx - TheRealWidth / 2
        TheField.value[i] = gauss(x, sigma, mu, max)
    return TheField


def TheGaussian1DdiffusionField(TheNumberOfPoints, TheRealWidth, sigma, mu, t):
    TheField = aScalar1Dfield()
    TheField.numberOfPoints = TheNumberOfPoints
    TheField.dx = TheRealWidth / (TheNumberOfPoints - 1)
    TheField.scalarType = "real"
    TheField.value = [0.0 for i in range(TheNumberOfPoints)]
    sigma_t = sigma * pow(t, 0.5)
    max = 1 / sigma_t
    for i in range(TheNumberOfPoints):
        x = i * TheField.dx - TheRealWidth / 2
        TheField.value[i] = gauss(x, sigma_t, mu, max)
    return TheField


def gauss(x, sigma, mu, max):  # the gaussian function
    return max * math.exp((mu - x) * (x - mu) / sigma / sigma / 2)


def TheGaussianQuantum1DwavePacket(TheNumberOfPoints, TheRealWidth, m, a, k0, xStart, t):
    TheField = aScalar1Dfield()
    TheField.numberOfPoints = TheNumberOfPoints
    TheField.dx = TheRealWidth / (TheNumberOfPoints - 1)
    TheField.scalarType = "complex"
    TheField.value = [0.0j for i in range(TheNumberOfPoints)]
    for i in range(TheNumberOfPoints):
        x = i * TheField.dx - TheRealWidth / 2 - xStart
        TheField.value[i] = gaussianQuantum1DwavePacket(m, a, k0, x, t)
    return TheField


def gaussianQuantum1DwavePacket(m, a, k0, x, t):
    # m is the mass, a the packet width, k0 the momentum, x the x coordinate and t the time.
    r1 = pow(pow(a, 4) + 4 * t * t / (m * m), 0.25)
    theta = 0.5 * math.atan(2 * t / (m * a * a))
    phi = -theta - k0 * k0 * t / (2 * m)
    c1 = pow(math.e, 1j * (phi + k0 * x) - pow(x - k0 * t / m, 2) / (a * a + 2j * t / m))
    return pow(2 * a * a / math.pi, 0.25) * c1 / r1

Scalar2DfieldPainter.py[edit]

# 2023, October 20. Version 231020

# A scalar 2D field painter paints each pixel of the canvas with a color
# and a brightness which represent the value of the field at the center
# of the pixel. The palette attributes a color to a real number, usually
# chosen between 0 and 1 (or -1 and 1).
# If the field is complex, the painting represents the phase of the field with
# a color, and its squared modulus with brightness.

import math, cmath
from Painter import aPainter

class aScalar2DfieldPainter(aPainter): # This means that a real 2D field painter
    # has all the properties and methods of a painter.

    def paintsAreal2Dfield(self, TheField, ThePixelNumberOfPoints, ThePalette):
        # The pixel number of points in the number of points of the field on the
        # side of a pixel.
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                self.canvas.pixel[i, j] = ThePalette(TheField.value[i*ThePixelNumberOfPoints][j*ThePixelNumberOfPoints])

    def paintsAcomplex2Dfield(self, TheField, ThePixelNumberOfPoints, maxField2, ThePalette):
        # The pixel number of points in the number of points of the field on the
        # side of a pixel.
        for i in range(self.canvas.width):
            for j in range(self.canvas.height):
                psi = TheField.value[i*ThePixelNumberOfPoints][j*ThePixelNumberOfPoints]
                psiPhase = cmath.phase(psi)
                colorNumber = (psiPhase + math.pi)/(2*math.pi)  # a real number between 0 and 1
                psiMod = abs(psi)
                brightness = psiMod*psiMod/maxField2  # a real number between 0 and 1
                self.canvas.pixel[i, j] = ThePalette(colorNumber, brightness)

Scalar2Dfield.py[edit]

# 2023, September 21. Version 230921.

import math

class aScalar2Dfield:
    width = 0 # number of points
    height = 0
    value = []
    dx = 0.001
    limitCondition = "periodic" # or "constant"
    scalarType = "real" # or "complex"

    def NextField(self, TheFlow, The_dt):
        TheNextField = self.plus(TheFlow(self).times(The_dt))
        if self.limitCondition == "constant":
            for i in range(self.width):
                TheNextField.value[i][0] = self.value[i][0]
                TheNextField.value[i][self.height-1] = self.value[i][self.height-1]
            for j in range(self.height):
                TheNextField.value[0][j] = self.value[0][j]
                TheNextField.value[self.width-1][j] = self.value[self.width-1][j]
        return TheNextField

    # Limit conditions:
    def ilc(self, i): # ilc(i) transforms i according to the limit conditions.
        # i is the column number of a point of the field.
        k = i
        if self.limitCondition == "periodic":
            if i<0:
                k += self.width
            else:
                if i>self.width - 1:
                    k -= self.width
        if self.limitCondition == "constant":
            if i<0:
                k = 0
            else:
                if i>self.width - 1:
                    k = self.width - 1
        return k

    def jlc(self, j): # jlc(j) transforms j according to the limit conditions.
        # j is the row number of a point of the field.
        l = j
        if self.limitCondition == "periodic":
            if j<0:
                l += self.height
            else:
                if j>self.height - 1 :
                    l -= self.height
        if self.limitCondition == "constant":
            if j<0:
                l = 0
            else:
                if j>self.height - 1:
                    l = self.height - 1
        return l

    def ddx(self): # ddx means d/dx, the partial x derivative of the field.
        TheDerivedField = aNew2DscalarField(self.width, self.height, self.dx,  self.scalarType, self.limitCondition)
        for i in range(self.width):
            for j in range(self.height):
                TheDerivedField.value[i][j] = (self.value[self.ilc(i+1)][j] - self.value[self.ilc(i-1)][j])/2/self.dx
        return TheDerivedField

    def ddy(self): # ddy means d/dy, the partial y derivative of the field.
        TheDerivedField = aNew2DscalarField(self.width, self.height, self.dx,  self.scalarType, self.limitCondition)
        for i in range(self.width):
            for j in range(self.height):
                TheDerivedField.value[i][j] = (self.value[i][self.jlc(j+1)] -
                                        self.value[i][self.jlc(j-1)])/2/self.dx
        return TheDerivedField

    def plus(self, TheField): # addition of two fields
        fieldSum = aNew2DscalarField(self.width, self.height, self.dx, self.scalarType, self.limitCondition)
        for i in range(self.width):
            for j in range(self.height):
                fieldSum.value[i][j] = self.value[i][j] + TheField.value[i][j]
        return fieldSum

    def times(self, TheScalar): # multiplication of a field by a scalar
        fieldProd = aNew2DscalarField(self.width, self.height, self.dx, self.scalarType, self.limitCondition)
        for i in range(self.width):
            for j in range(self.height):
                fieldProd.value[i][j] = self.value[i][j]*TheScalar
        return fieldProd

    def diffusion2Dflow(self, sigma):
        return self.ddx().ddx().plus(self.ddy().ddy()).times(sigma*sigma)

    def free2DquantumParticleFlow(self, m):  # Schrödinger equation for a free 2D particle:
        # d psi / dt = (i hbar / 2m) (d2 psi / dx2 + d2 psi / dy2)
        # Here hbar = 1
        return self.ddx().ddx().plus(self.ddy().ddy()).times(1j/(2*m))

def aNew2DscalarField(TheWidth, TheHeight, The_dx, TheScalarType, TheLimitCondition):
    TheField = aScalar2Dfield()
    TheField.width = TheWidth
    TheField.height = TheHeight
    TheField.dx = The_dx
    TheField.scalarType = TheScalarType
    TheField.limitCondition = TheLimitCondition
    TheField.value = [ [] for i in range(TheWidth)]
    for i in range(TheWidth):
        TheField.value[i] = [0.0 for j in range(TheHeight)]
    return TheField


def aNew2DgaussianDiffusionField(TheNumberOfPoints, TheRealWidth, sigma, t):
    TheField = aScalar2Dfield()
    TheField.width = TheNumberOfPoints
    TheField.height = TheNumberOfPoints
    TheField.dx = TheRealWidth/(TheNumberOfPoints-1)
    TheField.scalarType = "real"
    TheField.limitCondition = "periodic"
    TheField.value = [ [] for i in range(TheNumberOfPoints)]
    for i in range(TheNumberOfPoints):
        TheField.value[i] = [ 0.0 for j in range(TheNumberOfPoints)]
        for j in range(TheNumberOfPoints):
            TheField.value[i][j] = 1.2/t*math.exp(-(pow(i*TheField.dx - TheRealWidth/2, 2) +
                                                  pow(j*TheField.dx - TheRealWidth/2, 2))/(2*sigma*sigma*t))
    return TheField

def productField(field1, field2):
    if field1.scalarType == "complex" or field2.scalarType == "complex":
        TheScalarType = "complex"
    else:
        TheScalarType = "real"
    TheProduct = aNew2DscalarField(field1.width, field2.width, field1.dx,
                                   TheScalarType, field1.limitCondition)
    for i in range(TheProduct.width):
        for j in range(TheProduct.height):
            TheProduct.value[i][j] = field1.value[i]*field2.value[j]
    return TheProduct