Commit 06c5a73b authored by thomas.forbriger's avatar thomas.forbriger

libseife [TASK]: handle case of nstep=1 in C-version

parent ddf18bee
......@@ -7,8 +7,6 @@
* This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/
*
* $Id$
*
* ----
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -25,19 +23,18 @@
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
* ----
*
* \date 28/06/2005
* \date 22/06/2020
*
* remove tides (implementation)
*
* REVISIONS and CHANGES
* - 28/06/2005 V1.0 Thomas Forbriger
* - 22/06/2020 V1.1 handle case of nstep=1 explicitly
*
* ============================================================================
*/
#define TF_CSEIFE_TIDES_C_VERSION \
"TF_CSEIFE_TIDES_C V1.0 "
#define TF_CSEIFE_TIDES_C_CVSID \
"$Id$"
"TF_CSEIFE_TIDES_C V1.1"
#include <cseife.h>
#include <math.h>
......@@ -62,19 +59,24 @@ void seife_tid(double* x, int n, double dt, int nstep)
/* System generated locals */
int i__1, i__2, i__3, i__4;
/* Builtin functions */
// double atan(double), cos(double), sin(double);
/* Local variables */
double tdif, omeg[6];
int ndim;
double step, tint, omeg0, xgez1, xgez2, xgez3, a[169] /*
static double tdif, omeg[6];
static int ndim;
static double step, tint, omeg0, xgez1, xgez2, xgez3, a[169] /*
was [13][13] */, c__[13], d__[13], e[13], f[13];
int i__, j, k;
double t;
int nfreq, k2, k3;
double tstep2;
int jj;
double rs[13], sx, dth;
int nco;
double cor[6], dur;
static int i__, j, k;
static double t;
static int nfreq, k2, k3;
static double tstep2;
static int jj;
static double rs[13], sx, dth;
static int nco;
static double cor[6], dur;
// extern /* Subroutine */ int seife_gauss__(doublereal *, integer *,
// integer *, doublereal *, doublereal *);
/* remove tides. number of frequencies is automatically chosen according */
/* to the total length of the record. */
......@@ -92,6 +94,11 @@ void seife_tid(double* x, int n, double dt, int nstep)
step = (double) nstep;
tstep2 = (step - one) * dth;
tint = tstep2 + dth;
if (nstep == 1) {
step = 1.;
tstep2 = 0.;
tint = dth;
}
/* determine the number of frequencies required for a good fit */
dur = n * dt / 3600.;
nfreq = 6;
......@@ -124,84 +131,140 @@ void seife_tid(double* x, int n, double dt, int nstep)
a[i__ + k * 13 - 14] = 0.;
}
}
if (nstep == 1) {
/* ====================================================================== */
/* do not average */
/* ============== */
/* correction for averaging over nstep samples */
i__2 = nfreq;
for (j = 1; j <= i__2; ++j) {
cor[j - 1] = 1.;
}
/* set up system of linear equations */
c__[0] = one;
i__2 = n;
for (j = 1; j <= i__2; ++j) {
sx = zero;
i__1 = j;
for (jj = j; jj <= i__1; ++jj) {
sx += x[jj];
}
t = (j - 1) * dt;
i__1 = nfreq;
for (k = 1; k <= i__1; ++k) {
c__[(k << 1) - 1] = cos(omeg[k - 1] * t);
c__[k * 2] = sin(omeg[k - 1] * t);
}
i__1 = nco;
for (i__ = 1; i__ <= i__1; ++i__) {
rs[i__ - 1] += sx * c__[i__ - 1];
i__3 = nco;
for (k = 1; k <= i__3; ++k) {
a[i__ + k * 13 - 14] += c__[i__ - 1] * c__[k - 1];
}
}
}
/* solve for partial amplitudes */
seife_gauss(a, nco, ndim, rs, f);
i__2 = n;
for (j = 1; j <= i__2; ++j) {
t = (j - 1) * dt;
/* remove average and tides */
xgez1 = f[0];
i__1 = nfreq;
for (k = 1; k <= i__1; ++k) {
k2 = k << 1;
k3 = k2 + 1;
c__[k2 - 1] = cos(omeg[k - 1] * t);
c__[k3 - 1] = sin(omeg[k - 1] * t);
xgez1 = xgez1 + f[k2 - 1] * c__[k2 - 1] + f[k3 - 1] * c__[k3
- 1];
}
x[j] -= xgez1;
}
} else {
/* ====================================================================== */
/* average over nstep samples */
/* ========================== */
/* correction for averaging over nstep samples */
i__2 = nfreq;
for (j = 1; j <= i__2; ++j) {
i__2 = nfreq;
for (j = 1; j <= i__2; ++j) {
/* L102: */
cor[j - 1] = step * sin(omeg[j - 1] * dth) / sin(step * omeg[j - 1] *
dth);
}
cor[j - 1] = step * sin(omeg[j - 1] * dth) / sin(step * omeg[j -
1] * dth);
}
/* set up system of linear equations */
c__[0] = one;
i__2 = n - nstep + 1;
i__1 = nstep;
for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
sx = zero;
i__3 = j + nstep - 1;
for (jj = j; jj <= i__3; ++jj) {
c__[0] = one;
i__2 = n - nstep + 1;
i__1 = nstep;
for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
sx = zero;
i__3 = j + nstep - 1;
for (jj = j; jj <= i__3; ++jj) {
/* L12: */
sx += x[jj];
}
sx /= step;
t = (j - 1) * dt + tstep2;
i__3 = nfreq;
for (k = 1; k <= i__3; ++k) {
c__[(k << 1) - 1] = cos(omeg[k - 1] * t);
sx += x[jj];
}
sx /= step;
t = (j - 1) * dt + tstep2;
i__3 = nfreq;
for (k = 1; k <= i__3; ++k) {
c__[(k << 1) - 1] = cos(omeg[k - 1] * t);
/* L103: */
c__[k * 2] = sin(omeg[k - 1] * t);
}
i__3 = nco;
for (i__ = 1; i__ <= i__3; ++i__) {
rs[i__ - 1] += sx * c__[i__ - 1];
i__4 = nco;
for (k = 1; k <= i__4; ++k) {
c__[k * 2] = sin(omeg[k - 1] * t);
}
i__3 = nco;
for (i__ = 1; i__ <= i__3; ++i__) {
rs[i__ - 1] += sx * c__[i__ - 1];
i__4 = nco;
for (k = 1; k <= i__4; ++k) {
/* L2: */
a[i__ + k * 13 - 14] += c__[i__ - 1] * c__[k - 1];
a[i__ + k * 13 - 14] += c__[i__ - 1] * c__[k - 1];
}
}
}
}
/* solve for partial amplitudes */
seife_gauss(a, nco, ndim, rs, f);
i__4 = n - nstep + 1;
i__3 = nstep;
for (j = 1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
t = (j - 1) * dt + tstep2;
seife_gauss(a, nco, ndim, rs, f);
i__4 = n - nstep + 1;
i__3 = nstep;
for (j = 1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
t = (j - 1) * dt + tstep2;
/* remove average and tides */
i__1 = nfreq;
for (k = 1; k <= i__1; ++k) {
k2 = k << 1;
k3 = k2 + 1;
c__[k2 - 1] = cos(omeg[k - 1] * t);
c__[k3 - 1] = sin(omeg[k - 1] * t);
d__[k2 - 1] = -omeg[k - 1] * c__[k3 - 1];
d__[k3 - 1] = omeg[k - 1] * c__[k2 - 1];
e[k2 - 1] = -omeg[k - 1] * d__[k3 - 1];
i__1 = nfreq;
for (k = 1; k <= i__1; ++k) {
k2 = k << 1;
k3 = k2 + 1;
c__[k2 - 1] = cos(omeg[k - 1] * t);
c__[k3 - 1] = sin(omeg[k - 1] * t);
d__[k2 - 1] = -omeg[k - 1] * c__[k3 - 1];
d__[k3 - 1] = omeg[k - 1] * c__[k2 - 1];
e[k2 - 1] = -omeg[k - 1] * d__[k3 - 1];
/* L104: */
e[k3 - 1] = omeg[k - 1] * d__[k2 - 1];
}
xgez1 = f[0];
xgez2 = zero;
xgez3 = zero;
i__1 = nfreq;
for (k = 1; k <= i__1; ++k) {
k2 = k << 1;
k3 = k2 + 1;
xgez1 += cor[k - 1] * (f[k2 - 1] * c__[k2 - 1] + f[k3 - 1] * c__[
k3 - 1]);
xgez2 += cor[k - 1] * (f[k2 - 1] * d__[k2 - 1] + f[k3 - 1] * d__[
k3 - 1]);
e[k3 - 1] = omeg[k - 1] * d__[k2 - 1];
}
xgez1 = f[0];
xgez2 = zero;
xgez3 = zero;
i__1 = nfreq;
for (k = 1; k <= i__1; ++k) {
k2 = k << 1;
k3 = k2 + 1;
xgez1 += cor[k - 1] * (f[k2 - 1] * c__[k2 - 1] + f[k3 - 1] *
c__[k3 - 1]);
xgez2 += cor[k - 1] * (f[k2 - 1] * d__[k2 - 1] + f[k3 - 1] *
d__[k3 - 1]);
/* L105: */
xgez3 += cor[k - 1] * (f[k2 - 1] * e[k2 - 1] + f[k3 - 1] * e[k3 -
1]) / two;
}
if (j > n - (nstep << 1) + 1) {
nstep = n + 1 - j;
}
i__1 = j + nstep - 1;
for (jj = j; jj <= i__1; ++jj) {
tdif = (jj - j) * dt - tstep2;
xgez3 += cor[k - 1] * (f[k2 - 1] * e[k2 - 1] + f[k3 - 1] * e[
k3 - 1]) / two;
}
if (j > n - (nstep << 1) + 1) {
nstep = n + 1 - j;
}
i__1 = j + nstep - 1;
for (jj = j; jj <= i__1; ++jj) {
tdif = (jj - j) * dt - tstep2;
/* L3: */
x[jj] = x[jj] - xgez1 - xgez2 * tdif - xgez3 * tdif * tdif;
x[jj] = x[jj] - xgez1 - xgez2 * tdif - xgez3 * tdif * tdif;
}
}
}
} /* seife_tides__ */
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment