danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Removed tabs from c_functions.c refaactored quadgr

Browse Source
master
Per A Brodtkorb 9 years ago
parent 4366d3eeb0
commit c2a9247d00
2 changed files with 467 additions and 460 deletions
Show Stats Download Patch File Download Diff File

307
wafo/integrate.py
Unescape Escape View File

@ -1003,163 +1003,170 @@ def richardson(Q, k):
return R return R
def quadgr(fun, a, b, abseps=1e-5, max_iter=17): class _Quadgr(object):
'''
Gauss-Legendre quadrature with Richardson extrapolation. def __call__(self, fun, a, b, abseps=1e-5, max_iter=17):
'''
[Q,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function Gauss-Legendre quadrature with Richardson extrapolation.
FUN from A to B with an absolute error tolerance TOL. FUN is a function
handle and must accept vector arguments. TOL is 1e-6 by default. Q is [Q,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function
the integral approximation and ERR is an estimate of the absolute error. FUN from A to B with an absolute error tolerance TOL. FUN is a function
handle and must accept vector arguments. TOL is 1e-6 by default. Q is
QUADGR uses a 12-point Gauss-Legendre quadrature. The error estimate is the integral approximation and ERR is an estimate of the absolute
based on successive interval bisection. Richardson extrapolation error.
accelerates the convergence for some integrals, especially integrals
with endpoint singularities. QUADGR uses a 12-point Gauss-Legendre quadrature. The error estimate is
based on successive interval bisection. Richardson extrapolation
Examples accelerates the convergence for some integrals, especially integrals
-------- with endpoint singularities.
>>> import numpy as np
>>> Q, err = quadgr(np.log,0,1) Examples
>>> quadgr(np.exp,0,9999*1j*np.pi) --------
(-2.0000000000122662, 2.1933237448479304e-09) >>> import numpy as np
>>> Q, err = quadgr(np.log,0,1)
>>> quadgr(lambda x: np.sqrt(4-x**2),0,2,1e-12) >>> quadgr(np.exp,0,9999*1j*np.pi)
(3.1415926535897811, 1.5809575870662229e-13) (-2.0000000000122662, 2.1933237448479304e-09)
>>> quadgr(lambda x: x**-0.75,0,1) >>> quadgr(lambda x: np.sqrt(4-x**2),0,2,1e-12)
(4.0000000000000266, 5.6843418860808015e-14) (3.1415926535897811, 1.5809575870662229e-13)
>>> quadgr(lambda x: 1./np.sqrt(1-x**2),-1,1) >>> quadgr(lambda x: x**-0.75,0,1)
(3.141596056985029, 6.2146261559092864e-06) (4.0000000000000266, 5.6843418860808015e-14)
>>> quadgr(lambda x: np.exp(-x**2),-np.inf,np.inf,1e-9) #% sqrt(pi) >>> quadgr(lambda x: 1./np.sqrt(1-x**2),-1,1)
(1.7724538509055152, 1.9722334876348668e-11) (3.141596056985029, 6.2146261559092864e-06)
>>> quadgr(lambda x: np.cos(x)*np.exp(-x),0,np.inf,1e-9) >>> quadgr(lambda x: np.exp(-x**2),-np.inf,np.inf,1e-9) #% sqrt(pi)
(0.50000000000000044, 7.3296813063450372e-11) (1.7724538509055152, 1.9722334876348668e-11)
>>> quadgr(lambda x: np.cos(x)*np.exp(-x),0,np.inf,1e-9)
(0.50000000000000044, 7.3296813063450372e-11)
See also
--------
QUAD,
QUADGK
'''
# Author: jonas.lundgren@saabgroup.com, 2009. license BSD
# Order limits (required if infinite limits)
a = np.asarray(a)
b = np.asarray(b)
if a == b:
Q = b - a
err = b - a
return Q, err
elif np.real(a) > np.real(b):
reverse = True
a, b = b, a
else:
reverse = False
# Infinite limits
if np.isinf(a) | np.isinf(b):
# Check real limits
if ~ np.isreal(a) | ~np.isreal(b) | np.isnan(a) | np.isnan(b):
raise ValueError('Infinite intervals must be real.')
# Change of variable
if np.isfinite(a) & np.isinf(b):
# a to inf
[Q, err] = quadgr(lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2,
0, 1, abseps)
elif np.isinf(a) & np.isfinite(b):
# -inf to b
[Q, err] = quadgr(lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2,
-1, 0, abseps)
else: # -inf to inf
[Q1, err1] = quadgr(lambda t: fun(t / (1 - t)) / (1 - t) ** 2,
0, 1, abseps / 2)
[Q2, err2] = quadgr(lambda t: fun(t / (1 + t)) / (1 + t) ** 2,
-1, 0, abseps / 2)
Q = Q1 + Q2
err = err1 + err2
# Reverse direction
if reverse:
Q = -Q
return Q, err
# Gauss-Legendre quadrature (12-point)
xq = np.asarray(
[0.12523340851146894, 0.36783149899818018, 0.58731795428661748,
0.76990267419430469, 0.9041172563704748, 0.98156063424671924])
wq = np.asarray(
[0.24914704581340288, 0.23349253653835478, 0.20316742672306584,
0.16007832854334636, 0.10693932599531818, 0.047175336386511842])
xq = np.hstack((xq, -xq))
wq = np.hstack((wq, wq))
nq = len(xq)
dtype = np.result_type(fun(a), fun(b))
# Initiate vectors
Q0 = zeros(max_iter, dtype=dtype) # Quadrature
Q1 = zeros(max_iter, dtype=dtype) # First Richardson extrapolation
Q2 = zeros(max_iter, dtype=dtype) # Second Richardson extrapolation
# One interval
hh = (b - a) / 2 # Half interval length
x = (a + b) / 2 + hh * xq # Nodes
# Quadrature
Q0[0] = hh * np.sum(wq * fun(x), axis=0)
# Successive bisection of intervals
for k in range(1, max_iter):
# Interval bisection
hh = hh / 2
x = np.hstack([x + a, x + b]) / 2
# Quadrature
Q0[k] = hh * np.sum(wq * np.sum(np.reshape(fun(x), (-1, nq)),
axis=0),
axis=0)
# Richardson extrapolation
if k >= 5:
Q1[k] = richardson(Q0, k)
Q2[k] = richardson(Q1, k)
elif k >= 3:
Q1[k] = richardson(Q0, k)
# Estimate absolute error
if k >= 6:
Qv = np.hstack((Q0[k], Q1[k], Q2[k]))
Qw = np.hstack((Q0[k - 1], Q1[k - 1], Q2[k - 1]))
elif k >= 4:
Qv = np.hstack((Q0[k], Q1[k]))
Qw = np.hstack((Q0[k - 1], Q1[k - 1]))
else:
Qv = np.atleast_1d(Q0[k])
Qw = Q0[k - 1]
errors = np.atleast_1d(abs(Qv - Qw))
j = errors.argmin()
err = errors[j]
Q = Qv[j]
if k >= 2: # and not iscomplex:
_val, err1 = dea3(Q0[k - 2], Q0[k - 1], Q0[k])
# Convergence
if (err < abseps) | ~np.isfinite(Q):
break
else:
warnings.warn('Max number of iterations reached without ' +
'convergence.')
See also if ~ np.isfinite(Q):
-------- warnings.warn('Integral approximation is Infinite or NaN.')
QUAD,
QUADGK
'''
# Author: jonas.lundgren@saabgroup.com, 2009. license BSD
# Order limits (required if infinite limits)
a = np.asarray(a)
b = np.asarray(b)
if a == b:
Q = b - a
err = b - a
return Q, err
elif np.real(a) > np.real(b):
reverse = True
a, b = b, a
else:
reverse = False
# Infinite limits
if np.isinf(a) | np.isinf(b):
# Check real limits
if ~ np.isreal(a) | ~np.isreal(b) | np.isnan(a) | np.isnan(b):
raise ValueError('Infinite intervals must be real.')
# Change of variable
if np.isfinite(a) & np.isinf(b):
# a to inf
[Q, err] = quadgr(lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2,
0, 1, abseps)
elif np.isinf(a) & np.isfinite(b):
# -inf to b
[Q, err] = quadgr(lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2,
-1, 0, abseps)
else: # -inf to inf
[Q1, err1] = quadgr(lambda t: fun(t / (1 - t)) / (1 - t) ** 2,
0, 1, abseps / 2)
[Q2, err2] = quadgr(lambda t: fun(t / (1 + t)) / (1 + t) ** 2,
-1, 0, abseps / 2)
Q = Q1 + Q2
err = err1 + err2
# The error estimate should not be zero
err = err + 2 * np.finfo(Q).eps
# Reverse direction # Reverse direction
if reverse: if reverse:
Q = -Q Q = -Q
return Q, err
# Gauss-Legendre quadrature (12-point) return Q, err
xq = np.asarray(
[0.12523340851146894, 0.36783149899818018, 0.58731795428661748,
0.76990267419430469, 0.9041172563704748, 0.98156063424671924])
wq = np.asarray(
[0.24914704581340288, 0.23349253653835478, 0.20316742672306584,
0.16007832854334636, 0.10693932599531818, 0.047175336386511842])
xq = np.hstack((xq, -xq))
wq = np.hstack((wq, wq))
nq = len(xq)
dtype = np.result_type(fun(a), fun(b))
# Initiate vectors
Q0 = zeros(max_iter, dtype=dtype) # Quadrature
Q1 = zeros(max_iter, dtype=dtype) # First Richardson extrapolation
Q2 = zeros(max_iter, dtype=dtype) # Second Richardson extrapolation
# One interval
hh = (b - a) / 2 # Half interval length
x = (a + b) / 2 + hh * xq # Nodes
# Quadrature
Q0[0] = hh * np.sum(wq * fun(x), axis=0)
# Successive bisection of intervals
for k in range(1, max_iter):
# Interval bisection
hh = hh / 2
x = np.hstack([x + a, x + b]) / 2
# Quadrature
Q0[k] = hh * np.sum(wq * np.sum(np.reshape(fun(x), (-1, nq)), axis=0),
axis=0)
# Richardson extrapolation
if k >= 5:
Q1[k] = richardson(Q0, k)
Q2[k] = richardson(Q1, k)
elif k >= 3:
Q1[k] = richardson(Q0, k)
# Estimate absolute error
if k >= 6:
Qv = np.hstack((Q0[k], Q1[k], Q2[k]))
Qw = np.hstack((Q0[k - 1], Q1[k - 1], Q2[k - 1]))
elif k >= 4:
Qv = np.hstack((Q0[k], Q1[k]))
Qw = np.hstack((Q0[k - 1], Q1[k - 1]))
else:
Qv = np.atleast_1d(Q0[k])
Qw = Q0[k - 1]
errors = np.atleast_1d(abs(Qv - Qw))
j = errors.argmin()
err = errors[j]
Q = Qv[j]
if k >= 2: # and not iscomplex:
_val, err1 = dea3(Q0[k - 2], Q0[k - 1], Q0[k])
# Convergence
if (err < abseps) | ~np.isfinite(Q):
break
else:
warnings.warn('Max number of iterations reached without convergence.')
if ~ np.isfinite(Q):
warnings.warn('Integral approximation is Infinite or NaN.')
# The error estimate should not be zero
err = err + 2 * np.finfo(Q).eps
# Reverse direction
if reverse:
Q = -Q
return Q, err quadgr = _Quadgr()
def boole(y, x): def boole(y, x):

480
wafo/source/c_library/c_functions.c
Unescape Escape View File

@ -35,7 +35,7 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) {
if (( *(y+0) > *(y+1)) && ( *(y+1) > *(y+2)) ){ if (( *(y+0) > *(y+1)) && ( *(y+1) > *(y+2)) ){
info[0] = -1; info[0] = -1;
return; /*This is not a sequence of turningpoints, exit */ return; /*This is not a sequence of turningpoints, exit */
} }
if ((*(y+0) < *(y+1)) && (*(y+1)< *(y+2))){ if ((*(y+0) < *(y+1)) && (*(y+1)< *(y+2))){
@ -52,61 +52,61 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) {
xplus=*(y+2*i+2); xplus=*(y+2*i+2);
if(i!=0){ if(i!=0){
j=i-1; j=i-1;
while((j>=0) && (*(y+2*j+1)<=*(y+2*i+1))){ while((j>=0) && (*(y+2*j+1)<=*(y+2*i+1))){
if( (*(y+2*j)<xminus) ){ if( (*(y+2*j)<xminus) ){
xminus=*(y+2*j); xminus=*(y+2*j);
Tmi=Tstart+2*j; Tmi=Tstart+2*j;
} /*if */ } /*if */
j--; j--;
} /*while j*/ } /*while j*/
} /*if i */ } /*if i */
if ( xminus >= xplus){ if ( xminus >= xplus){
if ( (*(y+2*i+1)-xminus) >= hmin){ if ( (*(y+2*i+1)-xminus) >= hmin){
*(ind+ix)=Tmi; *(ind+ix)=Tmi;
ix++; ix++;
*(ind+ix)=(Tstart+2*i+1); *(ind+ix)=(Tstart+2*i+1);
ix++; ix++;
} /*if*/ } /*if*/
goto L180; goto L180;
} }
j=i+1; j=i+1;
while((j<NC) ) { while((j<NC) ) {
if (*(y+2*j+1) >= *(y+2*i+1)) goto L170; if (*(y+2*j+1) >= *(y+2*i+1)) goto L170;
if( (*(y+2*j+2) <= xplus) ){ if( (*(y+2*j+2) <= xplus) ){
xplus=*(y+2*j+2); xplus=*(y+2*j+2);
Tpl=(Tstart+2*j+2); Tpl=(Tstart+2*j+2);
}/*if*/ }/*if*/
j++; j++;
} /*while*/ } /*while*/
if ( (*(y+2*i+1)-xminus) >= hmin) { if ( (*(y+2*i+1)-xminus) >= hmin) {
*(ind+ix)=Tmi; *(ind+ix)=Tmi;
ix++; ix++;
*(ind+ix)=(Tstart+2*i+1); *(ind+ix)=(Tstart+2*i+1);
ix++; ix++;
} /*if*/ } /*if*/
goto L180; goto L180;
L170: L170:
if (xplus <= xminus ) { if (xplus <= xminus ) {
if ( (*(y+2*i+1)-xminus) >= hmin){ if ( (*(y+2*i+1)-xminus) >= hmin){
*(ind+ix)=Tmi; *(ind+ix)=Tmi;
ix++; ix++;
*(ind+ix)=(Tstart+2*i+1); *(ind+ix)=(Tstart+2*i+1);
ix++; ix++;
} /*if*/ } /*if*/
/*goto L180;*/ /*goto L180;*/
} }
else{ else{
if ( (*(y+2*i+1)-xplus) >= hmin) { if ( (*(y+2*i+1)-xplus) >= hmin) {
*(ind+ix)=(Tstart+2*i+1); *(ind+ix)=(Tstart+2*i+1);
ix++; ix++;
*(ind+ix)=Tpl; *(ind+ix)=Tpl;
ix++; ix++;
} /*if*/ } /*if*/
} /*elseif*/ } /*elseif*/
L180: L180:
iy=i; iy=i;
@ -128,42 +128,42 @@ void findrfc(double *y1,double hmin, int *ind, int n,int *info) {
void findcross(double *y, double v, int *ind, int n, int *info) void findcross(double *y, double v, int *ind, int n, int *info)
{ int i,start, ix=0,dcross=0; { int i,start, ix=0,dcross=0;
start=0; start=0;
if ( y[0]< v){ if ( y[0]< v){
dcross=-1; /* first is a up-crossing*/ dcross=-1; /* first is a up-crossing*/
} }
else if ( y[0]> v){ else if ( y[0]> v){
dcross=1; /* first is a down-crossing*/ dcross=1; /* first is a down-crossing*/
} }
else if ( y[0]== v){ else if ( y[0]== v){
/* Find out what type of crossing we have next time.. */ /* Find out what type of crossing we have next time.. */
for (i=1; i<n; i++) { for (i=1; i<n; i++) {
start=i; start=i;
if ( y[i]< v){ if ( y[i]< v){
ind[ix] = i-1; /* first crossing is a down crossing*/ ind[ix] = i-1; /* first crossing is a down crossing*/
ix++; ix++;
dcross=-1; /* The next crossing is a up-crossing*/ dcross=-1; /* The next crossing is a up-crossing*/
goto L120; goto L120;
} }
else if ( y[i]> v){ else if ( y[i]> v){
ind[ix] = i-1; /* first crossing is a up-crossing*/ ind[ix] = i-1; /* first crossing is a up-crossing*/
ix++; ix++;
dcross=1; /*The next crossing is a down-crossing*/ dcross=1; /*The next crossing is a down-crossing*/
goto L120; goto L120;
} }
} }
} }
L120: L120:
for (i=start; i<n-1; i++) { for (i=start; i<n-1; i++) {
if (( (dcross==-1) && (y[i]<=v) && (y[i+1] > v) ) || ((dcross==1 ) && (y[i]>=v) && (y[i+1] < v) ) ) { if (( (dcross==-1) && (y[i]<=v) && (y[i+1] > v) ) || ((dcross==1 ) && (y[i]>=v) && (y[i+1] < v) ) ) {
ind[ix] = i; ind[ix] = i;
ix++; ix++;
dcross=-dcross; dcross=-dcross;
} }
} }
info[0] = ix; info[0] = ix;
return; return;
} }
@ -198,11 +198,11 @@ void findcross(double *y, double v, int *ind, int n, int *info)
*/ */
void disufq(double *rvec, double *ivec, void disufq(double *rvec, double *ivec,
double *rA, double *iA, double *rA, double *iA,
double *w, double *kw, double *w, double *kw,
double h, double g, double h, double g,
int nmin, int nmax, int nmin, int nmax,
int m, int n) int m, int n)
{ {
double Epij, Edij; double Epij, Edij;
double tmp1, tmp2, tmp3, tmp4, kfact; double tmp1, tmp2, tmp3, tmp4, kfact;
@ -229,62 +229,62 @@ void disufq(double *rvec, double *ivec,
kw1 = kw[ix]; kw1 = kw[ix];
Epij = kw1; Epij = kw1;
for (i=0;i<m;i++,ixi++,iz1++) { for (i=0;i<m;i++,ixi++,iz1++) {
rrA = rA[ixi]*rA[ixi]; /// rrA = rA[ixi]*rA[ixi]; ///
iiA = iA[ixi]*iA[ixi]; /// iiA = iA[ixi]*iA[ixi]; ///
riA = rA[ixi]*iA[ixi]; /// riA = rA[ixi]*iA[ixi]; ///
/// Sum frequency effects along the diagonal /// Sum frequency effects along the diagonal
tmp1 = kfact*(rrA-iiA)*Epij; tmp1 = kfact*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*riA*Epij; tmp2 = kfact*2.0*riA*Epij;
rvec[iz1] += tmp1; rvec[iz1] += tmp1;
ivec[iz1] += tmp2; ivec[iz1] += tmp2;
//rvec[iz2] += tmp1; //rvec[iz2] += tmp1;
//ivec[iz2] -= tmp2; //ivec[iz2] -= tmp2;
//iz2++; //iz2++;
/// Difference frequency effects are zero along the diagonal /// Difference frequency effects are zero along the diagonal
/// and are thus not contributing to the mean. /// and are thus not contributing to the mean.
} }
for (jy = ix+1;jy<nmax;jy++){ for (jy = ix+1;jy<nmax;jy++){
kw2 = kw[jy]; kw2 = kw[jy];
Epij = 0.5*(kw2 + kw1); Epij = 0.5*(kw2 + kw1);
Edij = -0.5*(kw2 - kw1); Edij = -0.5*(kw2 - kw1);
//printf("Edij = %f Epij = %f \n", Edij,Epij); //printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix*m; ixi = ix*m;
jyi = jy*m; jyi = jy*m;
iz1 = ixi+jyi; iz1 = ixi+jyi;
iv1 = jyi-ixi; iv1 = jyi-ixi;
//iz2 = (n*m-iz1); //iz2 = (n*m-iz1);
//iv2 = (n*m-iv1); //iv2 = (n*m-iv1);
for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1; rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1;
ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2; ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2;
//rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1; //rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1;
//ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -= tmp2; //ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -= tmp2;
// iz2++; // iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1; rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2; ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
//rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1; //rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
//ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2; //ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
//iv2++; //iv2++;
} }
} }
} }
} }
@ -357,29 +357,29 @@ void disufq(double *rvec, double *ivec,
// iz2 = (n*m-iz1); // iz2 = (n*m-iz1);
// iv2 = n*m-iv1; // iv2 = n*m-iv1;
for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rvec[iz1] += tmp1;///rvec[i][jy+ix] += tmp1; rvec[iz1] += tmp1;///rvec[i][jy+ix] += tmp1;
ivec[iz1] += tmp2;///ivec[i][jy+ix] += tmp2; ivec[iz1] += tmp2;///ivec[i][jy+ix] += tmp2;
//rvec[iz2] += tmp1;///rvec[i][n*m-(jy+ix)] += tmp1; //rvec[iz2] += tmp1;///rvec[i][n*m-(jy+ix)] += tmp1;
//ivec[iz2] -= tmp2;///ivec[i][n*m-(jy+ix)] -= tmp2; //ivec[iz2] -= tmp2;///ivec[i][n*m-(jy+ix)] -= tmp2;
//iz2++; //iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1; rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
ivec[iv1] += tmp2;///ivec[i][jy-ix] -= tmp2; ivec[iv1] += tmp2;///ivec[i][jy-ix] -= tmp2;
//rvec[iv2] += tmp1; //rvec[iv2] += tmp1;
//ivec[iv2] -= tmp2; //ivec[iv2] -= tmp2;
//iv2++; //iv2++;
} }
} }
} }
@ -420,12 +420,12 @@ void disufq(double *rvec, double *ivec,
*/ */
void disufq2(double *rsvec, double *isvec, void disufq2(double *rsvec, double *isvec,
double *rdvec, double *idvec, double *rdvec, double *idvec,
double *rA, double *iA, double *rA, double *iA,
double *w, double *kw, double *w, double *kw,
double h, double g, double h, double g,
int nmin, int nmax, int nmin, int nmax,
int m, int n) int m, int n)
{ {
double Epij, Edij; double Epij, Edij;
double tmp1, tmp2, tmp3, tmp4, kfact; double tmp1, tmp2, tmp3, tmp4, kfact;
@ -455,62 +455,62 @@ void disufq2(double *rsvec, double *isvec,
kw1 = kw[ix]; kw1 = kw[ix];
Epij = kw1; Epij = kw1;
for (i=0;i<m;i++,ixi++,iz1++) { for (i=0;i<m;i++,ixi++,iz1++) {
rrA = rA[ixi]*rA[ixi]; /// rrA = rA[ixi]*rA[ixi]; ///
iiA = iA[ixi]*iA[ixi]; /// iiA = iA[ixi]*iA[ixi]; ///
riA = rA[ixi]*iA[ixi]; /// riA = rA[ixi]*iA[ixi]; ///
/// Sum frequency effects along the diagonal /// Sum frequency effects along the diagonal
tmp1 = kfact*(rrA-iiA)*Epij; tmp1 = kfact*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*riA*Epij; tmp2 = kfact*2.0*riA*Epij;
rsvec[iz1] += tmp1; rsvec[iz1] += tmp1;
isvec[iz1] += tmp2; isvec[iz1] += tmp2;
//rsvec[iz2] += tmp1; //rsvec[iz2] += tmp1;
//isvec[iz2] -= tmp2; //isvec[iz2] -= tmp2;
//iz2++; //iz2++;
/// Difference frequency effects are zero along the diagonal /// Difference frequency effects are zero along the diagonal
/// and are thus not contributing to the mean. /// and are thus not contributing to the mean.
} }
for (jy = ix+1;jy<nmax;jy++){ for (jy = ix+1;jy<nmax;jy++){
kw2 = kw[jy]; kw2 = kw[jy];
Epij = 0.5*(kw2 + kw1); Epij = 0.5*(kw2 + kw1);
Edij = -0.5*(kw2 - kw1); Edij = -0.5*(kw2 - kw1);
//printf("Edij = %f Epij = %f \n", Edij,Epij); //printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix*m; ixi = ix*m;
jyi = jy*m; jyi = jy*m;
iz1 = ixi+jyi; iz1 = ixi+jyi;
iv1 = jyi-ixi; iv1 = jyi-ixi;
//iz2 = (n*m-iz1); //iz2 = (n*m-iz1);
//iv2 = (n*m-iv1); //iv2 = (n*m-iv1);
for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rsvec[iz1] += tmp1; ///rvec[i][ix+jy] += tmp1; rsvec[iz1] += tmp1; ///rvec[i][ix+jy] += tmp1;
isvec[iz1] += tmp2; ///ivec[i][ix+jy] += tmp2; isvec[iz1] += tmp2; ///ivec[i][ix+jy] += tmp2;
//rsvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1; //rsvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1;
//isvec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] += tmp2; //isvec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] += tmp2;
//iz2++; //iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rdvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1; rdvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
idvec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2; idvec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
//rdvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1; //rdvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
//idvec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2; //idvec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
// iv2++; // iv2++;
} }
} }
} }
} }
@ -583,29 +583,29 @@ void disufq2(double *rsvec, double *isvec,
// iz2 = (n*m-iz1); // iz2 = (n*m-iz1);
// iv2 = (n*m-iv1); // iv2 = (n*m-iv1);
for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rsvec[iz1] += tmp1;///rsvec[i][jy+ix] += tmp1; rsvec[iz1] += tmp1;///rsvec[i][jy+ix] += tmp1;
isvec[iz1] += tmp2;///isvec[i][jy+ix] += tmp2; isvec[iz1] += tmp2;///isvec[i][jy+ix] += tmp2;
//rsvec[iz2] += tmp1;///rsvec[i][n*m-(jy+ix)] += tmp1; //rsvec[iz2] += tmp1;///rsvec[i][n*m-(jy+ix)] += tmp1;
//isvec[iz2] -= tmp2;///isvec[i][n*m-(jy-ix)] += tmp2; //isvec[iz2] -= tmp2;///isvec[i][n*m-(jy-ix)] += tmp2;
//iz2++; //iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rdvec[iv1] += tmp1; rdvec[iv1] += tmp1;
idvec[iv1] += tmp2; idvec[iv1] += tmp2;
//rdvec[iv2] += tmp1; //rdvec[iv2] += tmp1;
//idvec[iv2] -= tmp2; //idvec[iv2] -= tmp2;
// iv2++; // iv2++;
} }
} }
} }

Write Preview
Loading…
Cancel
Save