stf.c 9.64 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
/*-----------------------------------------------------------------------------------------
 * 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>.
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   inversion for source time function 
 *   31. August 2011 L. Rehor, T. Forbriger, M. Schaefer
 *  ----------------------------------------------------------------------*/

#include "fd.h"
#include "stfinv/stfinv.h"
#include "segy.h"

void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy_conv, float * source_time_function, int  **recpos, int  **recpos_loc, 
int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots, float FC){ 

	/* declaration of global variables */
	extern float DT, DH;
33
34
	extern int SEIS_FORMAT, MYID, NT, QUELLART, TIME_FILT, TIMEWIN;
	extern char  SEIS_FILE_VY[STRING_SIZE], SEIS_FILE_P[STRING_SIZE], PARA[STRING_SIZE], DATA_DIR[STRING_SIZE];
Tilman Steinweg's avatar
Tilman Steinweg committed
35
36
37
38
39
40
41
42
43
	extern int TRKILL_STF, NORMALIZE;
	extern char TRKILL_FILE_STF[STRING_SIZE];
	extern char SIGNAL_FILE[STRING_SIZE];
	
	/* declaration of variables for trace killing */
	int ** kill_tmp, *kill_vector, h, j;
	char trace_kill_file[STRING_SIZE];	
	FILE *ftracekill;
	
44
45
	float *picked_times=NULL;
	
Tilman Steinweg's avatar
Tilman Steinweg committed
46
47
48
49
50
	/* --------------- declaration of variables --------------- */
	unsigned int nrec, nsamp, i, npairs;
	float dt;
	float xr=0.0, yr=0.0;
	float XS=0.0, YS=0.0;
51
	char conv_y[STRING_SIZE], qw[STRING_SIZE], conv_y_tmp[STRING_SIZE], obs_y_tmp[STRING_SIZE], mod_y_tmp[STRING_SIZE];
Tilman Steinweg's avatar
Tilman Steinweg committed
52
53
54
55
56
	
	/* variables for wavelet */
	int nt, nts;
	float tshift, amp=0.0, fc, tau, t, ts, ag;
	float * wavelet, * stf_conv_wavelet, *psource=NULL;
57
	
Tilman Steinweg's avatar
Tilman Steinweg committed
58
59
60
	wavelet=vector(1,ns);
	stf_conv_wavelet=vector(1,ns);
	
61
62
	if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
	
Tilman Steinweg's avatar
Tilman Steinweg committed
63
64
65
66
67
	printf("\n================================================================================================\n\n");
	printf("\n ***** Inversion of Source Time Function - shot: %d - it: %d ***** \n\n",ishot,iter);
	
	/* if TRKILL_STF==1 a trace killing is applied */
	if(TRKILL_STF){
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
		kill_tmp = imatrix(1,ntr_glob,1,nshots);
		kill_vector = ivector(1,ntr_glob);
		
		ftracekill=fopen(TRKILL_FILE_STF,"r");
		
		if (ftracekill==NULL) err(" Trace kill file could not be opened!");
		
		for(i=1;i<=ntr_glob;i++){
			for(j=1;j<=nshots;j++){
				fscanf(ftracekill,"%d",&kill_tmp[i][j]);
			}
		}
		
		fclose(ftracekill);
		
		h=1;
		for(i=1;i<=ntr_glob;i++){
		kill_vector[h] = kill_tmp[i][ishot];
		h++;
Tilman Steinweg's avatar
Tilman Steinweg committed
87
88
89
90
91
92
93
94
95
96
97
98
99
		}
	} /* end if(TRKILL_STF)*/	
	
	if(TRKILL_STF){
		for(i=1;i<=ntr_glob;i++){
		if(i==1)printf("\n ***** Trace killing is applied for trace: ***** \n ***** \t");
			if(kill_vector[i]==1){
				printf("%d \t",i);
				for(j=1;j<=ns;j++){
				sectionvy[i][j]=0.0;
				sectionvy_obs[i][j]=0.0;
						}
				}	
100
		if(i==ntr_glob)printf(" ***** \n\n");
Tilman Steinweg's avatar
Tilman Steinweg committed
101
		}
102
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
103
104
	/* trace killing ends here */
	
105
106
107
108
109
	if(TIMEWIN==1){
		time_window(sectionvy, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
		time_window(sectionvy_obs, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
	}
	
Tilman Steinweg's avatar
Tilman Steinweg committed
110
111
112
113
	/* NORMALIZE TRACES */
	/*if(NORMALIZE==1){*/
	normalize_data(sectionvy,ntr_glob,ns);
	normalize_data(sectionvy_obs,ntr_glob,ns);
114
	/*}*/
Tilman Steinweg's avatar
Tilman Steinweg committed
115
116
117
118
119
120
121
122
123
124
	
	nrec=(unsigned int)ntr_glob;
	nsamp=(unsigned int)ns;
	dt=DT;
	npairs=nrec;
	
	/* source coordinates are written into trace header fields */
	XS=srcpos[1][ishot];	
	YS=srcpos[2][ishot];
	
125
	
Tilman Steinweg's avatar
Tilman Steinweg committed
126
	/* TF Software: see libstfinv */
127
	
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
151
152
153
154
155
156
157
158
159
	struct CTriples data;
	data.n=nrec;
	data.triples=(struct CWaveformTriple *)malloc(nrec*sizeof(struct CWaveformTriple));
	if (data.triples == NULL) {abort();}
	for (i=0;i<nrec;i++){
	
		xr=recpos[1][i+1]*DH;
		yr=recpos[2][i+1]*DH;
		
		data.triples[i].data=&sectionvy_obs[i+1][1];
		
		data.triples[i].synthetics=&sectionvy[i+1][1];
		
		data.triples[i].convolvedsynthetics=&sectionvy_conv[i+1][1];
		
		data.triples[i].header.sx=(unsigned int)iround(XS*1000.0);  /* X source coordinate */
		data.triples[i].header.sy=0.0;
		data.triples[i].header.sz=(unsigned int)iround(YS*1000.0);  /* source depth (positive) */
		data.triples[i].header.rx=(unsigned int)iround(xr*1000.0);  /* group coordinates */
		data.triples[i].header.ry=0.0;
		data.triples[i].header.rz=(unsigned int)iround(yr*1000.0);
		data.triples[i].header.sampling.n=nsamp;
		data.triples[i].header.sampling.dt=dt;
	}
	
	struct CWaveform stf;
	stf.series = &source_time_function[1];
	stf.sampling.n=nsamp;
	stf.sampling.dt=dt;
	
	
	initstfinvengine(data, stf, PARA);
160
	
Tilman Steinweg's avatar
Tilman Steinweg committed
161
162
163
164
165
	runstfinvengine();
	
	/* END TF Software */
	
	
166
	psource=vector(1,ns);
Tilman Steinweg's avatar
Tilman Steinweg committed
167
	
168
169
170
171
	if (QUELLART==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r"));
	if (QUELLART==7){
		inseis_source_wavelet(psource,ns,ishot);
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
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
	
	/* calculating wavelet SIN**3 for convoling with STF */
	tshift=srcpos[4][ishot];
	fc=srcpos[5][ishot];
	ts=1.0/fc;
	for (nt=1;nt<=ns;nt++){
		t=(float)nt*DT;
		switch (QUELLART){
					case 1 : 
						/* Old Ricker Wavelet */
						/* tau=PI*(t-ts-tshift)/(1.5*ts);
						amp=(((1.0-4.0*tau*tau)*exp(-2.0*tau*tau))); */ 

						/* New Ricker Wavelet, equal to SOFI2D */
						tau=PI*(t-1.5*ts-tshift)/(ts);
						amp=(((1.0-2.0*tau*tau)*exp(-tau*tau)));
					break;
					case 2 : 
						if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
						else amp=((sin(2.0*PI*(t-tshift)*fc) 
			    				-0.5*sin(4.0*PI*(t-tshift)*fc)));

						/*amp=((sin(2.0*PI*(t+tshift)*fc) 
			    			-0.5*sin(4.0*PI*(t+tshift)*fc)));*/
					break;
					case 3 : 
						if (nt<=nts) amp=psource[nt]; 
						else amp=0.0;
// 						amp=psource[nt];
					break;  /* source wavelet from file SOURCE_FILE */
					case 4 : 
						/*tau=PI*(t-ts-tshift)/(1.5*ts);*/ /* Ricker */
						/*amp=((t-ts-tshift)*exp(-2.0*tau*tau));*/
						
						if ((t<tshift) || (t>(tshift+ts))) amp=0.0;
						else amp=pow(sin(PI*(t+tshift)/ts),3.0);
						break; /* sinus raised to the power of three */
						
					break; 
					case 5 : 
				                /* first derivative of a Gaussian */
				                ts=1.2/fc;
					        ag  = PI*PI*fc*fc;
			                        amp = - 2.0 * ag * (t-ts) * exp(-ag*(t-ts)*(t-ts));
					break;					
					case 6 : 
					        /* Bandlimited Spike */
						amp=0.0;
						if(nt==1+iround(tshift/DT)){
						amp = 1.0;}
					break; 
					case 7 :
					        /* source wavelet from file SOURCE_FILE */ 
						amp=psource[nt];
					break;  
					case 8 :
					        /* integral of sinus raised to the power of three */
						if (t<tshift) {
							 amp=0.0;}
						if ((t>=tshift) && (t<=(tshift+ts))){
							amp=(ts/(0.75*PI))*(0.5-0.75*cos(PI*(t-tshift)/ts)+0.25*pow(cos(PI*(t-tshift)/ts),3.0));}
						if (t>(tshift+ts))
							{amp=ts/(0.75*PI);}
						break;                                                                                                                                           	
					default : 
						err("Which source-wavelet ? ");
					
					
		}/* end of switch (QUELLART)  */
		wavelet[nt]=amp;
		
	}/*  end of for (nt=1;nt<=ns;nt++) */	
244
245
	
	
Tilman Steinweg's avatar
Tilman Steinweg committed
246
247
248
249
	/* convolving wavelet with STF */
	conv_FD(wavelet,source_time_function,stf_conv_wavelet,ns);
	
	
250
251
252
253
254
255
256
257
258
259
260
// 	/* --------------- writing out the observed seismograms --------------- */
// 	sprintf(obs_y_tmp,"%s.shot%d.obs",SEIS_FILE_P,ishot);
// 	printf(" PE %d is writing %d observed seismograms (vy) for shot = %d to\n\t %s \n",MYID,ntr_glob,ishot,obs_y_tmp);
// 	outseis_glob(fp,fopen(obs_y_tmp,"w"),1,sectionvy_obs,recpos,recpos_loc,ntr_glob,srcpos,0,ns,SEIS_FORMAT,ishot,0);
// 	
// 	/* --------------- writing out the modelled seismograms --------------- */
// 	sprintf(mod_y_tmp,"%s.shot%d.mod",SEIS_FILE_P,ishot);
// 	printf(" PE %d is writing %d modelled seismograms (vy) for shot = %d to\n\t %s \n",MYID,ntr_glob,ishot,mod_y_tmp);
// 	outseis_glob(fp,fopen(mod_y_tmp,"w"),1,sectionvy,recpos,recpos_loc,ntr_glob,srcpos,0,ns,SEIS_FORMAT,ishot,0);
	
	
Tilman Steinweg's avatar
Tilman Steinweg committed
261
	/* --------------- writing out the source time function --------------- */
262
263
264
265
	if((TIME_FILT==1)||(TIME_FILT==2)){
		sprintf(qw,"%s.shot%d_%dHz",SIGNAL_FILE,ishot,(int)FC);
		printf(" PE %d is writing source time function for shot = %d to\n\t %s \n",MYID,ishot,qw);
		outseis_vector(fp,fopen(qw,"w"),1,stf_conv_wavelet,recpos,recpos_loc,ntr,srcpos,0,ns,SEIS_FORMAT,ishot,0);
Tilman Steinweg's avatar
Tilman Steinweg committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
	}
	
	sprintf(qw,"%s.shot%d",SIGNAL_FILE,ishot);
	printf(" PE %d is writing source time function for shot = %d to\n\t %s \n",MYID,ishot,qw);
	outseis_vector(fp,fopen(qw,"w"),1,stf_conv_wavelet,recpos,recpos_loc,ntr,srcpos,0,ns,SEIS_FORMAT,ishot,0);
	
	/*sprintf(conv_y_tmp,"%s.shot%d.forward",SEIS_FILE_VY,ishot);
	printf(" PE %d is writing %d seismograms (vy) for shot = %d to\n\t %s \n",MYID,ntr_glob,ishot,conv_y_tmp);
	outseis_glob(fp,fopen(conv_y_tmp,"w"),1,sectionvy,recpos,recpos_loc,ntr_glob,srcpos,0,ns,SEIS_FORMAT,ishot,0);*/
	
	/*freestfinvengine();
	free(data.triples);*/
	
	free_vector(wavelet,1,ns);
	free_vector(stf_conv_wavelet,1,ns);
	free_vector(psource,1,ns);
	
283
284
285
	/* free memory for time windowing and trace killing */
	if(TIMEWIN==1) free_vector(picked_times,1,ntr);

Tilman Steinweg's avatar
Tilman Steinweg committed
286
287
	/* free memory for trace killing */
	if(TRKILL_STF){
288
289
		free_imatrix(kill_tmp,1,ntr_glob,1,nshots);
		free_ivector(kill_vector,1,ntr_glob);
Tilman Steinweg's avatar
Tilman Steinweg committed
290
	}
291
}