Commit e1ad060d authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

first compiling version

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.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 1024
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent ea178b91
# this is <Makefile>
# ----------------------------------------------------------------------------
# $Id: Makefile,v 1.1 2002-11-08 20:33:35 forbrig Exp $
# $Id: Makefile,v 1.2 2002-11-11 15:27:00 forbrig Exp $
#
# Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
#
......@@ -16,16 +16,13 @@
all:
flist: Makefile
flist: Makefile $(wildcard *.f *.inc)
echo $^ | tr ' ' '\n' | sort > $@
.PHONY: edit
edit: flist; vim $<
.PHONY: clean
clean: ;
-find . -name \*.bak | xargs --no-run-if-empty /bin/rm -v
-/bin/rm -vf flist
ifdef CROSS_BASE
BINPREFIX=$(CROSS_BASE)/bin/dos-
......@@ -44,8 +41,9 @@ F2CFLAGS=-f -u
#CC=gcc
CFLAGS=-O2 -I${SERVERINCLUDEDIR} -I${LOCINCLUDEDIR}
LIBSRC=
LIBSRC=$(wildcard *.f)
LIBOBS=$(patsubst %.f,%.o,$(LIBSRC))
DOCS=libfourier.doc
docs: $(DOCS)
......@@ -58,9 +56,9 @@ docs: $(DOCS)
$(FC) $(FFLAGS) -c -o $@ $<
clean:
-/bin/rm *.o *.bak *.doc *.o77
-/bin/rm -fv *.o *.bak *.doc *.o77 flist
libfourier.doc: $(LIBSRC)
libfourier.doc: $(LIBSRC) $(wildcard *.inc)
makefdoc.pl $@ $^
......
c this is <filters.f>
c ----------------------------------------------------------------------------
c ($Id: filters.f,v 1.1 2002-11-11 15:27:00 forbrig Exp $)
c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c
c ready made filters
c
c REVISIONS and CHANGES
c 11/11/2002 V1.0 Thomas Forbriger
c
cS
c ============================================================================
c
subroutine foufil_clear
c
c initialize fourier filters
c
include 'filters.inc'
c
cE
c
call fou_clear
call foufil_frequency
c
return
end
c
cS
c----------------------------------------------------------------------
c
double precision function foufil_omega(in)
c
c calculate angular frequency
c in: will be treated as a frequency given in Hz, if fourier_frequency is true
c otherwise it will be treated as a period given in sec
c
include 'filters.inc'
c
double precision in
c
cE
c
double precision pi2, value
parameter(pi2=2.d0*3.14159265358979d0)
c
if (fourier_frequency) then
value=pi2*in
else
value=pi2/in
endif
c
foufil_omega=value
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine foufil_frequency
c
c set frequency mode
c
include 'filters.inc'
c
cE
fourier_frequency=.true.
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine foufil_period
c
c set period mode
c
include 'filters.inc'
c
cE
fourier_frequency=.false.
c
return
end
c
c======================================================================
c define filters
c
cS
c----------------------------------------------------------------------
c
subroutine foufil_int
c
c integrate once
c
cE
double complex value
value=(0.d0,0.d0)
call fou_pole(value)
value=(0.d0,1.d0)
call fou_cdenom(value)
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine foufil_dif
c
c differentiate once
c
cE
double complex value
value=(0.d0,0.d0)
call fou_zero(value)
value=(0.d0,1.d0)
call fou_cnumer(value)
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine foufil_hpb(in, ord)
c
c Butterworth high pass filter
c in: frequency or period (cf. foufil_omega)
c ord: order of filter
c
double precision in
integer ord
c
cE
double precision foufil_omega, omega, pi
integer i
double complex ime, pole, zero
parameter (ime=(0.d0,1.d0), pi=3.1415926535897931159d0)
c
if (ord.lt.1) stop 'ERROR (foufil_hpb): illegal order'
omega=foufil_omega(in)
c
zero=(0.d0,0.d0)
do i=1,ord
call fou_zero(zero)
pole=omega*exp(ime*pi*(2.d0*i-1.d0)/(2.d0*ord))
call fou_pole(pole)
enddo
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine foufil_lpb(in, ord)
c
c Butterworth low pass filter
c in: frequency or period (cf. foufil_omega)
c ord: order of filter
c
double precision in
integer ord
c
cE
double precision foufil_omega, omega, pi
integer i
double complex ime, pole, cfac
parameter (ime=(0.d0,1.d0), pi=3.1415926535897931159d0)
c
if (ord.lt.1) stop 'ERROR (foufil_hpb): illegal order'
omega=foufil_omega(in)
c
cfac=-ime*omega
do i=1,ord
call fou_cnumer(cfac)
pole=omega*exp(ime*pi*(2.d0*i-1.d0)/(2.d0*ord))
call fou_pole(pole)
enddo
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine foufil_lp2(in, h)
c
c second order low pass filter
c in: frequency or period (cf. foufil_omega)
c h: damping as a fraction of critical
c
double precision in, h
c
cE
double precision foufil_omega, omega
double complex ime, pole, cfac
parameter (ime=(0.d0,1.d0))
c
omega=foufil_omega(in)
c
cfac=-ime*omega
call fou_cnumer(cfac)
call fou_cnumer(cfac)
pole=omega*(ime*h+sqrt(1.d0-h*h))
call fou_pole(pole)
pole=omega*(ime*h-sqrt(1.d0-h*h))
call fou_pole(pole)
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine foufil_hp2(in, h)
c
c second order high pass filter
c in: frequency or period (cf. foufil_omega)
c h: damping as a fraction of critical
c
double precision in, h
c
cE
double precision foufil_omega, omega
double complex ime, pole, zero
parameter (ime=(0.d0,1.d0))
c
omega=foufil_omega(in)
c
zero=(0.d0,0.d0)
call fou_zero(zero)
call fou_zero(zero)
pole=omega*(ime*h+sqrt(1.d0-h*h))
call fou_pole(pole)
pole=omega*(ime*h-sqrt(1.d0-h*h))
call fou_pole(pole)
c
return
end
c
c
c ----- END OF filters.f -----
c this is <filters.inc>
c ----------------------------------------------------------------------------
c ($Id: filters.inc,v 1.1 2002-11-11 15:27:00 forbrig Exp $)
c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c
c some filter settings
c
c REVISIONS and CHANGES
c 11/11/2002 V1.0 Thomas Forbriger
c
cS
c ============================================================================
c
c define coefficients by frequency or by period
logical fourier_frequency
c
c common block
common /fourier_fil/ fourier_frequency
c
cE
c ----- END OF filters.inc -----
c this is <polesnzeroes.f>
c ----------------------------------------------------------------------------
c ($Id: polesnzeroes.f,v 1.1 2002-11-11 15:27:00 forbrig Exp $)
c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c
c pole and zero access subroutines
c
c REVISIONS and CHANGES
c 11/11/2002 V1.0 Thomas Forbriger
c
c ============================================================================
cS
c
subroutine fou_clear
c
c initialize poles and zeros database
c
include 'polesnzeros.inc'
c
cE
call fou_normal
fourier_numer=(1.d0,0.d0)
fourier_denom=(1.d0,0.d0)
fourier_npoles=0
fourier_nzeros=0
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine fou_inverse
c
c set inverse flag
c all following coefficients will be treated as inverse filters
c
include 'polesnzeros.inc'
c
cE
c
fourier_normal=.false.
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine fou_normal
c
c reset inverse flag
c all following coefficients will be treated as normal filters
c
include 'polesnzeros.inc'
c
cE
c
fourier_normal=.true.
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine fou_pole(pole)
c
c set one pole
c
double complex pole
c
include 'polesnzeros.inc'
c
cE
c
if (fourier_normal) then
if (fourier_npoles.lt.fourier_nmax) then
fourier_npoles=fourier_npoles+1
fourier_poles(fourier_npoles)=pole
else
stop 'ERROR (fourier_pole): too many poles'
endif
else
if (fourier_nzeros.lt.fourier_nmax) then
fourier_nzeros=fourier_nzeros+1
fourier_zeros(fourier_nzeros)=pole
else
stop 'ERROR (fourier_pole): too many zeros'
endif
endif
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine fou_zero(zero)
c
c set one zero
c
double complex zero
c
include 'polesnzeros.inc'
c
cE
c
if (fourier_normal) then
if (fourier_nzeros.lt.fourier_nmax) then
fourier_nzeros=fourier_nzeros+1
fourier_zeros(fourier_nzeros)=zero
else
stop 'ERROR (fourier_zero): too many zeros'
endif
else
if (fourier_npoles.lt.fourier_nmax) then
fourier_npoles=fourier_npoles+1
fourier_poles(fourier_npoles)=zero
else
stop 'ERROR (fourier_zero): too many poles'
endif
endif
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine fou_numer(factor)
c
c set one numerator factor
c
double precision factor
c
include 'polesnzeros.inc'
c
cE
c
if (fourier_normal) then
fourier_numer=factor*fourier_numer
else
fourier_denom=factor*fourier_denom
endif
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine fou_denom(factor)
c
c set one denominator factor
c
double precision factor
c
include 'polesnzeros.inc'
c
cE
c
if (fourier_normal) then
fourier_denom=factor*fourier_denom
else
fourier_numer=factor*fourier_numer
endif
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine fou_cnumer(factor)
c
c set one complex numerator factor
c
double complex factor
c
include 'polesnzeros.inc'
c
cE
c
if (fourier_normal) then
fourier_numer=factor*fourier_numer
else
fourier_denom=factor*fourier_denom
endif
c
return
end
c
cS
c----------------------------------------------------------------------
c
subroutine fou_cdenom(factor)
c
c set one complex denominator factor
c
double complex factor
c
include 'polesnzeros.inc'
c
cE
c
if (fourier_normal) then
fourier_denom=factor*fourier_denom
else
fourier_numer=factor*fourier_numer
endif
c
return
end
c
cS
c----------------------------------------------------------------------
c
double complex function fou_eval(omega)
c
c evaluate filter response coefficient at real frequency omega
c
double precision omega
c
include 'polesnzeros.inc'
c
cE
double complex value, thenumer, thedenom
integer i
c
thenumer=(1.d0,0.d0)
do i=1,fourier_nzeros
thenumer=thenumer*(fourier_zeros(i)-omega)
enddo
c
thedenom=(1.d0,0.d0)
do i=1,fourier_npoles
thedenom=thedenom*(fourier_poles(i)-omega)
enddo
c
value=thenumer*fourier_numer/(fourier_denom*thedenom)
c
fou_eval=value
c
return
end
c
c ----- END OF polesnzeroes.f -----
c this is <polesnzeros.inc>
cS
c ----------------------------------------------------------------------------
c ($Id: polesnzeros.inc,v 1.1 2002-11-11 15:27:00 forbrig Exp $)
c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c
c common containing pole and zero collection
c
c REVISIONS and CHANGES
c 11/11/2002 V1.0 Thomas Forbriger
c
c ============================================================================
c
c dimension
integer fourier_nmax
parameter(fourier_nmax=20)
c
c numerator
double complex fourier_numer
c
c denominator
double complex fourier_denom
c
c complex poles
double complex fourier_poles(fourier_nmax)
c
c complex zeros
double complex fourier_zeros(fourier_nmax)
c
c fill mode (normal or inverse)
logical fourier_normal
c
c counters
integer fourier_npoles, fourier_nzeros
c
c common block
common /fourier_pnz/ fourier_numer, fourier_denom,
& fourier_poles, fourier_zeros,
& fourier_npoles, fourier_nzeros, fourier_normal
c
cE
c ----- END OF polesnzeros.inc -----