NFFT  3.3.0
damp.c
1 /*
2  * Copyright (c) 2002, 2015 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 /* $Id: util.c 3483 2010-04-23 19:02:34Z keiner $ */
20 
21 #include "infft.h"
22 
27 R Y(modified_fejer)(const INT N, const INT kk)
28 {
29  return (K(2.0) / ((R) (N * N))
30  * (K(1.0) - FABS(K(2.0) * ((R) kk) + K(1.0) ) / ((R) N)));
31 }
32 
34 R Y(modified_jackson2)(const INT N, const INT kk)
35 {
36  INT kj;
37  const R n = (((R) N) / K(2.0) + K(1.0) ) / K(2.0);
38  R result, k;
39 
40  for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
41  {
42  k = (R)(ABS(kj));
43 
44  if (k / n < K(1.0) )
45  result += K(1.0)
46  - (K(3.0) * k + K(6.0) * n * POW(k, K(2.0) )
47  - K(3.0) * POW(k, K(3.0) ))
48  / (K(2.0) * n * (K(2.0) * POW(n, K(2.0) ) + K(1.0) ));
49  else
50  result += (K(2.0) * n - k) * (POW(2 * n - k, K(2.0) ) - K(1.0) )
51  / (K(2.0) * n * (K(2.0) * POW(n, K(2.0) ) + K(1.0) ));
52  }
53 
54  return result;
55 }
56 
58 R Y(modified_jackson4)(const INT N, const INT kk)
59 {
60  INT kj;
61  const R n = (((R) N) / K(2.0) + K(3.0) ) / K(4.0);
62  const R normalisation = (K(2416.0) * POW(n, K(7.0) )
63  + K(1120.0) * POW(n, K(5.0) ) + K(784.0) * POW(n, K(3.0) ) + K(720.0) * n);
64  R result, k;
65 
66  for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
67  {
68  k = (R)(ABS(kj));
69 
70  if (k / n < K(1.0) )
71  result += K(1.0)
72  - (K(1260.0) * k
73  + (K(1680.0) * POW(n, K(5.0) ) + K(2240.0) * POW(n, K(3.0) )
74  + K(2940.0) * n) * POW(k, K(2.0) )
75  - K(1715.0) * POW(k, K(3.0) )
76  - (K(560.0) * POW(n, K(3.0) ) + K(1400.0) * n) * POW(k, K(4.0) )
77  + K(490.0) * POW(k, K(5.0) ) + K(140.0) * n * POW(k, K(6.0) )
78  - K(35.0) * POW(k, K(7.0) )) / normalisation;
79 
80  if ((K(1.0) <= k / n) && (k / n < K(2.0) ))
81  result += ((K(2472.0) * POW(n, K(7.0) ) + K(336.0) * POW(n, K(5.0) )
82  + K(3528.0) * POW(n, K(3.0) ) - K(1296.0) * n)
83  - (K(392.0) * POW(n, K(6.0) ) - K(3920.0) * POW(n, K(4.0) )
84  + K(8232.0) * POW(n, K(2.0) ) - K(756.0) )*k
85  - (K(504.0)*POW(n, K(5.0)) + K(10080.0)*POW(n, K(3.0))
86  - K(5292.0)*n)*POW(k, K(2.0)) - (K(1960.0)*POW(n, K(4.0))
87  - K(7840.0)*POW(n, K(2.0)) + K(1029.0))*POW(k, K(3.0))
88  + (K(2520.0)*POW(n, K(3.0)) - K(2520.0)*n) * POW(k, K(4.0))
89  - (K(1176.0)*POW(n, K(2.0)) - K(294.0)) * POW(k, K(5.0))
90  + K(252.0)*n*POW(k, K(6.0)) - K(21.0)*POW(k, K(7.0)))/normalisation;
91 
92  if ((K(2.0) <= k / n) && (k / n < K(3.0) ))
93  result += (-(K(1112.0) * POW(n, K(7.0) ) - K(12880.0) * POW(n, K(5.0) )
94  + K(7448.0) * POW(n, K(3.0) ) - K(720.0) * n)
95  + (K(12152.0) * POW(n, K(6.0) ) - K(27440.0) * POW(n, K(4.0) )
96  + K(8232.0) * POW(n, K(2.0) ) - K(252.0) )*k
97  - (K(19320.0)*POW(n, K(5.0)) - K(21280.0)*POW(n, K(3.0))
98  + K(2940.0)*n)*POW(k, K(2.0)) + (K(13720.0)*POW(n, K(4.0))
99  - K(7840.0)*POW(n, K(2.0)) + K(343.0))*POW(k, K(3.0))
100  - (K(5320.0)*POW(n, K(3.0)) - K(1400.0)*n)*POW(k, K(4.0))
101  + (K(1176.0)*POW(n, K(2.0)) - K(98.0))*POW(k, K(5.0))
102  - K(140.0)*n*POW(k, K(6.0)) + K(7.0) * POW(k, K(7.0)))/normalisation;
103 
104  if ((K(3.0) <= k / n) && (k / n < K(4.0) ))
105  result += ((4 * n - k)
106  * (POW(4 * n - k, K(2.0) ) - K(1.0) )*(POW(4*n-k, K(2.0))
107  - K(4.0))*(POW(4*n-k, K(2.0)) - K(9.0)))/normalisation;
108  }
109 
110  return result;
111 }
112 
114 R Y(modified_sobolev)(const R mu, const INT kk)
115 {
116  R result;
117  INT kj, k;
118 
119  for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
120  {
121  k = ABS(kj);
122  if (k == 0)
123  result += K(1.0);
124  else
125  result += POW((R) k, -K(2.0) * mu);
126  }
127 
128  return result;
129 }
130 
132 R Y(modified_multiquadric)(const R mu, const R c, const INT kk)
133 {
134  R result;
135  INT kj, k;
136 
137  for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
138  {
139  k = ABS(kj);
140  result += POW((R)(k * k) + c * c, -mu);
141  }
142 
143  return result;
144 }