Openscad/math.scad

/***********************************CONSTANTES***********************************/
pi = 3.1415926535;

/***********************************ARITHMÉTIQUE***********************************/
function PGCD(a,b) = a > b ? PGCD(b , a-b) : a < b ? PGCD(a, b-a) : a;

/**********************************************************************************/
/***********************************GEOMETRIE***********************************/
/**********************************************************************************/

/***********************************FONCTIONS 2D***********************************/
//CONVERSION
function rad(angle = 0) = pi * angle / 180;
function degre(angle = 0) = angle * 180 / pi;

function calcule_distance_2d(x, y) = sqrt(pow(x, 2) + pow(y, 2));

function milieu( A, B) = [ (A[0] + B[0])/2, (A[1] + B[1])/2 ];

/***************************************************************************************
calcul de l'angle formé par la droite OA avec A[x, y] et l'axe des abscisses
pour connaitre le sens de l'angle il faut vérifier si y est positif ou négatif
pour éviter les divisions par 0 il faut vérifier si y == 0
***************************************************************************************/
//function calcule_angle_2d(x, y) = y < 0 ? 360-acos(x/calcule_distance_2d(x, y)) : y == 0 && x>=0 ? 0 : acos(x/calcule_distance_2d(x, y)); 
function calcule_angle_2d(x, y) = y < 0 ? 360-acos(x/calcule_distance_2d(x, y)) : y > 0 ? acos(x/calcule_distance_2d(x, y)) : x>=0 ? 0 : 180; 

/******************************DEVELOPPANTE DE CERCLE******************************/
function developpante_x(r = 1, t = 0) = r*(cos(t)+rad(t)*sin(t));
function developpante_y(r = 1, t = 0) = r*(sin(t)-rad(t)*cos(t));

//INTERSECTION
/*On cherche à calculer le point d'intersection
d'une développante de cercle et d'un cercle.
Pour se faire il faut chercher la valeur du paramètre a dans les équations :
x(a) = r*(cos(a) + a*sin(a))
y(a) = r*(sin(a) - a*cos(a))

En utilisant pythagore on conclue que le paramètre doit satisfaire :
tel que x(a)² + y(a)² = R²
avec R le rayon du cercle que coupe la développante.

[r*(cos(a) + a*sin(a))]² + [r*(sin(a) - a*cos(a))]² = R²
r² * (cos(a) + a*sin(a))² + r² * (sin(a) - a*cos(a))² = R²

On peut factoriser par r²
r² * [ (cos(a) + a*sin(a))² + (sin(a) - a*cos(a))² ] = R²
(cos(a) + a*sin(a))² + (sin(a) - a*cos(a))² = R²/r²

On utilise les identitées remarquables :
cos(a)² + 2a*cos(a)*sin(a) + a²*sin(a)² + sin(a)² - 2a*cos(a)*sin(a) + a²*cos(a)² = R²/r²
cos(a)² + a²*sin(a)² + sin(a)² + a²*cos(a)² = R²/r²
cos(a)² + sin(a)² + a²*sin(a)² + a²*cos(a)² = R²/r²
cos(a)² + sin(a)² + a²*( sin(a)² + cos(a)² ) = R²/r²

Sachant que cos² + sin² = 1 on obtient :
1 + a² = R²/r²
a² = R²/r² - 1
a = sqrt(R²/r² - 1)
*/

function intersect_rad(r=1, R=2) = sqrt( (R*R)/(r*r) - 1);
function intersect(r=1, R=2) = degre(intersect_rad(r, R));

/***********************************FONCTIONS 3D***********************************/
function projection_xy(vecteur) = [vecteur[0], vecteur[1], 0];

function conversion_cartesien_spherique(vecteur) =
[norm(vecteur),
norm(projection_xy(vecteur)) == 0 ? 0 : vecteur[1] < 0 ? 360-acos(vecteur[0]/norm(projection_xy(vecteur))) : acos(vecteur[0]/norm(projection_xy(vecteur))),
acos(vecteur[2]/norm(vecteur))];

/***********************************MODULES***********************************/

