File: slice.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 (64 lines) | stat: -rw-r--r-- 2,202 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
// slice.cmd --- Calculate length of intersection of plane with surface.
// Does not modify surface.
// plane eq: aa*x + bb*y + cc*z = dd
// output: prints length of slice, and area inside slice.
// Note all area inside slice is counted as positive!
// Try not to slice exactly through vertices!!

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

// Usage: Set plane coefficients, then do "slice".
// Results: Intersection length printed out, left in variable lensum.

slice_aa := 0; slice_bb := 0; slice_cc := 1; slice_dd := .1;  // set these for desired plane
slice := { 
       local any,xx1,yy1,zz1,xx2,yy2,zz2;
       local denom,lambda,zaa,zbb,darea;
       local xa,ya,za,xb,yb,zb,dx,dy,dz;

       lensum := 0; areasum := 0;
       foreach facet ff do 
       { any := 0; 
         foreach ff.edge ee do 
         {
           xx1 := ee.vertex[1].x; 
           yy1 := ee.vertex[1].y; 
           zz1 := ee.vertex[1].z; 
           xx2 := ee.vertex[2].x; 
           yy2 := ee.vertex[2].y; 
           zz2 := ee.vertex[2].z;
           denom := slice_aa*(xx1-xx2)+slice_bb*(yy1-yy2)+slice_cc*(zz1-zz2);
           if ( denom != 0.0 ) then 
           { 
             lambda := (slice_dd-slice_aa*xx2-slice_bb*yy2-slice_cc*zz2)/denom; 
             if ( (lambda >= 0) and (lambda <= 1) ) then 
             { 
               xa := xb; ya := yb; za := zb;
               xb := lambda*xx1+(1-lambda)*xx2; 
               yb := lambda*yy1+(1-lambda)*yy2;
               zb := lambda*zz1+(1-lambda)*zz2; 
               any := any+1; 
             } 
           }
        } ; 
        if any == 2 then
        { 
          dx := xa-xb; dy := ya-yb; dz := za - zb;
          lensum := lensum + sqrt(dx^2+dy^2+dz^2);
          zaa := za - slice_dd/slice_cc; zbb := zb - slice_dd/slice_cc;
          darea := sqrt((ya*zbb-yb*zaa)^2+(zaa*xb-zbb*xa)^2+(xa*yb-ya*xb)^2);
          areasum := areasum + darea/2;
        }
      };
      printf "Circumference: %18.15g  Area: %18.15g\n",lensum,areasum;
    } // end slice

// end slice.cmd

/* Usage:
   set  slice_aa, slice_bb, slice_cc, slice_dd
   call slice
*/