VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
Matrix.cpp
Go to the documentation of this file.
1 /**
2 * \file Matrix.cpp
3 *
4 * \author Developed under supervision of the PI Yuri Shprits
5 */
6 #ifndef matrix_array_MATRIX_CPP
7 #define matrix_array_MATRIX_CPP
8 
9 #include "Matrix.h"
10 
11 using namespace std;
12 
13 #if defined(_WIN32) || defined(_WIN64)
14 #define strncasecmp _strnicmp
15 #define strcasecmp _stricmp
16 #endif
17 
18 // Memory related functions
19 
20 /// Allocating memory for 1D matrix
21 /// \todo Move to Matrix.cpp
22 template<class T>inline T* matrix(long Rows)
23 {
24  T *m=new T[Rows];
25  // assert(m!=NULL);
26  if (m == NULL) {
27  throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", Rows*sizeof(T));
28  }
29  return m;
30 }
31 
32 /// Initilizing memory for 2D matrix
33 /// \todo Move to Matrix.cpp
34 template<class T>inline T** matrix(long Rows, long Columns)
35 {
36  // allocating memory for array of pinters
37  T **m=new T*[Rows];
38  // assert(m!=NULL);
39  if (m == NULL) {
40  throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", Rows * sizeof(T));
41  }
42  // allocating memory for data array
43  m[0] = new T[Rows * Columns];
44  // assert(m[0]!=NULL);
45  if (m[0] == NULL) {
46  throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", Rows * Columns * sizeof(T));
47  }
48  // assign pointers to data ranges
49  for(long i=1; i<Rows; i++) m[i] = m[i-1] + Columns;
50  return m;
51 }
52 
53 /// Initializing memory for 3D matrix
54 /// \todo Move to Matrix.cpp
55 template<class T>inline T*** matrix(int size_x, int size_y, int size_z)
56 {
57  // allocating memory for array of pointers to pointers
58  T ***m=new T**[size_x];
59  // assert(m!=NULL);
60  if (m == NULL) {
61  throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", size_x * sizeof(T));
62  }
63  for (int x = 0; x < size_x; x++) {
64  // for each pointer allocating memory for array of pointers
65  m[x] = new T*[size_y];
66  // assert(m[x]!=NULL);
67  if (m[x] == NULL) {
68  throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", size_y * sizeof(T));
69  }
70  }
71  //allocating memory for data array
72  m[0][0] = new T[size_x * size_y * size_z];
73  // assert(m[0][0]!=NULL);
74  if (m[0][0] == NULL) {
75  throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", size_x * size_y * size_z * sizeof(T));
76  }
77  for (int x = 0; x < size_x; x++) {
78  for (int y = 0; y < size_y; y++) {
79  // assign pointers to data ranges
80  m[x][y] = m[0][0] + (x*size_y + y)*size_z;
81  }
82  }
83 
84  return m;
85 }
86 
87 /// Freeing memory for 1D matrix
88 /// \todo Move to Matrix.cpp
89 template<class T>inline void free_matrix(T* m) {
90  delete m;
91 }
92 
93 /// Freeing memory for 2D matrix
94 /// \todo Move to Matrix.cpp
95 template<class T>inline void free_matrix(T** m) {
96  delete[](m[0]);
97  delete[](m);
98 }
99 
100 /// Freeing memory for 3D matrix
101 /// \todo Move to Matrix.cpp
102 template<class T>inline void free_matrix(T*** m, int size_x, int size_y) {
103  delete[](m[0][0]);
104  for (int x = 0; x < size_x; x++) {
105  delete[](m[x]);
106  }
107  delete[](m);
108 }
109 
110 //////////////////////////////////////////////////
111 //
112 // Matrix 1D
113 //
114 //////////////////////////////////////////////////
115 
116 /**
117 * Constructor.
118 *
119 * Runs allocating memory function and store matrix name.
120 *
121 * \param x_size - size of the matrix
122 * \param name - name of the matrix
123 *
124 */
125 template<class T>
126 Matrix1D<T>::Matrix1D( int x_size , string name) {
127  initialized = false;
128  this->name = name;
129  // allocating memory
130  AllocateMemory(x_size);
131 }
132 
133 /**
134 * Constructor.
135 *
136 * Runs allocating memory function.
137 *
138 * \param x_size - size of the matrix
139 *
140 */
141 template<class T>
142 Matrix1D<T>::Matrix1D( int x_size ) {
143  initialized = false;
144  // allocating memory
145  AllocateMemory(x_size);
146 }
147 
148 /**
149 * Constructor.
150 *
151 * Make new matrix equal to Matrix M.
152 *
153 * \param &M - matrix M
154 *
155 */
156 template<class T>
158  initialized = false;
159  this->operator = (M);
160 }
161 
162 /**
163 * Destructor. Destruct the class.
164 */
165 template<class T>
167  if (initialized) free_matrix<T>(matrix_array);
168 }
169 
170 /**
171 * Allocating memory
172 *
173 * \param x_size - size x
174 */
175 template<class T>
176 void Matrix1D<T>::AllocateMemory( int x_size ) {
177  this->size_x = x_size;
178  // using inline template for memory allocation
179  matrix_array = matrix<T>(x_size);
180  initialized = true;
181  // for (int i = 0; i < size_x; i++) matrix_array[i] = 2;
182 }
183 
184 
185 /**
186 * Operator [i], returns value of element i.
187 * If DEBUG_MODE defined, check if matrix has been initialized.
188 *
189 * \param i - number of element to return
190 */
191 template<class T>
192 inline T& Matrix1D<T>::operator[](int i) {
193 #ifdef DEBUG_MODE
194  if (!initialized) {
195  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
196  }
197 #endif
198  // return i-value
199  return matrix_array[i];
200 }
201 
202 /**
203 * Operator [i], returns value of element i, version returns 'const' value, can not be later modified.
204 * If DEBUG_MODE defined, check if matrix has been initialized.
205 *
206 * \param i - number of element to return
207 */
208 template<class T>
209 inline T& Matrix1D<T>::operator[](int i) const {
210 #ifdef DEBUG_MODE
211  if (!initialized) {
212  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
213  }
214 #endif
215  // return i-value
216  return matrix_array[i];
217 }
218 
219 /**
220 * Operator (x), returns value of element x.
221 * If DEBUG_MODE defined, check if matrix has been initialized.
222 * No dofference between [] and () operators for 1d-matrix class.
223 *
224 * \param x - number of element to return
225 */
226 template<class T>
227 inline T& Matrix1D<T>::operator()(int x) {
228 #ifdef DEBUG_MODE
229  if (!initialized) {
230  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
231  }
232 #endif
233  // return x-value
234  return matrix_array[x];
235 }
236 
237 /**
238 * Make matrix equal to matrix M
239 * Return the same instance of class Matrix.
240 *
241 * \param &M - matrix M
242 */
243 template<class T>
245  // Check, if both matrix aren't the same one
246  if (this!=&M) {
247  // check, if RHS matrix was initialized
248  if (M.initialized) {
249  // free LHS matrix, if it was initialized
250  if (initialized) free_matrix<T>(matrix_array);
251  // copy parameters
252  this->size_x = M.size_x;
253  this->name = M.name;
254  // allocating memory for LHS matrix
255  this->AllocateMemory(this->size_x);
256  // fast values copying as a memory range
257  memcpy( matrix_array, M.matrix_array, this->size_x * sizeof( T ) );
258  } else {
259  this->initialized = false;
260  }
261  }
262  return *this;
263 }
264 
265 /**
266 * Make matrix equal to value Val.
267 * Return the same instance of class Matrix.
268 *
269 * \param Val - value val
270 */
271 template<class T>
273  for (int x = 0; x < this->size_x; x++)
274  matrix_array[x] = Val;
275  return *this;
276 }
277 
278 /**
279 * Multiply a matrix to a value Val.
280 * Return new instance of class Matrix.
281 *
282 * \param Val - value Val
283 */
284 template<class T>
285 inline Matrix1D<T> Matrix1D<T>::operator* (const T Val) const {
286  Matrix1D<T> Tmp(*this);
287  for (int x = 0; x < this->size_x; x++)
288  Tmp[x] = matrix_array[x] * Val;
289  return Tmp;
290 }
291 
292 /**
293 * Divide a matrix to a value Val.
294 * Return new instance of class Matrix.
295 *
296 * \param Val - value Val
297 */
298 template<class T>
299 inline Matrix1D<T> Matrix1D<T>::operator/ (const T Val) const {
300  int x;
301  Matrix1D<T> Tmp(*this);
302  for (x = 0; x < this->size_x; x++)
303  Tmp[x] = matrix_array[x] / Val;
304  return Tmp;
305 }
306 
307 /**
308 * Multiply all values of matrix to values of matrix M.
309 *
310 * \param &M - matrix M.
311 */
312 template<class T>
313 inline Matrix1D<T> Matrix1D<T>::times (const Matrix1D<T> &M) const {
314  int x;
315  Matrix1D<T> Tmp(*this);
316  for (x = 0; x < this->size_x; x++)
317  Tmp[x] = matrix_array[x] * M.matrix_array[x];
318  return Tmp;
319 }
320 
321 /**
322 * Divide all values of matrix to values of matrix M.
323 *
324 * \param &M - matrix M.
325 */
326 template<class T>
327 inline Matrix1D<T> Matrix1D<T>::divide (const Matrix1D<T> &M) const {
328  int x;
329  Matrix1D<T> Tmp(*this);
330  for (x = 0; x < this->size_x; x++)
331  Tmp[x] = matrix_array[x] / M.matrix_array[x];
332  return Tmp;
333 }
334 
335 /**
336 * Norm of vector
337 */
338 template<class T>
339 inline T Matrix1D<T>::norm() const {
340  T res = 0;
341  int x;
342  for (x = 0; x < this->size_x; x++) {
343  res += matrix_array[x] * matrix_array[x];
344  }
345  return sqrt(res);
346 }
347 
348 /**
349 * Dot product
350 */
351 template<class T>
352 inline T Matrix1D<T>::dot( const Matrix1D<T> &W ) const {
353  if (this->size_x != W.size_x)
354  throw error_msg("DOT_PRODUCT", "Size is different");
355 
356  T res = 0;
357  int x;
358  for (x = 0; x < this->size_x; x++) {
359  res += matrix_array[x] * W[x];
360  }
361  return res;
362 }
363 
364 
365 /**
366 * Write matrix data to file.
367 */
368 template<class T>
369 void Matrix1D<T>::writeToFile(string filename) {
370  int x;
371  ofstream output(filename.c_str());
372  if (output==NULL && (filename.find("Debug") == string::npos)) {
373  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
374  }
375  output << "VARIABLES = \"" << ((this->name!="")?this->name:"function") << "\" "<< endl;
376  output << "ZONE T=\"" << filename << "\", I=" << size_x << endl;
377  output.setf(ios_base::scientific, ios_base::floatfield);
378  for (x = 0; x < size_x; x++) {
379  output << matrix_array[x] << endl;
380  }
381  output.close();
382 }
383 
384 /**
385 * Write matrix data to file with grid.
386 */
387 template<class T>
388 void Matrix1D<T>::writeToFile(string filename, Matrix1D<T> &grid_x) {
389  int x;
390  ofstream output(filename.c_str());
391  if (output==NULL && (filename.find("Debug") == string::npos)) {
392  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
393  }
394  output << "VARIABLES = \"" << ((grid_x.name!="")?grid_x.name:"x") << "\", \"" << ((this->name!="")?this->name:"function") << "\" "<< endl;
395  output << "ZONE T=\"" << filename << "\", I=" << size_x << endl;
396  output.setf(ios_base::scientific, ios_base::floatfield);
397  for (x = 0; x < size_x; x++) {
398  output << "\t" << grid_x[x] << "\t" << "\t" << matrix_array[x] << endl;
399  }
400  output.close();
401 }
402 
403 
404 /**
405 * Read matrix data from file.
406 */
407 template<class T>
408 void Matrix1D<T>::readFromFile(string filename) {
409  int x;
410  string inBuf;
411  if (!initialized) {
412  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
413  } else {
414  ifstream input(filename.c_str());
415  if (input != NULL && !input.eof()) {
416  // Skipping first two lines.
417  input >> inBuf;
418  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
419  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
420  input >> inBuf;
421  }
422  // read to the end of the line with 'zone'
423  input.ignore(9999, '\n');
424  for (x = 0; x < size_x; x++) {
425  input >> matrix_array[x];
426  }
427  } else {
428  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
429  //cout << "Error reading initial conditions input file " << filename << endl;
430  //getchar();
431  //exit(-1);
432  }
433  input.close();
434  }
435 }
436 
437 /**
438 * Read matrix data from file and check 'grid'
439 */
440 template<class T>
441 void Matrix1D<T>::readFromFile(string filename, Matrix1D<T> &grid_x) {
442  int x;
443  string inBuf;
444  double loaded_x;
445  double err = 1e-8;
446  if (!initialized) {
447  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
448  } else {
449  ifstream input(filename.c_str());
450  if (input != NULL && !input.eof()) {
451  // Skipping first two lines. Should serch for 'ZONE' and read from following line better - done.
452  /// \todo All matrix readings change to readFromFile() function
453  input >> inBuf;
454  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
455  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
456  input >> inBuf;
457  }
458  // read to the end of the line with 'zone'
459  input.ignore(9999, '\n');
460  for (x = 0; x < size_x; x++) {
461  input >> loaded_x;
462  // check if grid is the same
463  if (fabs(loaded_x - grid_x[x]) > err) {
464  throw error_msg("MATRIX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: %e\nGrid: %e\n", filename.c_str(), loaded_x, grid_x[x]);
465  } else {
466  input >> matrix_array[x];
467  }
468  }
469  } else {
470  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
471  //cout << "Error reading initial conditions input file " << filename << endl;
472  //getchar();
473  //exit(-1);
474  }
475  input.close();
476  }
477 }
478 
479 
480 //////////////////////////////////////////////////
481 //
482 // Matrix 2D
483 //
484 //////////////////////////////////////////////////
485 
486 /**
487 * Constructor.
488 * Allocate memory.
489 *
490 * \param x_size - x size
491 * \param y_size - y size
492 */
493 template<class T>
494 Matrix2D<T>::Matrix2D( int x_size, int y_size ) {
495  initialized = false;
496  // allocating memory
497  AllocateMemory(x_size, y_size);
498 }
499 
500 /**
501 * Constructor.
502 * Create new matrix from the Matrix M.
503 *
504 * \param &M - Matrix M
505 */
506 template<class T>
508  initialized = false;
509  this->operator = (M);
510 }
511 
512 /**
513 * Destructor.
514 */
515 template<class T>
517  if (initialized) free_matrix<T>(matrix_array);
518 }
519 
520 /**
521 * Allocate memory.
522 *
523 * \param x_size - x size
524 * \param y_size - y size
525 */
526 template<class T>
527 void Matrix2D<T>::AllocateMemory( int x_size, int y_size ) {
528  this->size_x = x_size;
529  this->size_y = y_size;
530  // using matrix inline template to allocate memory
531  matrix_array = matrix<T>(x_size, y_size);
532  initialized = true;
533  // for (int i = 0; i < size_x; i++)
534  // for (int j = 0; j < size_y; j++)
535  // matrix_array[i][j] = 0;
536 }
537 
538 /**
539 * Make matrix equal to Matrix M.
540 *
541 * \param &M - Matrix M.
542 */
543 template<class T>
545  // check, if both matrix are the same one
546  if (this!=&M) {
547  // check, if RHS matrix in initialized
548  if (M.initialized) {
549  // free LHS matrix if it was initialized
550  if (initialized && (size_x != M.size_x || size_y != M.size_y)) {
551  free_matrix<T>(matrix_array);
552  initialized = false;
553  }
554  if (!initialized) {
555  this->size_x = M.size_x;
556  this->size_y = M.size_y;
557  // allocating memory for LHS matrix (and creating the correct pointers M[0] - M[N] to the matrix rows)
558  this->AllocateMemory(this->size_x, this->size_y);
559  }
560  this->name = M.name;
561  // for (int x = 0; x < this->size_x; x++)
562  // memcpy( matrix_array[x], M.matrix_array[x], this->size_y * sizeof( T ) );
563  // fast values copying as a memory range
564  memcpy(matrix_array[0], M.matrix_array[0], this->size_x * this->size_y * sizeof(T) );
565  } else {
566  this->initialized = false;
567  }
568  }
569  return *this;
570 }
571 
572 /**
573 * Make matrix equal to value Val.
574 *
575 * \param val - value Val.
576 */
577 template<class T>
579  int x, y;
580  if (initialized) {
581  for (x = 0; x < this->size_x; x++)
582  for (y = 0; y < this->size_y; y++)
583  matrix_array[x][y] = val;
584  } else {
585  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
586  }
587 
588  return *this;
589 }
590 
591 /**
592 * Multiply matrix to Val.
593 *
594 * \param Val - value Val.
595 */
596 template<class T>
597 inline Matrix2D<T> Matrix2D<T>::operator* (const T Val) const {
598  int x, y;
599  Matrix2D<T> Tmp(*this);
600  for (x = 0; x < this->size_x; x++)
601  for (y = 0; y < this->size_y; y++)
602  Tmp[x][y] = matrix_array[x][y] * Val;
603  return Tmp;
604 }
605 
606 /**
607 * Divide matrix to Val.
608 *
609 * \param Val - value Val.
610 */
611 template<class T>
612 inline Matrix2D<T> Matrix2D<T>::operator/(const T Val) const {
613  int x, y;
614  Matrix2D<T> Tmp(*this);
615  for (x = 0; x < this->size_x; x++)
616  for (y = 0; y < this->size_y; y++)
617  Tmp[x][y] = matrix_array[x][y] / Val;
618  return Tmp;
619 }
620 
621 /**
622 * Divide all values of matrix to values of matrix M.
623 *
624 * \param &M - matrix M.
625 */
626 template<class T>
627 inline Matrix2D<T> Matrix2D<T>::divide (const Matrix2D<T> &M) const {
628  int x, y;
629  Matrix2D<T> Tmp(*this);
630  for (x = 0; x < this->size_x; x++)
631  for (y = 0; y < this->size_y; y++)
632  Tmp[x][y] = matrix_array[x][y] / M[x][y];
633  return Tmp;
634 }
635 
636 /**
637 * Multiply all values of matrix to values of matrix M.
638 *
639 * \param &M - matrix M.
640 */
641 template<class T>
642 inline Matrix2D<T> Matrix2D<T>::times (const Matrix2D<T> &M) const {
643  int x, y;
644  Matrix2D<T> Tmp(*this);
645  for (x = 0; x < this->size_x; x++)
646  for (y = 0; y < this->size_y; y++)
647  Tmp[x][y] = matrix_array[x][y] * M[x][y];
648  return Tmp;
649 }
650 
651 /**
652 * Divide all values of matrix to value Val.
653 *
654 * \param val - value to divide.
655 */
656 template<class T>
658  int x, y;
659  Matrix2D<T> Tmp(*this);
660  for (x = 0; x < this->size_x; x++)
661  for (y = 0; y < this->size_y; y++)
662  Tmp[x][y] = (matrix_array[x][y]>val)?matrix_array[x][y]:val;
663  return Tmp;
664 }
665 
666 /**
667 * Returns corresponding index of 1d array
668 */
669 template<class T>
670 inline int Matrix2D<T>::index1d(int x, int y) const {
671  return x*size_y + y;
672 }
673 
674 
675 
676 /**
677 * Writes the matrix to a file.
678 * File has two header lines.
679 *
680 * \todo Change the name 'writeToFile' in Matrix classes to 'write' or something.
681 *
682 * \param filename - file name
683 */
684 template<class T>
685 void Matrix2D<T>::writeToFile(string filename) {
686  int x, y;
687  ofstream output(filename.c_str());
688  if (output==NULL && (filename.find("Debug") == string::npos)) {
689  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
690  }
691  output << "VARIABLES = \""<< ((this->name!="")?this->name:"f") << "\" "<< endl;
692  output << "ZONE T=\"" << filename << "\", I=" << size_y << ", J= " << size_x << endl;
693  output.setf(ios_base::scientific, ios_base::floatfield);
694  for (x = 0; x < size_x; x++) {
695  for (y = 0; y < size_y; y++) {
696  output << matrix_array[x][y] << endl;
697  }
698  }
699  output.close();
700 }
701 
702 
703 /**
704 * Write the matrix to a file using other two matrixes as a grid.
705 * Simply that means - write all three matrixes to a file.
706 */
707 template<class T>
708 void Matrix2D<T>::writeToFile(string filename, Matrix2D<T> &grid_x, Matrix2D<T> &grid_y) {
709  int x, y;
710  ofstream output(filename.c_str());
711  output << "VARIABLES = \"" << ((grid_x.name!="")?grid_x.name:"x") << "\", \"" << ((grid_y.name!="")?grid_y.name:"y") << "\", \"" << ((this->name!="")?this->name:"f") << "\" "<< endl;
712  output << "ZONE T=\"" << filename << "\", I=" << size_y << ", J=" << size_x << endl;
713  output.setf(ios_base::scientific, ios_base::floatfield);
714  for (x = 0; x < size_x; x++) {
715  for (y = 0; y < size_y; y++) {
716  output << "\t" << grid_x[x][y] << "\t" << grid_y[x][y] << "\t" << matrix_array[x][y] << endl;
717  }
718  }
719  output.close();
720 }
721 
722 
723 /**
724 * Read matrix data from file.
725 */
726 template<class T>
727 void Matrix2D<T>::readFromFile(string filename) {
728  int x, y;
729  string inBuf;
730  if (!initialized) {
731  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
732  } else {
733  ifstream input(filename.c_str());
734  if (input != NULL && !input.eof()) {
735  // Skipping first two lines.
736  input >> inBuf;
737  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
738  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
739  input >> inBuf;
740  }
741  // read to the end of the line with 'zone'
742  input.ignore(9999, '\n');
743  for (x = 0; x < size_x; x++) {
744  for (y = 0; y < size_y; y++) {
745  input >> matrix_array[x][y];
746  }
747  }
748  } else {
749  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
750  //cout << "Error reading initial conditions input file " << filename << endl;
751  //getchar();
752  //exit(-1);
753  }
754  input.close();
755  }
756 }
757 
758 /**
759 * Read matrix data from file and check grid
760 */
761 template<class T>
762 void Matrix2D<T>::readFromFile(string filename, Matrix2D<T> &grid_x, Matrix2D<T> &grid_y) {
763  int x, y;
764  string inBuf;
765  double err = 1e-6;
766  double loaded_x, loaded_y;
767  if (!initialized) {
768  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
769  } else {
770  ifstream input(filename.c_str());
771  if (input != NULL && !input.eof()) {
772  // Skipping first two lines.
773  input >> inBuf;
774  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
775  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
776  input >> inBuf;
777  }
778  // read to the end of the line with 'zone'
779  input.ignore(9999, '\n');
780  for (x = 0; x < size_x; x++) {
781  for (y = 0; y < size_y; y++) {
782  input >> loaded_x >> loaded_y;
783  // check if grid is the same
784  if (fabs(loaded_x - grid_x[x][y]) > err || fabs(loaded_y - grid_y[x][y]) > err) {
785  throw error_msg("MATRIX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: %e, %e\nGrid: %e, %e\n", filename.c_str(), loaded_x, loaded_y, grid_x[x][y], grid_y[x][y]);
786  } else {
787  input >> matrix_array[x][y];
788  }
789  }
790  }
791  } else {
792  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
793  //cout << "Error reading initial conditions input file " << filename << endl;
794  //getchar();
795  //exit(-1);
796  }
797  input.close();
798  }
799 }
800 
801 
802 //////////////////////////////////////////////////
803 //
804 // Matrix 3D
805 //
806 //////////////////////////////////////////////////
807 
808 /**
809 * Constructor.
810 * Allocate memory.
811 */
812 template<class T>
813 Matrix3D<T>::Matrix3D( int x_size, int y_size, int z_size ) {
814  initialized = false;
815  // allocating memory
816  AllocateMemory(x_size, y_size, z_size);
817 }
818 
819 /**
820 * Constructor.
821 * Create matrix equal to Matrix M.
822 *
823 * \param &M - Matrix M.
824 */
825 template<class T>
827  initialized = false;
828  this->operator = (M);
829 }
830 
831 /**
832 * Destructor
833 */
834 template<class T>
836  if (initialized) free_matrix<T>(matrix_array, size_x, size_y);
837  // delete[] plane_array; - it's get deleted with matrix_array[0][0][0] or something?
838 }
839 
840 /**
841 * Allocating memory and filling it with zero-values.
842 */
843 template<class T>
844 void Matrix3D<T>::AllocateMemory( int x_size, int y_size, int z_size) {
845  this->size_x = x_size;
846  this->size_y = y_size;
847  this->size_z = z_size;
848  matrix_array = matrix<T>(x_size, y_size, z_size);
849  plane_array = matrix_array[0][0];
850  initialized = true;
851 #ifdef DEBUG_MODE
852  // should not initialize matrix with zeros, it can slow the code greatly in some cases
853  for (int i = 0; i < size_x; i++)
854  for (int j = 0; j < size_y; j++)
855  for (int k = 0; k < size_z; k++)
856  matrix_array[i][j][k] = 0;
857 #endif
858 }
859 
860 /**
861 * Operator [i], returns pointer to 2D array. Next [j][k] can be applied to return value.
862 * If DEBUG_MODE defined, check if matrix has been initialized.
863 *
864 * \param i - number of element to return
865 */
866 template<class T>
867 inline T** Matrix3D<T>::operator[] (int i) {
868 #ifdef DEBUG_MODE
869  if (!initialized) {
870  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
871  }
872 #endif
873  // return i-th pointer
874  return matrix_array[i];
875 }
876 
877 /**
878 * Operator (x, y, z), returns value of element [x][y][z].
879 * If DEBUG_MODE defined, check if matrix has been initialized.
880 *
881 */
882 template<class T>
883 inline T& Matrix3D<T>::operator() (int x, int y, int z) {
884 #ifdef DEBUG_MODE
885  if (!initialized) {
886  throw error_msg("MATRIX_ERROR", "Using not initialized matrix");
887  }
888  if ((x < 0 || x > size_x-1) || (y < 0 || y > size_y-1) || (z < 0 || z > size_z-1)) {
889  throw error_msg("MATRIX_ERROR", "Index is out of bound");
890  }
891 #endif
892  // return (x,y,z) value
893  return plane_array[(x*size_y + y)*size_z + z];
894 }
895 
896 /**
897 * Makes matrix equal to Matrix M.
898 *
899 * \param &M - Matrix M.
900 */
901 template<class T>
903  // check if not LHS and RHS matrix are the same
904  if (this!=&M) {
905  // check if RHS matrix is initialized
906  if (M.initialized) {
907  // free LHS matrix, if it is initialized
908  if (initialized && (size_x != M.size_x || size_y != M.size_y || size_z != M.size_z)) {
909  free_matrix<T>(matrix_array, this->size_x, this->size_y);
910  initialized = false;
911  }
912  if (!initialized) {
913  this->size_x = M.size_x;
914  this->size_y = M.size_y;
915  this->size_z = M.size_z;
916  // allocating memory for LHS matrix
917  this->AllocateMemory(this->size_x, this->size_y, this->size_z);
918  }
919  this->name = M.name;
920  // for (int x = 0; x < this->size_x; x++)
921  // for (int y = 0; y < this->size_y; y++)
922  // memcpy( matrix_array[x][y], M.matrix_array[x][y], this->size_z * sizeof( T ) );
923  // fast values copying as memory range
924  memcpy( matrix_array[0][0], M.matrix_array[0][0], this->size_x * this->size_y * this->size_z * sizeof( T ) );
925  } else {
926  this->initialized = false;
927  }
928  }
929  return *this;
930 }
931 
932 /**
933 * Makes 3D matrix from 2D matrix.
934 * The 3rd dimension makes equal to 1.
935 */
936 /* template<class T>
937 Matrix3D<T>& Matrix3D<T>::operator= (const Matrix2D<T> &M) {
938 double val;
939 if (M.initialized) {
940 if (initialized) free_matrix<T>(matrix_array, this->size_x, this->size_y);
941 this->size_x = 1;
942 this->size_y = M.size_x;
943 this->size_z = M.size_y;
944 //this->name = M.name;
945 this->AllocateMemory(this->size_x, this->size_y, this->size_z);
946 for (int y = 0; y < this->size_y; y++)
947 //memcpy( matrix_array[0][y], M.matrix_array[y], this->size_z * sizeof( T ) );
948 for (int z = 0; z < this->size_z; z++) {
949 this->matrix_array[0][y][z] = M[y][z];
950 }
951 } else {
952 this->initialized = false;
953 }
954 return *this;
955 }*/
956 
957 /**
958 * Makes Matrix equal to value Val.
959 */
960 template<class T>
962  // std::fill_n(matrix_array[0][0][0], size_x*size*y*size*z, Val); // might be faster
963  for (int x = 0; x < this->size_x; x++)
964  for (int y = 0; y < this->size_y; y++)
965  for (int z = 0; z < this->size_z; z++)
966  this->matrix_array[x][y][z] = Val;
967  return *this;
968 }
969 
970 
971 /**
972 * Matrix summation, result is stored into applied matrix (left hand side matrix)
973 */
974 template<class T>
976  int x, y, z;
977  for (x = 0; x < this->size_x; x++)
978  for (y = 0; y < this->size_y; y++)
979  for (z = 0; z < this->size_z; z++)
980  matrix_array[x][y][z] += M.matrix_array[x][y][z];
981  return *this;
982 }
983 
984 /**
985 * Matrix subtraction, result is stored into applied matrix (left hand side matrix)
986 */
987 template<class T>
989  int x, y, z;
990  for (x = 0; x < this->size_x; x++)
991  for (y = 0; y < this->size_y; y++)
992  for (z = 0; z < this->size_z; z++)
993  matrix_array[x][y][z] -= M.matrix_array[x][y][z];
994  return *this;
995 }
996 
997 /**
998 * Multiplication to a value. Result is stored into applied matrix (left hand side matrix)
999 */
1000 template<class T>
1002  int x, y, z;
1003  for (x = 0; x < this->size_x; x++)
1004  for (y = 0; y < this->size_y; y++)
1005  for (z = 0; z < this->size_z; z++)
1006  matrix_array[x][y][z] *= Val;
1007  return *this;
1008 }
1009 
1010 /**
1011 * Division by a value. Result is stored into applied matrix (left hand side matrix)
1012 */
1013 template<class T>
1015  int x, y, z;
1016  for (x = 0; x < this->size_x; x++)
1017  for (y = 0; y < this->size_y; y++)
1018  for (z = 0; z < this->size_z; z++)
1019  matrix_array[x][y][z] /= Val;
1020  return *this;
1021 }
1022 
1023 /**
1024 * Summation with a value. Result is stored into applied matrix (left hand side matrix)
1025 */
1026 template<class T>
1028  int x, y, z;
1029  for (x = 0; x < this->size_x; x++)
1030  for (y = 0; y < this->size_y; y++)
1031  for (z = 0; z < this->size_z; z++)
1032  matrix_array[x][y][z] += Val;
1033  return *this;
1034 }
1035 
1036 /**
1037 * Subtraction of a value. Result is stored into applied matrix (left hand side matrix)
1038 */
1039 template<class T>
1041  int x, y, z;
1042  for (x = 0; x < this->size_x; x++)
1043  for (y = 0; y < this->size_y; y++)
1044  for (z = 0; z < this->size_z; z++)
1045  matrix_array[x][y][z] -= Val;
1046  return *this;
1047 }
1048 
1049 /**
1050 * Multiplication between each element of the matrices (not a matrix multiplication). Result is stored into applied matrix (left hand side matrix)
1051 */
1052 template<class T>
1054  int x, y, z;
1055  for (x = 0; x < this->size_x; x++)
1056  for (y = 0; y < this->size_y; y++)
1057  for (z = 0; z < this->size_z; z++)
1058  matrix_array[x][y][z] *= M.matrix_array[x][y][z];
1059  return *this;
1060 }
1061 
1062 /**
1063 * Division of each element of one matrices to the element of another. Result is stored into applied matrix (left hand side matrix)
1064 */
1065 template<class T>
1067  int x, y, z;
1068  for (x = 0; x < this->size_x; x++)
1069  for (y = 0; y < this->size_y; y++)
1070  for (z = 0; z < this->size_z; z++)
1071  matrix_array[x][y][z] /= M.matrix_array[x][y][z];
1072  return *this;
1073 }
1074 
1075 /**
1076 * Add each element of the matrix to corresponds element of matrix M.
1077 */
1078 template<class T>
1080  int x, y, z;
1081  Matrix3D<T> Tmp(*this);
1082  for (x = 0; x < this->size_x; x++)
1083  for (y = 0; y < this->size_y; y++)
1084  for (z = 0; z < this->size_z; z++)
1085  Tmp[x][y][z] = matrix_array[x][y][z] + M.matrix_array[x][y][z];
1086  return Tmp;
1087 }
1088 
1089 /**
1090 * Substract each element of the matrix to corresponds element of matrix M.
1091 */
1092 template<class T>
1094  int x, y, z;
1095  Matrix3D<T> Tmp(*this);
1096  for (x = 0; x < this->size_x; x++)
1097  for (y = 0; y < this->size_y; y++)
1098  for (z = 0; z < this->size_z; z++)
1099  Tmp[x][y][z] = matrix_array[x][y][z] - M.matrix_array[x][y][z];
1100  return Tmp;
1101 }
1102 
1103 
1104 /**
1105 * Multiply each element of the matrix to Val, save result to a new matrix.
1106 */
1107 template<class T>
1108 inline Matrix3D<T> Matrix3D<T>::operator* (const T Val) const {
1109  int x, y, z;
1110  Matrix3D<T> Tmp(*this);
1111  for (x = 0; x < this->size_x; x++)
1112  for (y = 0; y < this->size_y; y++)
1113  for (z = 0; z < this->size_z; z++)
1114  Tmp[x][y][z] = matrix_array[x][y][z] * Val;
1115  return Tmp;
1116 }
1117 
1118 
1119 /**
1120 * Divide each element of the matrix to Val, save result to a new matrix.
1121 */
1122 template<class T>
1123 inline Matrix3D<T> Matrix3D<T>::operator/ (const T Val) const {
1124  int x, y, z;
1125  Matrix3D<T> Tmp(*this);
1126  for (x = 0; x < this->size_x; x++)
1127  for (y = 0; y < this->size_y; y++)
1128  for (z = 0; z < this->size_z; z++)
1129  Tmp[x][y][z] = matrix_array[x][y][z] / Val;
1130  return Tmp;
1131 }
1132 
1133 /**
1134 * Multiply each element of the matrix to corresponds element of matrix M.
1135 */
1136 template<class T>
1137 inline Matrix3D<T> Matrix3D<T>::times (const Matrix3D<T> &M) const {
1138  int x, y, z;
1139  Matrix3D<T> Tmp(*this);
1140  for (x = 0; x < this->size_x; x++)
1141  for (y = 0; y < this->size_y; y++)
1142  for (z = 0; z < this->size_z; z++)
1143  Tmp[x][y][z] = matrix_array[x][y][z] * M.matrix_array[x][y][z];
1144  return Tmp;
1145 }
1146 
1147 /**
1148 * Divide each element of the matrix to corresponds element of matrix M.
1149 */
1150 template<class T>
1152  int x, y, z;
1153  Matrix3D<T> Tmp(*this);
1154  for (x = 0; x < this->size_x; x++)
1155  for (y = 0; y < this->size_y; y++)
1156  for (z = 0; z < this->size_z; z++)
1157  Tmp[x][y][z] = matrix_array[x][y][z] / M.matrix_array[x][y][z];
1158  return Tmp;
1159 }
1160 
1161 
1162 /**
1163 * Write matrix to file.
1164 * File has two header lines.
1165 */
1166 template<class T>
1167 void Matrix3D<T>::writeToFile(string filename) {
1168  int x, y, z;
1169  ofstream output(filename.c_str());
1170  output << "VARIABLES = \""<< ((this->name!="")?this->name:"f") <<"\" "<< endl;
1171  output << "ZONE T=\"" << filename << "\", I=" << size_z << ", J=" << size_y << ", K=" << size_x << endl;
1172  output.setf(ios_base::scientific, ios_base::floatfield);
1173  for (x = 0; x < size_x; x++) {
1174  for (y = 0; y < size_y; y++) {
1175  for (z = 0; z < size_z; z++) {
1176  output << matrix_array[x][y][z] << endl;
1177  }
1178  }
1179  }
1180  output.close();
1181 }
1182 
1183 /**
1184 * Write matrix to file, using 3 other matrixes as a grid (simply - write all 4 matrixes to the file).
1185 * File has two header lines.
1186 */
1187 template<class T>
1188 void Matrix3D<T>::writeToFile(string filename, Matrix3D<T> &grid_x, Matrix3D<T> &grid_y, Matrix3D<T> &grid_z) {
1189  int x, y, z;
1190  ofstream output(filename.c_str());
1191  if (output==NULL && (filename.find("Debug") == string::npos)) {
1192  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
1193  }
1194  output << "VARIABLES = \"" << ((grid_x.name!="")?grid_x.name:"x") << "\", \"" << ((grid_y.name!="")?grid_y.name:"y") << "\", \"" << ((grid_z.name!="")?grid_z.name:"z") << "\", \"" << ((this->name!="")?this->name:"f") << "\" "<< endl;
1195  output << "ZONE T=\"" << filename << "\", I=" << size_z << ", J=" << size_y << ", K=" << size_x << endl;
1196  output.setf(ios_base::scientific, ios_base::floatfield);
1197  for (x = 0; x < size_x; x++) {
1198  for (y = 0; y < size_y; y++) {
1199  for (z = 0; z < size_z; z++) {
1200  output << "\t" << grid_x[x][y][z] << "\t" << grid_y[x][y][z] << "\t" << grid_z[x][y][z] << "\t" << matrix_array[x][y][z] << endl;
1201  }
1202  }
1203  }
1204  output.close();
1205 }
1206 
1207 /**
1208 * Read matrix data from file.
1209 */
1210 template<class T>
1211 void Matrix3D<T>::readFromFile(string filename) {
1212  int x, y, z;
1213  string inBuf;
1214  if (!initialized) {
1215  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
1216  } else {
1217  ifstream input(filename.c_str());
1218  if (input != NULL && !input.eof()) {
1219  // Skipping first two lines.
1220  input >> inBuf;
1221  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
1222  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
1223 
1224  input >> inBuf;
1225  }
1226  // read to the end of the line with 'zone'
1227  input.ignore(9999, '\n');
1228  for (x = 0; x < size_x; x++) {
1229  for (y = 0; y < size_y; y++) {
1230  for (z = 0; z < size_z; z++) {
1231  input >> matrix_array[x][y][z];
1232  }
1233  }
1234  }
1235  } else {
1236  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
1237  //cout << "Error reading initial conditions input file " << filename << endl;
1238  //getchar();
1239  //exit(-1);
1240  }
1241  input.close();
1242  }
1243 }
1244 
1245 /**
1246 * Read matrix data from file with grid
1247 */
1248 template<class T>
1249 void Matrix3D<T>::readFromFile(string filename, Matrix3D<T> &grid_x, Matrix3D<T> &grid_y, Matrix3D<T> &grid_z) {
1250  int x, y, z;
1251  string inBuf;
1252  if (!initialized) {
1253  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
1254  } else {
1255  ifstream input(filename.c_str());
1256  if (input != NULL && !input.eof()) {
1257  // Skipping first two lines.
1258  input >> inBuf;
1259  // !!! while (strcmp(_strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
1260  // !!! while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
1261  while (strcmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
1262  input >> inBuf;
1263  }
1264  // read to the end of the line with 'zone'
1265  input.ignore(9999, '\n');
1266  for (x = 0; x < size_x; x++) {
1267  for (y = 0; y < size_y; y++) {
1268  for (z = 0; z < size_z; z++) {
1269  input >> grid_x[x][y][z] >> grid_y[x][y][z] >> grid_z[x][y][z] >> matrix_array[x][y][z];
1270  // skip till the end of the line
1271  input.ignore(9999, '\n');
1272 
1273  // input >> loaded_x >> loaded_y >> loaded_z;
1274  // check if grid is the same
1275  // if (fabs(loaded_x - grid_x[x][y][z]) > err || fabs(loaded_y - grid_y[x][y][z]) > err || fabs(loaded_z - grid_z[x][y][z]) > err) {
1276  // throw error_msg("MATRIX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: %e, %e, %e\nGrid: %e, %e, %e\n", filename.c_str(), loaded_x, loaded_y, loaded_z, grid_x[x][y][z], grid_y[x][y][z], grid_z[x][y][z]);
1277  //} else {
1278  // input >> matrix_array[x][y][z];
1279  //}
1280  }
1281  }
1282  }
1283  } else {
1284  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
1285  //cout << "Error reading initial conditions input file " << filename << endl;
1286  //getchar();
1287  //exit(-1);
1288  }
1289  input.close();
1290  }
1291 }
1292 
1293 
1294 
1295 /**
1296 * Returns corresponding index of 1d array
1297 */
1298 template<class T>
1299 inline int Matrix3D<T>::index1d(int x, int y, int z) {
1300  return (x*size_y + y)*size_z + z;
1301 }
1302 
1303 
1304 /**
1305 * Return maximum value of the matrix.
1306 */
1307 template<class T>
1309  T tmp = 0;
1310  int x, y, z;
1311  for (x = 0; x < size_x; x++) {
1312  for (y = 0; y < size_y; y++) {
1313  for (z = 0; z < size_z; z++) {
1314  tmp = (tmp>matrix_array[x][y][z])?tmp:matrix_array[x][y][z];
1315  }
1316  }
1317  }
1318  return tmp;
1319 }
1320 
1321 /**
1322 * Return absolute maximum value of the matrix.
1323 */
1324 template<class T>
1326  T tmp = 0;
1327  int x, y, z;
1328  for (x = 0; x < size_x; x++) {
1329  for (y = 0; y < size_y; y++) {
1330  for (z = 0; z < size_z; z++) {
1331  tmp = (tmp>fabs((double)matrix_array[x][y][z]))?tmp:fabs((double)matrix_array[x][y][z]);
1332  }
1333  }
1334  }
1335  return tmp;
1336 }
1337 
1338 /**
1339 * Return absolute value of the matrix.
1340 */
1341 template<class T>
1343  Matrix3D<T> tmp(this->size_x, this->size_y, this->size_z);
1344  int x, y, z;
1345  for (x = 0; x < size_x; x++) {
1346  for (y = 0; y < size_y; y++) {
1347  for (z = 0; z < size_z; z++) {
1348  tmp[x][y][z] = (matrix_array[x][y][z]>0)?matrix_array[x][y][z]:-matrix_array[x][y][z];
1349  }
1350  }
1351  }
1352  return tmp;
1353 }
1354 
1355 
1356 /**
1357 * Make x-slice of 3d matrix - 2d matrix.
1358 */
1359 template<class T>
1361  int y, z;
1362  Matrix2D<T> tmp(this->size_y, this->size_z);
1363  tmp.name = this->name + "_slice";
1364  for (y = 0; y < size_y; y++) {
1365  for (z = 0; z < size_z; z++) {
1366  tmp[y][z] = matrix_array[p_x][y][z];
1367  }
1368  }
1369  return tmp;
1370 }
1371 
1372 /**
1373 * Make y-slice of 3d matrix - 2d matrix.
1374 */
1375 template<class T>
1377  int x, z;
1378  Matrix2D<T> tmp(this->size_x, this->size_z);
1379  tmp.name = this->name + "_slice";
1380  for (x = 0; x < size_x; x++) {
1381  for (z = 0; z < size_z; z++) {
1382  tmp[x][z] = matrix_array[x][p_y][z];
1383  }
1384  }
1385  return tmp;
1386 }
1387 
1388 /**
1389 * Make z-slice of 3d matrix - 2d matrix.
1390 */
1391 template<class T>
1393  int x, y;
1394  Matrix2D<T> tmp(this->size_x, this->size_y);
1395  tmp.name = this->name + "_slice";
1396  for (x = 0; x < size_x; x++) {
1397  for (y = 0; y < size_y; y++) {
1398  tmp[x][y] = matrix_array[x][y][p_z];
1399  }
1400  }
1401  return tmp;
1402 }
1403 
1404 
1405 /////////////////////////////////////////////////////////////////////////////
1406 //
1407 // Interpolations (mess)
1408 //
1409 //////////////////////////////////////////////////////////////////////////////
1410 
1411 
1412 /**
1413 * Linear interpolation.
1414 * Runs interpolation function from Numerical Recepies.
1415 */
1416 template<class T>
1417 void Matrix1D<T>::Interpolate(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
1418 
1419 
1420  //if (size_x != old_grid.size_x) {
1421  if (this->size_x != new_grid.size_x) {
1422  throw error_msg("INTERPOLATION_ERROR", "New grid and function size mismatch!\n");
1423  //cout << "new grid and function size mismatch!" << endl;
1424  //return;
1425  }
1426 
1427  //Matrix1D<T> res(new_grid.size_x);
1428  //Maths::Interpolation::Linear LinInterp(old_grid.size_x, old_grid.matrix_array, matrix_array);
1429  // initialize some class. It was like that when I downloaded the function.
1430  Maths::Interpolation::Linear LinInterp(old_grid.size_x, old_grid.matrix_array, old_function.matrix_array);
1431 
1432  int i;
1433  for (i = 0; i < new_grid.size_x; i++)
1434  //res[i] = LinInterp.getValue(new_grid.matrix_array[i]);
1435  matrix_array[i] = LinInterp.getValue(new_grid.matrix_array[i]);
1436 
1437  //return res;
1438 
1439 }
1440 
1441 /**
1442 * Spline interpolation.
1443 * Runs interpolation functions from Numerical Recepies.
1444 * Treats boundary conditions separete, i.e. use boundary values on baundary conditions and extrapolation, use linear interpolation for next-to-the-boundary points.
1445 * Has additional mechenism of smoothning the interpolation - makes it more linear according the parameters lin_spline_coef and max_second_der.
1446 * Spline interpolation is linear + additional terms from second derivatives (look Numerical Recepies). If second derivative > max_second_der the second derivatives are multiplid by lin_spline_coef.
1447 *
1448 * \param &old_function - old function
1449 * \param &old_grid - old grid
1450 * \param &new_grid - new grid
1451 * \param lb - low boundary value
1452 * \param ub - upper boundary value
1453 * \param lin_spline_coef - coefficient to multiply the second derivatives terms (makes interpolation more smooth)
1454 * \param max_second_der - maximum second derivatives, all greater derivatives gonna be multiplied by lin_spline_coef
1455 *
1456 */
1457 template<class T>
1458 void Matrix1D<T>::Spline(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub, double lin_spline_coef, double max_second_der) {
1459 
1460  //if (size_x != old_grid.size_x) {
1461  // check, if size of the grids are equal
1462  if (this->size_x != new_grid.size_x) {
1463  throw error_msg("SPLINE_INTERPOLATION", "New grid and function size mismatch!");
1464  }
1465  // if size == 1 - just copy the value
1466  if (this->size_x == 1) {
1467  this->matrix_array[0] = old_function.matrix_array[0];
1468  return;
1469  }
1470 
1471  // Idea was to find the intersection between grids, to interpolate insite the intersection and do not extrapolate outside. But is done later.
1472  int i=0, i_min=0, i_max=old_grid.size_x-1, interp_size=0;
1473  // for(i = 0; old_grid.matrix_array[i] <= new_grid.matrix_array[0] && i < old_grid.size_x-1; i++) {}
1474  // i_min = i;
1475  // for(i = old_grid.size_x-1; old_grid.matrix_array[i] >= new_grid.matrix_array[old_grid.size_x-1] && i > 0; i--) {}
1476  // i_max = i;
1477  //cout << i_min << " " << i_max << endl;
1478 
1479  // new range for spline interpolation - values close to the boundaryes are interpolated by linear interpolation
1480  i_min = 1; i_max = old_grid.size_x-2;
1481  interp_size = i_max-i_min+1;
1482 
1483  double yp0, ypn;
1484  //yp0 = 2*(old_function[1] - old_function[0])/(old_grid[1] - old_grid[0]) - (old_function[2] - old_function[1])/(old_grid[2] - old_grid[1]);
1485  //ypn = 2*(old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]) - (old_function[size_x-2] - old_function[size_x-3])/(old_grid[size_x-2] - old_grid[size_x-3]);
1486  // ypn = (old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]);
1487 
1488  // flag, means use smooth second derivative on the boundary
1489  yp0 = 1e99;
1490  //yp0 = -(old_function[2] - old_function[0])/(old_grid[2] - old_grid[0]);
1491  ypn = 1e99;//0;
1492 
1493  Matrix1D<T> y2(old_grid.size_x);
1494  //spline(old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, yp0, ypn, y2.matrix_array);
1495  //spline(&old_grid.matrix_array[1], &old_function.matrix_array[1], old_grid.size_x-2, yp0, ypn, &y2.matrix_array[1]);
1496  // calculating array of second derivatives
1497  spline(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], interp_size, yp0, ypn, &y2.matrix_array[i_min]);
1498 
1499  // make it more linear (lower second derivatives)
1500  for (i = 0; i < old_grid.size_x; i++) {
1501  // if second derivative value is too big
1502  if (fabs(y2[i]) >= max_second_der) {
1503  if (y2[i] > 0) {
1504  // decreasing by lin_spline_coef big second derivative values - makes interpolation more linear
1505  // see numerical recipes
1506  y2[i] = max_second_der + (y2[i] - max_second_der)*(1.0 - lin_spline_coef);
1507  } else {
1508  y2[i] = - max_second_der + (y2[i] + max_second_der)*(1.0 - lin_spline_coef);
1509  }
1510  }
1511  }
1512 
1513  // for all grid points - looking for new values
1514  double a; int k; // interpolation parameters
1515  for (i = 0; i < new_grid.size_x; i++) {
1516  if (i == 0) { // if we are on the boundary
1517  this->matrix_array[i] = lb;
1518  } else if (i == new_grid.size_x-1) { // if we are on the boundary
1519  this->matrix_array[i] = ub;
1520  } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[0] ) { // if it is extrapolation
1521  this->matrix_array[i] = lb;
1522  //this->matrix_array[i] = old_function.matrix_array[0];
1523  } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-1]) { // if it is extrapolation
1524  this->matrix_array[i] = ub;
1525  //this->matrix_array[i] = old_function.matrix_array[old_function.size_x-1];
1526  } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[1] ) { // if the new grid point is between two boundary points of the old grid
1527  // linear
1528  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[0])/(old_grid.matrix_array[1] - old_grid.matrix_array[0]);
1529  this->matrix_array[i] = old_function.matrix_array[0] + a*(old_function.matrix_array[1] - old_function.matrix_array[0]);
1530  } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-2]) { // if the new grid point is between two boundary points of the old grid
1531  // linear
1532  //double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-2])/(old_grid.matrix_array[old_grid.size_x-1] - old_grid.matrix_array[old_grid.size_x-2]);
1533  //this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-2] + a*(old_function.matrix_array[old_grid.size_x-1] - old_function.matrix_array[old_grid.size_x-2]);
1534  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-1])/(old_grid.matrix_array[old_grid.size_x-2] - old_grid.matrix_array[old_grid.size_x-1]);
1535  this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-1] + a*(old_function.matrix_array[old_grid.size_x-2] - old_function.matrix_array[old_grid.size_x-1]);
1536  /* } else if (i <= 1) { // !!!!!!!! always linear for the first two points ???
1537  int k = 1;
1538  while(old_grid.matrix_array[k] < new_grid.matrix_array[i]) {k++;}
1539  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[k-1])/(old_grid.matrix_array[k] - old_grid.matrix_array[k-1]);
1540  this->matrix_array[i] = old_function.matrix_array[k-1] + a*(old_function.matrix_array[k] - old_function.matrix_array[k-1]);
1541  } else if (i >= new_grid.size_x-2) { // !!!!!!!! always linear for the last two points ???
1542  int k = new_grid.size_x-1;
1543  while(old_grid.matrix_array[k] > new_grid.matrix_array[i]) {k--;}
1544  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[k])/(old_grid.matrix_array[k+1] - old_grid.matrix_array[k]);
1545  this->matrix_array[i] = old_function.matrix_array[k] + a*(old_function.matrix_array[k+1] - old_function.matrix_array[k]);*/
1546  } else {
1547  //splint(old_grid.matrix_array, old_function.matrix_array, y2.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i]);
1548  //splint(&old_grid.matrix_array[1], &old_function.matrix_array[1], &y2.matrix_array[1], old_grid.size_x-2, new_grid.matrix_array[i], &this->matrix_array[i]);
1549  // Shouldn't we pass whole matrix for the least argument????
1550  /// \todo Check the index of the last argument of the spline interpolation.
1551  splint(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], &y2.matrix_array[i_min], interp_size, new_grid.matrix_array[i], &this->matrix_array[i]);
1552  }
1553  }
1554 
1555 }
1556 
1557 /*template<class T>
1558 void Matrix1D<T>::Spline(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double ub, double lb) {
1559 
1560 //if (size_x != old_grid.size_x) {
1561 if (this->size_x != new_grid.size_x) {
1562 cout << "new grid and function size mismatch!" << endl;
1563 return;
1564 }
1565 
1566 double yp0, ypn;
1567 //yp0 = 2*(old_function[1] - old_function[0])/(old_grid[1] - old_grid[0]) - (old_function[2] - old_function[1])/(old_grid[2] - old_grid[1]);
1568 //ypn = 2*(old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]) - (old_function[size_x-2] - old_function[size_x-3])/(old_grid[size_x-2] - old_grid[size_x-3]);
1569 // ypn = (old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]);
1570 
1571 yp0 = 1e99;//-(old_function[1] - old_function[0])/(old_grid[1] - old_grid[0]);
1572 ypn = 1e99;//0;
1573 
1574 Matrix1D<double> y2(old_grid.size_x);
1575 //spline(old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, yp0, ypn, y2.matrix_array);
1576 spline(&old_grid.matrix_array[1], &old_function.matrix_array[1], old_grid.size_x-2, yp0, ypn, &y2.matrix_array[1]);
1577 
1578 int i;
1579 for (i = 0; i < new_grid.size_x; i++) {
1580 if (i == 0) {
1581 this->matrix_array[i] = lb;
1582 } else if (i == new_grid.size_x-1) {
1583 this->matrix_array[i] = ub;
1584 } else
1585 if (new_grid.matrix_array[i] <= old_grid.matrix_array[0] ) {
1586 this->matrix_array[i] = lb;//old_function.matrix_array[0];
1587 } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[1] ) {
1588 // linear
1589 double a = (new_grid.matrix_array[i] - old_grid.matrix_array[0])/(old_grid.matrix_array[1] - old_grid.matrix_array[0]);
1590 this->matrix_array[i] = old_function.matrix_array[0] + a*(old_function.matrix_array[1] - old_function.matrix_array[0]);
1591 } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-1]) {
1592 this->matrix_array[i] = ub;//1.e-21; //old_function.matrix_array[old_function.size_x-1];
1593 } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-2]) {
1594 // linear
1595 //double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-2])/(old_grid.matrix_array[old_grid.size_x-1] - old_grid.matrix_array[old_grid.size_x-2]);
1596 //this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-2] + a*(old_function.matrix_array[old_grid.size_x-1] - old_function.matrix_array[old_grid.size_x-2]);
1597 double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-1])/(old_grid.matrix_array[old_grid.size_x-2] - old_grid.matrix_array[old_grid.size_x-1]);
1598 this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-1] + a*(old_function.matrix_array[old_grid.size_x-2] - old_function.matrix_array[old_grid.size_x-1]);
1599 } else {
1600 //splint(old_grid.matrix_array, old_function.matrix_array, y2.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i]);
1601 splint(&old_grid.matrix_array[1], &old_function.matrix_array[1], &y2.matrix_array[1], old_grid.size_x-2, new_grid.matrix_array[i], &this->matrix_array[i]);
1602 }
1603 }
1604 
1605 }*/
1606 
1607 /**
1608 * Another spline interpolation. An attempt to make it more stable.
1609 * Runs interpolation functions from Numerical Recepies.
1610 * Treats boundary conditions separete, i.e. use boundary values on baundary conditions and extrapolation, use linear interpolation for next-to-the-boundary points.
1611 * Has additional mechenism of smoothning the interpolation - makes it more linear according the parameters lin_spline_coef and max_second_der.
1612 * Spline interpolation is linear + additional terms from second derivatives (look Numerical Recepies). If second derivative > max_second_der the second derivatives are multiplid by lin_spline_coef.
1613 *
1614 * \param &old_function - old function
1615 * \param &old_grid - old grid
1616 * \param &new_grid - new grid
1617 * \param lb - low boundary value
1618 * \param ub - upper boundary value
1619 * \param lin_spline_coef - coefficient to multiply the second derivatives terms (makes interpolation more smooth)
1620 * \param max_second_der - maximum second derivatives, all greater derivatives gonna be multiplied by lin_spline_coef
1621 *
1622 */
1623 template<class T>
1624 void Matrix1D<T>::Spline2(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub, double lin_spline_coef, double max_second_der) {
1625 
1626  //if (size_x != old_grid.size_x) {
1627  // check, if size of the grids are equal
1628  if (this->size_x != new_grid.size_x) {
1629  throw error_msg("SPLINE_INTERPOLATION", "New grid and function size mismatch!");
1630  }
1631  // if size == 1 - just copy the value
1632  if (this->size_x == 1) {
1633  this->matrix_array[0] = old_function.matrix_array[0];
1634  return;
1635  }
1636 
1637  // Idea was to find the intersection between grids, to interpolate inside the intersection and do not extrapolate outside. But is done later.
1638  int i=0, i_min=0, i_max=old_grid.size_x-1, interp_size=0;
1639  // for(i = 0; old_grid.matrix_array[i] <= new_grid.matrix_array[0] && i < old_grid.size_x-1; i++) {}
1640  // i_min = i;
1641  // for(i = old_grid.size_x-1; old_grid.matrix_array[i] >= new_grid.matrix_array[old_grid.size_x-1] && i > 0; i--) {}
1642  // i_max = i;
1643  //cout << i_min << " " << i_max << endl;
1644 
1645  // new range for spline interpolation - values close to the boundaries are interpolated by linear interpolation
1646  i_min = 0; i_max = old_grid.size_x-1;
1647  interp_size = i_max-i_min+1;
1648 
1649  double yp0, ypn;
1650  //yp0 = 2*(old_function[1] - old_function[0])/(old_grid[1] - old_grid[0]) - (old_function[2] - old_function[1])/(old_grid[2] - old_grid[1]);
1651  //ypn = 2*(old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]) - (old_function[size_x-2] - old_function[size_x-3])/(old_grid[size_x-2] - old_grid[size_x-3]);
1652  // ypn = (old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]);
1653 
1654  // flag, means use smooth second derivative on the boundary
1655  yp0 = 1e99;
1656  //yp0 = -(old_function[2] - old_function[0])/(old_grid[2] - old_grid[0]);
1657  ypn = 1e99;//0;
1658 
1659  Matrix1D<T> y2(old_grid.size_x);
1660  // calculating array of second derivatives
1661  spline(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], interp_size, yp0, ypn, &y2.matrix_array[i_min]);
1662 
1663  // for all grid points - looking for new values
1664  for (i = 0; i < new_grid.size_x; i++) {
1665  splint(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], &y2.matrix_array[i_min], interp_size, new_grid.matrix_array[i], &this->matrix_array[i]);
1666  }
1667 
1668 }
1669 
1670 
1671 /**
1672 * polinomial + lin
1673 * Some other wired interpolation attempts.
1674 */
1675 template<class T>
1676 void Matrix1D<T>::Polilinear(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub) {
1677 
1678  //if (size_x != old_grid.size_x) {
1679  if (this->size_x != new_grid.size_x) {
1680  throw error_msg("POLINOMIAL_INTERPOLATION", "New grid and function size mismatch!");
1681  }
1682 
1683  int i;
1684  for (i = 0; i < new_grid.size_x; i++)
1685  matrix_array[i] = polilinear( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], lb, ub);
1686 
1687 }
1688 
1689 /**
1690 * polinomial
1691 * Some other wired interpolation attempts.
1692 */
1693 template<class T>
1694 void Matrix1D<T>::Polint(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
1695 
1696  //if (size_x != old_grid.size_x) {
1697  if (this->size_x != new_grid.size_x) {
1698  throw error_msg("POLINT_INTERPOLATION", "New grid and function size mismatch!");
1699  }
1700 
1701  int i;
1702  double err;
1703  for (i = 0; i < new_grid.size_x; i++)
1704  polint( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i], &err);
1705 
1706 }
1707 
1708 // ratint
1709 template<class T>
1710 void Matrix1D<T>::Ratint(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
1711 
1712  //if (size_x != old_grid.size_x) {
1713  if (this->size_x != new_grid.size_x) {
1714  throw error_msg("RATIONAL_INTERPOLATION", "New grid and function size mismatch!");
1715  }
1716 
1717  int i;
1718  double err;
1719  for (i = 0; i < new_grid.size_x; i++)
1720  ratint( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i], &err);
1721 
1722 }
1723 
1724 //template<class T>
1725 //Matrix1D<T> Matrix1D<T>::Spline(Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
1726 //
1727 // if (this->size_x != old_grid.size_x) {
1728 // //if (this->size_x != new_grid.size_x) {
1729 // cout << "old grid and function size mismatch!" << endl;
1730 // return -1;
1731 // }
1732 //
1733 // double yp0, ypn;
1734 // yp0 = (matrix_array[0] - matrix_array[1])/(old_grid[0] - old_grid[1]);
1735 // ypn = (matrix_array[size_x-1] - matrix_array[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]);
1736 // Matrix1D<T> y2(old_grid.size_x);
1737 // spline(old_grid.matrix_array, this->matrix_array, old_grid.size_x, yp0, ypn, y2.matrix_array)
1738 //
1739 // int i;
1740 // Matrix1D<T> res(new_grid.size_x);
1741 // for (i = 0; i < new_grid.size_x; i++)
1742 // spline(old_grid.matrix_array, this->matrix_array, y2.matrix_array, new_grid.matrix_array[i], &res[i]);
1743 //
1744 // return res;
1745 //
1746 //}
1747 
1748 
1749 // interpolations
1750 
1751 /**
1752 * Akima interpolation
1753 */
1754 template<class T>
1755 void Matrix1D<T>::Akima_interpolation(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, int extrapolation_type, double lb, double ub) {
1756 
1757  /*
1758  for (i = 0; i < new_grid.size_x; i++) {
1759  if (i == 0) { // if we are on the boundary
1760  this->matrix_array[i] = lb;
1761  } else if (i == new_grid.size_x-1) { // if we are on the boundary
1762  this->matrix_array[i] = ub;
1763  } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[0] ) { // if it is extrapolation
1764  this->matrix_array[i] = lb;
1765  //this->matrix_array[i] = old_function.matrix_array[0];
1766  } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-1]) { // if it is extrapolation
1767  this->matrix_array[i] = ub;
1768  //this->matrix_array[i] = old_function.matrix_array[old_function.size_x-1];
1769  } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[1] ) { // if the new grid point is between two boundaryes poins of the old grid
1770  // linear
1771  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[0])/(old_grid.matrix_array[1] - old_grid.matrix_array[0]);
1772  this->matrix_array[i] = old_function.matrix_array[0] + a*(old_function.matrix_array[1] - old_function.matrix_array[0]);
1773  } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-2]) { // if the new grid point is between two boundaryes poins of the old grid
1774  // linear
1775  //double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-2])/(old_grid.matrix_array[old_grid.size_x-1] - old_grid.matrix_array[old_grid.size_x-2]);
1776  //this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-2] + a*(old_function.matrix_array[old_grid.size_x-1] - old_function.matrix_array[old_grid.size_x-2]);
1777  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-1])/(old_grid.matrix_array[old_grid.size_x-2] - old_grid.matrix_array[old_grid.size_x-1]);
1778  this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-1] + a*(old_function.matrix_array[old_grid.size_x-2] - old_function.matrix_array[old_grid.size_x-1]);
1779  } else if (i <= 1) { // always linear for the first point?
1780  int k = 1;
1781  // Wrong! We need to compare to grids, not indexes. Changed in ver 1.04.
1782  // while(k < i) {k++;}
1783  while(old_grid.matrix_array[k] < new_grid.matrix_array[i]) {k++;}
1784  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[k-1])/(old_grid.matrix_array[k] - old_grid.matrix_array[k-1]);
1785  this->matrix_array[i] = old_function.matrix_array[k-1] + a*(old_function.matrix_array[k] - old_function.matrix_array[k-1]);
1786  } else {
1787  akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, ub, lb);
1788  }
1789  }*/
1790 
1791 
1792  // add point
1793  /*int i_first = 0, i_last = new_grid.size_x-1;
1794  while(old_grid.matrix_array[i_first] <= new_grid.matrix_array[0]) {i_first++;}
1795  while(old_grid.matrix_array[i_last] >= new_grid.matrix_array[new_grid.size_x-1]) {i_last--;}
1796  Matrix1D<double> mod_grid(i_last-i_first+2);
1797  mod_grid[0] = lb;
1798  mod_grid[mod_grid.x_size-1] = ub;
1799 
1800  for (int i = 1; i < mod_grid.size_x-1; i++) {
1801  mod_grid[i] = old_grid.matrix_array[i_first + i - 1];
1802  }
1803 
1804  akima(mod_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, 1, ub, lb);*/
1805 
1806 
1807  /*int i_first = 0, i_last = new_grid.size_x-1;
1808  while(new_grid.matrix_array[i_first] <= old_grid.matrix_array[0]) {i_first++;}
1809  while(new_grid.matrix_array[i_last] >= old_grid.matrix_array[old_grid.size_x-1]) {i_last--;}
1810  int new_size_x = i_last - i_first + 1;
1811 
1812  akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x,
1813  &new_grid.matrix_array[i_first], &this->matrix_array[i_first], new_size_x,
1814  2,
1815  ub, lb);
1816 
1817  // boundaries (and extrapolation)
1818  int i;
1819  for (i = 0; i < i_first; i++) this->matrix_array[i] = lb;
1820  for (i = this->size_x-1; i > i_last; i--) this->matrix_array[i] = ub;*/
1821 
1822  // new grid is exactly one point wider than old grid
1823  /*int old_first = 0, old_last = old_grid.size_x-1;
1824  while(old_grid.matrix_array[old_first] <= new_grid.matrix_array[0]) {old_first++;}
1825  while(old_grid.matrix_array[old_last] >= new_grid.matrix_array[new_grid.size_x-4]) {old_last--;}
1826  int old_size_x = old_last - old_first + 1;
1827  int new_first = 0, new_last = new_grid.size_x-1;
1828  while(new_grid.matrix_array[new_first+1] <= old_grid.matrix_array[0]) {new_first++;}
1829  while(new_grid.matrix_array[new_last-1] >= old_grid.matrix_array[old_grid.size_x-4]) {new_last--;}
1830  int new_size_x = new_last - new_first + 1;
1831 
1832  akima(&old_grid.matrix_array[old_first], &old_function.matrix_array[old_first], old_size_x,
1833  &new_grid.matrix_array[new_first], &this->matrix_array[new_first], new_size_x, 1, ub, lb);
1834 
1835  // boundaries (and extrapolation)
1836  int i;
1837  for (i = 0; i < new_first; i++) this->matrix_array[i] = lb;
1838  for (i = this->size_x-1; i > new_last; i--) this->matrix_array[i] = ub;*/
1839 
1840  ///
1841  /* int i_first = 0, i_last = new_grid.size_x-1;
1842  while(old_grid.matrix_array[i_first] <= new_grid.matrix_array[0]) {i_first++;}
1843  while(old_grid.matrix_array[i_last] >= new_grid.matrix_array[new_grid.size_x-1]) {i_last--;}
1844  int new_size_x = i_last - i_first + 1;
1845 
1846  //akima(&old_grid.matrix_array[i_first], &old_function.matrix_array[i_first], new_size_x, new_grid.matrix_array, this->matrix_array, this->size_x, 1, ub, lb);
1847  // akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, 1, ub, lb);
1848  akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x,
1849  &new_grid.matrix_array[i_first], &this->matrix_array[i_first], new_size_x,
1850  2,
1851  ub, lb,
1852  max(this->matrix_array[0], old_grid.matrix_array[0]), min(this->matrix_array[this->x_size-1], old_grid.matrix_array[old_grid.x_size-1]));
1853 
1854  // boundaries (and extrapolation)
1855  int i;
1856  for (i = 0; i < i_first; i++) this->matrix_array[i] = lb;
1857  for (i = this->size_x-1; i > i_last; i--) this->matrix_array[i] = ub;*/
1858 
1859 
1860 
1861  akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, 1, ub, lb);
1862  matrix_array[0] = lb;
1863  matrix_array[size_x-1] = ub;
1864 }
1865 
1866 
1867 /**
1868 * Akima interpolation - second test version
1869 */
1870 template<class T>
1871 void Matrix1D<T>::Akima_interpolation2(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, int extrapolation_type, double lb, double ub) {
1872  akima2(old_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, ub, lb);
1873 }
1874 
1875 /**
1876 * Linear2D. Does not work probably or works only for regular grid.
1877 * Some other wired interpolation attempts.
1878 */
1879 template<class T>
1880 double Linear2D(double x, double y, T &old_grid_x, T &old_grid_y, T &f) {
1881 
1882  int i = 0, j = 0;
1883 
1884  while (x >= old_grid_x[++i][0] && i < old_grid_y.size_y-1);
1885  while (y >= old_grid_y[0][++j] && j < old_grid_x.size_x-1);
1886 
1887  double t = (x - old_grid_x[i-1][0])/(old_grid_x[i][0] - old_grid_x[i-1][0]);
1888  double u = (y - old_grid_y[0][j-1])/(old_grid_y[0][j] - old_grid_y[0][j-1]);
1889 
1890  return (1.0-t)*(1.0-u)*f[i-1][j-1] + t*(1.0-u)*f[i][j-1] + t*u*f[i][j] + (1.0-t)*u*f[i-1][j];
1891 }
1892 
1893 
1894 /**
1895 * Linear 2D interpolation.
1896 *
1897 * Not sure if it is working.
1898 */
1899 template<class T>
1900 void Matrix2D<T>::Interpolate(Matrix2D<T> &old_function, Matrix2D<T> &old_grid_x, Matrix2D<T> &old_grid_y, Matrix2D<T> &new_grid_x, Matrix2D<T> &new_grid_y) {
1901 
1902  Matrix2D<T> res(new_grid_x.size_x, new_grid_y.size_y);
1903 
1904  if (old_function.size_x != old_grid_x.size_x || old_function.size_y != old_grid_y.size_y) {
1905  throw error_msg("LINEAR_INTERPOLATION", "Old grid and function size mismatch!");
1906  }
1907 
1908  if (this->size_x != new_grid_x.size_x || this->size_y != new_grid_y.size_y) {
1909  throw error_msg("LINEAR_INTERPOLATION", "New grid and function size mismatch!");
1910  }
1911 
1912  int i, j;
1913  for (i = 0; i < new_grid_x.size_x; i++) {
1914  for (j = 0; j < new_grid_y.size_y; j++) {
1915  this->matrix_array[i][j] = Linear2D(new_grid_x[i][j], new_grid_y[i][j], old_grid_x, old_grid_y, old_function);
1916  }
1917  }
1918 
1919  return;
1920 
1921 }
1922 
1923 
1924 /**
1925 * Constructor for CalculationMatrix class.
1926 * \param x_size - x size
1927 * \param y_size - y size
1928 * \param z_size - z size
1929 * \param n_of_diags - NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL (main diagonal is not counted)
1930 */
1931 CalculationMatrix::CalculationMatrix(int x_size, int y_size, int z_size, int n_of_diags) {
1932  this->initialized = false;
1933  Initialize(x_size, y_size, z_size, n_of_diags);
1934 }
1935 
1936 /**
1937 * Allocating memory for CalculationMatrix
1938 */
1939 void CalculationMatrix::Initialize(int x_size, int y_size, int z_size, int n_of_diags) {
1940  this->initialized = false;
1941  this->x_size = x_size;
1942  this->y_size = y_size;
1943  this->total_size = x_size;
1944  if (y_size > 0) this->total_size = this->total_size * y_size;
1945  if (z_size > 0) this->total_size = this->total_size * z_size;
1946 
1947 
1948  // initializing diagonals
1949  int x, y, z, id;
1950  // !!! if (z_size > 0) { // means 3d
1951  if (z_size > 1) { // means 3d
1952  for (x = -n_of_diags; x <= n_of_diags; x++) {
1953  for (y = -n_of_diags; y <= n_of_diags; y++) {
1954  for (z = -n_of_diags; z <= n_of_diags; z++) {
1955  // calculating a diagonal number (id)
1956  id = this->index1d(x, y, z);
1957  // allocating memory for each diagonal
1958  (*this)[id] = Matrix1D<double>(this->total_size);
1959  (*this)[id] = 0;
1960  }
1961  }
1962  }
1963  // !!! } else if (y_size > 0) {
1964  } else if (y_size > 1) {
1965  for (x = -n_of_diags; x <= n_of_diags; x++) {
1966  for (y = -n_of_diags; y <= n_of_diags; y++) {
1967  // calculating a diagonal number (id)
1968  id = this->index1d(x, y);
1969  // allocating memory for each diagonal
1970  (*this)[id] = Matrix1D<double>(this->total_size);
1971  (*this)[id] = 0;
1972  }
1973  }
1974  // !!! } else if (x_size > 0) {
1975  } else if (x_size > 1) {
1976  for (x = -n_of_diags; x <= n_of_diags; x++) {
1977  // calculating a diagonal number (id)
1978  id = this->index1d(x);
1979  // allocating memory for each diagonal
1980  (*this)[id] = Matrix1D<double>(this->total_size);
1981  (*this)[id] = 0;
1982  }
1983  }
1984 
1985  // If im_size is not set, then it's 1d and ia_size should be zero for future
1986  if (y_size == 0) this->x_size = 0;
1987  // If il_size is not set, then it's 2d and im_size should be zero for future
1988  if (z_size == 0) this->y_size = 0;
1989 
1990  this->initialized = true;
1991 }
1992 
1993 /**
1994 * Function returns 1d index for 2D-3D arrays
1995 */
1996 int CalculationMatrix::index1d(int x, int y, int z) {
1997  return (z*y_size + y)*x_size + x;
1998 }
1999 
2000 /**
2001 * Save matrix to file.
2002 */
2003 void CalculationMatrix::writeToFile(string filename) {
2004  int in;
2005 
2006  ofstream output(filename.c_str());
2007  if (output==NULL && (filename.find("Debug") == string::npos)) {
2008  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
2009  }
2010  output << "VARIABLES = ";
2011  for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
2012  output << "\"" << it->first << "\", ";
2013  }
2014  output << endl;
2015  output << "ZONE T=\"" << "\", I=" << total_size << endl;
2016  output.setf(ios_base::scientific, ios_base::floatfield);
2017  for (in = 0; in < this->total_size; in++) {
2018  for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
2019  output << "\t" << it->second[in];
2020  }
2021  output << endl;
2022  }
2023  output.close();
2024 }
2025 
2026 
2027 //////////////////////////////////////////
2028 // Implementations
2029 //////////////////////////////////////////
2030 
2031 template class Matrix1D<double>;
2032 template class Matrix2D<double>;
2033 template class Matrix3D<double>;
2034 
2035 
2036 #endif