/***************************************************************************************
Modules pour créer une spirale

spirale_coord prend en arguments :
x1 l'abscisse du point de départ
y1 l'ordonnée du point de départ
x2 l'abscisse du point d'arrivé
y2 l'ordonnée du point d'arrivé
tour désigne le nombre de tour qui devront être effectué pour passer du point de départ au point d'arrivé
largeur désigne la largeur du bras de la spirale

spirale_point prend en arguments:
A un vecteur représentant le point de départ
B un vecteur représentant le point d'arrivé
tour désigne le nombre de tour qui devront être effectué pour passer du point de départ au point d'arrivé
largeur désigne la largeur du bras de la spirale

**********************************************************************************/
module spirale_coord(x1 = 10, y1 = 0, x2 = 10, y2 = 10, tours = 2, largeur = 1)
{
    //calcule des distances par rapport à l'origine
    d1 = calcule_distance_2d(x1, y1);
    d2 = calcule_distance_2d(x2, y2);
    
    //calcule des angles entre les vecteur OA et Ox
    theta1 = calcule_angle_2d(x1, y1);
    theta2 = calcule_angle_2d(x2, y2);
    
    //calcule de l'angle à balayer pour effectuer le nombre de tour voulu
    //pour lier le point de départ et le point d'arrivé
    angle = tours * 360 + (theta2 - theta1);
    
    //la spirale est de type linéaire, cela signifie que le rayon de la courbe va croitre de manière proportionnelle à l'angle
    coeff_dir = (d2 - d1)/angle;
    
    //création des points extérieur de la spirale du point de départ vers le point d'arrivé
    A = [for(i=[0:1:angle])
        [(coeff_dir * i + d1 + largeur/2) * cos(i + theta1),
         (coeff_dir * i + d1 + largeur/2) * sin(i + theta1)]];
    
    //création des points interieur de la spirale du point d'arrivé vers le point de départ
    B = [for(i=[angle:-1:0])
        [(coeff_dir * i + d1 - largeur/2) * cos(i + theta1),
         (coeff_dir * i + d1 - largeur/2) * sin(i + theta1)]];
    
    //combinaison des 2 séries de points
    C = concat(A,B);
    
    //tracé de la courbe finale
    polygon(C);
}

module spirale_coord2(x1 = 0, y1 = 0, x2 = 0, y2 = 10, tours = 2, largeur = 1)
{
    /*pour une spirale linéaire, ou d'Archimède, la relation est la suivante :
    r = a + b*theta
    cela signifie que les coordonnées cartésiennes sont :
    x(theta) = r*cos(theta) = (a + b*theta)*cos(theta)
    y(theta) = r*sin(theta) = (a + b*theta)*sin(theta)
    */
    function delta_x(t = 0) = cos(t) - rad(t)*sin(t);
    function delta_y(t = 0) = sin(t) + rad(t)*cos(t);
    
    //calcule des distances par rapport à l'origine
    d1 = calcule_distance_2d(x1, y1);
    d2 = calcule_distance_2d(x2, y2);
    
    //calcule des angles entre les vecteur OA et Ox
    theta1 = calcule_angle_2d(x1, y1);
    theta2 = calcule_angle_2d(x2, y2);
    
    //calcule de l'angle à balayer pour effectuer le nombre de tour voulu
    //pour lier le point de départ et le point d'arrivé
    angle = tours * 360 + abs(theta2 - theta1);
        
    //la spirale est de type linéaire, cela signifie que le rayon de la courbe va croitre de manière proportionnelle à l'angle
    coeff_dir = (d2 - d1)/angle;
    
    function vecteur(t = 0) = [[0, -1],[1,  0]]*[delta_x(t = t), delta_y(t = t)]/norm([delta_x(t = t), delta_y(t = t)]);
//    function vecteur(t = 0) = [largeur*cos(t)/2, largeur*sin(t)/2];
    
    //création des points extérieur de la spirale du point de départ vers le point d'arrivé
    A = [for(i=[0:1:angle])
        (coeff_dir * i + d1) * [cos(i + theta1), sin(i + theta1)] + vecteur(t = i + theta1) ];
    
    //création des points interieur de la spirale du point d'arrivé vers le point de départ
    B = [for(i=[angle:-1:0])
        (coeff_dir * i + d1) * [cos(i + theta1), sin(i + theta1)] - vecteur(t = i + theta1) ];
    
    //combinaison des 2 séries de points
    C = concat(A,B);
    
    //tracé de la courbe finale
    polygon(C);
}

module spirale_log_coord(x1 = 1, y1 = 0, x2 = 20, y2 = 0, tours = 2, largeur = 1)
{
    //calcule des distances par rapport à l'origine
    d1 = calcule_distance_2d(x1, y1);
    d2 = calcule_distance_2d(x2, y2);
    
    //calcule des angles entre les vecteur OA et Ox
    theta1 = calcule_angle_2d(x1, y1);
    theta2 = calcule_angle_2d(x2, y2);
	
