"Fossies" - the Fresh Open Source Software Archive

Member "mathmod-branches-r508-trunk/pariso/isosurface/internalfunctions.cpp" (8 Mar 2021, 5121 Bytes) of package /linux/misc/mathmod-11.0-source.zip:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "internalfunctions.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 10.1_vs_11.0.

    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 
   25 double maxim(double p1, double p2)
   26 {
   27     return p1 > p2 ? p1 : p2;
   28 }
   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 }
   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 }
   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 }
   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 }
  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 }