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

moved scratch space to gpiag12

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: 5184
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 88a887af
this is <COPYING>
============================================================================
libgrrefsub
------
$Id$
=========================================================================
Copyright 1990 by Joachim Ungerer (files: gr_mat.f, gr_refsub.f, grrtkc.f )
Copyright 1997, 2010 by Thomas Forbriger (IfG Stuttgart, BFO Schiltach)
The source code and other files in this directory are part of libgrrefsub
which compiles to the binary library libgrrefsub.a.
It contains the reflectivity code by J. Ungerer.
libgrrefsub is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
----- END OF COPYING -----
This is a legacy version of the repository. It may be incomplete as well as
inconsistent. See README.history for details. For the old stock of the
repository copyright and licence conditions apply as specified for versions
commited after 2015-03-01. Use recent versions as a base for new development.
The legacy version is only stored to keep a record of history.
# this is <Makefile>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2008 by Thomas Forbriger (BFO Schiltach)
#
# Makefile for libgrrefsub.a
#
# reflectivity code by J. Ungerer
#
# ----
# 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
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
# ----
#
# provides an API such that reflectivity can be used within programs like
# gremlin or syg
#
# REVISIONS and CHANGES
# 15/09/2008 V1.0 Thomas Forbriger
# use gfortran
# 08/10/2010 V1.1 migrate to new SVN environment
# discard f2c fallback option
# check environment
# 17/11/2010 V1.2 avoid incdep and makefdoc.tcl
# 11/01/2011 V1.3 avoid newlib
# 17/01/2011 V1.4 distinguish library creation and installation
#
# ============================================================================
#
.PHONY: all
all: install doc
.PHONY: doc
doc: libgrrefsub.doc
LIBRARIES=libgrrefsub.a
.PHONY: install
install: $(addprefix $(LOCLIBDIR)/,$(LIBRARIES))
$(LOCLIBDIR)/%: install-include %
mkdir -pv $(LOCLIBDIR)
/bin/mv -fv $(word 2,$^) $(LOCLIBDIR)
# install-include where no header files are to be installed
.PHONY: install-include
install-include:
.PHONY: reinstall
reinstall:
$(MAKE) clean
$(MAKE) install
#----------------------------------------------------------------------
CHECKVAR=$(if $($(1)),,$(error ERROR: missing variable $(1)))
CHECKVARS=$(foreach var,$(1),$(call CHECKVAR,$(var)))
$(call CHECKVARS,LOCINCLUDEDIR LOCLIBDIR LOCBINDIR)
FLAGS += $(MYFLAGS)
CFLAGS += -O2 $(FLAGS)
FFLAGS += -ff2c -Wall -ffixed-line-length-none $(FLAGS)
# use STATIC=-static to produce statically linked binaries
STATIC=
LDFLAGS=-L$(LOCLIBDIR) $(STATIC)
CPPFLAGS=-I$(LOCINCLUDEDIR) $(FLAGS)
LIBOBS=gr_refsub.o gr_mat.o gr_rtkc.o gr_setrheo.o gr_setsigma.o
#----------------------------------------------------------------------
# standard edit targets
clean:
-/bin/rm *.o *.bak *.d *.a *.doc
flist: $(wildcard *.f *.inc) Makefile
echo $^ | tr ' ' '\n' | sort > $@
.PHONY: edit
edit: flist; vim $<
#----------------------------------------------------------------------
libgrrefsub.doc: gr_refsub.f
/usr/bin/awk 'BEGIN{ hot=0; } \
/^cS/ { hot=1; \
print " c"; \
print FILENAME; \
print " c"; \
next; } \
/^cE/ { hot=0; } \
{ if (hot==1) { print " " $$0; } }' $^ > $@
#----------------------------------------------------------------------
%.d: %.f
echo $<: $(shell cat $< | egrep '^ +include' | cut -f 2 -d \' | sort | uniq) > $@
include $(patsubst %.f,%.d,$(wildcard *.f))
%.o: %.f; $(FC) -c -o $@ $< $(FFLAGS)
libgrrefsub.a: $(LIBOBS)
ar rcv libgrrefsub.a $(LIBOBS)
ranlib libgrrefsub.a
# ----- END OF Makefile -----
c this is <gr_dim.inc>
c------------------------------------------------------------------------------
c
c 18/11/97 by Thomas Forbriger (IfG Stuttgart)
c
c all parameter definitions
c
c REVISIONS and CHANGES
c 18/11/97 V1.0 Thomas Forbriger
c
c==============================================================================
c
c mod_mlay: maximum number of layers
c
integer mod_mlay
parameter (mod_mlay=60)
c
c ----- END OF gr_dim.inc -----
c this is <gr_mat.f>
c------------------------------------------------------------------------------
c
c These are matrix inversion and matrix multiplication routines written
c by Joachim Ungerer for his reflectivity program refseis
c
c This is part of a library version of the original code 'refseis.f'
c
c Copyright 1990 by Joachim Ungerer
c
c refseis.f is published as part of the diploma thesis of Joachim Ungerer
c without a specific license
c
c----------------------------------------------------------------------
c
c Copyright 1997, 2010 by Thomas Forbriger (IfG Stuttgart)
c
c ----
c This program is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 2 of the License, or
c (at your option) any later version.
c
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c GNU General Public License for more details.
c
c You should have received a copy of the GNU General Public License
c along with this program; if not, write to the Free Software
c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
c ----
c
c REVISIONS and CHANGES
c 24/07/97 V1.0 Thomas Forbriger
c
c==============================================================================
c
C23456789012345678901234567890123456789012345678901234567890123456789012
C Subroutine fuer Matrizeninversion
C Bem: DOUBLE COMPLEX muss implementiert sein
SUBROUTINE MATINV(A11,A12,A21,A22,B11,B12,B21,B22)
COMPLEX*16 A11,A12,A21,A22,B11,B12,B21,B22,D
D=1.D0/(A11*A22-A12*A21)
B11= A22*D
B12=-A12*D
B21=-A21*D
B22=A11*D
RETURN
END
C23456789012345678901234567890123456789012345678901234567890123456789012
C Subroutine fuer Matrizenmultiplikation
C Bem: DOUBLE COMPLEX muss implementiert sein
SUBROUTINE MATMUL(A11,A12,A21,A22,B11,B12,B21,B22,C11,C12,C21,C22)
COMPLEX*16 A11,A12,A21,A22,B11,B12,B21,B22,C11,C12,C21,C22
C11=A11*B11+A12*B21
C12=A11*B12+A12*B22
C21=A21*B11+A22*B21
C22=A21*B12+A22*B22
RETURN
END
c
c ----- END OF gr_mat.f -----
c this is <gr_model.inc>
c------------------------------------------------------------------------------
c
c 24/07/97 by Thomas Forbriger (IfG Stuttgart)
c
c This file contains a common block definition to hold information on the
c layered halfspace model
c
c REVISIONS and CHANGES
c 24/07/97 V1.0 Thomas Forbriger
c
c==============================================================================
c
c get parameter declarations
include 'gr_dim.inc'
c
c mod_nlay: number of layers described by model
c mod_z: top of layer in (kilometers)
c mod_t: thickness of layer in (kilometers)
c mod_a: Vp layer velocity (kilometers/second)
c mod_b: Vs layer velocity (kilometers/second)
c mod_r: layer density (g/ccm)
c mod_Qa: Vp quality factor
c mod_Qb: Vs quality factor
c mod_aC: complex Vp layer velocity
c mod_bC: complex Vs layer velocity
c
integer mod_nlay
double precision mod_z(0:mod_mlay), mod_t(0:mod_mlay)
double precision mod_a(0:mod_mlay), mod_b(0:mod_mlay)
double precision mod_r(0:mod_mlay)
double complex mod_aC(0:mod_mlay), mod_bC(0:mod_mlay)
double precision mod_Qa(0:mod_mlay), mod_Qb(0:mod_mlay)
common /gr_model/mod_aC, mod_bC, mod_z, mod_a, mod_b, mod_r,
& mod_Qa, mod_Qb, mod_t,
& mod_nlay
c
c ----- END OF gr_model.inc -----
This diff is collapsed.
c this is <gr_results.inc>
c------------------------------------------------------------------------------
c $Id$
c
c 29/07/97 by Thomas Forbriger (IfG Stuttgart)
c
c common block to hold intermediate results for one frequency
c
c IMPORTANT!
c include gr_model.inc first!
c
c REVISIONS and CHANGES
c 29/07/97 V1.0 Thomas Forbriger
c 26/03/02 V1.1 changed to use complex frequencies
c 28/03/02 V1.2 hold scaling factor for complex frequencies
c
c==============================================================================
c
c valid frequency
double complex rtc_omega
c
c real part of frequency
double precision real_omega
c
c slowness scaling factor when using complex frequencies
double complex cplx_scal, cplx_scalq
c
c reflectivities and transmissivities of layer stacks
double complex RP11, RP12, RP21, RP22
double complex RM11, RM12, RM21, RM22
double complex TP11, TP12, TP21, TP22
c
c source amplitudes
double complex Aq0,Bq0,Cq0,Dq0
c
c surface amplitudes
double complex B0,D0
c
c phase matrix
double complex E11(0:mod_mlay), E22(0:mod_mlay)
c
c the greens coefficient for the vertical displacement
double complex GREEN_Z
c the greens coefficient for the radial displacement
double complex GREEN_R
c
c here we hold the layer matrix stack
double complex
& RTD11(0:mod_mlay),RTD12(0:mod_mlay),
& RTD21(0:mod_mlay),RTD22(0:mod_mlay),
& RBU11(0:mod_mlay),RBU12(0:mod_mlay),
& RBU21(0:mod_mlay),RBU22(0:mod_mlay),
& RTU11(0:mod_mlay),RTU12(0:mod_mlay),
& RTU21(0:mod_mlay),RTU22(0:mod_mlay),
& TTUo11(0:mod_mlay),TTUo12(0:mod_mlay),
& TTUo21(0:mod_mlay),TTUo22(0:mod_mlay),
& S11(0:mod_mlay),S12(0:mod_mlay),
& S21(0:mod_mlay),S22(0:mod_mlay)
c
c common blocks
common /gr_intermediate/ RP11, RP12, RP21, RP22,
& RM11, RM12, RM21, RM22,
& TP11, TP12, TP21, TP22,
& B0,D0,
& Aq0,Bq0,Cq0,Dq0,
& GREEN_Z, E11, E22, GREEN_R,
& cplx_scal, cplx_scalq,
& rtc_omega, real_omega
c
common /gr_stack/ RTD11, RTD12, RTD21, RTD22,
& RBU11, RBU12, RBU21, RBU22,
& RTU11, RTU12, RTU21, RTU22,
& TTUo11, TTUo12, TTUo21, TTUo22,
& S11, S12, S21, S22
c
c ----- END OF gr_results.inc -----
c this is <gr_rtc.inc>
c------------------------------------------------------------------------------
c $Id$
c
c 24/07/97 by Thomas Forbriger (IfG Stuttgart)
c
c This file defines a common block to hold precalculated reflection and
c transmission coefficients for PSV waves.
c
c IMPORTANT!
c This file includes the model definition.
c
c REVISIONS and CHANGES
c 24/07/97 V1.0 Thomas Forbriger
c 26/03/02 V1.1 changed to complex version of slowness
c 27/03/02 V1.2 introduced correct definition of integration parameter
c
c==============================================================================
c
c include the model common block to know how large the dimension should be
include 'gr_model.inc'
c
c these are the precalculated coefficients:
double complex
& Rppd(mod_mlay),Rpsd(mod_mlay),Rssd(mod_mlay),Rspd(mod_mlay),
& Tppd(mod_mlay),Tpsd(mod_mlay),Tssd(mod_mlay),Tspd(mod_mlay),
& Rppu(mod_mlay),Rpsu(mod_mlay),Rssu(mod_mlay),Rspu(mod_mlay),
& Tppu(mod_mlay),Tpsu(mod_mlay),Tssu(mod_mlay),Tspu(mod_mlay)
c
c vertical slowness components
double complex vsa(0:mod_mlay), vsb(0:mod_mlay)
c
c and here wo hold the phase slowness (seconds/kilometer) related to them:
double complex rtc_u, rtc_uq
c
c real slowness
c this is defined to be
c
c real_p = rtc_u * (1 + ime*om_sigma/real_omega)
c
c where ime is the imaginary unit.
c real_p is used as integration parameter in the slowness integral. It is the
c parameter passed to gr_prep.
c
double precision real_p
c
c imaginary part of frequency and flag which tells that we use it
double precision om_sigma
logical freq_is_complex
c
c and now the common block:
common /gr_rtcoefficients/Rppd,Rpsd,Rssd,Rspd,Tppd,Tpsd,Tssd,Tspd,
& Rppu,Rpsu,Rssu,Rspu,Tppu,Tpsu,Tssu,Tspu,
& rtc_u, rtc_uq, vsa, vsb, om_sigma, real_p, freq_is_complex
c
c ----- END OF gr_rtc.inc -----
c this is <gr_rtkc.f>
c------------------------------------------------------------------------------
c $Id$
c
c This is the Reflection- Transmission-Coefficient Subroutine
c written by Joachim Ungerer to be used in his reflectivity program refseis.
c
c This is part of a library version of the original code 'refseis.f'
c
c Copyright 1990 by Joachim Ungerer
c
c refseis.f is published as part of the diploma thesis of Joachim Ungerer
c without a specific license
c
c----------------------------------------------------------------------
c
c Copyright 1997, 2010 by Thomas Forbriger (IfG Stuttgart)
c
c ----
c This program is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 2 of the License, or
c (at your option) any later version.
c
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c GNU General Public License for more details.
c
c You should have received a copy of the GNU General Public License
c along with this program; if not, write to the Free Software
c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
c ----
c
c REVISIONS and CHANGES
c 24/07/97 V1.0 Thomas Forbriger
c 29/07/97 V1.1 just calculating PSV-case
c 26/03/02 V1.2 prepare code to be used with complex frequencies
c
c==============================================================================
c
C23456789012345678901234567890123456789012345678901234567890123456789012
C******* Joachim Ungerer ********** Stuttgart *****17.3.90*************C
C C
C Subroutine RTKC C
C C
C Berechnet die Reflektions- und Transmissionskoeffizienten einer C
C ebenen Trennflaeche zwischen zwei homogenen Halbraeumen; fuer C
C komplexe Wellengeschwindigkeiten. C
C C
C Berechnung mit doppelter Genauigkeit: C
C DOUBLE COMPLEX muss implementiert sein C
C C
C Variablendokumentation: C
C C
C Uebergabevariablen: C
C u horizontale Langsamkeit C
C a1,b1 vert.Langsamkeiten des oberen Halbraumes C
C betaC1 S-Wellengeschw. des " " (komplex) C
C rho1 Dichte des " " C
C a2,b2 vert.Langsamkeiten des unteren Halbraumes C
C betaC2 S-Wellengeschw. des " " (komplex) C
C rho2 Dichte des " " C
C C
C Rueckgabevariablen: C
C Rppd,Rpsd,Tppd,Tpsd fuer P-Welle von oben C
C Rssd,Rspd,Tssd,Tspd fuer SV-Welle " C
C rd,td fuer SH-Welle " C
C Rppu,Rpsu,Tppu,Tpsu fuer P-Welle von unten C
C Rssu,Rspu,Tssu,Tspu fuer SV-Welle " C
C ru,tu fuer SH-Welle " C
C C
C Rechenvariablen: C
C ........ gehen aus der Deklaration und C
C ihrer Berechnung hervor C
C C
C**********************************************************************C
SUBROUTINE gr_rtkc(u,a1,b1,betaC1,rho1,a2,b2,betaC2,rho2,
& Rppd,Rpsd,Tppd,Tpsd,Rssd,Rspd,Tssd,Tspd,
& Rppu,Rpsu,Tppu,Tpsu,Rssu,Rspu,Tssu,Tspu)
c
c original (including SH)
c SUBROUTINE RTKC(u,a1,b1,betaC1,rho1,a2,b2,betaC2,rho2,
c & Rppd,Rpsd,Tppd,Tpsd,Rssd,Rspd,Tssd,Tspd,rd,td,
c & Rppu,Rpsu,Tppu,Tpsu,Rssu,Rspu,Tssu,Tspu,ru,tu)
C Deklaration der Variablen ********************************************
C Variablen fuer beide Faelle
COMPLEX*16 u,uQ
REAL*8 rho1,rho2,rho12
COMPLEX*16 Rppd,Rpsd,Tppd,Tpsd,Rssd,Rspd,Tssd,Tspd,
& Rppu,Rpsu,Tppu,Tpsu,Rssu,Rspu,Tssu,Tspu,
& betaC1,betaC2,a1,a2,b1,b2,a1b1,a2b2,a1b2,a2b1,ab,
& abm,rhoabm,C,C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,
& mue1,mue2
c COMPLEX*16 Rppd,Rpsd,Tppd,Tpsd,Rssd,Rspd,Tssd,Tspd,rd,td,
c & Rppu,Rpsu,Tppu,Tpsu,Rssu,Rspu,Tssu,Tspu,ru,tu,
c & betaC1,betaC2,a1,a2,b1,b2,a1b1,a2b2,a1b2,a2b1,ab,
c & abm,rhoabm,C,C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,
c & mue1,mue2,mue1b1,mue2b2,muebN,muebP
C nur fuer Einfall von oben
COMPLEX*16 D1D,D2D,DD,D1,D2,D3
C nur fuer Einfall von unten
COMPLEX*16 D1U,D2U,DU,U1,U2,U3
C Berechnungen ********************************************************C
C Berechnungen der Rechenvariablen fuer beide Faelle
uQ = u*u
mue1 = rho1*betaC1*betaC1
mue2 = rho2*betaC2*betaC2
c mue1b1= mue1*b1
c mue2b2= mue2*b2
c muebN = mue1b1-mue2b2
c muebP = mue1b1+mue2b2
rho12 = rho1*rho2
a1b1 = a1*b1
a2b2 = a2*b2
a1b2 = a1*b2
a2b1 = a2*b1
ab = a1b1*a2b2
abm = a1b2-a2b1
rhoabm= 2.D0*rho12*abm
C = 2.D0*(mue1-mue2)
C0 = C*uQ
C1 = C0-rho1
C2 = C0+rho2
C3 = C1+rho2
C4 = C2*a1-C1*a2
C5 = C2*b1-C1*b2
C6 = C3*C3*uQ
C7 = C*C0*ab
C8 = C1*C1*a2b2+rho12*a2b1
C9 = C2*C2*a1b1+rho12*a1b2
C10 = C3+C*a2b1
C11 = C3+C*a1b2
C Koeffizienten fuer Einfall von oben
C Rechenvariablen
D1D = C6+C8
D2D = C7+C9
DD = D1D+D2D
D1 = 2.D0*a1/DD
D2 = 2.D0*b1/DD
D3 = C3*C2+C*C1*a2b2
C Koeffizienten
Rppd = (D2D-D1D)/DD
Rpsd =-u*D1*D3
Tppd = rho1*D1*C5
Tpsd =-u*rho1*D1*C10
Rssd = Rppd-rhoabm/DD