	echo("<font color='#ffa400'>math.scad, module : spirale_log_coord</font>");
	echo("<font color='#ffa400'>theta1 = ", theta1, "</font>");
	echo("<font color='#ffa400'>theta2 = ", theta2, "</font>");
    
    //calcule de l'angle à balayer pour effectuer le nombre de tour voulu
    //pour lier le point de départ et le point d'arrivé
    delta_theta = theta2 - theta1 >= 0 ? theta2 - theta1 : 360 - (theta1 - theta2);
//    angle = tours * 360 - theta1 + theta2;
    angle = tours * 360 + delta_theta;
    theta2_ = angle + theta1;
    
	echo("<font color='#ffa400'>angle = ", angle, "</font>");

    /*une spirale logarithmique est décrite par : r = a*b^theta
    Si on possède 2 points :
    r1 = a*b^theta1
    r2 = a*b^theta2
    on peut extraire le coefficient a :
    ln(r1) = ln(a*b^theta1) = ln(a) + ln(b^theta1) = ln(a) + theta1*ln(b)
    ln(r2) = ln(a*b^theta2) = ln(a) + ln(b^theta2) = ln(a) + theta2*ln(b)
    
    ln(r2) - ln(r1) = ln(a) - ln(a) + theta2*ln(b) - theta1*ln(b)
    ln(r2/r1) = ln(b) * (theta2 - theta1)
    ln(r2/r1)/(theta2-theta1) = ln(b)
    <=> b = exp( ln(r2/r1)/(theta2-theta1) )
    
    */
    b = exp( ln(d2/d1)/(theta2_-theta1) );
    a = d1/(pow(b,theta1));
    
    //création des points extérieur de la spirale du point de départ vers le point d'arrivé
    A = [for(i=[0:1:angle])
        [(a*pow(b,(i + theta1)) + largeur/2) * cos(i + theta1),
         (a*pow(b,(i + theta1)) + largeur/2) * sin(i + theta1)]];

    //création des points interieur de la spirale du point d'arrivé vers le point de départ
    B = [for(i=[angle:-1:0])
        [(a*pow(b,(i + theta1)) - largeur/2) * cos(i + theta1),
         (a*pow(b,(i + theta1)) - largeur/2) * sin(i + theta1)]];
    
    //combinaison des 2 séries de points
    C = concat(A,B);
    
    //tracé de la courbe finale
    polygon(C);
}

//spirale_coord2(x1 = 2.5, y1 = 0, x2 = 20, y2 = 0, tours = 2, largeur = 1);

module spirale_point(A = [0, 0], B = [10, 10], tours = 1, largeur = 1)
{
    spirale(A[0], A[1], B[0], B[1], tours, largeur);
}

module polygone_regulier(cotes = 5, rayon = 10)
{
    polygon([for(i=[0: 360/cotes : 360])[rayon * cos(i), rayon * sin(i)]]);
}

module portion_disque(rayon = 10, angle = 180)
{
    A = [for(i=[0 : 1 : angle])[rayon * cos(i), rayon * sin(i)]];
    B = concat([[0, 0]], A);
    polygon(B);
}

module gaussienne(sigma = 0.5, mu = 0, interval = 20)
{
    polygon([   for(i=[-interval/2:0.1:interval/2])[i,  exp(-pow( (i - mu),2))/(2*pow(sigma,2))/(sigma * sqrt(2*pi))    ] ]);
}


module trompette(epaisseur = 2, hauteur = 10, r_min = 20, r_max = 40, precision = 100, facteur = 1)
{
	x1 = ln(hauteur/(exp( facteur * (r_max - r_min)) - 1))/facteur;
	x2 = r_max - r_min + x1;
	
	function calcule_xb(xa, l) = -sqrt(pow(l, 2)/(1 + facteur * facteur * exp(-2*xa * facteur))) + xa;
	function calcule_yb(xa, l) = sqrt(pow(l, 2)/(1 + facteur * facteur * exp(2*xa * facteur))) + exp(facteur * xa);
	
	A = [for(i=[x1:(x2-x1)/precision:x2])	[i ,exp(facteur *i)]]; 
	B = [for(i=[x2:-(x2-x1)/precision:x1])	[calcule_xb(i, epaisseur), calcule_yb(i, epaisseur)]];
		
	C = concat(A,B);
	
	
	rotate_extrude()
	{
		translate([  -x2 - r_min, -exp(facteur * x1), 0])
		{polygon(C);}
	}
}