1359 lines
33 KiB
C
1359 lines
33 KiB
C
|
|
#ifdef newc
|
|
#include <iostream>
|
|
#include <iomanip>
|
|
#include <fstream>
|
|
#include <sstream>
|
|
#include <strstream>
|
|
#include <cmath>
|
|
using namespace std;
|
|
#else
|
|
#include <iostream.h>
|
|
#include <iomanip.h>
|
|
#include <fstream.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#endif
|
|
#include <mpi.h>
|
|
|
|
#include "misc.h"
|
|
#include "macrodef.h"
|
|
#include "zbesh.h"
|
|
|
|
#define PI M_PI
|
|
|
|
void misc::tillherecheck(int myrank)
|
|
{
|
|
int atp = 1, tatp;
|
|
MPI_Allreduce(&atp, &tatp, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
if (myrank == 0)
|
|
cout << " here now: " << tatp << " processors." << endl;
|
|
}
|
|
void misc::tillherecheck(const char str[])
|
|
{
|
|
int myrank;
|
|
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
|
|
int atp = 1, tatp;
|
|
MPI_Allreduce(&atp, &tatp, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
if (myrank == 0)
|
|
{
|
|
cout << " here now: " << tatp << " processors." << endl;
|
|
cout << str << endl;
|
|
}
|
|
}
|
|
void misc::tillherecheck(MPI_Comm Comm_here, int out_rank, const char str[])
|
|
{
|
|
int myrank;
|
|
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
|
|
int atp = 1, tatp;
|
|
|
|
MPI_Allreduce(&atp, &tatp, 1, MPI_INT, MPI_SUM, Comm_here);
|
|
if (myrank == out_rank)
|
|
{
|
|
cout << " here now: " << tatp << " processors." << endl;
|
|
cout << str << endl;
|
|
}
|
|
}
|
|
void misc::tillherecheck(MPI_Comm Comm_here, int out_rank, const string str)
|
|
{
|
|
int myrank;
|
|
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
|
|
int atp = 1, tatp;
|
|
|
|
MPI_Allreduce(&atp, &tatp, 1, MPI_INT, MPI_SUM, Comm_here);
|
|
if (myrank == out_rank)
|
|
{
|
|
cout << " here now: " << tatp << " processors." << endl;
|
|
cout << str << endl;
|
|
}
|
|
}
|
|
// pick out value from input string
|
|
int misc::parse_parts(string str, string &sgrp, string &skey, string &sval, int &ind)
|
|
{
|
|
int pos1, pos2;
|
|
string s0;
|
|
|
|
ind = 0;
|
|
|
|
// remove comments
|
|
str = str.substr(0, str.find("#"));
|
|
if (rTrim(str).empty())
|
|
return 0; // continue;
|
|
|
|
// parse {group, key, val}
|
|
pos1 = str.find("::");
|
|
pos2 = str.find("=");
|
|
if (pos1 == string::npos || pos2 == string::npos)
|
|
return -1;
|
|
|
|
s0 = str.substr(0, pos1);
|
|
sgrp = lTrim(s0);
|
|
s0 = str.substr(pos1 + 2, pos2 - pos1 - 2);
|
|
skey = rTrim(s0);
|
|
s0 = str.substr(pos2 + 1);
|
|
sval = Trim(s0);
|
|
|
|
pos1 = sval.find("\"");
|
|
pos2 = sval.rfind("\"");
|
|
if (pos1 != string::npos)
|
|
{
|
|
sval = sval.substr(1, pos2 - 1);
|
|
}
|
|
|
|
pos1 = skey.find("[");
|
|
pos2 = skey.find("]");
|
|
if (pos1 != string::npos)
|
|
{
|
|
s0 = skey.substr(0, pos1);
|
|
ind = atoi(skey.substr(pos1 + 1, pos2 - pos1 - 1).c_str());
|
|
skey = s0;
|
|
}
|
|
|
|
return 1;
|
|
}
|
|
int misc::parse_parts(string str, string &sgrp, string &skey, string &sval, int &ind1, int &ind2)
|
|
{
|
|
int pos1, pos2;
|
|
string s0, s1;
|
|
|
|
ind1 = ind2 = 0;
|
|
|
|
// remove comments
|
|
str = str.substr(0, str.find("#"));
|
|
if (rTrim(str).empty())
|
|
return 0; // continue;
|
|
|
|
// parse {group, key, val}
|
|
pos1 = str.find("::");
|
|
pos2 = str.find("=");
|
|
if (pos1 == string::npos || pos2 == string::npos)
|
|
return -1;
|
|
|
|
s0 = str.substr(0, pos1);
|
|
sgrp = lTrim(s0);
|
|
s0 = str.substr(pos1 + 2, pos2 - pos1 - 2);
|
|
skey = rTrim(s0);
|
|
s0 = str.substr(pos2 + 1);
|
|
sval = Trim(s0);
|
|
|
|
pos1 = sval.find("\"");
|
|
pos2 = sval.rfind("\"");
|
|
if (pos1 != string::npos)
|
|
{
|
|
sval = sval.substr(1, pos2 - 1);
|
|
}
|
|
|
|
pos1 = skey.find("[");
|
|
pos2 = skey.find("]");
|
|
if (pos1 != string::npos)
|
|
{
|
|
s0 = skey.substr(0, pos1);
|
|
s1 = skey.substr(pos2 + 1);
|
|
ind1 = atoi(skey.substr(pos1 + 1, pos2 - pos1 - 1).c_str());
|
|
skey = s0;
|
|
}
|
|
|
|
pos1 = s1.find("[");
|
|
pos2 = s1.find("]");
|
|
if (pos1 != string::npos)
|
|
{
|
|
s0 = s1.substr(pos2 + 1);
|
|
ind2 = atoi(s1.substr(pos1 + 1, pos2 - pos1 - 1).c_str());
|
|
}
|
|
|
|
return 1;
|
|
}
|
|
int misc::parse_parts(string str, string &sgrp, string &skey, string &sval, int &ind1, int &ind2, int &ind3)
|
|
{
|
|
int pos1, pos2;
|
|
string s0, s1;
|
|
|
|
ind1 = ind2 = ind3 = 0;
|
|
|
|
// remove comments
|
|
str = str.substr(0, str.find("#"));
|
|
if (rTrim(str).empty())
|
|
return 0; // continue;
|
|
|
|
// parse {group, key, val}
|
|
pos1 = str.find("::");
|
|
pos2 = str.find("=");
|
|
if (pos1 == string::npos || pos2 == string::npos)
|
|
return -1;
|
|
|
|
s0 = str.substr(0, pos1);
|
|
sgrp = lTrim(s0);
|
|
s0 = str.substr(pos1 + 2, pos2 - pos1 - 2);
|
|
skey = rTrim(s0);
|
|
s0 = str.substr(pos2 + 1);
|
|
sval = Trim(s0);
|
|
|
|
pos1 = sval.find("\"");
|
|
pos2 = sval.rfind("\"");
|
|
if (pos1 != string::npos)
|
|
{
|
|
sval = sval.substr(1, pos2 - 1);
|
|
}
|
|
|
|
pos1 = skey.find("[");
|
|
pos2 = skey.find("]");
|
|
if (pos1 != string::npos)
|
|
{
|
|
s0 = skey.substr(0, pos1);
|
|
s1 = skey.substr(pos2 + 1);
|
|
ind1 = atoi(skey.substr(pos1 + 1, pos2 - pos1 - 1).c_str());
|
|
skey = s0;
|
|
}
|
|
|
|
pos1 = s1.find("[");
|
|
pos2 = s1.find("]");
|
|
if (pos1 != string::npos)
|
|
{
|
|
s0 = s1.substr(pos2 + 1);
|
|
ind2 = atoi(s1.substr(pos1 + 1, pos2 - pos1 - 1).c_str());
|
|
}
|
|
|
|
pos1 = s0.find("[");
|
|
pos2 = s0.find("]");
|
|
if (pos1 != string::npos)
|
|
{
|
|
ind3 = atoi(s0.substr(pos1 + 1, pos2 - pos1 - 1).c_str());
|
|
}
|
|
|
|
return 1;
|
|
}
|
|
// sent me from Roman Gold on 2010-10-8
|
|
void misc::gaulegf(double x1, double x2, double *x, double *w, int n)
|
|
{
|
|
int i, j, m;
|
|
double eps = 1.2E-16;
|
|
double p1, p2, p3, pp, xl, xm, z, z1;
|
|
|
|
m = (n + 1) / 2;
|
|
xm = 0.5 * (x2 + x1);
|
|
xl = 0.5 * (x2 - x1);
|
|
for (i = 0; i < m; i++)
|
|
{
|
|
z = cos(PI * ((double)i + 0.75) / ((double)n + 0.5));
|
|
do
|
|
{
|
|
p1 = 1.0;
|
|
p2 = 0.0;
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
p3 = p2;
|
|
p2 = p1;
|
|
p1 = ((2 * (double)j + 1) * z * p2 - (double)j * p3) / ((double)j + 1);
|
|
}
|
|
pp = n * (z * p1 - p2) / (z * z - 1.0);
|
|
z1 = z;
|
|
z = z1 - p1 / pp;
|
|
} while (fabs(z - z1) > eps);
|
|
x[i] = xm - xl * z;
|
|
x[n - 1 - i] = xm + xl * z;
|
|
w[i] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
|
|
w[n - 1 - i] = w[i];
|
|
}
|
|
} /* end gaulegf */
|
|
void misc::inversearray(double *aa, int NN)
|
|
{
|
|
int i, m;
|
|
m = (NN + 1) / 2;
|
|
double rr;
|
|
for (i = 0; i < m; i++)
|
|
{
|
|
rr = aa[i];
|
|
aa[i] = aa[NN - 1 - i];
|
|
aa[NN - 1 - i] = rr;
|
|
}
|
|
}
|
|
// Eq.(42) of PRD 77, 024027 (2008)
|
|
double misc::Wigner_d_function(int l, int m, int s, double costheta)
|
|
{
|
|
// we consider only theta in [0,pi]
|
|
int C1 = max(0, m - s), C2 = min(l + m, l - s);
|
|
|
|
double vv = 0;
|
|
double sinht = sqrt((1 - costheta) / 2.0), cosht = sqrt((1 + costheta) / 2.0);
|
|
if (C1 % 2 == 0)
|
|
{
|
|
for (int t = C1; t < C2 + 1; t += 2)
|
|
vv = vv + pow(cosht, 2 * l + m - s - 2 * t) * pow(sinht, 2 * t + s - m) /
|
|
(fact(l + m - t) * fact(l - s - t) * fact(t) * fact(t + s - m));
|
|
for (int t = C1 + 1; t < C2 + 1; t += 2)
|
|
vv = vv - pow(cosht, 2 * l + m - s - 2 * t) * pow(sinht, 2 * t + s - m) /
|
|
(fact(l + m - t) * fact(l - s - t) * fact(t) * fact(t + s - m));
|
|
}
|
|
else
|
|
{
|
|
for (int t = C1; t < C2 + 1; t += 2)
|
|
vv = vv - pow(cosht, 2 * l + m - s - 2 * t) * pow(sinht, 2 * t + s - m) /
|
|
(fact(l + m - t) * fact(l - s - t) * fact(t) * fact(t + s - m));
|
|
for (int t = C1 + 1; t < C2 + 1; t += 2)
|
|
vv = vv + pow(cosht, 2 * l + m - s - 2 * t) * pow(sinht, 2 * t + s - m) /
|
|
(fact(l + m - t) * fact(l - s - t) * fact(t) * fact(t + s - m));
|
|
}
|
|
return vv * sqrt(fact(l + m) * fact(l - m) * fact(l + s) * fact(l - s));
|
|
}
|
|
double misc::fact(int N)
|
|
{
|
|
if (N < 0)
|
|
cout << "error input for factorial." << endl;
|
|
double f;
|
|
if (N == 0)
|
|
f = 1;
|
|
else
|
|
f = N * fact(N - 1);
|
|
return f;
|
|
}
|
|
int misc::num_of_str(char *c)
|
|
{
|
|
int NN = 0, N1 = 0;
|
|
std::istringstream iss;
|
|
iss.str(c);
|
|
|
|
char c1[1000];
|
|
while (!iss.eof())
|
|
{
|
|
iss >> c1;
|
|
if (int(c1[0]) == 45 || int(c1[0]) == 46 || (int(c1[0]) > 47 && int(c1[0]) < 58))
|
|
NN++;
|
|
N1++;
|
|
}
|
|
|
|
char *c2 = c;
|
|
while (*(c2 + 1))
|
|
c2++;
|
|
if (int(*c2) == 32)
|
|
{
|
|
NN--;
|
|
N1--;
|
|
}
|
|
|
|
// cout<<"found "<<N1<<" generalized data including "<<NN<<" number type data"<<endl;
|
|
return NN;
|
|
}
|
|
// MNRAS 411, 2461 (2010)
|
|
// Eq.(20)-(22)
|
|
void misc::TVDrungekutta3(const int N, const double dT, double *f0,
|
|
double *f1, double *f_rhs, const int RK4)
|
|
{
|
|
const double F3o4 = 0.75, F1o4 = 0.25, F1o3 = 1.0 / 3, F2o3 = 2.0 / 3;
|
|
switch (RK4)
|
|
{
|
|
case 0:
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
f1[i] = f0[i] + dT * f_rhs[i];
|
|
f_rhs[i] = F1o4 * f1[i];
|
|
}
|
|
break;
|
|
case 1:
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
f1[i] = F3o4 * f0[i] + f_rhs[i] + F1o4 * dT * f1[i];
|
|
f_rhs[i] = F2o3 * f1[i];
|
|
}
|
|
break;
|
|
case 2:
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
f1[i] = F1o3 * f0[i] + f_rhs[i] + F2o3 * dT * f1[i];
|
|
}
|
|
break;
|
|
case 3:
|
|
break;
|
|
default:
|
|
cout << "misc::rungekutta4: something is wrong in RK4 counting!!" << endl;
|
|
}
|
|
}
|
|
void misc::rungekutta4(const int N, const double dT, double *f0,
|
|
double *f1, double *f_rhs, const int RK4)
|
|
{
|
|
const double F1o6 = 1.0 / 6, HLF = 0.5, TWO = 2;
|
|
switch (RK4)
|
|
{
|
|
case 0:
|
|
for (int i = 0; i < N; i++)
|
|
f1[i] = f0[i] + HLF * dT * f_rhs[i];
|
|
break;
|
|
case 1:
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
f_rhs[i] = f_rhs[i] + TWO * f1[i];
|
|
f1[i] = f0[i] + HLF * dT * f1[i];
|
|
}
|
|
break;
|
|
case 2:
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
f_rhs[i] = f_rhs[i] + TWO * f1[i];
|
|
f1[i] = f0[i] + dT * f1[i];
|
|
}
|
|
break;
|
|
case 3:
|
|
for (int i = 0; i < N; i++)
|
|
f1[i] = f0[i] + F1o6 * dT * (f1[i] + f_rhs[i]);
|
|
break;
|
|
default:
|
|
cout << "misc::rungekutta4: something is wrong in RK4 counting!!" << endl;
|
|
}
|
|
}
|
|
void misc::rungekutta4(const double dT, const std::vector<double> &f0,
|
|
std::vector<double> &f1, std::vector<double> &f_rhs, const int RK4)
|
|
{
|
|
const int N = f0.size();
|
|
const double F1o6 = 1.0 / 6, HLF = 0.5, TWO = 2;
|
|
switch (RK4)
|
|
{
|
|
case 0:
|
|
for (int i = 0; i < N; i++)
|
|
f1[i] = f0[i] + HLF * dT * f_rhs[i];
|
|
break;
|
|
case 1:
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
f_rhs[i] = f_rhs[i] + TWO * f1[i];
|
|
f1[i] = f0[i] + HLF * dT * f1[i];
|
|
}
|
|
break;
|
|
case 2:
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
f_rhs[i] = f_rhs[i] + TWO * f1[i];
|
|
f1[i] = f0[i] + dT * f1[i];
|
|
}
|
|
break;
|
|
case 3:
|
|
for (int i = 0; i < N; i++)
|
|
f1[i] = f0[i] + F1o6 * dT * (f1[i] + f_rhs[i]);
|
|
break;
|
|
default:
|
|
cout << "misc::rungekutta4: something is wrong in RK4 counting!!" << endl;
|
|
}
|
|
}
|
|
void misc::dividBlock(const int DIM, int *shape_here, double *bbox_here, const int pices, double *picef, int *shape_res, double *bbox_res,
|
|
const int min_width)
|
|
{
|
|
if (pices < 1)
|
|
{
|
|
cerr << "error in dividBlock: pices = " << pices << endl;
|
|
return;
|
|
}
|
|
if (pices == 1)
|
|
{
|
|
for (int i = 0; i < DIM; i++)
|
|
{
|
|
shape_res[i] = shape_here[i];
|
|
bbox_res[i] = bbox_here[i];
|
|
bbox_res[DIM + i] = bbox_here[DIM + i];
|
|
}
|
|
return;
|
|
}
|
|
|
|
double dd = picef[0];
|
|
for (int i = 1; i < pices; i++)
|
|
dd += picef[i];
|
|
|
|
if (feq(dd, 1, 1e-8))
|
|
{
|
|
int leg = shape_here[0];
|
|
int legi = 0;
|
|
for (int i = 1; i < DIM; i++)
|
|
{
|
|
if (leg < shape_here[i])
|
|
{
|
|
leg = shape_here[i];
|
|
legi = i;
|
|
}
|
|
}
|
|
|
|
int pic = 0;
|
|
|
|
for (int ip = 0; ip < pices; ip++)
|
|
{
|
|
for (int i = 0; i < DIM; i++)
|
|
{
|
|
if (i == legi)
|
|
{
|
|
if (ip == pices - 1)
|
|
shape_res[ip * DIM + i] = shape_here[i] - pic;
|
|
else
|
|
{
|
|
shape_res[ip * DIM + i] = shape_here[i] * picef[ip];
|
|
pic += shape_res[ip * DIM + i];
|
|
}
|
|
}
|
|
else
|
|
shape_res[ip * DIM + i] = shape_here[i];
|
|
}
|
|
}
|
|
|
|
for (int ip = 0; ip < pices; ip++)
|
|
{
|
|
for (int i = 0; i < DIM; i++)
|
|
{
|
|
#ifdef Vertex
|
|
#ifdef Cell
|
|
#error Both Cell and Vertex are defined
|
|
#endif
|
|
dd = (bbox_here[DIM + i] - bbox_here[i]) / (shape_here[i] - 1);
|
|
#else
|
|
#ifdef Cell
|
|
dd = (bbox_here[DIM + i] - bbox_here[i]) / shape_here[i];
|
|
#else
|
|
#error Not define Vertex nor Cell
|
|
#endif
|
|
#endif
|
|
|
|
if (i == legi)
|
|
{
|
|
if (shape_res[ip * DIM + i] < min_width)
|
|
{
|
|
cerr << "dividBlock: resulted too small shape, shapeo = " << shape_here[i] << ", shape = " << shape_res[ip * DIM + i] << ", min_width = " << min_width << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
|
|
if (ip == 0)
|
|
bbox_res[ip * 2 * DIM + i] = bbox_here[i];
|
|
#ifdef Vertex
|
|
#ifdef Cell
|
|
#error Both Cell and Vertex are defined
|
|
#endif
|
|
else
|
|
bbox_res[ip * 2 * DIM + i] = bbox_res[(ip - 1) * 2 * DIM + DIM + i] - ghost_width * dd + dd; // because for ip-1 we have already considered ghost points
|
|
#else
|
|
#ifdef Cell
|
|
else
|
|
bbox_res[ip * 2 * DIM + i] = bbox_res[(ip - 1) * 2 * DIM + DIM + i] - ghost_width * dd;
|
|
#else
|
|
#error Not define Vertex nor Cell
|
|
#endif
|
|
#endif
|
|
|
|
if (ip == pices - 1)
|
|
bbox_res[ip * 2 * DIM + DIM + i] = bbox_here[DIM + i];
|
|
#ifdef Vertex
|
|
#ifdef Cell
|
|
#error Both Cell and Vertex are defined
|
|
#endif
|
|
else
|
|
bbox_res[ip * 2 * DIM + DIM + i] = bbox_res[ip * 2 * DIM + i] + (shape_res[ip * DIM + i] - 1) * dd;
|
|
#else
|
|
#ifdef Cell
|
|
else
|
|
bbox_res[ip * 2 * DIM + DIM + i] = bbox_res[ip * 2 * DIM + i] + shape_res[ip * DIM + i] * dd;
|
|
#else
|
|
#error Not define Vertex nor Cell
|
|
#endif
|
|
#endif
|
|
|
|
if (ip > 0)
|
|
{
|
|
shape_res[ip * DIM + i] += ghost_width;
|
|
bbox_res[ip * 2 * DIM + i] -= ghost_width * dd;
|
|
}
|
|
if (ip < pices - 1)
|
|
{
|
|
shape_res[ip * DIM + i] += ghost_width;
|
|
bbox_res[ip * 2 * DIM + DIM + i] += ghost_width * dd;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
bbox_res[ip * 2 * DIM + i] = bbox_here[i];
|
|
bbox_res[ip * 2 * DIM + DIM + i] = bbox_here[DIM + i];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
cerr << "error in dividBlock: ";
|
|
for (int i = 0; i < pices; i++)
|
|
cerr << picef[i] << " ";
|
|
cerr << endl;
|
|
}
|
|
#if 0
|
|
// for check
|
|
int myrank;
|
|
MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
|
|
if(myrank == 0)
|
|
{
|
|
cerr<<"original one"<<endl;
|
|
cerr<<"shape: (";
|
|
for(int i=0;i<DIM;i++)
|
|
{
|
|
cerr<<shape_here[i];
|
|
if(i<DIM-1) cerr<<" ";
|
|
}
|
|
cerr<<")"<<endl;
|
|
cerr<<"range: (";
|
|
for(int i=0;i<DIM;i++)
|
|
{
|
|
cerr<<bbox_here[i];
|
|
if(i<DIM-1) cerr<<" ";
|
|
}
|
|
cerr<<") - (";
|
|
for(int i=0;i<DIM;i++)
|
|
{
|
|
cerr<<bbox_here[DIM+i];
|
|
if(i<DIM-1) cerr<<" ";
|
|
}
|
|
cerr<<")"<<endl;
|
|
for(int ip=0;ip<pices;ip++)
|
|
{
|
|
cerr<<"# "<<ip<<endl;
|
|
cerr<<"shape: (";
|
|
for(int i=0;i<DIM;i++)
|
|
{
|
|
cerr<<shape_res[ip*DIM+i];
|
|
if(i<DIM-1) cerr<<" ";
|
|
}
|
|
cerr<<")"<<endl;
|
|
cerr<<"range: (";
|
|
for(int i=0;i<DIM;i++)
|
|
{
|
|
cerr<<bbox_res[ip*2*DIM+i];
|
|
if(i<DIM-1) cerr<<" ";
|
|
}
|
|
cerr<<") - (";
|
|
for(int i=0;i<DIM;i++)
|
|
{
|
|
cerr<<bbox_res[ip*2*DIM+DIM+i];
|
|
if(i<DIM-1) cerr<<" ";
|
|
}
|
|
cerr<<")"<<endl;
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
void misc::swapvector(std::vector<double> &f0, std::vector<double> &f1)
|
|
{
|
|
const int N = f0.size();
|
|
double tt;
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
tt = f0[i];
|
|
f0[i] = f1[i];
|
|
f1[i] = tt;
|
|
}
|
|
}
|
|
complex<double> misc::complex_gamma(complex<double> z)
|
|
{
|
|
const double p[9] = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
|
|
771.32342877765313, -176.61502916214059, 12.507343278686905,
|
|
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
|
|
|
|
if (real(z) < 0.5)
|
|
{
|
|
return PI / (sin(PI * z) * complex_gamma(1.0 - z));
|
|
}
|
|
z -= 1.0;
|
|
complex<double> x = p[0];
|
|
for (int i = 1; i < 9; i++)
|
|
{
|
|
x += p[i] / (z + complex<double>(i, 0));
|
|
}
|
|
complex<double> t = z + (7 + 0.5);
|
|
t = sqrt(2 * PI) * pow(t, z + 0.5) * exp(-t) * x;
|
|
|
|
return t;
|
|
}
|
|
// also called Kummer function,
|
|
// Confluent hypergeometric function 1F1
|
|
#if 1
|
|
complex<double> misc::KummerComplex(const complex<double> a, const complex<double> b, complex<double> x)
|
|
{
|
|
// Default tolerance is tol = 1e-10. Feel free to change this as needed.
|
|
const double tol = 1e-10;
|
|
|
|
// Estimates the value by summing powers of the generalized hypergeometric
|
|
// series:
|
|
//
|
|
// sum(n=0-->Inf)[(a)_n*x^n/{(b)_n*n!}]
|
|
//
|
|
// until the specified tolerance is acheived.
|
|
|
|
complex<double> term = x * a / b;
|
|
complex<double> f = 1.0 + term;
|
|
int n = 1;
|
|
complex<double> an = a;
|
|
complex<double> bn = b;
|
|
int nmin = 100000;
|
|
|
|
while (n < nmin && (abs(term)) > tol)
|
|
{
|
|
n = n + 1;
|
|
an = an + 1.0;
|
|
bn = bn + 1.0;
|
|
term = x * term * an / bn / double(n);
|
|
f = f + term;
|
|
}
|
|
|
|
if ((abs(term)) > tol && n == nmin)
|
|
cout << "misc::KummerComplex has n > " << nmin << " with error " << abs(term) << endl
|
|
<< "a = " << a << " b = " << b << " x = " << x << endl;
|
|
|
|
return f;
|
|
}
|
|
// new code
|
|
#else
|
|
complex<double> misc::KummerComplex(const complex<double> a, const complex<double> b, complex<double> z)
|
|
{
|
|
// Default tolerance is tol = 1e-10. Feel free to change this as needed.
|
|
int precision = 15;
|
|
int m, j, k;
|
|
complex<double> cr, chg;
|
|
double cMax;
|
|
complex<double> g1, g2, g3;
|
|
complex<double> ba;
|
|
complex<double> cs1, cs2, cr1, cr2;
|
|
double c1Max, c2Max;
|
|
|
|
// Special cases
|
|
|
|
if (b.imag() == 0 && b.real() <= 0 && b.real() == int(b.real())) // b==-n;n=1,2,3,..
|
|
{
|
|
if (a.imag() == 0 && a.real() <= 0 && a.real() == int(a.real()) && abs(a) < abs(b)) // a==-m;m=1,2,..
|
|
{
|
|
m = int(-a.real());
|
|
cr = 1;
|
|
chg = 1;
|
|
|
|
cMax = abs(cr);
|
|
|
|
for (k = 1; k <= m; k++)
|
|
{
|
|
cr = cr * (k - 1.0 + a) / double(k) / (k - 1.0 + b) * z;
|
|
chg = chg + cr;
|
|
|
|
cMax = max(cMax, max(abs(cr), abs(chg)));
|
|
}
|
|
|
|
precision = 15 - int(log10(cMax / abs(chg)));
|
|
}
|
|
else if (a.imag() == 0 && a.real() <= 0 && a.real() == int(a.real()) && abs(a) == abs(b)) // a==b;
|
|
{
|
|
cout << "!!!Confluent hypergeometric function is indeterminate for input a = "
|
|
<< a << " b = " << b << " z = " << z << endl;
|
|
chg = 0;
|
|
}
|
|
else
|
|
{
|
|
cout << "!!!Confluent hypergeometric function error for input a = "
|
|
<< a << " b = " << b << " z = " << z << endl;
|
|
chg = 0;
|
|
}
|
|
}
|
|
else if (a == 0.0 || z == 0.0)
|
|
{
|
|
chg = 1;
|
|
}
|
|
else if (a == -1.0)
|
|
{
|
|
chg = 1.0 - z / b;
|
|
}
|
|
else if (a == b)
|
|
{
|
|
chg = exp(z);
|
|
}
|
|
else if ((a - b) == 1.0)
|
|
{
|
|
chg = (1.0 + z / b) * exp(z);
|
|
}
|
|
else if (a == 1.0 && b == 2.0)
|
|
{
|
|
chg = (exp(z) - 1.0) / z;
|
|
}
|
|
// finite number of elements in a row
|
|
else if (a.imag() == 0 && a.real() < 0 && a.real() == int(a.real()))
|
|
{
|
|
m = int(-a.real());
|
|
cr = 1;
|
|
chg = 1;
|
|
|
|
cMax = abs(cr);
|
|
|
|
for (k = 1; k <= m; k++)
|
|
{
|
|
cr = cr * (k - 1.0 + a) / double(k) / (k - 1.0 + b) * z;
|
|
chg = chg + cr;
|
|
|
|
cMax = max(cMax, max(abs(cr), abs(chg)));
|
|
}
|
|
|
|
precision = 15 - int(log10(cMax / abs(chg)));
|
|
}
|
|
else if (abs(z) > 10 * abs(a) && abs(z) > 10 * abs(b)) // Abramowitz Stegun 13.5.1
|
|
{
|
|
g1 = complex_gamma(a);
|
|
g2 = complex_gamma(b);
|
|
ba = b - a;
|
|
g3 = complex_gamma(ba);
|
|
|
|
cs1 = 1;
|
|
cs2 = 1;
|
|
cr1 = 1;
|
|
cr2 = 1;
|
|
|
|
c1Max = abs(cr1);
|
|
c2Max = abs(cr2);
|
|
|
|
for (j = 1; j <= 500; j++)
|
|
{
|
|
cr1 = -cr1 * (j - 1.0 + a) * (a - b + double(j)) / (z * double(j));
|
|
cr2 = cr2 * (j - 1.0 + b - a) * (double(j) - a) / (z * double(j));
|
|
cs1 = cs1 + cr1;
|
|
cs2 = cs2 + cr2;
|
|
|
|
c1Max = max(c1Max, max(abs(cr1), abs(cs1)));
|
|
c2Max = max(c2Max, max(abs(cr2), abs(cs2)));
|
|
|
|
if (abs(cr1) / abs(cs1) < 1e-15 && abs(cr2) / abs(cs2) < 1e-15)
|
|
break; // break j
|
|
|
|
if (j == 500)
|
|
{
|
|
cout << "Got to the " << j << " limit in the series of confluent hypergeometric function!" << endl;
|
|
chg = 0;
|
|
return chg;
|
|
}
|
|
}
|
|
|
|
precision = 15 - int(log10(max(c1Max / abs(cs1), c2Max / abs(cs2))));
|
|
|
|
double x = z.real();
|
|
double y = z.imag();
|
|
double phi;
|
|
complex<double> cfac, chg1, chg2;
|
|
int ns;
|
|
|
|
if (x == 0.0 && y >= 0.0)
|
|
phi = 0.5 * PI;
|
|
else if (x == 0.0 && y <= 0.0)
|
|
phi = -0.5 * PI;
|
|
else
|
|
phi = atan(y / x);
|
|
|
|
if (phi > -0.5 * PI && phi < 1.5 * PI)
|
|
ns = 1;
|
|
|
|
if (phi > -1.5 * PI && phi <= -0.5 * PI)
|
|
ns = -1;
|
|
|
|
cfac = exp(PI * ns * (complex<double>(0, 1)) * a);
|
|
|
|
if (y == 0)
|
|
cfac = cos(PI * a);
|
|
|
|
chg1 = g2 / g3 * pow(z, -a) * cfac * cs1;
|
|
chg2 = g2 / g1 * exp(z) * pow(z, a - b) * cs2;
|
|
chg = chg1 + chg2;
|
|
}
|
|
else // General case
|
|
{
|
|
chg = 1;
|
|
complex<double> crg = 1;
|
|
double cgMax = abs(crg);
|
|
|
|
for (j = 1; j <= 500; j++)
|
|
{
|
|
crg = crg * (j - 1.0 + a) / (double(j) * (j - 1.0 + b)) * z; // Abramowitz Stegun 13.1.2
|
|
chg = chg + crg;
|
|
|
|
cgMax = max(cgMax, max(abs(crg), abs(chg)));
|
|
|
|
if (abs(crg) / abs(chg) < 1e-15)
|
|
break; // break j
|
|
|
|
if (j == 500)
|
|
{
|
|
cout << "Got to the " << j << " limit in the series of confluent hypergeometric function!" << endl;
|
|
chg = 0;
|
|
return chg;
|
|
}
|
|
}
|
|
|
|
precision = 15 - int(log10(cgMax / abs(chg)));
|
|
}
|
|
|
|
if (precision <= 0)
|
|
{
|
|
precision = 0;
|
|
chg = 0;
|
|
}
|
|
|
|
if (precision < 10)
|
|
cout << "!!! Warning!!! Only about " << precision << " first digits are correct!!!" << endl;
|
|
|
|
return chg;
|
|
}
|
|
#endif
|
|
// Bessel function of the first kind: J_a
|
|
#if 0
|
|
//
|
|
// sum(m=0-->Inf)(-1)^m/m!/Gamma(m+a+1) (x/2)^{2 m+a}
|
|
//
|
|
complex<double> misc::First_Bessel(const complex<double> a,complex<double> x)
|
|
{
|
|
// Default tolerance is tol = 1e-10. Feel free to change this as needed.
|
|
const double tol = 1e-10;
|
|
|
|
x = x/2.0;
|
|
complex<double> term,term1=pow(x,a),term2=1.0/complex_gamma(a+1.0);
|
|
complex<double> f = term1*term2;
|
|
int m = 0;
|
|
const int mmax = 50;
|
|
|
|
term = f;
|
|
while(m < mmax && (abs(term)) > tol)
|
|
{
|
|
m++;
|
|
term1 = x*x*term1;
|
|
term2 = -term2/double(m*m);
|
|
term = term1*term2;
|
|
f = f + term;
|
|
}
|
|
|
|
if((abs(term)) > tol && m == mmax) cout<<"misc::First_Bessel has m > "<<mmax<<" with error "<<abs(term)<<" for x = "<<x<<endl;
|
|
|
|
return f;
|
|
}
|
|
#else
|
|
complex<double> misc::First_Bessel(double a, complex<double> x)
|
|
{
|
|
double xr = x.real(), xi = x.imag();
|
|
double yr, yi;
|
|
int IERR, NZ, N = 1, KODE = 1;
|
|
f_zbesj(xr, xi, a, KODE, N, yr, yi, NZ, IERR);
|
|
complex<double> f(yr, yi);
|
|
|
|
return f;
|
|
}
|
|
#endif
|
|
complex<double> misc::Rec_Int(const double xmin, const double xmax, complex<double> fun(double x))
|
|
{
|
|
// Default tolerance is tol = 1e-10. Feel free to change this as needed.
|
|
const double tol = 1e-8;
|
|
|
|
int N = int(xmax - xmin) * 10;
|
|
if (N < 1000)
|
|
N = 1000;
|
|
double dx = (xmax - xmin) / (N - 1);
|
|
complex<double> sum = 0, sum2 = 0;
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
sum2 += fun(xmin + i * dx);
|
|
}
|
|
sum2 = sum2 * dx;
|
|
|
|
int j = 1;
|
|
const int jmax = 10;
|
|
while (j < jmax && abs(sum2 - sum) > tol)
|
|
{
|
|
j++;
|
|
N = N * 2;
|
|
dx = (xmax - xmin) / (N - 1);
|
|
sum = sum2;
|
|
sum2 = 0;
|
|
for (int i = 0; i < N; i++)
|
|
{
|
|
sum2 += fun(xmin + i * dx);
|
|
}
|
|
sum2 = sum2 * dx;
|
|
|
|
// cout<<"j = "<<j<<" error = "<<abs(sum2-sum)<<endl;
|
|
}
|
|
|
|
if (j == jmax)
|
|
cout << "misc::Rec_Int has j > " << jmax << ", error = " << abs(sum2 - sum) << endl;
|
|
|
|
return sum2;
|
|
}
|
|
complex<double> misc::Simpson_Int(const double xmin, const double xmax, complex<double> fun(double x))
|
|
{
|
|
// Default tolerance is tol = 1e-10. Feel free to change this as needed.
|
|
const double tol = 1e-8;
|
|
|
|
int N = 1000;
|
|
double dx = (xmax - xmin) / (N - 1);
|
|
complex<double> sum = 0, sum2 = 0;
|
|
sum2 = 17.0 * fun(xmin) + 59.0 * fun(xmin + dx) + 43.0 * fun(xmin + 2 * dx) + 49.0 * fun(xmin + 3 * dx);
|
|
for (int i = 4; i < N - 3; i++)
|
|
{
|
|
sum2 += 48.0 * fun(xmin + i * dx);
|
|
}
|
|
sum2 = sum2 + 17.0 * fun(xmax) + 59.0 * fun(xmax - dx) + 43.0 * fun(xmax - 2 * dx) + 49.0 * fun(xmax - 3 * dx);
|
|
sum2 = sum2 * dx / 48.0;
|
|
|
|
int j = 1;
|
|
const int jmax = 50;
|
|
while (j < jmax && abs(sum2 - sum) > tol)
|
|
{
|
|
j++;
|
|
N = N * 2;
|
|
dx = (xmax - xmin) / (N - 1);
|
|
sum = sum2;
|
|
sum2 = 17.0 * fun(xmin) + 59.0 * fun(xmin + dx) + 43.0 * fun(xmin + 2 * dx) + 49.0 * fun(xmin + 3 * dx);
|
|
for (int i = 4; i < N - 3; i++)
|
|
{
|
|
sum2 += 48.0 * fun(xmin + i * dx);
|
|
}
|
|
sum2 = sum2 + 17.0 * fun(xmax) + 59.0 * fun(xmax - dx) + 43.0 * fun(xmax - 2 * dx) + 49.0 * fun(xmax - 3 * dx);
|
|
sum2 = sum2 * dx / 48.0;
|
|
|
|
// cout<<"j = "<<j<<" error = "<<abs(sum2-sum)<<endl;
|
|
}
|
|
|
|
if (j == jmax)
|
|
cout << "misc::Simpson_Int has j > " << jmax << ", error = " << abs(sum2 - sum) << endl;
|
|
|
|
return sum2;
|
|
}
|
|
complex<double> misc::Simpson3o8_Int(const double xmin, const double xmax, complex<double> fun(double x))
|
|
{
|
|
// Default tolerance is tol = 1e-10. Feel free to change this as needed.
|
|
const double tol = 1e-8;
|
|
|
|
int m = 300, N;
|
|
N = 3 * m + 2;
|
|
double dx = (xmax - xmin) / (N - 1);
|
|
complex<double> sum = 0, sum2;
|
|
sum2 = fun(xmin) + fun(xmax);
|
|
for (int i = 0; i < m; i++)
|
|
{
|
|
sum2 += 3.0 * (fun(xmin + (3 * i + 1) * dx) + fun(xmin + (3 * i + 2) * dx)) + 2.0 * fun(xmin + (3 * i + 3) * dx);
|
|
// cout<<sum2<<endl;
|
|
// cout<<fun(xmin+(3*i+1)*dx)<<endl;
|
|
// cout<<fun(xmin+(3*i+2)*dx)<<endl;
|
|
// cout<<fun(xmin+(3*i+3)*dx)<<endl;
|
|
// if(abs(sum2) != abs(sum2)) exit(0);
|
|
}
|
|
sum2 = sum2 * dx * 3.0 / 8.0;
|
|
|
|
int j = 1;
|
|
const int jmax = 10;
|
|
while (j < jmax && abs(sum2 - sum) > tol)
|
|
{
|
|
j++;
|
|
m = m * 2;
|
|
N = 3 * m + 2;
|
|
dx = (xmax - xmin) / (N - 1);
|
|
sum = sum2;
|
|
sum2 = fun(xmin) + fun(xmax);
|
|
for (int i = 0; i < m; i++)
|
|
{
|
|
sum2 += 3.0 * (fun(xmin + (3 * i + 1) * dx) + fun(xmin + (3 * i + 2) * dx)) + 2.0 * fun(xmin + (3 * i + 3) * dx);
|
|
}
|
|
sum2 = sum2 * dx * 3.0 / 8.0;
|
|
|
|
// cout<<"j = "<<j<<" error = "<<abs(sum2-sum)<<endl;
|
|
}
|
|
|
|
if (j == jmax)
|
|
cout << "misc::Simpson3o8_Int has j > " << jmax << ", error = " << abs(sum2 - sum) << endl;
|
|
|
|
return sum2;
|
|
}
|
|
#if 0
|
|
complex<double> misc::Gauss_Int(const double xmin,const double xmax,complex<double> fun(double x))
|
|
{
|
|
// Default tolerance is tol = 1e-10. Feel free to change this as needed.
|
|
const double tol = 1e-8;
|
|
|
|
int N=int(xmax-xmin)*10;
|
|
if(N<1000) N = 1000;
|
|
double *arcostheta,*wtcostheta;
|
|
// weight function cover all of [xmin,xmax]
|
|
arcostheta = new double[N];
|
|
wtcostheta = new double[N];
|
|
|
|
gaulegf(xmin,xmax,arcostheta,wtcostheta,N);
|
|
complex<double> sum=0,sum2=0;
|
|
for(int i =0;i<N;i++)
|
|
{
|
|
sum2 += fun(arcostheta[i])*wtcostheta[i];
|
|
// cout<<sum2<<endl;
|
|
// cout<<arcostheta[i]<<","<<fun(arcostheta[i])<<endl;
|
|
// if(abs(sum2) != abs(sum2)) exit(0);
|
|
}
|
|
delete[] arcostheta;
|
|
delete[] wtcostheta;
|
|
|
|
int j=1;
|
|
const int jmax = 10;
|
|
while(j < jmax && abs(sum2-sum) > tol)
|
|
{
|
|
j++;
|
|
N = N*2;
|
|
arcostheta = new double[N];
|
|
wtcostheta = new double[N];
|
|
|
|
gaulegf(xmin,xmax,arcostheta,wtcostheta,N);
|
|
sum=sum2;
|
|
sum2=0;
|
|
for(int i =0;i<N;i++)
|
|
{
|
|
sum2 += fun(arcostheta[i])*wtcostheta[i];
|
|
}
|
|
delete[] arcostheta;
|
|
delete[] wtcostheta;
|
|
|
|
// cout<<"j = "<<j<<" error = "<<abs(sum2-sum)<<endl;
|
|
}
|
|
|
|
if(j == jmax) cout<<"misc::Gauss_Int has j > "<<jmax<<", error = "<<abs(sum2-sum)<<endl;
|
|
|
|
return sum2;
|
|
}
|
|
#else
|
|
complex<double> misc::Gauss_Int(const double xmin, const double xmax, complex<double> fun(double x))
|
|
{
|
|
// Default tolerance is tol = 1e-10. Feel free to change this as needed.
|
|
const double tol = 1e-8;
|
|
|
|
// int N=int(xmax-xmin)*10;
|
|
// if(N<1000) N = 1000;
|
|
|
|
int N = 1000;
|
|
complex<double> sum = 0, sum2 = 0;
|
|
sum2 = gaulegf(xmin, xmax, N, fun);
|
|
|
|
int j = 1;
|
|
const int jmax = 10;
|
|
while (j < jmax && abs(sum2 - sum) > tol)
|
|
{
|
|
j++;
|
|
N = N * 2;
|
|
sum = sum2;
|
|
sum2 = gaulegf(xmin, xmax, N, fun);
|
|
|
|
cout << "j = " << j << " error = " << abs(sum2 - sum) << endl;
|
|
}
|
|
|
|
// if(j == jmax)
|
|
cout << "misc::Gauss_Int has j > " << jmax << ", error = " << abs(sum2 - sum) << endl;
|
|
|
|
return sum2;
|
|
}
|
|
#endif
|
|
complex<double> misc::gaulegf(double x1, double x2, int n, complex<double> fun(double x))
|
|
{
|
|
int i, j, m;
|
|
double eps = 1.2E-16;
|
|
double p1, p2, p3, pp, xl, xm, z, z1;
|
|
double w;
|
|
|
|
m = (n + 1) / 2;
|
|
xm = 0.5 * (x2 + x1);
|
|
xl = 0.5 * (x2 - x1);
|
|
complex<double> sum = 0;
|
|
for (i = 0; i < m; i++)
|
|
{
|
|
z = cos(PI * ((double)i + 0.75) / ((double)n + 0.5));
|
|
do
|
|
{
|
|
p1 = 1.0;
|
|
p2 = 0.0;
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
p3 = p2;
|
|
p2 = p1;
|
|
p1 = ((2 * (double)j + 1) * z * p2 - (double)j * p3) / ((double)j + 1);
|
|
}
|
|
pp = n * (z * p1 - p2) / (z * z - 1.0);
|
|
z1 = z;
|
|
z = z1 - p1 / pp;
|
|
// cout<<"here"<<endl;
|
|
} while (fabs(z - z1) > eps);
|
|
|
|
// cout<<"there"<<endl;
|
|
w = 2.0 * xl / ((1.0 - z * z) * pp * pp);
|
|
sum += w * (fun(xm - xl * z) + fun(xm + xl * z));
|
|
|
|
if (!isfinite(abs(sum)))
|
|
{
|
|
cout << xm - xl * z << "," << xm + xl * z << endl;
|
|
cout << fun(xm - xl * z) << "," << fun(xm + xl * z) << endl;
|
|
}
|
|
// cout<<i<<","<<m<<endl;
|
|
}
|
|
|
|
return sum;
|
|
}
|
|
/*
|
|
This computes an in-place (in data and out data are all stored in x and y) complex-to-complex FFT
|
|
x and y are the real and imaginary arrays of 2^m points.
|
|
dir = 1 gives forward transform
|
|
X(n) = 1/N sum_{k=0} ^{N-1} x(k) exp(-jk2 pi n/N) for n=0...N-1
|
|
dir = -1 gives reverse transform
|
|
x(n) = sum_{k=0} ^{N-1} X(k) exp( jk2 pi n/N) for n=0...N-1
|
|
*/
|
|
void misc::FFT(short int dir, long m, double *x, double *y)
|
|
{
|
|
long n, i, i1, j, k, i2, l, l1, l2;
|
|
double c1, c2, tx, ty, t1, t2, u1, u2, z;
|
|
|
|
/* Calculate the number of points */
|
|
n = 1;
|
|
for (i = 0; i < m; i++)
|
|
n *= 2;
|
|
|
|
/* Do the bit reversal */
|
|
i2 = n >> 1;
|
|
j = 0;
|
|
for (i = 0; i < n - 1; i++)
|
|
{
|
|
if (i < j)
|
|
{
|
|
tx = x[i];
|
|
ty = y[i];
|
|
x[i] = x[j];
|
|
y[i] = y[j];
|
|
x[j] = tx;
|
|
y[j] = ty;
|
|
}
|
|
k = i2;
|
|
while (k <= j)
|
|
{
|
|
j -= k;
|
|
k >>= 1;
|
|
}
|
|
j += k;
|
|
}
|
|
|
|
/* Compute the FFT */
|
|
c1 = -1.0;
|
|
c2 = 0.0;
|
|
l2 = 1;
|
|
for (l = 0; l < m; l++)
|
|
{
|
|
l1 = l2;
|
|
l2 <<= 1;
|
|
u1 = 1.0;
|
|
u2 = 0.0;
|
|
for (j = 0; j < l1; j++)
|
|
{
|
|
for (i = j; i < n; i += l2)
|
|
{
|
|
i1 = i + l1;
|
|
t1 = u1 * x[i1] - u2 * y[i1];
|
|
t2 = u1 * y[i1] + u2 * x[i1];
|
|
x[i1] = x[i] - t1;
|
|
y[i1] = y[i] - t2;
|
|
x[i] += t1;
|
|
y[i] += t2;
|
|
}
|
|
z = u1 * c1 - u2 * c2;
|
|
u2 = u1 * c2 + u2 * c1;
|
|
u1 = z;
|
|
}
|
|
c2 = sqrt((1.0 - c1) / 2.0);
|
|
if (dir == 1)
|
|
c2 = -c2;
|
|
c1 = sqrt((1.0 + c1) / 2.0);
|
|
}
|
|
|
|
/* Scaling for forward transform */
|
|
if (dir == 1)
|
|
{
|
|
for (i = 0; i < n; i++)
|
|
{
|
|
x[i] /= n;
|
|
y[i] /= n;
|
|
}
|
|
}
|
|
}
|
|
// assume a[0] a[1]......a[NN/2-1] a[NN/2] ...... a[NN-1]
|
|
// 0 df (NN/2-1)*df combine of \pm NN/2*df -df
|
|
// 0 1 2 3 4 5
|
|
// ^ ^ ^ o ^ ^
|
|
// 0 1 2 3
|
|
// ^ ^ o ^
|
|
void misc::Low_Pass_Filt(const int NN, double *a)
|
|
{
|
|
// we use 2/3 law, NN/2 * 2/3 = NN/3
|
|
for (int i = 0; i < NN / 3; i++)
|
|
{
|
|
a[NN / 2 + i] = 0;
|
|
a[NN / 2 - i] = 0;
|
|
}
|
|
}
|
|
void misc::polyinterp(double t, double &rr, double *ti, double *ri, const int ORD)
|
|
{
|
|
// (x -x_1)...(x -x_i-1)(x -x_i+1)...(x -x_N)
|
|
// ------------------------------------------------f_i
|
|
// (x_i-x_1)...(x_i-x_i-1)(x_i-x_i+1)...(x_i-x_N)
|
|
|
|
rr = 0;
|
|
for (int i = 0; i < ORD; i++)
|
|
{
|
|
double ss = 1, xx = 1;
|
|
for (int j = 0; j < ORD; j++)
|
|
{
|
|
if (j != i)
|
|
{
|
|
ss *= t - ti[j];
|
|
xx *= ti[i] - ti[j];
|
|
}
|
|
}
|
|
rr += ss / xx * ri[i];
|
|
}
|
|
#if 0
|
|
if(!isfinite(rr))
|
|
{
|
|
cout.setf(ios::scientific);
|
|
cout<<"misc::polyinterp: error at t = "<<setprecision(16)<<setw(15)<<t<<endl;
|
|
for(int i=0;i<ORD;i++) cout<<setprecision(16)<<setw(15)<<ti[i]<<","<<setprecision(16)<<setw(15)<<ri[i]<<endl;
|
|
}
|
|
#endif
|
|
}
|
|
void misc::polyinterp_d1(double t, double &rr, double *ti, double *ri, const int ORD)
|
|
{
|
|
// sum_{j != i}[(x -x_1)...(x -x_j-1)(x -x_j+1)...(x -x_i-1)(x -x_i+1)...(x -x_N)]
|
|
// ----------------------------------------------------------------------------------------f_i
|
|
// (x_i-x_1)...(x_i-x_i-1)(x_i-x_i+1)...(x_i-x_N)
|
|
|
|
rr = 0;
|
|
for (int i = 0; i < ORD; i++)
|
|
{
|
|
double ss = 0, xx = 1;
|
|
for (int j = 0; j < ORD; j++)
|
|
{
|
|
if (j != i)
|
|
{
|
|
double tt = 1;
|
|
for (int k = 0; k < ORD; k++)
|
|
{
|
|
if (k != j)
|
|
tt *= t - ti[k];
|
|
}
|
|
ss += tt;
|
|
xx *= ti[i] - ti[j];
|
|
}
|
|
}
|
|
rr += ss / xx * ri[i];
|
|
}
|
|
}
|
|
void misc::next2power(long int Nin, long int &Nout, int &M)
|
|
{
|
|
M = 0;
|
|
Nout = 1;
|
|
while (Nout < Nin)
|
|
{
|
|
M++;
|
|
Nout *= 2;
|
|
}
|
|
// return Nout=2^M > Nin
|
|
}
|
|
int misc::MYpow2(int i)
|
|
{
|
|
if (i == 0)
|
|
return 1;
|
|
else if (i > 0)
|
|
return 2 * MYpow2(i - 1);
|
|
else
|
|
return MYpow2(i + 1) / 2;
|
|
}
|