taper_grad_shot.c 10.3 KB
Newer Older
tilman.metz's avatar
tilman.metz committed
1
/*-----------------------------------------------------------------------------------------
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
tilman.metz's avatar
tilman.metz committed
3
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
4
 * This file is part of IFOS.
5
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
6
 * IFOS is free software: you can redistribute it and/or modify
tilman.metz's avatar
tilman.metz committed
7 8
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
9
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
10
 * IFOS is distributed in the hope that it will be useful,
tilman.metz's avatar
tilman.metz committed
11 12 13
 * 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.
14
 *
tilman.metz's avatar
tilman.metz committed
15
 * 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>.
17
 -----------------------------------------------------------------------------------------*/
tilman.metz's avatar
tilman.metz committed
18

19
/*------------------------------------------------------------------------
tilman.metz's avatar
tilman.metz committed
20 21 22 23
 * taper gradient with a gaussian frame to damp inversion artefacts near the sources and receivers
 sws == 1 vertical taper (for tomography geometry)
 sws == 2 horizontal taper (for reflection geometry)
 sws == 3 local radial taper at the source and receiver positions
24
 sws == 4
25
 ------------------------------------------------------------------------*/
tilman.metz's avatar
tilman.metz committed
26 27 28 29 30

#include "fd.h"

