-
Notifications
You must be signed in to change notification settings - Fork 2
/
cap_sub.c
83 lines (73 loc) · 2.2 KB
/
cap_sub.c
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
#include "cap.h"
void taper(float *aa, int n) {
int i, m;
float tt, pi1;
float taper_fraction = 0.3;//fraction of the waveform tapered from both sides
m = rint(taper_fraction*n);
pi1 = 3.1415926/m;
for (i=0; i<m; i++) {
tt = 0.5*(1.-cos(i*pi1));
aa[i] *= tt;
aa[n-i-1] *= tt;
}
}
int discard_bad_data(int nda, DATA *obs, SOLN sol, float sig, float rms_cut[]) {
int i, j, n;
float minCC = 50.;
COMP *spt;
n = 0;
for(i=0;i<nda;i++,obs++) {
spt = obs->com;
for(j=0; j<NCP; j++,spt++) {
if (sol.cfg[i][j]<minCC && sol.error[i][j]/sig > rms_cut[j]) {
spt->on_off = 0;
n++;
}
}
}
return(n);
}
//void principal_values(float a[]) {
// int i;
// float **b, **v;
// b = matrix(1,3,1,3);
// v = matrix(1,3,1,3);
// for(i=0;i<3;i++) b[i+1][i+1]=a[i];
// b[1][2]=b[2][1]=a[3];b[1][3]=b[3][1]=a[4];b[2][3]=b[3][2]=a[5];
// jacobi(b,3,a,v,&i);
// eigsrt(a,v,3);
// for(i=0;i<3;i++) {a[i] = a[i+1]; if (a[i]<0.0001*a[1]) a[i]=0.0001*a[1];}
// free_convert_matrix(b,1,3,1,3);
// free_matrix(v,1,3,1,3);
//}
float *cutTrace(float *trace, int npts, int offset, int n) {
int m;
float *cut;
cut = (float *) calloc(n, sizeof(float));
if (cut == NULL) return cut;
m = n+offset;
if (offset<0) {
if (m>npts) m = npts;
if (m>0) memcpy(cut-offset, trace, m*sizeof(float));
} else {
if (m>npts) n = npts-offset;
if (n>0) memcpy(cut, trace+offset, n*sizeof(float));
}
return cut;
}
// This function checks (1) if observed and predicted polarities match and (2) if
// abs(radiation amp) is greater than threshold (fm_thr). if (1) and (2) both
// true then solution is permissible. KEY: Letting fm_thr be negative allows for
// possible bad observations, eg. for stations near nodal lines.
// INPUT: moment tensor (mt[3][3]), set of polarity observations (fm),
// number of observations (n) and value of threshold(fm_thr).
// OUTPUT: -1 == bad solution or 0 == permissible solution
int check_first_motion(float mt[3][3], FM *fm, int n, float fm_thr) {
int i;
FM *pt;
for(pt=fm,i=0;i<n;i++,pt++) {
if (pt->type*radpmt(mt, pt->alpha, pt->az, abs(pt->type))/abs(pt->type)<fm_thr)
return -1;
}
return 0;
}