time_window_glob.c 5.97 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 /*-----------------------------------------------------------------------------------------
 * 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>.
-----------------------------------------------------------------------------------------*/

/*------------------------------------------------------------------------
20
 *   Apply time damping (after Brossier (2009)
21
22
23
24
25
26
27
28
 *  ----------------------------------------------------------------------*/
#include "fd.h"

void time_window_glob(float **sectiondata, int iter, int ntr_glob, int ns, int ishot){

	/* declaration of variables */
	extern float DT;
	extern float GAMMA, TWLENGTH_PLUS, TWLENGTH_MINUS;
29
	extern int TW_IND, USE_WORKFLOW, WORKFLOW_STAGE;
30
31
	extern char PICKS_FILE[STRING_SIZE];
	char pickfile_char[STRING_SIZE];
32
33
34
	float time, dumpa, dumpb, dumpc, dumpd;
	float dump1, dump2, dump3, dump4, dump5, dump6;
	float taper, taper1, taper2, taper3;
35
36
37
38
	int i, j;
	
	float **picked_times_m=NULL;
	float *picked_times=NULL;
39
    float **dummysection=NULL;
40
41
42
	
	FILE *fptime;
	
43
	if(TW_IND==1){
44
		picked_times_m = matrix(1,3,1,ntr_glob);
45
46
47
48
    }else if(TW_IND==2){
        picked_times_m = matrix(1,6,1,ntr_glob);
        dummysection = matrix(1,ntr_glob,1,ns);
    }else{
49
		picked_times = vector(1,ntr_glob);
50
51
    }
    
52
	/* read picked first arrival times */
53
54
55
56
57
58
59
    if(USE_WORKFLOW){
        sprintf(pickfile_char,"%s_%i_%i.dat",PICKS_FILE,ishot,WORKFLOW_STAGE);
        fptime=fopen(pickfile_char,"r");
        if (fptime == NULL) {
            sprintf(pickfile_char,"%s_%i.dat",PICKS_FILE,ishot);
            fptime=fopen(pickfile_char,"r");
            if (fptime == NULL) {
Florian Wittkamp's avatar
Florian Wittkamp committed
60
                declare_error(" picks_?.dat could not be opened !");
61
62
63
64
65
66
            }
        }
    }else{
        sprintf(pickfile_char,"%s_%i.dat",PICKS_FILE,ishot);
        fptime=fopen(pickfile_char,"r");
        if (fptime == NULL) {
Florian Wittkamp's avatar
Florian Wittkamp committed
67
            declare_error(" picks_?.dat could not be opened !");
68
69
        }
    }
70
	
71
	if(TW_IND==1){
72
		for(i=1;i<=ntr_glob;i++){
73
74
			fscanf(fptime,"%f%f%f",&dump1,&dump2,&dump3);
			picked_times_m[1][i] = dump1;
75
76
77
			picked_times_m[2][i] = dump2;
			picked_times_m[3][i] = dump3;
		}
78
79
80
81
82
83
84
85
86
87
88
	}if(TW_IND==2){
        for(i=1;i<=ntr_glob;i++){
            fscanf(fptime,"%f%f%f%f%f%f",&dump1,&dump2,&dump3,&dump4,&dump5,&dump6);
            picked_times_m[1][i] = dump1;
            picked_times_m[2][i] = dump2;
            picked_times_m[3][i] = dump3;
            picked_times_m[4][i] = dump4;
            picked_times_m[5][i] = dump5;
            picked_times_m[6][i] = dump6;
        }
    }else{
89
		for(i=1;i<=ntr_glob;i++){
90
91
			fscanf(fptime,"%f",&dump1);
			picked_times[i] = dump1;
92
93
94
95
96
		}
	}
	
	fclose(fptime);
	
97
	if(TW_IND==1){
98
99
100
101
102
		for(i=1;i<=ntr_glob;i++){
		for(j=2;j<=ns;j++){
		
			time = (float)(j * DT);
			
103
104
			dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
			taper = exp(-GAMMA*dumpa*dumpa);
105
			
106
107
			dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); 
			taper1 = exp(-GAMMA*dumpb*dumpb);
108
109
110
111
112
113
114
115
116
117
118
			
			if(time>=picked_times_m[1][i]+picked_times_m[2][i]){
			sectiondata[i][j] = sectiondata[i][j] * taper;}
			
			if(time<=picked_times_m[1][i]-picked_times_m[3][i]){
			sectiondata[i][j] = sectiondata[i][j] * taper1;}
			
			sectiondata[i][j] = sectiondata[i][j];
			
		}
		}
laura.gassner's avatar
laura.gassner committed
119
	}else if(TW_IND==2){
120
121
122
123
124
125
126
127
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
        for(i=1;i<=ntr_glob;i++){
        for(j=2;j<=ns;j++){
        
            time = (float)(j * DT);
            
            dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
            taper = exp(-GAMMA*dumpa*dumpa);
            
            dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); 
            taper1 = exp(-GAMMA*dumpb*dumpb);
            
            dummysection[i][j] = sectiondata[i][j];
            
            if(time>=picked_times_m[1][i]+picked_times_m[2][i]){
            dummysection[i][j] = sectiondata[i][j] * taper;}
            
            if(time<=picked_times_m[1][i]-picked_times_m[3][i]){
            dummysection[i][j] = sectiondata[i][j] * taper1;}
                        
            dumpc = (time-picked_times_m[4][i]-picked_times_m[5][i]);
            taper2 = exp(-GAMMA*dumpc*dumpc);
            
            dumpd = (time-picked_times_m[4][i]+picked_times_m[6][i]); 
            taper3 = exp(-GAMMA*dumpd*dumpd);
            
            if(time>=picked_times_m[4][i]+picked_times_m[5][i]){
            sectiondata[i][j] = sectiondata[i][j] * taper2;}
            
            if(time<=picked_times_m[4][i]-picked_times_m[6][i]){
            sectiondata[i][j] = sectiondata[i][j] * taper3;}
            
            sectiondata[i][j] = sectiondata[i][j] + dummysection[i][j];
            
        }
        }
    }else{
156
157
158
159
160
		for(i=1;i<=ntr_glob;i++){
		for(j=2;j<=ns;j++){
		
			time = (float)(j * DT);
			
161
162
			dumpa = (time-picked_times[i]-TWLENGTH_PLUS);
			taper = exp(-GAMMA*dumpa*dumpa);
163
			
164
165
			dumpb = (time-picked_times[i]+TWLENGTH_MINUS); 
			taper1 = exp(-GAMMA*dumpb*dumpb);
166
167
168
169
170
171
172
173
174
175
176
177
			
			if(time>=picked_times[i]+TWLENGTH_PLUS){
			sectiondata[i][j] = sectiondata[i][j] * taper;}
			
			if(time<=picked_times[i]-TWLENGTH_MINUS){
			sectiondata[i][j] = sectiondata[i][j] * taper1;}
			
			sectiondata[i][j] = sectiondata[i][j];
		}
		}
	}
	
178
	if(TW_IND==1){
179
		free_matrix(picked_times_m,1,3,1,ntr_glob);
180
181
182
183
    }else if(TW_IND==2){
        free_matrix(picked_times_m,1,6,1,ntr_glob);
        free_matrix(dummysection,1,ntr_glob,1,ns);
    }else{
184
		free_vector(picked_times,1,ntr_glob);
185
186
187
    }
    
} /* end of function time_window_glob.c */