Programmes, images

# -*- coding: utf-8 -*-

# Intersection de la feuille (du feuilletage linéaire dX/dt=X, dY/dt=mu.Y (Im(mu)>0))
# passant par un point donné à distance rayonSphere de l'origine avec la sphère S^3(rayonSphere)

# Entrées :
# - (u,v) \in S^3(rayonSphere)
# - nombre de points sur la bi-orbite
# Sortie : projection stéréographique dans R^3 de la bi-orbite de (u,v)


########################

import numpy as np
import subprocess
from povwriter import *
from multiprocessing import Manager, Process, current_process
import time
#import sys

I = complex(0,1)
J = -1./2+I*np.sqrt(3)/2
pi = np.pi
debut = time.time()

########################
########################

## Norme d'un vecteur
def normeVecteur(Z):
 return np.sqrt(abs(Z[0])**2+abs(Z[1])**2)


## Projection stéréographique de S^3(rayonSphere) dans C^2 sur R^3 comme plan équatorial
def projectionStereo(Z,rayonSphere):
 return [rayonSphere*Z[0].real/(rayonSphere-Z[1].imag),
 rayonSphere*Z[0].imag/(rayonSphere-Z[1].imag),
 rayonSphere*Z[1].real/(rayonSphere-Z[1].imag)]


## Rotation dans C^2
def rotationC2(Z):
 alpha = 1/np.sqrt(2)
 beta = 1/np.sqrt(2)
 R = np.array([[alpha,-beta.conjugate()],[beta,alpha.conjugate()]])
 return np.dot(R,Z)


## Symétrie d'ordre 3
def symetrie(Z):
 return np.array([J.conjugate()*Z[0],J*Z[1]])


########################
########################

## Calcul et dessin des deux orbites périodiques
def orbitesPeriodiques(rayonSphere) :
 discretisation = 10**3
 cercle1 = [] ; cercle2 = []
 for i in range(discretisation+1):
 cercle1.append(projectionStereo(rotationC2(np.array([rayonSphere*np.exp(2*I*pi*i/discretisation),0])),rayonSphere))
 cercle2.append(projectionStereo(rotationC2(np.array([0,rayonSphere*np.exp(2*I*pi*i/discretisation)])),rayonSphere))
 for i in range(discretisation):
 Cylinder(cercle1[i],cercle1[i+1],
 0.005,
 Texture(
 #Finish(
 # ambient = 0.0,
 # diffuse = 0.0,
 # reflection = 1,
 # specular = 1
 #),
 Pigment(color = (0,1,0)))
 ).write(image)
 Cylinder(cercle2[i],cercle2[i+1],
 0.005,
 Texture(
 #Finish(
 # ambient = 0.0,
 # diffuse = 0.0,
 # reflection = 1,
 # specular = 1
 #),
 Pigment(color = (0,0,1)))
 ).write(image)

########################
########################


def ecrirePOV(image,orbitePartielle,rayonSphere):
 for i in range(len(orbitePartielle)-1):
 Cylinder(projectionStereo(rotationC2(orbitePartielle[i]),rayonSphere),
 projectionStereo(rotationC2(orbitePartielle[i+1]),rayonSphere),
 0.005,
 Texture(
 #Finish(
 # ambient = 0.0,
 # diffuse = 0.0,
 # reflection = 0.85,
 # specular = 1
 #),
 Pigment(color = (1,1,1)))
 ).write(image)
 if (symCondInit == 1):
 Cylinder(projectionStereo(rotationC2(symetrie(orbitePartielle[i])),rayonSphere),
 projectionStereo(rotationC2(symetrie(orbitePartielle[i+1])),rayonSphere),
 0.005,
 Texture(
 #Finish(
 # ambient = 0.0,
 # diffuse = 0.0,
 # reflection = 0.85,
 # specular = 1
 #),
 Pigment(color = (0.8,0.8,0.8)))
 ).write(image)
 Cylinder(projectionStereo(rotationC2(symetrie(symetrie(orbitePartielle[i]))),rayonSphere),
 projectionStereo(rotationC2(symetrie(symetrie(orbitePartielle[i+1]))),rayonSphere),
 0.005,
 Texture(
 #Finish(
 # ambient = 0.0,
 # diffuse = 0.0,
 # reflection = 0.85,
 # specular = 1
 #),
 Pigment(color = (0.6,0.6,0.6)))
 ).write(image)


