GRASS GIS 7 Programmer's Manual  7.0.2(2015)-r00000
func2d.c
Go to the documentation of this file.
1 
35 #include <stdio.h>
36 #include <math.h>
37 #include <grass/gis.h>
38 #include <grass/interpf.h>
39 
40 
41 /* parameter description from DESCRIPTION.INTERP */
48 double IL_crst(double r,
49  double fi
50  )
51 {
52  double rfsta2 = fi * fi * r / 4.;
53 
54  static double c[4] = { 8.5733287401, 18.0590169730, 8.6347608925,
55  0.2677737343
56  };
57  static double b[4] = { 9.5733223454, 25.6329561486, 21.0996530827,
58  3.9584969228
59  };
60  double ce = 0.57721566;
61 
62  static double u[10] = { 1.e+00, -.25e+00,
63  .055555555555556e+00, -.010416666666667e+00, /*fixed bug 415.. repl. by 416.. */
64  .166666666666667e-02, -2.31481481481482e-04,
65  2.83446712018141e-05, -3.10019841269841e-06,
66  3.06192435822065e-07, -2.75573192239859e-08
67  };
68  double x = rfsta2;
69  double res;
70 
71  double e1, ea, eb;
72 
73 
74  if (x < 1.e+00) {
75  res = x * (u[0] + x * (u[1] + x * (u[2] + x * (u[3] + x * (u[4] + x *
76  (u[5] +
77  x *
78  (u[6] +
79  x *
80  (u[7] +
81  x *
82  (u[8] +
83  x *
84  u
85  [9])))))))));
86  return (res);
87  }
88 
89  if (x > 25.e+00)
90  e1 = 0.00;
91  else {
92  ea = c[3] + x * (c[2] + x * (c[1] + x * (c[0] + x)));
93  eb = b[3] + x * (b[2] + x * (b[1] + x * (b[0] + x)));
94  e1 = (ea / eb) / (x * exp(x));
95  }
96  res = e1 + ce + log(x);
97  return (res);
98 }
99 
100 
106 int IL_crstg(double r,
107  double fi,
108  double *gd1,
109  double *gd2
110  )
111 {
112  double r2 = r;
113  double rfsta2 = fi * fi * r / 4.;
114  double x, exm, oneme, hold;
115  double fsta2 = fi * fi / 2.;
116 
117  x = rfsta2;
118  if (x < 0.001) {
119  *gd1 = 1. - x / 2. + x * x / 6. - x * x * x / 24.;
120  *gd2 = fsta2 * (-.5 + x / 3. - x * x / 8. + x * x * x / 30.);
121  }
122  else {
123  if (x < 35.e+00) {
124  exm = exp(-x);
125  oneme = 1. - exm;
126  *gd1 = oneme / x;
127  hold = x * exm - oneme;
128  *gd2 = (hold + hold) / (r2 * x);
129  }
130  else {
131  *gd1 = 1. / x;
132  *gd2 = -2. / (x * r2);
133  }
134  }
135  return 1;
136 }
int IL_crstg(double r, double fi, double *gd1, double *gd2)
Definition: func2d.c:106
double IL_crst(double r, double fi)
Definition: func2d.c:48
double b
double r