1 /////////////////////////////////////////////////////////////////////////////
2 // einspline: a library for creating and evaluating B-splines //
3 // Copyright (C) 2007 Kenneth P. Esler, Jr. //
4 // //
5 // This program is free software; you can redistribute it and/or modify //
6 // it under the terms of the GNU General Public License as published by //
7 // the Free Software Foundation; either version 2 of the License, or //
8 // (at your option) any later version. //
9 // //
10 // This program is distributed in the hope that it will be useful, //
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of //
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
13 // GNU General Public License for more details. //
14 // //
15 // You should have received a copy of the GNU General Public License //
16 // along with this program; if not, write to the Free Software //
17 // Foundation, Inc., 51 Franklin Street, Fifth Floor, //
18 // Boston, MA 02110-1301 USA //
19 /////////////////////////////////////////////////////////////////////////////
20
21 #include "nugrid.h"
22 #include <cmath>
23 #include <stdlib.h>
24 #include <assert.h>
25
26 #include <boost/math/special_functions/log1p.hpp>
27 #include <boost/math/special_functions/expm1.hpp>
28 using namespace boost::math;
29
30
31 int
center_grid_reverse_map(void * gridptr,double x)32 center_grid_reverse_map (void* gridptr, double x)
33 {
34 center_grid *grid = (center_grid *)gridptr;
35
36 x -= grid->center;
37 double index =
38 copysign (log1p(fabs(x)*grid->aInv)*grid->bInv, x);
39 return (int)floor(grid->half_points + index - grid->even_half);
40 }
41
42 int
log_grid_reverse_map(void * gridptr,double x)43 log_grid_reverse_map (void *gridptr, double x)
44 {
45 log_grid *grid = (log_grid *)gridptr;
46
47 int index = (int) floor(grid->ainv*log(x*grid->startinv));
48
49 if (index < 0)
50 return 0;
51 else
52 return index;
53 }
54
55
56 int
general_grid_reverse_map(void * gridptr,double x)57 general_grid_reverse_map (void* gridptr, double x)
58 {
59 NUgrid* grid = (NUgrid*) gridptr;
60 int N = grid->num_points;
61 double *points = grid->points;
62 if (x <= points[0])
63 return (0);
64 else if (x >= points[N-1])
65 return (N-1);
66 else {
67 int hi = N-1;
68 int lo = 0;
69 bool done = false;
70 while (!done) {
71 int i = (hi+lo)>>1;
72 if (points[i] > x)
73 hi = i;
74 else
75 lo = i;
76 done = (hi-lo)<2;
77 }
78 return (lo);
79 }
80 }
81
82 NUgrid*
create_center_grid(double start,double end,double ratio,int num_points)83 create_center_grid (double start, double end, double ratio,
84 int num_points)
85 {
86 center_grid *grid = new center_grid;
87 if (grid != NULL) {
88 assert (ratio > 1.0);
89 grid->start = start;
90 grid->end = end;
91 grid->center = 0.5*(start + end);
92 grid->num_points = num_points;
93 grid->half_points = num_points/2;
94 grid->odd = ((num_points % 2) == 1);
95 grid->b = log(ratio) / (double)(grid->half_points-1);
96 grid->bInv = 1.0/grid->b;
97 grid->points = new double[num_points];
98 if (grid->odd) {
99 grid->even_half = 0.0;
100 grid->odd_one = 1;
101 grid->a = 0.5*(end-start)/expm1(grid->b*grid->half_points);
102 grid->aInv = 1.0/grid->a;
103 for (int i=-grid->half_points; i<=grid->half_points; i++) {
104 double sign;
105 if (i<0)
106 sign = -1.0;
107 else
108 sign = 1.0;
109 grid->points[i+grid->half_points] =
110 sign * grid->a*expm1(grid->b*abs(i))+grid->center;
111 }
112 }
113 else {
114 grid->even_half = 0.5;
115 grid->odd_one = 0;
116 grid->a =
117 0.5*(end-start)/expm1(grid->b*(-0.5+grid->half_points));
118 grid->aInv = 1.0/grid->a;
119 for (int i=-grid->half_points; i<grid->half_points; i++) {
120 double sign;
121 if (i<0) sign = -1.0;
122 else sign = 1.0;
123 grid->points[i+grid->half_points] =
124 sign * grid->a*expm1(grid->b*fabs(0.5+i)) + grid->center;
125 }
126 }
127 grid->reverse_map = center_grid_reverse_map;
128 grid->code = CENTER;
129 }
130 return (NUgrid*) grid;
131 }
132
133
134 NUgrid*
create_log_grid(double start,double end,int num_points)135 create_log_grid (double start, double end,
136 int num_points)
137 {
138 log_grid *grid = new log_grid;
139 grid->code = LOG;
140 grid->start = start;
141 grid->end = end;
142 grid->num_points = num_points;
143 grid->points = new double[num_points];
144 grid->a = 1.0/(double)(num_points-1)*log(end/start);
145 grid->ainv = 1.0/grid->a;
146 grid->startinv = 1.0/start;
147 for (int i=0; i<num_points; i++)
148 grid->points[i] = start*exp(grid->a*(double)i);
149 grid->reverse_map = log_grid_reverse_map;
150 return (NUgrid*) grid;
151 }
152
153
154 NUgrid*
create_general_grid(double * points,int num_points)155 create_general_grid (double *points, int num_points)
156 {
157 NUgrid* grid = new NUgrid;
158 if (grid != NULL) {
159 grid->reverse_map = general_grid_reverse_map;
160 grid->code = GENERAL;
161 grid->points = new double[num_points];
162 grid->start = points[0];
163 grid->end = points[num_points-1];
164 grid->num_points = num_points;
165 for (int i=0; i<num_points; i++)
166 grid->points[i] = points[i];
167 grid->code = GENERAL;
168 }
169 return grid;
170 }
171
172 void
destroy_grid(NUgrid * grid)173 destroy_grid (NUgrid *grid)
174 {
175 delete[] grid->points;
176 delete grid;
177 }
178