pup.c 7.75 KB
Newer Older
niklas.thiel's avatar
niklas.thiel committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
/*-----------------------------------------------------------------------------------------
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
 *
 * This file is part of IFOS.
 * 
 * IFOS 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.
 * 
 * IFOS 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 IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/

/*
 *-------------------------------------------------------------
 *
 * perform wavefield separation in FK domain and save it to file
 *  
 *
 *-------------------------------------------------------------
*/

#include "fd.h"

void	pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int  **recpos, int  **recpos_loc,int ntr, float ** srcpos, int ishot, int ns, int iter, int sw ) {
31
		extern float DT,DH;
niklas.thiel's avatar
niklas.thiel committed
32 33 34 35 36
		extern float VEL, DENS, ANGTAPI, ANGTAPO;
		extern char SEIS_FILE[STRING_SIZE];
		extern int SEIS_FORMAT;
		extern FILE *FP;
	
37 38
		int		i, j, k, nx_pot, ny_pot;
		float 	nyq_k, nyq_f, dk, df, angi, ango, ka, dtr;
niklas.thiel's avatar
niklas.thiel committed
39 40
		float	*w, *kx;
		float	**filtwk, **filtwk_full;
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
41
		float 	**data_pim,**data_vyim,**data_pup,**data_pupim,**data_pt,**data_vyt;
niklas.thiel's avatar
niklas.thiel committed
42 43 44
		char	pf[STRING_SIZE], pf1[STRING_SIZE], pf2[STRING_SIZE];
		FILE	*file;
		
45 46 47
		/* calculate receiver sampling */
		dtr=(recpos[1][ntr_glob]-recpos[1][1])/(ntr_glob-1)*DH;
		fprintf(FP,"\n Receiver sampling was calculatet to dtr: %f m \n ",dtr);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
48
					
49 50 51
		/* double the dimensions to avoid artefacts and recompute the dimensions to become values of 2^m */
		nx_pot = (int)(pow(2,(ceil(log(ntr_glob*2)/log(2)))));
		ny_pot = (int)(pow(2,(ceil(log(ns*2)/log(2)))));
niklas.thiel's avatar
niklas.thiel committed
52
		
53
		/* test: write p-data to file*/
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
54
		if (ishot==1 && iter==1) {
niklas.thiel's avatar
niklas.thiel committed
55
		sprintf(pf1,"%s_data_p_.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
56 57 58 59
			file = fopen(pf1,"wb");
			for (i=1;i<=ntr_glob;i++)
				for (j=1;j<=ns;j++)	fwrite(&data_p[i][j],sizeof(float),1,file);
			fclose(file);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
60
		}
61
		/* test: write vz-data to file*/
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
62
		if (ishot==1 && iter==1) {
niklas.thiel's avatar
niklas.thiel committed
63
		sprintf(pf1,"%s_data_vz_.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
64 65 66 67
			file = fopen(pf1,"wb");
			for (i=1;i<=ntr_glob;i++)
				for (j=1;j<=ns;j++)	fwrite(&data_vy[i][j],sizeof(float),1,file);
			fclose(file);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
68
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
69 70 71 72 73 74 75 76 77 78 79 80 81 82
		/* temporary array */
		filtwk=matrix(1,ny_pot/2,1,nx_pot/2);
		filtwk_full=matrix(1,ny_pot,1,nx_pot);
		data_pim=matrix(1,ny_pot,1,nx_pot);
		data_vyim=matrix(1,ny_pot,1,nx_pot);
		data_pup=matrix(1,ny_pot,1,nx_pot);
		data_pupim=matrix(1,ny_pot,1,nx_pot);
		data_pt=matrix(1,ny_pot,1,nx_pot);
		data_vyt=matrix(1,ny_pot,1,nx_pot);
		
		w=vector(1,ny_pot/2);
		kx=vector(1,nx_pot/2);

		/* calculate variable for wave field separation */
83
		nyq_k = 1/(2*dtr);	/* Nyquist of data in first dimension */
niklas.thiel's avatar
niklas.thiel committed
84
		nyq_f = 1/(2*DT);	/* Nyquist of data in time dimension */
85 86
		dk = 1/((nx_pot/2)*dtr);	/* Wavenumber increment */
		df = 1/((ny_pot/2)*DT);		/* Frequency increment */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
87
		
88 89
		for (j=1;j<=ny_pot/2;j++) w[j]=df*(j-1)*2*PI;	/* vector of angular frequencies*/
		for (j=1;j<=nx_pot/2;j++) kx[j]=dk*(j-1)*2*PI;	/* vector of angular wavenumbers*/
niklas.thiel's avatar
niklas.thiel committed
90
		  
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
91 92 93
		angi = ANGTAPI * PI/180.0;	/* angle at which taper begin */
		ango = ANGTAPO * PI/180.0;	/* angle at which taper ends */

94 95
		/* calculate Filter vor wavefield separation for first quadrant (after Kluever, 2008) */
		fprintf(FP,"Calculate Filter for wavefield separation... \n ");
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
96
		for (k=1;k<=ny_pot/2;k++) {
niklas.thiel's avatar
niklas.thiel committed
97
			ka = fabs(w[k])/VEL;
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
98
			for (i=1;i<=nx_pot/2;i++) {
niklas.thiel's avatar
niklas.thiel committed
99 100 101 102 103 104 105 106 107 108 109
				if(kx[i] < ka*sin(angi))	{
					filtwk[k][i] = (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i]));
				}else{
					if(kx[i] < ka*sin(ango) && kx[i] > ka*sin(angi)){
						filtwk[k][i] = 0.5*(1.0+cos(PI/((ka*sin(ango))-(ka*sin(angi))) * (kx[i]-(ka*sin(angi))))) * (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i]));
					}else {
						filtwk[k][i] = 0.0;
					}
				}
			}
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
110
			
111
		/* expand filter to all 4 quadrants */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
112 113 114 115 116 117 118 119 120 121
		for (k=1;k<=ny_pot;k++) { 
			for (i=1;i<=nx_pot;i++) {
				if(k<=ny_pot/2 && i<=nx_pot/2) {
					filtwk_full[k][i]=filtwk[(ny_pot/2)-(k-1)][(nx_pot/2)-(i-1)];
				}else if(k<=ny_pot/2 && i>nx_pot/2) {
					filtwk_full[k][i]=filtwk[(ny_pot/2)-(k-1)][i-(nx_pot/2)];
				}else if(k>ny_pot/2 && i<=nx_pot/2) {
					filtwk_full[k][i]=filtwk[k-(ny_pot/2)][(nx_pot/2)-(i-1)];
				}else if(k>ny_pot/2 && i>nx_pot/2) {
					filtwk_full[k][i]=filtwk[k-(ny_pot/2)][i-(nx_pot/2)];
niklas.thiel's avatar
niklas.thiel committed
122 123 124 125 126
				}
				
			}
		}
		
127
		/* write filter to file*/
niklas.thiel's avatar
niklas.thiel committed
128 129
		if (iter==1 && ishot==1) {
		sprintf(pf1,"%s_spectrum.filter.bin",SEIS_FILE);
niklas.thiel's avatar
niklas.thiel committed
130
			file = fopen(pf1,"wb");
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
131 132
			for (i=1;i<=nx_pot;i++)
				for (j=1;j<=ny_pot;j++)	fwrite(&filtwk_full[j][i],sizeof(float),1,file);
niklas.thiel's avatar
niklas.thiel committed
133
			fclose(file);
niklas.thiel's avatar
niklas.thiel committed
134 135
		}
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
136 137 138 139 140 141 142 143
		/* copy data to temp array */
		for(k=1;k<=ntr_glob;k++) {
			for (i=1;i<=ns;i++) {
			  data_pt[i][k]=data_p[k][i];
			  data_vyt[i][k]=data_vy[k][i];
			}
		}
		
144
		fprintf(FP,"2D FFT of p-data and vz-data... \n");
niklas.thiel's avatar
niklas.thiel committed
145
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
146 147 148
		/* transformation into fk domain */
		fft2(data_pt, data_pim, ny_pot, nx_pot,1);

niklas.thiel's avatar
niklas.thiel committed
149
		/* write spectrum to file */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
150
		if (ishot==1 && iter==1) {
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
151
			sprintf(pf2,"%s_spectrum.p.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
niklas.thiel's avatar
niklas.thiel committed
152
			file = fopen(pf2,"wb");
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
153 154
			for (i=1;i<=nx_pot;i++)
				for (j=1;j<=ny_pot;j++)	fwrite(&data_pt[j][i],sizeof(float),1,file);
niklas.thiel's avatar
niklas.thiel committed
155
			fclose(file);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
156
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
157 158 159 160
		
		fft2(data_vyt, data_vyim,  ny_pot, nx_pot,1);
		
			/* write spectrum to file */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
161
		if (ishot==1 && iter==1) {
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
162 163 164 165 166
			sprintf(pf2,"%s_spectrum.vz.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
			file = fopen(pf2,"wb");
			for (i=1;i<=nx_pot;i++)
				for (j=1;j<=ny_pot;j++)	fwrite(&data_vyt[j][i],sizeof(float),1,file);
			fclose(file);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
167
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
168 169
		
		fprintf(FP,"Perform wavefield separation... \n");
niklas.thiel's avatar
niklas.thiel committed
170 171
		
		/* perform wavefield separation */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
172 173 174 175
		for (k=1;k<=nx_pot;k++) {
			for (i=1;i<=ny_pot;i++) {
				data_pup[i][k]=0.5*(data_pt[i][k]-(filtwk_full[i][k]*data_vyt[i][k]));
				data_pupim[i][k]=0.5*(data_pim[i][k]-(filtwk_full[i][k]*data_vyim[i][k]));
niklas.thiel's avatar
niklas.thiel committed
176 177
			}
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
178
		/* write spectrum to file */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
179 180 181 182 183 184 185
		if (ishot==1 && iter==1) {
			sprintf(pf2,"%s_spectrum.pup.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
			file = fopen(pf2,"wb");
			for (i=1;i<=nx_pot;i++)
				for (j=1;j<=ny_pot;j++)	fwrite(&data_pup[j][i],sizeof(float),1,file);
			fclose(file);
		}
niklas.thiel's avatar
niklas.thiel committed
186
		
187
		fprintf(FP,"Reverse 2D FFT...\n ");
niklas.thiel's avatar
niklas.thiel committed
188 189
		
		/* reverse 2d fft */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
190 191
		fft2(data_pup, data_pupim, ny_pot, nx_pot,-1);
		
192
		fprintf(FP,"Write PUP data to file...\n ");
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
193 194
		
		/* write spectrum to file */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
195 196 197 198 199 200 201
		if (ishot==1 && iter==1) {
			sprintf(pf2,"%s_data_pup.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
			file = fopen(pf2,"wb");
			for (i=1;i<=nx_pot;i++)
				for (j=1;j<=ny_pot;j++)	fwrite(&data_pup[j][i],sizeof(float),1,file);
			fclose(file);
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
202
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
203
		fprintf(FP,"Write PUP-data into p-data array...\n ");
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
204 205 206 207 208 209
		/* overwrite p data */
		for (k=1;k<=ntr_glob;k++) {
			for (i=1;i<=ns;i++) {
				data_p[k][i]=data_pup[i][k];
			}
		}
niklas.thiel's avatar
niklas.thiel committed
210 211 212
		
		/* save pup*/
		if (sw == 1){
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
213
		  sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot);
niklas.thiel's avatar
niklas.thiel committed
214 215
		}
		if (sw == 2){ /* for step length calculation*/
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
216
		  sprintf(pf,"%s_pup.su.step.shot%i",SEIS_FILE,ishot);
niklas.thiel's avatar
niklas.thiel committed
217
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
218
		outseis_glob(fp,fopen(pf,"w"), 0, data_p,recpos,recpos_loc,ntr,srcpos,1,ns,SEIS_FORMAT,ishot,1);
niklas.thiel's avatar
niklas.thiel committed
219
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
220 221 222 223 224 225 226 227 228 229 230 231

						
		free_matrix(filtwk, 1,ny_pot/2,1,nx_pot/2);
		free_matrix(filtwk_full, 1,ny_pot,1,nx_pot);
		free_matrix(data_pim, 1,ny_pot,1,nx_pot);
		free_matrix(data_vyim, 1,ny_pot,1,nx_pot);
		free_matrix(data_pup, 1,ny_pot,1,nx_pot);
		free_matrix(data_pupim, 1,ny_pot,1,nx_pot);
		free_matrix(data_pt, 1,ny_pot,1,nx_pot);
		free_matrix(data_vyt, 1,ny_pot,1,nx_pot);
		free_vector(w,1,ny_pot/2);
		free_vector(kx,1,ny_pot/2);
niklas.thiel's avatar
niklas.thiel committed
232
}