RichardsFastSlow/backup/obsolete/geometry_spline.cpp

#include "geometry_spline.hpp"

void
GeometrySpline::initDefault(){
  lX=10;
  nX=100;
  nX_max=nmax_Qt;
  dX=lX/(nX-1);
  hsoil=new double[nmax_Qt];
  dhsoil=new double[nmax_Qt];
  hbot=new double[nmax_Qt];
  dhbot=new double[nmax_Qt];
  nZ=new size_t[nmax_Qt];
  dZ=new double[nmax_Qt];
  Z=new double*[nmax_Qt];
  for(size_t i=0;i<nmax_Qt;++i){
    Z[i]=new double[nmax_Qt];
  }
  for(size_t i=0;i<nX;++i){
    hsoil[i]=0.8;
    hbot[i]=0.2;
    dhsoil[i]=0;
    dhbot[i]=0;
  }
  depth=5;
  nZ_max=100;
  factor=5/0.6;
  initZ(depth/nZ_max,false);
}

void
GeometrySpline::update(double _lX,size_t _nX,double _depth,size_t _nZ_max,Spline& Ssoil,Spline& Sbot){
  lX=_lX;
  depth=_depth;
  nX=_nX;
  nZ_max=_nZ_max;
  dX=lX/(nX-1);
  double hs_max=-inf,hb_min=inf;
  double v;
  for(size_t k=0;k<nX;++k){
    double x=k*dX/lX;
    v=hsoil[k]=Ssoil(x);
    hs_max=max(v,hs_max);
    dhsoil[k]=Ssoil.derivate(x);
    v=hbot[k]=Sbot(x);
    hb_min=min(v,hb_min);
    dhbot[k]=Sbot.derivate(x);
  }
  double dZ_avg=depth/_nZ_max;
  factor=depth/(hs_max-hb_min);
  for(size_t k=0;k<nX;++k){
    hsoil[k]*=factor;
    hbot[k]*=factor;
  }
  initZ(dZ_avg,false);
}

GeometrySpline::~GeometrySpline(){
  if(hsoil!=nullptr){
    delete[] hsoil;
    delete[] dhsoil;
    delete[] hbot;
    delete[] dhbot;
    delete[] dZ;
    delete[] nZ;
    for(size_t i=0;i<nX_max;++i){
      delete[] Z[i];
    }
    delete[] Z;
  }
  hsoil=nullptr;
}