-
Notifications
You must be signed in to change notification settings - Fork 6
/
make_initial_model_dix.m
executable file
·526 lines (456 loc) · 17.2 KB
/
make_initial_model_dix.m
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
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PROGRAM:
% make_initial_model_dix.m
%
% PROGRAMMERS:
% Matt Haney and Victor Tsai
%
% Last revision date:
% 26 April 2017
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This code is distributed as part of the source-code package
% raylee_inversion_codes
% that accompanies Haney and Tsai (2017). The package can be downloaded
% from the Geophysics source-code archive at
% http://software.seg.org/2017/0003/index.html
% Use of this code is subject to acceptance of the terms and conditions
% that can be found at http://software.seg.org/disclaimer.txt
% Copyright 2017 by The Society of Exploration Geophysicists (SEG)
% Reference:
% Haney, M. M., Tsai, V. C. (2017) Perturbational and nonperturbational
% inversion of Rayleigh-wave velocities, Geophysics, 82(3), F15-F28.
% doi: 10.1190/geo2016-0397.1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Program make_initial_model_dix is a Matlab script to make an
% initial model for use in iterative Rayleigh/Scholte wave phase velocity
% inversion using the Dix method (Haney and Tsai, 2015).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Input: Listed below at beginning of script
%
% Output:
%
% files velocity_values.txt, velocity_values_errs.txt,
% frequency_values.txt, mode_values.txt, vtype_values.txt
% vp_init.txt, vs_init.txt, rho_init.txt, vpf.txt, rhof.txt,
% grid_values_solid.txt, grid_values_fluid.txt, input_params.txt
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is a script
clear all
% time this
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% begin input parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% velocity data, a row vector from low to high frequency
c = load('modx_phase_vels.ascii');
% frequency raster, a row vector from low to high frequency
fr = load('modx_freqs.ascii');
% estimates of high and low error bars
c_low = c*(1-.02);
c_hi = c*(1+.02);
% Flag for which inversion to do
% 1 = Rayleigh wave phase velocity, homogeneous form
% 2 = Rayleigh wave phase velocity, power-law form w/0.25 Poisson's ratio
% 3 = Love wave phase velocity, power-law form
% 4 = Rayleigh wave group velocity, homogeneous form
% 5 = Rayleigh wave phase velocity, power-law form w/0.3 Poisson's ratio
% 6 = Rayleigh wave group velocity, power-law form w/0.25 Poisson's ratio
% 7 = Love wave group velocity, power-law form
% 8 = Rayleigh wave group velocity, power-law form w/0.3 Poisson's ratio
inv_flag_ar = 1;
% multiplicative coefficient between wavelength and sensitivity depth
xcof = 0.5;
% number of layers per wavelength at the sensitivity depth
lrho = 20;
% range of correlation length factors to scan over
lcormultv = [10:10:1000];
% range of model variance factors to scan over
sigmscalev = [1:.1:20];
% low value of acceptable chi-squared
chiwinlo = 1;
% high value of acceptable chi-squared
chiwinhi = 1.5;
% multiple of minumum depth for model to extend to
zminfact = 1;
% multiple of maximum depth for model to extend to
zmaxfact = 4;
% Poisson's ratio (used for homogeneous model formulation and for bulding
% the output Vp model from the Vs model)
v = 0.45;
% Nominal power law exponent (used for group velocity power law forms)
alph = 0.3;
% These parameters are needed simply to pass output to nonlinear inversion:
% Nominal density value (used for building the output density model)
rhod = 2000; %kg/m^3
% Number of fluid layers to place above model
Nnf = 0;
% Thickness of fluid layers placed above model (if Nnf=0, not subsequently
% used)
hfval = 200; % m
% Density of fluid layers (if Nnf=0, not subsequently used)
rhovfal = 1000; % kg/m^3
% Velocity of fluid layers (if Nnf=0, not subsequently used)
vpvfal = 1500; % m/s
% Flag whether Vp/Vs ratio is to be fixed in nonlinear inversion
pratioflag = 1;
% Max number of nonlinear updates
nupdats = 15;
% Model standard deviation factor
mstdfact = 2;
% Smoothness scale, typically less than zmax
smscl = 10; % m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end input parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% build layered model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% min and max depths of exponentially increasing layer thicknesses
zmin = min(xcof*(c./fr))*zminfact;
zmax = max(xcof*(c./fr))*zmaxfact;
% maximum number of layers
Nmax = (xcof*lrho)*log(zmax/zmin);
% layer thickness assuming the last layer between zmin and zmax is smaller
% than it should be
thkso = diff([0 zmin*exp([0:1:Nmax]/(xcof*lrho)) zmax]);
% put small layers in the shallow part above zmin
thks = [ thkso(2)*ones(1,ceil(zmin/thkso(2))) thkso(2:end) ];
% the first layer is a little thinner
thks(1) = thks(1) + (zmin - thkso(2)*ceil(zmin/thkso(2)));
z = [0 cumsum(thks) inf];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% scan over two regularization parameters for acceptable models
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% size of data and model
dsz = length(fr);
msz = length(z)-1;
% Set parameter for inversion
inv_flag = inv_flag_ar;
% set some factors that are used to define kernel for homogeneous form
% f is the factor that multiplies shear wave velocity to get Rayleigh
% wave velocity
f = sqrt((8/3) - ((-16 + 56*v - 40*(v^2))/...
(3*(2^(2/3))*(-1 + v)*(-11 + 78*v - 123*(v^2) + 56*(v^3) + ...
3*sqrt(3)*sqrt(-5 + 36*v - 94*(v^2) + 148*(v^3) - 165*(v^4) + ...
112*(v^5) - 32*(v^6)))^(1/3))) + ((1/(3*(-1 + v)))*(2^(2/3))*(-11 + ...
78*v - 123*(v^2) + 56*(v^3) + 3*sqrt(3)*sqrt(-5 + 36*v - 94*(v^2) + ...
148*(v^3) - 165*(v^4) + 112*(v^5) - 32*(v^6)))^(1/3)));
% make sure there isn't a small imaginary part
f = real(f);
% this is the denominator of the Dix expression in Haney and Tsai (2015),
% multiplied by 8
g = -((((-2 + f^2)^3)/((1 - f^2)^(3/2))) - ...
(((2*sqrt(2)*sqrt((-2 + f^2 + 2*v - 2*(f^2)*v)/(-1 + v))*(4 - 4*v + ...
(f^2)*(-1 + 2*v))))/((2 - 2*v + (f^2)*(-1 + 2*v)))) + ((8*(-2 + ...
f^2)*(2 - 2*f^2 + sqrt(2 - 2*(f^2))*sqrt((-2 + f^2 + 2*v - ...
2*(f^2)*v)/(-1 + v))))/((-1 +f^2)*(2*sqrt(1 - f^2) + ...
sqrt(2)*sqrt((-2 + f^2 + 2*v - 2*(f^2)*v)/(-1 + v))))));
% these are factors and exponential terms used in the homogeneous
% formulation
fctr1 = -((f^2 - 2)^2)*(8 - 8*(f^2) + f^4)/((1 - f^2)^1.5)/g;
exp1 = 2*sqrt(1 - f^2);
fctr2 = -8*(f^2 - 2)*2*(1+sqrt(1 - f^2)*sqrt(1 + ((f^2)*((2*v - 1)/...
(2 - 2*v)))))/sqrt(1 - f^2)/g;
exp2 = sqrt(1 - f^2)+sqrt(1 + ((f^2)*((2*v - 1)/(2 - 2*v))));
fctr3 = 4*sqrt(1 + ((f^2)*((2*v - 1)/(2 - 2*v))))*(16*(v-1) + ...
(f^2)*(f^2 - 8)*(2*v - 1))/(2 - 2*v + (f^2)*(2*v - 1))/g;
exp3 = 2*sqrt(1 + ((f^2)*((2*v - 1)/(2 - 2*v))));
% Calculate kernel functions (Haney and Tsai, 2015)
% f_homo = Rayleigh wave phase velocity, homogeneous form
% f_homo_U = Rayleigh wave group velocity, homogeneous form
% f_ray1 = Rayleigh wave phase velocity, power law form, Poisson's ratio
% 0.25
% f_ray2 = Rayleigh wave phase velocity, power law form, Poisson's ratio
% 0.3
% f_love = Love wave phase velocity, power law form
% angular frequency, wavenumber, and wavenumber times depth
w = 2*pi*fr; k = w./c; kz = k'*z;
% Phase velocity, homogeneous form
f_homo = fctr3*exp(-exp3*kz)+fctr2*exp(-exp2*kz)+fctr1*exp(-exp1*kz);
f_homo(:,end) = 0;
% Group velocity, homogeneous form
f_homo_U = fctr3*(1-exp3*kz).*exp(-exp3*kz)+...
fctr2*(1-exp2*kz).*exp(-exp2*kz)+...
fctr1*(1-exp1*kz).*exp(-exp1*kz);
f_homo_U(:,end) = 0;
% Poisson ratio of 0.3
f_ray2 = -130.253*exp(-1.8362*kz) + 6.88812*exp(-1.7556*kz) + ...
271.540*exp(-1.7380*kz) - 12.1077*exp(-1.6859*kz) - ...
183.478*exp(-1.6750*kz) + 1.70238*exp(-1.6574*kz) - ...
142.361*exp(-1.6398*kz) + 341.126*exp(-1.6053*kz) + ...
4.76194*exp(-1.5877*kz) - 159.030*exp(-1.5356*kz);
% Poisson ratio of 0.25
f_ray1 = -103.14*exp(-1.8586*kz) + 6.1446*exp(-1.7714*kz) + ...
217.120*exp(-1.7555*kz) - 10.312*exp(-1.7012*kz) - ...
160.68*exp(-1.6842*kz) + 1.2856*exp(-1.6683*kz) - ...
115.00*exp(-1.6524*kz) + 294.66*exp(-1.6140*kz) + ...
4.0924*exp(-1.5981*kz) - 135.49*exp(-1.5438*kz);
% Love waves
f_love = -(1+0.85^2)*exp(-2*0.85*kz);
% Construct G matrix
if inv_flag == 1
G = diff(f_homo,1,2);
elseif inv_flag == 2
G = diff(f_ray1,1,2);
elseif inv_flag == 3
G = diff(f_love,1,2);
elseif inv_flag == 4
G = diff(f_homo_U,1,2);
elseif inv_flag == 5
G = diff(f_ray2,1,2);
elseif inv_flag == 6
G = ((1-alph)^2)*diff(f_ray1,1,2);
elseif inv_flag == 7
G = ((1-alph)^2)*diff(f_love,1,2);
else
G = ((1-alph)^2)*diff(f_ray2,1,2);
end
% Depths of layer tops
z_mid = z(1:(length(z)-1));
% Defines an array of data stdevs based on the hi and low values
sigma_d_ar = ((c_hi.^2-c_low.^2)/2);
% Data covariance matrix
Cd = diag(sigma_d_ar.^2);
% Create Xia model interp (0.5c/f) and midpoint distance matrix
% find best fit power law, find slope at shallow point, and
% do the extrapolation to the surface
% if power law gives negative velocity at surface, use linear function
% if linear function gives negative velocity at surface use shallowest
% velocity
% repeat similarly for deepest point in model
% Power law fit
[aap bbp] = polyfit(log(xcof*c./fr),log(c/0.88),1);
% Linear fit
[aal bbl] = polyfit(xcof*c./fr,c/0.88,1);
% extrapolation to surface
if (c(end)/.88 - exp(aap(2))*aap(1)*((xcof*c(end)/fr(end))^aap(1)) > 0)
czro = c(end)/.88 - exp(aap(2))*aap(1)*((xcof*c(end)/fr(end))^aap(1));
xiatrap = [c/0.88 czro];
elseif (c(end)/.88 - aal(1)*xcof*c(end)/fr(end) > 0)
czro = c(end)/.88 - aal(1)*xcof*c(end)/fr(end);
xiatrap = [c/0.88 czro];
else
czro = c(end)/.88;
xiatrap = [c/0.88 czro];
end
% extrapolation to deepest point in model
if (c(1)/.88 - exp(aap(2))*aap(1)*((xcof*c(1)/fr(1))^(aap(1)-1))*(xcof*c(1)/fr(1)-max(z_mid)) > 0)
cdeep = c(1)/.88 - exp(aap(2))*aap(1)*((xcof*c(1)/fr(1))^(aap(1)-1))*(xcof*c(1)/fr(1)-max(z_mid));
xiatrap = [cdeep xiatrap];
elseif (c(1)/.88 - aal(1)*(xcof*c(1)/fr(1)-max(z_mid)) > 0)
cdeep = c(1)/.88 - aal(1)*(xcof*c(1)/fr(1)-max(z_mid));
xiatrap = [cdeep xiatrap];
else
cdeep = c(1)/.88;
xiatrap = [cdeep xiatrap];
end
% interpolate onto model
xia_int = interp1([max(z_mid) xcof*c./fr 0],xiatrap,z_mid,'linear');
% Layer top distance matrix used in model covariance matrix
z_mid_mat = zeros(length(z_mid));
for i=1:length(z_mid)
for j=1:length(z_mid)
z_mid_mat(i,j)=abs(z_mid(i)-z_mid(j));
end
end
% counter for acceptable models
modcountr = 0;
% counter over correlation length factors
lcntr = 0;
% scan over correlation length factors
for lcormult=lcormultv
% increment counter
lcntr = lcntr + 1;
% counter over model variance factors
mcntr = 0;
% scan over model variance factors
for sigmscale=sigmscalev
% increment counter
mcntr = mcntr + 1;
% define model correlation length
% lcormult times the median layer thickness
l_cor = median(thks)*lcormult;
% Defines model stdev based on average data stdev
sigma_m = (mean(sigma_d_ar)*sigmscale);
% model covariance matrix
Cm = sigma_m^2*exp(-z_mid_mat/l_cor);
% Creates augmented data vector and G matrix
d_aug = [Cd^-0.5*c.^2'; Cm^-0.5*xia_int.^2'];
G_aug = [Cd^-0.5*G; Cm^-0.5];
% Does inversion using augmented d and G,
beta_sq_r2 = (G_aug'*G_aug)\G_aug'*d_aug;
% Chi-squared of model
chi_sqd(lcntr,mcntr) = transpose(c.^2' - G*beta_sq_r2)*...
(Cd^-1)*(c.^2' - G*beta_sq_r2)/dsz;
% Does model have negative values of squared shear velocity?
mod_neg(lcntr,mcntr) = min(beta_sq_r2>0);
% If the model is acceptable, save it in a matrix
if(chi_sqd(lcntr,mcntr) > chiwinlo & chi_sqd(lcntr,mcntr) < ...
chiwinhi & mod_neg(lcntr,mcntr) ~= 0)
modcountr = modcountr + 1;
beta_sq_r2_all(:,modcountr) = beta_sq_r2;
else
end
% End loop over model standard deviation
end
lcntr
% End loop over correlation length
end
% mean of acceptable vs models
if (modcountr > 0)
beta_sq_r2m = mean(beta_sq_r2_all')';
elseif (min(min(chi_sqd)) > chiwinhi)
% an error message if no models found
error('No acceptable models found: All chi-squareds above acceptable window. Consider expanding search grid over regularization parameters or expanding acceptable chi-squared range.');
elseif (max(max(chi_sqd)) < chiwinlo)
% an error message if no models found
error('No acceptable models found: All chi-squareds below acceptable window. Consider expanding search grid over regularization parameters or expanding acceptable chi-squared range.');
else
% an error message if no models found
error('No acceptable models found: All chi-squareds above or below acceptable window. Consider densifying search grid over regularization parameters')
end
% plot Chi-squared and range of acceptable models
figure
fsize = 14;
imagesc(sigmscalev,lcormultv,chi_sqd.*mod_neg); colormap('jet');
hh = colorbar('EastOutside','FontSize',fsize,...
'FontWeight','bold'); %,'Ytick',[-0.015:0.005:0.015]);
label = sprintf(' Chi-squared ');
set(get(hh,'YLabel'),'String',label,'FontSize',fsize,'FontWeight','bold')
hold on
contour(sigmscalev,lcormultv,chi_sqd.*mod_neg,[chiwinlo chiwinhi],'w--','LineWidth',4);
set(gca,'Fontsize',fsize,'FontWeight','bold');
title('Chi-squared misfit for Dix inversion')
xlabel('Model standard deviation factor','FontSize',fsize,'FontWeight','bold')
ylabel('Model correlation length factor','FontSize',fsize,...
'FontWeight','bold')
orient landscape
print(gcf,'-dpsc','chi_squared.ps');
% write out input files for nonlinear inversion
Nf = length(c);
Nn = length(thks);
vout = c;
voute = c_hi-c;
fks = fr;
% Mode number, 1=fundamental (only option implemented for Dix inversion)
modnv = ones(1,length(c));
% Type of velocity measurement, 0=phase, 1=group
if (inv_flag_ar == 1 | inv_flag_ar == 2 | inv_flag_ar == 3 | ...
inv_flag_ar == 5)
vtypv = zeros(1,length(c));
else
vtypv = ones(1,length(c));
end
% Output Vs velocity model
vsv = sqrt(beta_sq_r2m(1:Nn)');
% Output Vp velocity model
vpv = sqrt((2*v-2)/(2*v-1))*vsv;
% Output density model
rhov = rhod*ones(1,Nn);
vpvf = vpvfal*ones(1,Nnf);
rhovf = rhovfal*ones(1,Nnf);
h = thks;
hfv = hfval*ones(1,Nnf);
% write out velocities in single column format
fid = fopen('velocity_values.txt','w');
for ii=1:Nf
fprintf(fid,'%10.5f\n',vout(ii));
end
fclose(fid);
% write out velocity error bars in single column format
fid = fopen('velocity_values_errs.txt','w');
for ii=1:Nf
fprintf(fid,'%10.5f\n',voute(ii));
end
fclose(fid);
% write out frequencies in single column format
fid = fopen('frequency_values.txt','w');
for ii=1:Nf
fprintf(fid,'%10.5f\n',fks(ii));
end
fclose(fid);
% write out variable modnv in single column format
fid = fopen('mode_values.txt','w');
for ii=1:Nf
fprintf(fid,'%10.5f\n',modnv(ii));
end
fclose(fid);
% write out variable vtype in single column format
fid = fopen('vtype_values.txt','w');
for ii=1:Nf
fprintf(fid,'%10.5f\n',vtypv(ii));
end
fclose(fid);
if (pratioflag == 0)
% write out Vp model in single column format
fid = fopen('vp_init.txt','w');
for ii=1:length(vpv)
fprintf(fid,'%10.5f\n',vpv(ii));
end
fclose(fid);
elseif (pratioflag == 1)
fid = fopen('vp_init.txt','w');
fprintf(fid,'%10.5f\n',sqrt((2*v-2)/(2*v-1)));
fclose(fid);
else
end
% write out Vs model in single column format
fid = fopen('vs_init.txt','w');
for ii=1:length(vsv)
fprintf(fid,'%10.5f\n',vsv(ii));
end
fclose(fid);
% write out density model in single column format
fid = fopen('rho_init.txt','w');
for ii=1:length(rhov)
fprintf(fid,'%10.5f\n',rhov(ii));
end
fclose(fid);
% write out Vp model in fluid in single column format
fid = fopen('vpf.txt','w');
for ii=1:length(vpvf)
fprintf(fid,'%10.5f\n',vpvf(ii));
end
fclose(fid);
% write out density model in fluid in single column format
fid = fopen('rhof.txt','w');
for ii=1:length(rhovf)
fprintf(fid,'%10.5f\n',rhovf(ii));
end
fclose(fid);
% write out element thicknesses in solid in single column format
fid = fopen('grid_values_solid.txt','w');
for ii=1:Nn
fprintf(fid,'%10.5f\n',h(ii));
end
fclose(fid);
% write out element thicknesses in fluid in single column format
fid = fopen('grid_values_fluid.txt','w');
for ii=1:Nnf
fprintf(fid,'%10.5f\n',hfv(ii));
end
fclose(fid);
% write out input parameters file
fidt = fopen('input_params.txt','w');
fprintf(fidt,'%% input parameters for Rayleigh/Scholte wave inversion\n');
fprintf(fidt,'\n');
fprintf(fidt,'%i %% flag for fixed poisson''s ratio (0=no,1=yes)\n',pratioflag);
fprintf(fidt,'%10.5f %% smoothness scale (m)\n',smscl);
fprintf(fidt,'%10.5f %% a priori model standard deviation factor\n',mstdfact);
fprintf(fidt,'%i %% maximum number of updates (iterations)\n',nupdats);
fprintf(fidt,'%i %% number of measurements\n',Nf);
fprintf(fidt,'%i %% number of elements in solid part of model\n',Nn);
fprintf(fidt,'%i %% number of elements in fluid part of model\n',Nnf);
fprintf(fidt,'%10.5f %% lower chi squared window\n',chiwinlo);
fprintf(fidt,'%10.5f %% higher chi squared window\n',chiwinhi);
fclose(fidt);
% end the timer
toc