mirror of
https://github.com/HongjianFang/DSurfTomo.git
synced 2025-05-09 00:51:12 +08:00
Compare commits
No commits in common. "stable" and "v1.2" have entirely different histories.
7
.gitignore
vendored
7
.gitignore
vendored
@ -1,7 +0,0 @@
|
|||||||
*.o
|
|
||||||
*.mod
|
|
||||||
Applications
|
|
||||||
example_*
|
|
||||||
.gitignore
|
|
||||||
temp
|
|
||||||
DSurfTomo
|
|
||||||
21
README.md
21
README.md
@ -1,8 +1,15 @@
|
|||||||
DSurfTomo is a surface wave tomography program which inverts surface wave dispersion data directly to 3D shear wavespeed models without the intermediate step of constructing the phase or group velocity maps.
|
#DSurfTomo
|
||||||
The fast marching method (FMM) (Rawlinson et al., 2004) is used to compute, at each period, surface wave travel times and ray paths between sources and receivers. This avoids the assumption of great-circle propagation that is used in most surface wave tomographic studies, but which is not appropriate in complex media.
|
DSurfTomo is the surface wave inversion program which can directly invert surface wave
|
||||||
To show its usage, an example of Taipei Basin tomography is provided, including scripts that are used for data reformating and plotting results. Interested users are recommended to refer to Fang et al. (2015 , GJI) for the detail description of the method.
|
dispersion data to 3D shear wave speed without the intermediate step of constructing
|
||||||
|
the phase or group velocity maps.A fast march method (FMM) is used to compute, at each
|
||||||
Fang, H., Yao, H., Zhang, H., Huang, Y. C., & van der Hilst, R. D. (2015). Direct inversion of surface wave dispersion for three-dimensional shallow crustal structure based on ray tracing: methodology and application. Geophysical Journal International, 201(3), 1251-1263.
|
period, surface wave travel times and ray paths between sources and receivers. This
|
||||||
|
avoids the assumption of great-circle propagation that is used in most surface wave
|
||||||
Rawlinson, N. & Sambridge, M., 2004. Wave front evolution in strongly heterogeneous layered media using the fast marching method, Geophys. J. Int., 156(3), 631–647
|
tomographic studies, but which is not appropriate in complex media.
|
||||||
|
For detail description of the method, please refer to:
|
||||||
|
Fang, H., Yao, H., Zhang, H., Huang, Y. C., & van der Hilst, R. D. (2015).
|
||||||
|
Direct inversion of surface wave dispersion for three-dimensional shallow
|
||||||
|
crustal structure based on ray tracing: methodology and application. Geophysical Journal International, 201(3), 1251-1263.
|
||||||
|
|
||||||
|
This is the source code without users' manual and examples about how to use it. If you need a detail description and some examples,
|
||||||
|
feel free to send an email to Hongjian Fang (fanghj@mail.ustc.edu.cn) with a brief introduction about yourself and why you want to
|
||||||
|
use it. We will be happy to help.
|
||||||
|
|||||||
File diff suppressed because it is too large
Load Diff
@ -1 +0,0 @@
|
|||||||
../src/DSurfTomo
|
|
||||||
BIN
bin/DSurfTomo
Executable file
BIN
bin/DSurfTomo
Executable file
Binary file not shown.
13
configure
vendored
13
configure
vendored
@ -1,17 +1,16 @@
|
|||||||
#!/usr/bin/env python
|
#!/usr/bin/env python
|
||||||
# this is a script to install surf_tomo and surf_tomo_syn in your system
|
# this is a script to install surf_tomo and surf_tomo_syn in your system
|
||||||
# written by Hongjian Fang(fanghj@mail.ustc.edu.cn)
|
# written by Hongjian Fang(fanghj@mail.ustc.edu.cn)
|
||||||
|
|
||||||
import os
|
import os
|
||||||
if 'bin' in os.listdir('.'):
|
if 'bin' in os.listdir('.'):
|
||||||
print ('Installation beginning')
|
print 'installation beginning'
|
||||||
else:
|
else:
|
||||||
print ('Installation beginning')
|
print 'installation beginning'
|
||||||
os.mkdir('bin')
|
os.mkdir('bin')
|
||||||
os.chdir('src')
|
os.chdir('src')
|
||||||
os.system('make clean')
|
os.system('make clean')
|
||||||
os.system('make')
|
os.system('make')
|
||||||
os.system('cp DSurfTomo ../bin')
|
os.system('cp DSurfTomo ../bin')
|
||||||
print ('--------------------------------------')
|
print '--------------------------------------'
|
||||||
print ('DSurfTomo compiling Finished')
|
print 'Finishing DSurfTomo compiling'
|
||||||
print ('--------------------------------------')
|
print '--------------------------------------'
|
||||||
|
|||||||
BIN
doc/ManualDSurfTomo.pdf
Normal file
BIN
doc/ManualDSurfTomo.pdf
Normal file
Binary file not shown.
Binary file not shown.
Binary file not shown.
10
example/.gmtcommands4
Normal file
10
example/.gmtcommands4
Normal file
@ -0,0 +1,10 @@
|
|||||||
|
# GMT common arguments shelf
|
||||||
|
-Ba0.2f0.1:S velocity (km/s):
|
||||||
|
-JM2i
|
||||||
|
-R121.35/121.605/24.975/25.2
|
||||||
|
-x2.3i
|
||||||
|
-Y5i
|
||||||
|
-y-2.7i
|
||||||
|
-jM2i
|
||||||
|
EOF
|
||||||
|
|
||||||
107
example/.gmtdefaults4
Normal file
107
example/.gmtdefaults4
Normal file
@ -0,0 +1,107 @@
|
|||||||
|
#
|
||||||
|
# GMT-SYSTEM 4.5.15 [64-bit] Defaults file
|
||||||
|
#
|
||||||
|
#-------- Plot Media Parameters -------------
|
||||||
|
PAGE_COLOR = white
|
||||||
|
PAGE_ORIENTATION = landscape
|
||||||
|
PAPER_MEDIA = a4
|
||||||
|
#-------- Basemap Annotation Parameters ------
|
||||||
|
ANNOT_MIN_ANGLE = 20
|
||||||
|
ANNOT_MIN_SPACING = 0
|
||||||
|
ANNOT_FONT_PRIMARY = Helvetica
|
||||||
|
ANNOT_FONT_SIZE_PRIMARY = 6p
|
||||||
|
ANNOT_OFFSET_PRIMARY = 0.2c
|
||||||
|
ANNOT_FONT_SECONDARY = Helvetica
|
||||||
|
ANNOT_FONT_SIZE_SECONDARY = 16p
|
||||||
|
ANNOT_OFFSET_SECONDARY = 0.2c
|
||||||
|
DEGREE_SYMBOL = ring
|
||||||
|
HEADER_FONT = Helvetica
|
||||||
|
HEADER_FONT_SIZE = 36p
|
||||||
|
HEADER_OFFSET = 0.5c
|
||||||
|
LABEL_FONT = Helvetica
|
||||||
|
LABEL_FONT_SIZE = 5p
|
||||||
|
LABEL_OFFSET = 0.1c
|
||||||
|
OBLIQUE_ANNOTATION = 1
|
||||||
|
PLOT_CLOCK_FORMAT = hh:mm:ss
|
||||||
|
PLOT_DATE_FORMAT = yyyy-mm-dd
|
||||||
|
PLOT_DEGREE_FORMAT = ddd:mm:ss
|
||||||
|
Y_AXIS_TYPE = hor_text
|
||||||
|
#-------- Basemap Layout Parameters ---------
|
||||||
|
BASEMAP_AXES = WESN
|
||||||
|
BASEMAP_FRAME_RGB = black
|
||||||
|
BASEMAP_TYPE = fancy
|
||||||
|
FRAME_PEN = 1.25p
|
||||||
|
FRAME_WIDTH = 0.1c
|
||||||
|
GRID_CROSS_SIZE_PRIMARY = 0c
|
||||||
|
GRID_PEN_PRIMARY = 0.25p
|
||||||
|
GRID_CROSS_SIZE_SECONDARY = 0c
|
||||||
|
GRID_PEN_SECONDARY = 0.5p
|
||||||
|
MAP_SCALE_HEIGHT = 0.2c
|
||||||
|
POLAR_CAP = 85/90
|
||||||
|
TICK_LENGTH = 0.1c
|
||||||
|
TICK_PEN = 0.3p
|
||||||
|
X_AXIS_LENGTH = 25c
|
||||||
|
Y_AXIS_LENGTH = 15c
|
||||||
|
X_ORIGIN = 2.5c
|
||||||
|
Y_ORIGIN = 2.5c
|
||||||
|
UNIX_TIME = FALSE
|
||||||
|
UNIX_TIME_POS = BL/-2c/-2c
|
||||||
|
UNIX_TIME_FORMAT = %Y %b %d %H:%M:%S
|
||||||
|
#-------- Color System Parameters -----------
|
||||||
|
COLOR_BACKGROUND = black
|
||||||
|
COLOR_FOREGROUND = white
|
||||||
|
COLOR_NAN = 128
|
||||||
|
COLOR_IMAGE = adobe
|
||||||
|
COLOR_MODEL = rgb
|
||||||
|
HSV_MIN_SATURATION = 1
|
||||||
|
HSV_MAX_SATURATION = 0.1
|
||||||
|
HSV_MIN_VALUE = 0.3
|
||||||
|
HSV_MAX_VALUE = 1
|
||||||
|
#-------- PostScript Parameters -------------
|
||||||
|
CHAR_ENCODING = ISOLatin1+
|
||||||
|
DOTS_PR_INCH = 300
|
||||||
|
GLOBAL_X_SCALE = 1
|
||||||
|
GLOBAL_Y_SCALE = 1
|
||||||
|
N_COPIES = 1
|
||||||
|
PS_COLOR = rgb
|
||||||
|
PS_IMAGE_COMPRESS = lzw
|
||||||
|
PS_IMAGE_FORMAT = ascii
|
||||||
|
PS_LINE_CAP = butt
|
||||||
|
PS_LINE_JOIN = miter
|
||||||
|
PS_MITER_LIMIT = 35
|
||||||
|
PS_VERBOSE = FALSE
|
||||||
|
TRANSPARENCY = 0
|
||||||
|
#-------- I/O Format Parameters -------------
|
||||||
|
D_FORMAT = %.12lg
|
||||||
|
FIELD_DELIMITER = tab
|
||||||
|
GRIDFILE_FORMAT = nf
|
||||||
|
GRIDFILE_SHORTHAND = FALSE
|
||||||
|
INPUT_CLOCK_FORMAT = hh:mm:ss
|
||||||
|
INPUT_DATE_FORMAT = yyyy-mm-dd
|
||||||
|
IO_HEADER = FALSE
|
||||||
|
N_HEADER_RECS = 1
|
||||||
|
NAN_RECORDS = pass
|
||||||
|
OUTPUT_CLOCK_FORMAT = hh:mm:ss
|
||||||
|
OUTPUT_DATE_FORMAT = yyyy-mm-dd
|
||||||
|
OUTPUT_DEGREE_FORMAT = D
|
||||||
|
XY_TOGGLE = FALSE
|
||||||
|
#-------- Projection Parameters -------------
|
||||||
|
ELLIPSOID = WGS-84
|
||||||
|
MAP_SCALE_FACTOR = default
|
||||||
|
MEASURE_UNIT = cm
|
||||||
|
#-------- Calendar/Time Parameters ----------
|
||||||
|
TIME_FORMAT_PRIMARY = full
|
||||||
|
TIME_FORMAT_SECONDARY = full
|
||||||
|
TIME_EPOCH = 2000-01-01T12:00:00
|
||||||
|
TIME_IS_INTERVAL = OFF
|
||||||
|
TIME_INTERVAL_FRACTION = 0.5
|
||||||
|
TIME_LANGUAGE = us
|
||||||
|
TIME_UNIT = d
|
||||||
|
TIME_WEEK_START = Sunday
|
||||||
|
Y2K_OFFSET_YEAR = 1950
|
||||||
|
#-------- Miscellaneous Parameters ----------
|
||||||
|
HISTORY = TRUE
|
||||||
|
INTERPOLANT = akima
|
||||||
|
LINE_STEP = 0.025c
|
||||||
|
VECTOR_SHAPE = 0
|
||||||
|
VERBOSE = FALSE
|
||||||
@ -1,21 +1,21 @@
|
|||||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||||
c INPUT PARAMETERS
|
c INPUT PARAMETERS
|
||||||
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||||
surfdataTB.dat c: data file
|
surfdataTB.dat c: data file
|
||||||
18 18 9 c: nx ny nz (grid number in lat lon and depth direction)
|
18 18 9 c: nx ny nz (grid number in lat lon and depth direction)
|
||||||
25.2 121.35 c: goxd gozd (upper left point,[lat,lon])
|
25.2 121.35 c: goxd gozd (upper left point,[lat,lon])
|
||||||
0.015 0.017 c: dvxd dvzd (grid interval in lat and lon direction)
|
0.015 0.017 c: dvxd dvzd (grid interval in lat and lon direction)
|
||||||
20 c: max(sources, receivers)
|
25 c: nsrc*maxf
|
||||||
4.0 1.0 c: weight damp
|
4.0 0.0 c: weight damp
|
||||||
3 c: sablayers (for computing depth kernel, 2~5)
|
3 c: nsublayer (numbers of sublayers for each grid interval:grid --> layer)
|
||||||
0.5 2.8 c: minimum velocity, maximum velocity (a priori information)
|
0.5 2.8 c: minimum velocity, maximum velocity
|
||||||
10 c: maximum iteration
|
10 c: maxiter (iteration number)
|
||||||
0.2 c: sparsity fraction
|
0.2 c: sparsity fraction
|
||||||
26 c: kmaxRc (followed by periods)
|
26 c: kmaxRc (followed by periods)
|
||||||
0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0
|
0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0
|
||||||
0 c: kmaxRg
|
0 c: kmaxRg
|
||||||
0 c: kmaxLc
|
0 c: kmaxLc
|
||||||
0 c: kmaxLg
|
0 c: kmaxLg
|
||||||
0 c: synthetic flag(0:real data,1:synthetic)
|
0 c: synthetic flag(0:real data,1:synthetic)
|
||||||
0.02 c: noiselevel
|
0.02 c: noiselevel
|
||||||
3.0 c: threshold
|
2.5 c: threshold
|
||||||
49
example/plotslicenewVs.gmt
Normal file
49
example/plotslicenewVs.gmt
Normal file
@ -0,0 +1,49 @@
|
|||||||
|
#!/bin/csh
|
||||||
|
|
||||||
|
#-----------------------------------------------------------
|
||||||
|
# 2015-06-08 Hongjian Fang
|
||||||
|
# example:csh plotslicenewVs.gmt DSurfTomo.inMeasure.dat 0.2 0.4 0.8 1.1
|
||||||
|
#-----------------------------------------------------------
|
||||||
|
|
||||||
|
gmtset ANOT_FONT_SIZE 6
|
||||||
|
gmtset LABEL_FONT_SIZE 6
|
||||||
|
gmtset FRAME_WIDTH=0.1c
|
||||||
|
gmtset LABEL_OFFSET=0.1c
|
||||||
|
gmtset LABEL_FONT_SIZE=5p
|
||||||
|
gmtset TICK_LENGTH=0.1c
|
||||||
|
gmtset TICK_PEN=0.3p
|
||||||
|
set inp3D = $1
|
||||||
|
set ps = horizVs.ps
|
||||||
|
set J = -JM2i #size for plot
|
||||||
|
set cpt = slice.cpt
|
||||||
|
|
||||||
|
# start
|
||||||
|
makecpt -Cseis -T0.6/1.5/0.1 > $cpt #velocity boundary
|
||||||
|
set RMAP = `minmax -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
||||||
|
psbasemap $RMAP $J -Ba0.1f0.05WseN -P -K -Y5i >$ps
|
||||||
|
awk '{if($3==depth1) print $1,$2,$4}' depth1=$2 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||||
|
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||||
|
rm -rf tmp.grd
|
||||||
|
psscale -Cslice.cpt -Ba0.1f0.05:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
||||||
|
|
||||||
|
makecpt -Cseis -T0.8/1.8/0.1 > $cpt #velocity boundary
|
||||||
|
psbasemap $RMAP $J -Ba0.1f0.05NwsE -P -K -O -X2.3i >>$ps
|
||||||
|
awk '{if($3==depth2) print $1,$2,$4}' depth2=$3 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||||
|
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||||
|
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
||||||
|
|
||||||
|
makecpt -Cseis -T1.1/2.0/0.1 > $cpt #velocity boundary
|
||||||
|
psbasemap $RMAP $J -Ba0.1f0.05WneS -P -K -O -Y-2.7i -X-2.3i >>$ps
|
||||||
|
awk '{if($3==depth2) print $1,$2,$4}' depth2=$4 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||||
|
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||||
|
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -K -P >> $ps
|
||||||
|
|
||||||
|
makecpt -Cseis -T1.3/2.4/0.1 > $cpt #velocity boundary
|
||||||
|
psbasemap $RMAP $J -Ba0.1f0.05SwnE -P -K -O -X2.3i >>$ps
|
||||||
|
awk '{if($3==depth2) print $1,$2,$4}' depth2=$5 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||||
|
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||||
|
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -P >> $ps
|
||||||
|
|
||||||
|
rm *cpt *grd
|
||||||
|
ps2eps -f $ps
|
||||||
|
evince *eps &
|
||||||
@ -1,146 +0,0 @@
|
|||||||
#/* -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
|
|
||||||
#
|
|
||||||
# File Name : GenerateDSurfTomoInputFile.py
|
|
||||||
#
|
|
||||||
# Purpose : Generate the Input File for DSurfTomo using the data themselves
|
|
||||||
#
|
|
||||||
# Creation Date : 05-07-2023
|
|
||||||
#
|
|
||||||
# Last Modified : Tue 18 Jul 2023 08:10:22 PM CST
|
|
||||||
#
|
|
||||||
# Created By : Hongjian Fang: fanghj1990@gmail.com
|
|
||||||
#
|
|
||||||
#_._._._._._._._._._._._._._._._._._._._._.*/
|
|
||||||
|
|
||||||
import os
|
|
||||||
import numpy as np
|
|
||||||
import glob
|
|
||||||
import pandas as pd
|
|
||||||
|
|
||||||
#These are for generating the right format
|
|
||||||
DataPath = 'TaipeiRawData'
|
|
||||||
DataNames = 'CDisp*.dat'
|
|
||||||
wavetype = 2
|
|
||||||
phasetype = 0
|
|
||||||
|
|
||||||
# These following can be changed in DSurfTomo.in
|
|
||||||
nz = 9
|
|
||||||
sublayers = 3
|
|
||||||
minvel = 0.5
|
|
||||||
maxvel = 2.8
|
|
||||||
noiselevel = 0.02
|
|
||||||
datathreshold = 1.5
|
|
||||||
DataFiles = glob.glob(DataPath+'/'+DataNames)
|
|
||||||
|
|
||||||
sta1latall = []
|
|
||||||
sta1lonall = []
|
|
||||||
sta2latall = []
|
|
||||||
sta2lonall = []
|
|
||||||
periodsall = []
|
|
||||||
dispersionall = []
|
|
||||||
stationpairID = []
|
|
||||||
|
|
||||||
for ifile in DataFiles:
|
|
||||||
dispersion = np.loadtxt(ifile)
|
|
||||||
pairID = ifile.split('/')[-1]
|
|
||||||
sta1lon = dispersion[0,0]
|
|
||||||
sta1lat = dispersion[0,1]
|
|
||||||
sta2lon = dispersion[1,0]
|
|
||||||
sta2lat = dispersion[1,1]
|
|
||||||
nperiods,_ = dispersion.shape
|
|
||||||
nperiods = nperiods - 2
|
|
||||||
periods = np.zeros(nperiods,)
|
|
||||||
disper = np.zeros(nperiods,)
|
|
||||||
for ii in range(nperiods):
|
|
||||||
periods[ii] = dispersion[ii+2,0]
|
|
||||||
disper[ii] = dispersion[ii+2,1]
|
|
||||||
|
|
||||||
dispersionall = np.hstack([dispersionall,disper])
|
|
||||||
periodsall = np.hstack([periodsall,periods])
|
|
||||||
sta1latall = np.hstack([sta1latall,np.ones(nperiods,)*sta1lat])
|
|
||||||
sta1lonall = np.hstack([sta1lonall,np.ones(nperiods,)*sta1lon])
|
|
||||||
sta2latall = np.hstack([sta2latall,np.ones(nperiods,)*sta2lat])
|
|
||||||
sta2lonall = np.hstack([sta2lonall,np.ones(nperiods,)*sta2lon])
|
|
||||||
stationpairID = stationpairID + [pairID] * nperiods
|
|
||||||
|
|
||||||
dataall = pd.DataFrame({'sta1lat':sta1latall, 'sta1lon':sta1lonall, \
|
|
||||||
'sta2lat':sta2latall, 'sta2lon':sta2lonall, \
|
|
||||||
'periods':periodsall, 'dispersion': dispersionall, \
|
|
||||||
'pairid':stationpairID})
|
|
||||||
|
|
||||||
dataall = dataall.sort_values(by = ['periods', 'pairid'])
|
|
||||||
dataall = dataall.set_index('periods')
|
|
||||||
with open('surfdata.dat','w') as fout:
|
|
||||||
UniqPeriods = dataall.index.unique()
|
|
||||||
for iperiod,period in enumerate(UniqPeriods):
|
|
||||||
datasubperiod = dataall.loc[period]
|
|
||||||
if isinstance(datasubperiod,pd.Series):
|
|
||||||
continue
|
|
||||||
datasubperiod = datasubperiod.reset_index().set_index('sta1lat')
|
|
||||||
sta1lat = datasubperiod.index.unique()
|
|
||||||
for ista1 in sta1lat:
|
|
||||||
datasubstation = datasubperiod.loc[ista1]
|
|
||||||
if isinstance(datasubstation,pd.DataFrame):
|
|
||||||
fout.write(f'# {datasubstation.index[0]} {datasubstation["sta1lon"].iloc[0]} {iperiod+1} {wavetype} {phasetype}\n')
|
|
||||||
for ista2 in range(len(datasubstation)):
|
|
||||||
fout.write(f'{datasubstation["sta2lat"].iloc[ista2]} {datasubstation["sta2lon"].iloc[ista2]} {datasubstation["dispersion"].iloc[ista2]}\n')
|
|
||||||
else:
|
|
||||||
fout.write(f'# {datasubstation.name} {datasubstation["sta1lon"]} {iperiod+1} {wavetype} {phasetype}\n')
|
|
||||||
fout.write(f'{datasubstation["sta2lat"]} {datasubstation["sta2lon"]} {datasubstation["dispersion"]}\n')
|
|
||||||
|
|
||||||
|
|
||||||
print('Finish reformatting disperison data')
|
|
||||||
|
|
||||||
# Determin the grid interval, 1/3 of the smallest station distance
|
|
||||||
|
|
||||||
# Weight and Damp, you will have to run the code a couple of time to decide
|
|
||||||
# The default values is 5.0 and 1.0, respectively
|
|
||||||
weight = 1.0
|
|
||||||
damp = 1.0
|
|
||||||
|
|
||||||
# Sparsity fraction, this parameter represent the sparsity of the sensitivity matrix and is set to be 0.2 for safe keeping, but could be as low as 0.02.it is not important as long as your code does not throw you memory erros
|
|
||||||
sparsityfraction = 0.1
|
|
||||||
|
|
||||||
# Iteration number, the default is 10 time, most often you get stable results after 4-5 iterations.
|
|
||||||
maxiteration = 10
|
|
||||||
|
|
||||||
|
|
||||||
distall = np.sqrt((dataall['sta2lat']-dataall['sta1lat'])**2+(dataall['sta2lon']-dataall['sta1lon'])**2)
|
|
||||||
mindist = distall.min()
|
|
||||||
gridintval = int(np.ceil(mindist/3*1000))/1000
|
|
||||||
# Determine the origin, 1 grids
|
|
||||||
originLat = np.max([np.max(dataall.sta1lat.max()),np.max(dataall.sta2lat.max())]) + 1*gridintval
|
|
||||||
originLon = np.min([np.min(dataall.sta1lon.min()),np.min(dataall.sta2lon.min())]) - 1*gridintval
|
|
||||||
largestLat = np.max([dataall['sta1lat'].max()-dataall['sta1lat'].min(),dataall['sta2lat'].max()-dataall['sta2lat'].min()])
|
|
||||||
largestLon = np.max([dataall['sta1lon'].max()-dataall['sta1lon'].min(),dataall['sta2lon'].max()-dataall['sta2lon'].min()])
|
|
||||||
inputfile = open('DSurfTomo.in','w')
|
|
||||||
dummy = 'cccc'
|
|
||||||
inputfile.write(f'{dummy*10}\n')
|
|
||||||
inputfile.write(f'{dummy*10}\n')
|
|
||||||
inputfile.write(f'{dummy*10}\n')
|
|
||||||
inputfile.write('surfdata.dat\n')
|
|
||||||
nx = int(np.ceil(largestLat/gridintval))+5
|
|
||||||
ny = int(np.ceil(largestLon/gridintval))+5
|
|
||||||
nreceivers = len(dataall.sta1lat.unique())+1
|
|
||||||
inputfile.write(f'{nx} {ny} {nz}\n')
|
|
||||||
inputfile.write(f'{originLat:<7.3f} {originLon:7.3f}\n')
|
|
||||||
inputfile.write(f'{gridintval:<7.3f} {gridintval:7.3f}\n')
|
|
||||||
inputfile.write(f'{nreceivers}\n')
|
|
||||||
inputfile.write(f'{weight} {damp}\n')
|
|
||||||
inputfile.write(f'{sublayers}\n')
|
|
||||||
inputfile.write(f'{minvel} {maxvel}\n')
|
|
||||||
inputfile.write(f'{maxiteration}\n')
|
|
||||||
inputfile.write(f'{sparsityfraction}\n')
|
|
||||||
nperiods = len(UniqPeriods)
|
|
||||||
inputfile.write(f'{nperiods}\n')
|
|
||||||
#inputfile.write(f'{*UniqPeriods}\n')
|
|
||||||
inputfile.write(' '.join([str(iperiod) for iperiod in UniqPeriods])+'\n')
|
|
||||||
inputfile.write(f'0\n')
|
|
||||||
inputfile.write(f'0\n')
|
|
||||||
inputfile.write(f'0\n')
|
|
||||||
inputfile.write(f'0\n')
|
|
||||||
inputfile.write(f'{noiselevel}\n')
|
|
||||||
inputfile.write(f'{datathreshold}\n')
|
|
||||||
inputfile.close()
|
|
||||||
|
|
||||||
print('Finishing generating input file\n')
|
|
||||||
@ -6,20 +6,20 @@ import numpy as np
|
|||||||
|
|
||||||
#parameters need to be changed
|
#parameters need to be changed
|
||||||
#start
|
#start
|
||||||
nx=18
|
nx=75
|
||||||
ny=18
|
ny=96
|
||||||
|
nz=17
|
||||||
minvel=0.8
|
minvel=0.8
|
||||||
velgrad=0.5
|
velgrad=0.5
|
||||||
dep1=np.array([0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5])
|
dep1=1.5+np.array([-1.5, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,10.0,11.0,13.0,16.0,20.0,30.0])
|
||||||
#dep1=np.array([0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.1,1.3,1.5,1.8,2.1,2.5])
|
vel1=np.loadtxt('mod.1d')
|
||||||
nz=len(dep1)
|
|
||||||
#end
|
#end
|
||||||
vs1=np.zeros(nz)
|
vs1=np.zeros(nz)
|
||||||
mod=np.zeros((nz*ny,nx))
|
mod=np.zeros((nz*ny,nx))
|
||||||
for k in range(nz):
|
for k in range(nz):
|
||||||
for j in range(ny):
|
for j in range(ny):
|
||||||
for i in range(nx):
|
for i in range(nx):
|
||||||
mod[k*ny+j,i]= minvel+dep1[k]*velgrad
|
mod[k*ny+j,i]= vel1[k]/1.75#minvel+dep1[k]*velgrad
|
||||||
with open('MOD','w') as fp:
|
with open('MOD','w') as fp:
|
||||||
for i in range(nz):
|
for i in range(nz):
|
||||||
fp.write('%5.1f' % dep1[i])
|
fp.write('%5.1f' % dep1[i])
|
||||||
@ -30,4 +30,4 @@ with open('MOD','w') as fp:
|
|||||||
fp.write('%7.3f' % mod[k*ny+j,i])
|
fp.write('%7.3f' % mod[k*ny+j,i])
|
||||||
fp.write('\n')
|
fp.write('\n')
|
||||||
for i in range(nz):
|
for i in range(nz):
|
||||||
print (dep1[i]),
|
print dep1[i],
|
||||||
|
|||||||
@ -12,12 +12,11 @@ nz=9
|
|||||||
minvel=0.9
|
minvel=0.9
|
||||||
velgrad=0.6
|
velgrad=0.6
|
||||||
dep1=[0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5]
|
dep1=[0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5]
|
||||||
anosize=0.5
|
anosize=1.0
|
||||||
amplitude=0.4
|
amplitude=0.4
|
||||||
#end
|
#end
|
||||||
nmax = np.max([nx,ny])
|
x=range(1,nx+1)
|
||||||
x=range(1,nmax+1)
|
y=range(1,ny+1)
|
||||||
y=range(1,nmax+1)
|
|
||||||
z=range(1,nz+1)
|
z=range(1,nz+1)
|
||||||
z=np.ones(nz)
|
z=np.ones(nz)
|
||||||
bg=np.zeros((nz,nx,ny))
|
bg=np.zeros((nz,nx,ny))
|
||||||
@ -25,7 +24,7 @@ cross=np.zeros((nz,ny))
|
|||||||
vs1=np.zeros(nz)
|
vs1=np.zeros(nz)
|
||||||
xy=np.kron(np.sin(anosize*np.array(y)),np.sin(anosize*np.array(x)))
|
xy=np.kron(np.sin(anosize*np.array(y)),np.sin(anosize*np.array(x)))
|
||||||
xyz=np.kron(z,xy)
|
xyz=np.kron(z,xy)
|
||||||
pxy=xyz.reshape(nz,nmax,nmax)
|
pxy=xyz.reshape(nz,nx,ny)
|
||||||
for k in range(nz):
|
for k in range(nz):
|
||||||
for j in range(ny):
|
for j in range(ny):
|
||||||
for i in range(nx):
|
for i in range(nx):
|
||||||
@ -36,7 +35,7 @@ for k in range(nz):
|
|||||||
for i in range(nx):
|
for i in range(nx):
|
||||||
mod[(k)*ny+j,i]=bg[k,i,j]+pxy[k,i,j]*amplitude
|
mod[(k)*ny+j,i]=bg[k,i,j]+pxy[k,i,j]*amplitude
|
||||||
k=5
|
k=5
|
||||||
plt.imshow(mod[k*ny:(k+1)*ny,:],cmap='jet_r')
|
plt.imshow(mod[k*ny:(k+1)*ny,:],cmap='jet_r',interpolation='bicubic')
|
||||||
np.savetxt('MOD.true',mod,fmt='%4.4f')
|
np.savetxt('MOD.true',mod,fmt='%4.4f')
|
||||||
plt.colorbar()
|
plt.colorbar()
|
||||||
plt.show()
|
plt.show()
|
||||||
|
|||||||
@ -1,69 +0,0 @@
|
|||||||
#!/bin/csh
|
|
||||||
|
|
||||||
#-----------------------------------------------------------
|
|
||||||
# 2015-06-08 Hongjian Fang
|
|
||||||
# step 1: get velicoty along the track from each depth-layer
|
|
||||||
# step 2: combine tracks into 2D profile
|
|
||||||
#-----------------------------------------------------------
|
|
||||||
|
|
||||||
# change the following parameters according to your case
|
|
||||||
# start
|
|
||||||
set inp3D = DSurfTomo.inMeasure.dat
|
|
||||||
set pstart = 121.4/25.1 #start point for 2D profile
|
|
||||||
set pend = 121.55/25.1 #end point for 2D profile
|
|
||||||
set pstart2 = 121.5/25.0 #start point for 2D profile
|
|
||||||
set pend2 = 121.5/25.14 #end point for 2D profile
|
|
||||||
set ddint = 0.017/0.015 #grid interval
|
|
||||||
# over
|
|
||||||
set diInt = 0.5 #distance interval (km)
|
|
||||||
set J = -JX6.0i/-1.5i #size for plot
|
|
||||||
set cpt = slice.cpt
|
|
||||||
set ps = figcross.ps
|
|
||||||
gmt makecpt -Cjet -I -T0.5/2.5/0.1 > $cpt #velocity boundary
|
|
||||||
|
|
||||||
# start
|
|
||||||
set pro2d = tmp.profile.xzv
|
|
||||||
set RMAP = `gmt gmtinfo -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
|
||||||
rm -fr $pro2d
|
|
||||||
|
|
||||||
foreach layer(`awk '{print $3}' $inp3D| uniq`)
|
|
||||||
set llayer = `echo $layer | awk '{printf "%04d",$1*1000}'`
|
|
||||||
awk '{if($3==j) print $1,$2,$4}' j=$layer $inp3D |gmt surface -Gtmp.grd $RMAP -I$ddint -T0
|
|
||||||
# be careful -I: grid spacing for x/y
|
|
||||||
gmt project -C$pstart -E$pend -G$diInt -Q > track.dat
|
|
||||||
gmt grdtrack track.dat -Gtmp.grd | awk '{print $3,j,$4}' j=$layer > tmp.D$llayer
|
|
||||||
cat tmp.D$llayer >> $pro2d
|
|
||||||
end
|
|
||||||
#
|
|
||||||
set RMAP = `gmt gmtinfo -C $pro2d | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
|
||||||
gmt surface $pro2d -Gtmp.grd $RMAP -I$ddint -T0
|
|
||||||
#gmt grdfilter tmp.grd -D0 -Fg4 -Gtmpf.grd
|
|
||||||
gmt grdimage tmp.grd $J $RMAP -Bf2.5a5:'Distance (km)':/f0.5a1:'Depth (km)':SenW -C$cpt -K -Y5.2i> $ps
|
|
||||||
#gmt psscale -Cslice.cpt -Ba0.5f0.25:'Vs': -D2.5i/-0.7i/6.00/0.2h -O -K -P >> $ps
|
|
||||||
|
|
||||||
rm -fr tmp.D* tmp.grd tmp.profile.xzv track.dat
|
|
||||||
|
|
||||||
# start
|
|
||||||
set pro2d = tmp.profile.xzv
|
|
||||||
set RMAP = `gmt gmtinfo -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
|
||||||
rm -fr $pro2d
|
|
||||||
|
|
||||||
foreach layer(`awk '{print $3}' $inp3D| uniq`)
|
|
||||||
set llayer = `echo $layer | awk '{printf "%04d",$1*1000}'`
|
|
||||||
awk '{if($3==j) print $1,$2,$4}' j=$layer $inp3D |gmt surface -Gtmp.grd $RMAP -I$ddint -T0
|
|
||||||
# be careful -I: grid spacing for x/y
|
|
||||||
gmt project -C$pstart2 -E$pend2 -G$diInt -Q > track.dat
|
|
||||||
gmt grdtrack track.dat -Gtmp.grd | awk '{print $3,j,$4}' j=$layer > tmp.D$llayer
|
|
||||||
cat tmp.D$llayer >> $pro2d
|
|
||||||
end
|
|
||||||
#
|
|
||||||
set RMAP = `gmt gmtinfo -C $pro2d | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
|
||||||
gmt surface $pro2d -Gtmp.grd $RMAP -I$ddint -T0
|
|
||||||
#gmt grdfilter tmp.grd -D0 -Fg4 -Gtmpf.grd
|
|
||||||
gmt grdimage tmp.grd $J $RMAP -Bf2.5a5:'Distance (km)':/f0.5a1:'Depth (km)':SenW -C$cpt -K -O -Y-2.0i>> $ps
|
|
||||||
gmt psscale -Cslice.cpt -Ba0.5f0.25:'Vs': -D2.5i/-0.7i/6.00/0.2h -O -P >> $ps
|
|
||||||
|
|
||||||
rm -fr tmp.D* tmp.grd tmp.profile.xzv track.dat
|
|
||||||
p
|
|
||||||
ps2pdf $ps
|
|
||||||
|
|
||||||
@ -1,33 +0,0 @@
|
|||||||
#!/usr/bin/env python
|
|
||||||
|
|
||||||
import numpy as np
|
|
||||||
from matplotlib import pylab as plt
|
|
||||||
|
|
||||||
#parameters you need to change
|
|
||||||
#begin
|
|
||||||
numray = 3320 # how many rays do you have (all periods)
|
|
||||||
rb = 1 # ray index for some period (first one)
|
|
||||||
re = 50 # ray index for some period (last one)
|
|
||||||
maxseg = 500 #
|
|
||||||
#end
|
|
||||||
|
|
||||||
n = 0
|
|
||||||
numseg = np.zeros(numray,)
|
|
||||||
raylat = np.zeros((numray,maxseg))
|
|
||||||
raylon = np.zeros((numray,maxseg))
|
|
||||||
with open('raypath.out') as ray:
|
|
||||||
for line in ray:
|
|
||||||
linseg = line.split()
|
|
||||||
if linseg[0]=='#':
|
|
||||||
n = n+1
|
|
||||||
seg = 0
|
|
||||||
numseg[n-1] = float(linseg[1])
|
|
||||||
else:
|
|
||||||
seg = seg+1
|
|
||||||
raylat[n-1,seg-1] = float(linseg[0])
|
|
||||||
raylon[n-1,seg-1] = float(linseg[1])
|
|
||||||
|
|
||||||
for ii in range(rb,re):
|
|
||||||
plt.plot(raylon[ii,0:numseg[ii]],raylat[ii,0:numseg[ii]],'k-')
|
|
||||||
|
|
||||||
plt.show()
|
|
||||||
@ -7,42 +7,42 @@
|
|||||||
# csh plotslice.gmt SurfTomo.in.tvMeasure.dat 0.2 0.4 0.8 1.4
|
# csh plotslice.gmt SurfTomo.in.tvMeasure.dat 0.2 0.4 0.8 1.4
|
||||||
#-----------------------------------------------------------
|
#-----------------------------------------------------------
|
||||||
|
|
||||||
#gmt set ANOT_FONT_SIZE 6
|
gmtset ANOT_FONT_SIZE 6
|
||||||
#gmt set LABEL_FONT_SIZE 6
|
gmtset LABEL_FONT_SIZE 6
|
||||||
#gmt set FRAME_WIDTH=0.1c
|
gmtset FRAME_WIDTH=0.1c
|
||||||
#gmt set LABEL_OFFSET=0.1c
|
gmtset LABEL_OFFSET=0.1c
|
||||||
#gmt set LABEL_FONT_SIZE=5p
|
gmtset LABEL_FONT_SIZE=5p
|
||||||
#gmt set TICK_LENGTH=0.1c
|
gmtset TICK_LENGTH=0.1c
|
||||||
#gmt set TICK_PEN=0.3p
|
gmtset TICK_PEN=0.3p
|
||||||
set inp3D = $1
|
set inp3D = $1
|
||||||
set ps = horizVsnew.true.ps
|
set ps = horizVsnew.ps
|
||||||
set J = -JM2i #size for plot
|
set J = -JM2i #size for plot
|
||||||
set cpt = slice.cpt
|
set cpt = slice.cpt
|
||||||
|
|
||||||
# start
|
# start
|
||||||
gmt makecpt -Cseis -D -T0.5/2.5/0.1 > $cpt #velocity boundary
|
makecpt -Cseis -T0.6/1.5/0.1 > $cpt #velocity boundary
|
||||||
set RMAP = `gmt gmtinfo -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
set RMAP = `minmax -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
||||||
gmt psbasemap $RMAP $J -Ba0.1f0.05WseN -P -K -Y5i >$ps
|
psbasemap $RMAP $J -Ba0.1f0.05WseN -P -K -Y5i >$ps
|
||||||
awk '{if($3==depth1) print $1,$2,$4}' depth1=$2 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
awk '{if($3==depth1) print $1,$2,$4}' depth1=$2 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||||
gmt grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||||
rm -rf tmp.grd
|
rm -rf tmp.grd
|
||||||
gmt psscale -Cslice.cpt -Ba0.1f0.05:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
psscale -Cslice.cpt -Ba0.1f0.05:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
||||||
|
|
||||||
#gmt makecpt -Cseis -D -T0.8/1.8/0.1 > $cpt #velocity boundary
|
makecpt -Cseis -T0.8/1.8/0.1 > $cpt #velocity boundary
|
||||||
gmt psbasemap $RMAP $J -Ba0.1f0.05NwsE -P -K -O -X2.3i >>$ps
|
psbasemap $RMAP $J -Ba0.1f0.05NwsE -P -K -O -X2.3i >>$ps
|
||||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$3 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
awk '{if($3==depth2) print $1,$2,$4}' depth2=$3 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||||
gmt grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||||
gmt psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
||||||
|
|
||||||
#gmt makecpt -Cseis -D -T1.1/2.0/0.1 > $cpt #velocity boundary
|
makecpt -Cseis -T1.1/2.0/0.1 > $cpt #velocity boundary
|
||||||
gmt psbasemap $RMAP $J -Ba0.1f0.05WneS -P -K -O -Y-2.7i -X-2.3i >>$ps
|
psbasemap $RMAP $J -Ba0.1f0.05WneS -P -K -O -Y-2.7i -X-2.3i >>$ps
|
||||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$4 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
awk '{if($3==depth2) print $1,$2,$4}' depth2=$4 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||||
gmt grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||||
gmt psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -K -P >> $ps
|
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -K -P >> $ps
|
||||||
|
|
||||||
#gmt makecpt -Cseis -D -T1.3/2.4/0.1 > $cpt #velocity boundary
|
makecpt -Cseis -T1.3/2.4/0.1 > $cpt #velocity boundary
|
||||||
gmt psbasemap $RMAP $J -Ba0.1f0.05SwnE -P -K -O -X2.3i >>$ps
|
psbasemap $RMAP $J -Ba0.1f0.05SwnE -P -K -O -X2.3i >>$ps
|
||||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$5 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
awk '{if($3==depth2) print $1,$2,$4}' depth2=$5 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||||
gmt grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||||
gmt psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -P >> $ps
|
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -P >> $ps
|
||||||
|
|
||||||
|
|||||||
@ -92,17 +92,6 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, &
|
|||||||
dlncg_dlnvs(nn,i) = (cg2(nn)-cg1(nn))/(dlnVs*vsz(i))
|
dlncg_dlnvs(nn,i) = (cg2(nn)-cg1(nn))/(dlnVs*vsz(i))
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
!note here, no integral from z1 to z2 is needed. The grid based nodes
|
|
||||||
! have already taken into consideration of it. Layer based need integration.
|
|
||||||
! Test passed
|
|
||||||
!if (i == 1) then
|
|
||||||
! dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*thkm(1)/2.0
|
|
||||||
!else if (i == mmax) then
|
|
||||||
! dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*thkm(mmax-1)/2.0
|
|
||||||
!else
|
|
||||||
! dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0
|
|
||||||
!endif
|
|
||||||
|
|
||||||
|
|
||||||
vpm(i) = vpz(i) - 0.5*dlnVp*vpz(i)
|
vpm(i) = vpz(i) - 0.5*dlnVp*vpz(i)
|
||||||
call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,&
|
call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,&
|
||||||
@ -121,16 +110,6 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, &
|
|||||||
! cga = 0.5*(cg1(nn)+cg2(nn))
|
! cga = 0.5*(cg1(nn)+cg2(nn))
|
||||||
dlncg_dlnvp(nn,i) = (cg2(nn)-cg1(nn))/(dlnVp*vpz(i))
|
dlncg_dlnvp(nn,i) = (cg2(nn)-cg1(nn))/(dlnVp*vpz(i))
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
|
|
||||||
!if (i == 1) then
|
|
||||||
! dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*thkm(1)/2.0
|
|
||||||
!else if (i == mmax) then
|
|
||||||
! dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*thkm(mmax-1)/2.0
|
|
||||||
!else
|
|
||||||
! dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0
|
|
||||||
!endif
|
|
||||||
|
|
||||||
rhom(i) = rhoz(i) - 0.5*dlnrho*rhoz(i)
|
rhom(i) = rhoz(i) - 0.5*dlnrho*rhoz(i)
|
||||||
call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,&
|
call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,&
|
||||||
rmax,rdep,rvp,rvs,rrho,rthk)
|
rmax,rdep,rvp,rvs,rrho,rthk)
|
||||||
@ -148,15 +127,6 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, &
|
|||||||
! cga = 0.5*(cg1(nn)+cg2(nn))
|
! cga = 0.5*(cg1(nn)+cg2(nn))
|
||||||
dlncg_dlnrho(nn,i) = (cg2(nn)-cg1(nn))/(dlnrho*rhoz(i))
|
dlncg_dlnrho(nn,i) = (cg2(nn)-cg1(nn))/(dlnrho*rhoz(i))
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
! if (i == 1) then
|
|
||||||
! dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*thkm(1)/2.0
|
|
||||||
! else if (i == mmax) then
|
|
||||||
! dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*thkm(mmax-1)/2.0
|
|
||||||
! else
|
|
||||||
! dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0
|
|
||||||
! endif
|
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
sen_vsRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnvs(1:kmaxRc,1:mmax)
|
sen_vsRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnvs(1:kmaxRc,1:mmax)
|
||||||
sen_vpRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnvp(1:kmaxRc,1:mmax)
|
sen_vpRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnvp(1:kmaxRc,1:mmax)
|
||||||
@ -940,7 +910,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
|||||||
goxdf,gozdf,dvxdf,dvzdf,kmaxRc,kmaxRg,kmaxLc,kmaxLg, &
|
goxdf,gozdf,dvxdf,dvzdf,kmaxRc,kmaxRg,kmaxLc,kmaxLg, &
|
||||||
tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk, &
|
tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk, &
|
||||||
scxf,sczf,rcxf,rczf,nrc1,nsrcsurf1,kmax,nsrcsurf,nrcf, &
|
scxf,sczf,rcxf,rczf,nrc1,nsrcsurf1,kmax,nsrcsurf,nrcf, &
|
||||||
nar)
|
nar,writepath)
|
||||||
USE globalp
|
USE globalp
|
||||||
USE traveltime
|
USE traveltime
|
||||||
IMPLICIT NONE
|
IMPLICIT NONE
|
||||||
@ -1027,10 +997,11 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
|||||||
integer ii,jj,kk,nn,istep
|
integer ii,jj,kk,nn,istep
|
||||||
integer level,maxlevel,maxleveld,HorizonType,VerticalType,PorS
|
integer level,maxlevel,maxleveld,HorizonType,VerticalType,PorS
|
||||||
real,parameter::ftol=1e-4
|
real,parameter::ftol=1e-4
|
||||||
|
integer writepath
|
||||||
integer ig, igroup
|
integer ig, igroup
|
||||||
|
|
||||||
gdx=8
|
gdx=5
|
||||||
gdz=8
|
gdz=5
|
||||||
asgr=1
|
asgr=1
|
||||||
sgdl=8
|
sgdl=8
|
||||||
sgs=8
|
sgs=8
|
||||||
@ -1098,7 +1069,6 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
|||||||
if(kmaxRc.gt.0) then
|
if(kmaxRc.gt.0) then
|
||||||
iwave=2
|
iwave=2
|
||||||
igr=0
|
igr=0
|
||||||
print*,'Rayleigh wave phase velocity depth kernel'
|
|
||||||
call depthkernel(nx,ny,nz,vels,pvRc,sen_vsRc,sen_vpRc, &
|
call depthkernel(nx,ny,nz,vels,pvRc,sen_vsRc,sen_vpRc, &
|
||||||
sen_rhoRc,iwave,igr,kmaxRc,tRc,depz,minthk)
|
sen_rhoRc,iwave,igr,kmaxRc,tRc,depz,minthk)
|
||||||
endif
|
endif
|
||||||
@ -1108,9 +1078,8 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
|||||||
igr=0
|
igr=0
|
||||||
! print*,kmax
|
! print*,kmax
|
||||||
call caldespersion(nx,ny,nz,vels,pvRc, &
|
call caldespersion(nx,ny,nz,vels,pvRc, &
|
||||||
iwave,igr,kmaxRg,tRg,depz,minthk)
|
iwave,igr,kmax,tRg,depz,minthk)
|
||||||
igr=1
|
igr=1
|
||||||
print*,'Rayleigh wave group velocity depth kernel'
|
|
||||||
call depthkernel(nx,ny,nz,vels,pvRg,sen_vsRg,sen_vpRg, &
|
call depthkernel(nx,ny,nz,vels,pvRg,sen_vsRg,sen_vpRg, &
|
||||||
sen_rhoRg,iwave,igr,kmaxRg,tRg,depz,minthk)
|
sen_rhoRg,iwave,igr,kmaxRg,tRg,depz,minthk)
|
||||||
endif
|
endif
|
||||||
@ -1119,14 +1088,14 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
|||||||
iwave=1
|
iwave=1
|
||||||
igr=0
|
igr=0
|
||||||
call depthkernel(nx,ny,nz,vels,pvLc,sen_vsLc,sen_vpLc, &
|
call depthkernel(nx,ny,nz,vels,pvLc,sen_vsLc,sen_vpLc, &
|
||||||
sen_rhoLc,iwave,igr,kmaxLc,tLc,depz,minthk)
|
sen_rhoLc,iwave,igr,kmax,tLc,depz,minthk)
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if(kmaxLg.gt.0) then
|
if(kmaxLg.gt.0) then
|
||||||
iwave=1
|
iwave=1
|
||||||
igr=0
|
igr=0
|
||||||
call caldespersion(nx,ny,nz,vels,pvLc, &
|
call caldespersion(nx,ny,nz,vels,pvLc, &
|
||||||
iwave,igr,kmaxLg,tLg,depz,minthk)
|
iwave,igr,kmax,tLg,depz,minthk)
|
||||||
igr=1
|
igr=1
|
||||||
call depthkernel(nx,ny,nz,vels,pvLg,sen_vsLg,sen_vpLg, &
|
call depthkernel(nx,ny,nz,vels,pvLg,sen_vsLg,sen_vpLg, &
|
||||||
sen_rhoLg,iwave,igr,kmaxLg,tLg,depz,minthk)
|
sen_rhoLg,iwave,igr,kmaxLg,tLg,depz,minthk)
|
||||||
@ -1283,7 +1252,6 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
|||||||
ALLOCATE(ttnr(idm2,idm1))
|
ALLOCATE(ttnr(idm2,idm1))
|
||||||
ALLOCATE(nstsr(idm2,idm1))
|
ALLOCATE(nstsr(idm2,idm1))
|
||||||
ENDIF
|
ENDIF
|
||||||
!ttnr(1:nnzb,1:nnxb)=ttn(1:nnzb,1:nnxb)
|
|
||||||
ttnr=ttn
|
ttnr=ttn
|
||||||
nstsr=nsts
|
nstsr=nsts
|
||||||
ogx=vnl
|
ogx=vnl
|
||||||
@ -1378,11 +1346,13 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
|||||||
! if required.
|
! if required.
|
||||||
!
|
!
|
||||||
if (igrt(srcnum,knumi) == 0 .or. (ig == 2 .and. igrt(srcnum,knumi) == 1)) then
|
if (igrt(srcnum,knumi) == 0 .or. (ig == 2 .and. igrt(srcnum,knumi) == 1)) then
|
||||||
|
! a little stupid, remember to change latter
|
||||||
|
if (igrt(srcnum,knumi) == 1) then
|
||||||
|
call gridder(velf0)
|
||||||
|
endif
|
||||||
count11=count11+1
|
count11=count11+1
|
||||||
CALL rpaths(x,z,fdm,rcxf(istep,srcnum,knumi),rczf(istep,srcnum,knumi))
|
CALL rpaths(x,z,fdm,rcxf(istep,srcnum,knumi),rczf(istep,srcnum,knumi),writepath)
|
||||||
row(1:nparpi)=0.0
|
row(1:nparpi)=0.0
|
||||||
! change the Brocher relationship if depth is larger than 35 km
|
|
||||||
if (depz(nz-1)<35.0) then
|
|
||||||
do jj=1,nvz
|
do jj=1,nvz
|
||||||
do kk=1,nvx
|
do kk=1,nvx
|
||||||
if(abs(fdm(jj,kk)).ge.ftol) then
|
if(abs(fdm(jj,kk)).ge.ftol) then
|
||||||
@ -1400,28 +1370,6 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
|||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
! use a different formulation between Vp, Vs, and Rho, according to Luosong's
|
|
||||||
else
|
|
||||||
do jj=1,nvz
|
|
||||||
do kk=1,nvx
|
|
||||||
if(abs(fdm(jj,kk)).ge.ftol) then
|
|
||||||
coe_a=(2.2110-0.8984*2*vels(kk+1,jj+1,1:nz-1)+&
|
|
||||||
0.2786*3*vels(kk+1,jj+1,1:nz-1)**2-0.02412*4*vels(kk+1,jj+1,1:nz-1)**3)
|
|
||||||
vpft=0.9098 + 2.2110*vels(kk+1,jj+1,1:nz-1) - 0.8984*vels(kk+1,jj+1,1:nz-1)**2+ &
|
|
||||||
0.2786*vels(kk+1,jj+1,1:nz-1)**3 - 0.02412*vels(kk+1,jj+1,1:nz-1)**4
|
|
||||||
coe_rho=coe_a*(1.6612-0.4721*2*vpft+&
|
|
||||||
0.0671*3*vpft**2-0.0043*4*vpft**3+&
|
|
||||||
0.000106*5*vpft**4)
|
|
||||||
|
|
||||||
row((jj-1)*nvx+kk:(nz-2)*nvz*nvx+(jj-1)*nvx+kk:nvx*nvz)=&
|
|
||||||
(sen_vp(jj*(nvx+2)+kk+1,knumi,1:nz-1)*coe_a+&
|
|
||||||
sen_rho(jj*(nvx+2)+kk+1,knumi,1:nz-1)*coe_rho+&
|
|
||||||
sen_vs(jj*(nvx+2)+kk+1,knumi,1:nz-1))*fdm(jj,kk)
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
endif
|
|
||||||
|
|
||||||
do nn=1,nparpi
|
do nn=1,nparpi
|
||||||
if(abs(row(nn)).gt.ftol) then
|
if(abs(row(nn)).gt.ftol) then
|
||||||
nar=nar+1
|
nar=nar+1
|
||||||
@ -1768,7 +1716,7 @@ END SUBROUTINE srtimes
|
|||||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
!SUBROUTINE rpaths(wrgf,csid,cfd,scx,scz)
|
!SUBROUTINE rpaths(wrgf,csid,cfd,scx,scz)
|
||||||
!SUBROUTINE rpaths()
|
!SUBROUTINE rpaths()
|
||||||
SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz)
|
SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz,writepath)
|
||||||
USE globalp
|
USE globalp
|
||||||
IMPLICIT NONE
|
IMPLICIT NONE
|
||||||
INTEGER, PARAMETER :: i5=SELECTED_REAL_KIND(5,10)
|
INTEGER, PARAMETER :: i5=SELECTED_REAL_KIND(5,10)
|
||||||
@ -1789,6 +1737,7 @@ SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz)
|
|||||||
!fang!------------------------------------------------
|
!fang!------------------------------------------------
|
||||||
real fdm(0:nvz+1,0:nvx+1)
|
real fdm(0:nvz+1,0:nvx+1)
|
||||||
REAL(KIND=i10) surfrcx,surfrcz
|
REAL(KIND=i10) surfrcx,surfrcz
|
||||||
|
integer writepath
|
||||||
!fang!------------------------------------------------
|
!fang!------------------------------------------------
|
||||||
!
|
!
|
||||||
! ipx,ipz = Coordinates of cell containing current point
|
! ipx,ipz = Coordinates of cell containing current point
|
||||||
@ -2273,14 +2222,14 @@ SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz)
|
|||||||
! Write ray paths to output file
|
! Write ray paths to output file
|
||||||
!
|
!
|
||||||
!fang! IF(wrgf.EQ.csid.OR.wrgf.LT.0)THEN
|
!fang! IF(wrgf.EQ.csid.OR.wrgf.LT.0)THEN
|
||||||
!if(writepath == 1) then
|
if(writepath == 1) then
|
||||||
! WRITE(40,*)'#',nrp
|
WRITE(40,*)'#',nrp
|
||||||
! DO j=1,nrp
|
DO j=1,nrp
|
||||||
! rayx=(pi/2-rgx(j))*180.0/pi
|
rayx=(pi/2-rgx(j))*180.0/pi
|
||||||
! rayz=rgz(j)*180.0/pi
|
rayz=rgz(j)*180.0/pi
|
||||||
! WRITE(40,*)rayx,rayz
|
WRITE(40,*)rayx,rayz
|
||||||
! ENDDO
|
ENDDO
|
||||||
!endif
|
endif
|
||||||
!fang! ENDIF
|
!fang! ENDIF
|
||||||
!
|
!
|
||||||
! Write partial derivatives to output file
|
! Write partial derivatives to output file
|
||||||
|
|||||||
BIN
src/DSurfTomo
Executable file
BIN
src/DSurfTomo
Executable file
Binary file not shown.
12
src/Makefile
12
src/Makefile
@ -1,19 +1,17 @@
|
|||||||
CMD = DSurfTomo
|
CMD = DSurfTomo
|
||||||
FC = gfortran
|
FC = gfortran
|
||||||
FFLAGS = -O -ffixed-line-length-none -ffloat-store\
|
FFLAGS = -O3 -ffixed-line-length-none -ffloat-store\
|
||||||
-fbounds-check -m64 -mcmodel=medium
|
-W -fbounds-check -m64 -mcmodel=medium
|
||||||
F90SRCS = lsmrDataModule.f90 lsmrblasInterface.f90\
|
F90SRCS = lsmrDataModule.f90 lsmrblasInterface.f90\
|
||||||
lsmrblas.f90 lsmrModule.f90 delsph.f90\
|
lsmrblas.f90 lsmrModule.f90 delsph.f90\
|
||||||
aprod.f90 gaussian.f90 getpercentile.f90
|
aprod.f90 gaussian.f90 main.f90
|
||||||
FSRCS = surfdisp96.f slarnv.f slaruv.f
|
FSRCS = surfdisp96.f
|
||||||
OBJS = $(F90SRCS:%.f90=%.o) $(FSRCS:%.f=%.o) CalSurfG.o main.o
|
OBJS = $(F90SRCS:%.f90=%.o) $(FSRCS:%.f=%.o) CalSurfG.o
|
||||||
all:$(CMD)
|
all:$(CMD)
|
||||||
$(CMD):$(OBJS)
|
$(CMD):$(OBJS)
|
||||||
$(FC) -fopenmp $^ -o $@
|
$(FC) -fopenmp $^ -o $@
|
||||||
CalSurfG.o:CalSurfG.f90
|
CalSurfG.o:CalSurfG.f90
|
||||||
$(FC) -fopenmp $(FFLAGS) -c $< -o $@
|
$(FC) -fopenmp $(FFLAGS) -c $< -o $@
|
||||||
main.o:main.f90
|
|
||||||
$(FC) -fopenmp $(FFLAGS) -c $< -o $@
|
|
||||||
%.o: %.f90
|
%.o: %.f90
|
||||||
$(FC) $(FFLAGS) -c $(@F:.o=.f90) -o $@
|
$(FC) $(FFLAGS) -c $(@F:.o=.f90) -o $@
|
||||||
%.o: %.f
|
%.o: %.f
|
||||||
|
|||||||
@ -1,51 +0,0 @@
|
|||||||
SUBROUTINE getpercentile(N,array,q25,q75)
|
|
||||||
real q25
|
|
||||||
real q75
|
|
||||||
integer N
|
|
||||||
real,intent(in) :: array(*)
|
|
||||||
integer idx
|
|
||||||
|
|
||||||
real RA(N)
|
|
||||||
RA(1:N) = array(1:N)
|
|
||||||
L=N/2+1
|
|
||||||
IR=N
|
|
||||||
!The index L will be decremented from its initial value during the
|
|
||||||
!"hiring" (heap creation) phase. Once it reaches 1, the index IR
|
|
||||||
!will be decremented from its initial value down to 1 during the
|
|
||||||
!"retirement-and-promotion" (heap selection) phase.
|
|
||||||
10 continue
|
|
||||||
if(L > 1)then
|
|
||||||
L=L-1
|
|
||||||
RRA=RA(L)
|
|
||||||
else
|
|
||||||
RRA=RA(IR)
|
|
||||||
RA(IR)=RA(1)
|
|
||||||
IR=IR-1
|
|
||||||
if(IR.eq.1)then
|
|
||||||
RA(1)=RRA
|
|
||||||
idx = int(0.25*N)
|
|
||||||
q25 = RA(idx)
|
|
||||||
idx = int(0.75*N)
|
|
||||||
q75 = RA(idx)
|
|
||||||
return
|
|
||||||
end if
|
|
||||||
end if
|
|
||||||
I=L
|
|
||||||
J=L+L
|
|
||||||
20 if(J.le.IR)then
|
|
||||||
if(J < IR)then
|
|
||||||
if(RA(J) < RA(J+1)) J=J+1
|
|
||||||
end if
|
|
||||||
if(RRA < RA(J))then
|
|
||||||
RA(I)=RA(J)
|
|
||||||
I=J; J=J+J
|
|
||||||
else
|
|
||||||
J=IR+1
|
|
||||||
end if
|
|
||||||
goto 20
|
|
||||||
end if
|
|
||||||
RA(I)=RRA
|
|
||||||
goto 10
|
|
||||||
|
|
||||||
END
|
|
||||||
|
|
||||||
1233
src/main.f90
1233
src/main.f90
File diff suppressed because it is too large
Load Diff
178
src/slarnv.f
178
src/slarnv.f
@ -1,178 +0,0 @@
|
|||||||
*> \brief \b SLARNV returns a vector of random numbers from a uniform or normal distribution.
|
|
||||||
*
|
|
||||||
* =========== DOCUMENTATION ===========
|
|
||||||
*
|
|
||||||
* Online html documentation available at
|
|
||||||
* http://www.netlib.org/lapack/explore-html/
|
|
||||||
*
|
|
||||||
*> \htmlonly
|
|
||||||
*> Download SLARNV + dependencies
|
|
||||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarnv.f">
|
|
||||||
*> [TGZ]</a>
|
|
||||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarnv.f">
|
|
||||||
*> [ZIP]</a>
|
|
||||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarnv.f">
|
|
||||||
*> [TXT]</a>
|
|
||||||
*> \endhtmlonly
|
|
||||||
*
|
|
||||||
* Definition:
|
|
||||||
* ===========
|
|
||||||
*
|
|
||||||
* SUBROUTINE SLARNV( IDIST, ISEED, N, X )
|
|
||||||
*
|
|
||||||
* .. Scalar Arguments ..
|
|
||||||
* INTEGER IDIST, N
|
|
||||||
* ..
|
|
||||||
* .. Array Arguments ..
|
|
||||||
* INTEGER ISEED( 4 )
|
|
||||||
* REAL X( * )
|
|
||||||
* ..
|
|
||||||
*
|
|
||||||
*
|
|
||||||
*> \par Purpose:
|
|
||||||
* =============
|
|
||||||
*>
|
|
||||||
*> \verbatim
|
|
||||||
*>
|
|
||||||
*> SLARNV returns a vector of n random real numbers from a uniform or
|
|
||||||
*> normal distribution.
|
|
||||||
*> \endverbatim
|
|
||||||
*
|
|
||||||
* Arguments:
|
|
||||||
* ==========
|
|
||||||
*
|
|
||||||
*> \param[in] IDIST
|
|
||||||
*> \verbatim
|
|
||||||
*> IDIST is INTEGER
|
|
||||||
*> Specifies the distribution of the random numbers:
|
|
||||||
*> = 1: uniform (0,1)
|
|
||||||
*> = 2: uniform (-1,1)
|
|
||||||
*> = 3: normal (0,1)
|
|
||||||
*> \endverbatim
|
|
||||||
*>
|
|
||||||
*> \param[in,out] ISEED
|
|
||||||
*> \verbatim
|
|
||||||
*> ISEED is INTEGER array, dimension (4)
|
|
||||||
*> On entry, the seed of the random number generator; the array
|
|
||||||
*> elements must be between 0 and 4095, and ISEED(4) must be
|
|
||||||
*> odd.
|
|
||||||
*> On exit, the seed is updated.
|
|
||||||
*> \endverbatim
|
|
||||||
*>
|
|
||||||
*> \param[in] N
|
|
||||||
*> \verbatim
|
|
||||||
*> N is INTEGER
|
|
||||||
*> The number of random numbers to be generated.
|
|
||||||
*> \endverbatim
|
|
||||||
*>
|
|
||||||
*> \param[out] X
|
|
||||||
*> \verbatim
|
|
||||||
*> X is REAL array, dimension (N)
|
|
||||||
*> The generated random numbers.
|
|
||||||
*> \endverbatim
|
|
||||||
*
|
|
||||||
* Authors:
|
|
||||||
* ========
|
|
||||||
*
|
|
||||||
*> \author Univ. of Tennessee
|
|
||||||
*> \author Univ. of California Berkeley
|
|
||||||
*> \author Univ. of Colorado Denver
|
|
||||||
*> \author NAG Ltd.
|
|
||||||
*
|
|
||||||
*> \date December 2016
|
|
||||||
*
|
|
||||||
*> \ingroup OTHERauxiliary
|
|
||||||
*
|
|
||||||
*> \par Further Details:
|
|
||||||
* =====================
|
|
||||||
*>
|
|
||||||
*> \verbatim
|
|
||||||
*>
|
|
||||||
*> This routine calls the auxiliary routine SLARUV to generate random
|
|
||||||
*> real numbers from a uniform (0,1) distribution, in batches of up to
|
|
||||||
*> 128 using vectorisable code. The Box-Muller method is used to
|
|
||||||
*> transform numbers from a uniform to a normal distribution.
|
|
||||||
*> \endverbatim
|
|
||||||
*>
|
|
||||||
* =====================================================================
|
|
||||||
SUBROUTINE SLARNV( IDIST, ISEED, N, X )
|
|
||||||
*
|
|
||||||
* -- LAPACK auxiliary routine (version 3.7.0) --
|
|
||||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
|
||||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
|
||||||
* December 2016
|
|
||||||
*
|
|
||||||
* .. Scalar Arguments ..
|
|
||||||
INTEGER IDIST, N
|
|
||||||
* ..
|
|
||||||
* .. Array Arguments ..
|
|
||||||
INTEGER ISEED( 4 )
|
|
||||||
REAL X( * )
|
|
||||||
* ..
|
|
||||||
*
|
|
||||||
* =====================================================================
|
|
||||||
*
|
|
||||||
* .. Parameters ..
|
|
||||||
REAL ONE, TWO
|
|
||||||
PARAMETER ( ONE = 1.0E+0, TWO = 2.0E+0 )
|
|
||||||
INTEGER LV
|
|
||||||
PARAMETER ( LV = 128 )
|
|
||||||
REAL TWOPI
|
|
||||||
PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
|
|
||||||
* ..
|
|
||||||
* .. Local Scalars ..
|
|
||||||
INTEGER I, IL, IL2, IV
|
|
||||||
* ..
|
|
||||||
* .. Local Arrays ..
|
|
||||||
REAL U( LV )
|
|
||||||
* ..
|
|
||||||
* .. Intrinsic Functions ..
|
|
||||||
INTRINSIC COS, LOG, MIN, SQRT
|
|
||||||
* ..
|
|
||||||
* .. External Subroutines ..
|
|
||||||
EXTERNAL SLARUV
|
|
||||||
* ..
|
|
||||||
* .. Executable Statements ..
|
|
||||||
*
|
|
||||||
DO 40 IV = 1, N, LV / 2
|
|
||||||
IL = MIN( LV / 2, N-IV+1 )
|
|
||||||
IF( IDIST.EQ.3 ) THEN
|
|
||||||
IL2 = 2*IL
|
|
||||||
ELSE
|
|
||||||
IL2 = IL
|
|
||||||
END IF
|
|
||||||
*
|
|
||||||
* Call SLARUV to generate IL2 numbers from a uniform (0,1)
|
|
||||||
* distribution (IL2 <= LV)
|
|
||||||
*
|
|
||||||
CALL SLARUV( ISEED, IL2, U )
|
|
||||||
*
|
|
||||||
IF( IDIST.EQ.1 ) THEN
|
|
||||||
*
|
|
||||||
* Copy generated numbers
|
|
||||||
*
|
|
||||||
DO 10 I = 1, IL
|
|
||||||
X( IV+I-1 ) = U( I )
|
|
||||||
10 CONTINUE
|
|
||||||
ELSE IF( IDIST.EQ.2 ) THEN
|
|
||||||
*
|
|
||||||
* Convert generated numbers to uniform (-1,1) distribution
|
|
||||||
*
|
|
||||||
DO 20 I = 1, IL
|
|
||||||
X( IV+I-1 ) = TWO*U( I ) - ONE
|
|
||||||
20 CONTINUE
|
|
||||||
ELSE IF( IDIST.EQ.3 ) THEN
|
|
||||||
*
|
|
||||||
* Convert generated numbers to normal (0,1) distribution
|
|
||||||
*
|
|
||||||
DO 30 I = 1, IL
|
|
||||||
X( IV+I-1 ) = SQRT( -TWO*LOG( U( 2*I-1 ) ) )*
|
|
||||||
$ COS( TWOPI*U( 2*I ) )
|
|
||||||
30 CONTINUE
|
|
||||||
END IF
|
|
||||||
40 CONTINUE
|
|
||||||
RETURN
|
|
||||||
*
|
|
||||||
* End of SLARNV
|
|
||||||
*
|
|
||||||
END
|
|
||||||
447
src/slaruv.f
447
src/slaruv.f
@ -1,447 +0,0 @@
|
|||||||
*> \brief \b SLARUV returns a vector of n random real numbers from a uniform distribution.
|
|
||||||
*
|
|
||||||
* =========== DOCUMENTATION ===========
|
|
||||||
*
|
|
||||||
* Online html documentation available at
|
|
||||||
* http://www.netlib.org/lapack/explore-html/
|
|
||||||
*
|
|
||||||
*> \htmlonly
|
|
||||||
*> Download SLARUV + dependencies
|
|
||||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaruv.f">
|
|
||||||
*> [TGZ]</a>
|
|
||||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaruv.f">
|
|
||||||
*> [ZIP]</a>
|
|
||||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaruv.f">
|
|
||||||
*> [TXT]</a>
|
|
||||||
*> \endhtmlonly
|
|
||||||
*
|
|
||||||
* Definition:
|
|
||||||
* ===========
|
|
||||||
*
|
|
||||||
* SUBROUTINE SLARUV( ISEED, N, X )
|
|
||||||
*
|
|
||||||
* .. Scalar Arguments ..
|
|
||||||
* INTEGER N
|
|
||||||
* ..
|
|
||||||
* .. Array Arguments ..
|
|
||||||
* INTEGER ISEED( 4 )
|
|
||||||
* REAL X( N )
|
|
||||||
* ..
|
|
||||||
*
|
|
||||||
*
|
|
||||||
*> \par Purpose:
|
|
||||||
* =============
|
|
||||||
*>
|
|
||||||
*> \verbatim
|
|
||||||
*>
|
|
||||||
*> SLARUV returns a vector of n random real numbers from a uniform (0,1)
|
|
||||||
*> distribution (n <= 128).
|
|
||||||
*>
|
|
||||||
*> This is an auxiliary routine called by SLARNV and CLARNV.
|
|
||||||
*> \endverbatim
|
|
||||||
*
|
|
||||||
* Arguments:
|
|
||||||
* ==========
|
|
||||||
*
|
|
||||||
*> \param[in,out] ISEED
|
|
||||||
*> \verbatim
|
|
||||||
*> ISEED is INTEGER array, dimension (4)
|
|
||||||
*> On entry, the seed of the random number generator; the array
|
|
||||||
*> elements must be between 0 and 4095, and ISEED(4) must be
|
|
||||||
*> odd.
|
|
||||||
*> On exit, the seed is updated.
|
|
||||||
*> \endverbatim
|
|
||||||
*>
|
|
||||||
*> \param[in] N
|
|
||||||
*> \verbatim
|
|
||||||
*> N is INTEGER
|
|
||||||
*> The number of random numbers to be generated. N <= 128.
|
|
||||||
*> \endverbatim
|
|
||||||
*>
|
|
||||||
*> \param[out] X
|
|
||||||
*> \verbatim
|
|
||||||
*> X is REAL array, dimension (N)
|
|
||||||
*> The generated random numbers.
|
|
||||||
*> \endverbatim
|
|
||||||
*
|
|
||||||
* Authors:
|
|
||||||
* ========
|
|
||||||
*
|
|
||||||
*> \author Univ. of Tennessee
|
|
||||||
*> \author Univ. of California Berkeley
|
|
||||||
*> \author Univ. of Colorado Denver
|
|
||||||
*> \author NAG Ltd.
|
|
||||||
*
|
|
||||||
*> \date December 2016
|
|
||||||
*
|
|
||||||
*> \ingroup OTHERauxiliary
|
|
||||||
*
|
|
||||||
*> \par Further Details:
|
|
||||||
* =====================
|
|
||||||
*>
|
|
||||||
*> \verbatim
|
|
||||||
*>
|
|
||||||
*> This routine uses a multiplicative congruential method with modulus
|
|
||||||
*> 2**48 and multiplier 33952834046453 (see G.S.Fishman,
|
|
||||||
*> 'Multiplicative congruential random number generators with modulus
|
|
||||||
*> 2**b: an exhaustive analysis for b = 32 and a partial analysis for
|
|
||||||
*> b = 48', Math. Comp. 189, pp 331-344, 1990).
|
|
||||||
*>
|
|
||||||
*> 48-bit integers are stored in 4 integer array elements with 12 bits
|
|
||||||
*> per element. Hence the routine is portable across machines with
|
|
||||||
*> integers of 32 bits or more.
|
|
||||||
*> \endverbatim
|
|
||||||
*>
|
|
||||||
* =====================================================================
|
|
||||||
SUBROUTINE SLARUV( ISEED, N, X )
|
|
||||||
*
|
|
||||||
* -- LAPACK auxiliary routine (version 3.7.0) --
|
|
||||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
|
||||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
|
||||||
* December 2016
|
|
||||||
*
|
|
||||||
* .. Scalar Arguments ..
|
|
||||||
INTEGER N
|
|
||||||
* ..
|
|
||||||
* .. Array Arguments ..
|
|
||||||
INTEGER ISEED( 4 )
|
|
||||||
REAL X( N )
|
|
||||||
* ..
|
|
||||||
*
|
|
||||||
* =====================================================================
|
|
||||||
*
|
|
||||||
* .. Parameters ..
|
|
||||||
REAL ONE
|
|
||||||
PARAMETER ( ONE = 1.0E0 )
|
|
||||||
INTEGER LV, IPW2
|
|
||||||
REAL R
|
|
||||||
PARAMETER ( LV = 128, IPW2 = 4096, R = ONE / IPW2 )
|
|
||||||
* ..
|
|
||||||
* .. Local Scalars ..
|
|
||||||
INTEGER I, I1, I2, I3, I4, IT1, IT2, IT3, IT4, J
|
|
||||||
* ..
|
|
||||||
* .. Local Arrays ..
|
|
||||||
INTEGER MM( LV, 4 )
|
|
||||||
* ..
|
|
||||||
* .. Intrinsic Functions ..
|
|
||||||
INTRINSIC MIN, MOD, REAL
|
|
||||||
* ..
|
|
||||||
* .. Data statements ..
|
|
||||||
DATA ( MM( 1, J ), J = 1, 4 ) / 494, 322, 2508,
|
|
||||||
$ 2549 /
|
|
||||||
DATA ( MM( 2, J ), J = 1, 4 ) / 2637, 789, 3754,
|
|
||||||
$ 1145 /
|
|
||||||
DATA ( MM( 3, J ), J = 1, 4 ) / 255, 1440, 1766,
|
|
||||||
$ 2253 /
|
|
||||||
DATA ( MM( 4, J ), J = 1, 4 ) / 2008, 752, 3572,
|
|
||||||
$ 305 /
|
|
||||||
DATA ( MM( 5, J ), J = 1, 4 ) / 1253, 2859, 2893,
|
|
||||||
$ 3301 /
|
|
||||||
DATA ( MM( 6, J ), J = 1, 4 ) / 3344, 123, 307,
|
|
||||||
$ 1065 /
|
|
||||||
DATA ( MM( 7, J ), J = 1, 4 ) / 4084, 1848, 1297,
|
|
||||||
$ 3133 /
|
|
||||||
DATA ( MM( 8, J ), J = 1, 4 ) / 1739, 643, 3966,
|
|
||||||
$ 2913 /
|
|
||||||
DATA ( MM( 9, J ), J = 1, 4 ) / 3143, 2405, 758,
|
|
||||||
$ 3285 /
|
|
||||||
DATA ( MM( 10, J ), J = 1, 4 ) / 3468, 2638, 2598,
|
|
||||||
$ 1241 /
|
|
||||||
DATA ( MM( 11, J ), J = 1, 4 ) / 688, 2344, 3406,
|
|
||||||
$ 1197 /
|
|
||||||
DATA ( MM( 12, J ), J = 1, 4 ) / 1657, 46, 2922,
|
|
||||||
$ 3729 /
|
|
||||||
DATA ( MM( 13, J ), J = 1, 4 ) / 1238, 3814, 1038,
|
|
||||||
$ 2501 /
|
|
||||||
DATA ( MM( 14, J ), J = 1, 4 ) / 3166, 913, 2934,
|
|
||||||
$ 1673 /
|
|
||||||
DATA ( MM( 15, J ), J = 1, 4 ) / 1292, 3649, 2091,
|
|
||||||
$ 541 /
|
|
||||||
DATA ( MM( 16, J ), J = 1, 4 ) / 3422, 339, 2451,
|
|
||||||
$ 2753 /
|
|
||||||
DATA ( MM( 17, J ), J = 1, 4 ) / 1270, 3808, 1580,
|
|
||||||
$ 949 /
|
|
||||||
DATA ( MM( 18, J ), J = 1, 4 ) / 2016, 822, 1958,
|
|
||||||
$ 2361 /
|
|
||||||
DATA ( MM( 19, J ), J = 1, 4 ) / 154, 2832, 2055,
|
|
||||||
$ 1165 /
|
|
||||||
DATA ( MM( 20, J ), J = 1, 4 ) / 2862, 3078, 1507,
|
|
||||||
$ 4081 /
|
|
||||||
DATA ( MM( 21, J ), J = 1, 4 ) / 697, 3633, 1078,
|
|
||||||
$ 2725 /
|
|
||||||
DATA ( MM( 22, J ), J = 1, 4 ) / 1706, 2970, 3273,
|
|
||||||
$ 3305 /
|
|
||||||
DATA ( MM( 23, J ), J = 1, 4 ) / 491, 637, 17,
|
|
||||||
$ 3069 /
|
|
||||||
DATA ( MM( 24, J ), J = 1, 4 ) / 931, 2249, 854,
|
|
||||||
$ 3617 /
|
|
||||||
DATA ( MM( 25, J ), J = 1, 4 ) / 1444, 2081, 2916,
|
|
||||||
$ 3733 /
|
|
||||||
DATA ( MM( 26, J ), J = 1, 4 ) / 444, 4019, 3971,
|
|
||||||
$ 409 /
|
|
||||||
DATA ( MM( 27, J ), J = 1, 4 ) / 3577, 1478, 2889,
|
|
||||||
$ 2157 /
|
|
||||||
DATA ( MM( 28, J ), J = 1, 4 ) / 3944, 242, 3831,
|
|
||||||
$ 1361 /
|
|
||||||
DATA ( MM( 29, J ), J = 1, 4 ) / 2184, 481, 2621,
|
|
||||||
$ 3973 /
|
|
||||||
DATA ( MM( 30, J ), J = 1, 4 ) / 1661, 2075, 1541,
|
|
||||||
$ 1865 /
|
|
||||||
DATA ( MM( 31, J ), J = 1, 4 ) / 3482, 4058, 893,
|
|
||||||
$ 2525 /
|
|
||||||
DATA ( MM( 32, J ), J = 1, 4 ) / 657, 622, 736,
|
|
||||||
$ 1409 /
|
|
||||||
DATA ( MM( 33, J ), J = 1, 4 ) / 3023, 3376, 3992,
|
|
||||||
$ 3445 /
|
|
||||||
DATA ( MM( 34, J ), J = 1, 4 ) / 3618, 812, 787,
|
|
||||||
$ 3577 /
|
|
||||||
DATA ( MM( 35, J ), J = 1, 4 ) / 1267, 234, 2125,
|
|
||||||
$ 77 /
|
|
||||||
DATA ( MM( 36, J ), J = 1, 4 ) / 1828, 641, 2364,
|
|
||||||
$ 3761 /
|
|
||||||
DATA ( MM( 37, J ), J = 1, 4 ) / 164, 4005, 2460,
|
|
||||||
$ 2149 /
|
|
||||||
DATA ( MM( 38, J ), J = 1, 4 ) / 3798, 1122, 257,
|
|
||||||
$ 1449 /
|
|
||||||
DATA ( MM( 39, J ), J = 1, 4 ) / 3087, 3135, 1574,
|
|
||||||
$ 3005 /
|
|
||||||
DATA ( MM( 40, J ), J = 1, 4 ) / 2400, 2640, 3912,
|
|
||||||
$ 225 /
|
|
||||||
DATA ( MM( 41, J ), J = 1, 4 ) / 2870, 2302, 1216,
|
|
||||||
$ 85 /
|
|
||||||
DATA ( MM( 42, J ), J = 1, 4 ) / 3876, 40, 3248,
|
|
||||||
$ 3673 /
|
|
||||||
DATA ( MM( 43, J ), J = 1, 4 ) / 1905, 1832, 3401,
|
|
||||||
$ 3117 /
|
|
||||||
DATA ( MM( 44, J ), J = 1, 4 ) / 1593, 2247, 2124,
|
|
||||||
$ 3089 /
|
|
||||||
DATA ( MM( 45, J ), J = 1, 4 ) / 1797, 2034, 2762,
|
|
||||||
$ 1349 /
|
|
||||||
DATA ( MM( 46, J ), J = 1, 4 ) / 1234, 2637, 149,
|
|
||||||
$ 2057 /
|
|
||||||
DATA ( MM( 47, J ), J = 1, 4 ) / 3460, 1287, 2245,
|
|
||||||
$ 413 /
|
|
||||||
DATA ( MM( 48, J ), J = 1, 4 ) / 328, 1691, 166,
|
|
||||||
$ 65 /
|
|
||||||
DATA ( MM( 49, J ), J = 1, 4 ) / 2861, 496, 466,
|
|
||||||
$ 1845 /
|
|
||||||
DATA ( MM( 50, J ), J = 1, 4 ) / 1950, 1597, 4018,
|
|
||||||
$ 697 /
|
|
||||||
DATA ( MM( 51, J ), J = 1, 4 ) / 617, 2394, 1399,
|
|
||||||
$ 3085 /
|
|
||||||
DATA ( MM( 52, J ), J = 1, 4 ) / 2070, 2584, 190,
|
|
||||||
$ 3441 /
|
|
||||||
DATA ( MM( 53, J ), J = 1, 4 ) / 3331, 1843, 2879,
|
|
||||||
$ 1573 /
|
|
||||||
DATA ( MM( 54, J ), J = 1, 4 ) / 769, 336, 153,
|
|
||||||
$ 3689 /
|
|
||||||
DATA ( MM( 55, J ), J = 1, 4 ) / 1558, 1472, 2320,
|
|
||||||
$ 2941 /
|
|
||||||
DATA ( MM( 56, J ), J = 1, 4 ) / 2412, 2407, 18,
|
|
||||||
$ 929 /
|
|
||||||
DATA ( MM( 57, J ), J = 1, 4 ) / 2800, 433, 712,
|
|
||||||
$ 533 /
|
|
||||||
DATA ( MM( 58, J ), J = 1, 4 ) / 189, 2096, 2159,
|
|
||||||
$ 2841 /
|
|
||||||
DATA ( MM( 59, J ), J = 1, 4 ) / 287, 1761, 2318,
|
|
||||||
$ 4077 /
|
|
||||||
DATA ( MM( 60, J ), J = 1, 4 ) / 2045, 2810, 2091,
|
|
||||||
$ 721 /
|
|
||||||
DATA ( MM( 61, J ), J = 1, 4 ) / 1227, 566, 3443,
|
|
||||||
$ 2821 /
|
|
||||||
DATA ( MM( 62, J ), J = 1, 4 ) / 2838, 442, 1510,
|
|
||||||
$ 2249 /
|
|
||||||
DATA ( MM( 63, J ), J = 1, 4 ) / 209, 41, 449,
|
|
||||||
$ 2397 /
|
|
||||||
DATA ( MM( 64, J ), J = 1, 4 ) / 2770, 1238, 1956,
|
|
||||||
$ 2817 /
|
|
||||||
DATA ( MM( 65, J ), J = 1, 4 ) / 3654, 1086, 2201,
|
|
||||||
$ 245 /
|
|
||||||
DATA ( MM( 66, J ), J = 1, 4 ) / 3993, 603, 3137,
|
|
||||||
$ 1913 /
|
|
||||||
DATA ( MM( 67, J ), J = 1, 4 ) / 192, 840, 3399,
|
|
||||||
$ 1997 /
|
|
||||||
DATA ( MM( 68, J ), J = 1, 4 ) / 2253, 3168, 1321,
|
|
||||||
$ 3121 /
|
|
||||||
DATA ( MM( 69, J ), J = 1, 4 ) / 3491, 1499, 2271,
|
|
||||||
$ 997 /
|
|
||||||
DATA ( MM( 70, J ), J = 1, 4 ) / 2889, 1084, 3667,
|
|
||||||
$ 1833 /
|
|
||||||
DATA ( MM( 71, J ), J = 1, 4 ) / 2857, 3438, 2703,
|
|
||||||
$ 2877 /
|
|
||||||
DATA ( MM( 72, J ), J = 1, 4 ) / 2094, 2408, 629,
|
|
||||||
$ 1633 /
|
|
||||||
DATA ( MM( 73, J ), J = 1, 4 ) / 1818, 1589, 2365,
|
|
||||||
$ 981 /
|
|
||||||
DATA ( MM( 74, J ), J = 1, 4 ) / 688, 2391, 2431,
|
|
||||||
$ 2009 /
|
|
||||||
DATA ( MM( 75, J ), J = 1, 4 ) / 1407, 288, 1113,
|
|
||||||
$ 941 /
|
|
||||||
DATA ( MM( 76, J ), J = 1, 4 ) / 634, 26, 3922,
|
|
||||||
$ 2449 /
|
|
||||||
DATA ( MM( 77, J ), J = 1, 4 ) / 3231, 512, 2554,
|
|
||||||
$ 197 /
|
|
||||||
DATA ( MM( 78, J ), J = 1, 4 ) / 815, 1456, 184,
|
|
||||||
$ 2441 /
|
|
||||||
DATA ( MM( 79, J ), J = 1, 4 ) / 3524, 171, 2099,
|
|
||||||
$ 285 /
|
|
||||||
DATA ( MM( 80, J ), J = 1, 4 ) / 1914, 1677, 3228,
|
|
||||||
$ 1473 /
|
|
||||||
DATA ( MM( 81, J ), J = 1, 4 ) / 516, 2657, 4012,
|
|
||||||
$ 2741 /
|
|
||||||
DATA ( MM( 82, J ), J = 1, 4 ) / 164, 2270, 1921,
|
|
||||||
$ 3129 /
|
|
||||||
DATA ( MM( 83, J ), J = 1, 4 ) / 303, 2587, 3452,
|
|
||||||
$ 909 /
|
|
||||||
DATA ( MM( 84, J ), J = 1, 4 ) / 2144, 2961, 3901,
|
|
||||||
$ 2801 /
|
|
||||||
DATA ( MM( 85, J ), J = 1, 4 ) / 3480, 1970, 572,
|
|
||||||
$ 421 /
|
|
||||||
DATA ( MM( 86, J ), J = 1, 4 ) / 119, 1817, 3309,
|
|
||||||
$ 4073 /
|
|
||||||
DATA ( MM( 87, J ), J = 1, 4 ) / 3357, 676, 3171,
|
|
||||||
$ 2813 /
|
|
||||||
DATA ( MM( 88, J ), J = 1, 4 ) / 837, 1410, 817,
|
|
||||||
$ 2337 /
|
|
||||||
DATA ( MM( 89, J ), J = 1, 4 ) / 2826, 3723, 3039,
|
|
||||||
$ 1429 /
|
|
||||||
DATA ( MM( 90, J ), J = 1, 4 ) / 2332, 2803, 1696,
|
|
||||||
$ 1177 /
|
|
||||||
DATA ( MM( 91, J ), J = 1, 4 ) / 2089, 3185, 1256,
|
|
||||||
$ 1901 /
|
|
||||||
DATA ( MM( 92, J ), J = 1, 4 ) / 3780, 184, 3715,
|
|
||||||
$ 81 /
|
|
||||||
DATA ( MM( 93, J ), J = 1, 4 ) / 1700, 663, 2077,
|
|
||||||
$ 1669 /
|
|
||||||
DATA ( MM( 94, J ), J = 1, 4 ) / 3712, 499, 3019,
|
|
||||||
$ 2633 /
|
|
||||||
DATA ( MM( 95, J ), J = 1, 4 ) / 150, 3784, 1497,
|
|
||||||
$ 2269 /
|
|
||||||
DATA ( MM( 96, J ), J = 1, 4 ) / 2000, 1631, 1101,
|
|
||||||
$ 129 /
|
|
||||||
DATA ( MM( 97, J ), J = 1, 4 ) / 3375, 1925, 717,
|
|
||||||
$ 1141 /
|
|
||||||
DATA ( MM( 98, J ), J = 1, 4 ) / 1621, 3912, 51,
|
|
||||||
$ 249 /
|
|
||||||
DATA ( MM( 99, J ), J = 1, 4 ) / 3090, 1398, 981,
|
|
||||||
$ 3917 /
|
|
||||||
DATA ( MM( 100, J ), J = 1, 4 ) / 3765, 1349, 1978,
|
|
||||||
$ 2481 /
|
|
||||||
DATA ( MM( 101, J ), J = 1, 4 ) / 1149, 1441, 1813,
|
|
||||||
$ 3941 /
|
|
||||||
DATA ( MM( 102, J ), J = 1, 4 ) / 3146, 2224, 3881,
|
|
||||||
$ 2217 /
|
|
||||||
DATA ( MM( 103, J ), J = 1, 4 ) / 33, 2411, 76,
|
|
||||||
$ 2749 /
|
|
||||||
DATA ( MM( 104, J ), J = 1, 4 ) / 3082, 1907, 3846,
|
|
||||||
$ 3041 /
|
|
||||||
DATA ( MM( 105, J ), J = 1, 4 ) / 2741, 3192, 3694,
|
|
||||||
$ 1877 /
|
|
||||||
DATA ( MM( 106, J ), J = 1, 4 ) / 359, 2786, 1682,
|
|
||||||
$ 345 /
|
|
||||||
DATA ( MM( 107, J ), J = 1, 4 ) / 3316, 382, 124,
|
|
||||||
$ 2861 /
|
|
||||||
DATA ( MM( 108, J ), J = 1, 4 ) / 1749, 37, 1660,
|
|
||||||
$ 1809 /
|
|
||||||
DATA ( MM( 109, J ), J = 1, 4 ) / 185, 759, 3997,
|
|
||||||
$ 3141 /
|
|
||||||
DATA ( MM( 110, J ), J = 1, 4 ) / 2784, 2948, 479,
|
|
||||||
$ 2825 /
|
|
||||||
DATA ( MM( 111, J ), J = 1, 4 ) / 2202, 1862, 1141,
|
|
||||||
$ 157 /
|
|
||||||
DATA ( MM( 112, J ), J = 1, 4 ) / 2199, 3802, 886,
|
|
||||||
$ 2881 /
|
|
||||||
DATA ( MM( 113, J ), J = 1, 4 ) / 1364, 2423, 3514,
|
|
||||||
$ 3637 /
|
|
||||||
DATA ( MM( 114, J ), J = 1, 4 ) / 1244, 2051, 1301,
|
|
||||||
$ 1465 /
|
|
||||||
DATA ( MM( 115, J ), J = 1, 4 ) / 2020, 2295, 3604,
|
|
||||||
$ 2829 /
|
|
||||||
DATA ( MM( 116, J ), J = 1, 4 ) / 3160, 1332, 1888,
|
|
||||||
$ 2161 /
|
|
||||||
DATA ( MM( 117, J ), J = 1, 4 ) / 2785, 1832, 1836,
|
|
||||||
$ 3365 /
|
|
||||||
DATA ( MM( 118, J ), J = 1, 4 ) / 2772, 2405, 1990,
|
|
||||||
$ 361 /
|
|
||||||
DATA ( MM( 119, J ), J = 1, 4 ) / 1217, 3638, 2058,
|
|
||||||
$ 2685 /
|
|
||||||
DATA ( MM( 120, J ), J = 1, 4 ) / 1822, 3661, 692,
|
|
||||||
$ 3745 /
|
|
||||||
DATA ( MM( 121, J ), J = 1, 4 ) / 1245, 327, 1194,
|
|
||||||
$ 2325 /
|
|
||||||
DATA ( MM( 122, J ), J = 1, 4 ) / 2252, 3660, 20,
|
|
||||||
$ 3609 /
|
|
||||||
DATA ( MM( 123, J ), J = 1, 4 ) / 3904, 716, 3285,
|
|
||||||
$ 3821 /
|
|
||||||
DATA ( MM( 124, J ), J = 1, 4 ) / 2774, 1842, 2046,
|
|
||||||
$ 3537 /
|
|
||||||
DATA ( MM( 125, J ), J = 1, 4 ) / 997, 3987, 2107,
|
|
||||||
$ 517 /
|
|
||||||
DATA ( MM( 126, J ), J = 1, 4 ) / 2573, 1368, 3508,
|
|
||||||
$ 3017 /
|
|
||||||
DATA ( MM( 127, J ), J = 1, 4 ) / 1148, 1848, 3525,
|
|
||||||
$ 2141 /
|
|
||||||
DATA ( MM( 128, J ), J = 1, 4 ) / 545, 2366, 3801,
|
|
||||||
$ 1537 /
|
|
||||||
* ..
|
|
||||||
* .. Executable Statements ..
|
|
||||||
*
|
|
||||||
I1 = ISEED( 1 )
|
|
||||||
I2 = ISEED( 2 )
|
|
||||||
I3 = ISEED( 3 )
|
|
||||||
I4 = ISEED( 4 )
|
|
||||||
*
|
|
||||||
DO 10 I = 1, MIN( N, LV )
|
|
||||||
*
|
|
||||||
20 CONTINUE
|
|
||||||
*
|
|
||||||
* Multiply the seed by i-th power of the multiplier modulo 2**48
|
|
||||||
*
|
|
||||||
IT4 = I4*MM( I, 4 )
|
|
||||||
IT3 = IT4 / IPW2
|
|
||||||
IT4 = IT4 - IPW2*IT3
|
|
||||||
IT3 = IT3 + I3*MM( I, 4 ) + I4*MM( I, 3 )
|
|
||||||
IT2 = IT3 / IPW2
|
|
||||||
IT3 = IT3 - IPW2*IT2
|
|
||||||
IT2 = IT2 + I2*MM( I, 4 ) + I3*MM( I, 3 ) + I4*MM( I, 2 )
|
|
||||||
IT1 = IT2 / IPW2
|
|
||||||
IT2 = IT2 - IPW2*IT1
|
|
||||||
IT1 = IT1 + I1*MM( I, 4 ) + I2*MM( I, 3 ) + I3*MM( I, 2 ) +
|
|
||||||
$ I4*MM( I, 1 )
|
|
||||||
IT1 = MOD( IT1, IPW2 )
|
|
||||||
*
|
|
||||||
* Convert 48-bit integer to a real number in the interval (0,1)
|
|
||||||
*
|
|
||||||
X( I ) = R*( REAL( IT1 )+R*( REAL( IT2 )+R*( REAL( IT3 )+R*
|
|
||||||
$ REAL( IT4 ) ) ) )
|
|
||||||
*
|
|
||||||
IF (X( I ).EQ.1.0) THEN
|
|
||||||
* If a real number has n bits of precision, and the first
|
|
||||||
* n bits of the 48-bit integer above happen to be all 1 (which
|
|
||||||
* will occur about once every 2**n calls), then X( I ) will
|
|
||||||
* be rounded to exactly 1.0. In IEEE single precision arithmetic,
|
|
||||||
* this will happen relatively often since n = 24.
|
|
||||||
* Since X( I ) is not supposed to return exactly 0.0 or 1.0,
|
|
||||||
* the statistically correct thing to do in this situation is
|
|
||||||
* simply to iterate again.
|
|
||||||
* N.B. the case X( I ) = 0.0 should not be possible.
|
|
||||||
I1 = I1 + 2
|
|
||||||
I2 = I2 + 2
|
|
||||||
I3 = I3 + 2
|
|
||||||
I4 = I4 + 2
|
|
||||||
GOTO 20
|
|
||||||
END IF
|
|
||||||
*
|
|
||||||
10 CONTINUE
|
|
||||||
*
|
|
||||||
* Return final value of seed
|
|
||||||
*
|
|
||||||
ISEED( 1 ) = IT1
|
|
||||||
ISEED( 2 ) = IT2
|
|
||||||
ISEED( 3 ) = IT3
|
|
||||||
ISEED( 4 ) = IT4
|
|
||||||
RETURN
|
|
||||||
*
|
|
||||||
* End of SLARUV
|
|
||||||
*
|
|
||||||
END
|
|
||||||
@ -56,7 +56,7 @@ c-----
|
|||||||
integer NL, NL2, NLAY
|
integer NL, NL2, NLAY
|
||||||
parameter(NL=200,NLAY=200,NL2=NL+NL)
|
parameter(NL=200,NLAY=200,NL2=NL+NL)
|
||||||
integer NP
|
integer NP
|
||||||
parameter (NP=80)
|
parameter (NP=60)
|
||||||
|
|
||||||
c-----
|
c-----
|
||||||
c LIN - unit for FORTRAN read from terminal
|
c LIN - unit for FORTRAN read from terminal
|
||||||
@ -308,7 +308,7 @@ c ----- print *, itst,iq,t(k),t1a,t1b,cc0,cc1,gvel
|
|||||||
if(iverb(ifunc).eq.0)then
|
if(iverb(ifunc).eq.0)then
|
||||||
iverb(ifunc) = 1
|
iverb(ifunc) = 1
|
||||||
write(LOT,*)'improper initial value in disper - no zero found'
|
write(LOT,*)'improper initial value in disper - no zero found'
|
||||||
!write(*,*)'WARNING:improper initial value in disper - no zero found'
|
write(*,*)'WARNING:improper initial value in disper - no zero found'
|
||||||
write(LOT,*)'in fundamental mode '
|
write(LOT,*)'in fundamental mode '
|
||||||
write(LOT,*)'This may be due to low velocity zone '
|
write(LOT,*)'This may be due to low velocity zone '
|
||||||
write(LOT,*)'causing reverse phase velocity dispersion, '
|
write(LOT,*)'causing reverse phase velocity dispersion, '
|
||||||
@ -500,7 +500,7 @@ c 2 - Rayleigh Wave
|
|||||||
c iflag I*4 0 - Initialize
|
c iflag I*4 0 - Initialize
|
||||||
c 1 - Make model for Love or Rayleigh Wave
|
c 1 - Make model for Love or Rayleigh Wave
|
||||||
c-----
|
c-----
|
||||||
parameter(NL=200,NP=80)
|
parameter(NL=200,NP=60)
|
||||||
real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
|
real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
|
||||||
integer mmax,llw
|
integer mmax,llw
|
||||||
c common/modl/ d,a,b,rho,rtp,dtp,btp
|
c common/modl/ d,a,b,rho,rtp,dtp,btp
|
||||||
@ -563,7 +563,7 @@ c nev = 0 force interval halving
|
|||||||
c nev = 1 permit neville iteration if conditions are proper
|
c nev = 1 permit neville iteration if conditions are proper
|
||||||
c nev = 2 neville iteration is being used
|
c nev = 2 neville iteration is being used
|
||||||
c-----
|
c-----
|
||||||
parameter (NL=200,NP=80)
|
parameter (NL=200,NP=60)
|
||||||
implicit double precision (a-h,o-z)
|
implicit double precision (a-h,o-z)
|
||||||
real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
|
real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
|
||||||
dimension x(20),y(20)
|
dimension x(20),y(20)
|
||||||
@ -704,7 +704,7 @@ c
|
|||||||
function dltar1(wvno,omega,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi)
|
function dltar1(wvno,omega,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi)
|
||||||
c find SH dispersion values.
|
c find SH dispersion values.
|
||||||
c
|
c
|
||||||
parameter (NL=200,NP=80)
|
parameter (NL=200,NP=60)
|
||||||
implicit double precision (a-h,o-z)
|
implicit double precision (a-h,o-z)
|
||||||
real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
|
real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
|
||||||
integer llw,mmax
|
integer llw,mmax
|
||||||
@ -768,7 +768,7 @@ c
|
|||||||
c & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz)
|
c & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz)
|
||||||
c find P-SV dispersion values.
|
c find P-SV dispersion values.
|
||||||
c
|
c
|
||||||
parameter (NL=200,NP=80)
|
parameter (NL=200,NP=60)
|
||||||
implicit double precision (a-h,o-z)
|
implicit double precision (a-h,o-z)
|
||||||
dimension e(5),ee(5),ca(5,5)
|
dimension e(5),ee(5),ca(5,5)
|
||||||
real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
|
real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user