RichardsFastSlow/src/geometry.cpp

#include "geometry.hpp"

double inf = numeric_limits<double>::infinity();

Geometry::Geometry(double _lX,size_t _nX,size_t _nZ,Func _hsoil,Func _dhsoil,Func _hbot,Func _dhbot){
  lX=_lX;
  nX=_nX;
  dX=lX/nX;
  hsoil=new double[nX];
  dhsoil=new double[nX];
  hbot=new double[nX];
  dhbot=new double[nX];
  double hs_max=-inf,hb_min=inf;
  double v;
  for(size_t k=0;k<nX;++k){
    double x=k*dX;
    v=hsoil[k]=_hsoil(x);
    hs_max=max(v,hs_max);
    dhsoil[k]=_dhsoil(x);
    v=hbot[k]=_hbot(x);
    hb_min=min(v,hb_min);
    dhbot[k]=_dhbot(x);
  }
  double dZ_avg=(hs_max-hb_min)/_nZ;
  nZ=new size_t[nX];
  dZ=new double[nX];
  Z=new double*[nX];
  for(size_t k=0;k<nX;++k){
    double d=hsoil[k]-hbot[k];
    size_t n=d/dZ_avg;
    double dz=d/n;
    dZ[k]=dz;
    nZ[k]=n;
    Z[k]=new double[n];
    for(size_t j=0;j<n;++j){
      Z[k][j]=hbot[k]+j*dz;
    }
  }
}