c ********************************************************************
C                                                                      C
C Programm REFSEIS (Version 02.10.1990 by Joachim Ungerer)             C
C Developed for SUN 4/110 workstation, written in SUN-FORTRAN          C
C                                                                      C 
C Part of the thesis of Joachim Ungerer                                C
C Work performed at Institute of Geophysics, Stuttgart University      C
C                                                                      C
C Functions and limitations of REFSEIS.F:                              C
C The program REFSEIS.F generates complete (body and surface wave)     C
C synthetic near- and far-field seismograms for a generalized point    C
C source of Brustle-Muller model or an impulse excitation which is     C
C located within horizontal layered medium. Theoretical basis of the   C
C program is presented in the thesis [Ungerer, 1990].                  C
C Medium can be arbitrarily layered. As a special case, full space can C
C be modeled. Layer parameters include velocities of P and S waves,    C
C density, Q of P and S waves. In the framework of acausal damping     C
C model, strong damping effect can be satisfactorily taken into        C
C account which depends on the station distribution and layer          C
C parameters. Elastic case can be approximated by very high Q.         C
C REFSEIS.F calculates near- and far-field seismograms separately      C
C [Ungerer, 1990]. Either displacement (Typ=1) or velocity (Typ=2) or  C
C acceleration (Typ=3) is calculated. The source excitation function   C
C is implemented in the first part of the program. In addition to the  C
C seismograms, their corresponding spectra can also be obtained in the C
C output.                                                              C
C Seismograms are represented in a cylindrical coordinate system with  C
C origin at the epicenter as three components -- radial(r),            C
C azimuthal(phi) and vertical(z).                                      C
C The number of receivers can be arbitrary in principle. their         C
C locations with respect to the source is described by epicentral      C
C distance and azimuth.                                                C
C There is only one limitation on the location of the source. It can   C
C not be put at the surface, since the analytic expression on which    C
C the program is based would no longer be defined [Ungerer, 1990].     C
c ***** The source must not be vertically under a station *****        c
C The source mechanism and its orientation are described by the        C
C Cartician components of the source moment tensor.                    C
C Difficulties and numerical disturbances arise mainly in the slowness C
C integration. The integral is exact only when the upper bound is      C
C unlimited. Numerically we need to truncate it. Either a small        C
C integration range or a rough sampling can lead to considerable       C
C disturbances or the truncation phase in the resultant seismograms    C
C which are dependent on the source depth, the epicentral distance of  C
C the receiver and the Q factor. These disturbances can never be       C
C completely eliminated. But with proper selection of integration      C
C range, sampling interval, source depth and Q, useful seismograms can C
C be produced by REFSEIS which allows the investication of the near-   C
C and far-field seismograms from a particular event.                   C
C                                                                      C
C *********************************************************************C
c
C Program by Joachim Unterer, 1990
c Corrections by Karl Koch, 1991, 1995
c Diese Version enthaelt die Besselfunktions-Subroutinen nach
c "Numerical Recipes in Fortran", nicht optimal da nicht 
c echt doppeltgenau, E.W. Maerz 2004
c Kleinere Anpassungen an den NAS FortranPlus Compiler, E.W. Maerz 04
c minor updates October 2009, E.W.
c
c *******************************************************************
c
c DEFINITIONS (translation of the German text below, by EW)
c The author's definitions are not quite clear. According to a few
c test cases, it appears that the following definitions are used:
c The moment tensor is defined as in Aki and Richard's book (first
c edition, p. 117, Box 4.4):
c strike: cw from North, fault plane dipping to the right
c dip: measured against the horizontal
c rake: direction of displacement of the upper block, measured in the
c fault plane from the horizontal upward (ccw as seen from above)
c
c Azimuth: from epicenter to station, cw from North
c Radial component: outward from the epicenter
c Transverse component: cw from radial
c Vertical: up
c
c *******************************************************************
c
c Bem. von EW, August 2004:
c Ungerers Definitionen des Momenten-Tensors, des Azimuts und der
c Seismogrammkomponenten koennten missverstanden werden.
c
c Anhand einiger Testfaelle halte ich folgende Def. fuer zutreffend:
c => Momenten-Tensor wie in Aki & Richards (1. Auflage) S. 117 Box 4.4
c    Die Formeln von A&R sind im Matlab-Skript "Momentens" programmiert.
c    Dieses Skript erzeugt auch ein Abstrahldiagramm fuer P-Wellen.
c => Stations-Azimut vom Epizentrum zur Station, von Nord ueber Ost (cw)
c => Radialkomponente weg vom Epizentrum
c => Transversalkomponente rechts (90 Grad cw) von der radialen
c => Vertikalkomponente nach oben (also alles wie ueblich)
c
c Definitionen der Herdwinkel in A&R und im Matlab-Skript:
c => strike analog zum Azimut, Orientierung so dass Flaeche nach
c           rechts einfaellt
c => dip    gegen die Horizontale, <= 90 Grad
c => rake   Verschiebungsrichtung des Hangenden (also des Blocks
c           rechts von der Bruchflaeche). Der Winkel wird in der 
c           Bruchflaeche von der Horizontalen ausgehend ccw gemessen
c           (also nach oben positiv)
c           0 Grad: Horizontalverschiebung, 90 Grad: Aufschiebung
c
c *********************************************************************c
C                                                                      C
C Brief description of the program:                                    C
C REFSEIS.F calculates the seismogram components of a generalized      C
C point source with strength M0 (seismic moment) at depth ZQ within a  C
C horizontally layered medium according to reflectivity method.        C
C Type and orientation of the source are defined by the components of  C
C its moment tensor. Excitation of the source is Brustle-Muller model  C
C or a delta function.                                                 C
C Seismogram components are calculated separately as:                  C
C      Far-field components   SFr,SFphi,SFz                            C
C      Near-field components  SNr,SNphi,SNz                            C
C They are represented in a cylindrical coordinate system with the     C
C origin at the epicenter.                                             C
C The seismogram can be reduced to receiver dependent first arrival    C
C time by reduction velocity Vred. They can also be shifted recerver-  C
C independently to the left by Tli amount or to the right by Tre       C
C amount.                                                              C
C Different seismogram types can be generated:                         C
C      -Displacement seimogram   (Typ=1)                               C
C      -Velocity seismogram   (Typ=2)                                  C
C      -Acceleration seismogram    (Typ=3)                             C
C                                                                      C
C NOTE!!!! IF THE SOURCE IS AN IMPULSE (Thd = 0), ONLY DISPLACEMENT IS C
C          GENERATED NO MATTER WHAT TYPE IS SET !!!!!!!!!!!!!!!!!!!    C
C                                                                      C
C Input is read from the file:   refseis.par                           C
C Protocol is written to the file:  refseis.out                        C
C Spectra are written to refspec.plt                                   C
C Seismograms are written to refseis.sig                               C
C                                                                      C
C The files refseis.spe and refseis.sig are optionally generated.      C
C                                                                      C
C The file refseis.out is a protocol. The kind of                      C
C information contained is determined through the switches.            C
C  (1) vertical slownesses of layers and reflection and transmission   C
C      coefficients of the interfaces.                                 C
C  (2) reflectivities and transmissivities of the source layer         C
C          Reflectivity matrix      RM(inus) = (RM11,RM12,RM21,RM22)   C
C                    "              RP(lus)  = (RP11,RP12,RP21,RP22)   C
C          Reflectivity scalar      rm(inus) = rm                      C
C                    "              rp(lus)  = rp                      C
C          Transmissivity matrix    TP(lus)  = (TP11,TP12,TP21,TP22)   C
C          Transmissivity scalar    tp(lus)  = tp                      C
C  (3) surface amplitutes         B0(0,1,2),D0(0,1,2),F0(1,2)          C
C  (4) normalized impules spectra (complex) VF,VN in [s*s/(g*cm)]*1E-20C
C  (5) spectrum moduli    BSF,BSN of source strength   M0              C
C  (7) far- and near-field seismograms  in [cm,cm/s,cm/(s*s)]          C
C                                                                      C
C                                                                      C
C subroutines:                                                         C
C       RTKC      calculates reflection and transmission coefficients  C
C       DFORK     fast Fourier transform                               C
C       MATINV    Matrix invertion                                     C
C       MAMULT    Matrix multiplication                                C
C All subroutines are at the end of the program.                       C
C                                                                      C
C Deviations from the FORTRAN standard:                                C
C special functions:                                                   C
C       d_j0,d_j1,d_jn  Bessel functions of 0-th, 1-th and n-th order  C
C                       (when called in SUN FORTRAN, should be         C
C                        declared as DOUBLE PRECISION)                 C
C       dtime(tarray)   determination of calculation time              C
C                                                                      C
C extended command and functions:                                      C
C      Die Berechnungen finden mit doppelter Genauigkeit statt.        C
C      Daher laeuft dieses Programm nur, wenn in FORTRAN folgende      C
C      Befehle implementiert sind:  DOUBLE COMPLEX bzw COMPLEX*16,     C
C                                   CDABS,CDEXP,CDSQRT,DCMPLX,DCONJG,  C
C                                   DREAL                              C
C      fuer NAS-Fortran statt cdabs usw.: abs, exp, sqrt               C
C                                                                      C
C Note1: As the Fourier transform of the displacement has a pole at    C
C        frequency zero, we define fr(1)=1.D-30 for displacement       C
C        calculation.                                                  C
C                                                                      C
C Note2: The occurrence of the truncation phase may be reduced by      C
C        the extention of the slowness integration range or the        C
C        increase of the sampling points Nu.                           C
C                                                                      C
C                                                                      C
C Description of variables:                                            C
C (The names of the variables are chosen mostly in accordance to the   C
C thesis)                                                              C
C                                                                      C
C    Eingabevariablen:                                                 C
C          Die Eingabe des Modells sollte nach der Vorlage des Ein-    C
C          gabefiles refseis.par erfolgen.                             C
C          Bis auf die Schaltervariable werden alle Eingabevariablen   C
C          formatiert eingelesen, sie sind deshalb rechtsbuendig in    C
C          ihre jeweiligen Eingabefelder zu setzen.                    C
C          text0         Kopfzeile fuer Text                           C
C      I.  Schalter (0=nein, 1=ja, 2=ja und Ende)                      C
C          Sch(1)        vert.Langsamkeiten,Ref.u.Transmissionskoeff.  C
C          Sch(2)        Reflektivitaeten und Transmissivitaeten       C
C          Sch(3)        Oberflaechenamplituden                        C
C          Sch(4)        normierte Impuls-Spektren(komplex)            C
C          Sch(5)        endgueltige Spektren(Betraege)                C
C          Sch(6)        Plot-File fuer die endgueltigen Spektren      C
C          Sch(7)        Endgueltige Nah- und Fernfeldseismogramme     C
C          Sch(8)        Plot-File fuer Seismogramme                   C
C          Sch(9)        noch frei                                     C
C          Sch(10)       noch frei                                     C
C      II. Schichtdaten                                                C
C          z(i)          Tiefe der Schichtdecke der i-ten Schicht [km] C
C                        (die "nullte" Schicht repraesentiert den obe- C
C                        ren Halbraum,z.B. Atmosphaere,ihre Tiefe ist  C
C                        definitionsgemaess Null)                      C
C          alpha(i)      P-Wellengeschw.[km/s] in Schicht i            C
C          beta(i)       S-   "                "                       C
C          rho(i)        Dichte[g/ccm] der Schicht i                   C
C          Qa(i)         Q-Faktor fuer P-Wellen                        C
C          Qb(i)            "       "  S-Wellen                        C
C          Bem: Grosse Schichtdicken koennen zum Ueberlauf durch Ex-   C
C          ponentialterme fuehren. Abhilfe durch Einfuegen virtueller  C
C          Trennflaechen.                                              C
C      III.Herddaten                                                   C
C          ZQ            Tiefe[km] der ErdbebenQuelle                  C
C          M0            Staerke des Momententensors in [dyn*cm]       C
C          Mxx,Myy,Mzz   normierter Momententensor  M=(ni*fj+nj*fi)    C
C          Mxy,Mxz,Myz   (Bem:seismischer Momententensor MS=M0*M)      C
C          Bem1: Ein grosser Abstand des Herdes zur Schichtdecke kann  C
C          zum Ueberlauf durch die Exponentialterme ea,eb fuehren. Ab- C
C          hilfe durch Einfuegen einer virtuellen Trennflaeche nahe    C
C          ueber dem Herd.                                             C
C          Bem2: Die Umrechnung von Streichen(S),Fallen(F) und Dislo-  C
C          kationswinkel(D) beim Double-Couple in die Vektoren f und n C
C          erfolgt durch: f=(cosD*cosS+cosF*sinD*sinS,cosD*sinS-cosF*  C
C          sinD*cosS,-sinD*sinF); n=(sinF*sinS,-sinF*cosS,cosF)        C
C          [siehe hierzu auch Aki/Richards S.106 ff]                   C
C      IV. Daten fuer Normierte Anregungsfunktion                      C
C          typ      1=Verschiebung,2=Geschwindigkeit,3=Beschleunigung  C
C          The      Einsatzzeit des Herdes                             C
C          Thd      Dauer des Zeitverlaufs von The bis Endwert Null    C
C                   (Momentenanstiegszeit)                             C
C          Bem: Impulsanregung fuer The=Thd=0                          C
C      V.  Reduktionsdaten                                             C
C          Vred     Reduktionsgeschwindigkeit[km/s] (bei Vred=0 keine  C
C                   empfaengerabhaengige Reduktion),die zugehoerige    C
C                   Reduktionszeit ist hypozentralbezogen.             C
C          Tli      Zeitverschiebung[s] des Seismogramms nach rechts   C
C          Tre      Zeitverschiebung[s] des Seismogramms nach links    C
C      VI. Empfaengerdaten                                             C
C          r(E)      Radialentf.[km] des Empfaengers E vom Epizentrum  C
C          phi(E)    Winkel[Grad] des Empfaengers E zum Epizentrum     C
C     VII. Langsamkeitsdaten                                           C
C          umin,umax  Langsamkeitsbereich [s/km]                       C
C          uwil,uwir  Langsamkeits-Window links,rechts [s/km]          C
C                     (Ueblich sind umin=uwil=0, uwir sollte reziprok  C
C                     der kleinsten Wellengeschwindigkeit der Herd-    C
C                     schicht oder groesser sein)                      C
C          Nu         Zahl der gleichabstaendigen Langsamkeit.(NU muss C
C                     mindestens 2 sein)                               C
C          Bem:Langsamkeitsbereich und Integrationsgenauigkeit Nu sind C
C          die masgebenden Groessen um moeglichst stoerungsfreie Seis- C
C          mogramme zu erzeugen. Das rechteckigen Langsamkeitswindow   C
C          [uwil,uwir] wird beidseitig mit einem Cos-Taper bis zu den  C
C          jeweiligen Grenzen [umin,umax] auf Null heruntergezogen ->  C
C          evt. Trennung von Raum und Oberflaechenwellen durch Wahl    C
C          geeigneter Langsamkeitsbereiche.                            C
C     VIII.Frequenz- und Zeitdaten                                     C
C          fmin,fmax     Frequenzbereich [Hz](physikalisch)            C
C          fwil,fwir     Frequenz-Window links,rechts [Hz]             C 
C          Dt            zeitliches Abtastintervall [sec]              C
C          SL            Seismogrammlaenge (Zweierpotenz)=Anzahl der   C
C                        Stuetzstellen des Seismogrammes(maximal=MSL)  C
C          Bem: das rechteckige Frequenzwindow [fwil,fwir] wird beid-  C
C          seitig mit einem Cos-Taper bis zu den Grenzen [fmin,fmax]   C
C          auf Null heruntergezogen -> Frequenzfilter                  C
C                                                                      C
C    Laufvariablen                                                     C
C          E       fuer Empfaenger                                     C
C          f       fuer Frequenzen, teilweise auch fuer Zeiten         C
C          i       fuer Schichten, und als Zaehlindex                  C
C          k       Zaehlindex                                          C
C          l       fuer Langsamkeiten                                  C
C                                                                      C
C    Parameter(koennen nach Belieben im Programm eingestellt werden)   C
C          IME     Imaginaere Einheit                                  C
C          ME      Max.Zahl der Empfaenger                             C
C          Mf          "        physikalischen Frequenzen(=MSL/2+1)    C
C                               bis zur Nyquistfrequenz Fny            C
C          MS          "        Schichten                              C
C          MSL         "        Stuetzstellen (Max.Seismogrammlaenge)  C
C          hin,rueck  -1.,+1. zur Steuerung der Fouriertransformation  C
C          pi      =3.14159265358979                                   C
C                                                                      C
C    Charaktervariablen                                                C
C          typ0    ='Verschiebung','Geschwindigkeit' oder              C
C                   'Beschleunigung', je nach Schalter(9)              C
C          impuls0 ='Impulsanregung'                                   C
C                                                                      C
C    Variablen fuer Rechenzeit                                         C
C          dtime    Funktion zur Bestimmung der Rechenzeit             C
C          tarray   ihre zweidim. Variable(User-Zeit,System-Zeit)      C
C          time     Aufnahmevariable fuer Rechenzeit                   C
C                                                                      C
C    Hilfsvariablen zur Zwischenspeicherung                            C
C          help[reell],hilf,gew                                        C
C          H11,H12,H21,H22                                             C
C          I11,I12,I21,I22                                             C
C                                                                      C
C    Rechenvariablen(fest)                                             C
C          ap   =Wurzel(SL)*Dt Aperiodizitaetsfaktor fuer Fourierkoeff.C
C          alphC(i)  komplexe P-Wellengeschwindigkeit                  C
C          betaC(i)     "     S-Wellen     "                           C
C          d(i)      Dicke der Schicht i in [km]                       C
C          Df        Schrittweite der Frequenz,Frequenzintervall       C
C          Du        Schrittweite der Langsamkeit u                    C
C          DDu       getaperte und gewichtete Integrationsschrittweite C
C          FL        Frequenzlaenge des vollstaendigen Spektrums       C
C          fmi,fma   Index der Frequenzen fmin,fmax                    C
C          Fny       Nyquistfrequenz (fmax muss < Fny sein)            C
C          fr(f)     f-te Frequenz in [Hz]                             C
C          fwl,fwr   Index der Frequenzen fwil,fwir                    C
C          h         Schicht in der der Herd sich befindet             C
C          M0dim     =M0/1.D20 zur Dimensionierung der Seismogramme    C
C          n         Zahl der Schichten                                C
C          NE        Zahl der Empfaenger                               C
C          Nf        Zahl der phys. Frequenzen im Bereich [fmax,fmin]  C
C          Nff        "    "    "       "      im Fenster [fwil,fwir]  C
C          Ny        Index der Nyquistfrequenz Fny (max.phys.Freqenz)  C
C          phiB(E)   Winkel des Empfaengers E in Bogenmass             C
C          t(f)      f-te Zeit in [s], d.h. Zeitwerte des Seismogramms C
C          TL        Zeitlaenge des Seismogramms (Periode)             C
C          tred(f)   reduzierte Zeitskala  / modified 8-2-91 K.Koch /  C
C                                                                      C
C          Variablen der normierten Herdfunktion                       C
C          hfkt(tvar,T1,T2) Definiton der Herdfunktion                 C
C          tvar           ihre Zeitvariable                            C
C          T1,T2          ihre Begrenzungsparameter(T1=The,T2=The+Thd) C
C          g(f)           Feld in dem hfkt je Frequenz abgelegt wird   C
C          Die Dimension der Herdfunktion richtet sich nach ihrer In-  C
C          terpretation:als Verschiebung [ ]  ->  Bodenverschiebung    C
C                       als 1.Ableitung [1/s] ->  Bodengeschwindigkeit C
C                       als 2.Ableitung [1/(s*s)]->  "  beschleunigung C
C                                                                      C
C    Rechenvariablen(veraenderlich)                                    C
C          a(i)     vert.Langsamkeit der P-Wellen in Schicht i [s/km]  C
C          b(i)      "        "      der S-Wellen       "              C
C          E11(i),E22(i)     Phasenmatrix der Schicht i (Diagonalel.)  C
C          ea         =exp[IME*w*a(h)*(ZQ-z(h))]                       C
C          eb         =   "      b(h)    "                             C
C          earg       Argument fuer Exponentialfunktionen              C
C          ftap(f)    Fensterfunktion der Frequenz mit Cos-Taper zur   C
C                     Verwendung als Frequenzfilter,das Rechteckwindow C
C                     [fwil,fwir] wird beidseitig mit einem Cos-Taper  C
C                     bis zu den Frequenzen fmax bzw. fmin auf Null    C
C                     heruntergezogen.                                 C
C          red        Reduktionsfaktor:red=exp[i*w*(R(E)/Vred+Tli-Tre] C
C                     mit R(E)=Wurzel[ZQ*ZQ+r(E)*r(E)]=Hypozentralentf.C
C          u          horizontale Langsamkeit in [s/km]                C
C          uQ         ihr Quadrat                                      C
C          utap       Fensterfunktion der Langsamkeit mit Cos-Taper    C
C                     zum Glaetten der Stoerphasen                     C
C          w(f)       Kreisfrequenz Omega [1/s]                        C
C                                                                      C
C          Reflektions und Transmissionskoeffizienten d.Trennflaeche i C
C          Rppd(i),Rpsd(i),Tppd(i),Tpsd(i)  fuer  P-Wellen von oben    C
C          Rssd(i),Rspd(i),Tssd(i),Tspd(i)  fuer SV-Wellen      "      C
C          rd(i),td(i)                      fuer SH-Wellen      "      C
C          Rppu(i),Rpsu(i),Tppu(i),Tpsu(i)  fuer  P-Wellen von unten   C
C          Rssu(i),Rspu(i),Tssu(i),Tspu(i)  fuer SV-Wellen      "      C
C          ru(i),tu(i)                      fuer SH-Wellen      "      C
C                                                                      C
C          Reflektivitaets- und Transmissivitaetsmatrizen              C
C          RTD11(i),RTD12(i) Reflektivitaetsmatrix f.Einfall von oben, C
C        & RTD21(i),RTD22(i) fuer unteren Stapel ab Decke Schicht i    C 
C          RBU11(i),RBU12(i) Reflektivitaetsmatrix f.Einfall v. unten, C
C        & RBU21(i),RBU21(i) fuer oberen Stapel ab Boden Schicht i     C
C          RTU11(i),RTU12(i) Reflektivitaetsmatrix f.Einfall v. unten, C
C        & RTU21(i),RTU22(i) fuer oberen Stapel ab Decke Schicht i     C
C          TTUo11(i),TTUo12(i) Transmissivitaetsmatrix f.Einf.v.unten, C
C        & TTUo21(i),TTUo22(i) ab Decke Schicht i bis Decke Schicht 0  C
C          S11(i),S12(i)     Schichtmatrix von Schicht i,zur Berech-   C
C        & S21(i),S22(i)     nung der Transmissivitaetsmatrix TTUo     C
C                                                                      C
C          Reflektivitaets- und Transmissivitaetsskalare               C
C          rtd(i)   SH-Reflektivitaet ab Decke Schicht i abwaerts      C
C          rbu(i)   SH-Reflektivitaet ab Boden Schicht i aufwaerts     C
C          rtu(i)   SH-Reflektivitaet ab Decke Schicht i    "          C
C          ttuo(i)  SH-Transmissivitaet ab Decke Schicht i  "          C
C          s(i)     Schichtskalar von Schicht i,fuer Skalar TTU        C  
C                                                                      C
C          Reflektivitaeten und Transmissivitaeten der Herdschicht     C
C          RM11,RM12,RM21,RM22   Reflektivitaetsmatrix RM(inus)        C
C          RP11,RP12,RP21,RP22   Reflektivitaetsmatrix RP(lus)         C
C          TP11,TP12,TP21,TP22   Transmissivitaetsmatrix TP(lus)       C
C          rm                    SH-Reflektivitaet rm(inus)            C
C          rp                    SH-Reflektivitaet rp(lus)             C
C          tp                    SH-Transmissivitaet tp(lus)           C
C                                                                      C
C          Amplituden der Wellenfelder                                 C
C          AQ(0,1,2)  Quell-Amplituden des abwaertslaufenden           C
C          CQ(0,1,2)  Wellenfeldes der Quelle                          C
C          EQ(1,2)         "                                           C
C          BQ(0,1,2)  Quell-Amplituden des aufwaertslaufenden          C
C          DQ(0,1,2)  Wellenfeldes der Quelle                          C
C          FQ(1,1)         "                                           C
C          B0(0,1,2)  Oberflaechenamplituden                           C
C          D0(0,1,2)           "                                       C
C          F0(1,2)             "                                       C
C                                                                      C
C          Empfaengerabhaengige Variablen                              C
C          K0(E),K1(E)  Faktoren des Winkelanteils,enthalten auch      C
C          K2(E),K3(E)  die Komponenten des Momententensors            C
C          L1(E),L2(E)                                                 C
C          jarg         Argument fuer Besselfunktionen                 C
C          J0,J1,J2     zur Aufnahme der Werte der Besselfunktionen    C
C                       (d_j0,d_j1,d_jn, wobei n=2)                    C
C                                                                      C
C          Integrationsvariablen                                       C
C          In1...In6  Hilfsvariablen zur Berechnung der Integranden    C
C          FFI        =w*w/(4pi*rho(h))  Fernfeld-Integrandenfaktor    C
C          NFI        =w/(r*4pi*rho(h))  Nahfeld-    "                 C
C          Integrandenglieder je Langsamkeit fuer f-te Frequenz        C
C          des E-ten Empfaengers:  /modified 8-2-91 k.koch /           C
C          IDFr       Int.glied der radialen    Fernfeldkomponente     C
C          IDNr          "             "        Nahfeldkomponente      C
C          IDFphi        "        transversalen Fernfeldkomponente     C
C          IDNphi        "             "        Nahfeldkomponente      C
C          IDFz          "        vertikalen    Fernfeldkomponente     C
C          IDNz          "             "        Nahfeldkomponente      C
C                                                                      C
C          normierte Teil-Impulsspektren fuer f-te Frequenz des E-ten  C
C          Empfaengers, fuer Verschiebung in [s*s/(g*cm)]* 1E-20:      C
C          VFr(f,E)   Radiales      Fernfeld-Impulsspektrum            C
C          VNr(f,E)      "           Nahfeld-      "                   C
C          VFphi(f,E) Transversales Fernfeld-      "                   C
C          VNphi(f,E)    "           Nahfeld-      "                   C
C          VFz(f,E)   Vertikales    Fernfeld-      "                   C
C          VNz(f,E)      "           Nahfeld-      "                   C
C                                                                      C
C          End-Spektren(Betraege) fuer Staerke M0:                     C
C          BSFr(f),BSFphi(f),BSFz(f)       / modified 8-2-91 K.Koch /  C
C          BSNr(f),BSNphi(f),BSNz(f)       / modified 8-2-91 K.Koch /  C
C          in [cm*s],[cm],[cm/s]                                       C
C                                                                      C
C                                                                      C
C    Ergebnissvariablen                                                C
C          Vollstaendige Seismogramme  bzw. Spektren                   C
C          SFr(f),SFphi(f),SFz(f)   Fernfeld  /modified 8-2-91 K.Koch/ C
C          SNr(f),SNphi(f),SNz(f)   Nahfeld   /modified 8-2-91 K.Koch/ C
C          fuer Bodenverschiebung in     [cm]         bzw. [cm*s]      C
C          fuer Bodengeschwindigkeit in  [cm/s]       bzw. [cm]        C
C          fuer Bodenbeschleunigung in   [cm/(s*s)]   bzw. [cm/s]      C
C**********************************************************************C