Commit 775f627d authored by thomas.forbriger's avatar thomas.forbriger
Browse files

libs/libseife [FEATURE]: add M3 to tide removal in C-version

revision of tide removal code:
- both versions (C and Fortran): do case of nstep=1 explicitly
  seife by default averages over serval samples and uses the tidal
  hamonic as well as its derivatives for corrections; this appears to
  have saved computation time in the very early times but is of no
  relevance nowadays; by not averaging the code (and the computation)
  becomes much more straight-forward and less error-prone
- C-version: add M3 frequency
  The number of frequencies was fixed to 6 in the original version; to
  allow for the addition of frequencies, the dimensions of arrays now
  are det by C preprocessor macros; the function seife_gauss is modified
  accordingly to allow for a larger system of equations
parents e73b9e98 e69dcd9c
this is <COPYING> this is <COPYING>
============================================================================ ============================================================================
libseife libseife
--------
$Id$
============================================================================ ============================================================================
The source code in this directory is part of libseife which compiles to The source code in this directory is part of libseife which compiles to
......
# this is <Makefile> # this is <Makefile>
# ---------------------------------------------------------------------------- # ----------------------------------------------------------------------------
# ($Id$)
# #
# 25/10/2000 by Thomas Forbriger (IMGF Frankfurt) # 25/10/2000 by Thomas Forbriger (IMGF Frankfurt)
# #
......
...@@ -3,8 +3,6 @@ ...@@ -3,8 +3,6 @@
* *
* ---------------------------------------------------------------------------- * ----------------------------------------------------------------------------
* *
* $Id$
*
* ---- * ----
* This program is free software; you can redistribute it and/or modify * 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 * it under the terms of the GNU General Public License as published by
...@@ -37,8 +35,6 @@ ...@@ -37,8 +35,6 @@
*/ */
#define TF_CSEIFE_C_VERSION \ #define TF_CSEIFE_C_VERSION \
"TF_CSEIFE_C V1.2" "TF_CSEIFE_C V1.2"
#define TF_CSEIFE_C_CVSID \
"$Id$"
#include <cseife.h> #include <cseife.h>
#include <stdio.h> #include <stdio.h>
......
...@@ -3,7 +3,6 @@ ...@@ -3,7 +3,6 @@
* *
* ---------------------------------------------------------------------------- * ----------------------------------------------------------------------------
* *
* $Id$
* \author Thomas Forbriger * \author Thomas Forbriger
* \date 14/01/2005 * \date 14/01/2005
* *
...@@ -40,8 +39,6 @@ ...@@ -40,8 +39,6 @@
#define TF_CSEIFE_H_VERSION \ #define TF_CSEIFE_H_VERSION \
"TF_CSEIFE_H V1.2" "TF_CSEIFE_H V1.2"
#define TF_CSEIFE_H_CVSID \
"$Id$"
#define SEIFE_EXIT_FAILURE 1 #define SEIFE_EXIT_FAILURE 1
......
...@@ -7,8 +7,6 @@ ...@@ -7,8 +7,6 @@
* This algorithm was part of seife.f. A current version of seife.f can be * This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/ * obtained from http://www.software-for-seismometry.de/
* *
* $Id$
*
* ---- * ----
* This program is free software; you can redistribute it and/or modify * 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 * it under the terms of the GNU General Public License as published by
...@@ -36,8 +34,6 @@ ...@@ -36,8 +34,6 @@
*/ */
#define TF_CSEIFE_DERIV_C_VERSION \ #define TF_CSEIFE_DERIV_C_VERSION \
"TF_CSEIFE_DERIV_C V1.0 " "TF_CSEIFE_DERIV_C V1.0 "
#define TF_CSEIFE_DERIV_C_CVSID \
"$Id$"
#include <cseife.h> #include <cseife.h>
......
...@@ -7,8 +7,6 @@ ...@@ -7,8 +7,6 @@
* This algorithm was part of seife.f. A current version of seife.f can be * This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/ * obtained from http://www.software-for-seismometry.de/
* *
* $Id$
*
* ---- * ----
* This program is free software; you can redistribute it and/or modify * 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 * it under the terms of the GNU General Public License as published by
...@@ -32,13 +30,13 @@ ...@@ -32,13 +30,13 @@
* REVISIONS and CHANGES * REVISIONS and CHANGES
* - 28/06/2005 V1.0 Thomas Forbriger * - 28/06/2005 V1.0 Thomas Forbriger
* - 15/11/2010 V1.1 do not use tfmacros.h * - 15/11/2010 V1.1 do not use tfmacros.h
* - 22/06/2020 V1.2 increase array size to allow for more tidal
* frequencies
* *
* ============================================================================ * ============================================================================
*/ */
#define TF_CSEIFE_GAUSS_C_VERSION \ #define TF_CSEIFE_GAUSS_C_VERSION \
"TF_CSEIFE_GAUSS_C V1.1" "TF_CSEIFE_GAUSS_C V1.2"
#define TF_CSEIFE_GAUSS_C_CVSID \
"$Id$"
#include <cseife.h> #include <cseife.h>
#include <stdio.h> #include <stdio.h>
...@@ -51,6 +49,7 @@ ...@@ -51,6 +49,7 @@
the code was derived through f2c, but modified thereafter the code was derived through f2c, but modified thereafter
*/ */
#define C_MSIZE 20
void seife_gauss(double *aik, int m, int n, double* rs, double* f) void seife_gauss(double *aik, int m, int n, double* rs, double* f)
{ {
/* System generated locals */ /* System generated locals */
...@@ -58,15 +57,15 @@ void seife_gauss(double *aik, int m, int n, double* rs, double* f) ...@@ -58,15 +57,15 @@ void seife_gauss(double *aik, int m, int n, double* rs, double* f)
double d__1; double d__1;
/* Local variables */ /* Local variables */
int imax[14]; int imax[(C_MSIZE+1)];
double h__[15]; double h__[(C_MSIZE+2)];
int j, k, l; int j, k, l;
double q; double q;
int index; int index;
double aikmax; double aikmax;
SEIFE_CHECKERROR( m>13 , "seife_gauss", "matrix is too large" ) SEIFE_CHECKERROR( m>C_MSIZE , "seife_gauss", "matrix is too large" )
SEIFE_CHECKERROR( n>13 , "seife_gauss", "matrix is too large" ) SEIFE_CHECKERROR( n>C_MSIZE , "seife_gauss", "matrix is too large" )
/* solve linear equations */ /* solve linear equations */
/* Parameter adjustments */ /* Parameter adjustments */
--f; --f;
......
...@@ -7,8 +7,6 @@ ...@@ -7,8 +7,6 @@
* This algorithm was part of seife.f. A current version of seife.f can be * This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/ * obtained from http://www.software-for-seismometry.de/
* *
* $Id$
*
* ---- * ----
* This program is free software; you can redistribute it and/or modify * 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 * it under the terms of the GNU General Public License as published by
...@@ -37,8 +35,6 @@ ...@@ -37,8 +35,6 @@
*/ */
#define TF_CSEIFE_REKFL_C_VERSION \ #define TF_CSEIFE_REKFL_C_VERSION \
"TF_CSEIFE_REKFL_C V1.1" "TF_CSEIFE_REKFL_C V1.1"
#define TF_CSEIFE_REKFL_C_CVSID \
"$Id$"
#include <stdio.h> #include <stdio.h>
#include <cseife.h> #include <cseife.h>
......
...@@ -7,8 +7,6 @@ ...@@ -7,8 +7,6 @@
* This algorithm was part of seife.f. A current version of seife.f can be * This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/ * obtained from http://www.software-for-seismometry.de/
* *
* $Id$
*
* ---- * ----
* This program is free software; you can redistribute it and/or modify * 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 * it under the terms of the GNU General Public License as published by
...@@ -36,8 +34,6 @@ ...@@ -36,8 +34,6 @@
*/ */
#define TF_CSEIFE_RFK_C_VERSION \ #define TF_CSEIFE_RFK_C_VERSION \
"TF_CSEIFE_RFK_C V1.0 " "TF_CSEIFE_RFK_C V1.0 "
#define TF_CSEIFE_RFK_C_CVSID \
"$Id$"
#include <cseife.h> #include <cseife.h>
#include <math.h> #include <math.h>
......
...@@ -7,8 +7,6 @@ ...@@ -7,8 +7,6 @@
* This algorithm was part of seife.f. A current version of seife.f can be * This algorithm was part of seife.f. A current version of seife.f can be
* obtained from http://www.software-for-seismometry.de/ * obtained from http://www.software-for-seismometry.de/
* *
* $Id$
*
* ---- * ----
* This program is free software; you can redistribute it and/or modify * 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 * it under the terms of the GNU General Public License as published by
...@@ -25,19 +23,19 @@ ...@@ -25,19 +23,19 @@
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
* ---- * ----
* *
* \date 28/06/2005 * \date 22/06/2020
* *
* remove tides (implementation) * remove tides (implementation)
* *
* REVISIONS and CHANGES * REVISIONS and CHANGES
* - 28/06/2005 V1.0 Thomas Forbriger * - 28/06/2005 V1.0 Thomas Forbriger
* - 22/06/2020 V1.1 handle case of nstep=1 explicitly
* V1.2 make array size flexible and add M3 frequency
* *
* ============================================================================ * ============================================================================
*/ */
#define TF_CSEIFE_TIDES_C_VERSION \ #define TF_CSEIFE_TIDES_C_VERSION \
"TF_CSEIFE_TIDES_C V1.0 " "TF_CSEIFE_TIDES_C V1.1"
#define TF_CSEIFE_TIDES_C_CVSID \
"$Id$"
#include <cseife.h> #include <cseife.h>
#include <math.h> #include <math.h>
...@@ -54,7 +52,10 @@ void seife_tid(double* x, int n, double dt, int nstep) ...@@ -54,7 +52,10 @@ void seife_tid(double* x, int n, double dt, int nstep)
{ {
/* Initialized data */ /* Initialized data */
static double omega[6] = { 1.93227,.92954,2.,1.00274,1.89567,.89293 }; #define C_MFREQ 7
#define C_MDIM (C_MFREQ*2+1)
static double omega[C_MFREQ] = { 2.898450424,
1.93227,.92954,2.,1.00274,1.89567,.89293 };
static double zero = 0.; static double zero = 0.;
static double one = 1.; static double one = 1.;
static double two = 2.; static double two = 2.;
...@@ -62,19 +63,25 @@ void seife_tid(double* x, int n, double dt, int nstep) ...@@ -62,19 +63,25 @@ void seife_tid(double* x, int n, double dt, int nstep)
/* System generated locals */ /* System generated locals */
int i__1, i__2, i__3, i__4; int i__1, i__2, i__3, i__4;
/* Builtin functions */
// double atan(double), cos(double), sin(double);
/* Local variables */ /* Local variables */
double tdif, omeg[6]; static double tdif, omeg[C_MFREQ];
int ndim; static int ndim;
double step, tint, omeg0, xgez1, xgez2, xgez3, a[169] /* static double step, tint, omeg0, xgez1, xgez2, xgez3,
was [13][13] */, c__[13], d__[13], e[13], f[13]; a[C_MDIM*C_MDIM], c__[C_MDIM], d__[C_MDIM],
int i__, j, k; e[C_MDIM], f[C_MDIM];
double t; static int i__, j, k;
int nfreq, k2, k3; static double t;
double tstep2; static int nfreq, k2, k3;
int jj; static double tstep2;
double rs[13], sx, dth; static int jj;
int nco; static double rs[C_MDIM], sx, dth;
double cor[6], dur; static int nco;
static double cor[C_MFREQ], dur;
// extern /* Subroutine */ int seife_gauss__(doublereal *, integer *,
// integer *, doublereal *, doublereal *);
/* remove tides. number of frequencies is automatically chosen according */ /* remove tides. number of frequencies is automatically chosen according */
/* to the total length of the record. */ /* to the total length of the record. */
...@@ -82,7 +89,7 @@ void seife_tid(double* x, int n, double dt, int nstep) ...@@ -82,7 +89,7 @@ void seife_tid(double* x, int n, double dt, int nstep)
--x; --x;
/* Function Body */ /* Function Body */
ndim = 13; ndim = C_MDIM;
dth = dt / two; dth = dt / two;
if (nstep == 0) { if (nstep == 0) {
...@@ -92,20 +99,25 @@ void seife_tid(double* x, int n, double dt, int nstep) ...@@ -92,20 +99,25 @@ void seife_tid(double* x, int n, double dt, int nstep)
step = (double) nstep; step = (double) nstep;
tstep2 = (step - one) * dth; tstep2 = (step - one) * dth;
tint = tstep2 + dth; tint = tstep2 + dth;
if (nstep == 1) {
step = 1.;
tstep2 = 0.;
tint = dth;
}
/* determine the number of frequencies required for a good fit */ /* determine the number of frequencies required for a good fit */
dur = n * dt / 3600.; dur = n * dt / 3600.;
nfreq = 6; nfreq = C_MFREQ;
if (dur < 35.) { if (dur < 35.) {
nfreq = 5; nfreq = 6;
} }
if (dur < 18.) { if (dur < 18.) {
nfreq = 4; nfreq = 5;
} }
if (dur < 14.) { if (dur < 14.) {
nfreq = 3; nfreq = 4;
} }
if (dur < 5.) { if (dur < 5.) {
nfreq = 2; nfreq = 3;
} }
nco = (nfreq << 1) + 1; nco = (nfreq << 1) + 1;
omeg0 = atan(one) * 8. / 86400.; omeg0 = atan(one) * 8. / 86400.;
...@@ -121,87 +133,143 @@ void seife_tid(double* x, int n, double dt, int nstep) ...@@ -121,87 +133,143 @@ void seife_tid(double* x, int n, double dt, int nstep)
i__2 = nco; i__2 = nco;
for (k = 1; k <= i__2; ++k) { for (k = 1; k <= i__2; ++k) {
/* L1: */ /* L1: */
a[i__ + k * 13 - 14] = 0.; a[i__ + k * C_MDIM - (C_MDIM+1)] = 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 * C_MDIM - (C_MDIM+1)] += 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 */ /* correction for averaging over nstep samples */
i__2 = nfreq; i__2 = nfreq;
for (j = 1; j <= i__2; ++j) { for (j = 1; j <= i__2; ++j) {
/* L102: */ /* L102: */
cor[j - 1] = step * sin(omeg[j - 1] * dth) / sin(step * omeg[j - 1] * cor[j - 1] = step * sin(omeg[j - 1] * dth) / sin(step * omeg[j -
dth); 1] * dth);
} }
/* set up system of linear equations */ /* set up system of linear equations */
c__[0] = one; c__[0] = one;
i__2 = n - nstep + 1; i__2 = n - nstep + 1;
i__1 = nstep; i__1 = nstep;
for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
sx = zero; sx = zero;
i__3 = j + nstep - 1; i__3 = j + nstep - 1;
for (jj = j; jj <= i__3; ++jj) { for (jj = j; jj <= i__3; ++jj) {
/* L12: */ /* L12: */
sx += x[jj]; sx += x[jj];
} }
sx /= step; sx /= step;
t = (j - 1) * dt + tstep2; t = (j - 1) * dt + tstep2;
i__3 = nfreq; i__3 = nfreq;
for (k = 1; k <= i__3; ++k) { for (k = 1; k <= i__3; ++k) {
c__[(k << 1) - 1] = cos(omeg[k - 1] * t); c__[(k << 1) - 1] = cos(omeg[k - 1] * t);
/* L103: */ /* L103: */
c__[k * 2] = sin(omeg[k - 1] * t); c__[k * 2] = sin(omeg[k - 1] * t);
} }
i__3 = nco; i__3 = nco;
for (i__ = 1; i__ <= i__3; ++i__) { for (i__ = 1; i__ <= i__3; ++i__) {
rs[i__ - 1] += sx * c__[i__ - 1]; rs[i__ - 1] += sx * c__[i__ - 1];
i__4 = nco; i__4 = nco;
for (k = 1; k <= i__4; ++k) { for (k = 1; k <= i__4; ++k) {
/* L2: */ /* L2: */
a[i__ + k * 13 - 14] += c__[i__ - 1] * c__[k - 1]; a[i__ + k * C_MDIM - (C_MDIM+1)] += c__[i__ - 1] * c__[k - 1];
}
} }
} }
}
/* solve for partial amplitudes */ /* solve for partial amplitudes */
seife_gauss(a, nco, ndim, rs, f); seife_gauss(a, nco, ndim, rs, f);
i__4 = n - nstep + 1; i__4 = n - nstep + 1;
i__3 = nstep; i__3 = nstep;
for (j = 1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) { for (j = 1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
t = (j - 1) * dt + tstep2; t = (j - 1) * dt + tstep2;
/* remove average and tides */ /* remove average and tides */
i__1 = nfreq; i__1 = nfreq;
for (k = 1; k <= i__1; ++k) { for (k = 1; k <= i__1; ++k) {
k2 = k << 1; k2 = k << 1;
k3 = k2 + 1; k3 = k2 + 1;
c__[k2 - 1] = cos(omeg[k - 1] * t); c__[k2 - 1] = cos(omeg[k - 1] * t);
c__[k3 - 1] = sin(omeg[k - 1] * t); c__[k3 - 1] = sin(omeg[k - 1] * t);
d__[k2 - 1] = -omeg[k - 1] * c__[k3 - 1]; d__[k2 - 1] = -omeg[k - 1] * c__[k3 - 1];
d__[k3 - 1] = omeg[k - 1] * c__[k2 - 1]; d__[k3 - 1] = omeg[k - 1] * c__[k2 - 1];
e[k2 - 1] = -omeg[k - 1] * d__[k3 - 1]; e[k2 - 1] = -omeg[k - 1] * d__[k3 - 1];
/* L104: */ /* L104: */
e[k3 - 1] = omeg[k - 1] * d__[k2 - 1]; e[k3 - 1] = omeg[k - 1] * d__[k2 - 1];
} }
xgez1 = f[0]; xgez1 = f[0];
xgez2 = zero; xgez2 = zero;
xgez3 = zero; xgez3 = zero;
i__1 = nfreq; i__1 = nfreq;
for (k = 1; k <= i__1; ++k) { for (k = 1; k <= i__1; ++k) {
k2 = k << 1; k2 = k << 1;
k3 = k2 + 1; k3 = k2 + 1;
xgez1 += cor[k - 1] * (f[k2 - 1] * c__[k2 - 1] + f[k3 - 1] * c__[ xgez1 += cor[k - 1] * (f[k2 - 1] * c__[k2 - 1] + f[k3 - 1] *
k3 - 1]); c__[k3 - 1]);
xgez2 += cor[k - 1] * (f[k2 - 1] * d__[k2 - 1] + f[k3 - 1] * d__[ xgez2 += cor[k - 1] * (f[k2 - 1] * d__[k2 - 1] + f[k3 - 1] *
k3 - 1]); d__[k3 - 1]);
/* L105: */ /* L105: */
xgez3 += cor[k - 1] * (f[k2 - 1] * e[k2 - 1] + f[k3 - 1] * e[k3 - xgez3 += cor[k - 1] * (f[k2 - 1] * e[k2 - 1] + f[k3 - 1] * e[
1]) / two; k3 - 1]) / two;
} }
if (j > n - (nstep << 1) + 1) { if (j > n - (nstep << 1) + 1) {
nstep = n + 1 - j; nstep = n + 1 - j;
} }
i__1 = j + nstep - 1; i__1 = j + nstep - 1;