14 #define PI 3.141592653589793
15 #define speedlight 299792458.0
20 inline double Lagrange(
const int n,
const int i,
const double x)
23 double xn = x+n/2.-0.5;
25 for(
int j=0;
j<
n;
j++)
26 if(
j!=i) c *= (xn-
j)/(i-
j);
37 template<
class DataType_t>
51 q[i] = p[i] + (x--)*(p[i+1]-p[
i]);
56 q[0] += xm*(x--)*(q[1]-q[0]);
58 q[1] += xm*(x--)*(q[2]-q[1]);
60 q[2] += xm*(x--)*(q[3]-q[2]);
62 q[3] += xm*(x--)*(q[4]-q[3]);
64 q[4] += xm*(x--)*(q[5]-q[4]);
66 q[5] += xm*(x--)*(q[6]-q[5]);
70 q[i] += xm*(x--)*(q[i+1]-q[
i]);
79 inline double signPDF(
const size_t m,
const size_t k)
84 double rho = (double(m)-double(k))/n;
88 for(i=1; i<k+1; i++){ pdf -=
log(
double(m+i)/
double(i)); }
89 pdf -=
log(
double(n))-(n+1)*
log(2.);
90 pdf -=
log(sqrt(2.*
PI/n));
92 else pdf = n*rho*rho/2.;
106 for(
int k=1;
k<
n;
k++){
108 if((y*=Y/(
k+1))>1.e290)
break;
114 return x>=n ? x-n + (2+(n-1)/sqrt(2.))/(n+1)-(n-sqrt(n/4.)-1)*
log(x/n) :
115 pow(x/n,sqrt(n))*(2+(n-1)/sqrt(2.))/(n+1);
126 inline double logNormCL(
double Y,
double p=0.,
double s=1.,
double a=0.0001){
128 double z =
a*sqrt(
log(4.));
132 return erfc(z/sqrt(2.))/2.;
143 inline double logNormArg(
double CL,
double p=0.,
double s=1.,
double a=0.){
144 double z =
fabs(sqrt(-2.14*
log(2.*CL))-0.506);
145 if(
a <= 0.)
return p+z*
s;
146 double f =
a*sqrt(
log(4.));
147 z = (exp((z+
a)*
a)-1)/a;
148 return p+z*s*f/sinh(f);
161 for(i=0; i<
n; i++) c[i]=0.;
164 if(abs(n/2-i)<=m/2)
continue;
167 if(abs(n/2-j)<=m/2)
continue;
168 if(j!=i) c[
i] *= double(n/2-j)/(i-
j);
200 double p0 = 1.000000000190015;
201 double p1 = 76.18009172947146;
202 double p2 = -86.50532032941677;
203 double p3 = 24.01409824083091;
204 double p4 = -1.231739572450155;
205 double p5 = 1.208650973866179e-3;
206 double p6 = -5.395239384953e-6;
207 double g = p0+p1/(r+1)+p2/(r+2)+p3/(r+3)+p4/(r+4)+p5/(r+5)+p6/(r+6);
208 return g*pow(r+5.5,r+0.5)*exp(-(r+5.5))*sqrt(2*
PI)/
r;
215 double a = r>0 ? 1/r : 0;
219 for(
int i=1;
i<
n;
i++) {
232 while(P>p) { x *= 0.999; P=
Gamma1G(r,x); }
double gammaCL(double x, double n)
static double g(double e)
wavearray< double > a(hp.size())
double iGamma(double r, double p)
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
double logNormCL(double Y, double p=0., double s=1., double a=0.0001)
double Nevill(const double x0, int n, DataType_t *p, double *q)
watplot p1(const_cast< char * >("TFMap1"))
double signPDF(const size_t m, const size_t k)
double Lagrange(const int n, const int i, const double x)
double logNormArg(double CL, double p=0., double s=1., double a=0.)
void fLagrange(int n, int m, double *c)
double gammaCLa(double Y, int n)
watplot p2(const_cast< char * >("TFMap2"))
double fabs(const Complex &x)
double iGamma1G(double r, double p)