########################
########################


## Intersection du champ de vecteurs avec S^3(rayonSphere)
def champS3(X,sens):
 mu = (4*J+5)/(4*J-1)
 gamma = (abs(X[0])**2+mu.conjugate()*abs(X[1])**2)*I
 return(sens*gamma*np.array([X[0],mu*X[1]]))


## Schéma de Runge-Kutta
def rungeKutta(X,sens):
 dt = 2**(-10)
 k1 = champS3(X,sens)
 k2 = champS3(X+k1*dt/2,sens)
 k3 = champS3(X+k2*dt/2,sens)
 k4 = champS3(X+k3*dt,sens)
 return X+(k1+2*k2+2*k3+k4)*dt/6


## Calcul de la bi-orbite d'un point dans C^2 (en fait sur S^3(rayonSphere), cf. définition champS3)
def orbiteIntersectionSphere(_,nbPointsOrbite,condInit,sens):
 trajectoire = [condInit]
 for point in range(nbPointsOrbite):
 trajectoire.append(rungeKutta(trajectoire[-1],sens))
 processusCalcul[int(current_process().name.split("-")[1])-2] = trajectoire



########################
########################
########################


#symCondInit = input("Conditions initiales symétriques ? ")
#nbPointsOrbite = input("Nombre de points calculés : ")
#nbPointsOrbite = int(sys.argv[1])
symCondInit = 1
nbPointsOrbite = 6*10**3
C1 = np.array([(np.sqrt(3)+I)/(2*np.sqrt(2)),np.exp(I*pi/4)/np.sqrt(2)])
rayonSphere = normeVecteur(C1)


if __name__ == '__main__':
 image = POVFile("singularite-lineaire.pov")
 Camera(location = (0, 1, -5.5), look_at = (0, 0, 0)).write(image)
 LightSource1)10, 10, -10), color = (1, 1, 1.write(image)
 LightSource2)15, 15, -10), color = (1, 1, 1 .write(image)
 LightSource3)-15, 15, -10), color = (1, 1, 1.write(image)
 LightSource4)15, -15, -10), color = (1, 1, 1.write(image)
 LightSource5)10, 10, -10), color = (0, 0, 0.3.write(image)
 LightSource6)-10, 10, -10), color = (0, 0.3, 0.write(image)
 LightSource7)10, -10, -10), color = (0.3, 0, 0.write(image)
 
 manager = Manager()
 processusCalcul = manager.list(range(2))
 p = [Process(target=orbiteIntersectionSphere,args=(processusCalcul,nbPointsOrbite,C1,-1)),
 Process(target=orbiteIntersectionSphere,args=(processusCalcul,nbPointsOrbite,C1,1))]
 for each in p: each.start()
 for each in p: each.join()
 
 orbitesPeriodiques(rayonSphere)
 #il est bon de faire une petite rotation pour éviter qu'une des deux orbites périodiques ne rencontre le pôle de projection
 ecrirePOV(image,processusCalcul[0],rayonSphere)
 ecrirePOV(image,processusCalcul[1],rayonSphere)
 
 fin = time.time()
 print(fin-debut)

subprocess.Popen(["povray", "+A +FN +P +w1920 +h1080 /Users/Emilie/Seafile/Mathrice-Travail/Feuilletages-holomorphes-Jouanolou/Calculs-numeriques/singularite-lineaire/python-pov/singularite-lineaire.pov"])