Mercurial > hg > octave-lyh
annotate libinterp/corefcn/__contourc__.cc @ 15220:61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
* __contourc__.cc, __qp__.cc, quadcc.cc, qz.cc, eigs.cc, mex.cc,
graphics.cc, CMatrix.cc, CollocWt.cc, DASPK-opts.in, DASRT-opts.in,
DASSL-opts.in, LSODE-opts.in, Quad-opts.in, Range.cc, dMatrix.cc,
eigs-base.cc, fCMatrix.cc, fMatrix.cc: Replace DBL_EPSILON and
FLT_EPSILON with std::numeric_limits<T>::epsilon.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 23 Aug 2012 12:13:59 -0400 |
| parents | 2fc554ffbc28 |
| children | 704e15f8fecd |
| rev | line source |
|---|---|
| 6257 | 1 /* Contour lines for function evaluated on a grid. |
| 2 | |
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
12805
diff
changeset
|
3 Copyright (C) 2007-2012 Kai Habel |
| 7017 | 4 Copyright (C) 2004, 2007 Shai Ayal |
| 6257 | 5 |
| 6 Adapted to an oct file from the stand alone contourl by Victro Munoz | |
| 7 Copyright (C) 2004 Victor Munoz | |
| 8 | |
| 9 Based on contour plot routine (plcont.c) in PLPlot package | |
| 10 http://plplot.org/ | |
| 11 | |
| 12 Copyright (C) 1995, 2000, 2001 Maurice LeBrun | |
| 13 Copyright (C) 2000, 2002 Joao Cardoso | |
| 14 Copyright (C) 2000, 2001, 2002, 2004 Alan W. Irwin | |
| 15 Copyright (C) 2004 Andrew Ross | |
| 16 | |
| 17 This file is part of Octave. | |
| 18 | |
| 19 Octave is free software; you can redistribute it and/or modify it | |
| 20 under the terms of the GNU General Public License as published by the | |
| 7016 | 21 Free Software Foundation; either version 3 of the License, or (at your |
| 22 option) any later version. | |
| 6257 | 23 |
| 24 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 25 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 26 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 27 for more details. | |
| 28 | |
| 29 You should have received a copy of the GNU General Public License | |
| 7016 | 30 along with Octave; see the file COPYING. If not, see |
| 31 <http://www.gnu.org/licenses/>. | |
| 6257 | 32 |
| 33 */ | |
| 34 | |
| 35 #ifdef HAVE_CONFIG_H | |
| 36 #include <config.h> | |
| 37 #endif | |
| 38 | |
| 7042 | 39 #include <cfloat> |
| 40 | |
| 6257 | 41 #include "quit.h" |
| 42 | |
|
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14501
diff
changeset
|
43 #include "defun.h" |
| 6257 | 44 #include "error.h" |
| 45 #include "oct-obj.h" | |
| 46 | |
| 47 static Matrix this_contour; | |
| 48 static Matrix contourc; | |
| 49 static int elem; | |
| 50 | |
| 7042 | 51 // This is the quanta in which we increase this_contour. |
| 6257 | 52 #define CONTOUR_QUANT 50 |
| 53 | |
| 7042 | 54 // Add a coordinate point (x,y) to this_contour. |
| 6257 | 55 |
| 56 static void | |
| 7042 | 57 add_point (double x, double y) |
| 6257 | 58 { |
| 59 if (elem % CONTOUR_QUANT == 0) | |
| 60 this_contour = this_contour.append (Matrix (2, CONTOUR_QUANT, 0)); | |
| 61 | |
| 62 this_contour (0, elem) = x; | |
| 63 this_contour (1, elem) = y; | |
| 64 elem++; | |
| 65 } | |
| 66 | |
| 7042 | 67 // Add contents of current contour to contourc. |
| 6257 | 68 // this_contour.cols () - 1; |
| 69 | |
| 70 static void | |
| 7042 | 71 end_contour (void) |
| 6257 | 72 { |
| 73 if (elem > 2) | |
| 74 { | |
| 75 this_contour (1, 0) = elem - 1; | |
| 76 contourc = contourc.append (this_contour.extract_n (0, 0, 2, elem)); | |
| 77 } | |
| 7042 | 78 |
| 6257 | 79 this_contour = Matrix (); |
| 80 elem = 0; | |
| 81 } | |
| 82 | |
| 7042 | 83 // Start a new contour, and add contents of current one to contourc. |
| 6257 | 84 |
| 85 static void | |
| 7042 | 86 start_contour (double lvl, double x, double y) |
| 6257 | 87 { |
| 7042 | 88 end_contour (); |
| 6257 | 89 this_contour.resize (2, 0); |
| 7042 | 90 add_point (lvl, 0); |
| 91 add_point (x, y); | |
| 6257 | 92 } |
| 93 | |
| 94 static void | |
| 7042 | 95 drawcn (const RowVector& X, const RowVector& Y, const Matrix& Z, |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
96 double lvl, int r, int c, double ct_x, double ct_y, |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
97 unsigned int start_edge, bool first, charMatrix& mark) |
| 6257 | 98 { |
| 7042 | 99 double px[4], py[4], pz[4], tmp; |
| 7070 | 100 unsigned int stop_edge, next_edge, pt[2]; |
| 7042 | 101 int next_r, next_c; |
| 6257 | 102 |
| 7042 | 103 //get x, y, and z - lvl for current facet |
| 104 px[0] = px[3] = X(c); | |
| 105 px[1] = px[2] = X(c+1); | |
| 6257 | 106 |
| 7042 | 107 py[0] = py[1] = Y(r); |
| 108 py[2] = py[3] = Y(r+1); | |
| 6257 | 109 |
| 7042 | 110 pz[3] = Z(r+1, c) - lvl; |
| 111 pz[2] = Z(r+1, c + 1) - lvl; | |
| 112 pz[1] = Z(r, c+1) - lvl; | |
| 113 pz[0] = Z(r, c) - lvl; | |
| 6257 | 114 |
| 7042 | 115 // Facet edge and point naming assignment. |
| 116 // | |
| 117 // 0-----1 .-0-. | |
| 118 // | | | | | |
| 119 // | | 3 1 | |
| 120 // | | | | | |
| 121 // 3-----2 .-2-. | |
| 6257 | 122 |
| 7042 | 123 // Get mark value of current facet. |
| 124 char id = static_cast<char> (mark(r, c)); | |
| 125 | |
| 126 // Check startedge s. | |
| 127 if (start_edge == 255) | |
| 6257 | 128 { |
| 7042 | 129 // Find start edge. |
| 7070 | 130 for (unsigned int k = 0; k < 4; k++) |
| 131 if (static_cast<char> (1 << k) & id) | |
| 7042 | 132 start_edge = k; |
| 6257 | 133 } |
| 134 | |
| 7042 | 135 if (start_edge == 255) |
| 136 return; | |
| 137 | |
| 138 // Decrease mark value of current facet for start edge. | |
| 7070 | 139 mark(r, c) -= static_cast<char> (1 << start_edge); |
| 6257 | 140 |
| 7042 | 141 // Next point (clockwise). |
| 142 pt[0] = start_edge; | |
| 143 pt[1] = (pt[0] + 1) % 4; | |
| 144 | |
| 145 // Calculate contour segment start if first of contour. | |
| 146 if (first) | |
| 147 { | |
| 148 tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]); | |
| 149 | |
| 150 if (xisnan (tmp)) | |
| 151 ct_x = ct_y = 0.5; | |
| 152 else | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
153 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
154 ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
155 ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
156 } |
| 7042 | 157 |
| 158 start_contour (lvl, ct_x, ct_y); | |
| 159 } | |
| 160 | |
| 8333 | 161 // Find stop edge. |
| 162 // FIXME -- perhaps this should use a while loop? | |
| 7070 | 163 for (unsigned int k = 1; k <= 4; k++) |
| 7042 | 164 { |
| 165 if (start_edge == 0 || start_edge == 2) | |
| 166 stop_edge = (start_edge + k) % 4; | |
| 167 else | |
| 168 stop_edge = (start_edge - k) % 4; | |
| 169 | |
| 7070 | 170 if (static_cast<char> (1 << stop_edge) & id) |
| 7042 | 171 break; |
| 172 } | |
| 173 | |
| 174 pt[0] = stop_edge; | |
| 175 pt[1] = (pt[0] + 1) % 4; | |
| 176 tmp = fabs (pz[pt[1]]) / fabs (pz[pt[0]]); | |
| 6257 | 177 |
| 7042 | 178 if (xisnan (tmp)) |
| 179 ct_x = ct_y = 0.5; | |
| 180 else | |
| 181 { | |
| 182 ct_x = px[pt[0]] + (px[pt[1]] - px[pt[0]])/(1 + tmp); | |
| 183 ct_y = py[pt[0]] + (py[pt[1]] - py[pt[0]])/(1 + tmp); | |
| 184 } | |
| 185 | |
| 186 // Add point to contour. | |
| 187 add_point (ct_x, ct_y); | |
| 188 | |
| 189 // Decrease id value of current facet for start edge. | |
| 7070 | 190 mark(r, c) -= static_cast<char> (1 << stop_edge); |
| 7042 | 191 |
| 192 // Find next facet. | |
| 193 next_c = c; | |
| 194 next_r = r; | |
| 195 | |
| 196 if (stop_edge == 0) | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
197 next_r--; |
| 7042 | 198 else if (stop_edge == 1) |
| 199 next_c++; | |
| 200 else if (stop_edge == 2) | |
| 201 next_r++; | |
| 202 else if (stop_edge == 3) | |
| 203 next_c--; | |
| 204 | |
| 205 // Check if next facet is not done yet. | |
| 206 // Go to next facet. | |
| 207 if (next_r >= 0 && next_c >= 0 && next_r < mark.rows () | |
| 208 && next_c < mark.cols () && mark(next_r, next_c) > 0) | |
| 209 { | |
| 210 next_edge = (stop_edge + 2) % 4; | |
| 211 drawcn (X, Y, Z, lvl, next_r, next_c, ct_x, ct_y, next_edge, false, mark); | |
| 6257 | 212 } |
| 213 } | |
| 214 | |
| 215 static void | |
| 7042 | 216 mark_facets (const Matrix& Z, charMatrix& mark, double lvl) |
| 6257 | 217 { |
| 7070 | 218 unsigned int nr = mark.rows (); |
| 219 unsigned int nc = mark.cols (); | |
| 7042 | 220 |
| 221 double f[4]; | |
| 222 | |
| 7070 | 223 for (unsigned int c = 0; c < nc; c++) |
| 224 for (unsigned int r = 0; r < nr; r++) | |
| 7042 | 225 { |
| 226 f[0] = Z(r, c) - lvl; | |
| 227 f[1] = Z(r, c+1) - lvl; | |
| 228 f[3] = Z(r+1, c) - lvl; | |
| 229 f[2] = Z(r+1, c+1) - lvl; | |
| 230 | |
| 7070 | 231 for (unsigned int i = 0; i < 4; i++) |
|
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
232 if (fabs(f[i]) < std::numeric_limits<double>::epsilon ()) |
|
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
233 f[i] = std::numeric_limits<double>::epsilon (); |
| 7042 | 234 |
| 235 if (f[1] * f[2] < 0) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
236 mark(r, c) += 2; |
| 7042 | 237 |
| 238 if (f[0] * f[3] < 0) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
239 mark(r, c) += 8; |
| 7042 | 240 } |
| 6257 | 241 |
| 7070 | 242 for (unsigned int r = 0; r < nr; r++) |
| 243 for (unsigned int c = 0; c < nc; c++) | |
| 7042 | 244 { |
| 245 f[0] = Z(r, c) - lvl; | |
| 246 f[1] = Z(r, c+1) - lvl; | |
| 247 f[3] = Z(r+1, c) - lvl; | |
| 248 f[2] = Z(r+1, c+1) - lvl; | |
| 249 | |
| 7070 | 250 for (unsigned int i = 0; i < 4; i++) |
|
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
251 if (fabs(f[i]) < std::numeric_limits<double>::epsilon ()) |
|
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
252 f[i] = std::numeric_limits<double>::epsilon (); |
| 7042 | 253 |
| 254 if (f[0] * f[1] < 0) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
255 mark(r, c) += 1; |
| 7042 | 256 |
| 257 if (f[2] * f[3] < 0) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
258 mark(r, c) += 4; |
| 7042 | 259 } |
| 260 } | |
| 261 | |
| 262 static void | |
| 263 cntr (const RowVector& X, const RowVector& Y, const Matrix& Z, double lvl) | |
| 264 { | |
| 7070 | 265 unsigned int nr = Z.rows (); |
| 266 unsigned int nc = Z.cols (); | |
| 7042 | 267 |
| 268 charMatrix mark (nr - 1, nc - 1, 0); | |
| 269 | |
| 270 mark_facets (Z, mark, lvl); | |
| 271 | |
| 272 // Find contours that start at a domain edge. | |
| 273 | |
| 7070 | 274 for (unsigned int c = 0; c < nc - 1; c++) |
| 6257 | 275 { |
| 7042 | 276 // Top. |
| 277 if (mark(0, c) & 1) | |
| 278 drawcn (X, Y, Z, lvl, 0, c, 0.0, 0.0, 0, true, mark); | |
| 279 | |
| 280 // Bottom. | |
| 281 if (mark(nr - 2, c) & 4) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
282 drawcn (X, Y, Z, lvl, nr - 2, c, 0.0, 0.0, 2, true, mark); |
| 6257 | 283 } |
| 7042 | 284 |
| 7070 | 285 for (unsigned int r = 0; r < nr - 1; r++) |
| 7042 | 286 { |
| 287 // Left. | |
| 288 if (mark(r, 0) & 8) | |
| 289 drawcn (X, Y, Z, lvl, r, 0, 0.0, 0.0, 3, true, mark); | |
| 290 | |
| 291 // Right. | |
| 292 if (mark(r, nc - 2) & 2) | |
| 293 drawcn (X, Y, Z, lvl, r, nc - 2, 0.0, 0.0, 1, true, mark); | |
| 294 } | |
| 295 | |
| 7070 | 296 for (unsigned int r = 0; r < nr - 1; r++) |
| 297 for (unsigned int c = 0; c < nc - 1; c++) | |
| 7042 | 298 if (mark (r, c) > 0) |
| 299 drawcn (X, Y, Z, lvl, r, c, 0.0, 0.0, 255, true, mark); | |
| 6257 | 300 } |
| 301 | |
|
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14501
diff
changeset
|
302 DEFUN (__contourc__, args, , |
| 6257 | 303 "-*- texinfo -*-\n\ |
|
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14501
diff
changeset
|
304 @deftypefn {Built-in Function} {} __contourc__ (@var{x}, @var{y}, @var{z}, @var{levels})\n\ |
| 6945 | 305 Undocumented internal function.\n\ |
| 6257 | 306 @end deftypefn") |
| 307 { | |
| 308 octave_value retval; | |
| 309 | |
| 310 if (args.length () == 4) | |
| 311 { | |
| 312 RowVector X = args (0).row_vector_value (); | |
| 313 RowVector Y = args (1).row_vector_value (); | |
| 7042 | 314 Matrix Z = args (2).matrix_value (); |
| 6257 | 315 RowVector L = args (3).row_vector_value (); |
| 316 | |
| 317 if (! error_state) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
318 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
319 contourc.resize (2, 0); |
| 6257 | 320 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
321 for (int i = 0; i < L.length (); i++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
322 cntr (X, Y, Z, L (i)); |
| 6257 | 323 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
324 end_contour (); |
| 6257 | 325 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
326 retval = contourc; |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
327 } |
| 6257 | 328 else |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
329 error ("__contourc__: invalid argument values"); |
| 6257 | 330 } |
| 331 else | |
| 332 print_usage (); | |
| 333 | |
| 334 return retval; | |
| 335 } | |
|
12805
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
336 |
|
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
337 /* |
|
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
338 ## No test needed for internal helper function. |
|
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
339 %!assert (1) |
|
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
340 */ |