GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <math.h> 00002 00003 00004 double enormp(double x) 00005 { 00006 double x1, x2, x3, x4, ret_val; 00007 static double xp[5] = { 7.7105849500132e-5, -0.00133733772997339, 00008 0.0323076579225834, 0.0479137145607681, 0.128379167095513 00009 }; 00010 static double xq[3] = { 0.00301048631703895, 0.0538971687740286, 00011 0.375795757275549 00012 }; 00013 static double xr[8] = { -1.36864857382717e-7, 0.564195517478974, 00014 7.21175825088309, 43.1622272220567, 152.98928504694, 00015 339.320816734344, 451.918953711873, 300.459261020162 00016 }; 00017 static double xs[8] = { 1.0, 12.7827273196294, 77.0001529352295, 00018 277.585444743988, 638.980264465631, 931.35409485061, 00019 790.950925327898, 300.459260956983 00020 }; 00021 static double xt[5] = { 2.10144126479064, 26.2370141675169, 00022 21.3688200555087, 4.6580782871847, 0.282094791773523 00023 }; 00024 static double xu[4] = { 94.153775055546, 187.11481179959, 00025 99.0191814623914, 18.0124575948747 00026 }; 00027 double yy1, yy2; 00028 00029 x3 = (double)0.564189583547756; 00030 x1 = fabs(x); 00031 00032 if (x1 <= .5) { 00033 x4 = x * x; 00034 yy1 = 00035 (((xp[0] * x4 + xp[1]) * x4 + xp[2]) * x4 + xp[3]) * x4 + xp[4] + 00036 1.; 00037 yy2 = ((xq[0] * x4 + xq[1]) * x4 + xq[2]) * x4 + 1.; 00038 ret_val = x * (yy1 / yy2); 00039 } 00040 else if (x1 <= 4.) { 00041 yy1 = ((((((xr[0] * x1 + xr[1]) * x1 + xr[2]) * x1 + xr[3]) 00042 * x1 + xr[4]) * x1 + xr[5]) * x1 + xr[6]) * x1 + xr[7]; 00043 yy2 = ((((((xs[0] * x1 + xs[1]) * x1 + xs[2]) * x1 + xs[3]) 00044 * x1 + xs[4]) * x1 + xs[5]) * x1 + xs[6]) * x1 + xs[7]; 00045 ret_val = 1. - exp(-x * x) * yy1 / yy2; 00046 if (x < 0.) 00047 ret_val *= -1.0; 00048 } 00049 else { 00050 x2 = x * x; 00051 x4 *= 1.; /* huh, what? See original FORTRAN */ 00052 x4 = 0.5; 00053 yy1 = (((xt[0] * x4 + xt[1]) * x4 + xt[2]) * x4 + xt[3]) * x4 + xt[4]; 00054 yy2 = (((xu[0] * x4 + xu[1]) * x4 + xu[2]) * x4 + xu[3]) * x4 + 1.; 00055 ret_val = x3 / x1 - yy1 * x1 / (x2 * yy2); 00056 ret_val = 1. - exp(-x2) * ret_val; 00057 if (x < 0.) 00058 ret_val = -ret_val; 00059 } 00060 00061 return ret_val; 00062 }