calc_res.c 8.22 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
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
/*-----------------------------------------------------------------------------------------
 * Copyright (C) 2013  For the list of authors, see file AUTHORS.
 *
 * This file is part of DENISE.
 * 
 * DENISE is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
 * 
 * DENISE is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   Calculate Data Residuals                                  
 *   last update 29/03/08, D.Koehn
 *  ----------------------------------------------------------------------*/
#include "fd.h"

double calc_res(float **sectiondata, float **section, float **sectiondiff, float **sectiondiffold, int ntr, int ns, int LNORM, float L2, int itest, int sws, int swstestshot, int ntr_glob, int **recpos_loc, int nsrc_glob, int ishot, int iter){

/* declaration of variables */
extern float DT, WATERLEVEL_LNORM8;
extern int REC1, REC2, MYID, ACOUSTIC;
extern int TRKILL, NORMALIZE, FC, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY;
float RMS, signL1, intseis;
int Lcount,i,j,invtime,k,h, umax=0;
float l2;
float abs_section, abs_sectiondata, sectiondata_mult_section;
float intseis_s, intseis_sd;
float **intseis_section=NULL, **intseis_sectiondata=NULL;
float **intseis_sectiondata_envelope=NULL, **intseis_section_envelope=NULL, **intseis_section_hilbert=NULL, **dummy_1=NULL, **dummy_2=NULL;
40

Tilman Steinweg's avatar
Tilman Steinweg committed
41
42
43
44
if(LNORM==8){
	intseis_section_envelope = matrix(1,ntr,1,ns);
	intseis_sectiondata_envelope = matrix(1,ntr,1,ns);
	intseis_section_hilbert = matrix(1,ntr,1,ns);
45
	dummy_1 = matrix(1,ntr,1,ns);
Tilman Steinweg's avatar
Tilman Steinweg committed
46
	dummy_2 = matrix(1,ntr,1,ns);
47
48
}

Tilman Steinweg's avatar
Tilman Steinweg committed
49
50
51
52
53
54
55
56
57
58
59
intseis_section = matrix(1,ntr,1,ns);  /* declaration of variables for integration */
intseis_sectiondata = matrix(1,ntr,1,ns);

/* sectiondiff will be set to zero */
umax=ntr*ns;
zero(&sectiondiff[1][1],umax);

/* TRACE KILLING */
int ** kill_tmp, *kill_vector;	/* declaration of variables for trace killing */
char trace_kill_file[STRING_SIZE];	
FILE *ftracekill;
60

Tilman Steinweg's avatar
Tilman Steinweg committed
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
if(TRKILL==1){
	/* sectiondiff will be set to zero */
	umax=ntr*ns;
	zero(&sectiondiff[1][1],umax);
	
	kill_tmp = imatrix(1,ntr_glob,1,nsrc_glob);
	kill_vector = ivector(1,ntr);

	ftracekill=fopen(TRKILL_FILE,"r");

	if (ftracekill==NULL) err(" Trace kill file could not be opened!");

	for(i=1;i<=ntr_glob;i++){
		for(j=1;j<=nsrc_glob;j++){
			fscanf(ftracekill,"%d",&kill_tmp[i][j]);
		}
	}
78
	
Tilman Steinweg's avatar
Tilman Steinweg committed
79
	fclose(ftracekill);
80
	
Tilman Steinweg's avatar
Tilman Steinweg committed
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
	h=1;
	for(i=1;i<=ntr;i++){
	   kill_vector[h] = kill_tmp[recpos_loc[3][i]][ishot];
	   h++;
	}
} /* end if(TRKILL)*/

RMS=0.0;
Lcount=1;  

/* Integration of measured and synthetic data  */
for(i=1;i<=ntr;i++){
	intseis_s=0;
	intseis_sd=0;
	if (VELOCITY==0){	/* only integtration if displacement seismograms are inverted */
		for(j=1;j<=ns;j++){
			intseis_s  += section[i][j];
			intseis_section[i][j]     = intseis_s*DT;
			intseis_sd += sectiondata[i][j];
			intseis_sectiondata[i][j] = intseis_sd*DT;
		}
	}else{
		for(j=1;j<=ns;j++){
			intseis_section[i][j]     = section[i][j];
			intseis_sectiondata[i][j] = sectiondata[i][j];
		}
	}
}

/* TIME WINDOWING */
if(TIMEWIN==1){
112
113
	time_window(intseis_section, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
	time_window(intseis_sectiondata, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
Tilman Steinweg's avatar
Tilman Steinweg committed
114
115
116
117
}

/* NORMALIZE TRACES */
if(NORMALIZE==1){
118
119
	normalize_data(intseis_section,ntr,ns);
	normalize_data(intseis_sectiondata,ntr,ns);
Tilman Steinweg's avatar
Tilman Steinweg committed
120
121
122
123
124
125
126
}

/* calculate RMS */
/*for(i=1;i<=ntr;i++){
      for(j=1;j<=ns;j++){*/
      /*RMS += (sectionpdata[i][j]-sectionvx[i][j]) * (sectionpdata[i][j]-sectionvx[i][j]);*/
      /*RMS += (sectionpdata[i][j]) * (sectionpdata[i][j]);*/
127
      /*      Lcount++;
Tilman Steinweg's avatar
Tilman Steinweg committed
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
      }
} */

/*RMS=sqrt(RMS/Lcount);*/
RMS=1.0;

/* calculate weighted data residuals and reverse time direction */
/* calculate kind of "energy" */

/* calculate envelope and hilbert transform for LNORM==8 */	
if (LNORM==8){
	calc_envelope(intseis_sectiondata,intseis_sectiondata_envelope,ns,ntr);
	calc_envelope(intseis_section,intseis_section_envelope,ns,ntr);	
	calc_hilbert(intseis_section,intseis_section_hilbert,ns,ntr);
	for(i=1;i<=ntr;i++){
		for(j=1;j<=ns;j++){
			dummy_1[i][j] = (intseis_section_envelope[i][j] - intseis_sectiondata_envelope[i][j]) / (intseis_section_envelope[i][j]+WATERLEVEL_LNORM8) * intseis_section_hilbert[i][j];
			}
	}
	calc_hilbert(dummy_1,dummy_2,ns,ntr);
}
/* end of calculate envelope and hilbert transform for LNORM==8 */

151
152
153
for(i=1;i<=ntr;i++){
	if((TRKILL==1)&&(kill_vector[i]==1))
	continue;
Tilman Steinweg's avatar
Tilman Steinweg committed
154
155
	
	if ((LNORM==5) || (LNORM==7)){
156
		abs_sectiondata=0.0;
Tilman Steinweg's avatar
Tilman Steinweg committed
157
158
159
160
161
162
163
164
165
166
167
168
169
		abs_section=0.0;
		sectiondata_mult_section=0.0;
	
		for(j=1;j<=ns;j++){
			abs_sectiondata+=intseis_sectiondata[i][j]*intseis_sectiondata[i][j];
			abs_section+=intseis_section[i][j]*intseis_section[i][j];
			sectiondata_mult_section+=intseis_sectiondata[i][j]*intseis_section[i][j]; /* calculation of dot product for measured (section) and synthetic (sectiondata) data*/
		}
		abs_sectiondata=sqrt(abs_sectiondata);
		abs_section=sqrt(abs_section);
	}
	/* calculate residual seismograms and norm */	
	
170
171
	if((TRKILL==1)&&(kill_vector[i]==1))
	continue;
Tilman Steinweg's avatar
Tilman Steinweg committed
172
173
174
175
176

	/*reverse time direction */
	invtime=ns;	
	
	for(j=1;j<=ns;j++){
177
178
179
		/*printf("%d \t %d \t %e \t %e \n",i,j,sectionpdata[i][j],sectionp[i][j]);*/
		/* calculate L1 residuals */
		if(LNORM==1){
Tilman Steinweg's avatar
Tilman Steinweg committed
180
181
182
183
184
			if(((sectiondata[i][j]-section[i][j])/RMS)>0){signL1=1.0;}
			if(((sectiondata[i][j]-section[i][j])/RMS)<0){signL1=-1.0;}
			if(((sectiondata[i][j]-section[i][j])/RMS)==0){signL1=0.0;}
			
			sectiondiff[i][invtime]=signL1/RMS;
185
186
187
188
		}
		
		/* calculate L2 residuals */
		if(LNORM==2){
Tilman Steinweg's avatar
Tilman Steinweg committed
189
			sectiondiff[i][invtime] = intseis_section[i][j]-intseis_sectiondata[i][j];
190
191
192
193
		}
		
		/* calculate Cauchy residuals */  /* NOT UP TO DATE */
		if(LNORM==3){
Tilman Steinweg's avatar
Tilman Steinweg committed
194
			sectiondiff[i][invtime]=((sectiondata[i][j]-section[i][j])/RMS)/(1+(((sectiondata[i][j]-section[i][j])/RMS)*((sectiondata[i][j]-section[i][j])/RMS)))/RMS; 
195
196
197
198
		}
		
		/* calculate sech residuals */   /* NOT UP TO DATE */
		if(LNORM==4){
Tilman Steinweg's avatar
Tilman Steinweg committed
199
			sectiondiff[i][invtime]=(tanh((sectiondata[i][j]-section[i][j])/RMS))/RMS;
200
201
202
203
		}
		
		/* calculate LNORM 5 or LNORM 7 residuals */
		if((LNORM==5) || (LNORM==7)){
Tilman Steinweg's avatar
Tilman Steinweg committed
204
			sectiondiff[i][invtime]=((intseis_section[i][j]*sectiondata_mult_section)/(abs_section*abs_section*abs_section*abs_sectiondata)) - (intseis_sectiondata[i][j]/(abs_section*abs_sectiondata));
205
206
207
208
		}
		
		/* calculate LNORM 8 residuals */
		if (LNORM==8){
Tilman Steinweg's avatar
Tilman Steinweg committed
209
			sectiondiff[i][invtime]=intseis_section[i][j]/(intseis_section_envelope[i][j]+WATERLEVEL_LNORM8)*(intseis_section_envelope[i][j]-intseis_sectiondata_envelope[i][j])-dummy_2[i][j];
210
211
212
213
214
215
		}
		
		/*sectionpdiff[i][invtime]=sectionp[i][j];*/
		
		/* calculate norm for different LNORM */
		if((LNORM==2) && (swstestshot==1)){
Tilman Steinweg's avatar
Tilman Steinweg committed
216
			L2+=sectiondiff[i][invtime]*sectiondiff[i][invtime];
217
218
219
		}
		
		if(((LNORM==5)||(LNORM==7)) && (swstestshot==1)){
Tilman Steinweg's avatar
Tilman Steinweg committed
220
			L2-=(intseis_sectiondata[i][j]*intseis_section[i][j])/(abs_sectiondata*abs_section);
221
222
223
		}
		
		if((LNORM==8) && (swstestshot==1)){
Tilman Steinweg's avatar
Tilman Steinweg committed
224
			L2+=(intseis_section_envelope[i][invtime]-intseis_sectiondata_envelope[i][invtime]) * (intseis_section_envelope[i][invtime]-intseis_sectiondata_envelope[i][invtime]);
225
226
227
		}
		
		if((sws==2)&&(swstestshot==1)){
Tilman Steinweg's avatar
Tilman Steinweg committed
228
			L2+=fabs(sectiondiff[i][invtime])*fabs(sectiondiffold[i][j]);
229
230
231
232
233
234
		}
			
		/*L2+=sectiondiff[i][invtime];*/
		
		invtime--;          /* reverse time direction */
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
235
236
237
238
239
240
241
242
243
}

l2=L2;

/* free memory for integrated seismograms */
free_matrix(intseis_section,1,ntr,1,ns);
free_matrix(intseis_sectiondata,1,ntr,1,ns);

if(TRKILL==1){
244
245
246
247
	free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob);
	free_ivector(kill_vector,1,ntr);
}

Tilman Steinweg's avatar
Tilman Steinweg committed
248
249
250
251
252
253
254
255
256
if(LNORM==8){
	free_matrix(intseis_sectiondata_envelope,1,ntr,1,ns);
	free_matrix(intseis_section_envelope,1,ntr,1,ns);
	free_matrix(intseis_section_hilbert,1,ntr,1,ns);
	free_matrix(dummy_1,1,ntr,1,ns);
	free_matrix(dummy_2,1,ntr,1,ns);
	}
return l2;
} /* end of function */