TclVec.C

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

Generated on Tue Sep 19 01:17:14 2017 for NAMD by  doxygen 1.4.7