00001 /*************************************************************************** 00002 *cr 00003 *cr (C) Copyright 1995-2019 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.46 $ $Date: 2020年11月04日 20:57:00 $ 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 00240 int tcl_get_vector(const char *s, float *val, Tcl_Interp *interp) { 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 00267 int tcl_get_vecarray(const char *s, int &num, float *&val, Tcl_Interp *interp) { 00268 const char **pos; 00269 int cnt; 00270 if (Tcl_SplitList(interp, s, &cnt, &pos) != TCL_OK) { 00271 Tcl_SetResult(interp, (char *) "need numeric values", TCL_STATIC); 00272 return TCL_ERROR; 00273 } 00274 00275 num = 3*cnt; 00276 val = new float[num]; 00277 for (int i=0; i<cnt; i++) { 00278 if (tcl_get_vector(pos[i], val+i*3, interp) != TCL_OK) { 00279 num = 0; 00280 delete [] val; 00281 val = NULL; 00282 ckfree((char *) pos); // free of tcl data 00283 return TCL_ERROR; 00284 } 00285 } 00286 00287 ckfree((char *) pos); // free of tcl data 00288 return TCL_OK; 00289 } 00290 00291 00292 int tcl_get_array(const char *s, int &num, float *&val, Tcl_Interp *interp) { 00293 const char **pos; 00294 if (Tcl_SplitList(interp, s, &num, &pos) != TCL_OK) { 00295 Tcl_SetResult(interp, (char *) "need numeric values", TCL_STATIC); 00296 return TCL_ERROR; 00297 } 00298 00299 val = new float[num]; 00300 for (int i=0; i<num; i++) { 00301 double a; 00302 if (Tcl_GetDouble(interp, pos[i], &a) != TCL_OK) { 00303 num = 0; 00304 delete [] val; 00305 val = NULL; 00306 ckfree((char *) pos); // free of tcl data 00307 return TCL_ERROR; 00308 } 00309 val[i] = (float) a; 00310 } 00311 00312 ckfree((char *) pos); // free of tcl data 00313 return TCL_OK; 00314 } 00315 00316 00317 int tcl_get_intarray(const char *s, int &num, int *&val, Tcl_Interp *interp) { 00318 const char **pos; 00319 if (Tcl_SplitList(interp, s, &num, &pos) != TCL_OK) { 00320 Tcl_SetResult(interp, (char *) "need numeric values", TCL_STATIC); 00321 return TCL_ERROR; 00322 } 00323 00324 val = new int[num]; 00325 for (int i=0; i<num; i++) { 00326 int a; 00327 if (Tcl_GetInt(interp, pos[i], &a) != TCL_OK) { 00328 num = 0; 00329 delete [] val; 00330 val = NULL; 00331 ckfree((char *) pos); // free of tcl data 00332 return TCL_ERROR; 00333 } 00334 val[i] = (int) a; 00335 } 00336 00337 ckfree((char *) pos); // free of tcl data 00338 return TCL_OK; 00339 } 00340 00341 00342 // append the matrix into the Tcl result 00343 void tcl_append_matrix(Tcl_Interp *interp, const float *mat) { 00344 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL); 00345 for (int i=0; i<4; i++) { 00346 Tcl_Obj *m = Tcl_NewListObj(0, NULL); 00347 for (int j=0; j<4; j++) 00348 Tcl_ListObjAppendElement(interp, m, Tcl_NewDoubleObj(mat[4*j+i])); 00349 Tcl_ListObjAppendElement(interp, tcl_result, m); 00350 } 00351 Tcl_SetObjResult(interp, tcl_result); 00352 } 00353 00354 00355 // speed up the matrix * vector routines -- DIFFERENT ERROR MESSAGES 00356 // THAN THE TCL VERSION 00357 // speedup is nearly 25 fold 00358 static int obj_vectrans(ClientData, Tcl_Interp *interp, int argc, 00359 Tcl_Obj * const objv[]) { 00360 if (argc != 3) { 00361 Tcl_WrongNumArgs(interp, 1, objv, (char *)"?matrix? ?vector?"); 00362 return TCL_ERROR; 00363 } 00364 00365 // get the matrix data 00366 float mat[16]; 00367 if (tcl_get_matrix( 00368 Tcl_GetStringFromObj(objv[0],NULL), interp, objv[1], mat) != TCL_OK) { 00369 return TCL_ERROR; 00370 } 00371 00372 // for the vector 00373 Tcl_Obj **vec; 00374 int vec_size; 00375 if (Tcl_ListObjGetElements(interp, objv[2], &vec_size, &vec) != TCL_OK) 00376 return TCL_ERROR; 00377 00378 if (vec_size != 3 && vec_size != 4) { 00379 Tcl_SetResult(interp, (char *) "vectrans: vector must be of size 3 or 4", 00380 TCL_STATIC); 00381 return TCL_ERROR; 00382 } 00383 00384 float opoint[4]; 00385 opoint[3] = 0; 00386 for (int i=0; i<vec_size; i++) { 00387 double tmp; 00388 if (Tcl_GetDoubleFromObj(interp, vec[i], &tmp) != TCL_OK) { 00389 Tcl_SetResult(interp, (char *) "vectrans: non-numeric in vector", TCL_STATIC); 00390 return TCL_ERROR; 00391 } 00392 opoint[i] = (float)tmp; 00393 } 00394 // vector data is in vec_data 00395 float npoint[4]; 00396 00397 npoint[0]=opoint[0]*mat[0]+opoint[1]*mat[4]+opoint[2]*mat[8]+opoint[3]*mat[12] 00398 ; 00399 npoint[1]=opoint[0]*mat[1]+opoint[1]*mat[5]+opoint[2]*mat[9]+opoint[3]*mat[13] 00400 ; 00401 npoint[2]=opoint[0]*mat[2]+opoint[1]*mat[6]+opoint[2]*mat[10]+opoint[3]*mat[14 00402 ]; 00403 npoint[3]=opoint[0]*mat[3]+opoint[1]*mat[7]+opoint[2]*mat[11]+opoint[3]*mat[15 00404 ]; 00405 // return it 00406 00407 { 00408 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL); 00409 for (int i=0; i<vec_size; i++) 00410 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(npoint[i])); 00411 Tcl_SetObjResult(interp, tcl_result); 00412 } 00413 return TCL_OK; 00414 } 00415 00416 00417 // Function: transmult m1 m2 ... mn 00418 // Returns: the product of the matricies 00419 // speedup is 136347 / 1316 = factor of 104 00420 static int obj_transmult(ClientData, Tcl_Interp *interp, int argc, 00421 Tcl_Obj * const objv[]) { 00422 // make there there are at least two values 00423 if (argc < 3) { 00424 Tcl_WrongNumArgs(interp, 1, objv, (char *)"mx my ?m1? ?m2? ..."); 00425 return TCL_ERROR; 00426 } 00427 // Get the first matrix 00428 float mult[16]; 00429 if (tcl_get_matrix("transmult: ", interp, objv[1], mult) != TCL_OK) { 00430 return TCL_ERROR; 00431 } 00432 int i = 2; 00433 float pre[16]; 00434 while (i < argc) { 00435 if (tcl_get_matrix("transmult: ", interp, objv[i], pre) != TCL_OK) { 00436 return TCL_ERROR; 00437 } 00438 // premultiply mult by tmp 00439 float tmp[4]; 00440 for (int k=0; k<4; k++) { 00441 tmp[0] = mult[k]; 00442 tmp[1] = mult[4+k]; 00443 tmp[2] = mult[8+k]; 00444 tmp[3] = mult[12+k]; 00445 for (int j=0; j<4; j++) { 00446 mult[4*j+k] = pre[4*j]*tmp[0] + pre[4*j+1]*tmp[1] + 00447 pre[4*j+2]*tmp[2] + pre[4*j+3]*tmp[3]; 00448 } 00449 } 00450 i++; 00451 } 00452 tcl_append_matrix(interp, mult); 00453 return TCL_OK; 00454 } 00455 00456 static int obj_transvec(ClientData, Tcl_Interp *interp, int argc, 00457 Tcl_Obj * const objv[]) { 00458 if (argc != 2) { 00459 Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?"); 00460 return TCL_ERROR; 00461 } 00462 00463 int num; 00464 Tcl_Obj **data; 00465 if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 00466 return TCL_ERROR; 00467 if (num != 3) { 00468 Tcl_AppendResult(interp, "transvec: vector must have three elements",NULL); 00469 return TCL_ERROR; 00470 } 00471 double x,y,z; 00472 if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK || 00473 Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK || 00474 Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) { 00475 Tcl_SetResult(interp, (char *)"transvec: non-numeric in vector", TCL_STATIC); 00476 return TCL_ERROR; 00477 } 00478 Matrix4 mat; 00479 mat.transvec((float) x,(float) y,(float) z); 00480 tcl_append_matrix(interp, mat.mat); 00481 return TCL_OK; 00482 } 00483 00484 static int obj_transvecinv(ClientData, Tcl_Interp *interp, int argc, 00485 Tcl_Obj * const objv[]) { 00486 if (argc != 2) { 00487 Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?"); 00488 return TCL_ERROR; 00489 } 00490 00491 int num; 00492 Tcl_Obj **data; 00493 if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 00494 return TCL_ERROR; 00495 if (num != 3) { 00496 Tcl_AppendResult(interp, "transvecinv: vector must have three elements",NULL); 00497 return TCL_ERROR; 00498 } 00499 double x,y,z; 00500 if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK || 00501 Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK || 00502 Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) { 00503 Tcl_SetResult(interp, (char *)"transvecinv: non-numeric in vector", TCL_STATIC); 00504 return TCL_ERROR; 00505 } 00506 Matrix4 mat; 00507 mat.transvecinv((float) x,(float) y,(float) z); 00508 tcl_append_matrix(interp, mat.mat); 00509 return TCL_OK; 00510 } 00511 00512 // Returns the transformation matrix needed to rotate by a certain 00513 // angle around a given axis. 00514 // Tcl syntax: 00515 // transabout v amount [deg|rad|pi] 00516 // The increase in speed from Tcl to C++ is 15 fold 00517 static int obj_transabout(ClientData, Tcl_Interp *interp, int argc, 00518 Tcl_Obj * const objv[]) { 00519 if (argc != 3 && argc != 4) { 00520 Tcl_WrongNumArgs(interp, 1, objv, (char *)"axis amount [deg|rad|pi]"); 00521 return TCL_ERROR; 00522 } 00523 00524 int num; 00525 Tcl_Obj **data; 00526 // get the axis 00527 if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 00528 return TCL_ERROR; 00529 if (num != 3) { 00530 Tcl_AppendResult(interp, "transabout: vector must have three elements",NULL); 00531 return TCL_ERROR; 00532 } 00533 double x,y,z; 00534 if (Tcl_GetDoubleFromObj(interp, data[0], &x) != TCL_OK || 00535 Tcl_GetDoubleFromObj(interp, data[1], &y) != TCL_OK || 00536 Tcl_GetDoubleFromObj(interp, data[2], &z) != TCL_OK) { 00537 Tcl_SetResult(interp, (char *)"transabout: non-numeric in vector", TCL_STATIC); 00538 return TCL_ERROR; 00539 } 00540 00541 // get the amount 00542 double amount; 00543 if (Tcl_GetDoubleFromObj(interp, objv[2], &amount) != TCL_OK) { 00544 Tcl_SetResult(interp, (char *)"transabout: non-numeric angle", TCL_STATIC); 00545 return TCL_ERROR; 00546 } 00547 00548 // get units 00549 if (argc == 4) { 00550 if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "deg")) { 00551 amount = DEGTORAD(amount); 00552 } else if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "rad")) { 00553 // amount = amount; 00554 } else if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "pi")) { 00555 amount = amount*VMD_PI; 00556 } else { 00557 Tcl_AppendResult(interp, "transabout: unit must be deg|rad|pi",NULL); 00558 return TCL_ERROR; 00559 } 00560 } else { 00561 // If no unit was specified assume that we have degrees 00562 amount = DEGTORAD(amount); 00563 } 00564 00565 float axis[3]; 00566 axis[0] = (float) x; 00567 axis[1] = (float) y; 00568 axis[2] = (float) z; 00569 00570 Matrix4 mat; 00571 mat.rotate_axis(axis, (float) amount); 00572 tcl_append_matrix(interp, mat.mat); 00573 return TCL_OK; 00574 } 00575 00576 static int obj_veclength(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) { 00577 00578 if (argc != 2) { 00579 Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?"); 00580 return TCL_ERROR; 00581 } 00582 00583 int num; 00584 Tcl_Obj **data; 00585 if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 00586 return TCL_ERROR; 00587 00588 double length = 0.; 00589 for (int i=0; i<num; i++) { 00590 double tmp; 00591 if (Tcl_GetDoubleFromObj(interp, data[i], &tmp) != TCL_OK) { 00592 Tcl_SetResult(interp, (char *) "veclength: non-numeric in vector", TCL_STATIC); 00593 return TCL_ERROR; 00594 } else { 00595 length += tmp*tmp; 00596 } 00597 } 00598 00599 length = sqrt(length); 00600 Tcl_Obj *tcl_result = Tcl_GetObjResult(interp); 00601 Tcl_SetDoubleObj(tcl_result, length); 00602 return TCL_OK; 00603 } 00604 00605 00606 static double* obj_getdoublearray(Tcl_Interp *interp, Tcl_Obj *const objv[], int *len) { 00607 int num; 00608 00609 Tcl_Obj **data; 00610 if (Tcl_ListObjGetElements(interp, objv[1], &num, &data) != TCL_OK) 00611 return NULL; 00612 00613 double *list = (double*) malloc(num*sizeof(double)); 00614 if (list == NULL) 00615 return NULL; 00616 00617 for (int i=0; i<num; i++) { 00618 double tmp; 00619 if (Tcl_GetDoubleFromObj(interp, data[i], &tmp) != TCL_OK) { 00620 Tcl_SetResult(interp, (char *) "veclength: non-numeric in vector", TCL_STATIC); 00621 free(list); 00622 return NULL; 00623 } 00624 list[i] = tmp; 00625 } 00626 00627 *len = num; 00628 00629 return list; 00630 } 00631 00632 00633 static int obj_vecsum(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) { 00634 if (argc != 2) { 00635 Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?"); 00636 return TCL_ERROR; 00637 } 00638 00639 int num; 00640 double *list = obj_getdoublearray(interp, objv, &num); 00641 if (list == NULL) 00642 return TCL_ERROR; 00643 00644 double sum = 0.; 00645 for (int i=0; i<num; i++) { 00646 sum += list[i]; 00647 } 00648 free(list); 00649 00650 Tcl_Obj *tcl_result = Tcl_GetObjResult(interp); 00651 Tcl_SetDoubleObj(tcl_result, sum); 00652 return TCL_OK; 00653 } 00654 00655 00656 static int obj_vecmean(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) { 00657 if (argc != 2) { 00658 Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?"); 00659 return TCL_ERROR; 00660 } 00661 00662 int num; 00663 double *list = obj_getdoublearray(interp, objv, &num); 00664 if (list == NULL) 00665 return TCL_ERROR; 00666 00667 double sum = 0.; 00668 for (int i=0; i<num; i++) { 00669 sum += list[i]; 00670 } 00671 sum /= (double) num; 00672 free(list); 00673 00674 Tcl_Obj *tcl_result = Tcl_GetObjResult(interp); 00675 Tcl_SetDoubleObj(tcl_result, sum); 00676 return TCL_OK; 00677 } 00678 00679 00680 static int obj_vecstddev(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const objv[]) { 00681 if (argc != 2) { 00682 Tcl_WrongNumArgs(interp, 1, objv, (char *)"?vector?"); 00683 return TCL_ERROR; 00684 } 00685 00686 int i, num; 00687 double* list = obj_getdoublearray(interp, objv, &num); 00688 if (list == NULL) 00689 return TCL_ERROR; 00690 00691 double mean = 0.; 00692 for (i=0; i<num; i++) { 00693 mean += list[i]; 00694 } 00695 mean /= (double) num; 00696 00697 double stddev = 0.; 00698 for (i=0; i<num; i++) { 00699 double tmp = list[i] - mean; 00700 stddev += tmp * tmp; 00701 } 00702 stddev /= (double) num; 00703 stddev = sqrt(stddev); 00704 free(list); 00705 00706 Tcl_Obj *tcl_result = Tcl_GetObjResult(interp); 00707 Tcl_SetDoubleObj(tcl_result, stddev); 00708 return TCL_OK; 00709 } 00710 00711 00712 int Vec_Init(Tcl_Interp *interp) { 00713 Tcl_CreateObjCommand(interp, "vecadd", obj_vecadd, 00714 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00715 Tcl_CreateObjCommand(interp, "vecsub", obj_vecsub, 00716 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00717 Tcl_CreateObjCommand(interp, "vecscale", obj_vecscale, 00718 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00719 Tcl_CreateObjCommand(interp, "transmult", obj_transmult, 00720 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00721 Tcl_CreateObjCommand(interp, "vectrans", obj_vectrans, 00722 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00723 Tcl_CreateObjCommand(interp, "veclength", obj_veclength, 00724 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00725 Tcl_CreateObjCommand(interp, "vecmean", obj_vecmean, 00726 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00727 Tcl_CreateObjCommand(interp, "vecstddev", obj_vecstddev, 00728 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00729 Tcl_CreateObjCommand(interp, "vecsum", obj_vecsum, 00730 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00731 Tcl_CreateObjCommand(interp, "transvec", obj_transvec, 00732 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00733 Tcl_CreateObjCommand(interp, "transvecinv", obj_transvecinv, 00734 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00735 Tcl_CreateObjCommand(interp, "transabout", obj_transabout, 00736 (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL); 00737 return TCL_OK; 00738 } 00739