pup.c 7.74 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 31 32 33 34 35 36
/*-----------------------------------------------------------------------------------------
 * 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 ) {
		extern float DT, DH;
		extern float VEL, DENS, ANGTAPI, ANGTAPO;
		extern char SEIS_FILE[STRING_SIZE];
		extern int SEIS_FORMAT;
		extern FILE *FP;
	
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
37
		int		i, j, k, nx_pot, ny_pot,nxy;
niklas.thiel's avatar
niklas.thiel committed
38 39 40
		float 	nyq_k, nyq_f, dk, df, angi, ango, ka;
		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;
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
45 46 47 48 49 50
					
		/* recompute the dimensions to become values of 2^m */
		nx_pot = (int)(pow(2,(ceil(log(ntr_glob)/log(2)))));
		ny_pot = (int)(pow(2,(ceil(log(ns)/log(2)))));
		nxy = max(nx_pot,ny_pot);
		nx_pot *=2;		/* double columns to avaoid filtering artefacts*/
niklas.thiel's avatar
niklas.thiel committed
51
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
52 53 54 55 56 57 58 59 60 61 62 63 64
		/* test: write filter to file*/
		sprintf(pf1,"%s_data_p_output.bin",SEIS_FILE);
			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);
			
		/* test: write filter to file*/
		sprintf(pf1,"%s_data_vz_output.bin",SEIS_FILE);
			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
niklas.thiel committed
65
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
66 67 68 69 70 71 72 73 74 75 76 77 78 79
		/* 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 */
niklas.thiel's avatar
niklas.thiel committed
80 81
		nyq_k = 1/(2*DH);	/* Nyquist of data in first dimension */
		nyq_f = 1/(2*DT);	/* Nyquist of data in time dimension */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
82 83 84 85 86
		dk = 1/(nx_pot*DH);	/* Wavenumber increment */
		df = 1/(ny_pot*DT);		/* Frequency increment */
		
		for (j=1;j<=ny_pot/2;j++) w[j]=df*(j-1);	/* vector of angular frequencies*/
		for (j=1;j<=nx_pot/2;j++) kx[j]=dk*(j-1);	/*vector of angular wavenumbers*/
niklas.thiel's avatar
niklas.thiel committed
87
		  
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
88 89
		angi = ANGTAPI * PI/180.0;	/* angle at which taper begin */
		ango = ANGTAPO * PI/180.0;	/* angle at which taper ends */
niklas.thiel's avatar
niklas.thiel committed
90 91
		
// 		fprintf(FP,"\n Testposition 2: angi: %f; ango: %f.\n ",angi, ango);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
92 93 94 95 96 97
//  	fprintf(FP,"\n Testposition 2: ns: %i; DT: %f.\n ",ns, DT);
//  	fprintf(FP,"\n Testposition 2: ntr_glob: %i; DH: %f.\n ",ntr_glob, DH);
// 		fprintf(FP,"\n Testposition 2: DENS: %f; VEL: %f.\n ",DENS, VEL);
// 		fprintf(FP,"\n Testposition 3: angi: %f; ango: %f.\n ",angi, ango);

		fprintf(FP,"\n Calculate Filter for wavefield separation...\n ",angi, ango);
niklas.thiel's avatar
niklas.thiel committed
98
		/* calculate Filter vor wavefield separation (after Kluever, 2008) */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
99
		for (k=1;k<=ny_pot/2;k++) {
niklas.thiel's avatar
niklas.thiel committed
100
			ka = fabs(w[k])/VEL;
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
101
			for (i=1;i<=nx_pot/2;i++) {
niklas.thiel's avatar
niklas.thiel committed
102 103 104 105 106 107 108 109 110 111 112
				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
113
			
niklas.thiel's avatar
niklas.thiel committed
114 115
				
		/* expand filter */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
116 117 118 119 120 121 122 123 124 125
		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
126 127 128 129 130
				}
				
			}
		}
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
131 132
		/* test: write filter to file*/
		sprintf(pf1,"%s_spectrum.filt.shot%i.it%i.bin",SEIS_FILE,ishot,iter);
niklas.thiel's avatar
niklas.thiel committed
133
			file = fopen(pf1,"wb");
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
134 135
			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
136
			fclose(file);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
137 138 139 140 141 142 143 144 145 146
					
		/* 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];
			}
		}
		
			fprintf(FP,"2D FFT of p-data and vz-data... \n");
niklas.thiel's avatar
niklas.thiel committed
147
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
148 149 150
		/* transformation into fk domain */
		fft2(data_pt, data_pim, ny_pot, nx_pot,1);

niklas.thiel's avatar
niklas.thiel committed
151
		/* write spectrum to file */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
152
			sprintf(pf2,"%s_spectrum.p.real.shot%i.it%i.bin",SEIS_FILE,ishot, iter);
niklas.thiel's avatar
niklas.thiel committed
153
			file = fopen(pf2,"wb");
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
154 155
			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
156
			fclose(file);
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
157 158 159 160 161 162 163 164 165 166 167 168
		
		fft2(data_vyt, data_vyim,  ny_pot, nx_pot,1);
		
			/* write spectrum to file */
			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);

		
		fprintf(FP,"Perform wavefield separation... \n");
niklas.thiel's avatar
niklas.thiel committed
169 170
		
		/* perform wavefield separation */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
171 172 173 174
		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
175 176
			}
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
177 178 179 180 181 182 183
		/* write spectrum to file */
		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
184
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
185
		fprintf(FP,"\n Reverse 2D FFT...\n ");
niklas.thiel's avatar
niklas.thiel committed
186 187
		
		/* reverse 2d fft */
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
		fft2(data_pup, data_pupim, ny_pot, nx_pot,-1);
		
		fprintf(FP,"\n Write PUP data to file...\n ");
		
		/* write spectrum to file */
		sprintf(pf2,"%s_pup.hand.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);
		
		/* 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
205 206 207 208 209 210 211 212
		
		/* save pup*/
		if (sw == 1){
		  sprintf(pf,"%s_pup.su.shot%i.it%i",SEIS_FILE,ishot,iter);
		}
		if (sw == 2){ /* for step length calculation*/
		  sprintf(pf,"%s_pup.su.step.shot%i.it%i",SEIS_FILE,ishot,iter);
		}
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
213
		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
214
		
niklas.thiel's avatar
BUGFIX  
niklas.thiel committed
215 216 217 218 219 220 221 222 223 224 225 226 227
		fprintf(FP,"\n Write PUP-data into p-data array...\n ");

						
		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
228
}