/
measure_windows_xcorr.f90
98 lines (72 loc) · 3.02 KB
/
measure_windows_xcorr.f90
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
!
! $Id:$
!
!----------------------------------------------------------------------
subroutine measure_windows_xcorr
use xcorr_constants
use seismo_variables
use measurement_variables
implicit none
double precision, dimension(npt) :: syn_lp_local, obs_lp_local
double precision :: dtau_local, dlnA_local
double precision :: CC_after_local, dlnA_after_local, tshift_after_local
double precision, dimension(npt) :: obs_win_local, syn_win_local, obs_rec_local
integer :: iwin, nlen
double precision :: lp_min, lp, win_length
logical :: quality_ok
double precision, dimension(npt) :: fp_local, fq_local
! Loop over windows
! ---------------------------------------------------------------------
write(*,*) "Cross-correlation measurements"
do iwin = 1, num_win
win_length = win_end(iwin) - win_start(iwin)
nlen=win_length/dt
! initialise global measurements
! ---------------------------------------------------------------------
n_freq(iwin) = 0
fr(:,iwin) = 0
dtau_w(:,iwin) = 0; dlnA_w(:,iwin) = 0
F1_after(iwin) = 0; F2_after(iwin) = 0
Tshift_after(iwin) = 0; CC_after(iwin) = 0 ; dlnA_after = 0
! minimum lowpass corner frequency for measurement = WIN_LP_FREQ
! -------------------------------------------------------------------
lp_min = 1.0/WIN_LP_PERIOD
lp=lp_min
quality_ok=.true.
! initialise all local measurements and erros to zero
! ---------------------------------------------------------------------
dtau_local = 0 ; dlnA_local = 0
obs_win_local(:) = 0 ; syn_win_local(:) = 0; obs_rec_local(:) = 0
fp_local(:) = 0 ; fq_local(:) = 0
obs_lp_local(:)=0; syn_lp_local(:)=0
obs_lp_local(1:npts)=obs_lp(1:npts)
syn_lp_local(1:npts)=synt_lp(1:npts)
! make the xcorr measurement
! -------------------------------------------------------------------
call xcorr_measure(win_start(iwin),win_end(iwin),&
syn_lp_local,obs_lp_local,npts,dt,b,&
tshift_after_local,CC_after_local,dlnA_after_local,&
dtau_local,dlnA_local,&
obs_win_local,syn_win_local,obs_rec_local,&
CALC_ADJ,fp_local,fq_local,DEBUG)
! ---------------------------------------------------------------------
! measurements
n_freq(iwin) = 1
fr(1,iwin) = lp
dtau_w(1,iwin) = dtau_local
dlnA_w(1,iwin) = dlnA_local
Tshift_after(iwin) = tshift_after_local
CC_after(iwin)=CC_after_local
dlnA_after(iwin)=dlnA_after_local
syn_win(1:nlen,iwin) = syn_win_local(1:nlen)
obs_win(1:nlen,iwin) = obs_win_local(1:nlen)
obs_rec(1:nlen,iwin) = obs_rec_local(1:nlen)
! version without frequency iteration
write(*,'("Window ",i2," : dTau = ",f5.2," dlnA = ",f5.3," CC = ",f5.3," lp at ",f6.2," s ")') iwin, &
dtau_local, dlnA_local, CC_after_local, 1/(lp)
if(calc_adj) then
fp_adj_win(1:nlen,iwin) = fp_local(1:nlen)
fq_adj_win(1:nlen,iwin) = fq_local(1:nlen)
endif
enddo
end subroutine measure_windows_xcorr