void taper_grad_shot(float ** waveconv,float ** taper_coeff, float **srcpos, int nshots, int **recpos, int ntr, int ishot, int sws)
{
31 32 33 34
    
    /* extern variables */
    
    extern float DH;
tilman.metz's avatar
tilman.metz committed
35
	extern int FREE_SURF, NX, NY, NXG, NYG;
36 37
    extern int NPROCX, NPROCY, MYID, POS[3];
    extern FILE *FP;
38
    extern char TAPER_FILE_NAME[STRING_SIZE];
39
    extern int VERBOSE;
40 41
    extern int USE_WORKFLOW, WORKFLOW_STAGE;
    
42
    /* local variables */
43 44
    int i, j, ii, jj, n;
    int ijc, iy, ix, xx, yy, srctaper_gridpt, i1, j1;
45 46
    
    /*extern int GRADT1, GRADT2, GRADT3, GRADT4;*/
47 48
    float  a, grad_tap, **waveconvtmp;
    char  taper_file[STRING_SIZE];
49 50 51 52 53
    
    extern float SRTRADIUS;
    extern int SRTSHAPE, FILTSIZE;
    float **m, **edgemat, **mm, **msum, minm, maxm, x, y, rad, **taper_coeff_glob;
    float maxrad;
54
    FILE *fp_taper = NULL;
55 56 57 58
    
    
    if(sws==1){
        
tilman.metz's avatar
tilman.metz committed
59 60
        /* =================================== */
        /* taper source and receiver positions */
61 62
        /* =================================== */
        
tilman.metz's avatar
tilman.metz committed
63 64 65
        /* Convert from meters to gridpoints -> minimum 5x5 gridpoints */
        srctaper_gridpt = (int)(ceil(2.0*SRTRADIUS/DH));
        if (srctaper_gridpt<5)  srctaper_gridpt = 5;
66
        
tilman.metz's avatar
tilman.metz committed
67 68 69 70 71
        m               = matrix(1,srctaper_gridpt,1,srctaper_gridpt);
        edgemat         = matrix(1,4,1,1);
        mm              = matrix(1,NYG,1,NXG);
        msum            = matrix(1,NYG,1,NXG);
        taper_coeff_glob= matrix(1,NYG,1,NXG);
72 73
        waveconvtmp     = matrix(0,NY+1,0,NX+1);
        
tilman.metz's avatar
tilman.metz committed
74
        for (iy=1;iy<=NYG;iy++)
75 76
            for (ix=1;ix<=NXG;ix++)  msum[iy][ix] = 1.0;
        
tilman.metz's avatar
tilman.metz committed
77
        MPI_Barrier(MPI_COMM_WORLD);
78
        
tilman.metz's avatar
tilman.metz committed
79 80 81
        /*****************************/
        /* Taper at source positions */
        /*****************************/
82 83
        
        a = 1.0;
tilman.metz's avatar
tilman.metz committed
84 85
        maxrad = sqrt(2.0*SRTRADIUS*SRTRADIUS);
        for (j=1;j<=srctaper_gridpt;j++) {
86 87 88 89 90 91 92 93 94 95 96 97 98
            for (i=1;i<=srctaper_gridpt;i++) {
                x = ((float)i-((float)srctaper_gridpt)/2.0-0.5)*DH;
                y = ((float)j-((float)srctaper_gridpt)/2.0-0.5)*DH;
                rad = sqrt(x*x+y*y);
                
                switch (SRTSHAPE) {
                    case 1:
                        m[j][i] = erf(a*rad/maxrad);
                        break;
                    case 2:
                        if (rad>0)      m[j][i] = log(rad);
                        else            m[j][i] = 0.0;
                        break;
tilman.metz's avatar
tilman.metz committed
99
                }
100
            }
tilman.metz's avatar
tilman.metz committed
101
        }
102
        
tilman.metz's avatar
tilman.metz committed
103 104 105
        /* generate local taper matrix */
        minm = minimum_m(m,srctaper_gridpt,srctaper_gridpt);
        for (j=1;j<=srctaper_gridpt;j++)
106 107
            for (i=1;i<=srctaper_gridpt;i++)  m[j][i] -= minm;
        
tilman.metz's avatar
tilman.metz committed
108 109 110 111 112 113 114 115
        /* normalize taper matrix to max of values at the centre of all 4 taper area edges,     */
        /* not the global maximum, which is located at the corners                              */
        edgemat[1][1] = m[1][srctaper_gridpt/2];
        edgemat[2][1] = m[srctaper_gridpt/2][1];
        edgemat[3][1] = m[srctaper_gridpt/2][srctaper_gridpt];
        edgemat[4][1] = m[srctaper_gridpt][srctaper_gridpt/2];
        maxm = maximum_m(edgemat,1,4);
        for (j=1;j<=srctaper_gridpt;j++)
116 117 118 119
            for (i=1;i<=srctaper_gridpt;i++) {
                m[j][i] /= maxm;
                if (m[j][i]>1.0)  m[j][i] = 1.0;
            }
tilman.metz's avatar
tilman.metz committed
120 121
        /* get central position within the taper */
        ijc = (int)(ceil((float)srctaper_gridpt/2));
122
        
tilman.metz's avatar
tilman.metz committed
123 124 125
        /*********************/
        /* loop over sources */
        /*for (n=1;n<=nshots;n++) {*/
126
        n=ishot;
tilman.metz's avatar
tilman.metz committed
127
        for (iy=1;iy<=NYG;iy++)
128 129
            for (ix=1;ix<=NXG;ix++)  mm[iy][ix] = 1.0;
        
tilman.metz's avatar
tilman.metz committed
130 131 132
        i = iround(srcpos[1][n]/DH);
        j = iround(srcpos[2][n]/DH);
        for (iy=1;iy<=srctaper_gridpt;iy++) {
133 134 135 136 137 138
            for (ix=1;ix<=srctaper_gridpt;ix++) {
                xx = i + ix - ijc;
                yy = j + iy - ijc;
                if ((xx<1) || (xx>NXG) || (yy<1) || (yy>NYG))  continue;
                mm[yy][xx] = m[iy][ix];
            }
tilman.metz's avatar
tilman.metz committed
139
        }
140
        
tilman.metz's avatar
tilman.metz committed
141
        for (iy=1;iy<=NYG;iy++)
142 143 144
            for (ix=1;ix<=NXG;ix++)
                if (msum[iy][ix] > mm[iy][ix])
                    msum[iy][ix] = mm[iy][ix];
145
                
146
        
tilman.metz's avatar
tilman.metz committed
147 148
        minm = minimum_m(msum,NXG,NYG);
        for (iy=1;iy<=NYG;iy++)
149 150
            for (ix=1;ix<=NXG;ix++)  msum[iy][ix] -= minm;
        
tilman.metz's avatar
tilman.metz committed
151 152
        maxm = maximum_m(msum,NXG,NYG);
        for (iy=1;iy<=NYG;iy++)
153 154 155 156 157 158 159 160
            for (ix=1;ix<=NXG;ix++) {
                taper_coeff_glob[iy][ix] = msum[iy][ix]/maxm;
                if ((POS[1]==((ix-1)/NX)) && (POS[2]==((iy-1)/NY))){
                    ii = ix-POS[1]*NX;
                    jj = iy-POS[2]*NY;
                    /*Diese Zeile wurde von Daniel auskommentiert. taper_coeff[jj][ii] = taper_coeff_glob[iy][ix] * ((float)(iy*DH));*/
                    /*taper_coeff[jj][ii] = ((float)(iy*DH)) * ((float)(iy*DH)) * ((float)(iy*DH));*/
                    taper_coeff[jj][ii]=msum[iy][ix]/maxm;
tilman.metz's avatar
tilman.metz committed
161
                }
162 163 164 165
            }
        
        
        /* apply taper on local gradient */
tilman.metz's avatar
tilman.metz committed
166
        for (j=1;j<=NY;j++){
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192
            for (i=1;i<=NX;i++){
                waveconv[j][i]*=taper_coeff[j][i];
                waveconvtmp[j][i] = waveconv[j][i];
            }
        }
        
        
        /* apply filter at shot and receiver points */
        n=ishot;
        i1 = iround(srcpos[1][n]/DH);
        j1 = iround(srcpos[2][n]/DH);
        
        for (i=i1-FILTSIZE;i<=i1+FILTSIZE;i++){
            for (j=j1-FILTSIZE;j<=j1+FILTSIZE;j++){
                if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
                    ii = i-POS[1]*NX;
                    jj = j-POS[2]*NY;
                    if (jj>0){
                        waveconvtmp[jj][ii] = 0.0;
                        taper_coeff[jj][ii] = 0.0;
                    }
                }
            }
        }
        
        
tilman.metz's avatar
tilman.metz committed
193 194
        /* apply taper on local gradient */
        for (j=1;j<=NY;j++){
195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
            for (i=1;i<=NX;i++){
                waveconv[j][i] = waveconvtmp[j][i];
            }
        }
        
        
        free_matrix(m,1,srctaper_gridpt,1,srctaper_gridpt);
        free_matrix(edgemat,1,4,1,1);
        free_matrix(mm,1,NYG,1,NXG);
        free_matrix(msum,1,NYG,1,NXG);
        free_matrix(taper_coeff_glob,1,NYG,1,NXG);
        free_matrix(waveconvtmp,0,NX+1,0,NY+1);
        
        
    }	/* end of sws==1 */
    
