-
Notifications
You must be signed in to change notification settings - Fork 3
/
crustal_velocities_full.py
executable file
·210 lines (201 loc) · 18.9 KB
/
crustal_velocities_full.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
import numpy as np
import math
length=input('Enter Profile length: ')
reso=input('Enter resolution along length: ')
data=np.loadtxt('post_processing_output.dat',usecols=(0,1,2,3,4,5,6))
vp_vs=input('Enter Vp/Vs value: ')
#########################################################
#### Anelastic attenuation model from Jackson et al., 2010
########################################################
AA1=816
alfa3=0.36
energi=293.0E03
volexp=1.20E-05
rgas=8.314472
pi=3.1415926
Qp=[]
Qs=[]
dsize=10
oscill=50
### data format
# X Z T P Vp Vs density Material mk??
## caluclating Vp in crust according to Brocher 2005
f=0
for i in range(len(data)):
if data[i,5]==0:
dens=data[i,6]/1000
# From emprical relations
data[i,4]= 39.128*dens - 63.064*dens**2 + 37.083*dens**3 - 9.1819*dens**4 + 0.8215*dens**5
#f=vp_vs - data[i,1]/1000
data[i,5]=data[i,4]/vp_vs
### Accroding to equation 6 Brocher BSSA 2005
#data[i,5]= 0.7858 -1.2344*data[i,4] + 0.7949*data[i,4]**2 - 0.1238*data[i,4]**3 + 0.0064*data[i,4]**4
### Accroding to equation 7 Brocher BSSA 2005
#data[i,5]= (data[i,4]-1.36)/1.16 # does not work
### Accroding to equation 8 Brocher BSSA 2005
#data[i,5]= 2.88 + 0.52*(data[i,4]-5.25) # works but not perfect
## Constant values
#data[i,4]=6.52
#data[i,5]=3.0
Qp.append(1000)
Qs.append(200)
else:
parexp=math.exp((-(energi+(volexp*data[i,3])))/(rgas*(data[i,2]+273.0E0)))
sqatt50=AA1*(((oscill*(1.0E0/(dsize*1000.0E0)))*parexp))**alfa3
Qp.append((1/sqatt50)*(9/4))
Qs.append(1/sqatt50)
f=open("post_processing_output_crust_mantle.dat","w")
for i in range(len(data)):
f.writelines(" %f %f %f %f %f %f %f %f %f \n " % (data[i,0],data[i,1],data[i,6]/1000, data[i,4],data[i,5],data[i,2],data[i,3],Qp[i],Qs[i]))
f.close()
data=np.loadtxt("post_processing_output_crust_mantle.dat",)
x_nodes=length/reso
Qp=[]
Qs=[]
for i in range(x_nodes):
s=str(str(i*reso)+"_vel.dat")
f=open(s,"w")
f.write("MODEL\n\
TEST MODEL \n\
ISOTROPIC \n\
KGS \n\
FLAT EARTH \n\
1-D \n\
CONSTANT VELOCITY \n\
LINE08 \n\
LINE09 \n\
LINE10 \n\
LINE11 \n \
HR VP VS RHO QP QS ETAP ETAS FREFP FREFS \n")
for j in range(i+ (94*i),i+ (94*i) + 94):
#print i,j
#Parameters
f.writelines("%f \t %f \t %f \t %f \t %s \t %s \t %s \t %s \t %s \t %s \n"
%(-(data[j+1,1]+data[j,1]) , data[j,4], data[j,5], data[j,6], data[j,7], data[j,8], "0.0", "0.0", "1.0", "1.0"))
f.writelines("0.000000 9.030200 4.870200 3.506800 377.930000 146.570000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 9.360100 5.080600 3.931700 413.660000 162.500000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 9.528000 5.186400 3.927300 417.320000 164.870000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 9.696200 5.292200 3.923300 419.940000 166.800000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 9.864000 5.398900 3.921800 422.550000 168.780000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 10.032000 5.504700 3.920600 425.510000 170.820000 0.0 0.0 1.0 1.0 \n")
f.writelines("0.000000 10.200000 5.610400 3.920100 428.690000 172.930000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 10.790900 5.960700 4.238700 1350.540000 549.450000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 10.922200 6.089800 4.298600 1311.170000 543.480000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.055300 6.210000 4.356500 1277.930000 537.630000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.135500 6.242400 4.411800 1269.440000 531.910000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.222800 6.279900 4.465000 1260.680000 526.320000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.306800 6.316400 4.516200 1251.690000 520.830000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.389700 6.351900 4.565400 1243.020000 515.460000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.470400 6.386000 4.592600 1234.540000 510.200000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.549300 6.418200 4.619800 1226.520000 505.050000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.626500 6.451400 4.646700 1217.910000 500.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.702000 6.482200 4.673500 1210.020000 495.050000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.776800 6.513100 4.700100 1202.040000 490.200000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.849100 6.543100 4.726600 1193.990000 485.440000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.920800 6.572800 4.752800 1186.060000 480.770000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 11.989100 6.600900 4.779000 1178.190000 476.190000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.057100 6.628500 4.805000 1170.530000 471.700000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.124700 6.655400 4.830700 1163.160000 467.290000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.191200 6.681300 4.856200 1156.040000 462.960000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.255800 6.707000 4.881700 1148.760000 458.720000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.318100 6.732300 4.906900 1141.320000 454.550000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.381300 6.757900 4.932100 1134.010000 450.450000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.442700 6.782000 4.957000 1127.020000 446.430000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.503000 6.805600 4.981700 1120.090000 442.480000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.563800 6.828900 5.006200 1108.580000 436.680000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.622600 6.851700 5.030600 1097.160000 431.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.680700 6.874300 5.054800 1085.970000 425.530000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.738400 6.897200 5.078900 1070.380000 418.410000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.795600 6.919400 5.102700 1064.230000 414.940000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.852400 6.941600 5.126400 1058.030000 411.520000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.909300 6.962500 5.149900 1048.090000 406.500000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 12.966300 6.985200 5.173200 1042.070000 403.230000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 13.022600 7.006900 5.196300 1032.140000 398.410000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 13.078600 7.028600 5.219200 1018.380000 392.160000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 13.133700 7.050400 5.242000 1008.790000 387.600000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 13.189500 7.072200 5.264600 999.440000 383.140000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 13.246500 7.093200 5.287000 990.770000 378.790000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 13.301700 7.114400 5.309200 985.630000 375.940000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 13.358400 7.136800 5.331300 976.810000 371.750000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.500000 13.415600 7.158400 5.353100 968.460000 367.650000 0.0 0.0 1.0 1.0 \n")
f.writelines("48.500000 13.474100 7.180400 5.374800 960.360000 363.640000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 13.531100 7.203100 5.396200 952.000000 359.710000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.000000 13.589900 7.225300 5.417600 940.880000 354.610000 0.0 0.0 1.0 1.0 \n")
f.writelines("0.000000 13.649800 7.248500 5.438700 933.210000 350.880000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.670000 13.649800 7.248500 5.693400 722.730000 271.740000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.660000 13.653300 7.259300 5.719600 726.870000 273.970000 0.0 0.0 1.0 1.0 \n")
'''
f.writelines("52.170000 13.657000 7.270000 5.745800 725.110000 273.970000 0.0 0.0 1.0 1.0 \n")
f.writelines("0.000000 13.660100 7.281700 5.772100 723.120000 273.970000 0.0 0.0 1.0 1.0 \n")
f.writelines("47.830000 8.000000 0.000000 9.914500 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.038200 0.000000 9.994200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.128300 0.000000 10.072200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.221300 0.000000 10.148500 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.340000 8.312200 0.000000 10.223300 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.400100 0.000000 10.296400 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.486100 0.000000 10.367900 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.569200 0.000000 10.437800 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.649600 0.000000 10.506200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.728300 0.000000 10.573100 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.803600 0.000000 10.638500 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.876100 0.000000 10.702300 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 8.946100 0.000000 10.764700 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.340000 9.013800 0.000000 10.825700 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.079200 0.000000 10.885200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.142600 0.000000 10.943400 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.204200 0.000000 11.000100 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.263400 0.000000 11.055500 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.320500 0.000000 11.109500 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.376000 0.000000 11.162300 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.429700 0.000000 11.213700 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.481400 0.000000 11.263900 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.340000 9.530600 0.000000 11.312700 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.577700 0.000000 11.360400 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.623200 0.000000 11.406900 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.667300 0.000000 11.452100 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.710000 0.000000 11.496200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.751300 0.000000 11.539100 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.791400 0.000000 11.580900 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.830400 0.000000 11.621600 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.868200 0.000000 11.661200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.340000 9.905100 0.000000 11.699800 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.941000 0.000000 11.737300 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 9.976100 0.000000 11.773700 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.010300 0.000000 11.809200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.043900 0.000000 11.843700 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.076800 0.000000 11.877200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.109500 0.000000 11.909800 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.141500 0.000000 11.941400 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.173900 0.000000 11.972200 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.340000 10.204900 0.000000 12.000100 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.232900 0.000000 12.031100 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.256500 0.000000 12.059300 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.330000 10.274500 0.000000 12.086700 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("49.930000 10.285400 0.000000 12.113300 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("0.000000 10.289000 0.000000 12.139100 57822.000000 0.000000 0.0 0.0 1.0 1.0 \n")
f.writelines("51.110000 11.042700 3.504300 12.703700 633.260000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.058500 3.518700 12.728900 629.890000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.720000 11.071800 3.531400 12.753000 626.870000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.085000 3.543500 12.776000 624.080000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.098300 3.555100 12.798000 621.500000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.116600 3.566100 12.818800 619.710000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.720000 11.131600 3.576500 12.838700 617.780000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.145700 3.586400 12.857400 615.930000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.159000 3.595700 12.875100 614.210000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.171500 3.604400 12.891700 612.620000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.720000 11.183200 3.612600 12.907200 611.120000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.194100 3.620200 12.921700 609.740000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.204100 3.627200 12.935100 608.480000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.213400 3.633700 12.947400 607.310000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.720000 11.221900 3.639600 12.958600 606.260000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.229500 3.645000 12.968800 605.280000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.236400 3.649800 12.977900 604.440000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.242400 3.654000 12.985900 603.690000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.720000 11.247700 3.657700 12.992900 603.040000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.252100 3.660800 12.998800 602.490000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.255700 3.663300 13.003600 602.050000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.258600 3.665300 13.007400 601.700000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.720000 11.260600 3.666700 13.010000 601.460000 85.030000 0.0 0.0 1.0 1.0 \n")
f.writelines("50.710000 11.261800 3.667500 13.011700 601.320000 85.030000 0.0 0.0 1.0 1.0 \n")
'''
f.close()