File: interp.cc

package info (click to toggle)
gri 2.4.2-1
  • links: PTS
  • area: main
  • in suites: potato
  • size: 4,540 kB
  • ctags: 1,966
  • sloc: cpp: 32,542; lisp: 3,243; perl: 806; makefile: 548; sh: 253
file content (127 lines) | stat: -rw-r--r-- 3,570 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <stddef.h>
#include	"extern.hh"
#include        "files.hh"
#include        "superus.hh"

// `interpolate x grid to .left. .right. .inc.|{/.cols.}'
// `interpolate y grid to .bottom. .top. .inc.|{/.rows.}'
bool
interpolateCmd()
{
    bool is_x = false;
    if (word_is(1, "x")) 
	is_x = true;
    else if (!word_is(1, "y")) {
	err("Second word must be `x' or `y'");
    }
    int num;
    double start, stop, increment;
    if (is_x) {
	Require(getdnum(_word[4], &start), READ_WORD_ERROR(".left."));
	Require(getdnum(_word[5], &stop),  READ_WORD_ERROR(".right."));
    } else {
	Require(getdnum(_word[4], &start), READ_WORD_ERROR(".bottom."));
	Require(getdnum(_word[5], &stop),  READ_WORD_ERROR(".top."));
    }
    switch (_nword) {
    case 7:
	if (is_x) {
	    Require(getdnum(_word[6], &increment), 
		    READ_WORD_ERROR(".cols."));
	} else {
	    Require(getdnum(_word[6], &increment), 
		    READ_WORD_ERROR(".rows."));
	}
	num = int(1 + (stop - start) / increment);
	break;
    case 8:
	Require(word_is(6, "/"), err("6-th word must be `/'"));
	Require(getinum(_word[7], &num), READ_WORD_ERROR(".inc."));
	if (num < 2) {
	    err("Increment .inc. must exceed 1");
	    return false;
	}
	increment = (stop - start) / (num - 1);
	break;
    }
    if (num < 0) {
	err("Sign of increment disagrees with start and stop values");
	return false;
    }
    //
    // The grid must already exist
    if (!_grid_exists) {
	err("Cannot `convert grid to image' since no grid data exist yet");
	return false;
    }
    if (!_xgrid_exists) {
	err("Cannot `convert grid to image' since x-grid doesn't exist yet");
	return false;
    }
    if (!_ygrid_exists) {
	err("Cannot `convert grid to image' since y-grid doesn't exist yet");
	return false;
    }
    increment *= 0.9999999;	// to avoid overflow
    double znew;
    GriMatrix<double> _f_xy_new;
    if (is_x) {
	vector<double> ygrid((size_t)_num_ymatrix_data, 0.0);
	_f_xy_new.set_size(num, _num_ymatrix_data);
	unsigned int i, j;
	j = _num_ymatrix_data - 1;
	do {
	    ygrid[j] = _ymatrix[j];
	    double xnew = start;
	    for (i = 0; i < (unsigned int)num; i++) {
		grid_interp(xnew, _ymatrix[j], &znew);
		_f_xy_new(i, j) = znew;
		xnew += increment;
	    }
	} while (j-- != 0);
	// Dump into official storage
	allocate_grid_storage(num, _num_ymatrix_data);
	allocate_xmatrix_storage(num);
	for (i = 0; i < (unsigned int)num; i++) 
	    _xmatrix[i] = start + i * increment;
	for (j = 0; j < _num_ymatrix_data; j++)
	    _ymatrix[j] = ygrid[j];
	for (j = 0; j < _num_ymatrix_data; j++) {
	    for (i = 0; i < (unsigned int)num; i++) {
		_f_xy(i, j) = _f_xy_new(i, j);
		_legit_xy(i, j) = true;
	    }
	}
    } else {
	vector<double> xgrid((size_t)_num_xmatrix_data, 0.0);
	_f_xy_new.set_size(_num_xmatrix_data, num);
	unsigned int i, j;
	for (i = 0; i < _num_xmatrix_data; i++) {
	    xgrid[i] = _xmatrix[i];
	    double ynew = start;
	    for (j = 0; j < (unsigned int)num; j++) {
		grid_interp(_xmatrix[i], ynew, &znew);
		_f_xy_new(i, j) = znew;
		ynew += increment;
	    }
	}
	// Dump into official storage
	allocate_grid_storage(_num_xmatrix_data, num);
	allocate_ymatrix_storage(num);
	for (i = 0; i < _num_xmatrix_data; i++) 
	    _xmatrix[i] = xgrid[i];
	for (j = 0; j < (unsigned int) num; j++)
	    _ymatrix[j] = start + j * increment;
	for (i = 0; i < _num_xmatrix_data; i++) {
	    for (j = 0; j < (unsigned int)num; j++) {
		_f_xy(i, j) = _f_xy_new(i, j);
		_legit_xy(i, j) = true;
	    }
	}
    }
    return true;
}