time_window_glob.c 5.93 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;
		}
laura.gassner's avatar
laura.gassner committed
78
	}else if(TW_IND==2){
79 80 81 82 83 84 85 86 87 88
        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
			dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
104
			taper = exp(-GAMMA*dumpa);
105
			
106
			dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); 
107
			taper1 = exp(-GAMMA*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
        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]);
126
            taper = exp(-GAMMA*dumpa);
127 128
            
            dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); 
129
            taper1 = exp(-GAMMA*dumpb);
130 131 132 133 134 135 136 137 138 139
            
            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]);
140
            taper2 = exp(-GAMMA*dumpc);
141 142
            
            dumpd = (time-picked_times_m[4][i]+picked_times_m[6][i]); 
143
            taper3 = exp(-GAMMA*dumpd);
144 145 146 147 148 149 150 151 152 153 154 155
            
            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
			dumpa = (time-picked_times[i]-TWLENGTH_PLUS);
162
			taper = exp(-GAMMA*dumpa);
163
			
164
			dumpb = (time-picked_times[i]+TWLENGTH_MINUS); 
165
			taper1 = exp(-GAMMA*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 */