Commit f06fddf2 authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

official qssp2010 version from GFZ web page

parent 90835227
......@@ -9,4 +9,4 @@ Makefile
config.log
config.status
*.o
src/qssp
src/fomosto_qssp2010
AC_INIT([qssp], [1.3])
AC_INIT([fomosto-qssp], [2010])
AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AC_PROG_F77
AC_CONFIG_FILES([
......
# This is the input file of FORTRAN77 program "qssp2010" for calculating
# synthetic seismograms of a self-gravitating, spherically symmetric,
# isotropic and viscoelastic earth.
#
# by
# Rongjiang Wang <wang@gfz-potsdam.de>
# Helmholtz-Centre Potsdam
# GFZ German Reseach Centre for Geosciences
# Telegrafenberg, D-14473 Potsdam, Germany
#
# Last modified: Potsdam, July, 2010
#
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# If not specified, SI Unit System is used overall!
#
# Coordinate systems:
# spherical (r,t,p) with r = radial,
# t = co-latitude,
# p = east longitude.
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
#
# UNIFORM RECEIVER DEPTH
# ======================
# 1. uniform receiver depth [km]
#-------------------------------------------------------------------------------------------
0.00
#-------------------------------------------------------------------------------------------
#
# TIME (FREQUENCY) SAMPLING
# =========================
# 1. time window [sec], sampling interval [sec]
# 2. max. frequency [Hz] of Green's functions
# 3. max. slowness [s/km] of Green's functions
# Note: if the near-field static displacement is desired, the maximum slowness should not
# be smaller than the S wave slowness in the receiver layer
# 4. anti-aliasing factor (> 0 & < 1), if it is <= 0 or >= 1/e (~ 0.4), then
# default value of 1/e is used (e.g., 0.1 = alias phases will be suppressed
# to 10% of their original amplitude)
#
# Note: The computation effort increases linearly the time window and
# quadratically with the cut-off frequency.
#-------------------------------------------------------------------------------------------
200.0 0.25
0.50
0.40
0.25
#-------------------------------------------------------------------------------------------
#
# SELF-GRAVITATING EFFECT
# =======================
# 1. the critical frequency [Hz] and the critical harmonic degree, below which
# the self-gravitating effect should be included
#-------------------------------------------------------------------------------------------
0.0 0
#-------------------------------------------------------------------------------------------
#
# WAVE TYPES
# ==========
# 1. selection (1/0 = yes/no) of speroidal modes (P-SV waves), selection of toroidal modes
# (SH waves), minimum and maximum cutoff harmonic degrees
# Note: if the near-field static displacement is desired, the minimum cutoff harmonic
# degree should not be smaller than, e.g., 2000.
#-------------------------------------------------------------------------------------------
1 1 5000 25000
#-------------------------------------------------------------------------------------------
# GREEN'S FUNCTION FILES
# ======================
# 1. number of discrete source depths, estimated radius of each source patch [km] and
# directory for Green's functions
# 2. list of the source depths [km], the respective file names of the Green's
# functions (spectra) and the switch number (0/1) (0 = do not calculate
# this Green's function because it exists already, 1 = calculate or update
# this Green's function. Note: update is required if any of the above
# parameters is changed)
#-------------------------------------------------------------------------------------------
4 2.5 './GreenFunction/'
2.25 'Green_02km' 0
7.25 'Green_07km' 0
12.25 'Green_12km' 0
17.25 'Green_17km' 0
#-------------------------------------------------------------------------------------------
#
# MULTI-EVENT SOURCE PARAMETERS
# =============================
# 1. number of discrete point sources and selection of the source data format
# (1 or 2)
# 2. list of the multi-event sources
# Format 1:
# M-Unit Mrr Mtt Mpp Mrt Mrp Mtp Lat Lon Depth T_origin T_rise
# [Nm] [deg] [deg] [km] [sec] [sec]
# Format 2:
# Moment Strike Dip Rake Lat Lon Depth T_origin T_rise
# [Nm] [deg] [deg] [deg] [deg] [deg] [km] [sec] [sec]
#-------------------------------------------------------------------------------------------
48 2
1.000E+18 84.7914 90.0000 180.0000 40.7923 27.4302 2.5000 19.5863 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.7923 27.4302 7.5000 19.2511 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.7923 27.4302 12.5000 19.1380 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.7923 27.4302 17.5000 19.2511 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.7968 27.4905 2.5000 17.5144 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.7968 27.4905 7.5000 17.1386 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.7968 27.4905 12.5000 17.0115 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.7968 27.4905 17.5000 17.1386 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8013 27.5508 2.5000 15.4573 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8013 27.5508 7.5000 15.0302 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8013 27.5508 12.5000 14.8851 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8013 27.5508 17.5000 15.0302 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8057 27.6112 2.5000 13.4218 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8057 27.6112 7.5000 12.9276 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8057 27.6112 12.5000 12.7587 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8057 27.6112 17.5000 12.9276 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8275 27.9130 2.5000 4.6779 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8275 27.9130 7.5000 2.9769 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8275 27.9130 12.5000 2.1265 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8275 27.9130 17.5000 2.9769 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8318 27.9734 2.5000 4.1667 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8318 27.9734 7.5000 2.0833 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8318 27.9734 12.5000 0.0000 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8318 27.9734 17.5000 2.0833 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8360 28.0338 2.5000 4.6779 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8360 28.0338 7.5000 2.9769 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8360 28.0338 12.5000 2.1264 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8360 28.0338 17.5000 2.9769 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8402 28.0942 2.5000 5.9538 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8402 28.0942 7.5000 4.7357 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8402 28.0942 12.5000 4.2529 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8402 28.0942 17.5000 4.7357 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8607 28.3964 2.5000 15.4573 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8607 28.3964 7.5000 15.0302 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8607 28.3964 12.5000 14.8851 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8607 28.3964 17.5000 15.0302 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8647 28.4569 2.5000 17.5144 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8647 28.4569 7.5000 17.1386 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8647 28.4569 12.5000 17.0115 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8647 28.4569 17.5000 17.1386 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8687 28.5173 2.5000 19.5863 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8687 28.5173 7.5000 19.2510 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8687 28.5173 12.5000 19.1380 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8687 28.5173 17.5000 19.2510 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8726 28.5778 2.5000 21.6688 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8726 28.5778 7.5000 21.3662 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8726 28.5778 12.5000 21.2644 2.0000
1.000E+18 84.7914 90.0000 180.0000 40.8726 28.5778 17.5000 21.3662 2.0000
#-------------------------------------------------------------------------------------------
#
# RECEIVER PARAMETERS
# ===================
# 1. output file name and and selection of output format:
# 1 = cartesian: vertical(z)/north(n)/east(e);
# 2 = spherical: radial(r)/theta(t)/phi(p)
# (Note: if output format 2 is selected, the epicenter (T_origin = 0)
# 2. output time window [sec] (<= Green's function time window)
# 3. selection of order of Butterworth bandpass filter (if <= 0, then no filtering), lower
# and upper corner frequencies (smaller than the cut-off frequency defined above)
# 4. lower and upper slowness cut-off [s/km] (slowness band-pass filter)
# 5. number of receiver
# 6. list of the station parameters
# Format:
# Lat Lon Name Time_reduction
# [deg] [deg] [sec]
# (Note: Time_reduction = start time of the time window)
#-------------------------------------------------------------------------------------------
'marmara' 1
200.0
3 0.0000 0.40
0.000 0.400
55
40.8455 29.9624 'UCG2' -20.0
40.7867 29.4507 'TUBI' -20.0
40.5655 29.3719 'DUM2' -20.0
40.0975 29.1314 'ULUT' -20.0
40.8521 29.1179 'BAD1' -20.0
40.8197 29.1127 'YANT' -20.0
41.0608 29.0614 'KANT' -20.0
40.5344 28.7820 'BOZT' -20.0
40.9887 28.7239 'AVCT' -20.0
40.2653 28.3326 'KART' -20.0
40.9669 27.9617 'MER1' -20.0
40.3932 27.8079 'ERDT' -20.0
40.0199 27.7150 'CHMT' -20.0
40.1864 27.6958 'ALAT' -20.0
40.6114 27.5869 'MADT' -20.0
40.0830 27.5633 'ATCT' -20.0
40.9507 26.9985 'KRDT' -20.0
40.4683 26.5873 'YENT' -20.0
40.1257 26.5236 'ATHT' -20.0
40.3841 26.4857 'TYF1' -20.0
40.8658 28.9909 'YSST' -20.0
40.8747 28.9735 'SVRT' -20.0
40.9902 27.9810 'BOTS' -20.0
40.8799 29.0688 'BRGA' -20.0
40.8746 29.1295 'BUYA' -20.0
39.3360 26.6663 'CUND' -20.0
39.6041 28.6267 'DST' -20.0
40.6070 26.1406 'ENZZ' -20.0
39.5212 30.8492 'ESKT' -20.0
40.1997 25.9183 'GOKC' -20.0
40.0466 27.6860 'GONE' -20.0
40.2850 30.3167 'GPA' -20.0
39.4552 26.1329 'GPNR' -20.0
40.9908 28.6678 'HVHR' -20.0
40.8780 29.0941 'HYBA' -20.0
40.3368 29.4728 'IZI' -20.0
41.0628 29.0608 'KAVA' -20.0
41.4844 28.3000 'KIYI' -20.0
40.7598 29.3554 'LAFA' -20.0
40.0452 28.8914 'ORLT' -20.0
40.6010 29.7000 'OSMT' -20.0
40.7401 30.3271 'SAUA' -20.0
40.7042 29.1517 'SBT1' -20.0
40.8783 28.5133 'SBT2' -20.0
40.8850 27.9783 'SBT3' -20.0
40.8283 27.5350 'SBT4' -20.0
40.6311 28.8804 'SBT5' -20.0
39.0231 29.2220 'SHAP' -20.0
40.9998 28.5392 'SINB' -20.0
39.1005 28.9805 'SMAA' -20.0
40.9902 27.5357 'TKR' -20.0
40.8128 29.2657 'TUZL' -20.0
40.9908 28.6678 'YKBL' -20.0
40.6953 29.3703 'YLVH' -20.0
41.1288 26.6323 'UKOP' -20.0
#-------------------------------------------------------------------------------------------
#
# LAYERED EARTH MODEL (IASP91)
# ============================
# 1. number of data lines of the layered model and selection for including
# the physical dispersion according Kamamori & Anderson (1977)
#-------------------------------------------------------------------------------------------
29 0
#-------------------------------------------------------------------------------------------
#
# MULTILAYERED MODEL PARAMETERS (source site)
# ===========================================
# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs
#-------------------------------------------------------------------------------------------
1 0.0 5.80000 3.20000 2.60000 1456.0 600.0
2 3.0 5.80000 3.20000 2.60000 1456.0 600.0
3 3.0 5.80000 3.20000 2.60000 1456.0 600.0
4 15.0 5.80000 3.20000 2.60000 1456.0 600.0
5 15.0 6.80000 3.90000 2.90000 1350.0 600.0
6 24.4 6.80000 3.90000 2.90000 1350.0 600.0
7 24.4 8.11061 4.49094 3.38076 1446.0 600.0
8 40.0 8.10119 4.48486 3.37906 1446.0 600.0
9 60.0 8.08907 4.47715 3.37688 1447.0 600.0
10 80.0 8.07689 4.46954 3.37471 1447.0 600.0
11 80.0 8.07689 4.46954 3.37471 195.0 80.0
12 115.0 8.05540 4.45643 3.37091 195.0 80.0
13 150.0 8.03370 4.44361 3.36710 195.0 80.0
14 185.0 8.01180 4.43180 3.36330 195.0 80.0
15 220.0 7.98970 4.41885 3.35950 195.0 80.0
16 220.0 8.55896 4.64391 3.43578 362.0 143.0
17 265.0 8.64552 4.67540 3.46264 365.0 143.0
18 310.0 8.73209 4.70690 3.48951 367.0 143.0
19 355.0 8.81867 4.73840 3.51639 370.0 143.0
20 400.0 8.90522 4.76989 3.54325 372.0 143.0
21 400.0 9.13397 4.93259 3.72378 366.0 143.0
22 450.0 9.38990 5.07842 3.78678 365.0 143.0
23 500.0 9.64588 5.22428 3.84980 364.0 143.0
24 550.0 9.90185 5.37014 3.91282 363.0 143.0
25 600.0 10.15782 5.51601 3.97584 362.0 143.0
26 600.0 10.15782 5.51601 3.97584 362.0 143.0
27 635.0 10.21023 5.54311 3.98399 362.0 143.0
28 670.0 10.26622 5.57020 3.99214 362.0 143.0
29 670.0 10.75131 5.94508 4.38071 759.0 312.0
##############
30 721.0 10.91005 6.09418 4.41241 744.0 312.0
31 771.0 11.06557 6.24046 4.44317 730.0 312.0
32 771.0 11.06557 6.24046 4.44317 730.0 312.0
33 871.0 11.24490 6.31091 4.50372 737.0 312.0
34 971.0 11.41560 6.37813 4.56307 743.0 312.0
35 1071.0 11.57828 6.44232 4.62129 750.0 312.0
36 1171.0 11.73357 6.50370 4.67844 755.0 312.0
37 1271.0 11.88209 6.56250 4.73460 761.0 312.0
38 1371.0 12.02445 6.61891 4.78983 766.0 312.0
39 1471.0 12.16126 6.67317 4.84422 770.0 312.0
40 1571.0 12.29316 6.72548 4.89783 775.0 312.0
41 1671.0 12.42075 6.77606 4.95073 779.0 312.0
42 1771.0 12.54466 6.82512 5.00299 784.0 312.0
43 1871.0 12.66550 6.87289 5.05469 788.0 312.0
44 1971.0 12.78389 6.91957 5.10590 792.0 312.0
45 2071.0 12.90045 6.96538 5.15669 795.0 312.0
46 2171.0 13.01579 7.01053 5.20713 799.0 312.0
47 2271.0 13.13055 7.05525 5.25729 803.0 312.0
48 2371.0 13.24532 7.09974 5.30724 807.0 312.0
49 2471.0 13.36074 7.14423 5.35706 811.0 312.0
50 2571.0 13.47742 7.18892 5.40681 815.0 312.0
51 2671.0 13.59597 7.23403 5.45657 819.0 312.0
52 2741.0 13.68041 7.26597 5.49145 822.0 312.0
53 2771.0 13.68753 7.26575 5.50642 823.0 312.0
54 2871.0 13.71168 7.26486 5.55641 826.0 312.0
55 2891.0 13.71660 7.26466 5.56645 826.0 312.0
56 2891.0 8.06482 0.00000 9.90349 57822.0 0.0
57 2971.0 8.19939 0.00000 10.02940 57822.0 0.0
58 3071.0 8.36019 0.00000 10.18134 57822.0 0.0
59 3171.0 8.51298 0.00000 10.32726 57822.0 0.0
60 3271.0 8.65805 0.00000 10.46727 57822.0 0.0
61 3371.0 8.79573 0.00000 10.60152 57822.0 0.0
62 3471.0 8.92632 0.00000 10.73012 57822.0 0.0
63 3571.0 9.05015 0.00000 10.85321 57822.0 0.0
64 3671.0 9.16752 0.00000 10.97091 57822.0 0.0
65 3771.0 9.27876 0.00000 11.08335 57822.0 0.0
66 3871.0 9.38418 0.00000 11.19067 57822.0 0.0
67 3971.0 9.48409 0.00000 11.29298 57822.0 0.0
68 4071.0 9.57881 0.00000 11.39042 57822.0 0.0
69 4171.0 9.66865 0.00000 11.48311 57822.0 0.0
70 4271.0 9.75393 0.00000 11.57119 57822.0 0.0
71 4371.0 9.83496 0.00000 11.65478 57822.0 0.0
72 4471.0 9.91206 0.00000 11.73401 57822.0 0.0
73 4571.0 9.98554 0.00000 11.80900 57822.0 0.0
74 4671.0 10.05572 0.00000 11.87990 57822.0 0.0
75 4771.0 10.12291 0.00000 11.94682 57822.0 0.0
76 4871.0 10.18743 0.00000 12.00989 57822.0 0.0
77 4971.0 10.24959 0.00000 12.06924 57822.0 0.0
78 5071.0 10.30971 0.00000 12.12500 57822.0 0.0
79 5149.5 10.35568 0.00000 12.16634 57822.0 0.0
80 5149.5 11.02827 3.50432 12.76360 625.0 85.0
81 5171.0 11.03643 3.51002 12.77493 624.0 85.0
82 5271.0 11.07249 3.53522 12.82501 620.0 85.0
83 5371.0 11.10542 3.55823 12.87073 615.0 85.0
84 5471.0 11.13521 3.57905 12.91211 611.0 85.0
85 5571.0 11.16186 3.59767 12.94912 608.0 85.0
86 5671.0 11.18538 3.61411 12.98178 605.0 85.0
87 5771.0 11.20576 3.62835 13.01009 603.0 85.0
88 5871.0 11.22301 3.64041 13.03404 600.0 85.0
89 5971.0 11.23712 3.65027 13.05364 599.0 85.0
90 6071.0 11.24809 3.65794 13.06888 597.0 85.0
91 6171.0 11.25593 3.66342 13.07977 596.0 85.0
92 6271.0 11.26064 3.66670 13.08630 596.0 85.0
93 6321.0 11.26142 3.66725 13.08739 596.0 85.0
94 6371.0 11.26220 3.66780 13.08848 596.0 85.0
#---------------------------------end of all inputs-----------------------------------------
This diff is collapsed.
bin_PROGRAMS = qssp
qssp_SOURCES = brunewvlet.f butterworth.f caxcb.f cdsvd500.f cmemcpy.f disazi.f four1.f moments.f pemain.f qpdifmat2.f qpdifmat4.f qpdifmat6.f qpdlegendre.f qpfftinv.f qpgetinp.f qpgrnspec.f qplegendre.f qpmain.f qppsvkern.f qppsvkerng.f qpqmodel.f qpshkern.f qpsmat0.f qpsmatc.f qpsmat.f qpsprop0.f qpsprop1.f qpsprop.f qpspropg0.f qpspropg1.f qpspropg.f qpstart0.f qpstart4.f qpstart6.f qpsublayer.f qptmat.f qptprop1.f qptprop.f qpwvint.f ruku.f spbdphj.f spbdphy.f spbjh.f spbphj.f spbphy.f spbpsj.f spbpsy.f taper.f qpglobal.h skip_comments.f
bin_PROGRAMS = fomosto_qssp2010
fomosto_qssp2010_SOURCES = brunewvlet.f butterworth.f caxcb.f cdsvd500.f cmemcpy.f \
disazi.f four1.f moments.f pemain.f qpdifmat2.f qpdifmat4.f \
qpdifmat6.f qpdlegendre.f qpfftinv.f qpgetinp.f qpgrnspec.f \
qplegendre.f qpmain.f qppsvkern.f qppsvkerng.f qpqmodel.f \
qpshkern.f qpsmat0.f qpsmatc.f qpsmat.f qpsprop0.f qpsprop1.f \
qpsprop.f qpspropg0.f qpspropg1.f qpspropg.f qpstart0.f \
qpstart4.f qpstart6.f qpsublayer.f qptmat.f qptprop1.f \
qptprop.f qpwvint.f ruku.f spbdphj.f spbdphy.f spbjh.f \
spbphj.f spbphy.f spbpsj.f spbpsy.f taper.f banddpasss.f \
qptprop1st.f qpdifmat0.f qpdifmatl.f qpdifmats.f swavelet.f \
qpglobal.h skip_comments.f
subroutine bandpass(n,f1,f2,df,nf,h)
implicit none
integer n,nf
double precision f1,f2,df
double complex h(nf)
c
integer l,k
double precision w,w1,w2,delt,arg
c
double precision PI
data PI/3.14159265358979d0/
c
if(n.le.0.or.f1.ge.f2)then
do l=1,nf
h(l)=(1.d0,0.d0)
enddo
else if(f1.le.0.d0)then
w2=2.d0*PI*f2
do l=1,nf
w=2.d0*PI*dble(l-1)*df
delt=w2-w1
h(l)=(1.d0,0.d0)
do k=1,n
arg=PI*dble(2*k-1)/dble(2*n)
h(l)=h(l)*dcmplx(0.d0,-delt)
& /dcmplx(w-delt*dcos(arg),-delt*dsin(arg))
enddo
enddo
else
w1=2.d0*PI*f1
w2=2.d0*PI*f2
do l=1,nf
w=2.d0*PI*dble(l-1)*df
delt=w*(w2-w1)
h(l)=(1.d0,0.d0)
do k=1,n
arg=PI*dble(2*k-1)/dble(2*n)
h(l)=h(l)*dcmplx(0.d0,-delt)
& /dcmplx(w*w-w1*w2-delt*dcos(arg),-delt*dsin(arg))
enddo
enddo
endif
return
end
subroutine d2dfit(n,dis,p,b,ndeg,disfit)
implicit none
integer n,ndeg
double precision dis(n),disfit(n),p(n,0:ndeg),b(0:ndeg)
c
integer i,ideg
double precision x,dx
c
dx=2.d0/dble(n-1)
do i=1,n
x=-1.d0+dble(i-1)*dx
p(i,0)=1.d0
p(i,1)=x
do ideg=2,ndeg
p(i,ideg)=(dble(2*ideg-1)*x*p(i,ideg-1)
& -dble(ideg-1)*p(i,ideg-2))/dble(ideg)
enddo
enddo
c
do ideg=0,ndeg
b(ideg)=0.5d0*(dis(1)*p(1,ideg)+dis(n)*p(n,ideg))
do i=2,n-1
b(ideg)=b(ideg)+dis(i)*p(i,ideg)
enddo
b(ideg)=b(ideg)*dx*0.5d0*dble(2*ideg+1)
enddo
c
do i=1,n
x=-1.+dble(i-1)*dx
disfit(i)=0.d0
do ideg=0,ndeg
disfit(i)=disfit(i)+b(ideg)*p(i,ideg)
enddo
enddo
return
end
\ No newline at end of file
subroutine qpdifmat0(ly,ldeg,rr,mat)
implicit none
integer ly,ldeg
double precision rr
double complex mat(6,6)
c
include 'qpglobal.h'
c
c 3x3 coefficient matrix for spheroidal mode l = 0
c
double precision mass,rorr,beta,dro,ro1
double complex cup,clw,crr,crorr,cxirr,clarr,cmurr,cgrrr,cvprr
double complex cgarr,gamma
c
double complex c1,c2,c4
data c1,c2,c4/(1.d0,0.d0),(2.d0,0.d0),(4.d0,0.d0)/
c
crr=dcmplx(rr,0.d0)
cup=(crr-crrlw(ly))/(crrup(ly)-crrlw(ly))
clw=c1-cup
clarr=cup*claup(ly)+clw*clalw(ly)
cmurr=cup*cmuup(ly)+clw*cmulw(ly)
cvprr=cup*cvpup(ly)+clw*cvplw(ly)
cxirr=clarr+c2*cmurr
beta=dlog(roup(ly)/rolw(ly))/(rrup(ly)-rrlw(ly))
c
if(rr.le.rrlw(ly))then
mass=0.d0
crorr=crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
else if(ly.ge.lyos)then
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
dro=(rorr-rolw(ly))/(rr-rrlw(ly))
ro1=rorr-dro*rrlw(ly)
mass=PI*(rr-rrlw(ly))*((4.d0/3.d0)*ro1
& *(rr**2+rr*rrlw(ly)+rrlw(ly)**2)
& +dro*(rr**3+rr**2*rrlw(ly)
& +rr*rrlw(ly)**2+rrlw(ly)**3))
else
rorr=rolw(ly)*dexp(dlog(roup(ly)/rolw(ly))
& *(rr-rrlw(ly))/(rrup(ly)-rrlw(ly)))
crorr=dcmplx(rorr,0.d0)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
mass=2.d0*PI2/beta**3
& *(rorr*(beta*rr*(beta*rr-2.d0)+2.d0)
& -rolw(ly)*(beta*rrlw(ly)*(beta*rrlw(ly)-2.d0)+2.d0))
endif
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG*mass/rr**2,0.d0)
c
c stabilization at the static limit
c
if(ly.lt.lyos)then
if(cdabs(comi)/PI2.lt.fbvatm)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvatm))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
else if(ly.ge.lyos.and.ly.lt.lyob)then
if(cdabs(comi)/PI2.lt.fbvocean)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvocean))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
else if(ly.ge.lycm.and.ly.lt.lycc)then
if(cdabs(comi)/PI2.lt.fbvcore)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvcore))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
endif
c
mat(1,1)=(c1-c2*clarr/cxirr)/crr
mat(1,2)=c1/(crorr*cvprr**2)/crr
mat(1,3)=(0.d0,0.d0)
c
mat(2,1)=c4*(cmurr*(c1+c2*clarr/cxirr)/crr-crorr*cgrrr)
& -crorr*comi2*crr
mat(2,2)=c2*clarr/cxirr/crr
mat(2,3)=(0.d0,0.d0)
c
mat(3,1)=cgarr/crr
mat(3,2)=(0.d0,0.d0)
mat(3,3)=(0.d0,0.d0)
c
return
end
\ No newline at end of file
subroutine qpdifmatl(ly,ldeg,rr,mat)
implicit none
integer ly,ldeg
double precision rr
double complex mat(6,6)
c
include 'qpglobal.h'
c
c 4x4 coefficient matrix for spheroidal mode l > 0 in liquid media
c
double precision mass,rorr,beta,dro,ro1
double complex cldeg,clp1,cll1
double complex cup,clw,crr,crorr,cgrrr,cvprr
double complex cgarr,gamma,cn2rr
logical bv
c
double complex c1,c2,c4
data c1,c2,c4/(1.d0,0.d0),(2.d0,0.d0),(4.d0,0.d0)/
c
cldeg=dcmplx(dble(ldeg),0.d0)
clp1=cldeg+c1
cll1=cldeg*clp1
c
crr=dcmplx(rr,0.d0)
cup=(crr-crrlw(ly))/(crrup(ly)-crrlw(ly))
clw=c1-cup
cvprr=cup*cvpup(ly)+clw*cvplw(ly)
beta=dlog(roup(ly)/rolw(ly))/(rrup(ly)-rrlw(ly))
c
if(rr.le.rrlw(ly))then
mass=0.d0
crorr=crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
else if(ly.ge.lyos)then
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
dro=(rorr-rolw(ly))/(rr-rrlw(ly))
ro1=rorr-dro*rrlw(ly)
mass=PI*(rr-rrlw(ly))*((4.d0/3.d0)*ro1
& *(rr**2+rr*rrlw(ly)+rrlw(ly)**2)
& +dro*(rr**3+rr**2*rrlw(ly)
& +rr*rrlw(ly)**2+rrlw(ly)**3))
else
rorr=rolw(ly)*dexp(dlog(roup(ly)/rolw(ly))
& *(rr-rrlw(ly))/(rrup(ly)-rrlw(ly)))
crorr=dcmplx(rorr,0.d0)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
mass=2.d0*PI2/beta**3
& *(rorr*(beta*rr*(beta*rr-2.d0)+2.d0)
& -rolw(ly)*(beta*rrlw(ly)*(beta*rrlw(ly)-2.d0)+2.d0))
endif
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG*mass/rr**2,0.d0)
c
c stabilization at the static limit
c
if(ly.lt.lyos)then
bv=fbvatm.gt.0.d0
if(cdabs(comi)/PI2.lt.fbvatm)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvatm))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma