stf.c 12.2 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
/*-----------------------------------------------------------------------------------------
Florian Wittkamp's avatar
Florian Wittkamp committed
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
Tilman Steinweg's avatar
Tilman Steinweg committed
3
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
4
 * This file is part of IFOS.
Tilman Steinweg's avatar
Tilman Steinweg committed
5
 * 
Florian Wittkamp's avatar
Florian Wittkamp committed
6
 * IFOS is free software: you can redistribute it and/or modify
Tilman Steinweg's avatar
Tilman Steinweg committed
7 8 9
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
 * 
Florian Wittkamp's avatar
Florian Wittkamp committed
10
 * IFOS is distributed in the hope that it will be useful,
Tilman Steinweg's avatar
Tilman Steinweg committed
11 12 13 14 15
 * 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
Florian Wittkamp's avatar
Florian Wittkamp committed
16
 * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
Tilman Steinweg's avatar
Tilman Steinweg committed
17 18 19 20 21 22 23 24 25 26
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
 *   inversion for source time function 
 *  ----------------------------------------------------------------------*/

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

27
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, int SH,int nsrc_glob){
Tilman Steinweg's avatar
Tilman Steinweg committed
28 29 30

	/* declaration of global variables */
	extern float DT, DH;
31
	extern int SEIS_FORMAT, MYID, NT, SOURCE_SHAPE, TIME_FILT, TIMEWIN, TAPER_STF;
32
	extern char  PARA[STRING_SIZE], DATA_DIR[STRING_SIZE];
33
	extern int TRKILL_STF, NORMALIZE, USE_WORKFLOW, WORKFLOW_STAGE;
Tilman Steinweg's avatar
Tilman Steinweg committed
34 35
	extern char TRKILL_FILE_STF[STRING_SIZE];
	extern char SIGNAL_FILE[STRING_SIZE];
Florian Wittkamp's avatar
Florian Wittkamp committed
36 37
    extern char SIGNAL_FILE_SH[STRING_SIZE];
    
38 39 40 41 42 43 44
    extern int TRKILL_STF_OFFSET;
    extern float TRKILL_STF_OFFSET_LOWER;
    extern float TRKILL_STF_OFFSET_UPPER;
    
    extern int USE_WORKFLOW;
    extern int WORKFLOW_STAGE;
    
Tilman Steinweg's avatar
Tilman Steinweg committed
45
	/* declaration of variables for trace killing */
Florian Wittkamp's avatar
Florian Wittkamp committed
46
	int ** kill_tmp = NULL, *kill_vector = NULL, h, j;
47
	char trace_kill_file[STRING_SIZE];
Tilman Steinweg's avatar
Tilman Steinweg committed
48 49 50 51 52 53 54
	FILE *ftracekill;
	
	/* --------------- 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;
Florian Wittkamp's avatar
Florian Wittkamp committed
55
	char qw[STRING_SIZE];
Tilman Steinweg's avatar
Tilman Steinweg committed
56 57
	
	/* variables for wavelet */
Florian Wittkamp's avatar
Florian Wittkamp committed
58
	int nt, nts = 0;
Tilman Steinweg's avatar
Tilman Steinweg committed
59 60
	float tshift, amp=0.0, fc, tau, t, ts, ag;
	float * wavelet, * stf_conv_wavelet, *psource=NULL;
61
	
Tilman Steinweg's avatar
Tilman Steinweg committed
62 63 64 65 66 67
	wavelet=vector(1,ns);
	stf_conv_wavelet=vector(1,ns);
	
	printf("\n================================================================================================\n\n");
	printf("\n ***** Inversion of Source Time Function - shot: %d - it: %d ***** \n\n",ishot,iter);
	
68 69 70 71 72
    /* if TRKILL_STF==1 a trace killing is applied */
    if(TRKILL_STF){
        kill_tmp = imatrix(1,ntr_glob,1,nshots);
        kill_vector = ivector(1,ntr_glob);
        
73 74 75 76 77 78 79 80 81 82 83
        
        /*------------------*/
        /* clear kill table */
        /*------------------*/
        for(i=1;i<=nsrc_glob;i++){
            for (j=1; j<=ntr_glob; j++) {
                kill_tmp[j][i]=0;
            }
        }
        
        if(TRKILL_STF_OFFSET==1) {
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
            
            if(MYID==0) {
                printf("Automatic offset based TraceKill for STF\n");
            }
            
            /* Generate TraceKill file on the fly */
            create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,-100,TRKILL_STF_OFFSET_LOWER,TRKILL_STF_OFFSET_UPPER);
            
        } else {
            if(USE_WORKFLOW){
                sprintf(trace_kill_file,"%s_%i.dat",TRKILL_FILE_STF,WORKFLOW_STAGE);
                ftracekill=fopen(trace_kill_file,"r");
                if (ftracekill==NULL){
                    sprintf(trace_kill_file,"%s.dat",TRKILL_FILE_STF);
                    ftracekill=fopen(trace_kill_file,"r");
                    if (ftracekill==NULL){
                        declare_error(" Trace kill file could not be opened!");
                    }
                }
            }else{
104 105 106
                sprintf(trace_kill_file,"%s.dat",TRKILL_FILE_STF);
                ftracekill=fopen(trace_kill_file,"r");
                if (ftracekill==NULL){
Florian Wittkamp's avatar
Florian Wittkamp committed
107
                    declare_error(" Trace kill file could not be opened!");
108 109
                }
            }
110 111 112 113 114
            
            for(i=1;i<=ntr_glob;i++){
                for(j=1;j<=nshots;j++){
                    fscanf(ftracekill,"%d",&kill_tmp[i][j]);
                }
115
            }
116 117
            
            fclose(ftracekill);
118 119 120 121 122
            
            if(TRKILL_STF_OFFSET==2) {
                /* Generate TraceKill file on the fly */
                create_trkill_table(kill_tmp,ntr_glob,recpos,nsrc_glob,srcpos,-100,TRKILL_STF_OFFSET_LOWER,TRKILL_STF_OFFSET_UPPER);
            }
123
        }
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
        h=1;
        for(i=1;i<=ntr_glob;i++){
            kill_vector[h] = kill_tmp[i][ishot];
            h++;
        }
    } /* 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;
                }
            }
            if(i==ntr_glob)printf(" ***** \n\n");
        }
    }
Tilman Steinweg's avatar
Tilman Steinweg committed
144 145
	/* trace killing ends here */
	
