USML
boundary_grid.h
1 
6 #ifndef USML_OCEAN_BOUNDARY_GRID_H
7 #define USML_OCEAN_BOUNDARY_GRID_H
8 
9 #include <usml/ocean/boundary_model.h>
10 #include <usml/ocean/reflect_loss_rayleigh.h>
11 
12 namespace usml {
13 namespace ocean {
14 
32 template< class DATA_TYPE, int NUM_DIMS > class boundary_grid
33  : public boundary_model
34 {
35  //**************************************************
36  // height model
37 
38 protected:
39 
42 
43 public:
44 
54  virtual void height(const wposition& location, matrix<double>* rho,
55  wvector* normal = NULL, bool quick_interp = false) {
56  switch (NUM_DIMS) {
57 
58  //***************
59  // 1-D grids
60 
61  case 1:
62  if(quick_interp) { this->_height->interp_type(0, GRID_INTERP_LINEAR); }
63  else{ this->_height->interp_type(0, GRID_INTERP_PCHIP); }
64  if (normal) {
65 
66  matrix<double> gtheta(location.size1(), location.size2());
67  matrix<double> t(location.size1(), location.size2());
68  this->_height->interpolate(location.theta(), rho, &gtheta);
69  t = min(element_div(gtheta,*rho),1.0); // slope = tan(angle)
70  normal->theta( // normal = -sin(angle)
71  element_div( -t, sqrt(1.0+abs2(t)) ));
72  normal->phi(scalar_matrix<double>(location.size1(),location.size2(),0.0));
73  normal->rho( sqrt(1.0-abs2(normal->theta())) ) ; // r=sqrt(1-t^2)
74  } else {
75  this->_height->interpolate(location.theta(), rho);
76  }
77  break;
78 
79  //***************
80  // 2-D grids
81 
82  case 2:
83  if(quick_interp) {
84  this->_height->interp_type(0, GRID_INTERP_LINEAR);
85  this->_height->interp_type(1, GRID_INTERP_LINEAR);
86  } else{
87  this->_height->interp_type(0, GRID_INTERP_PCHIP);
88  this->_height->interp_type(1, GRID_INTERP_PCHIP);
89  }
90  if (normal) {
91  matrix<double> gtheta(location.size1(), location.size2());
92  matrix<double> gphi(location.size1(), location.size2());
93  matrix<double> t(location.size1(), location.size2());
94  matrix<double> p(location.size1(), location.size2());
95  this->_height->interpolate(location.theta(), location.phi(),
96  rho, &gtheta, &gphi);
97 
98  t = element_div(gtheta, *rho); // slope = tan(angle)
99  p = element_div(gphi, element_prod(*rho, sin(location.theta())));
100  normal->theta( // normal = -sin(angle)
101  element_div( -t, sqrt(1.0+abs2(t)) ));
102  normal->phi(
103  element_div( -p, sqrt(1.0+abs2(p)) ));
104  normal->rho(sqrt( // r=sqrt(1-t^2-p^2)
105  1.0 - abs2(normal->theta()) - abs2(normal->phi()) ));
106  } else {
107  this->_height->interpolate(location.theta(), location.phi(),
108  rho);
109  }
110  break;
111 
112  //***************
113  // error
114 
115  default:
116  throw std::invalid_argument("dataset must be 1-D or 2-D");
117  break;
118  }
119  }
120 
130  virtual void height(const wposition1& location, double* rho,
131  wvector1* normal = NULL, bool quick_interp = false) {
132  switch (NUM_DIMS) {
133 
134  //***************
135  // 1-D grids
136 
137  case 1:
138  if(quick_interp) { this->_height->interp_type(0, GRID_INTERP_LINEAR); }
139  else{ this->_height->interp_type(0, GRID_INTERP_PCHIP); }
140  if (normal) {
141  double theta = location.theta();
142  double gtheta;
143  *rho = this->_height->interpolate(&theta, &gtheta);
144  const double t = gtheta / (*rho); // slope = tan(angle)
145  normal->theta(-t / sqrt(1.0 + t * t)); // normal = -sin(angle)
146  normal->phi(0.0);
147  const double N = normal->theta() * normal->theta() ;
148  normal->rho( sqrt(1.0-N) ); // r=sqrt(1-t^2)
149  } else {
150  double theta = location.theta();
151  *rho = this->_height->interpolate(&theta);
152  }
153  break;
154 
155  //***************
156  // 2-D grids
157 
158  case 2:
159  if(quick_interp) {
160  this->_height->interp_type(0, GRID_INTERP_LINEAR);
161  this->_height->interp_type(1, GRID_INTERP_LINEAR);
162  } else{
163  this->_height->interp_type(0, GRID_INTERP_PCHIP);
164  this->_height->interp_type(1, GRID_INTERP_PCHIP);
165  }
166  if (normal) {
167  double loc[2] = { location.theta(), location.phi() };
168  double grad[2];
169  *rho = this->_height->interpolate(loc, grad);
170  const double t = grad[0] / (*rho); // slope = tan(angle)
171  const double p = grad[1] / ((*rho) * sin(location.theta()));
172  normal->theta(-t / sqrt(1.0 + t * t)); // normal = -sin(angle)
173  normal->phi(-p / sqrt(1.0 + p * p));
174  const double N = normal->theta() * normal->theta()
175  + normal->phi() * normal->phi();
176  normal->rho( sqrt(1.0-N) ); // r=sqrt(1-t^2-p^2)
177  } else {
178  double loc[2] = { location.theta(), location.phi() };
179  *rho = this->_height->interpolate(loc);
180  }
181  break;
182 
183  //***************
184  // error
185 
186  default:
187  throw std::invalid_argument("bathymetry must be 1-D or 2-D");
188  break;
189  }
190  }
191 
192  //**************************************************
193  // initialization
194 
209  this->_height->interp_type(0,GRID_INTERP_PCHIP);
210  this->_height->interp_type(1,GRID_INTERP_PCHIP);
211  this->_height->edge_limit(0,true);
212  this->_height->edge_limit(1,true);
213  if (_reflect_loss_model == NULL) {
216  }
217  }
218 
222  virtual ~boundary_grid() {
223  delete _height;
224  }
225 
226 };
227 
228 } // end of namespace ocean
229 } // end of namespace usml
230 
231 #endif
unsigned size2() const
Number of columns in each coordinate.
Definition: wvector.h:264
Definition: reflect_loss_rayleigh.h:96
World location in geodetic earth coordinates (latitude, longitude, and altitude). ...
Definition: wposition.h:39
BOOST_UBLAS_INLINE matrix_unary1_traits< E, scalar_abs2< typename E::value_type > >::result_type abs2(const matrix_expression< E > &e)
Magnitude squared of a complex matrix.
Definition: matrix_math.h:270
unsigned size1() const
Number of rows in each coordinate.
Definition: wvector.h:256
Definition: data_grid.h:23
BOOST_UBLAS_INLINE boost::enable_if< boost::is_convertible< T2, typename E1::value_type >, typename matrix_binary_scalar2_traits< E1, const T2, scalar_min< typename E1::value_type, T2 > >::result_type >::type min(const matrix_expression< E1 > &e1, const T2 &e2)
Minimum between a matrix and a scalar.
Definition: matrix_math.h:158
A reflection loss model computes the changes in amplitude and phase that result from the reflection o...
Definition: reflect_loss_model.h:30
Bottom model constructed from a 1-D or 2-D data grid.
Definition: boundary_grid.h:32
Individual world vector in spherical earth coordinates.
Definition: wvector1.h:27
double theta() const
Retrieves the colatitude component of the spherical earth coordinate system.
Definition: wvector1.h:91
virtual void height(const wposition &location, matrix< double > *rho, wvector *normal=NULL, bool quick_interp=false)
Compute the height of the boundary and it's surface normal at a series of locations.
Definition: boundary_grid.h:54
reflect_loss_model * _reflect_loss_model
Reference to the reflection loss model.
Definition: boundary_model.h:95
double phi() const
Retrieves the longitude component of the spherical earth coordinate system.
Definition: wvector1.h:122
data_grid< DATA_TYPE, NUM_DIMS > * _height
Boundary for all locations.
Definition: boundary_grid.h:41
void reflect_loss(reflect_loss_model *reflect_loss)
Define a new reflection loss model.
Definition: boundary_model.h:104
const matrix< double > & phi() const
Retrieves the longitude component of the spherical earth coordinate system.
Definition: wvector.h:204
const matrix< double > & theta() const
Retrieves the colatitude component of the spherical earth coordinate system.
Definition: wvector.h:140
BOOST_UBLAS_INLINE matrix_unary1_traits< E, scalar_sin< typename E::value_type > >::result_type sin(const matrix_expression< E > &e)
Sine of a matrix.
Definition: matrix_math.h:354
World location in geodetic earth coordinates (latitude, longitude, and altitude). ...
Definition: wposition1.h:23
N-dimensional data set and its associated axes.
Definition: data_grid.h:64
virtual ~boundary_grid()
Delete boundary grid.
Definition: boundary_grid.h:222
Models plane wave reflection loss from a flat fluid-solid interface.
Definition: reflect_loss_rayleigh.h:86
boundary_grid(data_grid< DATA_TYPE, NUM_DIMS > *height, reflect_loss_model *reflect_loss=NULL)
Initialize depth and reflection loss components for a boundary.
Definition: boundary_grid.h:206
Definition: data_grid.h:24
World vector in spherical earth coordinates.
Definition: wvector.h:42
A "boundary model" computes the environmental parameters of the ocean's surface or bottom...
Definition: boundary_model.h:58
BOOST_UBLAS_INLINE matrix_unary1_traits< E, scalar_sqrt< typename E::value_type > >::result_type sqrt(const matrix_expression< E > &e)
Square root of a matrix.
Definition: matrix_math.h:296
virtual void height(const wposition1 &location, double *rho, wvector1 *normal=NULL, bool quick_interp=false)
Compute the height of the boundary and it's surface normal at a single location.
Definition: boundary_grid.h:130