Skip to content

IFOS2D
IFOS2D
Commit f6697382 authored by Florian Wittkamp's avatar Florian Wittkamp
Browse files
Options

Taper: Check if EOF is reached 

Changed output of L2 to STDOUT if WAVETYPE != 3
parent 0d135133
Hide whitespace changes
Inline Side-by-side
Showing with 509 additions and 501 deletions
src/denise.c
View file @ f6697382
......@@ -2802,7 +2802,7 @@ int main(int argc, char **argv){
energy_sum_all_shots = 0.0;
MPI_Allreduce(&energy_all_shots,&energy_sum_all_shots,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
if(MYID==0) printf("\n\n PSV: L2=%f",L2sum_all_shots/energy_sum_all_shots);
if(MYID==0&&(WAVETYPE==3)) printf("\n\n PSV: L2=%f",L2sum_all_shots/energy_sum_all_shots);
}
if(WAVETYPE==2||WAVETYPE==3){
L2sum_SH = 0.0;
......@@ -2814,7 +2814,7 @@ int main(int argc, char **argv){
energy_sum_all_shots_SH = 0.0;
MPI_Allreduce(&energy_all_shots_SH,&energy_sum_all_shots_SH,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
if(MYID==0) printf("\n\n SH: L2=%f",L2sum_all_shots_SH/energy_sum_all_shots_SH);
if(MYID==0&&(WAVETYPE==3)) printf("\n SH: L2=%f",L2sum_all_shots_SH/energy_sum_all_shots_SH);
}
sum_killed_traces=0;
MPI_Allreduce(&killed_traces,&sum_killed_traces,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
......@@ -2835,7 +2835,7 @@ int main(int argc, char **argv){
L2t[1]+=L2sum_SH/energy_sum_SH;
L2t[4]+=L2sum_all_shots_SH/energy_sum_all_shots_SH;
}
if(MYID==0) printf("\n\n Sum: L2=%f",L2t[4]);
if(MYID==0&&(WAVETYPE==3)) printf("\n Sum: L2=%f",L2t[4]);
break;
case 7:
......
src/taper_grad.c
View file @ f6697382
......@@ -2,19 +2,19 @@
* Copyright (C) 2016 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>.
-----------------------------------------------------------------------------------------*/
-----------------------------------------------------------------------------------------*/
/*
* taper gradient with a gaussian frame to damp inversion artefacts near the sources and receivers
......@@ -28,574 +28,582 @@
void taper_grad(float ** waveconv,float ** taper_coeff, float **srcpos, int nshots, int **recpos, int ntr, int sws)
{
/* extern variables */
extern float DH;
extern int FREE_SURF, NX, NY, NXG, NYG;
extern int NPROCX, NPROCY, MYID, POS[3];
extern FILE *FP;
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
/* local variables */
int i, j, h, ifw, ii, jj, n, xb, yb, xe, ye, taperlength,taperlength2, VTON, SRTON;
int ijc, iy, ix, iii, jjj, xx, yy, srctaper_gridpt, i1, j1;
extern int GRADT1, GRADT2, GRADT3, GRADT4;
float amp, a, *window, grad_tap, **waveconvtmp;
char modfile[STRING_SIZE];
extern float SRTRADIUS;
extern int SRTSHAPE, FILTSIZE;
float **m, **edgemat, **mm, **msum, minm, maxm, x, y, rad, **taper_coeff_glob;
float maxrad;
char taper_file[STRING_SIZE];
/* extern variables */
extern float DH;
extern int FREE_SURF, NX, NY, NXG, NYG;
extern int NPROCX, NPROCY, MYID, POS[3];
extern FILE *FP;
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_U[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
/* local variables */
int i, j, h, ifw, ii, jj, n, xb, yb, xe, ye, taperlength,taperlength2, VTON, SRTON;
int ijc, iy, ix, iii, jjj, xx, yy, srctaper_gridpt, i1, j1;
extern int GRADT1, GRADT2, GRADT3, GRADT4;
float amp, a, *window, grad_tap, **waveconvtmp;
char modfile[STRING_SIZE];
extern float SRTRADIUS;
extern int SRTSHAPE, FILTSIZE;
float **m, **edgemat, **mm, **msum, minm, maxm, x, y, rad, **taper_coeff_glob;
float maxrad;
char taper_file[STRING_SIZE];
extern int VERBOSE;
FILE *fp_taper;
/*SRTSHAPE=2;
SRTRADIUS=25.0;
filtsize=2;*/
FILE *fp_taper;
/*SRTSHAPE=2;
SRTRADIUS=25.0;
filtsize=2;*/
/* =============== */
/* Vertical taper */
/* =============== */
if(sws==1){
/* define taper geometry */
taperlength=GRADT2-GRADT1;
taperlength2=GRADT4-GRADT3;
ifw = GRADT4-GRADT1;
/*printf("%d \t %d \t %d \t %d \n",GRADT1, GRADT2, GRADT3, GRADT4);
printf("%d \t %d \t %d \n",taperlength, taperlength2,ifw);*/
if (MYID==0&&VERBOSE)
{
fprintf(FP,"\n **Message from taper_grid (printed by PE %d):\n",MYID);
fprintf(FP," Coefficients for gradient taper are now calculated.\n");
}
waveconvtmp = matrix(0,NY+1,0,NX+1);
window=vector(1,ifw);
/* Gaussian window */
a=3; /* damping coefficient */
for (i=1;i<=ifw;i++){
window[i] = 1.0;
if(i<=taperlength){
window[i] = exp(-(1.0/2.0)*(a*(i-taperlength)/(taperlength/2.0))*(a*(i-taperlength)/(taperlength/2.0)));
}
if(i>=(ifw-taperlength2)){
window[i] = exp(-(1.0/2.0)*(a*(i-(ifw-taperlength2))/(taperlength2/2.0))*(a*(i-(ifw-taperlength2))/(taperlength2/2.0)));
}
}
/* loop over global grid */
for (j=1;j<=NYG;j++){
h=1;
for (i=1;i<=NXG;i++){
grad_tap=0.0;
if((i>GRADT1)&&(i<GRADT4)){
grad_tap=window[h];
h++;
}
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
taper_coeff[jj][ii]=grad_tap;
}
}
}
/* apply taper on local gradient */
for (j=1;j<=NY;j++){
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 */
for (n=1;n<=nshots;n++)
{
i = iround(srcpos[1][n]/DH);
j = iround(srcpos[2][n]/DH);
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
waveconvtmp[jj][ii] = 0.0;
}
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
}
for (n=1;n<=ntr;n++)
{
i = recpos[1][n];
j = recpos[2][n];
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
waveconvtmp[jj][ii] = 0.0;
}
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
}
sprintf(modfile,"taper_coeff_vert.bin");
writemod(modfile,taper_coeff,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
free_vector(window,1,ifw);
free_matrix(waveconvtmp,0,NX+1,0,NY+1);
}
if(sws==1){
/* ======================== */
/* Horizontal Taper */
/* ======================== */
if(sws==2){
/* define taper geometry */
taperlength=GRADT2-GRADT1;
taperlength2=GRADT4-GRADT3;
ifw = GRADT2-GRADT1+1;
if (MYID==0&&VERBOSE){
printf("%d \t %d \t %d \t %d \n",GRADT1, GRADT2, GRADT3, GRADT4);
printf("%d \t %d \t %d \n",taperlength, taperlength2,ifw);
ifw = GRADT4-GRADT1;
/*printf("%d \t %d \t %d \t %d \n",GRADT1, GRADT2, GRADT3, GRADT4);
printf("%d \t %d \t %d \n",taperlength, taperlength2,ifw);*/
if (MYID==0&&VERBOSE)
{
fprintf(FP,"\n **Message from taper_grid (printed by PE %d):\n",MYID);
fprintf(FP," Coefficients for gradient taper are now calculated.\n");
}
waveconvtmp = matrix(0,NY+1,0,NX+1);
window=vector(1,ifw);
/* Gaussian window */
/* Gaussian window */
a=3; /* damping coefficient */
for (i=1;i<=ifw;i++){
window[i] = 1.0;
window[i] = 1.0;
if(i<=taperlength){
window[i] = exp(-(1.0/2.0)*(a*(i-taperlength)/(taperlength/2.0))*(a*(i-taperlength)/(taperlength/2.0)));
window[i] = exp(-(1.0/2.0)*(a*(i-taperlength)/(taperlength/2.0))*(a*(i-taperlength)/(taperlength/2.0)));
}
if(i>=(ifw-taperlength2)){
window[i] = exp(-(1.0/2.0)*(a*(i-(ifw-taperlength2))/(taperlength2/2.0))*(a*(i-(ifw-taperlength2))/(taperlength2/2.0)));
}
}
/* loop over global grid */
h=1;
for (j=1;j<=NYG;j++){
for (i=1;i<=NXG;i++){
grad_tap=0.0;
if((j>=GRADT1)&&(j<=GRADT2)){
grad_tap=window[h];
}
if(j>GRADT2){
/*grad_tap=((float)(j*DH))*((float)(j*DH))*((float)(j*DH));*/
/* grad_tap=((float)(j*DH))*((float)(j*DH));*/
/*grad_tap=((float)(j*DH));*/ /*Version of Daniel */
grad_tap=1.0;
}
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
taper_coeff[jj][ii]=grad_tap;
}
h=1;
for (i=1;i<=NXG;i++){
grad_tap=0.0;
if((i>GRADT1)&&(i<GRADT4)){
grad_tap=window[h];
h++;
}
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
taper_coeff[jj][ii]=grad_tap;
}
if(j>=GRADT1){h++;}
}
}
/* apply taper on local gradient */
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i]*=taper_coeff[j][i];
waveconvtmp[j][i] = waveconv[j][i];
}
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 */
/*for (n=1;n<=nshots;n++)
for (n=1;n<=nshots;n++)
{
i = iround(srcpos[1][n]/DH);
j = iround(srcpos[2][n]/DH);
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;
}
i = iround(srcpos[1][n]/DH);
j = iround(srcpos[2][n]/DH);
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
waveconvtmp[jj][ii] = 0.0;
}
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
}
for (n=1;n<=ntr;n++)
{
i = recpos[1][n];
j = recpos[2][n];
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;
}
i = recpos[1][n];
j = recpos[2][n];
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
/*waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;*/
waveconvtmp[jj][ii] = 0.0;
}
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
}
sprintf(modfile,"taper_coeff_vert.bin");
writemod(modfile,taper_coeff,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
free_vector(window,1,ifw);
free_matrix(waveconvtmp,0,NX+1,0,NY+1);
}
/* ======================== */
/* Horizontal Taper */
/* ======================== */
if(sws==2){
/* define taper geometry */
taperlength=GRADT2-GRADT1;
taperlength2=GRADT4-GRADT3;
ifw = GRADT2-GRADT1+1;
if (MYID==0&&VERBOSE){
printf("%d \t %d \t %d \t %d \n",GRADT1, GRADT2, GRADT3, GRADT4);
printf("%d \t %d \t %d \n",taperlength, taperlength2,ifw);
fprintf(FP,"\n **Message from taper_grid (printed by PE %d):\n",MYID);
fprintf(FP," Coefficients for gradient taper are now calculated.\n");
}
waveconvtmp = matrix(0,NY+1,0,NX+1);
window=vector(1,ifw);
/* Gaussian window */
a=3; /* damping coefficient */
for (i=1;i<=ifw;i++){
window[i] = 1.0;
if(i<=taperlength){
window[i] = exp(-(1.0/2.0)*(a*(i-taperlength)/(taperlength/2.0))*(a*(i-taperlength)/(taperlength/2.0)));
}
}
/* loop over global grid */
h=1;
for (j=1;j<=NYG;j++){
for (i=1;i<=NXG;i++){
grad_tap=0.0;
if((j>=GRADT1)&&(j<=GRADT2)){
grad_tap=window[h];
}
}*/
if(j>GRADT2){
/*grad_tap=((float)(j*DH))*((float)(j*DH))*((float)(j*DH));*/
/* grad_tap=((float)(j*DH))*((float)(j*DH));*/
/*grad_tap=((float)(j*DH));*/ /*Version of Daniel */
grad_tap=1.0;
}
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
taper_coeff[jj][ii]=grad_tap;
}
}
if(j>=GRADT1){h++;}
}
/* apply taper on local gradient */
for (j=1;j<=NY;j++){
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 */
/*for (n=1;n<=nshots;n++)
{
i = iround(srcpos[1][n]/DH);
j = iround(srcpos[2][n]/DH);
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;
}
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
}
for (n=1;n<=ntr;n++)
{
i = recpos[1][n];
j = recpos[2][n];
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii = i-POS[1]*NX;
jj = j-POS[2]*NY;
waveconvtmp[jj][ii] = 1*waveconv[jj][ii]
+ 8*(waveconv[jj-1][ii] + waveconv[jj+1][ii] + waveconv[jj][ii-1] + waveconv[jj][ii+1])
+ 4*(waveconv[jj-1][ii-1] + waveconv[jj-1][ii+1] + waveconv[jj+1][ii+1] + waveconv[jj+1][ii-1]);
waveconvtmp[jj][ii] = waveconvtmp[jj][ii]/49;
}
}
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
waveconv[j][i] = waveconvtmp[j][i];
}
}*/
sprintf(modfile,"taper_coeff_hor.bin");
writemod(modfile,taper_coeff,3);
writemod(modfile,taper_coeff,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
if (MYID==0) mergemod(modfile,3);
free_vector(window,1,ifw);
free_matrix(waveconvtmp,0,NX+1,0,NY+1);
} /* end of sws==2 */
/* =================================== */
/* taper source and receiver positions */
/* =================================== */
if(sws==3) {
/* Convert from meters to gridpoints -> minimum 5x5 gridpoints */
srctaper_gridpt = (int)(ceil(2.0*SRTRADIUS/DH));
if (srctaper_gridpt<5) srctaper_gridpt = 5;
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);
waveconvtmp = matrix(0,NY+1,0,NX+1);
for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++) msum[iy][ix] = 1.0;
MPI_Barrier(MPI_COMM_WORLD);
/*****************************/
/* Taper at source positions */
/*****************************/
a = 1.0;
maxrad = sqrt(2.0*SRTRADIUS*SRTRADIUS);
for (j=1;j<=srctaper_gridpt;j++) {
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;
}
}
} /* end of sws==2 */
/* =================================== */
/* taper source and receiver positions */
/* =================================== */
if(sws==3) {
/* Convert from meters to gridpoints -> minimum 5x5 gridpoints */
srctaper_gridpt = (int)(ceil(2.0*SRTRADIUS/DH));
if (srctaper_gridpt<5) srctaper_gridpt = 5;
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);
waveconvtmp = matrix(0,NY+1,0,NX+1);
for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++) msum[iy][ix] = 1.0;
MPI_Barrier(MPI_COMM_WORLD);
/*****************************/
/* Taper at source positions */
/*****************************/
a = 1.0;
maxrad = sqrt(2.0*SRTRADIUS*SRTRADIUS);
for (j=1;j<=srctaper_gridpt;j++) {
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;
}
/* generate local taper matrix */
minm = minimum_m(m,srctaper_gridpt,srctaper_gridpt);
for (j=1;j<=srctaper_gridpt;j++)
for (i=1;i<=srctaper_gridpt;i++) m[j][i] -= minm;
/* 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++)
for (i=1;i<=srctaper_gridpt;i++) {
m[j][i] /= maxm;
if (m[j][i]>1.0) m[j][i] = 1.0;
}
/* get central position within the taper */
ijc = (int)(ceil((float)srctaper_gridpt/2));
/*********************/
/* loop over sources */
for (n=1;n<=nshots;n++) {
for (iy=1;iy<=NYG;iy++)
for (ix=1;ix<=NXG;ix++) mm[iy][ix] = 1.0;
i = iround(srcpos[1][n]/DH);
j = iround(srcpos[2][n]/DH);
for (iy=1;iy<=srctaper_gridpt;iy++) {
for (ix=1;ix<=srctaper_gridpt;ix++) {
xx = i + ix - ijc;