File: gaussmap.cmd

package info (click to toggle)
evolver 2.70+ds-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 17,148 kB
  • sloc: ansic: 127,395; makefile: 209; sh: 98
file content (71 lines) | stat: -rw-r--r-- 2,461 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
// gaussmap.cmd
// Converts each vertex coordinate to unit normal.
// No vertices should be on constraints or fixed.
// Saves original coordinates so "revert" restores original surface.
// Before running "gaussmap", vertices should be freed from all
// constraints and boundaries.

// Programmer: Ken Brakke, brakke@susqu.edu, http://www.susqu.edu/brakke

// Usage: gaussmap

define vertex attribute gaussx real[3]
define vertex attribute oldx real[3]
gaussmap := { 
   foreach vertex vv do  
   { local connum;
     connum := 1;
     do { unset vv constraint connum; connum += 1; }
     while ( connum < 30 );
     vv.oldx[1] := vv.x[1];
     vv.oldx[2] := vv.x[2];
     vv.oldx[3] := vv.x[3];
     vv.gaussx[1] := vv.vertexnormal[1];
     vv.gaussx[2] := vv.vertexnormal[2];
     vv.gaussx[3] := vv.vertexnormal[3];
   };
   set vertex x gaussx[1];
   set vertex y gaussx[2];
   set vertex z gaussx[3];
}

revert := {
   set vertex x oldx[1];
   set vertex y oldx[2];
   set vertex z oldx[3];
}

// Command to find spherical area of gauss map, Lagrange model
calc_gaussarea := {
   local jnx, jnxx;
   if quadratic then lagrange 2; 
   gaussarea := 0;
   jnx := 2;
   jnxx := 3;
   foreach facet ff do
   { local angle_a, angle_b, angle_c;
     angle_a := asin(sqrt(
        (ff.vertex[1].vertexnormal[1]-ff.vertex[jnx].vertexnormal[1])^2
       +(ff.vertex[1].vertexnormal[2]-ff.vertex[jnx].vertexnormal[2])^2
       +(ff.vertex[1].vertexnormal[3]-ff.vertex[jnx].vertexnormal[3])^2)/2);
     angle_b := asin(sqrt(
        (ff.vertex[jnxx].vertexnormal[1]-ff.vertex[jnx].vertexnormal[1])^2
       +(ff.vertex[jnxx].vertexnormal[2]-ff.vertex[jnx].vertexnormal[2])^2
       +(ff.vertex[jnxx].vertexnormal[3]-ff.vertex[jnx].vertexnormal[3])^2)/2);
     angle_c := asin(sqrt(
        (ff.vertex[jnxx].vertexnormal[1]-ff.vertex[1].vertexnormal[1])^2
       +(ff.vertex[jnxx].vertexnormal[2]-ff.vertex[1].vertexnormal[2])^2
       +(ff.vertex[jnxx].vertexnormal[3]-ff.vertex[1].vertexnormal[3])^2)/2);
     gaussarea += 4*atan(sqrt(tan((angle_a+angle_b+angle_c)/2)
                         *tan((angle_a+angle_b-angle_c)/2)
                         *tan((angle_a-angle_b+angle_c)/2)
                         *tan((-angle_a+angle_b+angle_c)/2)));
   };
   printf "Gauss map area: %18.15f or %18.15f*pi\n",gaussarea,gaussarea/pi;
}
     
// End gaussmap.cmd

// Usage: gaussmap
//        revert