Commit 785dc80a authored by thomas.forbriger's avatar thomas.forbriger

[MERGE] (newgr|master): merge master into branch newgr

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.

The master branch has developped significantly since the laste merge.
parents 7250918f ff082569
......@@ -53,10 +53,23 @@ junk*
*.TT
*.su
*.sff
#
# temporary files
*.tmp
*.htmp
#
# automatic prototypes
*.P
#
# secondary LaTeX output
*.aux
*.dvi
*.log
*.toc
#
# ignore all source code files, which can automatically be created
# and which store user specific settings
# ----------------------------------------------------------------
src/green/gremlin1/libs/glq_dim.inc
#
#
# ----- END OF .gitignore -----
#!/bin/sh
# this is <DL1sum.sh>
# ----------------------------------------------------------------------------
VERSION=2014-05-09
VERSION=2015-01-07
#
# Copyright (c) 2014 by Thomas Forbriger (BFO Schiltach)
#
......@@ -26,6 +26,7 @@ VERSION=2014-05-09
#
# REVISIONS and CHANGES
# 09/05/2014 V1.0 Thomas Forbriger
# 07/01/2015 implement option -y|--year
#
# ============================================================================
#
......@@ -80,8 +81,8 @@ verbose() {
# ----
# read command line options
# -------------------------
TEMP=$(getopt -o hb:t:Dv --long \
help,basepath:,tmpfile:,verbose,debug \
TEMP=$(getopt -o hb:t:y:Dv --long \
help,year:,basepath:,tmpfile:,verbose,debug \
-n $(basename $0) -- "$@") || {
echo >&2
echo >&2 "ERROR: command line parameters are incorrect!"
......@@ -114,6 +115,7 @@ while true; do
--help|-h) usage; echo; longusage; exit 1;;
--) shift; break;;
-b|--basepath) BASEPATH=$2; shift;;
-y|--year) YEAR=$2; shift;;
-t|--tmpfile) TMPFILE=$2; shift;;
-v|--verbose) VERBOSE=1;;
-D) set -x ;;
......
this is <COPYING>
============================================================================
green/grepg
---------
-----------
$Id$
============================================================================
Display Fourier Bessel expansion coefficients
Copyright 1997,2010 by Thomas Forbriger (IfG Stuttgart)
Copyright (c) 1997, 2010 by Thomas Forbriger (IfG Stuttgart)
Copyright (c) 1983-2001 by the California Institute of Technology
The source code and other files in this directory are part of grepg.
Most components of grepg are provided under terms of the
GNU General Public License
See README.PGPLOT.copyright.notice for conditions applying to
grepg_phasewedg.f
----
This program is free software; you can redistribute it and/or modify
......
......@@ -6,22 +6,18 @@ Display Fourier Bessel expansion coefficients
$Id$
============================================================================
This directory contains the following programs:
This directory contains the following program:
grepg: display Fourier Bessel expansion coefficients
like created by greda or syg
libraries/packages needed to compile the code
f2c: Fortran to C converter from netlib
PGPLOT: libpgplot.a and g77 version libpgplot77.a
they are called libf2cpgplot52.a and libpgplot52.a here
source code can be obtained from
PGPLOT: source code can be obtained from
www.astro.caltech.edu/~tjp/pgplot/
libtf: package containing libtf.a
you can obtain the latter package from where you obtained the present code
============================================================================
Installation:
......@@ -42,8 +38,6 @@ libraries and library header files are expected:
LOCLIBDIR defines location of binary libraries
LOCINCLUDEDIR defines location of C/C++ header files (prototypes)
LOCBINDIR defines location of binary executables
============================================================================
......
......@@ -6,7 +6,7 @@ c
c Copyright (c) 1983-2001 by the California Institute of Technology
c Copyright (c) 1999 by Thomas Forbriger (IfG Stuttgart)
c
c This is a copy of code written by Tim Pearson for the PGPLOT library
c This is a modfied copy of code written by Tim Pearson for the PGPLOT library
c For the license conditions see
c http://www.astro.caltech.edu/~tjp/pgplot
c or README.PGPLOT.copyright.notice
......@@ -22,7 +22,7 @@ c I: plot wedge for full color phase and amplitude plot by using the HUE
c value to represent the phase angle
c G: use selfmade color map to represent angle values from 0. to 360. deg.
c
c They main purpose of this routine is to:
c The main purpose of this modification of the original subroutine is to:
c 1. provide the full color phase angle wedge
c 2. to create a wedge with axis labels suitable for the range from 0. to
c 360. (i.e. stepwidth of major tick marks is 45.)
......
......@@ -33,6 +33,18 @@
* - 31/10/2012 V1.1 distinguish warning and critical warning; the latter
* for duplicate sample values being at variance
* - 07/11/2012 V1.2 add sample values for duplicate entries
* - 07/01/2015 V1.3
* modifications applied to File::put() function:
* - accept duplicate time stamps as part of normal
* operation; this is consistent with the behaviour
* of DL1logger as implemented on 20/03/2014; it is
* further consistent with the tested behaviour of
* the ThiesDL1 logger (see note below).
* - add sample values for duplicate time stamps
* - tr-mode (input stream tolerates redundant samples)
* now suppresses a notice-message which otherwise
* would be output to the terminal in cases of
* duplicate time stamps
*
* ============================================================================
* A statement of Volker König at ThiesClima regarding the duplicate entries
......@@ -48,9 +60,9 @@
* ============================================================================
*/
#define DATRW_THIESDL1FILE_CC_VERSION \
"DATRW_THIESDL1FILE_CC V1.2"
"DATRW_THIESDL1FILE_CC V1.3"
#define DATRW_THIESDL1FILE_CC_CVSID \
"$Id$"
"2015-01-07"
#include <datrwxx/thiesdl1file.h>
#include <datrwxx/thiesdl1line.h>
......@@ -348,47 +360,20 @@ namespace datrw {
}
if (i>=0 && i<Mnsamples)
{
int countvalue=line.counts();
if (Mfilled(i))
{
oss.clear();
oss.str("");
oss << "WARNING: duplicate sample (index "
oss << "NOTICE: duplicate sample time (index "
<< i << "): " << line.line();
std::cerr << oss.str() << std::endl;
Mtracefree.append(oss.str());
if (countvalue != Miseries(i))
{
oss.clear();
oss.str("");
oss << "CRITICAL WARNING: counts for duplicate sample differ: ";
std::cerr << oss.str() << std::endl;
Mtracefree.append(oss.str());
}
if ((countvalue > 0) || (countvalue != Miseries(i)))
if (!this->Mbetolerantagainstredundant)
{
oss.clear();
oss.str("");
oss << " previous value: " << Miseries(i)
<< " this value: " << countvalue;
if (countvalue > 0)
{
oss << " (adding this to previous)";
}
std::cerr << oss.str() << std::endl;
Mtracefree.append(oss.str());
}
}
if (!this->Mbetolerantagainstredundant)
{
DL1_rcassert (!Mfilled(i),
"duplicate sample",
this->Mheader.earliestdate,
this->Mheader.latestdate,
line.line());
Mtracefree.append(oss.str());
}
Mfilled(i)=true;
Miseries(i)+=countvalue;
Miseries(i) += line.counts();
}
else
{
......
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libemod
# -------------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libfouries
# ----------------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libgrrefsub
# -----------------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libsff
# ------------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libsffu
# -------------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libstfinv
# ---------------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libtf
# -----------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libtime
# -------------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libts
# -----------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
# this is <.gitignore>
# ============================================================================
# files to be ignored by git
# --------------------------
# files in src/libs/libttsynt
# ---------------------------
#
# source code documentation
*.doc
#
# ----- END OF .gitignore -----
......@@ -35,6 +35,7 @@ c
c REVISIONS and CHANGES
c 06/06/2005 V1.0 Thomas Forbriger
c 26/06/2013 V1.1 check definition of PSD values
c 04/12/2014 V1.2 fix: declare return type of function fnlnm
c
c ============================================================================
c
......@@ -42,14 +43,14 @@ c
c
character*(*) version
parameter(version=
& 'NLNMTAB V1.1 write table of NLNM values')
& 'NLNMTAB V1.2 write table of NLNM values')
character*(*) NLNMTAB_CVS_ID
parameter(NLNMTAB_CVS_ID=
& '$Id$')
c
c
logical useperiod, uselog
double precision f,p,f1,f2,df
double precision f,p,f1,f2,df,fnlnm
integer k,n
character*80 outstring
c commandline
......
......@@ -48,6 +48,11 @@
* and segment (analysis) series
* - 10/01/2011 V1.8b corrected Fourier transformation formula in
* doxygen documentation
* - 03/12/2014 V1.9 provide output file format selector
* - 08/01/2015 V1.10 FIX: when using -scalerbw take bandwidth from
* opt.scaledecades rather than opt.decades
* - 29/01/2015 V1.11 FEATURE: provide calculation of n-th derivative
* of time series
*
* \note 08/01/2010:
* Scaling for foutra power spectrum was tested against theory.
......@@ -65,7 +70,7 @@
* ============================================================================
*/
#define FOUTRA_VERSION \
"FOUTRA V1.8 Fourier transforms"
"FOUTRA V1.10 Fourier transforms"
#define FOUTRA_CVSID \
"$Id$"
......@@ -89,8 +94,8 @@
#include <tfxx/xcmdline.h>
#include <tfxx/misc.h>
#include <tfxx/handle.h>
#include <sffostream.h>
#include <datrwxx/readany.h>
#include <datrwxx/writeany.h>
#include <fourier/fftwaff.h>
#include <sstream>
......@@ -100,12 +105,13 @@ using std::endl;
struct Options {
bool verbose, overwrite, debug, asciiout, logascii;
std::string inputformat, asciibase;
std::string inputformat, asciibase, outputformat;
bool amplitudespectrum, powerspectrum, boxcartaper;
bool avgconstbw, avgrelbw, avgasciionly;
bool demean, detrend, scalerbw, adjustdivisor;
bool derivative;
int ndemean, ndetrend, divisor;
double decades, scaledecades, asciidecades;
double decades, scaledecades, asciidecades, nderivative;
int avgsamples;
bool reportrms, harmonicsignal, padzeroes;
int padfactor, nsegments;
......@@ -134,11 +140,12 @@ int main(int iargc, char* argv[])
char usage_text[]=
{
FOUTRA_VERSION "\n"
"usage: foutra [-v] [-o] [-type type] [-D] [-ASCII[=base]]" "\n"
"usage: foutra [-v] [-o] [-type type] [-Type type] [-D] [-ASCII[=base]]\n"
" [-amplitude] [-power] [-logascii[=n]]" "\n"
" [-boxcar] [-avg[=n]] [-rbw[=n]] [-avgascii]" "\n"
" [-demean[=n]] [-dtrend[=n]] [-scalerbw[=n]]" "\n"
" [-divisor[=n]] [-rms] [-harmonic] [-pad n]" "\n"
" [-derivative n]" "\n"
" [-nsegments n]" "\n"
" outfile infile [t:T] [infile [t:T] ...]" "\n"
" or: foutra --help|-h" "\n"
......@@ -165,11 +172,12 @@ int main(int iargc, char* argv[])
"-amplitude calculate amplitude spectrum" "\n"
"-power calculate power spectrum" "\n"
"-type type select input file type" "\n"
"-type type select input file type" "\n"
"-Type type select output file type" "\n"
"-avg[=n] smooth power spectrum by averaging over n samples" "\n"
"-rbw[=n] smooth power spectrum by averaging over n decades" "\n"
"-demean[=n] remove average (determined from n samples)" "\n"
"-detrend[=n] remove trend (determined from n samples)" "\n"
"-derivative[=n] take n-th derivative of time series" "\n"
"-scalerbw[=n] scale to mean value in n decades" "\n"
"-divisor[=n] FFT becomes very inefficient if the factorization" "\n"
" of the number of samples includes large prime numbers." "\n"
......@@ -285,6 +293,10 @@ int main(int iargc, char* argv[])
{"nsegments",arg_yes,"1"},
// 21: extended help
{"xhelp",arg_no,"-"},
// 22: input format
{"Type",arg_yes,"sff"},
// 23: input format
{"derivative",arg_opt,"1"},
{NULL}
};
......@@ -308,7 +320,7 @@ int main(int iargc, char* argv[])
{
cerr << usage_text << endl;
cerr << help_text << endl;
datrw::supported_input_data_types(cerr);
datrw::supported_data_types(cerr);
if (cmdline.optset(21)) { datrw::online_help(cerr); }
exit(0);
}
......@@ -344,6 +356,9 @@ int main(int iargc, char* argv[])
opt.padzeroes=cmdline.optset(19);
opt.padfactor=cmdline.int_arg(19);
opt.nsegments=cmdline.int_arg(20);
opt.outputformat=cmdline.string_arg(22);
opt.derivative=cmdline.optset(23);
opt.nderivative=cmdline.double_arg(23);
TFXX_assert((opt.divisor > 0), "illegal value for argument divisor");
TFXX_assert((opt.nsegments > 0), "illegal value for argument nsegments");
......@@ -428,6 +443,15 @@ int main(int iargc, char* argv[])
if (opt.harmonicsignal || opt.amplitudespectrum || opt.powerspectrum)
{
processfree.append("An appropriately scaled FFT (libdrfftw) is applied");
if (opt.derivative)
{
processfree.append("Calculate values for derivative with respect to");
processfree.append(" time by multiplication of the Fourier coefficients");
std::ostringstream freeline;
freeline << " with the angular frequency to the power of "
<< opt.nderivative;
processfree.append(freeline.str());
}
}
if (opt.amplitudespectrum)
{
......@@ -520,7 +544,7 @@ int main(int iargc, char* argv[])
TFXX_assert((!file.good()),"ERROR: output file exists!");
}
std::ofstream ofs(outfile.c_str());
sff::SFFostream<Tseries> os(ofs, opt.debug);
datrw::oanystream os(ofs, opt.outputformat);
// prepare file FREE block
sff::FREE filefree;
......@@ -739,6 +763,19 @@ int main(int iargc, char* argv[])
TFXX_debug(opt.debug, "main", "create iterators");
aff::Iterator<Tseries> S(series);
aff::Browser<Tfft::Tspectrum> C(coeff);
// take derivative if selected
if (opt.derivative)
{
for (int i=coeff.f()+1; i<=coeff.l(); ++i)
{
double frequency=df*(i-coeff.f())*2*M_PI;
double thepower=opt.nderivative;
double factor=std::pow(frequency,thepower);
coeff(i) *= factor;
}
}
TFXX_debug(opt.debug, "main", "calculate square of modulus and copy");
while(S.valid() && C.valid())
{
......@@ -797,8 +834,8 @@ int main(int iargc, char* argv[])
{
TFXX_debug(opt.debug, "main",
"scale to average in relative bandwidth");
double bwfactor=(std::pow(10.,opt.decades)-1.)/
std::pow(10.,opt.decades/2.);
double bwfactor=(std::pow(10.,opt.scaledecades)-1.)/
std::pow(10.,opt.scaledecades/2.);
// cout << "bwfactor: " << bwfactor << endl;
// cout << "index range: "
// << series.f() << "-" << series.l() << endl;
......@@ -908,10 +945,10 @@ int main(int iargc, char* argv[])
/*----------------------------------------------------------------------*/
// adjust sampling interval to frequency interval
wid2.dt=df;
if (opt.amplitudespectrum || opt.powerspectrum || opt.harmonicsignal)
{ wid2.dt=df; }
wid2.nsamples=series.size();
os << series;
os << wid2;
TFXX_debug(opt.debug, "main", " series and WID are written");
if (is.hasinfo()) { sff::INFO info; is >> info; os << info; }
......@@ -924,6 +961,7 @@ int main(int iargc, char* argv[])
tracefree.append(processfree);
os << tracefree;
}
os << series;
// write ASCII output if requested
if (opt.asciiout)
......
......@@ -44,7 +44,7 @@ edit: flist; vim $<
clean: ;
-find . -name \*.bak | xargs --no-run-if-empty /bin/rm -v
-/bin/rm -vf flist *.xxx *.xxx.* *.sff *.TZ *.TR *.grn *.su
-/bin/rm -vf *.ps *.pdf
-/bin/rm -vf *.ps *.pdf *.asc *.ascii
#----------------------------------------------------------------------
......@@ -65,6 +65,9 @@ soutifutests:
noise.xxx.ps: Makefile.foutra
/bin/rm -fv noise.xxx.*
$(MAKE) -f $< $@
noise.sig.xxx.ps: Makefile.foutra
/bin/rm -fv noise.xxx.*
$(MAKE) -f $< $@
harmspec.xxx.ps: Makefile.foutra
/bin/rm -fv sine?.xxx.sff
$(MAKE) -f $< $@
......
......@@ -42,7 +42,7 @@ NOISERMS=70.
NOISEDT=0.01
NOISET=500.
NOISEOPT=-logascii -avgascii -rbw=0.5
NOISEPROG=siggen
NOISEPROG=siggenx
NOISETYPE=14
NOISEDEBUG=
# NOISEOPT=-pad 5 -nsegments 5
......@@ -56,14 +56,16 @@ noise.xxx.fnycalc:
echo -e "0.5\n$(NOISEDT)\n/\np" > $@
noise.xxx.fny: noise.xxx.fnycalc
dc --file=$< > $@
noise.xxx.sff:
noise.xxx.ascii:
$(NOISEPROG) $(NOISETYPE) $@ $(NOISEDEBUG)\
-v -o -T $(NOISET) -d $(NOISEDT) -a $(NOISERMS)
noise.xxx.001.asc: noise.xxx.sff
-v -o -T $(NOISET) -d $(NOISEDT) -a $(NOISERMS) -ot ascii
noise.xxx.001.asc: noise.xxx.ascii
foutra -v -o -rms -ASCII=$(patsubst %.001.asc,%,$@) $(NOISEOPT) \
-power junk.xxx.sff $< \
-type ascii -power junk.xxx.sff $< \
| tee noise.xxx.log | grep rms-value > noise.xxx.rms
cat noise.xxx.log
noise.tap.xxx.ascii: noise.xxx.ascii
foutra -v -o -type ascii -Type ascii $@ $<
noise.xxx.gpt: noise.xxx.001.asc noise.xxx.fny
echo "set output 'noise.xxx.ps'" > $@
echo "set terminal postscript enhanced" >> $@
......@@ -80,6 +82,9 @@ noise.xxx.gpt: noise.xxx.001.asc noise.xxx.fny
echo "PSD t 'expected PSD'" >> $@
noise.xxx.ps: noise.xxx.gpt
gnuplot $<
noise.sig.xxx.ps: noise.xxx.ascii noise.tap.xxx.ascii
stuplox -d $@/ps -g -l 3,2,1 -h 1.2,1.2,1.5,1.2 -ty ascii \
-u counts -s x -i -Y 0.8 -c fin $^
#----------------------------------------------------------------------
# test foutra scaling for harmonic signals
......
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