![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00022 #include "config.h" 00023 00024 #include <stdio.h> 00025 #include <math.h> 00026 #include <float.h> 00027 #ifdef HAVE_COMPLEX_H 00028 #include <complex.h> 00029 #endif 00030 00031 #include "kernels.h" 00032 00038 C gaussian(R x, int der, const R *param) /* K(x)=EXP(-x^2/c^2) */ 00039 { 00040 R c = param[0]; 00041 R value = K(0.0); 00042 00043 switch (der) 00044 { 00045 case 0 : value = EXP(-x*x/(c*c)); break; 00046 case 1 : value = -K(2.0)*x/(c*c)*EXP(-x*x/(c*c)); break; 00047 case 2 : value = K(2.0)*EXP(-x*x/(c*c))*(-c*c+K(2.0)*x*x)/(c*c*c*c); break; 00048 case 3 : value = -K(4.0)*x*EXP(-x*x/(c*c))*(-K(3.0)*c*c+K(2.0)*x*x)/(c*c*c*c*c*c); break; 00049 case 4 : value = K(4.0)*EXP(-x*x/(c*c))*(K(3.0)*c*c*c*c-K(12.0)*c*c*x*x+K(4.0)*x*x*x*x)/(c*c*c*c*c*c*c*c); break; 00050 case 5 : value = -K(8.0)*x*EXP(-x*x/(c*c))*(K(15.0)*c*c*c*c-K(20.0)*c*c*x*x+K(4.0)*x*x*x*x)/POW(c,K(10.0)); break; 00051 case 6 : value = K(8.0)*EXP(-x*x/(c*c))*(-K(15.0)*c*c*c*c*c*c+K(90.0)*x*x*c*c*c*c-K(60.0)*x*x*x*x*c*c+K(8.0)*x*x*x*x*x*x)/POW(c,K(12.0)); break; 00052 case 7 : value = -K(16.0)*x*EXP(-x*x/(c*c))*(-K(105.0)*c*c*c*c*c*c+K(210.0)*x*x*c*c*c*c-K(84.0)*x*x*x*x*c*c+K(8.0)*x*x*x*x*x*x)/POW(c,K(14.0)); break; 00053 case 8 : value = K(16.0)*EXP(-x*x/(c*c))*(K(105.0)*c*c*c*c*c*c*c*c-K(840.0)*x*x*c*c*c*c*c*c+K(840.0)*x*x*x*x*c*c*c*c-K(224.0)*x*x*x*x*x*x*c*c+K(16.0)*x*x*x*x*x*x*x*x)/POW(c,K(16.0)); break; 00054 case 9 : value = -K(32.0)*x*EXP(-x*x/(c*c))*(K(945.0)*c*c*c*c*c*c*c*c-K(2520.0)*x*x*c*c*c*c*c*c+K(1512.0)*x*x*x*x*c*c*c*c-K(288.0)*x*x*x*x*x*x*c*c+K(16.0)*x*x*x*x*x*x*x*x)/POW(c,K(18.0)); break; 00055 case 10 : value = K(32.0)*EXP(-x*x/(c*c))*(-K(945.0)*POW(c,K(10.0))+K(9450.0)*x*x*c*c*c*c*c*c*c*c-K(12600.0)*x*x*x*x*c*c*c*c*c*c+K(5040.0)*x*x*x*x*x*x*c*c*c*c-K(720.0)*x*x*x*x*x*x*x*x*c*c+K(32.0)*POW(x,K(10.0)))/POW(c,K(20.0)); break; 00056 case 11 : value = -K(64.0)*x*EXP(-x*x/(c*c))*(-K(10395.0)*POW(c,K(10.0))+K(34650.0)*x*x*c*c*c*c*c*c*c*c-K(27720.0)*x*x*x*x*c*c*c*c*c*c+K(7920.0)*x*x*x*x*x*x*c*c*c*c-K(880.0)*x*x*x*x*x*x*x*x*c*c+K(32.0)*POW(x,K(10.0)))/POW(c,K(22.0)); break; 00057 case 12 : value = K(64.0)*EXP(-x*x/(c*c))*(K(10395.0)*POW(c,K(12.0))-K(124740.0)*x*x*POW(c,K(10.0))+K(207900.0)*x*x*x*x*c*c*c*c*c*c*c*c-K(110880.0)*x*x*x*x*x*x*c*c*c*c*c*c+K(23760.0)*x*x*x*x*x*x*x*x*c*c*c*c-K(2112.0)*POW(x,K(10.0))*c*c+K(64.0)*POW(x,K(12.0)))/POW(c,K(24.0)); break; 00058 default : value = K(0.0); 00059 } 00060 00061 return value; 00062 } 00063 00064 C multiquadric(R x, int der, const R *param) /* K(x)=SQRT(x^2+c^2) */ 00065 { 00066 R c=param[0]; 00067 R value=K(0.0); 00068 00069 switch (der) 00070 { 00071 case 0 : value=SQRT(x*x+c*c); break; 00072 case 1 : value=K(1.0)/(SQRT(x*x+c*c))*x; break; 00073 case 2 : value=c*c/SQRT(POW(x*x+c*c,K(3.0))); break; 00074 case 3 : value=-K(3.0)*x*c*c/SQRT(POW(x*x+c*c,K(5.0))); break; 00075 case 4 : value=K(3.0)*c*c*(K(4.0)*x*x-c*c)/SQRT(POW(x*x+c*c,K(7.0))); break; 00076 case 5 : value=-K(15.0)*x*c*c*(K(4.0)*x*x-K(3.0)*c*c)/SQRT(POW(x*x+c*c,K(9.0))); break; 00077 case 6 : value=K(45.0)*c*c*(K(8.0)*x*x*x*x-K(12.0)*x*x*c*c+c*c*c*c)/SQRT(POW(x*x+c*c,K(11.0))); break; 00078 case 7 : value=-K(315.0)*x*c*c*(K(8.0)*x*x*x*x-K(20.0)*x*x*c*c+K(5.0)*c*c*c*c)/SQRT(POW(x*x+c*c,K(13.0))); break; 00079 case 8 : value=K(315.0)*c*c*(K(64.0)*x*x*x*x*x*x-K(240.0)*x*x*x*x*c*c+K(120.0)*x*x*c*c*c*c-K(5.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(15.0))); break; 00080 case 9 : value=-K(2835.0)*x*c*c*(K(64.0)*x*x*x*x*x*x-K(336.0)*x*x*x*x*c*c+K(280.0)*x*x*c*c*c*c-K(35.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(17.0))); break; 00081 case 10 : value=K(14175.0)*c*c*(K(128.0)*x*x*x*x*x*x*x*x-K(896.0)*x*x*x*x*x*x*c*c+K(1120.0)*x*x*x*x*c*c*c*c-K(280.0)*x*x*c*c*c*c*c*c+K(7.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(19.0))); break; 00082 case 11 : value=-K(155925.0)*x*c*c*(K(128.0)*x*x*x*x*x*x*x*x-K(1152.0)*x*x*x*x*x*x*c*c+K(2016.0)*x*x*x*x*c*c*c*c-K(840.0)*x*x*c*c*c*c*c*c+K(63.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(21.0))); break; 00083 case 12 : value=K(467775.0)*c*c*(K(1260.0)*x*x*c*c*c*c*c*c*c*c-K(21.0)*POW(c,K(10.0))+K(512.0)*POW(x,K(10.0))-K(5760.0)*x*x*x*x*x*x*x*x*c*c+K(13440.0)*x*x*x*x*x*x*c*c*c*c-K(8400.0)*x*x*x*x*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(23.0))); break; 00084 default : value=K(0.0); 00085 } 00086 00087 return value; 00088 } 00089 00090 C inverse_multiquadric(R x, int der, const R *param) /* K(x)=1/SQRT(x^2+c^2) */ 00091 { 00092 R c=param[0]; 00093 R value=K(0.0); 00094 00095 switch (der) 00096 { 00097 case 0 : value=K(1.0)/SQRT(x*x+c*c); break; 00098 case 1 : value=-K(1.0)/(SQRT(POW(x*x+c*c,K(3.0))))*x; break; 00099 case 2 : value=(K(2.0)*x*x-c*c)/SQRT(POW(x*x+c*c,K(5.0))); break; 00100 case 3 : value=-K(3.0)*x*(K(2.0)*x*x-K(3.0)*c*c)/SQRT(POW(x*x+c*c,K(7.0))); break; 00101 case 4 : value=K(3.0)*(K(8.0)*x*x*x*x-K(24.0)*x*x*c*c+K(3.0)*c*c*c*c)/SQRT(POW(x*x+c*c,K(9.0))); break; 00102 case 5 : value=-K(15.0)*x*(K(8.0)*x*x*x*x-K(40.0)*x*x*c*c+K(15.0)*c*c*c*c)/SQRT(POW(x*x+c*c,K(11.0))); break; 00103 case 6 : value=K(45.0)*(K(16.0)*x*x*x*x*x*x-K(120.0)*x*x*x*x*c*c+K(90.0)*x*x*c*c*c*c-K(5.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(13.0))); break; 00104 case 7 : value=-K(315.0)*x*(K(16.0)*x*x*x*x*x*x-K(168.0)*x*x*x*x*c*c+K(210.0)*x*x*c*c*c*c-K(35.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(15.0))); break; 00105 case 8 : value=K(315.0)*(K(128.0)*x*x*x*x*x*x*x*x-K(1792.0)*x*x*x*x*x*x*c*c+K(3360.0)*x*x*x*x*c*c*c*c-K(1120.0)*x*x*c*c*c*c*c*c+K(35.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(17.0))); break; 00106 case 9 : value=-K(2835.0)*x*(K(128.0)*x*x*x*x*x*x*x*x-K(2304.0)*x*x*x*x*x*x*c*c+K(6048.0)*x*x*x*x*c*c*c*c-K(3360.0)*x*x*c*c*c*c*c*c+K(315.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(19.0))); break; 00107 case 10 : value=K(14175.0)*(K(256.0)*POW(x,K(10.0))-K(5760.0)*x*x*x*x*x*x*x*x*c*c+K(20160.0)*x*x*x*x*x*x*c*c*c*c-K(16800.0)*x*x*x*x*c*c*c*c*c*c+K(3150.0)*x*x*c*c*c*c*c*c*c*c-K(63.0)*POW(c,K(10.0)))/SQRT(POW(x*x+c*c,K(21.0))); break; 00108 case 11 : value=-K(155925.0)*x*(K(256.0)*POW(x,K(10.0))-K(7040.0)*x*x*x*x*x*x*x*x*c*c+K(31680.0)*x*x*x*x*x*x*c*c*c*c-K(36960.0)*x*x*x*x*c*c*c*c*c*c+K(11550.0)*x*x*c*c*c*c*c*c*c*c-K(693.0)*POW(c,K(10.0)))/SQRT(POW(x*x+c*c,K(23.0))); break; 00109 case 12 : value=K(467775.0)*(K(231.0)*POW(c,K(12.0))+K(190080.0)*x*x*x*x*x*x*x*x*c*c*c*c-K(16632.0)*x*x*POW(c,K(10.0))-K(295680.0)*x*x*x*x*x*x*c*c*c*c*c*c+K(138600.0)*x*x*x*x*c*c*c*c*c*c*c*c+K(1024.0)*POW(x,K(12.0))-K(33792.0)*POW(x,K(10.0))*c*c)/SQRT(POW(x*x+c*c,K(25.0))); break; 00110 default : value=K(0.0); 00111 } 00112 00113 return value; 00114 } 00115 00116 C logarithm(R x, int der, const R *param) /* K(x)=LOG |x| */ 00117 { 00118 R value=K(0.0); 00119 00120 (void)param; 00121 00122 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00123 else switch (der) 00124 { 00125 case 0 : value=LOG(FABS(x)); break; 00126 case 1 : value=(x<0 ? -1 : 1)/FABS(x); break; 00127 case 2 : value=-1/(x*x); break; 00128 case 3 : value=K(2.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(3.0)); break; 00129 case 4 : value=-K(6.0)/(x*x*x*x); break; 00130 case 5 : value=K(24.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(5.0)); break; 00131 case 6 : value=-K(120.0)/(x*x*x*x*x*x); break; 00132 case 7 : value=K(720.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(7.0)); break; 00133 case 8 : value=-K(5040.0)/(x*x*x*x*x*x*x*x); break; 00134 case 9 : value=K(40320.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(9.0)); break; 00135 case 10 : value=-K(362880.0)/POW(x,K(10.0)); break; 00136 case 11 : value=K(3628800.0)*(x<0 ? -1 : 1)/POW(FABS(x),K(11.0)); break; 00137 case 12 : value=-K(39916800.0)/POW(x,K(12.0)); break; 00138 case 13 : value=K(479001600.0)/POW(x,K(13.0)); break; 00139 case 14 : value=-K(6227020800.0)/POW(x,K(14.0)); break; 00140 case 15 : value=K(87178291200.0)/POW(x,K(15.0)); break; 00141 case 16 : value=-K(1307674368000.0)/POW(x,K(16.0)); break; 00142 case 17 : value=K(20922789888000.0)/POW(x,K(17.0)); break; 00143 default : value=K(0.0); 00144 } 00145 00146 return value; 00147 } 00148 00149 C thinplate_spline(R x, int der, const R *param) /* K(x) = x^2 LOG |x| */ 00150 { 00151 R value=K(0.0); 00152 00153 (void)param; 00154 00155 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00156 else switch (der) 00157 { 00158 case 0 : value=x*x*LOG(FABS(x)); break; 00159 case 1 : value=K(2.0)*x*LOG(FABS(x))+x; break; 00160 case 2 : value=K(2.0)*LOG(FABS(x))+K(3.0); break; 00161 case 3 : value=K(2.0)/x; break; 00162 case 4 : value=-K(2.0)/(x*x); break; 00163 case 5 : value=K(4.0)/(x*x*x); break; 00164 case 6 : value=-K(12.0)/(x*x*x*x); break; 00165 case 7 : value=K(48.0)/(x*x*x*x*x); break; 00166 case 8 : value=-K(240.0)/(x*x*x*x*x*x); break; 00167 case 9 : value=K(1440.0)/(x*x*x*x*x*x*x); break; 00168 case 10 : value=-K(10080.0)/(x*x*x*x*x*x*x*x); break; 00169 case 11 : value=K(80640.0)/(x*x*x*x*x*x*x*x*x); break; 00170 case 12 : value=-K(725760.0)/POW(x,K(10.0)); break; 00171 default : value=K(0.0); 00172 } 00173 00174 return value; 00175 } 00176 00177 C one_over_square(R x, int der, const R *param) /* K(x) = 1/x^2 */ 00178 { 00179 R value=K(0.0); 00180 00181 (void)param; 00182 00183 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00184 else switch (der) 00185 { 00186 case 0 : value=K(1.0)/(x*x); break; 00187 case 1 : value=-K(2.0)/(x*x*x); break; 00188 case 2 : value=K(6.0)/(x*x*x*x); break; 00189 case 3 : value=-K(24.0)/(x*x*x*x*x); break; 00190 case 4 : value=K(120.0)/(x*x*x*x*x*x); break; 00191 case 5 : value=-K(720.0)/(x*x*x*x*x*x*x); break; 00192 case 6 : value=K(5040.0)/(x*x*x*x*x*x*x*x); break; 00193 case 7 : value=-K(40320.0)/(x*x*x*x*x*x*x*x*x); break; 00194 case 8 : value=K(362880.0)/POW(x,K(10.0)); break; 00195 case 9 : value=-K(3628800.0)/POW(x,K(11.0)); break; 00196 case 10 : value=K(39916800.0)/POW(x,K(12.0)); break; 00197 case 11 : value=-K(479001600.0)/POW(x,K(13.0)); break; 00198 case 12 : value=K(6227020800.0)/POW(x,K(14.0)); break; 00199 default : value=K(0.0); 00200 } 00201 00202 return value; 00203 } 00204 00205 C one_over_modulus(R x, int der, const R *param) /* K(x) = 1/|x| */ 00206 { 00207 R value=K(0.0); 00208 00209 (void)param; 00210 00211 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00212 else switch (der) 00213 { 00214 case 0 : value=K(1.0)/FABS(x); break; 00215 case 1 : value=-1/x/FABS(x); break; 00216 case 2 : value=K(2.0)/POW(FABS(x),K(3.0)); break; 00217 case 3 : value=-K(6.0)/(x*x*x)/FABS(x); break; 00218 case 4 : value=K(24.0)/POW(FABS(x),K(5.0)); break; 00219 case 5 : value=-K(120.0)/(x*x*x*x*x)/FABS(x); break; 00220 case 6 : value=K(720.0)/POW(FABS(x),K(7.0)); break; 00221 case 7 : value=-K(5040.0)/(x*x*x*x*x*x*x)/FABS(x); break; 00222 case 8 : value=K(40320.0)/POW(FABS(x),K(9.0)); break; 00223 case 9 : value=-K(362880.0)/(x*x*x*x*x*x*x*x*x)/FABS(x); break; 00224 case 10 : value=K(3628800.0)/POW(FABS(x),K(11.0)); break; 00225 case 11 : value=-K(39916800.0)/POW(x,K(11.0))/FABS(x); break; 00226 case 12 : value=K(479001600.0)/POW(FABS(x),K(13.0)); break; 00227 default : value=K(0.0); 00228 } 00229 00230 return value; 00231 } 00232 00233 C one_over_x(R x, int der, const R *param) /* K(x) = 1/x */ 00234 { 00235 R value=K(0.0); 00236 00237 (void)param; 00238 00239 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00240 else switch (der) 00241 { 00242 case 0 : value=K(1.0)/x; break; 00243 case 1 : value=-K(1.0)/(x*x); break; 00244 case 2 : value=K(2.0)/(x*x*x); break; 00245 case 3 : value=-K(6.0)/(x*x*x*x); break; 00246 case 4 : value=K(24.0)/(x*x*x*x*x); break; 00247 case 5 : value=-K(120.0)/(x*x*x*x*x*x); break; 00248 case 6 : value=K(720.0)/(x*x*x*x*x*x*x); break; 00249 case 7 : value=-K(5040.0)/(x*x*x*x*x*x*x*x); break; 00250 case 8 : value=K(40320.0)/(x*x*x*x*x*x*x*x*x); break; 00251 case 9 : value=-K(362880.0)/POW(x,K(10.0)); break; 00252 case 10 : value=K(3628800.0)/POW(x,K(11.0)); break; 00253 case 11 : value=-K(39916800.0)/POW(x,K(12.0)); break; 00254 case 12 : value=K(479001600.0)/POW(x,K(13.0)); break; 00255 default : value=K(0.0); 00256 } 00257 00258 return value; 00259 } 00260 00261 C inverse_multiquadric3(R x, int der, const R *param) /* K(x) = 1/SQRT(x^2+c^2)^3 */ 00262 { 00263 R c=param[0]; 00264 R value=K(0.0); 00265 00266 switch (der) 00267 { 00268 case 0 : value=K(1.0)/(SQRT(POW(x*x+c*c,K(3.0)))); break; 00269 case 1 : value=-K(3.0)/SQRT(POW(x*x+c*c,K(5.0)))*x; break; 00270 case 2 : value=K(3.0)*(K(4.0)*x*x-c*c)/SQRT(POW(x*x+c*c,K(7.0))); break; 00271 case 3 : value=-K(15.0)*x*(K(4.0)*x*x-K(3.0)*c*c)/SQRT(POW(x*x+c*c,K(9.0))); break; 00272 case 4 : value=K(45.0)*(K(8.0)*x*x*x*x-K(12.0)*x*x*c*c+c*c*c*c)/SQRT(POW(x*x+c*c,K(11.0))); break; 00273 case 5 : value=-K(315.0)*x*(K(8.0)*x*x*x*x-K(20.0)*x*x*c*c+K(5.0)*c*c*c*c)/SQRT(POW(x*x+c*c,K(13.0))); break; 00274 case 6 : value=K(315.0)*(K(64.0)*x*x*x*x*x*x-K(240.0)*x*x*x*x*c*c+K(120.0)*x*x*c*c*c*c-K(5.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(15.0))); break; 00275 case 7 : value=-K(2835.0)*x*(K(64.0)*x*x*x*x*x*x-K(336.0)*x*x*x*x*c*c+K(280.0)*x*x*c*c*c*c-K(35.0)*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(17.0))); break; 00276 case 8 : value=K(14175.0)*(K(128.0)*x*x*x*x*x*x*x*x-K(896.0)*x*x*x*x*x*x*c*c+K(1120.0)*x*x*x*x*c*c*c*c-K(280.0)*x*x*c*c*c*c*c*c+K(7.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(19.0))); break; 00277 case 9 : value=-K(155925.0)*x*(K(128.0)*x*x*x*x*x*x*x*x-K(1152.0)*x*x*x*x*x*x*c*c+K(2016.0)*x*x*x*x*c*c*c*c-K(840.0)*x*x*c*c*c*c*c*c+K(63.0)*c*c*c*c*c*c*c*c)/SQRT(POW(x*x+c*c,K(21.0))); break; 00278 case 10 : value=K(467775.0)*(K(512.0)*POW(x,K(10.0))-K(5760.0)*x*x*x*x*x*x*x*x*c*c+K(13440.0)*x*x*x*x*x*x*c*c*c*c-K(8400.0)*x*x*x*x*c*c*c*c*c*c+K(1260.0)*x*x*c*c*c*c*c*c*c*c-K(21.0)*POW(c,K(10.0)))/SQRT(POW(x*x+c*c,K(23.0))); break; 00279 case 11 : value=-K(6081075.0)*x*(K(512.0)*POW(x,K(10.0))-K(7040.0)*x*x*x*x*x*x*x*x*c*c+K(21120.0)*x*x*x*x*x*x*c*c*c*c-K(18480.0)*x*x*x*x*c*c*c*c*c*c+K(4620.0)*x*x*c*c*c*c*c*c*c*c-K(231.0)*POW(c,K(10.0)))/SQRT(POW(x*x+c*c,K(25.0))); break; 00280 case 12 : value=K(42567525.0)*(K(1024.0)*POW(x,K(12.0))+K(27720.0)*x*x*x*x*c*c*c*c*c*c*c*c+K(33.0)*POW(c,K(12.0))-K(2772.0)*x*x*POW(c,K(10.0))-K(73920.0)*x*x*x*x*x*x*c*c*c*c*c*c+K(63360.0)*x*x*x*x*x*x*x*x*c*c*c*c-K(16896.0)*POW(x,K(10.0))*c*c)/SQRT(POW(x*x+c*c,K(27.0))); break; 00281 default : value=K(0.0); 00282 } 00283 00284 return value; 00285 } 00286 00287 C sinc_kernel(R x, int der, const R *param) /* K(x) = SIN(cx)/x */ 00288 { 00289 R c=param[0]; 00290 R value=K(0.0); 00291 00292 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00293 else switch (der) 00294 { 00295 case 0 : value=SIN(c*x)/x; break; 00296 case 1 : value=(COS(c*x)*c*x-SIN(c*x))/(x*x); break; 00297 case 2 : value=-(SIN(c*x)*c*c*x*x+K(2.0)*COS(c*x)*c*x-K(2.0)*SIN(c*x))/(x*x*x); break; 00298 case 3 : value=-(COS(c*x)*c*c*c*x*x*x-K(3.0)*SIN(c*x)*c*c*x*x-K(6.0)*COS(c*x)*c*x+K(6.0)*SIN(c*x))/(x*x*x*x); break; 00299 case 4 : value=(SIN(c*x)*c*c*c*c*x*x*x*x+K(4.0)*COS(c*x)*c*c*c*x*x*x-K(12.0)*SIN(c*x)*c*c*x*x-K(24.0)*COS(c*x)*c*x+K(24.0)*SIN(c*x))/(x*x*x*x*x); break; 00300 case 5 : value=(COS(c*x)*c*c*c*c*c*x*x*x*x*x-K(5.0)*SIN(c*x)*c*c*c*c*x*x*x*x-K(20.0)*COS(c*x)*c*c*c*x*x*x+K(60.0)*SIN(c*x)*c*c*x*x+K(120.0)*COS(c*x)*c*x-K(120.0)*SIN(c*x))/(x*x*x*x*x*x); break; 00301 case 6 : value=-(SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(6.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x-K(30.0)*SIN(c*x)*c*c*c*c*x*x*x*x-K(120.0)*COS(c*x)*c*c*c*x*x*x+K(360.0)*SIN(c*x)*c*c*x*x+K(720.0)*COS(c*x)*c*x-K(720.0)*SIN(c*x))/(x*x*x*x*x*x*x); break; 00302 case 7 : value=-(COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(7.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(42.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x+K(210.0)*SIN(c*x)*c*c*c*c*x*x*x*x+K(840.0)*COS(c*x)*c*c*c*x*x*x-K(2520.0)*SIN(c*x)*c*c*x*x-K(5040.0)*COS(c*x)*c*x+K(5040.0)*SIN(c*x))/(x*x*x*x*x*x*x*x); break; 00303 case 8 : value=(SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(8.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(56.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(336.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x+K(1680.0)*SIN(c*x)*c*c*c*c*x*x*x*x+K(6720.0)*COS(c*x)*c*c*c*x*x*x-K(20160.0)*SIN(c*x)*c*c*x*x-K(40320.0)*COS(c*x)*c*x+K(40320.0)*SIN(c*x))/(x*x*x*x*x*x*x*x*x); break; 00304 case 9 : value=(COS(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-K(9.0)*SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(72.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(504.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(3024.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x-K(15120.0)*SIN(c*x)*c*c*c*c*x*x*x*x-K(60480.0)*COS(c*x)*c*c*c*x*x*x+K(181440.0)*SIN(c*x)*c*c*x*x+K(362880.0)*COS(c*x)*c*x-K(362880.0)*SIN(c*x))/POW(x,K(10.0)); break; 00305 case 10 : value=-(SIN(c*x)*POW(c,K(10.0))*POW(x,K(10.0))+K(10.0)*COS(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-K(90.0)*SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(720.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(5040.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(30240.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x-K(151200.0)*SIN(c*x)*c*c*c*c*x*x*x*x-K(604800.0)*COS(c*x)*c*c*c*x*x*x+K(1814400.0)*SIN(c*x)*c*c*x*x+K(3628800.0)*COS(c*x)*c*x-K(3628800.0)*SIN(c*x))/POW(x,K(11.0)); break; 00306 case 11 : value=-(COS(c*x)*POW(c,K(11.0))*POW(x,K(11.0))-K(11.0)*SIN(c*x)*POW(c,K(10.0))*POW(x,K(10.0))-K(110.0)*COS(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+K(990.0)*SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(7920.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(55440.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(332640.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x+K(1663200.0)*SIN(c*x)*c*c*c*c*x*x*x*x+K(6652800.0)*COS(c*x)*c*c*c*x*x*x-K(19958400.0)*SIN(c*x)*c*c*x*x-K(39916800.0)*COS(c*x)*c*x+K(39916800.0)*SIN(c*x))/POW(x,K(12.0)); break; 00307 case 12 : value=(SIN(c*x)*POW(c,K(12.0))*POW(x,K(12.0))+K(12.0)*COS(c*x)*POW(c,K(11.0))*POW(x,K(11.0))-K(132.0)*SIN(c*x)*POW(c,K(10.0))*POW(x,K(10.0))-K(1320.0)*COS(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+K(11880.0)*SIN(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(95040.0)*COS(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(665280.0)*SIN(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(3991680.0)*COS(c*x)*c*c*c*c*c*x*x*x*x*x+K(19958400.0)*SIN(c*x)*c*c*c*c*x*x*x*x+K(79833600.0)*COS(c*x)*c*c*c*x*x*x-K(239500800.0)*SIN(c*x)*c*c*x*x-K(479001600.0)*COS(c*x)*c*x+K(479001600.0)*SIN(c*x))/POW(x,K(13.0)); break; 00308 default : value=K(0.0); 00309 } 00310 00311 return value; 00312 } 00313 00314 C cosc(R x, int der, const R *param) /* K(x) = COS(cx)/x */ 00315 { 00316 R c=param[0]; 00317 R value=K(0.0); 00318 R sign; 00319 00320 if (x<0) sign=-K(1.0); else sign=K(1.0); 00321 x=FABS(x); 00322 00323 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00324 else switch (der) 00325 { 00326 case 0 : value=COS(c*x)/x; break; 00327 case 1 : value=-(SIN(c*x)*c*x+COS(c*x))/(x*x); break; 00328 case 2 : value=(-COS(c*x)*c*c*x*x+K(2.0)*SIN(c*x)*c*x+K(2.0)*COS(c*x))/(x*x*x); break; 00329 case 3 : value=(SIN(c*x)*c*c*c*x*x*x+K(3.0)*COS(c*x)*c*c*x*x-K(6.0)*SIN(c*x)*c*x-K(6.0)*COS(c*x))/(x*x*x*x); break; 00330 case 4 : value=(COS(c*x)*c*c*c*c*x*x*x*x-K(4.0)*SIN(c*x)*c*c*c*x*x*x-K(12.0)*COS(c*x)*c*c*x*x+K(24.0)*SIN(c*x)*c*x+K(24.0)*COS(c*x))/(x*x*x*x*x); break; 00331 case 5 : value=-(SIN(c*x)*c*c*c*c*c*x*x*x*x*x+K(5.0)*COS(c*x)*c*c*c*c*x*x*x*x-K(20.0)*SIN(c*x)*c*c*c*x*x*x-K(60.0)*COS(c*x)*c*c*x*x+K(120.0)*SIN(c*x)*c*x+K(120.0)*COS(c*x))/(x*x*x*x*x*x); break; 00332 case 6 : value=-(COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(6.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x-K(30.0)*COS(c*x)*c*c*c*c*x*x*x*x+K(120.0)*SIN(c*x)*c*c*c*x*x*x+K(360.0)*COS(c*x)*c*c*x*x-K(720.0)*SIN(c*x)*c*x-K(720.0)*COS(c*x))/(x*x*x*x*x*x*x); break; 00333 case 7 : value=(SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(7.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(42.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x-K(210.0)*COS(c*x)*c*c*c*c*x*x*x*x+K(840.0)*SIN(c*x)*c*c*c*x*x*x+K(2520.0)*COS(c*x)*c*c*x*x-K(5040.0)*SIN(c*x)*c*x-K(5040.0)*COS(c*x))/(x*x*x*x*x*x*x*x); break; 00334 case 8 : value=(COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(8.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(56.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(336.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x+K(1680.0)*COS(c*x)*c*c*c*c*x*x*x*x-K(6720.0)*SIN(c*x)*c*c*c*x*x*x-K(20160.0)*COS(c*x)*c*c*x*x+K(40320.0)*SIN(c*x)*c*x+K(40320.0)*COS(c*x))/(x*x*x*x*x*x*x*x*x); break; 00335 case 9 : value=-(SIN(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+K(9.0)*COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(72.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(504.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(3024.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x+K(15120.0)*COS(c*x)*c*c*c*c*x*x*x*x-K(60480.0)*SIN(c*x)*c*c*c*x*x*x-K(181440.0)*COS(c*x)*c*c*x*x+K(362880.0)*SIN(c*x)*c*x+K(362880.0)*COS(c*x))/POW(x,K(10.0)); break; 00336 case 10 : value=-(COS(c*x)*POW(c,K(10.0))*POW(x,K(10.0))-K(10.0)*SIN(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-K(90.0)*COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(720.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(5040.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(30240.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x-K(151200.0)*COS(c*x)*c*c*c*c*x*x*x*x+K(604800.0)*SIN(c*x)*c*c*c*x*x*x+K(1814400.0)*COS(c*x)*c*c*x*x-K(3628800.0)*SIN(c*x)*c*x-K(3628800.0)*COS(c*x))/POW(x,K(11.0)); break; 00337 case 11 : value=(SIN(c*x)*POW(c,K(11.0))*POW(x,K(11.0))+K(11.0)*COS(c*x)*POW(c,K(10.0))*POW(x,K(10.0))-K(110.0)*SIN(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x-K(990.0)*COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x+K(7920.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x+K(55440.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x-K(332640.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x-K(1663200.0)*COS(c*x)*c*c*c*c*x*x*x*x+K(6652800.0)*SIN(c*x)*c*c*c*x*x*x+K(19958400.0)*COS(c*x)*c*c*x*x-K(39916800.0)*SIN(c*x)*c*x-K(39916800.0)*COS(c*x))/POW(x,K(12.0)); break; 00338 case 12 : value=(COS(c*x)*POW(c,K(12.0))*POW(x,K(12.0))-K(12.0)*SIN(c*x)*POW(c,K(11.0))*POW(x,K(11.0))-K(132.0)*COS(c*x)*POW(c,K(10.0))*POW(x,K(10.0))+K(1320.0)*SIN(c*x)*c*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x*x+K(11880.0)*COS(c*x)*c*c*c*c*c*c*c*c*x*x*x*x*x*x*x*x-K(95040.0)*SIN(c*x)*c*c*c*c*c*c*c*x*x*x*x*x*x*x-K(665280.0)*COS(c*x)*c*c*c*c*c*c*x*x*x*x*x*x+K(3991680.0)*SIN(c*x)*c*c*c*c*c*x*x*x*x*x+K(19958400.0)*COS(c*x)*c*c*c*c*x*x*x*x-K(79833600.0)*SIN(c*x)*c*c*c*x*x*x-K(239500800.0)*COS(c*x)*c*c*x*x+K(479001600.0)*SIN(c*x)*c*x+K(479001600.0)*COS(c*x))/POW(x,K(13.0)); break; 00339 default : value=K(0.0); 00340 } 00341 value *= POW(sign, (R)(der)); 00342 00343 return value; 00344 } 00345 00346 C kcot(R x, int der, const R *param) /* K(x) = cot(cx) */ 00347 { 00348 R c=param[0]; 00349 R value=K(0.0); 00350 00351 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00352 else switch (der) 00353 { 00354 case 0 : value = K(1.0)/TAN(c * x); break; 00355 case 1 : value = -(K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * c; break; 00356 case 2 : value = K(2.0) * K(1.0)/TAN(c * x) * (K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * c * c; break; 00357 case 3 : value = -K(2.0) * (K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * POW(c, K(3.0)) * (K(1.0) + K(3.0) * POW(K(1.0)/TAN(c * x), K(2.0))); break; 00358 case 4 : value = K(8.0) * (K(1.0) + POW(K(1.0)/TAN(c * x), K(2.0))) * POW(c, K(4.0)) * K(1.0)/TAN(c * x) * (K(2.0) + K(3.0) * POW(K(1.0)/TAN(c * x), K(2.0))); break; 00359 case 5 : value = -K(0.8e1) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.5e1)) * (K(0.15e2) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.15e2) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.2e1)); break; 00360 case 6 : value = K(0.16e2) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.6e1)) * K(1.0)/TAN(c * x) * (K(0.60e2) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.45e2) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.17e2)); break; 00361 case 7 : value = -K(0.16e2) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.7e1)) * (K(0.525e3) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.315e3) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.231e3) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.17e2)); break; 00362 case 8 : value = K(0.128e3) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.8e1)) * K(1.0)/TAN(c * x) * (K(0.630e3) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.315e3) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.378e3) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.62e2)); break; 00363 case 9 : value = -K(0.128e3) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.9e1)) * (K(0.6615e4) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.2835e4) * POW(K(1.0)/TAN(c * x), K(0.8e1)) + K(0.5040e4) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.1320e4) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.62e2)); break; 00364 case 10 : value = K(0.256e3) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.10e2)) * K(1.0)/TAN(c * x) * (K(0.37800e5) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.14175e5) * POW(K(1.0)/TAN(c * x), K(0.8e1)) + K(0.34965e5) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.12720e5) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.1382e4)); break; 00365 case 11 : value = -K(0.256e3) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.11e2)) * (K(0.467775e6) * POW(K(1.0)/TAN(c * x), K(0.8e1)) + K(0.155925e6) * POW(K(1.0)/TAN(c * x), K(0.10e2)) + K(0.509355e6) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.238425e6) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.42306e5) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.1382e4)); break; 00366 case 12 : value = K(0.1024e4) * (K(0.1e1) + POW(K(1.0)/TAN(c * x), K(0.2e1))) * POW(c, K(0.12e2)) * K(1.0)/TAN(c * x) * (K(0.1559250e7) * POW(K(1.0)/TAN(c * x), K(0.8e1)) + K(0.467775e6) * POW(K(1.0)/TAN(c * x), K(0.10e2)) + K(0.1954260e7) * POW(K(1.0)/TAN(c * x), K(0.6e1)) + K(0.1121670e7) * POW(K(1.0)/TAN(c * x), K(0.4e1)) + K(0.280731e6) * POW(K(1.0)/TAN(c * x), K(0.2e1)) + K(0.21844e5)); break; 00367 default : value=K(0.0); 00368 } 00369 00370 return value; 00371 } 00372 00373 00374 C one_over_cube(R x, int der, const R *param) 00375 { 00376 R value=K(0.0); 00377 UNUSED(param); 00378 00379 if (FABS(x)<DBL_EPSILON) value=K(0.0); 00380 else switch (der) 00381 { 00382 case 0 : value = K(1.0)/(x*x*x); break; 00383 case 1 : value = -K(3.0)/(x*x*x*x); break; 00384 case 2 : value = K(12.0)/(x*x*x*x*x); break; 00385 case 3 : value = -K(60.0)/(x*x*x*x*x*x); break; 00386 case 4 : value = K(360.0)/(x*x*x*x*x*x*x); break; 00387 case 5 : value = -K(2520.0)/(x*x*x*x*x*x*x*x); break; 00388 case 6 : value = K(20160.0)/POW(x, K(9.0)); break; 00389 case 7 : value = -K(181440.0)/POW(x, K(10.0)); break; 00390 case 8 : value = K(1814400.0)/POW(x, K(11.0)); break; 00391 case 9 : value = -K(19958400.0)/POW(x, K(12.0)); break; 00392 case 10 : value = K(239500800.0)/POW(x, K(13.0)); break; 00393 case 11 : value = -K(3113510400.0)/POW(x, K(14.0)); break; 00394 case 12 : value = K(43589145600.0)/POW(x, K(15.0)); break; 00395 default : value=K(0.0); 00396 } 00397 00398 return value; 00399 } 00400 00401 00402 /* \} */ 00403 00404 /* kernels.c */