~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/geom/mat2.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**@(#)Mat2 handling
  2 **/
  3 
  4 #include <math.h>
  5 #include <stdio.h>
  6 #include <string.h>
  7 #include <tina/sys.h>
  8 #include <tina/sysfuncs.h>
  9 #include <tina/math.h>
 10 #include <tina/mathfuncs.h>
 11 
 12 static Mat2 mat2_0 = {Mat2_id};
 13 
 14 void   *mat2_alloc(void)
 15 {
 16     Mat2   *m = ts_ralloc(Mat2);
 17 
 18     m->el[0][0] = 0.0;
 19     m->el[0][1] = 0.0;
 20     m->el[1][0] = 0.0;
 21     m->el[1][1] = 0.0;
 22     return ((void *) m);
 23 }
 24 
 25 void   *mat2_make(Mat2 n)
 26 {
 27     Mat2   *m = ts_ralloc(Mat2);
 28 
 29     *m = n;
 30     return ((void *) m);
 31 }
 32 
 33 void    mat2_free(void *m)
 34 {
 35     rfree((void *) m);
 36 }
 37 
 38 Mat2    mat2(double mxx, double mxy, double myx, double myy)
 39 {
 40     Mat2    m = {Mat2_id};
 41 
 42     m.el[0][0] = mxx;
 43     m.el[0][1] = mxy;
 44     m.el[1][0] = myx;
 45     m.el[1][1] = myy;
 46     return (m);
 47 }
 48 
 49 Mat2    mat2_unit(void)
 50 {
 51     return (mat2(1.0, 0.0,
 52                  0.0, 1.0));
 53 }
 54 
 55 Mat2    mat2_zero(void)
 56 {
 57     return (mat2_0);
 58 }
 59 
 60 void    mat2_comps(Mat2 m, float *mxx, float *mxy, float *myx, float *myy)
 61 {
 62     *mxx = m.el[0][0];
 63     *mxy = m.el[0][1];
 64     *myx = m.el[1][0];
 65     *myy = m.el[1][1];
 66 }
 67 
 68 Vec2    mat2_rowx(Mat2 m)
 69 {
 70     return (vec2(m.el[0][0], m.el[0][1]));
 71 }
 72 
 73 Vec2    mat2_rowy(Mat2 m)
 74 {
 75     return (vec2(m.el[1][0], m.el[1][1]));
 76 }
 77 
 78 Vec2    mat2_colx(Mat2 m)
 79 {
 80     return (vec2(m.el[0][0], m.el[1][0]));
 81 }
 82 
 83 Vec2    mat2_coly(Mat2 m)
 84 {
 85     return (vec2(m.el[0][1], m.el[1][1]));
 86 }
 87 
 88 Mat2    mat2_of_rows(Vec2 rx, Vec2 ry)
 89 {
 90     Mat2    m = {Mat2_id};
 91 
 92     m.el[0][0] = rx.el[0];
 93     m.el[1][0] = rx.el[1];
 94     m.el[0][1] = ry.el[0];
 95     m.el[1][1] = ry.el[1];
 96     return (m);
 97 }
 98 
 99 Mat2    mat2_of_cols(Vec2 cx, Vec2 cy)
100 {
101     Mat2    m = {Mat2_id};
102 
103     m.el[0][0] = cx.el[0];
104     m.el[0][1] = cx.el[1];
105     m.el[1][0] = cy.el[0];
106     m.el[1][1] = cy.el[1];
107     return (m);
108 }
109 
110 Mat2    mat2_sum(Mat2 m, Mat2 n)
111 {
112     Mat2    sum = {Mat2_id};
113     int     i, j;
114 
115     for (i = 0; i < 2; ++i)
116         for (j = 0; j < 2; ++j)
117             sum.el[i][j] = m.el[i][j] + n.el[i][j];
118     return (sum);
119 }
120 
121 Mat2    mat2_diff(Mat2 m, Mat2 n)
122 {
123     Mat2    diff = {Mat2_id};
124     int     i, j;
125 
126     for (i = 0; i < 2; ++i)
127         for (j = 0; j < 2; ++j)
128             diff.el[i][j] = m.el[i][j] - n.el[i][j];
129     return (diff);
130 }
131 
132 Mat2    mat2_minus(Mat2 m)
133 {
134     Mat2    minus = {Mat2_id};
135     int     i, j;
136 
137     for (i = 0; i < 2; ++i)
138         for (j = 0; j < 2; ++j)
139             minus.el[i][j] = -m.el[i][j];
140     return (minus);
141 }
142 
143 Mat2    mat2_times(double k, Mat2 m)
144 {
145     Mat2    prod = {Mat2_id};
146     int     i, j;
147 
148     for (i = 0; i < 2; ++i)
149         for (j = 0; j < 2; ++j)
150             prod.el[i][j] = k * m.el[i][j];
151     return (prod);
152 }
153 
154 Mat2    mat2_prod(Mat2 m, Mat2 n)
155 {
156     Mat2    prod = {Mat2_id};
157     int     i, j, k;
158     double  sum;
159 
160     for (i = 0; i < 2; ++i)
161         for (j = 0; j < 2; ++j)
162         {
163             sum = 0.0;
164             for (k = 0; k < 2; ++k)
165                 sum += m.el[i][k] * n.el[k][j];
166             prod.el[i][j] = sum;
167         }
168     return (prod);
169 }
170 
171 Mat2    mat2_inverse(Mat2 m)
172 {
173     double  det, k;
174     float   mxx, mxy, myx, myy;
175 
176     mat2_comps(m, &mxx, &mxy, &myx, &myy);
177     det = mxx * myy - mxy * myx;
178     if (det == 0.0)
179         return (mat2_unit());
180     k = 1.0 / det;
181     return (mat2(k * myy, -k * mxy,
182                  -k * myx, k * mxx));
183 }
184 
185 Mat2    mat2_transpose(Mat2 m)
186 {
187     Mat2    trans = {Mat2_id};
188     int     i, j;
189 
190     for (i = 0; i < 2; ++i)
191         for (j = 0; j < 2; ++j)
192             trans.el[i][j] = m.el[j][i];
193     return (trans);
194 }
195 
196 Mat2    mat2_sym(Mat2 m)
197 {
198     double  mxx = mat2_xx(m);
199     double  mxy = mat2_xy(m);
200     double  myx = mat2_yx(m);
201     double  myy = mat2_yy(m);
202 
203     mxy = 0.5 * (mxy + myx);
204     return (mat2(mxx, mxy, mxy, myy));
205 }
206 
207 double  mat2_trace(Mat2 m)
208 {
209     return (m.el[0][0] + m.el[1][1]);
210 }
211 
212 double  mat2_det(Mat2 m)
213 {
214     return (mat2_xx(m) * mat2_yy(m) - mat2_xy(m) * mat2_yx(m));
215 }
216 
217 Bool    mat2_posdef(Mat2 m)
218 {
219     if (mat2_xx(m) <= 0.0)
220         return (false);
221     else if (mat2_det(m) <= 0.0)
222         return (false);
223     else
224         return (true);
225 }
226 
227 Vec2    mat2_vprod(Mat2 m, Vec2 v)
228 {
229     Vec2    prod = {Vec2_id};
230     int     i, j;
231     double  sum;
232 
233     for (i = 0; i < 2; ++i)
234     {
235         sum = 0.0;
236         for (j = 0; j < 2; ++j)
237             sum += m.el[i][j] * v.el[j];
238         prod.el[i] = sum;
239     }
240     return (prod);
241 }
242 
243 double  mat2_sprod(Vec2 v, Mat2 m, Vec2 w)
244 {
245     double  sum = 0.0;
246     int     i, j;
247 
248     for (i = 0; i < 2; ++i)
249         for (j = 0; j < 2; ++j)
250             sum += v.el[i] * m.el[i][j] * w.el[j];
251     return (sum);
252 }
253 
254 Mat2    mat2_tensor(Vec2 v, Vec2 w)
255 {
256     Mat2    prod = {Mat2_id};
257     int     i, j;
258 
259     for (i = 0; i < 2; ++i)
260         for (j = 0; j < 2; ++j)
261             prod.el[i][j] = v.el[i] * w.el[j];
262     return (prod);
263 }
264 
265 Mat2    mat2_read(FILE * fp)
266 {
267     Mat2    temp = {Mat2_id};
268     int     i;
269 
270     for (i = 0; i < 2; ++i)
271         (void) fscanf(fp, "%f %f", &temp.el[i][0], &temp.el[i][1]);
272     return (temp);
273 }
274 
275 void    mat2_print(FILE * fp, Mat2 m)
276 {
277     int     i;
278 
279     for (i = 0; i < 2; ++i)
280         (void) fprintf(fp, "%f %f\n", m.el[i][0], m.el[i][1]);
281 }
282 
283 void    mat2_pprint(FILE * fp, char *msg, Mat2 m)
284 {
285     int     i, n;
286 
287     n = strlen(msg);
288     (void) fprintf(fp, "%s", msg);
289     (void) fprintf(fp, "|%15.6f%15.6f|\n", m.el[0][0], m.el[0][1]);
290     for (i = 0; i < n; ++i)
291         (void) fputc(' ', stderr);
292     (void) fprintf(fp, "|%15.6f%15.6f|\n", m.el[1][0], m.el[1][1]);
293 }
294 
295 void    mat2_format(Mat2 m)
296 {
297     int     i;
298 
299     for (i = 0; i < 2; ++i)
300         format("%f %f \n", m.el[i][0], m.el[i][1]);
301 }
302 
303 /**
304 eigendirection and eigenvalues of symmetric 2by2 matrix
305 |a b|
306 |b c|
307 **/
308 
309 void    mat2_sym_eigen(Mat2 m, double *theta, double *lambda1, double *lambda2)
310 {
311     double  a, b, c, p, q, r;
312 
313     a = mat2_xx(m);
314     b = 0.5 * (mat2_xy(m) + mat2_yx(m));
315     c = mat2_yy(m);
316     p = (a + c) / 2.0;
317     q = (a - c) / 2.0;
318     r = sqrt(q * q + b * b);
319     if (b == 0.0 && q == 0.0)
320         *theta = 0.0;
321     else
322         *theta = 0.5 * atan2(b, q);
323     *lambda1 = p + r;
324     *lambda2 = p - r;
325 }
326 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.