211 212 213 214
    /* ======================== */
    /* Read Taper from file     */  
    /* ======================== */
           
215
    if((sws>=2)&&(sws<=4)){
216 217
        
        if (MYID==0&&VERBOSE)
218 219 220
        {
            fprintf(FP,"\n **Message from taper_grad_shot (printed by PE %d):\n",MYID);
            fprintf(FP," Coefficients for gradient tapers are now read in.\n");
tilman.metz's avatar
tilman.metz committed
221 222 223
        }
        
        if(sws==2){
224
        if(USE_WORKFLOW){
225
            sprintf(taper_file,"%s.shot%d_%i.vp",TAPER_FILE_NAME,ishot,WORKFLOW_STAGE);
226 227
            fp_taper=fopen(taper_file,"r");
            if(fp_taper==NULL){
228
                sprintf(taper_file,"%s.shot%d.vp",TAPER_FILE_NAME,ishot);
229 230 231
                fp_taper=fopen(taper_file,"r");
            }
        }else{
232
            sprintf(taper_file,"%s.shot%d.vp",TAPER_FILE_NAME,ishot);
233 234 235 236 237
            fp_taper=fopen(taper_file,"r");
        }
        }
        if(sws==3){   
        if(USE_WORKFLOW){
238
            sprintf(taper_file,"%s.shot%d_%i.vs",TAPER_FILE_NAME,ishot,WORKFLOW_STAGE);
239 240
            fp_taper=fopen(taper_file,"r");
            if(fp_taper==NULL){
241
                sprintf(taper_file,"%s.shot%d.vs",TAPER_FILE_NAME,ishot);
242 243 244
                fp_taper=fopen(taper_file,"r");
            }
        }else{
245
            sprintf(taper_file,"%s.shot%d.vs",TAPER_FILE_NAME,ishot);
246 247 248
            fp_taper=fopen(taper_file,"r");
        }
        }
tilman.metz's avatar
tilman.metz committed
249
        if(sws==4){
250
        if(USE_WORKFLOW){
251
            sprintf(taper_file,"%s.shot%d_%i.rho",TAPER_FILE_NAME,ishot,WORKFLOW_STAGE);
252 253
            fp_taper=fopen(taper_file,"r");
            if(fp_taper==NULL){
254
                sprintf(taper_file,"%s.shot%d.rho",TAPER_FILE_NAME,ishot);
255 256 257
                fp_taper=fopen(taper_file,"r");
            }
        }else{
258
            sprintf(taper_file,"%s.shot%d.rho",TAPER_FILE_NAME,ishot);
259 260 261
            fp_taper=fopen(taper_file,"r");
        }
        }
262
        
263
        if(fp_taper==NULL) {
264
            declare_error("Taper file could not be opened");
265
        }
266
        
tilman.metz's avatar
tilman.metz committed
267 268 269
        /* loop over global grid */
        for (i=1;i<=NXG;i++){
            for (j=1;j<=NYG;j++){
270
                
tilman.metz's avatar
tilman.metz committed
271
                fread(&grad_tap, sizeof(float), 1, fp_taper);
272
                
tilman.metz's avatar
tilman.metz committed
273
                if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
274 275 276 277 278
                    ii=i-POS[1]*NX;
                    jj=j-POS[2]*NY;
                    
                    taper_coeff[jj][ii]=grad_tap;
                    
tilman.metz's avatar
tilman.metz committed
279 280 281
                }
            }
        }
282 283 284 285 286 287 288
        
        for (j=1;j<=NY;j++){
            for (i=1;i<=NX;i++){
                waveconv[j][i]*=taper_coeff[j][i];
            }
        }
        
tilman.metz's avatar
tilman.metz committed
289 290
        fclose(fp_taper);
        
291 292
    }
}