146
	if(TIMEWIN==1){
147 148
		time_window_glob(sectionvy, iter, ntr_glob, ns, ishot);
		time_window_glob(sectionvy_obs, iter, ntr_glob, ns, ishot);
149 150
	}
	
Tilman Steinweg's avatar
Tilman Steinweg committed
151
	/* NORMALIZE TRACES */
152 153 154 155
	if(NORMALIZE==1){
		normalize_data(sectionvy,ntr_glob,ns);
		normalize_data(sectionvy_obs,ntr_glob,ns);
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
156 157 158 159 160 161 162 163 164 165
	
	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];
	
166
	
Tilman Steinweg's avatar
Tilman Steinweg committed
167
	/* TF Software: see libstfinv */
168
	
Tilman Steinweg's avatar
Tilman Steinweg committed
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
	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);
201
	
Tilman Steinweg's avatar
Tilman Steinweg committed
202 203 204 205 206
	runstfinvengine();
	
	/* END TF Software */
	
	
207
	psource=vector(1,ns);
Tilman Steinweg's avatar
Tilman Steinweg committed
208
	
209 210
	if (SOURCE_SHAPE==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r"));
	if (SOURCE_SHAPE==7){
Florian Wittkamp's avatar
Florian Wittkamp committed
211
		inseis_source_wavelet(psource,ns,ishot,SH);
212
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
213 214 215 216 217 218 219
	
	/* 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;
220
		switch (SOURCE_SHAPE){
Tilman Steinweg's avatar
Tilman Steinweg committed
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
					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 : 
Florian Wittkamp's avatar
Florian Wittkamp committed
278
						declare_error("Which source-wavelet ? ");
Tilman Steinweg's avatar
Tilman Steinweg committed
279 280
					
					
281
		}/* end of switch (SOURCE_SHAPE)  */
Tilman Steinweg's avatar
Tilman Steinweg committed
282 283 284
		wavelet[nt]=amp;
		
	}/*  end of for (nt=1;nt<=ns;nt++) */	
285 286
	
	
Tilman Steinweg's avatar
Tilman Steinweg committed
287 288 289
	/* convolving wavelet with STF */
	conv_FD(wavelet,source_time_function,stf_conv_wavelet,ns);
	
290
	if(TAPER_STF)
291
		taper(stf_conv_wavelet, ns, fc);
Tilman Steinweg's avatar
Tilman Steinweg committed
292
	
293 294 295 296 297 298 299 300 301 302
// 	/* --------------- 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
303
	/* --------------- writing out the source time function --------------- */
304
	if((TIME_FILT==1)||(TIME_FILT==2)){
Florian Wittkamp's avatar
Florian Wittkamp committed
305
        if(SH==0) {
306 307 308 309 310
            if(USE_WORKFLOW){
                sprintf(qw,"%s.stage%d.shot%d_%dHz.su",SIGNAL_FILE,WORKFLOW_STAGE,ishot,(int)FC);
            } else {
                sprintf(qw,"%s.shot%d_%dHz.su",SIGNAL_FILE,ishot,(int)FC);
            }
Florian Wittkamp's avatar
Florian Wittkamp committed
311
        } else {
312 313 314 315 316
            if(USE_WORKFLOW){
                sprintf(qw,"%s.stage%d.shot%d_%dHz.su",SIGNAL_FILE_SH,WORKFLOW_STAGE,ishot,(int)FC);
            } else {
                sprintf(qw,"%s.shot%d_%dHz.su",SIGNAL_FILE_SH,ishot,(int)FC);
            }
Florian Wittkamp's avatar
Florian Wittkamp committed
317
        }
318 319
		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
320
	}
Florian Wittkamp's avatar
Florian Wittkamp committed
321 322
    
    if(SH==0) {
323 324 325 326 327
        if(USE_WORKFLOW){
            sprintf(qw,"%s.stage%d.shot%d.su",SIGNAL_FILE,WORKFLOW_STAGE,ishot);
        } else {
            sprintf(qw,"%s.shot%d.su",SIGNAL_FILE,ishot);
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
328
    } else {
329 330 331 332 333
        if(USE_WORKFLOW){
            sprintf(qw,"%s.stage%d.shot%d.su",SIGNAL_FILE_SH,WORKFLOW_STAGE,ishot);
        } else {
            sprintf(qw,"%s.shot%d.su",SIGNAL_FILE_SH,ishot);
        }
Florian Wittkamp's avatar
Florian Wittkamp committed
334
    }
Tilman Steinweg's avatar
Tilman Steinweg committed
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
	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);
	
	/* free memory for trace killing */
	if(TRKILL_STF){
351 352
		free_imatrix(kill_tmp,1,ntr_glob,1,nshots);
		free_ivector(kill_vector,1,ntr_glob);
Tilman Steinweg's avatar
Tilman Steinweg committed
353
	}
354
}