00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00022
00024
00044
00045
00047
00048 #include <time.h>
00049 #include <math.h>
00050 #include <stdexcept>
00051 #include "arts.h"
00052 #include "math_funcs.h"
00053 #include "make_vector.h"
00054 #include "array.h"
00055
00056
00057
00061
00070 Numeric first( ConstVectorView x )
00071 {
00072 return x[0];
00073 }
00074
00083 Numeric last( ConstVectorView x )
00084 {
00085 return x[x.nelem()-1];
00086 }
00087
00088
00089
00091
00093
00095
00104 bool any( const ArrayOfIndex& x )
00105 {
00106 for ( Index i=0; i<x.nelem(); i++ ) {
00107 if ( x[i] )
00108 return true;
00109 }
00110 return false;
00111 }
00112
00113
00114
00115
00117
00119
00121
00139 void linspace(
00140 Vector& x,
00141 const Numeric start,
00142 const Numeric stop,
00143 const Numeric step )
00144 {
00145 Index n = (Index) floor( (stop-start)/step ) + 1;
00146 if ( n<1 )
00147 n=1;
00148 x.resize(n);
00149 for ( Index i=0; i<n; i++ )
00150 x[i] = start + i*step;
00151 }
00152
00168 void nlinspace(
00169 Vector& x,
00170 const Numeric start,
00171 const Numeric stop,
00172 const Index n )
00173 {
00174 assert( 1<n );
00175 x.resize(n);
00176 Numeric step = (stop-start)/(n-1) ;
00177 for ( Index i=0; i<n; i++ )
00178 x[i] = start + i*step;
00179 }
00180
00182
00198 void nlogspace(
00199 Vector& x,
00200 const Numeric start,
00201 const Numeric stop,
00202 const Index n )
00203 {
00204
00205 assert( 1<n );
00206
00207 assert( 0<start );
00208 assert( 0<stop );
00209
00210 x.resize(n);
00211 Numeric a = log(start);
00212 Numeric step = (log(stop)-a)/(n-1);
00213 x[0] = start;
00214 for ( Index i=1; i<n-1; i++ )
00215 x[i] = exp(a + i*step);
00216 x[n-1] = stop;
00217 }
00218
00234 Vector nlogspace(
00235 const Numeric start,
00236 const Numeric stop,
00237 const Index n )
00238 {
00239 Vector x;
00240 nlogspace( x, start, stop, n );
00241 return x;
00242 }
00243
00244
00245
00246
00248
00250
00259 Index interp_check( ConstVectorView x,
00260 ConstVectorView xi,
00261 const Index n_y )
00262 {
00263 const Index n = x.nelem();
00264 const Index ni = xi.nelem();
00265 Index order=1;
00266
00267
00268 if ( x[0] > x[n-1] )
00269 order = -1;
00270
00271 if ( n < 2 )
00272 throw runtime_error("Vector length for interpolation must be >= 2.");
00273
00274 if ( n != n_y )
00275 throw runtime_error("Sizes of input data to interpolation do not match.");
00276
00277
00278
00279 const Numeric lower_bound = x[0] - 0.5*(x[1]-x[0]);
00280 const Numeric upper_bound = x[n-1] + 0.5*(x[n-1]-x[n-2]);
00281
00282 if ( (order*xi[0]<order*lower_bound) || (order*xi[ni-1]>order*upper_bound) )
00283 {
00284 ostringstream os;
00285 os << "Interpolation points must be not more than\n"
00286 << "half a grid spacing outside the original range.\n"
00287 << "Int.: xi[0] = " << xi[0] << ", xi[ni-1] = " << xi[ni-1] << '\n'
00288 << "Orig.: x[0] = " << x[0] << ", x[n-1] = " << x[n-1];
00289 throw runtime_error(os.str());
00290 }
00291
00292 for ( Index i=0; i<n-1; i++ )
00293 {
00294 if ( order*x[i+1] < order*x[i] )
00295 throw runtime_error("Original interpolation grid must be ordered");
00296 }
00297
00298 for ( Index i=0; i<ni-1; i++ )
00299 {
00300 if ( order*xi[i+1] < order*xi[i] )
00301 throw runtime_error("Interpolation points must be ordered");
00302 }
00303
00304 return order;
00305 }
00306
00328 void interp_lin_vector( VectorView yi,
00329 ConstVectorView x,
00330 ConstVectorView y,
00331 ConstVectorView xi )
00332 {
00333
00334 Index order = interp_check( x, xi, y.nelem() );
00335
00336 Index i, j=0;
00337 const Index n=xi.nelem();
00338 const Index nx=x.nelem();
00339 Numeric w;
00340
00341
00342 assert( n==yi.nelem() );
00343
00344 for ( i=0; i<n; i++ )
00345 {
00346
00347 while ( j<nx-2 && order*x[j+1]<order*xi[i] )
00348 {
00349 ++j;
00350 }
00351
00352
00353
00354
00355
00356 assert( x[j+1]!=x[j] );
00357 w = (xi[i]-x[j]) / (x[j+1]-x[j]);
00358
00359
00360
00361 yi[i] = y[j] + w * (y[j+1]-y[j]);
00362 }
00363 }
00364
00382 void interp_lin_matrix(
00383 MatrixView Yi,
00384 ConstVectorView x,
00385 ConstMatrixView Y,
00386 ConstVectorView xi )
00387 {
00388
00389 Index order = interp_check( x, xi, Y.ncols() );
00390
00391 Index j=0;
00392 const Index n=xi.nelem(), nrow=Y.nrows(), nx=x.nelem();
00393 Numeric w;
00394
00395 assert( nrow == Yi.nrows() );
00396 assert( n == Yi.ncols() );
00397
00398 for (Index i=0; i<n; i++ )
00399 {
00400
00401 while ( j<nx-2 && order*x[j+1]<order*xi[i] )
00402 {
00403 ++j;
00404 }
00405
00406
00407
00408
00409
00410 assert( x[j+1]!=x[j] );
00411 w = (xi[i]-x[j]) / (x[j+1]-x[j]);
00412
00413
00414
00415 for( Index k=0; k<nrow; k++ )
00416 {
00417 Yi(k,i) = Y(k,j) + w * (Y(k,j+1)-Y(k,j));
00418 }
00419 }
00420 }
00421
00434 Numeric interp_lin(
00435 ConstVectorView x,
00436 ConstVectorView y,
00437 const Numeric xi )
00438 {
00439 Vector Yi(1);
00440 MakeVector Xi(xi);
00441
00442
00443
00444 interp_lin_vector( Yi, x, y, Xi );
00445 return Yi[0];
00446 }
00447
00448
00449
00450
00452
00454
00456
00467 void check_if_bool( const Index& x, const String& x_name )
00468 {
00469 if ( !(x==0 || x==1) )
00470 {
00471 ostringstream os;
00472 os << "The boolean *" << x_name << "* must either be 0 or 1.\n"
00473 << "The present value of *"<< x_name << "* is " << x << ".";
00474 throw runtime_error( os.str() );
00475 }
00476 }
00477
00478
00479
00481
00494 void check_if_in_range(
00495 const Numeric& x_low,
00496 const Numeric& x_high,
00497 const Numeric& x,
00498 const String& x_name )
00499 {
00500 if ( (x<x_low) || (x>x_high) )
00501 {
00502 ostringstream os;
00503 os << "The variable *" << x_name << "* must fulfill:\n"
00504 << " " << x_low << " <= " << x_name << " <= " << x_high << "\n"
00505 << "The present value of *"<< x_name << "* is " << x << ".";
00506 throw runtime_error( os.str() );
00507 }
00508 }
00509
00510
00511
00513
00526 void check_lengths( const Vector& x1, const String& x1_name,
00527 const Vector& x2, const String& x2_name )
00528 {
00529 if ( x1.nelem() != x2.nelem() )
00530 {
00531 ostringstream os;
00532 os << "The vectors *" << x1_name << "* and *" << x2_name << "*\n"
00533 << "must have the same lengths. \nThe lengths are: \n"
00534 << x1_name << ": " << x1.nelem() << "\n"
00535 << x2_name << ": " << x2.nelem();
00536 throw runtime_error( os.str() );
00537 }
00538 }
00539
00540
00541
00543
00558 void check_length_nrow( const Vector& x, const String& x_name,
00559 const Matrix& A, const String& A_name )
00560 {
00561 if ( x.nelem() != A.nrows() )
00562 {
00563 ostringstream os;
00564 os << "The length of vector *" << x_name << "* must be the\n"
00565 << "same as the number of rows of *" << A_name << "*.\n"
00566 << "The length of *" << x_name << "* is " << x.nelem() << ".\n"
00567 << "The number of rows of *" << A_name << "* is " << A.nrows()
00568 << ".";
00569 throw runtime_error( os.str() );
00570 }
00571 }
00572
00573
00574
00576
00591 void check_length_ncol( const Vector& x, const String& x_name,
00592 const Matrix& A, const String& A_name )
00593 {
00594 if ( x.nelem() != A.ncols() )
00595 {
00596 ostringstream os;
00597 os << "The length of vector *" << x_name << "* must be the\n"
00598 << "same as the number of columns of *" << A_name << "*.\n"
00599 << "The length of *" << x_name << "* is " << x.nelem() << ".\n"
00600 << "The number of columns of *" << A_name << "* is " << A.ncols()
00601 << ".";
00602 throw runtime_error( os.str() );
00603 }
00604 }
00605
00607
00621 void check_ncol_nrow( const Matrix& A, const String& A_name,
00622 const Matrix& B, const String& B_name )
00623 {
00624 if ( A.ncols() != B.nrows() )
00625 {
00626 ostringstream os;
00627 os << "The number of columns of *" << A_name << "* must be the\n"
00628 << "same as the number of rows of *" << B_name << "*."
00629 << "The number of columns of *" << A_name << "* is " << A.ncols()
00630 << ".\n"
00631 << "The number of rows of *" << B_name << "* is " << B.nrows()
00632 << ".";
00633 throw runtime_error( os.str() );
00634 }
00635 }