NAMD
TclVec.C
Go to the documentation of this file.
1 /***************************************************************************
2  *cr
3  *cr (C) Copyright 1995-2008 The Board of Trustees of the
4  *cr University of Illinois
5  *cr All Rights Reserved
6  *cr
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * RCS INFORMATION:
11  *
12  * $RCSfile: TclVec.C,v $
13  * $Author: jim $ $Locker: $ $State: Exp $
14  * $Revision: 1.2 $ $Date: 2008/09/17 16:19:54 $
15  *
16  ***************************************************************************
17  * DESCRIPTION:
18  * A C-based implementation of some performance-critical Tcl callable
19  * routines in VMD. The C implementation outperforms a raw Tcl version
20  * by a factor of three or so. The performance advantage helps
21  * significantly when doing analysis in VMD.
22  ***************************************************************************/
23 
24 #include <tcl.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #include <math.h>
28 // #include "TclCommands.h"
29 // #include "Matrix4.h"
30 // #include "utilities.h"
31 
32 #define VMD_PI 3.14159265358979323846
33 #define VMD_TWOPI (2.0 * VMD_PI)
34 #define DEGTORAD(a) (a*VMD_PI/180.0)
35 #define RADTODEG(a) (a*180.0/VMD_PI)
36 
37 /***************** override some of the vector routines for speed ******/
38 /* These should be the exact C equivalent to the corresponding Tcl */
39 /* vector commands */
40 
41 // Function: vecadd v1 v2 {v3 ...}
42 // Returns: the sum of vectors; must all be the same length
43 // The increase in speed from Tcl to C++ is 4561 / 255 == 18 fold
44 static int obj_vecadd(ClientData, Tcl_Interp *interp, int argc,
45  Tcl_Obj * const objv[]) {
46  if (argc < 3) {
47  Tcl_WrongNumArgs(interp, 1, objv, (char *)"vec1 vec2 ?vec3? ?vec4? ...");
48  return TCL_ERROR;
49  }
50  int num;
51  Tcl_Obj **data;
52  if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) {
53  return TCL_ERROR;
54  }
55  double *sum = new double[num];
56  int i;
57  for (i=0; i<num; i++) {
58  if (Tcl_GetDoubleFromObj(interp, data[i], sum+i) != TCL_OK) {
59  delete [] sum;
60  return TCL_ERROR;
61  }
62  }
63  // do the sums on the rest
64  int num2;
65  for (int term=2; term < argc; term++) {
66  if (Tcl_ListObjGetElements(interp, objv[term], &num2, &data) != TCL_OK) {
67  delete [] sum;
68  return TCL_ERROR;
69  }
70  if (num != num2) {
71  Tcl_SetResult(interp, (char *) "vecadd: two vectors don't have the same size", TCL_STATIC);
72  delete [] sum;
73  return TCL_ERROR;
74  }
75  for (i=0; i<num; i++) {
76  double df;
77  if (Tcl_GetDoubleFromObj(interp, data[i], &df) != TCL_OK) {
78  delete [] sum;
79  return TCL_ERROR;
80  }
81  sum[i] += df;
82  }
83  }
84 
85 
86  // and return the result
87  Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
88  for (i=0; i<num; i++) {
89  Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(sum[i]));
90  }
91  Tcl_SetObjResult(interp, tcl_result);
92  delete [] sum;
93  return TCL_OK;
94 }
95 
96 // Function: vecsub v1 v2
97 // Returns: v1 - v2
98 
99 
100 static int obj_vecsub(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
101 {
102  if (argc != 3) {
103  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?x? ?y?");
104  return TCL_ERROR;
105  }
106  int num1=0, num2=0;
107  Tcl_Obj **data1, **data2;
108  if (Tcl_ListObjGetElements(interp, objv[1], &num1, &data1) != TCL_OK)
109  return TCL_ERROR;
110  if (Tcl_ListObjGetElements(interp, objv[2], &num2, &data2) != TCL_OK)
111  return TCL_ERROR;
112 
113  if (num1 != num2) {
114  Tcl_SetResult(interp, (char *)"vecsub: two vectors don't have the same size", TCL_STATIC);
115  return TCL_ERROR;
116  }
117 
118  Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
119  for (int i=0; i<num1; i++) {
120  double d1=0, d2=0;
121  if (Tcl_GetDoubleFromObj(interp, data1[i], &d1) != TCL_OK) {
122  Tcl_SetResult(interp, (char *)"vecsub: non-numeric in first argument", TCL_STATIC);
123  return TCL_ERROR;
124  }
125  if (Tcl_GetDoubleFromObj(interp, data2[i], &d2) != TCL_OK) {
126  Tcl_SetResult(interp, (char *)"vecsub: non-numeric in second argument", TCL_STATIC);
127  return TCL_ERROR;
128  }
129  Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(d1-d2));
130  }
131  Tcl_SetObjResult(interp, tcl_result);
132  return TCL_OK;
133 }
134 
135 
136 // Function: vecscale
137 // Returns: scalar * vector or vector * scalar
138 // speedup is 1228/225 = 5.5 fold
139 static int obj_vecscale(ClientData, Tcl_Interp *interp, int argc,
140  Tcl_Obj * const objv[]) {
141  if (argc != 3) {
142  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?c? ?v?");
143  return TCL_ERROR;
144  }
145 
146  int num1, num2;
147  Tcl_Obj **data1, **data2;
148  if (Tcl_ListObjGetElements(interp, objv[1], &num1, &data1) != TCL_OK) {
149  return TCL_ERROR;
150  }
151  if (Tcl_ListObjGetElements(interp, objv[2], &num2, &data2) != TCL_OK) {
152  return TCL_ERROR;
153  }
154  if (num1 == 0 || num2 == 0) {
155  Tcl_SetResult(interp, (char *) "vecscale: parameters must have data", TCL_STATIC);
156  return TCL_ERROR;
157  } else if (num1 != 1 && num2 != 1) {
158  Tcl_SetResult(interp, (char *) "vecscale: one parameter must be a scalar value", TCL_STATIC);
159  return TCL_ERROR;
160  }
161 
162  int num;
163  Tcl_Obj *scalarobj, **vector;
164  if (num1 == 1) {
165  scalarobj = data1[0];
166  vector = data2;
167  num = num2;
168  } else {
169  scalarobj = data2[0];
170  vector = data1;
171  num = num1;
172  }
173 
174  double scalar;
175  if (Tcl_GetDoubleFromObj(interp, scalarobj, &scalar) != TCL_OK)
176  return TCL_ERROR;
177 
178  Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
179  for (int i=0; i<num; i++) {
180  double val;
181  if (Tcl_GetDoubleFromObj(interp, vector[i], &val) != TCL_OK) {
182  Tcl_SetResult(interp, (char *) "vecscale: non-numeric in vector", TCL_STATIC);
183  return TCL_ERROR;
184  }
185  val *= scalar;
186  Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(val));
187  }
188  Tcl_SetObjResult(interp, tcl_result);
189  return TCL_OK;
190 }
191 
193 // returns TCL_OK if good
194 // If bad, returns TCL_ERROR and sets the Tcl result to the error message
195 // The name of the function should be passed in 'fctn' so the error message
196 // can be constructed correctly
197 int tcl_get_matrix(const char *fctn, Tcl_Interp *interp,
198  Tcl_Obj *s, double *mat)
199 {
200  int num_rows;
201  Tcl_Obj **data_rows;
202  if (Tcl_ListObjGetElements(interp, s, &num_rows, &data_rows) != TCL_OK) {
203  char tmpstring[1024];
204  sprintf(tmpstring, "%s: badly formed matrix", fctn);
205  Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
206  return TCL_ERROR;
207  }
208  if (num_rows != 4) {
209  char tmpstring[1024];
210  sprintf(tmpstring, "%s: need a 4x4 matrix", fctn);
211  Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
212  return TCL_ERROR;
213  }
214  int num_row[4];
215  Tcl_Obj **data_row[4];
216  if (Tcl_ListObjGetElements(interp, data_rows[0], num_row+0, data_row+0) != TCL_OK ||
217  num_row[0] != 4 ||
218  Tcl_ListObjGetElements(interp, data_rows[1], num_row+1, data_row+1) != TCL_OK ||
219  num_row[1] != 4 ||
220  Tcl_ListObjGetElements(interp, data_rows[2], num_row+2, data_row+2) != TCL_OK ||
221  num_row[2] != 4 ||
222  Tcl_ListObjGetElements(interp, data_rows[3], num_row+3, data_row+3) != TCL_OK ||
223  num_row[3] != 4) {
224  Tcl_AppendResult(interp, fctn, ": poorly formed matrix", NULL);
225  return TCL_ERROR;
226  }
227  // now get the numbers
228  for (int i=0; i<4; i++) {
229  for (int j=0; j<4; j++) {
230  double tmp;
231  if (Tcl_GetDoubleFromObj(interp, data_row[i][j], &tmp) != TCL_OK) {
232  char tmpstring[1024];
233  sprintf(tmpstring, "%s: non-numeric in matrix", fctn);
234  Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
235  return TCL_ERROR;
236  } else {
237  mat[4*j+i] = (double) tmp; // Matrix4 is transpose of Tcl's matrix
238  }
239  }
240  }
241  return TCL_OK;
242 }
243 
244 #if 0
245 int tcl_get_vector(const char *s, double *val, Tcl_Interp *interp)
246 {
247  int num;
248  const char **pos;
249  if (Tcl_SplitList(interp, s, &num, &pos) != TCL_OK) {
250  Tcl_SetResult(interp, (char *) "need three data elements for a vector",
251  TCL_STATIC);
252  return TCL_ERROR;
253  }
254  if (num != 3) {
255  Tcl_SetResult(interp, (char *) "need three numbers for a vector", TCL_STATIC);
256  return TCL_ERROR;
257  }
258  double a[3];
259  if (Tcl_GetDouble(interp, pos[0], a+0) != TCL_OK ||
260  Tcl_GetDouble(interp, pos[1], a+1) != TCL_OK ||
261  Tcl_GetDouble(interp, pos[2], a+2) != TCL_OK) {
262  ckfree((char *) pos); // free of tcl data
263  return TCL_ERROR;
264  }
265  val[0] = (double) a[0];
266  val[1] = (double) a[1];
267  val[2] = (double) a[2];
268  ckfree((char *) pos); // free of tcl data
269  return TCL_OK;
270 }
271 #endif
272 
273 // append the matrix into the Tcl result
274 void tcl_append_matrix(Tcl_Interp *interp, const double *mat) {
275  Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
276  for (int i=0; i<4; i++) {
277  Tcl_Obj *m = Tcl_NewListObj(0, NULL);
278  for (int j=0; j<4; j++)
279  Tcl_ListObjAppendElement(interp, m, Tcl_NewDoubleObj(mat[4*j+i]));
280  Tcl_ListObjAppendElement(interp, tcl_result, m);
281  }
282  Tcl_SetObjResult(interp, tcl_result);
283 }
284 
285 // speed up the matrix * vector routines -- DIFFERENT ERROR MESSAGES
286 // THAN THE TCL VERSION
287 // speedup is nearly 25 fold
288 static int obj_vectrans(ClientData, Tcl_Interp *interp, int argc,
289  Tcl_Obj * const objv[])
290 {
291  if (argc != 3) {
292  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?matrix? ?vector?");
293  return TCL_ERROR;
294  }
295 
296  // get the matrix data
297  double mat[16];
298  if (tcl_get_matrix(
299  Tcl_GetStringFromObj(objv[0],NULL), interp, objv[1], mat) != TCL_OK) {
300  return TCL_ERROR;
301  }
302 
303  // for the vector
304  Tcl_Obj **vec;
305  int vec_size;
306  if (Tcl_ListObjGetElements(interp, objv[2], &vec_size, &vec) != TCL_OK)
307  return TCL_ERROR;
308 
309  if (vec_size != 3 && vec_size != 4) {
310  Tcl_SetResult(interp, (char *) "vectrans: vector must be of size 3 or 4",
311  TCL_STATIC);
312  return TCL_ERROR;
313  }
314 
315  double opoint[4];
316  opoint[3] = 0;
317  for (int i=0; i<vec_size; i++) {
318  double tmp;
319  if (Tcl_GetDoubleFromObj(interp, vec[i], &tmp) != TCL_OK) {
320  Tcl_SetResult(interp, (char *) "vectrans: non-numeric in vector", TCL_STATIC);
321  return TCL_ERROR;
322  }
323  opoint[i] = (double)tmp;
324  }
325  // vector data is in vec_data
326  double npoint[4];
327 
328  npoint[0]=opoint[0]*mat[0]+opoint[1]*mat[4]+opoint[2]*mat[8]+opoint[3]*mat[12]
329 ;
330  npoint[1]=opoint[0]*mat[1]+opoint[1]*mat[5]+opoint[2]*mat[9]+opoint[3]*mat[13]
331 ;
332  npoint[2]=opoint[0]*mat[2]+opoint[1]*mat[6]+opoint[2]*mat[10]+opoint[3]*mat[14
333 ];
334  npoint[3]=opoint[0]*mat[3]+opoint[1]*mat[7]+opoint[2]*mat[11]+opoint[3]*mat[15
335 ];
336  // return it
337 
338  {
339  Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
340  for (int i=0; i<vec_size; i++)
341  Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(npoint[i]));
342  Tcl_SetObjResult(interp, tcl_result);
343  }
344  return TCL_OK;
345 }
346 
347 
348 // Function: transmult m1 m2 ... mn
349 // Returns: the product of the matricies
350 // speedup is 136347 / 1316 = factor of 104
351 static int obj_transmult(ClientData, Tcl_Interp *interp, int argc,
352  Tcl_Obj * const objv[]) {
353  // make there there are at least two values
354  if (argc < 3) {
355  Tcl_WrongNumArgs(interp, 1, objv, (char *)"mx my ?m1? ?m2? ...");
356  return TCL_ERROR;
357  }
358  // Get the first matrix
359  double mult[16];
360  if (tcl_get_matrix("transmult: ", interp, objv[1], mult) != TCL_OK) {
361  return TCL_ERROR;
362  }
363  int i = 2;
364  double pre[16];
365  while (i < argc) {
366  if (tcl_get_matrix("transmult: ", interp, objv[i], pre) != TCL_OK) {
367  return TCL_ERROR;
368  }
369  // premultiply mult by tmp
370  double tmp[4];
371  for (int k=0; k<4; k++) {
372  tmp[0] = mult[k];
373  tmp[1] = mult[4+k];
374  tmp[2] = mult[8+k];
375  tmp[3] = mult[12+k];
376  for (int j=0; j<4; j++) {
377  mult[4*j+k] = pre[4*j]*tmp[0] + pre[4*j+1]*tmp[1] +
378  pre[4*j+2]*tmp[2] + pre[4*j+3]*tmp[3];
379  }
380  }
381  i++;
382  }
383  tcl_append_matrix(interp, mult);
384  return TCL_OK;
385 }
386 
387 static int obj_transvec(ClientData, Tcl_Interp *interp, int argc,
388  Tcl_Obj * const objv[]) {
389  if (argc != 2) {
390  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
391  return TCL_ERROR;
392  }
393 
394  int num;
395  Tcl_Obj **data;
396  if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK)
397  return TCL_ERROR;
398  if (num != 3) {
399  Tcl_AppendResult(interp, "transvec: vector must have three elements",NULL);
400  return TCL_ERROR;
401  }
402  double x,y,z;
403  if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK ||
404  Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK ||
405  Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) {
406  Tcl_SetResult(interp, (char *)"transvec: non-numeric in vector", TCL_STATIC);
407  return TCL_ERROR;
408  }
409  Matrix4 mat;
410  mat.transvec((double) x,(double) y,(double) z);
411  tcl_append_matrix(interp, mat.mat);
412  return TCL_OK;
413 }
414 
415 static int obj_transvecinv(ClientData, Tcl_Interp *interp, int argc,
416  Tcl_Obj * const objv[]) {
417  if (argc != 2) {
418  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
419  return TCL_ERROR;
420  }
421 
422  int num;
423  Tcl_Obj **data;
424  if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK)
425  return TCL_ERROR;
426  if (num != 3) {
427  Tcl_AppendResult(interp, "transvecinv: vector must have three elements",NULL);
428  return TCL_ERROR;
429  }
430  double x,y,z;
431  if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK ||
432  Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK ||
433  Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) {
434  Tcl_SetResult(interp, (char *)"transvecinv: non-numeric in vector", TCL_STATIC);
435  return TCL_ERROR;
436  }
437  Matrix4 mat;
438  mat.transvecinv((double) x,(double) y,(double) z);
439  tcl_append_matrix(interp, mat.mat);
440  return TCL_OK;
441 }
442 
443 // Returns the transformation matrix needed to rotate by a certain
444 // angle around a given axis.
445 // Tcl syntax:
446 // transabout v amount [deg|rad|pi]
447 // The increase in speed from Tcl to C++ is 15 fold
448 static int obj_transabout(ClientData, Tcl_Interp *interp, int argc,
449  Tcl_Obj * const objv[]) {
450  if (argc != 3 && argc != 4) {
451  Tcl_WrongNumArgs(interp, 1, objv, (char *)"axis amount [deg|rad|pi]");
452  return TCL_ERROR;
453  }
454 
455  int num;
456  Tcl_Obj **data;
457  // get the axis
458  if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK)
459  return TCL_ERROR;
460  if (num != 3) {
461  Tcl_AppendResult(interp, "transabout: vector must have three elements",NULL);
462  return TCL_ERROR;
463  }
464  double x,y,z;
465  if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK ||
466  Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK ||
467  Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) {
468  Tcl_SetResult(interp, (char *)"transabout: non-numeric in vector", TCL_STATIC);
469  return TCL_ERROR;
470  }
471 
472  // get the amount
473  double amount;
474  if (Tcl_GetDoubleFromObj(interp, objv[2], &amount) != TCL_OK) {
475  Tcl_SetResult(interp, (char *)"transabout: non-numeric angle", TCL_STATIC);
476  return TCL_ERROR;
477  }
478 
479  // get units
480  if (argc == 4) {
481  if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "deg")) {
482  amount = DEGTORAD(amount);
483  } else if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "rad")) {
484  // amount = amount;
485  } else if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "pi")) {
486  amount = amount*VMD_PI;
487  } else {
488  Tcl_AppendResult(interp, "transabout: unit must be deg|rad|pi",NULL);
489  return TCL_ERROR;
490  }
491  } else {
492  // If no unit was specified assume that we have degrees
493  amount = DEGTORAD(amount);
494  }
495 
496  double axis[3];
497  axis[0] = (double) x;
498  axis[1] = (double) y;
499  axis[2] = (double) z;
500 
501  Matrix4 mat;
502  mat.rotate_axis(axis, amount);
503  tcl_append_matrix(interp, mat.mat);
504  return TCL_OK;
505 }
506 
507 static int obj_veclength(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) {
508 
509  if (argc != 2) {
510  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
511  return TCL_ERROR;
512  }
513 
514  int num;
515  Tcl_Obj **data;
516  if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK)
517  return TCL_ERROR;
518 
519  double length = 0.;
520  for (int i=0; i<num; i++) {
521  double tmp;
522  if (Tcl_GetDoubleFromObj(interp, data[i], &tmp) != TCL_OK) {
523  Tcl_SetResult(interp, (char *) "veclength: non-numeric in vector", TCL_STATIC);
524  return TCL_ERROR;
525  } else {
526  length += tmp*tmp;
527  }
528  }
529 
530  length = sqrt(length);
531  Tcl_Obj *tcl_result = Tcl_GetObjResult(interp);
532  Tcl_SetDoubleObj(tcl_result, length);
533  return TCL_OK;
534 }
535 
536 
537 static double* obj_getdoublearray(Tcl_Interp *interp, Tcl_Obj *const objv[], int *len) {
538  int num;
539 
540  Tcl_Obj **data;
541  if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK)
542  return NULL;
543 
544  double *list = (double*) malloc(num*sizeof(double));
545  if (list == NULL)
546  return NULL;
547 
548  for (int i=0; i<num; i++) {
549  double tmp;
550  if (Tcl_GetDoubleFromObj(interp, data[i], &tmp) != TCL_OK) {
551  Tcl_SetResult(interp, (char *) "veclength: non-numeric in vector", TCL_STATIC);
552  free(list);
553  return NULL;
554  }
555  list[i] = tmp;
556  }
557 
558  *len = num;
559 
560  return list;
561 }
562 
563 
564 static int obj_vecsum(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) {
565  if (argc != 2) {
566  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
567  return TCL_ERROR;
568  }
569 
570  int num;
571  double *list = obj_getdoublearray(interp, objv, &num);
572  if (list == NULL)
573  return TCL_ERROR;
574 
575  double sum = 0.;
576  for (int i=0; i<num; i++) {
577  sum += list[i];
578  }
579  free(list);
580 
581  Tcl_Obj *tcl_result = Tcl_GetObjResult(interp);
582  Tcl_SetDoubleObj(tcl_result, sum);
583  return TCL_OK;
584 }
585 
586 
587 static int obj_vecmean(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) {
588  if (argc != 2) {
589  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
590  return TCL_ERROR;
591  }
592 
593  int num;
594  double *list = obj_getdoublearray(interp, objv, &num);
595  if (list == NULL)
596  return TCL_ERROR;
597 
598  double sum = 0.;
599  for (int i=0; i<num; i++) {
600  sum += list[i];
601  }
602  sum /= (double) num;
603  free(list);
604 
605  Tcl_Obj *tcl_result = Tcl_GetObjResult(interp);
606  Tcl_SetDoubleObj(tcl_result, sum);
607  return TCL_OK;
608 }
609 
610 
611 static int obj_vecstddev(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) {
612  if (argc != 2) {
613  Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
614  return TCL_ERROR;
615  }
616 
617  int i, num;
618  double* list = obj_getdoublearray(interp, objv, &num);
619  if (list == NULL)
620  return TCL_ERROR;
621 
622  double mean = 0.;
623  for (i=0; i<num; i++) {
624  mean += list[i];
625  }
626  mean /= (double) num;
627 
628  double stddev = 0.;
629  for (i=0; i<num; i++) {
630  double tmp = list[i] - mean;
631  stddev += tmp * tmp;
632  }
633  stddev /= (double) num;
634  stddev = sqrt(stddev);
635  free(list);
636 
637  Tcl_Obj *tcl_result = Tcl_GetObjResult(interp);
638  Tcl_SetDoubleObj(tcl_result, stddev);
639  return TCL_OK;
640 }
641 
642 
643 int Vec_Init(Tcl_Interp *interp) {
644  Tcl_CreateObjCommand(interp, "vecadd", obj_vecadd,
645  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
646  Tcl_CreateObjCommand(interp, "vecsub", obj_vecsub,
647  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
648  Tcl_CreateObjCommand(interp, "vecscale", obj_vecscale,
649  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
650  Tcl_CreateObjCommand(interp, "transmult", obj_transmult,
651  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
652  Tcl_CreateObjCommand(interp, "vectrans", obj_vectrans,
653  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
654  Tcl_CreateObjCommand(interp, "veclength", obj_veclength,
655  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
656  Tcl_CreateObjCommand(interp, "vecmean", obj_vecmean,
657  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
658  Tcl_CreateObjCommand(interp, "vecstddev", obj_vecstddev,
659  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
660  Tcl_CreateObjCommand(interp, "vecsum", obj_vecsum,
661  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
662  Tcl_CreateObjCommand(interp, "transvec", obj_transvec,
663  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
664  Tcl_CreateObjCommand(interp, "transvecinv", obj_transvecinv,
665  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
666  Tcl_CreateObjCommand(interp, "transabout", obj_transabout,
667  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
668  return TCL_OK;
669 }
670 
int Vec_Init(Tcl_Interp *interp)
Definition: TclVec.C:643
static int obj_vectrans(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:288
void transvecinv(double x, double y, double z)
apply a rotation such that the given vector is brought along &#39;x&#39;.
Definition: Matrix4.C:261
void tcl_append_matrix(Tcl_Interp *interp, const double *mat)
Definition: TclVec.C:274
static int obj_vecmean(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:587
static int obj_transvecinv(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:415
#define DEGTORAD(a)
Definition: TclVec.C:34
static int obj_transmult(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:351
4x4 matrix class with numerous operators, conversions, etc.
Definition: Matrix4.h:26
static int obj_vecscale(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:139
static int obj_vecsum(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:564
gridSize z
static int obj_vecadd(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:44
static int obj_vecstddev(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:611
#define VMD_PI
Definition: TclVec.C:32
int tcl_get_matrix(const char *fctn, Tcl_Interp *interp, Tcl_Obj *s, double *mat)
Given a string with a matrix in it, return the matrix.
Definition: TclVec.C:197
static int obj_transvec(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:387
void rotate_axis(const double axis[3], double angle)
apply a rotation around the given vector; angle in radians.
Definition: Matrix4.C:245
static int obj_transabout(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:448
double mat[16]
the matrix itself
Definition: Matrix4.h:33
void transvec(double x, double y, double z)
apply a rotation such that &#39;x&#39; is brought along the given vector.
Definition: Matrix4.C:252
static int obj_vecsub(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:100
static int obj_veclength(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
Definition: TclVec.C:507
gridSize y
gridSize x
static double * obj_getdoublearray(Tcl_Interp *interp, Tcl_Obj *const objv[], int *len)
Definition: TclVec.C:537