GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /*- 00003 * 00004 * Original program and various modifications: 00005 * Lubos Mitas 00006 * 00007 * GRASS4.1 version of the program and GRASS4.2 modifications: 00008 * H. Mitasova, 00009 * I. Kosinovsky, D. Gerdes 00010 * D. McCauley 00011 * 00012 * Copyright 1993, 1995: 00013 * L. Mitas , 00014 * H. Mitasova , 00015 * I. Kosinovsky, 00016 * D.Gerdes 00017 * D. McCauley 00018 * 00019 * modified by McCauley in August 1995 00020 * modified by Mitasova in August 1995, Nov. 1996 00021 * 00022 */ 00023 00024 #include <stdio.h> 00025 #include <math.h> 00026 #include <grass/gis.h> 00027 #include <grass/interpf.h> 00028 00029 double IL_crst(double r, double fi) 00030 /* 00031 * Radial basis function - completely regularized spline with 00032 * tension (d=2) 00033 */ 00034 { 00035 double rfsta2 = fi * fi * r / 4.; 00036 00037 static double c[4] = { 8.5733287401, 18.0590169730, 8.6347608925, 00038 0.2677737343 00039 }; 00040 static double b[4] = { 9.5733223454, 25.6329561486, 21.0996530827, 00041 3.9584969228 00042 }; 00043 double ce = 0.57721566; 00044 00045 static double u[10] = { 1.e+00, -.25e+00, 00046 .055555555555556e+00, -.010416666666667e+00, /*fixed bug 415.. repl. by 416.. */ 00047 .166666666666667e-02, -2.31481481481482e-04, 00048 2.83446712018141e-05, -3.10019841269841e-06, 00049 3.06192435822065e-07, -2.75573192239859e-08 00050 }; 00051 double x = rfsta2; 00052 double res; 00053 00054 double e1, ea, eb; 00055 00056 00057 if (x < 1.e+00) { 00058 res = x * (u[0] + x * (u[1] + x * (u[2] + x * (u[3] + x * (u[4] + x * 00059 (u[5] + 00060 x * 00061 (u[6] + 00062 x * 00063 (u[7] + 00064 x * 00065 (u[8] + 00066 x * 00067 u 00068 [9]))))))))); 00069 return (res); 00070 } 00071 00072 if (x > 25.e+00) 00073 e1 = 0.00; 00074 else { 00075 ea = c[3] + x * (c[2] + x * (c[1] + x * (c[0] + x))); 00076 eb = b[3] + x * (b[2] + x * (b[1] + x * (b[0] + x))); 00077 e1 = (ea / eb) / (x * exp(x)); 00078 } 00079 res = e1 + ce + log(x); 00080 return (res); 00081 } 00082 00083 00084 00085 int IL_crstg(double r, double fi, double *gd1, /* G1(r) */ 00086 double *gd2 /* G2(r) */ 00087 ) 00088 00089 00090 /* 00091 * Function for calculating derivatives (d=2) 00092 */ 00093 { 00094 double r2 = r; 00095 double rfsta2 = fi * fi * r / 4.; 00096 double x, exm, oneme, hold; 00097 double fsta2 = fi * fi / 2.; 00098 00099 x = rfsta2; 00100 if (x < 0.001) { 00101 *gd1 = 1. - x / 2. + x * x / 6. - x * x * x / 24.; 00102 *gd2 = fsta2 * (-.5 + x / 3. - x * x / 8. + x * x * x / 30.); 00103 } 00104 else { 00105 if (x < 35.e+00) { 00106 exm = exp(-x); 00107 oneme = 1. - exm; 00108 *gd1 = oneme / x; 00109 hold = x * exm - oneme; 00110 *gd2 = (hold + hold) / (r2 * x); 00111 } 00112 else { 00113 *gd1 = 1. / x; 00114 *gd2 = -2. / (x * r2); 00115 } 00116 } 00117 return 1; 00118 }