diff --git a/contrib/FourierModel/ContinuationPatch.cpp b/contrib/FourierModel/ContinuationPatch.cpp
index 49b62c9cb549423e28eadad32c58c5f4745735ee..62b6d01f5bcafb47710777ffe14e41bef3ba06b8 100644
--- a/contrib/FourierModel/ContinuationPatch.cpp
+++ b/contrib/FourierModel/ContinuationPatch.cpp
@@ -452,7 +452,7 @@ void ContinuationPatch::_ReprocessSeriesCoeff()
delete [] dataU;
delete [] dataV;
}
- else {
+ else if (IsVPeriodic()){
std::vector<double> u(_uM + 1), v(_vM);
for (int j = 0; j < _uM + 1; j++)
u[j] = (double)j / (double)_uM;
@@ -496,6 +496,52 @@ void ContinuationPatch::_ReprocessSeriesCoeff()
delete [] dataV;
delete [] dataU;
}
+ else {
+ std::vector<double> u(_uM + 1), v(_vM + 1);
+ for (int j = 0; j < _uM + 1; j++)
+ u[j] = (double)j / (double)_uM;
+ for (int j = 0; j < _vM + 1; j++)
+ v[j] = (double)j / (double)_vM;
+
+ std::complex<double> **dataU = new std::complex<double> *[_vM];
+ for (int k = 0; k < _vM; k++)
+ dataU[k] = new std::complex<double> [2*_uM + 1];
+
+ std::complex<double> *dataV = new std::complex<double>[2*_vM + 1];
+ for (int j = 0; j < _uM + 1; j++) {
+ for (int k = 0; k < _vM + 1; k++) {
+ dataV[k] = _Interpolate(0.5 * cos(M_PI * u[j]) + 0.5,
+ 0.5 * cos(M_PI * v[k]) + 0.5);
+ }
+ for (int k = 1; k < _vM+1; k++)
+ dataV[_vM + k] = dataV[_vM - k];
+ _BackwardFft(2*_vM + 1, dataV);
+ dataU[0][j] = 0.5 * dataV[0] / (double)_vM;
+ for (int k=1; k<_vM-1; k++)
+ dataU[k][j] = dataV[k] / (double)_vM;
+ dataU[_vM-1][j] = 0.5 * dataV[_vM-1] / (double)_vM;
+ }
+
+ for (int k = 0; k < _vM; k++) {
+ for (int j = 1; j < _uM+1; j++) {
+ dataU[k][_uM + j] = dataU[k][_uM - k];
+ }
+ _BackwardFft(2*_uM + 1, dataU[k]);
+ }
+
+ for (int j = 0; j < _uM; j++) {
+ for (int k = 0; k < _vM; k++) {
+ if ((j == 0) || (j == _uM - 1))
+ _coeffData[j][k] = 0.5 * dataU[k][j]/ (double)_uM;
+ else
+ _coeffData[j][k] = dataU[k][j]/ (double)_uM;
+ }
+ }
+ for (int k = 0; k < _vM; k++)
+ delete [] dataU[k];
+ delete [] dataU;
+ delete [] dataV;
+ }
}
std::complex<double> ContinuationPatch::
diff --git a/contrib/FourierModel/CylindricalProjectionSurface.cpp b/contrib/FourierModel/CylindricalProjectionSurface.cpp
index d3efbd70c59c2786166b41768f240ec42f224cdc..9ea0ed6a8637db7c4fc6c33a26d0f7598913a5f9 100755
--- a/contrib/FourierModel/CylindricalProjectionSurface.cpp
+++ b/contrib/FourierModel/CylindricalProjectionSurface.cpp
@@ -4,7 +4,7 @@ CylindricalProjectionSurface::CylindricalProjectionSurface
(int tag) : ProjectionSurface(1.)
{
SetTag(tag);
- SetName(std::string("Cylinder"));
+ SetName(std::string("cylinder"));
twoPi_ = 2 * M_PI;
diff --git a/contrib/FourierModel/FM_Reader.cpp b/contrib/FourierModel/FM_Reader.cpp
index 5c745756e9ba0008f1a3063a96a80d63ea64ca78..11666b29cb6e151764786c4da32f19c89acbc4b9 100644
--- a/contrib/FourierModel/FM_Reader.cpp
+++ b/contrib/FourierModel/FM_Reader.cpp
@@ -44,12 +44,9 @@ FM_Reader::FM_Reader(const char* fn)
else if (!strcmp(psName,cylinder))
_ps[i] = new CylindricalProjectionSurface
(psTag,origin,E0,E1,scale,psParams[0],psParams[1]);
- else if (!strcmp(psName,paraboloid)) {
- double K[2];
- K[0] = psParams[0]; K[1] = psParams[1];
+ else if (!strcmp(psName,paraboloid))
_ps[i] = new ParaboloidProjectionSurface
- (psTag,origin,E0,E1,scale,K);
- }
+ (psTag,origin,E0,E1,scale);
else if (!strcmp(psName,revolvedParabola)) {
double R = psParams[0];
double K[2];
diff --git a/contrib/FourierModel/FPatch.cpp b/contrib/FourierModel/FPatch.cpp
index 8c897275d1bf16e69bc54e0833fe50eab7effcbc..914713ff08b37de2def16ace375612b3424ab833 100644
--- a/contrib/FourierModel/FPatch.cpp
+++ b/contrib/FourierModel/FPatch.cpp
@@ -497,7 +497,7 @@ void FPatch::_ReprocessSeriesCoeff()
delete [] dataU;
delete [] dataV;
}
- else {
+ else if (_ps->IsVPeriodic()) {
std::vector<double> u(_uM + 1), v(_vM);
for (int j = 0; j < _uM + 1; j++)
u[j] = (double)j / (double)_uM;
@@ -541,6 +541,52 @@ void FPatch::_ReprocessSeriesCoeff()
delete [] dataV;
delete [] dataU;
}
+ else {
+ std::vector<double> u(_uM + 1), v(_vM + 1);
+ for (int j = 0; j < _uM + 1; j++)
+ u[j] = (double)j / (double)_uM;
+ for (int j = 0; j < _vM + 1; j++)
+ v[j] = (double)j / (double)_vM;
+
+ std::complex<double> **dataU = new std::complex<double> *[_vM];
+ for (int k = 0; k < _vM; k++)
+ dataU[k] = new std::complex<double> [2*_uM + 1];
+
+ std::complex<double> *dataV = new std::complex<double>[2*_vM + 1];
+ for (int j = 0; j < _uM + 1; j++) {
+ for (int k = 0; k < _vM + 1; k++) {
+ dataV[k] = _Interpolate(0.5 * cos(M_PI * u[j]) + 0.5,
+ 0.5 * cos(M_PI * v[k]) + 0.5);
+ }
+ for (int k = 1; k < _vM+1; k++)
+ dataV[_vM + k] = dataV[_vM - k];
+ _BackwardFft(2*_vM + 1, dataV);
+ dataU[0][j] = 0.5 * dataV[0] / (double)_vM;
+ for (int k=1; k<_vM-1; k++)
+ dataU[k][j] = dataV[k] / (double)_vM;
+ dataU[_vM-1][j] = 0.5 * dataV[_vM-1] / (double)_vM;
+ }
+
+ for (int k = 0; k < _vM; k++) {
+ for (int j = 1; j < _uM+1; j++) {
+ dataU[k][_uM + j] = dataU[k][_uM - k];
+ }
+ _BackwardFft(2*_uM + 1, dataU[k]);
+ }
+
+ for (int j = 0; j < _uM; j++) {
+ for (int k = 0; k < _vM; k++) {
+ if ((j == 0) || (j == _uM - 1))
+ _coeffData[j][k] = 0.5 * dataU[k][j]/ (double)_uM;
+ else
+ _coeffData[j][k] = dataU[k][j]/ (double)_uM;
+ }
+ }
+ for (int k = 0; k < _vM; k++)
+ delete [] dataU[k];
+ delete [] dataU;
+ delete [] dataV;
+ }
}
std::complex<double> FPatch::
diff --git a/contrib/FourierModel/ParaboloidProjectionSurface.cpp b/contrib/FourierModel/ParaboloidProjectionSurface.cpp
index 5407cca67755f8275205a6b73aafcaeaf486e5ce..26cb5647ee8d9288a768ace42a40557546136357 100755
--- a/contrib/FourierModel/ParaboloidProjectionSurface.cpp
+++ b/contrib/FourierModel/ParaboloidProjectionSurface.cpp
@@ -1,15 +1,12 @@
#include "ParaboloidProjectionSurface.h"
ParaboloidProjectionSurface::ParaboloidProjectionSurface
-(int tag) : ProjectionSurface(1.)
+(int tag) : ProjectionSurface()
{
SetTag(tag);
SetName(std::string("paraboloid"));
- twoPi_ = 2 * M_PI;
- K_[0] = K_[1] = 0.5;
-
- numParameters_ = 2;
+ numParameters_ = 0;
O_[0] = O_[1] = O_[2] = 0.;
@@ -21,15 +18,12 @@ ParaboloidProjectionSurface::ParaboloidProjectionSurface
}
ParaboloidProjectionSurface::ParaboloidProjectionSurface
-(int tag, double O[3], double E0[3], double E1[3], double scale[3],
- double K[2]) : ProjectionSurface(1.)
+(int tag, double O[3], double E0[3], double E1[3], double scale[3])
+ : ProjectionSurface()
{
SetTag(tag);
SetName(std::string("Paraboloid"));
- twoPi_ = 2 * M_PI;
- K_[0] = K[0]; K_[1] = K[1];
-
O_[0] = O[0]; O_[1] = O[1]; O_[2] = O[2];
E0_[0] = E0[0]; E0_[1] = E0[1]; E0_[2] = E0[2];
@@ -45,15 +39,15 @@ ParaboloidProjectionSurface::ParaboloidProjectionSurface
void ParaboloidProjectionSurface::F
(double u, double v, double &x, double &y, double &z)
{
- x = O_[0] + K_[0] * v * (scale_[1] * E1_[0] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[0] * sin(twoPi_ * (u - 0.5))) +
- K_[1] * v * v * scale_[0] * E0_[0];
- y = O_[1] + K_[0] * v * (scale_[1] * E1_[1] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[1] * sin(twoPi_ * (u - 0.5))) +
- K_[1] * v * v * scale_[0] * E0_[1];
- z = O_[2] + K_[0] * v * (scale_[1] * E1_[2] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[2] * sin(twoPi_ * (u - 0.5))) +
- K_[1] * v * v * scale_[0] * E0_[2];
+ x = O_[0] + (u - 0.5) * scale_[1] * E1_[0] + (v - 0.5) * scale_[2] * E2_[0] +
+ ((u - 0.5) * scale_[1] * (u - 0.5) * scale_[1] +
+ (v - 0.5) * scale_[2] * (v - 0.5) * scale_[2]) * scale_[0] * E0_[0];
+ y = O_[1] + (u - 0.5) * scale_[1] * E1_[1] + (v - 0.5) * scale_[2] * E2_[1] +
+ ((u - 0.5) * scale_[1] * (u - 0.5) * scale_[1] +
+ (v - 0.5) * scale_[2] * (v - 0.5) * scale_[2]) * scale_[0] * E0_[1];
+ z = O_[0] + (u - 0.5) * scale_[1] * E1_[2] + (v - 0.5) * scale_[2] * E2_[2] +
+ ((u - 0.5) * scale_[1] * (u - 0.5) * scale_[1] +
+ (v - 0.5) * scale_[2] * (v - 0.5) * scale_[2]) * scale_[0] * E0_[2];
}
bool ParaboloidProjectionSurface::Inverse
@@ -73,98 +67,61 @@ bool ParaboloidProjectionSurface::Inverse
RdotT /= scale_[1];
RdotNcT /= scale_[2];
- u = atan2(RdotNcT,RdotT);
- u /= twoPi_;
- u += 0.5;
-
- v = sqrt(RdotT * RdotT + RdotNcT * RdotNcT);
- v /= K_[0];
+ u = RdotT + 0.5;
+ v = RdotNcT + 0.5;
}
void ParaboloidProjectionSurface::Dfdu
(double u, double v, double &x, double &y, double &z)
{
- x = twoPi_ * K_[0] * v * (- scale_[1] * E1_[0] * sin(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[0] * cos(twoPi_ * (u - 0.5)));
- y = twoPi_ * K_[0] * v * (- scale_[1] * E1_[1] * sin(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[1] * cos(twoPi_ * (u - 0.5)));
- z = twoPi_ * K_[0] * v * (- scale_[1] * E1_[2] * sin(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[2] * cos(twoPi_ * (u - 0.5)));
+ x = scale_[1] * E1_[0] +
+ 2 * scale_[0] * scale_[1] * scale_[1] * (u - 0.5) * E0_[0];
+ y = scale_[1] * E1_[1] +
+ 2 * scale_[0] * scale_[1] * scale_[1] * (u - 0.5) * E0_[1];
+ z = scale_[1] * E1_[2] +
+ 2 * scale_[0] * scale_[1] * scale_[1] * (u - 0.5) * E0_[2];
}
void ParaboloidProjectionSurface::Dfdv
(double u, double v, double &x, double &y, double &z)
{
- x = K_[0] * (scale_[1] * E1_[0] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[0] * sin(twoPi_ * (u - 0.5))) +
- 2. * K_[1] * v * scale_[0] * E0_[0];
- y = K_[0] * (scale_[1] * E1_[1] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[1] * sin(twoPi_ * (u - 0.5))) +
- 2. * K_[1] * v * scale_[0] * E0_[1];
- z = K_[0] * (scale_[1] * E1_[2] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[2] * sin(twoPi_ * (u - 0.5))) +
- 2. * K_[1] * v * scale_[0] * E0_[2];
+ x = scale_[2] * E2_[0] + 2 * scale_[0] * scale_[2] * (v - 0.5) * E0_[0];
+ y = scale_[2] * E2_[1] + 2 * scale_[0] * scale_[2] * (v - 0.5) * E0_[1];
+ z = scale_[2] * E2_[2] + 2 * scale_[0] * scale_[2] * (v - 0.5) * E0_[2];
}
void ParaboloidProjectionSurface::Dfdfdudu
(double u,double v, double &x, double &y, double &z)
{
- x = - twoPi_ * twoPi_ * K_[0] * v *
- (scale_[1] * E1_[0] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[0] * sin(twoPi_ * (u - 0.5)));
- y = - twoPi_ * twoPi_ * K_[0] * v *
- (scale_[1] * E1_[1] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[1] * sin(twoPi_ * (u - 0.5)));
- z = - twoPi_ * twoPi_ * K_[0] * v *
- (scale_[1] * E1_[2] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[2] * sin(twoPi_ * (u - 0.5)));
+ x = 2 * scale_[0] * scale_[1] * E0_[0];
+ y = 2 * scale_[0] * scale_[1] * E0_[1];
+ z = 2 * scale_[0] * scale_[1] * E0_[2];
}
void ParaboloidProjectionSurface::Dfdfdudv
(double u, double v, double &x, double &y, double &z)
{
- x = K_[0] * twoPi_ * (- scale_[1] * E1_[0] * sin(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[0] * cos(twoPi_ * (u - 0.5)));
- y = K_[0] * twoPi_ * (- scale_[1] * E1_[1] * sin(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[1] * cos(twoPi_ * (u - 0.5)));
- z = K_[0] * twoPi_ * (- scale_[1] * E1_[2] * sin(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[2] * cos(twoPi_ * (u - 0.5)));
+ x = y = z = 0.;
}
void ParaboloidProjectionSurface::Dfdfdvdv
(double u, double v, double &x, double &y, double &z)
{
- x = 2. * K_[1] * scale_[0] * E0_[0];
- y = 2. * K_[1] * scale_[0] * E0_[1];
- z = 2. * K_[1] * scale_[0] * E0_[2];
+ x = 2. * scale_[0] * scale_[2] * E0_[0];
+ y = 2. * scale_[0] * scale_[2] * E0_[1];
+ z = 2. * scale_[0] * scale_[2] * E0_[2];
}
void ParaboloidProjectionSurface::Dfdfdfdududu
(double u,double v,double &x,double &y,double &z)
{
- x = twoPi_ * twoPi_ * twoPi_ * K_[0] * v *
- (scale_[1] * E1_[0] * sin(twoPi_ * (u - 0.5)) -
- scale_[2] * E2_[0] * cos(twoPi_ * (u - 0.5)));
- y = twoPi_ * twoPi_ * twoPi_ * K_[0] * v *
- (scale_[1] * E1_[1] * sin(twoPi_ * (u - 0.5)) -
- scale_[2] * E2_[1] * cos(twoPi_ * (u - 0.5)));
- z = twoPi_ * twoPi_ * twoPi_ * K_[0] * v *
- (scale_[1] * E1_[2] * sin(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[2] * cos(twoPi_ * (u - 0.5)));
+ x = y = z = 0.;
}
void ParaboloidProjectionSurface::Dfdfdfdududv
(double u,double v,double &x,double &y,double &z)
{
- x = - twoPi_ * twoPi_ * K_[0] *
- (scale_[1] * E1_[0] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[0] * sin(twoPi_ * (u - 0.5)));
- y = - twoPi_ * twoPi_ * K_[0] *
- (scale_[1] * E1_[1] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[1] * sin(twoPi_ * (u - 0.5)));
- z = - twoPi_ * twoPi_ * K_[0] *
- (scale_[1] * E1_[2] * cos(twoPi_ * (u - 0.5)) +
- scale_[2] * E2_[2] * sin(twoPi_ * (u - 0.5)));
+ x = y = z = 0.;
}
void ParaboloidProjectionSurface::Dfdfdfdudvdv
@@ -182,6 +139,7 @@ void ParaboloidProjectionSurface::Dfdfdfdvdvdv
bool ParaboloidProjectionSurface::OrthoProjectionOnSurface
(double x, double y, double z, double& u, double& v)
{
+ /*
double R[3];
R[0] = x - O_[0];
R[1] = y - O_[1];
@@ -219,18 +177,7 @@ bool ParaboloidProjectionSurface::OrthoProjectionOnSurface
std::vector<double> root = SolveCubic(a,b,c);
if (root.size()) {
- /*
- double xP,yP,zP;
- double minDist = 1.e12;
- for (int i=0;i<root.size();i++) {
- F(u,root[i],xP,yP,zP);
- double dist = sqrt((x-xP)*(x-xP)+(y-yP)*(y-yP)+(z-zP)*(z-zP));
- if (dist < minDist) {
- minDist = dist;
- v = root[i];
- }
- }
- */
+
for (int i=0;i<root.size();i++) {
if (root[i] >= 0.) {
v = root[i];
@@ -245,37 +192,22 @@ bool ParaboloidProjectionSurface::OrthoProjectionOnSurface
}
else
return false;
+*/
}
void ParaboloidProjectionSurface::
SetParameter(int i, double x)
{
- switch (i) {
- case 0:
- K_[0] = x;
- case 1:
- K_[1] = x;
- }
}
double ParaboloidProjectionSurface::
GetParameter(int i)
{
- switch (i) {
- case 0:
- return K_[0];
- case 1:
- return K_[1];
- }
+ return 0.;
}
std::string ParaboloidProjectionSurface::
GetLabel(int i)
{
- switch (i) {
- case 0:
- return std::string("Parabola Scale 1");
- case 1:
- return std::string("Parabola Scale 2");
- }
+ return std::string(" ");
}
diff --git a/contrib/FourierModel/ParaboloidProjectionSurface.h b/contrib/FourierModel/ParaboloidProjectionSurface.h
index 0d0a5f58ae93f8c7853b94e715f2ec3379b96460..dcefbe038e3a3e24db416eb51af8e3e969b94402 100755
--- a/contrib/FourierModel/ParaboloidProjectionSurface.h
+++ b/contrib/FourierModel/ParaboloidProjectionSurface.h
@@ -7,21 +7,14 @@
#include "ProjectionSurface.h"
class ParaboloidProjectionSurface : public ProjectionSurface {
- private:
- double twoPi_;
- double K_[2];
public:
ParaboloidProjectionSurface
(int tag);
ParaboloidProjectionSurface
- (int tag, double O[3], double E0[3], double E1[3], double scale[3],
- double K[2]);
+ (int tag, double O[3], double E0[3], double E1[3], double scale[3]);
ParaboloidProjectionSurface(ParaboloidProjectionSurface *ps)
: ProjectionSurface(ps)
{
- twoPi_ = ps->twoPi_;
- K_[0] = ps->K_[0];
- K_[1] = ps->K_[1];
}
virtual ~ParaboloidProjectionSurface
diff --git a/contrib/FourierModel/Utils.cpp b/contrib/FourierModel/Utils.cpp
index 344a4e209312084abbe8802c418245e5b676f87a..9d951238bf8ebb8ba10351b380e38e4179cbd559 100755
--- a/contrib/FourierModel/Utils.cpp
+++ b/contrib/FourierModel/Utils.cpp
@@ -84,6 +84,91 @@ std::vector<double> SolveCubic(double a, double b, double c)
return root;
}
+std::vector<double> SolveCubic(double a, double b, double c, double d)
+{
+ std::vector<double> root;
+ // find real roots of polynomial a*x^3 + b*x^2 + c*x + d=0
+
+ double tol=pow(10.0, -12.);
+
+ double error=1.;
+ int maxiter=100; int iter=0;
+
+ std::vector<double> initial_guess(5);
+ initial_guess[0]=0.;
+ initial_guess[1]=1.;
+ initial_guess[2]=-1.;
+ initial_guess[3]=0.5;
+ initial_guess[4]=-0.5;
+
+ bool find=0; int init=0;
+ double x0=0., x1, f, df;
+
+ while (find==0 && init< 5){
+
+ init=init+1;
+ x0=initial_guess[init-1];
+ error=1.; iter=0;
+
+ while (error>tol && iter<maxiter){
+
+ f=a*x0*x0*x0+b*x0*x0+c*x0+d;
+ df=3.*a*x0*x0+2.*b*x0+c;
+
+ x1=x0-f/df;
+
+ error=std::abs(x1-x0);
+ x0=x1;
+ iter=iter+1;
+ }
+
+ if (iter<maxiter){
+ find=1;
+ }
+ }
+
+ if (find) {
+ root.push_back(x0);
+
+ //%polynomial a*x^3 + b*x^2 + c*x + d = (x-x0)*(a*x^2+bb*x+cc)
+ //%where bb=a*x0, cc=-c/x0 (=b-a*x^2);
+
+ double aa=a, bb=b+a*x0, cc=c+b*x0+a*x0*x0;
+ double D=bb*bb-4.*aa*cc;
+
+ if (D>=0.){
+ root.push_back((-bb+sqrt(D))/2./aa);
+ root.push_back((-bb-sqrt(D))/2./aa);
+ }
+
+ /*
+ // it seems that the solution obtained by quadratic formula
+ //is not so accurate --> used as initial guess.
+
+ if (root.size() == 3){
+ for (int i=0 ; i < 2 ; i++){
+ x0=root[i+1];
+ error=1.; iter=0;
+ while (error>tol && iter<maxiter){
+
+ f=a*x0*x0*x0+b*x0+c;
+ df=a*3.*x0*x0+b;
+
+ x1=x0-f/df;
+
+ error=std::abs(x1-x0);
+ x0=x1;
+ iter=iter+1;
+ }
+ root[i+1]=x0;
+ }
+ }
+ */
+ }
+
+ return root;
+}
+
void find(std::vector<int> &a, int length, std::vector<int> &q, int &num){
int i,j=0;
diff --git a/contrib/FourierModel/Utils.h b/contrib/FourierModel/Utils.h
index 7fbfabbcd50d327fc91e045e691d5a3756ea6cc0..51c1191341ec81ade7848508b1da09fa5442e144 100755
--- a/contrib/FourierModel/Utils.h
+++ b/contrib/FourierModel/Utils.h
@@ -7,6 +7,7 @@
#include <iostream>
std::vector<double> SolveCubic(double a, double b, double c);
+std::vector<double> SolveCubic(double a, double b, double c, double d);
void find(std::vector<int> &a, int length, std::vector<int> &q, int &num);
int minVec(std::vector<int> &a,int n);
int maxVec(std::vector<int> &a,int n);