GRASS Programmer's Manual  6.4.2(2012)
func2d.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines