time_window.c 7.7 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
/*-----------------------------------------------------------------------------------------
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
 *   Apply time damping (after Brossier (2009))                                 
Tilman Steinweg's avatar
Tilman Steinweg committed
21 22 23
 *  ----------------------------------------------------------------------*/
#include "fd.h"

24
void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc, int ntr, int ns, int ishot){
Tilman Steinweg's avatar
Tilman Steinweg committed
25

26 27 28
    /* 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 35
    float time, dumpa, dumpb, dumpc, dumpd;
    float dump1, dump2, dump3, dump4, dump5, dump6;
    float taper, taper1, taper2, taper3;
    float *pick_tmp, **pick_tmp_m, **dummysection=NULL;
36
    int i, j, h;
37

38 39 40
    float *picked_times=NULL, **picked_times_m=NULL;
    
    FILE *fptime;
41

42
    if(TW_IND==1){
43 44
        picked_times_m = matrix(1,3,1,ntr_glob);
        pick_tmp_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);
        pick_tmp_m = matrix(1,6,1,ntr_glob);
        dummysection = matrix(1,ntr,1,ns);
49 50 51 52 53
    }else{
        picked_times = vector(1,ntr_glob);
        pick_tmp = vector(1,ntr_glob);
    }
    
54 55 56 57 58 59 60 61
    /* read picked first arrival times */
    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) {
62
                declare_error(" picks_?.dat could not be opened !");
63 64 65 66 67 68
            }
        }
    }else{
        sprintf(pickfile_char,"%s_%i.dat",PICKS_FILE,ishot);
        fptime=fopen(pickfile_char,"r");
        if (fptime == NULL) {
69
            declare_error(" picks_?.dat could not be opened !");
70
        }
71 72
    }
    
73
    if(TW_IND==1){
74
        for(i=1;i<=ntr_glob;i++){
75 76
            fscanf(fptime,"%f%f%f",&dump1,&dump2,&dump3);
            pick_tmp_m[1][i] = dump1;
77 78 79
            pick_tmp_m[2][i] = dump2;
            pick_tmp_m[3][i] = dump3;
        }
80 81 82 83 84 85 86 87 88 89
    }else 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);
            pick_tmp_m[1][i] = dump1;
            pick_tmp_m[2][i] = dump2;
            pick_tmp_m[3][i] = dump3;
            pick_tmp_m[4][i] = dump4;
            pick_tmp_m[5][i] = dump5;
            pick_tmp_m[6][i] = dump6;
        }
90 91
    }else{
        for(i=1;i<=ntr_glob;i++){
92 93
            fscanf(fptime,"%f",&dump1);
            pick_tmp[i] = dump1;
94 95 96 97 98 99 100
        }
    }
    
    fclose(fptime);
    
    /* distribute picks on CPUs */
    h=1;
101
    if(TW_IND==1){
102 103 104 105 106
        for(i=1;i<=ntr;i++){
            picked_times_m[1][h] = pick_tmp_m[1][recpos_loc[3][i]];
            picked_times_m[2][h] = pick_tmp_m[2][recpos_loc[3][i]];
            picked_times_m[3][h] = pick_tmp_m[3][recpos_loc[3][i]];
            
107 108 109 110 111 112 113 114 115 116 117
            h++;
        }
    }else if(TW_IND==2){
        for(i=1;i<=ntr;i++){
            picked_times_m[1][h] = pick_tmp_m[1][recpos_loc[3][i]];
            picked_times_m[2][h] = pick_tmp_m[2][recpos_loc[3][i]];
            picked_times_m[3][h] = pick_tmp_m[3][recpos_loc[3][i]];
            picked_times_m[4][h] = pick_tmp_m[4][recpos_loc[3][i]];
            picked_times_m[5][h] = pick_tmp_m[5][recpos_loc[3][i]];
            picked_times_m[6][h] = pick_tmp_m[6][recpos_loc[3][i]];
            
118 119 120 121 122 123 124 125 126 127
            h++;
        }
    }else{
        for(i=1;i<=ntr;i++){
            picked_times[h] = pick_tmp[recpos_loc[3][i]];
            
            h++;
        }
    }
    
128
    if(TW_IND==1)
129
        free_matrix(pick_tmp_m,1,3,1,ntr_glob);
130 131
    else if(TW_IND==2)
        free_matrix(pick_tmp_m,1,6,1,ntr_glob);
132 133 134
    else
        free_vector(pick_tmp,1,ntr_glob);
    
135
    if(TW_IND==1){
136 137 138 139 140
        for(i=1;i<=ntr;i++){
        for(j=2;j<=ns;j++){
        
            time = (float)(j * DT);
            
141
            dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
142
            taper = exp(-GAMMA*dumpa);
143
            
144
            dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); 
145
            taper1 = exp(-GAMMA*dumpb);
146 147 148 149 150 151 152 153
            
            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];
154 155 156 157 158 159 160 161 162
        }
        }
    }else if(TW_IND==2){
        for(i=1;i<=ntr;i++){
        for(j=2;j<=ns;j++){
        
            time = (float)(j * DT);
            
            dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
163
            taper = exp(-GAMMA*dumpa);
164 165
            
            dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]); 
166
            taper1 = exp(-GAMMA*dumpb);
167 168 169 170 171 172 173 174 175 176
            
            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]);
177
            taper2 = exp(-GAMMA*dumpc);
178 179
            
            dumpd = (time-picked_times_m[4][i]+picked_times_m[6][i]); 
180
            taper3 = exp(-GAMMA*dumpd);
181 182 183 184 185 186 187 188
            
            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];
189 190 191 192 193 194 195 196 197
            
        }
        }
    }else{
        for(i=1;i<=ntr;i++){
        for(j=2;j<=ns;j++){
        
            time = (float)(j * DT);
            
198
            dumpa = (time-picked_times[i]-TWLENGTH_PLUS);
199
            taper = exp(-GAMMA*dumpa);
200
            
201
            dumpb = (time-picked_times[i]+TWLENGTH_MINUS); 
202
            taper1 = exp(-GAMMA*dumpb);
203 204 205 206 207 208 209 210 211 212 213 214
            
            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];
        }     
        }
    }
    
215
    if(TW_IND==1){
216
        free_matrix(picked_times_m,1,3,1,ntr_glob);
217 218 219 220
    }else if(TW_IND==2){
        free_matrix(picked_times_m,1,6,1,ntr_glob);
        free_matrix(dummysection,1,ntr,1,ns);
    }else{
221
        free_vector(picked_times,1,ntr_glob);
222
    }
223
    
224
} /* end of function time_window.c */