NFFT 3.5.3alpha
damp.c
1/*
2 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
3 *
4 * This program is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU General Public License as published by the Free Software
6 * Foundation; either version 2 of the License, or (at your option) any later
7 * version.
8 *
9 * This program is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU General Public License along with
15 * this program; if not, write to the Free Software Foundation, Inc., 51
16 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18
19#include "infft.h"
20
25R Y(modified_fejer)(const INT N, const INT kk)
26{
27 return (K(2.0) / ((R) (N * N))
28 * (K(1.0) - FABS(K(2.0) * ((R) kk) + K(1.0) ) / ((R) N)));
29}
30
32R Y(modified_jackson2)(const INT N, const INT kk)
33{
34 INT kj;
35 const R n = (((R) N) / K(2.0) + K(1.0) ) / K(2.0);
36 R result, k;
37
38 for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
39 {
40 k = (R)(ABS(kj));
41
42 if (k / n < K(1.0) )
43 result += K(1.0)
44 - (K(3.0) * k + K(6.0) * n * POW(k, K(2.0) )
45 - K(3.0) * POW(k, K(3.0) ))
46 / (K(2.0) * n * (K(2.0) * POW(n, K(2.0) ) + K(1.0) ));
47 else
48 result += (K(2.0) * n - k) * (POW(2 * n - k, K(2.0) ) - K(1.0) )
49 / (K(2.0) * n * (K(2.0) * POW(n, K(2.0) ) + K(1.0) ));
50 }
51
52 return result;
53}
54
56R Y(modified_jackson4)(const INT N, const INT kk)
57{
58 INT kj;
59 const R n = (((R) N) / K(2.0) + K(3.0) ) / K(4.0);
60 const R normalisation = (K(2416.0) * POW(n, K(7.0) )
61 + K(1120.0) * POW(n, K(5.0) ) + K(784.0) * POW(n, K(3.0) ) + K(720.0) * n);
62 R result, k;
63
64 for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
65 {
66 k = (R)(ABS(kj));
67
68 if (k / n < K(1.0) )
69 result += K(1.0)
70 - (K(1260.0) * k
71 + (K(1680.0) * POW(n, K(5.0) ) + K(2240.0) * POW(n, K(3.0) )
72 + K(2940.0) * n) * POW(k, K(2.0) )
73 - K(1715.0) * POW(k, K(3.0) )
74 - (K(560.0) * POW(n, K(3.0) ) + K(1400.0) * n) * POW(k, K(4.0) )
75 + K(490.0) * POW(k, K(5.0) ) + K(140.0) * n * POW(k, K(6.0) )
76 - K(35.0) * POW(k, K(7.0) )) / normalisation;
77
78 if ((K(1.0) <= k / n) && (k / n < K(2.0) ))
79 result += ((K(2472.0) * POW(n, K(7.0) ) + K(336.0) * POW(n, K(5.0) )
80 + K(3528.0) * POW(n, K(3.0) ) - K(1296.0) * n)
81 - (K(392.0) * POW(n, K(6.0) ) - K(3920.0) * POW(n, K(4.0) )
82 + K(8232.0) * POW(n, K(2.0) ) - K(756.0) )*k
83 - (K(504.0)*POW(n, K(5.0)) + K(10080.0)*POW(n, K(3.0))
84 - K(5292.0)*n)*POW(k, K(2.0)) - (K(1960.0)*POW(n, K(4.0))
85 - K(7840.0)*POW(n, K(2.0)) + K(1029.0))*POW(k, K(3.0))
86 + (K(2520.0)*POW(n, K(3.0)) - K(2520.0)*n) * POW(k, K(4.0))
87 - (K(1176.0)*POW(n, K(2.0)) - K(294.0)) * POW(k, K(5.0))
88 + K(252.0)*n*POW(k, K(6.0)) - K(21.0)*POW(k, K(7.0)))/normalisation;
89
90 if ((K(2.0) <= k / n) && (k / n < K(3.0) ))
91 result += (-(K(1112.0) * POW(n, K(7.0) ) - K(12880.0) * POW(n, K(5.0) )
92 + K(7448.0) * POW(n, K(3.0) ) - K(720.0) * n)
93 + (K(12152.0) * POW(n, K(6.0) ) - K(27440.0) * POW(n, K(4.0) )
94 + K(8232.0) * POW(n, K(2.0) ) - K(252.0) )*k
95 - (K(19320.0)*POW(n, K(5.0)) - K(21280.0)*POW(n, K(3.0))
96 + K(2940.0)*n)*POW(k, K(2.0)) + (K(13720.0)*POW(n, K(4.0))
97 - K(7840.0)*POW(n, K(2.0)) + K(343.0))*POW(k, K(3.0))
98 - (K(5320.0)*POW(n, K(3.0)) - K(1400.0)*n)*POW(k, K(4.0))
99 + (K(1176.0)*POW(n, K(2.0)) - K(98.0))*POW(k, K(5.0))
100 - K(140.0)*n*POW(k, K(6.0)) + K(7.0) * POW(k, K(7.0)))/normalisation;
101
102 if ((K(3.0) <= k / n) && (k / n < K(4.0) ))
103 result += ((4 * n - k)
104 * (POW(4 * n - k, K(2.0) ) - K(1.0) )*(POW(4*n-k, K(2.0))
105 - K(4.0))*(POW(4*n-k, K(2.0)) - K(9.0)))/normalisation;
106 }
107
108 return result;
109}
110
112R Y(modified_sobolev)(const R mu, const INT kk)
113{
114 R result;
115 INT kj, k;
116
117 for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
118 {
119 k = ABS(kj);
120 if (k == 0)
121 result += K(1.0);
122 else
123 result += POW((R) k, -K(2.0) * mu);
124 }
125
126 return result;
127}
128
130R Y(modified_multiquadric)(const R mu, const R c, const INT kk)
131{
132 R result;
133 INT kj, k;
134
135 for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
136 {
137 k = ABS(kj);
138 result += POW((R)(k * k) + c * c, -mu);
139 }
140
141 return result;
142}
Internal header file for auxiliary definitions and functions.