← Back to branch summary

~ubuntu-branches/ubuntu/vivid/libctl/vivid

« back to all changes in this revision

Viewing changes to utils/geom.c

  • Committer: Package Import Robot
  • Author(s): Thorsten Alteholz
  • Date: 2012-07-24 18:00:26 UTC
  • mto: (11.1.1 experimental) (1.3.1)
  • mto: This revision was merged to the branch mainline in revision 12.
  • Revision ID: package-import@ubuntu.com-20120724180026-d5yjvlshxjrvu1kk
Tags: upstream-3.2
Import upstream version 3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
utils/geom.c
1
1
/* libctl: flexible Guile-based control files for scientific software 
2
 
 * Copyright (C) 1998-2009, Steven G. Johnson
 
2
 * Copyright (C) 1998-2012 Massachusetts Institute of Technology and Steven G. Johnson
3
3
 *
4
4
 * This library is free software; you can redistribute it and/or
5
5
 * modify it under the terms of the GNU Lesser General Public
53
53
#  define FREE1(p) free(p)
54
54
#endif
55
55
 
 
56
#define K_PI 3.14159265358979323846
 
57
 
56
58
/**************************************************************************/
57
59
 
58
60
/* If v is a vector in the lattice basis, normalize v so that
67
69
          *v);
68
70
}
69
71
 
 
72
static vector3 lattice_to_cartesian(vector3 v)
 
73
{
 
74
     return matrix3x3_vector3_mult(geometry_lattice.basis, v);
 
75
}
 
76
 
 
77
static vector3 cartesian_to_lattice(vector3 v)
 
78
{
 
79
     return matrix3x3_vector3_mult(matrix3x3_inverse(geometry_lattice.basis), 
 
80
                                   v);
 
81
}
 
82
 
70
83
/* "Fix" the parameters of the given object to account for the
71
84
   geometry_lattice basis, which may be non-orthogonal.  In particular,
72
85
   this means that the normalization of several unit vectors, such
80
93
     switch(o.which_subclass) {
81
94
         case GEOM CYLINDER:
82
95
              lattice_normalize(&o.subclass.cylinder_data->axis);
 
96
              if (o.subclass.cylinder_data->which_subclass == CYL WEDGE) {
 
97
                   vector3 a = o.subclass.cylinder_data->axis;
 
98
                   vector3 s = o.subclass.cylinder_data->subclass.wedge_data->wedge_start;
 
99
                   double p = vector3_dot(s, matrix3x3_vector3_mult(geometry_lattice.metric, a));
 
100
                   o.subclass.cylinder_data->subclass.wedge_data->e1 = 
 
101
                        vector3_minus(s, vector3_scale(p, a));
 
102
                   lattice_normalize(&o.subclass.cylinder_data->subclass.wedge_data->e1);
 
103
                   o.subclass.cylinder_data->subclass.wedge_data->e2 = 
 
104
                        cartesian_to_lattice(
 
105
                             vector3_cross(lattice_to_cartesian(o.subclass.cylinder_data->axis),
 
106
                                           lattice_to_cartesian(o.subclass.cylinder_data->subclass.wedge_data->e1)));
 
107
              }
83
108
              break;
84
109
         case GEOM BLOCK:
85
110
         {
224
249
             radius += (proj/height + 0.5) *
225
250
                  (o->subclass.cylinder_data->subclass.cone_data->radius2
226
251
                   - radius);
 
252
        else if (o->subclass.cylinder_data->which_subclass == CYL WEDGE) {
 
253
             number x = vector3_dot(rm, o->subclass.cylinder_data->subclass.wedge_data->e1);
 
254
             number y = vector3_dot(rm, o->subclass.cylinder_data->subclass.wedge_data->e2);
 
255
             number theta = atan2(y, x);
 
256
             number wedge_angle = o->subclass.cylinder_data->subclass.wedge_data->wedge_angle;
 
257
             if (wedge_angle > 0) {
 
258
                  if (theta < 0) theta = theta + 2 * K_PI;
 
259
                  if (theta > wedge_angle) return 0;
 
260
             }
 
261
             else {
 
262
                  if (theta > 0) theta = theta - 2 * K_PI;
 
263
                  if (theta < wedge_angle) return 0;
 
264
             }
 
265
        }
227
266
        return(radius != 0.0 && vector3_dot(r,rm) - proj*proj <= radius*radius);
228
267
      }
229
268
      else
389
428
          double d3 = fabs(fabs(proj.z) - 0.5 * size.z);
390
429
          if (d1 < d2 && d1 < d3)
391
430
               return matrix3x3_row1(o.subclass.block_data->projection_matrix);
392
 
          else if (d2 < d1 && d2 < d3)
 
431
          else if (d2 < d3)
393
432
               return matrix3x3_row2(o.subclass.block_data->projection_matrix);
394
433
          else
395
434
               return matrix3x3_row3(o.subclass.block_data->projection_matrix);
573
612
     switch (o.which_subclass) {
574
613
         case GEOM CYLINDER:
575
614
              switch (o.subclass.cylinder_data->which_subclass) {
 
615
                  case CYL WEDGE:
 
616
                       printf("wedge");
 
617
                       break;
576
618
                  case CYL CONE:
577
619
                       printf("cone");
578
620
                       break;
614
656
              if (o.subclass.cylinder_data->which_subclass == CYL CONE)
615
657
                   printf("%*s     radius2 %g\n", indentby, "",
616
658
                        o.subclass.cylinder_data->subclass.cone_data->radius2);
 
659
              else if (o.subclass.cylinder_data->which_subclass == CYL WEDGE)
 
660
                   printf("%*s     wedge-theta %g\n", indentby, "",
 
661
                        o.subclass.cylinder_data->subclass.wedge_data->wedge_angle);
617
662
              break;
618
663
         case GEOM SPHERE:
619
664
              printf("%*s     radius %g\n", indentby, "", 
1182
1227
     return 0.0;
1183
1228
}
1184
1229
 
1185
 
#define K_PI 3.14159265358979323846
1186
 
 
1187
1230
number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o,
1188
1231
                           number tol, integer maxeval)
1189
1232
{
1904
1947
     return o;
1905
1948
}
1906
1949
 
 
1950
geometric_object make_wedge(material_type material, vector3 center,
 
1951
                            number radius, number height, vector3 axis,
 
1952
                            number wedge_angle, vector3 wedge_start)
 
1953
{
 
1954
     geometric_object o = make_cylinder(material, center, radius,height, axis);
 
1955
     o.subclass.cylinder_data->which_subclass = CYL WEDGE;
 
1956
     o.subclass.cylinder_data->subclass.wedge_data = MALLOC1(wedge);
 
1957
     CHECK(o.subclass.cylinder_data->subclass.wedge_data, "out of memory");
 
1958
     o.subclass.cylinder_data->subclass.wedge_data->wedge_angle = wedge_angle;
 
1959
     o.subclass.cylinder_data->subclass.wedge_data->wedge_start = wedge_start;
 
1960
     geom_fix_object(o);
 
1961
     return o;
 
1962
}
 
1963
 
1907
1964
geometric_object make_sphere(material_type material, vector3 center,
1908
1965
                             number radius)
1909
1966
{