[GRASS-SVN] r36020 -
grass/branches/releasebranch_6_4/vector/v.vol.rst
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Feb 21 12:20:49 EST 2009
Author: neteler
Date: 2009-02-21 12:20:49 -0500 (Sat, 21 Feb 2009)
New Revision: 36020
Modified:
grass/branches/releasebranch_6_4/vector/v.vol.rst/user.h
grass/branches/releasebranch_6_4/vector/v.vol.rst/user3.c
grass/branches/releasebranch_6_4/vector/v.vol.rst/user4.c
Log:
Antonio Galea: reduced computational time to 30-40% by skipping unnecessary computations (merge from trunk, r36018)
Modified: grass/branches/releasebranch_6_4/vector/v.vol.rst/user.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.vol.rst/user.h 2009-02-21 17:19:47 UTC (rev 36019)
+++ grass/branches/releasebranch_6_4/vector/v.vol.rst/user.h 2009-02-21 17:20:49 UTC (rev 36020)
@@ -14,6 +14,7 @@
double crsd();
double crsdr2();
double crsdd();
+void crs_full();
int COGRR1();
int POINT();
int LINEQS();
Modified: grass/branches/releasebranch_6_4/vector/v.vol.rst/user3.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.vol.rst/user3.c 2009-02-21 17:19:47 UTC (rev 36019)
+++ grass/branches/releasebranch_6_4/vector/v.vol.rst/user3.c 2009-02-21 17:20:49 UTC (rev 36020)
@@ -235,6 +235,7 @@
double wm, dx, dy, dz, dxx, dyy, dxy, dxz, dyz, dzz, h, bmgd1,
bmgd2, etar, zcon, r, ww, wz, r2, hcell, zzcell2,
etarcell, rcell, wwcell, zzcell;
+ double x_crs,x_crsd,x_crsdd,x_crsdr2;
int n1, k1, k2, k, i1, l, l1, n4, n5, m, i;
int NGST, LSIZE, ngstc, nszc, ngstr, nszr, ngstl, nszl;
int POINT();
@@ -244,6 +245,8 @@
int bmask = 1;
static FCELL *cell = NULL;
+ int cond1 = (gradient != NULL) || (aspect1 != NULL) || (aspect2 != NULL);
+ int cond2 = (ncurv != NULL) || (gcurv != NULL) || (mcurv != NULL);
#define CEULER .57721566
/*
@@ -450,25 +453,31 @@
r = sqrt(r2);
etar = (fi * r) / 2.;
- h = h + b[m] * crs(etar);
- /* partial derivatives */
- bmgd1 = b[m] * crsd(etar, fi);
+ crs_full(
+ etar,fi,
+ &x_crs,
+ cond1?&x_crsd:NULL,
+ cond2?&x_crsdr2:NULL,
+ cond2?&x_crsdd:NULL
+ );
+ h = h + b[m] * x_crs;
+ if(cond1)
+ {
+ bmgd1 = b[m] * x_crsd;
dx = dx + bmgd1 * xx;
dy = dy + bmgd1 * w[m];
- /* dz = dz + bmgd1 * sqrt(wz2[m]); */
dz = dz + bmgd1 * wz1[m];
- bmgd2 = b[m] * crsdd(etar, fi);
- bmgd1 = b[m] * crsdr2(etar, fi);
- dxx = dxx + bmgd2 * xx2 + bmgd1 * xx2;
+ }
+ if(cond2)
+ {
+ bmgd2 = b[m] * x_crsdd;
+ bmgd1 = b[m] * x_crsdr2;
dyy = dyy + bmgd2 * w2[m] + bmgd1 * w2[m];
dzz = dzz + bmgd2 * wz2[m] + bmgd1 * wz2[m];
dxy = dxy + bmgd2 * xx * w[m] + bmgd1 * xx * w[m];
- dxz =
- dxz + bmgd2 * xx * wz1[m] +
- bmgd1 * xx * wz1[m];
- dyz =
- dyz + bmgd2 * w[m] * wz1[m] +
- bmgd1 * w[m] * wz1[m];
+ dxz = dxz + bmgd2 * xx * wz1[m] + bmgd1 * xx * wz1[m];
+ dyz = dyz + bmgd2 * w[m] * wz1[m] + bmgd1 * w[m] * wz1[m];
+ }
}
ww = h + wmin;
if ((cellinp != NULL) && (cellout != NULL) &&
Modified: grass/branches/releasebranch_6_4/vector/v.vol.rst/user4.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.vol.rst/user4.c 2009-02-21 17:19:47 UTC (rev 36019)
+++ grass/branches/releasebranch_6_4/vector/v.vol.rst/user4.c 2009-02-21 17:20:49 UTC (rev 36020)
@@ -145,7 +145,7 @@
static double a[5] = { 0.254829592, -0.284496736, 1.421413741,
-1.453152027, 1.061405429
};
- double p = 0.3275911;
+ static double p = 0.3275911;
double erf, t;
@@ -167,21 +167,21 @@
*/
{
- static double a[10];
- double tp = 1.1283791671 / 2.;
+ static double a[10] = {
+ 1.,
+ -1. / 3.,
+ 1. / 10.,
+ -1. / 42.,
+ 1. / (24. * 9.),
+ -1. / (120. * 11.),
+ 1. / (720. * 13.),
+ -1. / (5040. * 15.),
+ 1. / (40320. * 17.),
+ -1. / (362880. * 19.),
+ };
+ static double tp = 1.1283791671 / 2.;
double xx, res;
- a[0] = 1.;
- a[1] = -1. / 3.;
- a[2] = 1. / 10.;
- a[3] = -1. / 42.;
- a[4] = 1. / (24. * 9.);
- a[5] = -1. / (120. * 11.);
- a[6] = 1. / (720. * 13.);
- a[7] = -1. / (5040. * 15.);
- a[8] = 1. / (40320. * 17.);
- a[9] = -1. / (362880. * 19.);
-
if (x < XA) {
xx = x * x;
res =
@@ -204,37 +204,24 @@
}
+void crs_full(double x, double fi, double *crs, double *crsd, double *crsdr2, double *crsdd){
+ static double a[10] = { 1., -1. / 3., 1. / 10., -1. / 42., 1. / (24. * 9.), -1. / (120. * 11.),
+ 1. / (720. * 13.), -1. / (5040. * 15.), 1. / (40320. * 17.), -1. / (362880. * 19.) };
+ static double b[10] = { 0, -2. / 3., 4. / 10., -6. / 42., 8. / (24. * 9.), -10. / (120. * 11.),
+ 12. / (720. * 13.), -14. / (5040. * 15.), 16. / (40320. * 17.), -18. / (362880. * 19.) };
+ static double c[10] = { 0, 0, 8. / 10., -24. / 42., 48. / (24. * 9.), -80. / (120. * 11.),
+ 120. / (720. * 13.), -168. / (5040. * 15.), 16. * 14. / (40320. * 17.), -18. * 16. / (362880. * 19.) };
+ static double tp = 1.1283791671 / 2.;
+ double xx, r, r2, fi2, fi4, fi8, tmp1, tmp2;
-
-double crsd(double x, double fi)
-/*
- function for first derivatives
- */
-{
-
- static double a[10];
- double tp = 1.1283791671 / 2.;
- double xx, r, r2, fi2, fi4, res;
-
- a[1] = -2. / 3.;
- a[2] = 4. / 10.;
- a[3] = -6. / 42.;
- a[4] = 8. / (24. * 9.);
- a[5] = -10. / (120. * 11.);
- a[6] = 12. / (720. * 13.);
- a[7] = -14. / (5040. * 15.);
- a[8] = 16. / (40320. * 17.);
- a[9] = -18. / (362880. * 19.);
-
- r = 2. * x / fi;
- r2 = r * r;
fi2 = fi / 2.;
fi4 = fi2 * fi2;
+ fi8 = fi4 * fi4;
if (x < XA) {
xx = x * x;
- res =
- tp * fi4 * (a[1] +
+ *crs =
+ tp * xx * (a[1] +
xx * (a[2] +
xx * (a[3] +
xx * (a[4] +
@@ -243,111 +230,52 @@
xx * (a[7] +
xx * (a[8] +
xx *
- a
- [9]))))))));
- }
- else {
- res = (2. * tp * exp(-x * x) - (erf(x) / x)) / (2. * r2);
- }
- return (res);
-}
-
-double crsdr2(double x, double fi)
-/*
- function for first derivatives
- */
-{
-
- static double a[10];
- double tp = 1.1283791671 / 2.;
- double xx, r, r2, fi2, fi4, res;
-
- a[1] = -2. / 3.;
- a[2] = 4. / 10.;
- a[3] = -6. / 42.;
- a[4] = 8. / (24. * 9.);
- a[5] = -10. / (120. * 11.);
- a[6] = 12. / (720. * 13.);
- a[7] = -14. / (5040. * 15.);
- a[8] = 16. / (40320. * 17.);
- a[9] = -18. / (362880. * 19.);
-
- r = 2. * x / fi;
- r2 = r * r;
- fi2 = fi / 2.;
- fi4 = fi2 * fi2;
-
- if (x < XA) {
- xx = x * x;
- res =
- tp * fi4 * (a[1] +
- xx * (a[2] +
- xx * (a[3] +
- xx * (a[4] +
- xx * (a[5] +
- xx * (a[6] +
- xx * (a[7] +
- xx * (a[8] +
+ a[9]))))))));
+ if(crsd!=NULL){
+ *crsd =
+ tp * fi4 * (b[1] +
+ xx * (b[2] +
+ xx * (b[3] +
+ xx * (b[4] +
+ xx * (b[5] +
+ xx * (b[6] +
+ xx * (b[7] +
+ xx * (b[8] +
xx *
- a
- [9]))))))));
+ b[9]))))))));
}
- else {
- res = (2. * tp * exp(-x * x) - (erf(x) / x)) / (2. * r2 * r2);
+ if(crsdr2!=NULL){
+ *crsdr2 = *crsd;
}
- return (res);
-}
-
-double crsdd(double x, double fi)
-/*
- function for second derivatives
- */
-{
-
- static double a[10];
- double tp = 1.1283791671 / 2.;
- double xx, r, r2, fi8, res;
-
- a[2] = 8. / 10.;
- a[3] = -24. / 42.;
- a[4] = 48. / (24. * 9.);
- a[5] = -80. / (120. * 11.);
- a[6] = 120. / (720. * 13.);
- a[7] = -168. / (5040. * 15.);
- a[8] = 16. * 14. / (40320. * 17.);
- a[9] = -18. * 16. / (362880. * 19.);
-
+ if(crsdd!=NULL){
+ *crsdd =
+ tp * fi8 * (c[2] +
+ xx * (c[3] +
+ xx * (c[4] +
+ xx * (c[5] +
+ xx * (c[6] +
+ xx * (c[7] +
+ xx * (c[8] +
+ xx * c[9])))))));
+ }
+ } else {
+ tmp1 = erf(x) / x;
+ *crs = tmp1 / 2. - tp;
+ if((crsd!=NULL)||(crsdr2!=NULL)||(crsdd!=NULL)){
r = 2. * x / fi;
r2 = r * r;
- fi8 = fi * fi * fi * fi / 16.;
-
- if (x < XA) {
- xx = x * x;
- res =
- tp * fi8 * (a[2] +
- xx * (a[3] +
- xx * (a[4] +
- xx * (a[5] +
- xx * (a[6] +
- xx * (a[7] +
- xx * (a[8] +
- xx * a[9])))))));
+ tmp2 = exp(-x * x);
}
- else {
+ if(crsd!=NULL){ *crsd = (2. * tp * tmp2 - tmp1) / (2. * r2); }
+ if(crsdr2!=NULL){ *crsdr2 = *crsd / r2; }
/* Jarov vzorec */
- res =
- (erf(x) / (x * r2) -
- exp(-x * x) * (2 / r2 + fi * fi / 2) * tp) / (r2 * r2);;
+ if(crsdd!=NULL){ *crsdd = (tmp1 / r2 - tmp2 * (2 / r2 + fi * fi / 2) * tp) / (r2 * r2); }
}
- return (res);
+ return;
}
-
-
-
-
/*********solution of system of lin. equations*********/
int LINEQS(int DIM1, int N1, int N2, int *NERROR, double *DETERM)
More information about the grass-commit
mailing list