Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Compound Members   File Members   Related Pages  

TclVec.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr                                                                       
00003  *cr            (C) Copyright 1995-2011 The Board of Trustees of the           
00004  *cr                        University of Illinois                       
00005  *cr                         All Rights Reserved                        
00006  *cr                                                                   
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: TclVec.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.42 $      $Date: 2010/12/16 04:08:43 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   A C-based implementation of some performance-critical Tcl callable
00019  *   routines in VMD.  The C implementation outperforms a raw Tcl version
00020  *   by a factor of three or so.  The performance advantage helps 
00021  *   significantly when doing analysis in VMD.
00022  ***************************************************************************/
00023 
00024 #include <tcl.h>
00025 #include <stdlib.h>
00026 #include <string.h>
00027 #include <math.h>
00028 #include "TclCommands.h"
00029 #include "Matrix4.h"
00030 #include "utilities.h"
00031 
00032 /***************** override some of the vector routines for speed ******/
00033 /* These should be the exact C equivalent to the corresponding Tcl    */
00034 /* vector commands */
00035 
00036 // Function:  vecadd v1 v2 {v3 ...}
00037 //  Returns: the sum of vectors; must all be the same length
00038 //  The increase in speed from Tcl to C++ is 4561 / 255 == 18 fold
00039 static int obj_vecadd(ClientData, Tcl_Interp *interp, int argc, 
00040                        Tcl_Obj * const objv[]) {
00041   if (argc < 3) {
00042     Tcl_WrongNumArgs(interp, 1, objv, (char *)"vec1 vec2 ?vec3? ?vec4? ...");
00043     return TCL_ERROR;
00044   }
00045   int num;
00046   Tcl_Obj **data;
00047   if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) {
00048     return TCL_ERROR;
00049   }
00050   double *sum = new double[num];
00051   int i;
00052   for (i=0; i<num; i++) {
00053     if (Tcl_GetDoubleFromObj(interp, data[i], sum+i) != TCL_OK) {
00054       delete [] sum;
00055       return TCL_ERROR;
00056     }
00057   }
00058   // do the sums on the rest
00059   int num2;
00060   for (int term=2; term < argc; term++) {
00061     if (Tcl_ListObjGetElements(interp, objv[term], &num2, &data) != TCL_OK) {
00062       delete [] sum;
00063       return TCL_ERROR;
00064     }
00065     if (num != num2) {
00066       Tcl_SetResult(interp, (char *) "vecadd: two vectors don't have the same size", TCL_STATIC);
00067       delete [] sum;
00068       return TCL_ERROR;
00069     }
00070     for (i=0; i<num; i++) {
00071       double df;
00072       if (Tcl_GetDoubleFromObj(interp, data[i], &df) != TCL_OK) {
00073         delete [] sum;
00074         return TCL_ERROR;
00075       }
00076       sum[i] += df;
00077     }
00078   }
00079 
00080   
00081   // and return the result
00082   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00083   for (i=0; i<num; i++) {
00084     Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(sum[i]));
00085   }
00086   Tcl_SetObjResult(interp, tcl_result);
00087   delete [] sum;
00088   return TCL_OK;
00089 }
00090 
00091 // Function:  vecsub  v1 v2
00092 //  Returns:   v1 - v2
00093 
00094 
00095 static int obj_vecsub(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[])
00096 {
00097   if (argc != 3) {
00098     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?x? ?y?");
00099     return TCL_ERROR;
00100   }
00101   int num1=0, num2=0;
00102   Tcl_Obj **data1, **data2;
00103   if (Tcl_ListObjGetElements(interp, objv[1], &num1, &data1) != TCL_OK)
00104     return TCL_ERROR;
00105   if (Tcl_ListObjGetElements(interp, objv[2], &num2, &data2) != TCL_OK)
00106     return TCL_ERROR;
00107 
00108   if (num1 != num2) {
00109     Tcl_SetResult(interp, (char *)"vecsub: two vectors don't have the same size", TCL_STATIC);
00110     return TCL_ERROR;
00111   }
00112 
00113   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00114   for (int i=0; i<num1; i++) {
00115     double d1=0, d2=0;
00116     if (Tcl_GetDoubleFromObj(interp, data1[i], &d1) != TCL_OK) {
00117       Tcl_SetResult(interp, (char *)"vecsub: non-numeric in first argument", TCL_STATIC);
00118       return TCL_ERROR; 
00119     }
00120     if (Tcl_GetDoubleFromObj(interp, data2[i], &d2) != TCL_OK) {
00121       Tcl_SetResult(interp, (char *)"vecsub: non-numeric in second argument", TCL_STATIC);
00122       return TCL_ERROR; 
00123     }
00124     Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(d1-d2));
00125   }
00126   Tcl_SetObjResult(interp, tcl_result);
00127   return TCL_OK;
00128 }
00129 
00130 
00131 // Function: vecscale
00132 //  Returns: scalar * vector or vector * scalar
00133 // speedup is 1228/225 = 5.5 fold
00134 static int obj_vecscale(ClientData, Tcl_Interp *interp, int argc, 
00135                        Tcl_Obj * const objv[]) {
00136   if (argc != 3) {
00137     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?c? ?v?");
00138     return TCL_ERROR;
00139   }
00140     
00141   int num1, num2;
00142   Tcl_Obj **data1, **data2;
00143   if (Tcl_ListObjGetElements(interp, objv[1], &num1, &data1) != TCL_OK) {
00144     return TCL_ERROR;
00145   }
00146   if (Tcl_ListObjGetElements(interp, objv[2], &num2, &data2) != TCL_OK) {
00147     return TCL_ERROR;
00148   }
00149   if (num1 == 0 || num2 == 0) {
00150     Tcl_SetResult(interp, (char *) "vecscale: parameters must have data", TCL_STATIC);
00151     return TCL_ERROR;
00152   } else if (num1 != 1 && num2 != 1) {
00153     Tcl_SetResult(interp, (char *) "vecscale: one parameter must be a scalar value", TCL_STATIC);
00154     return TCL_ERROR;
00155   }
00156   
00157   int num;
00158   Tcl_Obj *scalarobj, **vector;
00159   if (num1 == 1) {
00160     scalarobj = data1[0];
00161     vector = data2;
00162     num = num2;
00163   } else {
00164     scalarobj = data2[0];
00165     vector = data1;
00166     num = num1;
00167   }
00168  
00169   double scalar;
00170   if (Tcl_GetDoubleFromObj(interp, scalarobj, &scalar) != TCL_OK)
00171     return TCL_ERROR;
00172 
00173   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00174   for (int i=0; i<num; i++) {
00175     double val;
00176     if (Tcl_GetDoubleFromObj(interp, vector[i], &val) != TCL_OK) {
00177       Tcl_SetResult(interp, (char *) "vecscale: non-numeric in vector", TCL_STATIC);
00178       return TCL_ERROR;
00179     }
00180     val *= scalar;
00181     Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(val));
00182   }
00183   Tcl_SetObjResult(interp, tcl_result);
00184   return TCL_OK;
00185 }
00186 
00188 // returns TCL_OK if good
00189 // If bad, returns TCL_ERROR and sets the Tcl result to the error message
00190 // The name of the function should be passed in 'fctn' so the error message
00191 // can be constructed correctly
00192 int tcl_get_matrix(const char *fctn, Tcl_Interp *interp, 
00193                           Tcl_Obj *s, float *mat)
00194 { 
00195   int num_rows;
00196   Tcl_Obj **data_rows;
00197   if (Tcl_ListObjGetElements(interp, s, &num_rows, &data_rows) != TCL_OK) {
00198     char tmpstring[1024];
00199     sprintf(tmpstring, "%s: badly formed matrix", fctn);
00200     Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
00201     return TCL_ERROR;
00202   }
00203   if (num_rows != 4) {
00204     char tmpstring[1024];
00205     sprintf(tmpstring, "%s: need a 4x4 matrix", fctn);
00206     Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
00207     return TCL_ERROR;
00208   }
00209   int num_row[4];
00210   Tcl_Obj **data_row[4];
00211   if (Tcl_ListObjGetElements(interp, data_rows[0], num_row+0, data_row+0) != TCL_OK ||
00212       num_row[0] != 4 ||
00213       Tcl_ListObjGetElements(interp, data_rows[1], num_row+1, data_row+1) != TCL_OK ||
00214       num_row[1] != 4 ||
00215       Tcl_ListObjGetElements(interp, data_rows[2], num_row+2, data_row+2) != TCL_OK ||
00216       num_row[2] != 4 ||
00217       Tcl_ListObjGetElements(interp, data_rows[3], num_row+3, data_row+3) != TCL_OK ||
00218       num_row[3] != 4) {
00219     Tcl_AppendResult(interp, fctn, ": poorly formed matrix", NULL);
00220     return TCL_ERROR;
00221   }
00222   // now get the numbers
00223   for (int i=0; i<4; i++) {
00224     for (int j=0; j<4; j++) {
00225       double tmp;
00226       if (Tcl_GetDoubleFromObj(interp, data_row[i][j], &tmp) != TCL_OK) {
00227         char tmpstring[1024];
00228         sprintf(tmpstring, "%s: non-numeric in matrix", fctn);
00229         Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
00230         return TCL_ERROR;
00231       } else {
00232         mat[4*j+i] = (float) tmp;  // Matrix4 is transpose of Tcl's matrix
00233       }
00234     }
00235   }
00236   return TCL_OK;
00237 }
00238 
00239 int tcl_get_vector(const char *s, float *val, Tcl_Interp *interp)
00240 {
00241   int num;
00242   const char **pos;
00243   if (Tcl_SplitList(interp, s, &num, &pos) != TCL_OK) {
00244     Tcl_SetResult(interp, (char *) "need three data elements for a vector", 
00245                   TCL_STATIC);
00246     return TCL_ERROR;
00247   }
00248   if (num != 3) {
00249     Tcl_SetResult(interp, (char *) "need three numbers for a vector", TCL_STATIC);
00250     return TCL_ERROR;
00251   }
00252   double a[3];
00253   if (Tcl_GetDouble(interp, pos[0], a+0) != TCL_OK ||
00254       Tcl_GetDouble(interp, pos[1], a+1) != TCL_OK ||
00255       Tcl_GetDouble(interp, pos[2], a+2) != TCL_OK) {
00256     ckfree((char *) pos); // free of tcl data
00257     return TCL_ERROR;
00258   }
00259   val[0] = (float) a[0];
00260   val[1] = (float) a[1];
00261   val[2] = (float) a[2];
00262   ckfree((char *) pos); // free of tcl data
00263   return TCL_OK;
00264 }
00265 
00266 // append the matrix into the Tcl result
00267 void tcl_append_matrix(Tcl_Interp *interp, const float *mat) {
00268   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00269   for (int i=0; i<4; i++) {
00270     Tcl_Obj *m = Tcl_NewListObj(0, NULL);
00271     for (int j=0; j<4; j++) 
00272       Tcl_ListObjAppendElement(interp, m, Tcl_NewDoubleObj(mat[4*j+i]));
00273     Tcl_ListObjAppendElement(interp, tcl_result, m);
00274   }
00275   Tcl_SetObjResult(interp, tcl_result);
00276 }
00277 
00278 // speed up the matrix * vector routines -- DIFFERENT ERROR MESSAGES
00279 // THAN THE TCL VERSION
00280 // speedup is nearly 25 fold
00281 static int obj_vectrans(ClientData, Tcl_Interp *interp, int argc, 
00282                   Tcl_Obj * const objv[])
00283 {
00284   if (argc != 3) {
00285     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?matrix? ?vector?");
00286     return TCL_ERROR;
00287   }
00288 
00289   // get the matrix data
00290   float mat[16];
00291   if (tcl_get_matrix(
00292     Tcl_GetStringFromObj(objv[0],NULL), interp, objv[1], mat) != TCL_OK) {
00293     return TCL_ERROR;
00294   }
00295   
00296   // for the vector
00297   Tcl_Obj **vec;
00298   int vec_size;
00299   if (Tcl_ListObjGetElements(interp, objv[2], &vec_size, &vec) != TCL_OK)
00300     return TCL_ERROR;
00301 
00302   if (vec_size != 3 && vec_size != 4) {
00303     Tcl_SetResult(interp, (char *) "vectrans: vector must be of size 3 or 4",
00304                   TCL_STATIC);
00305     return TCL_ERROR;
00306   }
00307 
00308   float opoint[4];
00309   opoint[3] = 0;
00310   for (int i=0; i<vec_size; i++) {
00311     double tmp;
00312     if (Tcl_GetDoubleFromObj(interp, vec[i], &tmp) != TCL_OK) {
00313       Tcl_SetResult(interp, (char *) "vectrans: non-numeric in vector", TCL_STATIC);
00314       return TCL_ERROR;
00315     }
00316     opoint[i] = (float)tmp;
00317   }
00318   // vector data is in vec_data
00319   float npoint[4];
00320  
00321   npoint[0]=opoint[0]*mat[0]+opoint[1]*mat[4]+opoint[2]*mat[8]+opoint[3]*mat[12]
00322 ;
00323   npoint[1]=opoint[0]*mat[1]+opoint[1]*mat[5]+opoint[2]*mat[9]+opoint[3]*mat[13]
00324 ;
00325   npoint[2]=opoint[0]*mat[2]+opoint[1]*mat[6]+opoint[2]*mat[10]+opoint[3]*mat[14
00326 ];
00327   npoint[3]=opoint[0]*mat[3]+opoint[1]*mat[7]+opoint[2]*mat[11]+opoint[3]*mat[15
00328 ];
00329   // return it
00330 
00331   {
00332   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00333   for (int i=0; i<vec_size; i++) 
00334     Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(npoint[i]));
00335   Tcl_SetObjResult(interp, tcl_result);
00336   }
00337   return TCL_OK;
00338 }
00339 
00340 
00341 // Function: transmult m1 m2 ... mn
00342 //  Returns: the product of the matricies
00343 // speedup is 136347 / 1316 = factor of 104
00344 static int obj_transmult(ClientData, Tcl_Interp *interp, int argc, 
00345                    Tcl_Obj * const objv[]) {
00346   // make there there are at least two values
00347   if (argc < 3) {
00348     Tcl_WrongNumArgs(interp, 1, objv, (char *)"mx my ?m1? ?m2? ...");
00349     return TCL_ERROR;
00350   }
00351   // Get the first matrix
00352   float mult[16];
00353   if (tcl_get_matrix("transmult: ", interp, objv[1], mult) != TCL_OK) {
00354     return TCL_ERROR;
00355   }
00356   int i = 2;
00357   float pre[16];
00358   while (i < argc) {
00359     if (tcl_get_matrix("transmult: ", interp, objv[i], pre) != TCL_OK) {
00360       return TCL_ERROR;
00361     }
00362     // premultiply mult by tmp
00363     float tmp[4];
00364     for (int k=0; k<4; k++) {
00365       tmp[0] = mult[k];
00366       tmp[1] = mult[4+k];
00367       tmp[2] = mult[8+k];
00368       tmp[3] = mult[12+k];
00369       for (int j=0; j<4; j++) {
00370         mult[4*j+k] = pre[4*j]*tmp[0] + pre[4*j+1]*tmp[1] +
00371           pre[4*j+2]*tmp[2] + pre[4*j+3]*tmp[3];
00372       }
00373     }
00374     i++;
00375   }
00376   tcl_append_matrix(interp, mult);
00377   return TCL_OK;
00378 }
00379 
00380 static int obj_transvec(ClientData, Tcl_Interp *interp, int argc, 
00381                    Tcl_Obj * const objv[]) {
00382   if (argc != 2) {
00383     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
00384     return TCL_ERROR;
00385   }
00386   
00387   int num;
00388   Tcl_Obj **data;
00389   if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 
00390     return TCL_ERROR;
00391   if (num != 3) {
00392     Tcl_AppendResult(interp, "transvec: vector must have three elements",NULL);
00393     return TCL_ERROR;
00394   }
00395   double x,y,z;
00396   if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK ||
00397       Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK ||
00398       Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) {
00399     Tcl_SetResult(interp, (char *)"transvec: non-numeric in vector", TCL_STATIC);
00400     return TCL_ERROR;
00401   }
00402   Matrix4 mat;
00403   mat.transvec((float) x,(float) y,(float) z);
00404   tcl_append_matrix(interp, mat.mat);
00405   return TCL_OK;
00406 }
00407 
00408 static int obj_transvecinv(ClientData, Tcl_Interp *interp, int argc, 
00409                    Tcl_Obj * const objv[]) {
00410   if (argc != 2) {
00411     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
00412     return TCL_ERROR;
00413   }
00414   
00415   int num;
00416   Tcl_Obj **data;
00417   if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 
00418     return TCL_ERROR;
00419   if (num != 3) {
00420     Tcl_AppendResult(interp, "transvecinv: vector must have three elements",NULL);
00421     return TCL_ERROR;
00422   }
00423   double x,y,z;
00424   if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK ||
00425       Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK ||
00426       Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) {
00427     Tcl_SetResult(interp, (char *)"transvecinv: non-numeric in vector", TCL_STATIC);
00428     return TCL_ERROR;
00429   }
00430   Matrix4 mat;
00431   mat.transvecinv((float) x,(float) y,(float) z);
00432   tcl_append_matrix(interp, mat.mat);
00433   return TCL_OK;
00434 }
00435 
00436 // Returns the transformation matrix needed to rotate by a certain
00437 // angle around a given axis.
00438 // Tcl syntax:
00439 // transabout v amount [deg|rad|pi]
00440 // The increase in speed from Tcl to C++ is 15 fold
00441 static int obj_transabout(ClientData, Tcl_Interp *interp, int argc, 
00442                    Tcl_Obj * const objv[]) {
00443   if (argc != 3 && argc != 4) {
00444     Tcl_WrongNumArgs(interp, 1, objv, (char *)"axis amount [deg|rad|pi]");
00445     return TCL_ERROR;
00446   }
00447   
00448   int num;
00449   Tcl_Obj **data;
00450   // get the axis
00451   if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 
00452     return TCL_ERROR;
00453   if (num != 3) {
00454     Tcl_AppendResult(interp, "transabout: vector must have three elements",NULL);
00455     return TCL_ERROR;
00456   }
00457   double x,y,z;
00458   if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK ||
00459       Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK ||
00460       Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) {
00461     Tcl_SetResult(interp, (char *)"transabout: non-numeric in vector", TCL_STATIC);
00462     return TCL_ERROR;
00463   }
00464 
00465   // get the amount
00466   double amount;
00467   if (Tcl_GetDoubleFromObj(interp, objv[2], &amount) != TCL_OK) {
00468     Tcl_SetResult(interp, (char *)"transabout: non-numeric angle", TCL_STATIC);
00469     return TCL_ERROR;
00470   }
00471 
00472   // get units
00473   if (argc == 4) {
00474     if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "deg")) {
00475       amount = DEGTORAD(amount);
00476     } else if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "rad")) {
00477       // amount = amount; 
00478     } else if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "pi")) {
00479       amount = amount*VMD_PI;
00480     } else {
00481       Tcl_AppendResult(interp, "transabout: unit must be deg|rad|pi",NULL);
00482       return TCL_ERROR;
00483     }
00484   } else {
00485     // If no unit was specified assume that we have degrees
00486     amount = DEGTORAD(amount);
00487   }
00488 
00489   float axis[3];
00490   axis[0] = (float) x;
00491   axis[1] = (float) y;
00492   axis[2] = (float) z;
00493 
00494   Matrix4 mat;
00495   mat.rotate_axis(axis, (float) amount);
00496   tcl_append_matrix(interp, mat.mat);
00497   return TCL_OK;
00498 }
00499 
00500 static int obj_veclength(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) {
00501 
00502   if (argc != 2) {
00503     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
00504     return TCL_ERROR;
00505   }
00506 
00507   int num;
00508   Tcl_Obj **data;
00509   if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 
00510     return TCL_ERROR;
00511 
00512   double length = 0.;
00513   for (int i=0; i<num; i++) {
00514     double tmp;
00515     if (Tcl_GetDoubleFromObj(interp, data[i], &tmp) != TCL_OK) {
00516       Tcl_SetResult(interp, (char *) "veclength: non-numeric in vector", TCL_STATIC);
00517       return TCL_ERROR;
00518     } else {
00519       length += tmp*tmp;
00520     }
00521   }
00522 
00523   length = sqrt(length);
00524   Tcl_Obj *tcl_result = Tcl_GetObjResult(interp);
00525   Tcl_SetDoubleObj(tcl_result, length);
00526   return TCL_OK; 
00527 }
00528 
00529 
00530 static double* obj_getdoublearray(Tcl_Interp *interp, Tcl_Obj *const objv[], int *len) {
00531   int num;
00532 
00533   Tcl_Obj **data;
00534   if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK)
00535     return NULL;
00536  
00537   double *list = (double*) malloc(num*sizeof(double));
00538   if (list == NULL)
00539     return NULL;
00540 
00541   for (int i=0; i<num; i++) {
00542     double tmp;
00543     if (Tcl_GetDoubleFromObj(interp, data[i], &tmp) != TCL_OK) {
00544       Tcl_SetResult(interp, (char *) "veclength: non-numeric in vector", TCL_STATIC);
00545       free(list);
00546       return NULL;
00547     }
00548     list[i] = tmp;
00549   }
00550 
00551   *len = num;
00552 
00553   return list;
00554 }
00555 
00556 
00557 static int obj_vecsum(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) {
00558   if (argc != 2) {
00559     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
00560     return TCL_ERROR;
00561   }
00562 
00563   int num;
00564   double *list = obj_getdoublearray(interp, objv, &num);
00565   if (list == NULL) 
00566     return TCL_ERROR;
00567 
00568   double sum = 0.;
00569   for (int i=0; i<num; i++) {
00570     sum += list[i];
00571   }
00572   free(list);
00573 
00574   Tcl_Obj *tcl_result = Tcl_GetObjResult(interp);
00575   Tcl_SetDoubleObj(tcl_result, sum);
00576   return TCL_OK;
00577 }
00578 
00579 
00580 static int obj_vecmean(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) {
00581   if (argc != 2) {
00582     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
00583     return TCL_ERROR;
00584   }
00585 
00586   int num;
00587   double *list = obj_getdoublearray(interp, objv, &num);
00588   if (list == NULL) 
00589     return TCL_ERROR;
00590 
00591   double sum = 0.;
00592   for (int i=0; i<num; i++) {
00593     sum += list[i];
00594   }
00595   sum /= (double) num;
00596   free(list);
00597 
00598   Tcl_Obj *tcl_result = Tcl_GetObjResult(interp);
00599   Tcl_SetDoubleObj(tcl_result, sum);
00600   return TCL_OK;
00601 }
00602 
00603 
00604 static int obj_vecstddev(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) {
00605   if (argc != 2) {
00606     Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?");
00607     return TCL_ERROR;
00608   }
00609 
00610   int i, num;
00611   double* list = obj_getdoublearray(interp, objv, &num);
00612   if (list == NULL) 
00613     return TCL_ERROR;
00614 
00615   double mean = 0.;
00616   for (i=0; i<num; i++) {
00617     mean += list[i];
00618   }
00619   mean /= (double) num;
00620 
00621   double stddev = 0.;
00622   for (i=0; i<num; i++) {
00623     double tmp = list[i] - mean;
00624     stddev += tmp * tmp; 
00625   }
00626   stddev /= (double) num;
00627   stddev = sqrt(stddev);
00628   free(list);
00629 
00630   Tcl_Obj *tcl_result = Tcl_GetObjResult(interp);
00631   Tcl_SetDoubleObj(tcl_result, stddev);
00632   return TCL_OK;
00633 }
00634 
00635 
00636 int Vec_Init(Tcl_Interp *interp) {
00637   Tcl_CreateObjCommand(interp, "vecadd", obj_vecadd,
00638                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00639   Tcl_CreateObjCommand(interp, "vecsub", obj_vecsub,
00640                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00641   Tcl_CreateObjCommand(interp, "vecscale", obj_vecscale,
00642                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00643   Tcl_CreateObjCommand(interp, "transmult", obj_transmult,
00644                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00645   Tcl_CreateObjCommand(interp, "vectrans", obj_vectrans,
00646                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00647   Tcl_CreateObjCommand(interp, "veclength", obj_veclength,
00648                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00649   Tcl_CreateObjCommand(interp, "vecmean", obj_vecmean,
00650                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00651   Tcl_CreateObjCommand(interp, "vecstddev", obj_vecstddev,
00652                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00653   Tcl_CreateObjCommand(interp, "vecsum", obj_vecsum,
00654                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00655   Tcl_CreateObjCommand(interp, "transvec", obj_transvec,
00656                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00657   Tcl_CreateObjCommand(interp, "transvecinv", obj_transvecinv,
00658                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00659   Tcl_CreateObjCommand(interp, "transabout", obj_transabout,
00660                     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00661   return TCL_OK;
00662 }
00663  

Generated on Wed May 22 01:49:08 2013 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002