1 /***************************************************************************
2 * Copyright (C) 2021 by Abderrahman Taha *
3 * *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU General Public License as published by *
7 * the Free Software Foundation; either version 2 of the License, or *
8 * (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU General Public License *
16 * along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
19 ***************************************************************************/
20 #define PI1 (double (314159265)/double (100000000))
21 #include <vector>
22 #include <math.h>
23 unsigned int ThreadsNumber =4;
24
maxim(double p1,double p2)25 double maxim(double p1, double p2)
26 {
27 return p1 > p2 ? p1 : p2;
28 }
p_skeletal_int(const double * pp)29 double p_skeletal_int(const double* pp)
30 {
31 double cx, cy,cz;
32 cx=cos(pp[0]);
33 cy=cos(pp[1]);
34 cz=cos(pp[2]);
35 return(cx+cy+cz+0.51*(cx*cy+cy*cz+cz*cx)+1.2);
36 }
f_hex_y(const double * pp)37 double f_hex_y(const double* pp)
38 {
39 double x1,y1,x2,y2, th;
40 double p[10];
41 for(int i=0; i<4; i++)
42 p[i] = pp[i];
43 x1=fabs(fmod(fabs(p[0]), sqrt(3.0))-sqrt(3.0)/2);
44 y1=fabs(fmod(fabs(p[1]), 3)-1.5);
45 x2=sqrt(3.0)/2-x1;
46 y2=1.5-y1;
47 if ((x1*x1+y1*y1)>(x2*x2+y2*y2))
48 {
49 x1=x2;
50 y1=y2;
51 }
52 if ((x1==0.0)&&(y1==0.0))
53 p[0]=0.000001;
54 th=atan2(y1,x1);
55 if (th<PI1/6)
56 return(y1);
57 else
58 {
59 y1=-sin(PI1/3)*x1+cos(PI1/3)*y1;
60 return(fabs(y1));
61 }
62 }
fmesh(const double * pp)63 double fmesh(const double* pp) // 40
64 {
65 double th, ph, r, r2, temp;
66 double p[10];
67
68 for(int i=0; i<10; i++)
69 p[i] = pp[i];
70 th = PI1 / p[3];
71 ph = PI1/ p[4];
72 r = fmod(p[0], p[3] * 2);
73 if (r < 0)
74 r += p[3] * 2;
75 r = fabs(r - p[3]) * p[5];
76 r2 = (p[1] - cos(p[2] * ph) * p[6]) * p[7];
77 temp = -sqrt(r2 * r2 + r * r);
78 r = fmod(p[0] - p[3], p[3] * 2);
79 if (r < 0)
80 r += p[3] * 2;
81 r = fabs(r - p[3]) * p[5];
82 r2 = (p[1] + cos(p[2] * ph) * p[6]) * p[7];
83 temp = maxim(-sqrt(r2 * r2 + r * r), temp);
84 r = fmod(p[2], p[4] * 2);
85 if (r < 0)
86 r += p[4] * 2;
87 r = fabs(r - p[4]) * p[5];
88 r2 = (p[1] + cos(p[0] * th) * p[6]) * p[7];
89 temp = maxim(-sqrt(r2 * r2 + r * r), temp);
90 r = fmod(p[2] - p[4], p[4] * 2);
91 if (r < 0)
92 r += p[4] * 2;
93 r = fabs(r - p[4]) * p[5];
94 r2 = (p[1] - cos(p[0] * th) * p[6]) * p[7];
95 return (-maxim(-sqrt(r2 * r2 + r * r), temp));
96 }
fhelix1(const double * pp)97 double fhelix1(const double* pp)
98 {
99 double r, r2, r3, temp, th, ph, x2;
100 double p[10];
101 for(int i=0; i<10; i++)
102 p[i] = pp[i];
103 r = sqrt(p[0] * p[0] + p[2] * p[2]);
104 if ((p[0] == 0.0) && (p[2] == 0.0))
105 p[0] = 0.000001;
106 th = atan2(p[2], p[0]);
107 th = fmod(th * p[3] + p[1] * p[4] * p[3], 2*PI1);
108 if (th < 0)
109 th += 2*PI1;
110 p[2] = (th - PI1) / p[7] / (p[4] * p[3]);
111 p[0] = r - p[6];
112 if (p[8] == 1.0)
113 r2 = sqrt(p[0] * p[0] + p[2] * p[2]);
114 else
115 {
116 if (p[9] != 0.0)
117 {
118 th = cos(p[9] * PI1/180);
119 ph = sin(p[9] * PI1/180);
120 x2 = p[0] * th - p[2] * ph;
121 p[2] = p[0] * ph + p[2] * th;
122 p[0] = x2;
123 }
124 if (p[8] != 0.0)
125 {
126 temp = 2. / p[8];
127 r2 = pow((pow(fabs(p[0]), temp) + pow(fabs(p[2]), temp)), p[8] *.5);
128 }
129 else
130 fabs(p[0]) > fabs(p[2]) ? r2 = fabs(p[0]) : r2 = fabs(p[2]);
131 }
132 (p[6] + r) < r2 ? r3 = (p[6] + r) : r3 = r2;
133 return (-p[5] + r3);
134 }
fhelix2(const double * pp)135 double fhelix2(const double* pp) // 26
136 {
137 double th, ph, x2, z2, r2, temp;
138 double p[10];
139 for(int i=0; i<10; i++)
140 p[i] = pp[i];
141 /* helical shape for (minor radius>major radius *
142 * cross section p[5] same as NFUNCTION = 6 */
143 th = p[1] * p[4];
144 ph = cos(th);
145 th = sin(th);
146 x2 = p[0] - p[6] * ph;
147 z2 = p[2] - p[6] * th;
148 p[0] = x2 * ph + z2 * th;
149 p[2] = (-x2 * th + z2 * ph);
150 if (p[8] == 1.0)
151 return (sqrt(p[0] * p[0] + p[2] * p[2]) - p[5]);
152 if (p[8] != 0.0)
153 {
154 temp = 2. / p[8];
155 r2 = pow((pow(fabs(p[0]), temp) + pow(fabs(p[2]), temp)), p[8] *0.5);
156 }
157 else
158 r2 = maxim(fabs(p[0]), fabs(p[2]));
159
160 return (r2 - p[5]);
161 }
162