$search

uspline.cc

Go to the documentation of this file.

/*
  James R. Diebel
  Stanford University

  Started: 22 November 2004
  Last revised: 2004

  uspline.cc - implements class defined in uspline.hh

  Depends on:
  - USpline Class (uspline.hh)
  - Matrix inversion routine (matrix.hh)
*/

#include <string.h>
#include <iostream>
#include "bmtk/uspline.hh"

using namespace std;

namespace bmtk {

USpline::USpline(int n_) {
  n = n_;
  b0 = 0; bn = 0;
  xmin = 0; xmax = 1; span = xmax - xmin;
  interval = span/float(n-1);
  dpdx = float(n-1)/span;
  a = new float[4*n];
  y = a; b = a+n; c = b+n; d = c+n;
  A = new float *[n];
  for (int i=0;i<n;i++) A[i] = new float[n];
}

USpline::USpline(int n_,float xmin_, float xmax_, float b0_, float bn_) {
  n = n_;
  b0 = b0_; bn = bn_;
  xmin = xmin_; xmax = xmax_; span = xmax - xmin;
  interval = span/float(n-1);
  dpdx = float(n-1)/span;
  a = new float[4*n];
  y = a; b = a+n; c = b+n; d = c+n;
  A = new float *[n];
  for (int i=0;i<n;i++) A[i] = new float[n];
}

USpline::USpline(float* y_, int n_) {
  n = n_;
  b0 = 0; bn = 0;
  xmin = 0; xmax = 1; span = xmax - xmin;
  interval = span/float(n-1);
  dpdx = float(n-1)/span;
  a = new float[n];
  y = a; b = a+n; c = b+n; d = c+n;
  A = new float *[n];
  for (int i=0;i<n;i++) A[i] = new float[n];

  for (int i=0;i<n;i++) {
    y[i] = y_[i];
  }
}

USpline::~USpline() {
  delete [] a;
  for (int i=0;i<n;i++) delete [] A[i];
  delete [] A;
}


void USpline::update() {
  for (int i=0;i<n;i++) {
    if (i==0) c[i] = b0/dpdx;
    else if (i==(n-1)) c[i] = bn/dpdx;
    else c[i] = 3*(y[i+1] - y[i-1]);
    memset(A[i],0,n*sizeof(float));
    for (int j=i-1;j<=i+1;j++) if (j>=0 && j<n) {
      if (i==j) if (i==0 || i==(n-1)) A[i][j] = 1; else A[i][j] = 4;
      else if (i!=0 && i!=(n-1)) A[i][j] = 1;
    }
  }
  gelimd2(A,c,b,n);
  for (int i=0;i<n;i++) {
    if (i<(n-1)) {
      c[i] = 3*(a[i+1] - a[i]) - 2*b[i] - b[i+1];
      d[i] = 2*(a[i] - a[i+1]) + b[i] + b[i+1];
    } else {
      c[i] = 3*(-a[i]) - 2*b[i];
      d[i] = 2*(a[i]) + b[i];
    }
  }
}

float USpline::evalf(float x) {
  float p = (x-xmin)/span;
  if (p<0) p = 0;
  if (p>1) p = 1;
  p *= (n-1);
  int i = int(p);
  p -= float(i);
  float pp = p*p;
  return a[i] + b[i]*p + c[i]*pp + d[i]*pp*p;
}

float USpline::evaldf(float x) {
  float p = (x-xmin)/span;
  if (p<0) p = 0;
  if (p>1) p = 1;
  p *= (n-1);
  int i = int(p);
  p -= float(i);
  return dpdx*(b[i] + 2*c[i]*p + 3*d[i]*p*p);
}

void USpline::print() {
  cout << n << " " << xmin << " " << xmax << " " << span << " " << interval
       << endl << endl;
  for (int i=0;i<n;i++) {
    for (int j=0;j<n;j++) {
      cout << A[i][j] << " ";
    }
    cout << "   " << a[i] << " " << b[i] << " " << c[i] << " " << d[i] << endl;
  }
}

} // namespace bmtk

---

bmtk
Author(s): Benjamin Pitzer
autogenerated on Mon Mar 4 11:05:51 2013