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