suea2df.c

su 的源代码库

	/* get source type, angle */	
	if (!getparstring("stype",&stype)) stype="p";
	if (!getparfloat("sang",&sang)) sang = 0;

	/* get wavelet type, source time length, average frequency */
	if (!getparstring("wtype",&wtype)) wtype="ri";
	if (!getparfloat("ts",&ts)) ts = 0.05;
	if (!getparfloat("favg",&favg)) favg = 50;

	if (!getparint("qsw",&qsw)) qsw = 0;
	if (!getparint("aniso",&aniso)) aniso = 0;
	if (!getparint("tsw",&tsw)) tsw = 0;
	if (!getparint("verbose",&verbose)) verbose = 0;
	if (!getparint("mode",&mode)) mode = 0;

	/***...begin......... shapshot and source files....... */
	/* get name of output snapshot, source files */ 
	getparstring("snfile",&snfile);
	getparstring("sofile",&sofile);

	/* determine snapshots times and sample numbers */
	nsnap = countparval("snaptime");
	if (nsnap!=0){
		snaptime = alloc1float(nsnap);
		isnap = alloc1int(nsnap);
		getparfloat("snaptime",snaptime);
		warn("nsnap=%d \n",nsnap);
		for (i=0; i<nsnap; i++) isnap[i]=NINT(snaptime[i]/dt);
	}
	/* if requested, open file for snapshots */
	if (*snfile!='\0') {
		if((sneisfp=fopen(snfile,"w"))==NULL)
			err("cannot open snfile=%s\n",snfile);
	}

	/* if requested, open file for source */
	if (*sofile!='\0') {
		if((soeisfp=fopen(sofile,"w"))==NULL)
			err("cannot open sofile=%s\n",sofile);
	}
	/***...end......... shapshot and source files....... */


	/**** begin ....... stiffnesses  .... */
	/* Must have stiffnesses c11 and c55 */
	if (!getparstring("c11file",&c11file)) 	c11file="c11_file";
	if (!getparstring("c55file",&c55file)) 	c55file="c55_file";

	/*.. open file for C11 */
	if((c11filefp=fopen(c11file,"r"))==NULL)
			err("cannot open c11file=%s\n",c11file);
	/*.. open file for C55 */
	if((c55filefp=fopen(c55file,"r"))==NULL)
			err("cannot open c55file=%s\n",c55file);


	if (aniso==1) { 
		/* Must have these additional stiffnesses for anisotropy */
		if (!getparstring("c13file",&c13file))	c13file="c13_file";
		if (!getparstring("c35file",&c35file))	c35file="c35_file";
		if (!getparstring("c33file",&c33file))	c33file="c33_file";
		if (!getparstring("c15file",&c15file))	c15file="c15_file";

		/*.. open file for C13 */
		if((c13filefp=fopen(c13file,"r"))==NULL)
				err("cannot open c13file=%s\n",c13file);

		/*.. open file for C35 */
		if((c35filefp=fopen(c35file,"r"))==NULL)
				err("cannot open c35file=%s\n",c35file);

		/*.. open file for C33 */
		if((c33filefp=fopen(c33file,"r"))==NULL)
				err("cannot open c33file=%s\n",c33file);

		/*.. open file for C15 */
		if((c15filefp=fopen(c15file,"r"))==NULL)
				err("cannot open c15file=%s\n",c15file);
	}
	/**** end ....... stiffnesses .... */
	
	/**** begin ..... boundary conditions .... */
	/* get boundary conditions */
	nbc = countparval("bc");
	if (nbc==4) {
		getparint("bc", bc);	
	} else {
		bc[0] = 10;
		bc[1] = 10;
		bc[2] = 10;
		bc[3] = 10;

		if (!((nbc==4) || (nbc==0)) ) 
			warn("Number of bc %d, using bc=10,10,10,10",bc);
	}
	if (!getparfloat("bc_a",&bc_a)) bc_a = 0.95;
	if (!getparfloat("bc_r",&bc_r)) bc_r = 0;
	warn("%f bc_a  %f bc_r\n",bc_a,bc_r);

	for(i=0; i<4; i++) {
		wbc[i]=2;
		if(bc[i]>2)wbc[i]=bc[i];
	}

	bc0=alloc1float(bc[0]);
	bc0[0]=bc_a;
	bc1=alloc1float(bc[1]);
	bc1[0]=bc_a;
	bc2=alloc1float(bc[2]);
	bc2[0]=bc_a;
	bc3=alloc1float(bc[3]);
	bc3[0]=bc_a;

	for(j=1; j<bc[0]; j++) bc0[j]=bc0[j-1]*bc_a*pow(j,-bc_r);
	for(j=1; j<bc[1]; j++) bc1[j]=bc1[j-1]*bc_a*pow(j,-bc_r);
	for(j=1; j<bc[2]; j++) bc2[j]=bc2[j-1]*bc_a*pow(j,-bc_r);
	for(j=1; j<bc[3]; j++) bc3[j]=bc3[j-1]*bc_a*pow(j,-bc_r);

		/*warn("%d %f\n",j,bc0[j]);*/

 	for(j=0; j<bc[3]; j++) warn(" %d %f\n",j,bc3[j]);
	
	warn("bc=%d %d %d %d",bc[0],bc[1],bc[2],bc[3]);
	warn("wbc=%d %d %d %d",wbc[0],wbc[1],wbc[2],wbc[3]);
	/**** end ..... boundary conditions .... */
	
	/***** begin .... seismogram files .... */
	/* determine seismogram recording lines */
	if (!getparfloat("hsz",&hsz)) hsz = 0.0;
	hs1 = NINT( (hsz - fz)/dz )+wbc[0];		/*horizontal line*/
	if (!getparfloat("vsx",&vsx)) vsx = 0.0;
	vs2 = NINT((vsx - fx)/dx )+wbc[1];		/*vertical line*/
	if (!getparint("verbose",&verbose)) verbose = 0;
	if (!getparstring("hsfile",&hsfile)) hsfile="hs.su";
	if (!getparstring("vrslfile",&vrslfile)) vrslfile="vsp.su";
	warn("hs1=%d vs2=%d \n",hs1, vs2);
	

	/* open files for seismograms */
	if((hseisfp=fopen(hsfile,"w+"))==NULL) 
		err("cannot open hsfile=%s\n",hsfile);
	if((vseisfp=fopen(vrslfile,"w+"))==NULL)
		err("cannot open vrslfile=%s\n",vrslfile);
	/***** end .... seismogram files .... */
		
	/* pad the nx and nz dimensions for boundary condition layers */
	nxpadded=nx+wbc[1]+wbc[3];		
	nzpadded=nz+wbc[0]+wbc[2];

	nt=NINT((lt-ft)/dt)+1;
	warn("dx=%f dz=%f dt=%f",dx,dz,dt);
	warn("nxpadded=%d nzpadded=%d nt=%d",nxpadded,nzpadded,nt);
	warn("array sizes");

	/* set grid parameters*/
	limits = alloc1int(4);
	 	dtx=dt/dx; dtz=dt/dz;
	 	energy=0;efac=1000000000000000000./((nxpadded-4)*(nzpadded-4));
	 	pfac=1./(dx*dz);
	 	is=NINT((-fx+sx)/dx)+wbc[1]; js=NINT((-fz+sz)/dz)+wbc[0];
 	warn("is=%d js=%d\n",is,js);

	/* check to se if source is inside the model */
 	if((sx<fx) || (sx>xmax))
		err("sx=%f source must lie within model \n",sx);
 	if((sz<fz) || (sz>zmax))
		err("sz=%f source must lie within model \n",sz);
		
	/* allocate space for model*/
	c11 = alloc2float(nxpadded,nzpadded);
	c11temp = alloc2float(nx,nz);
	c55 = alloc2float(nxpadded,nzpadded);
	c55temp = alloc2float(nx,nz);
	rho = alloc2float(nxpadded,nzpadded);
	rhotemp = alloc2float(nx,nz);
	ngrids=3;

	/* if attenuation is desired allocate space for q values */
	if(qsw==1) {
		q = alloc2float(nxpadded,nzpadded);
		ngrids+=1;
	}

	/* calculate the number of bytes in grids */
	nbytes=ngrids*nxpadded*nzpadded*FSIZE;

	/* if anisotropy is desired allocate space for anisotropy params */
	if(aniso==1) {
		c33 = alloc2float(nxpadded,nzpadded);
		c33temp = alloc2float(nx,nz);
		c13 = alloc2float(nxpadded,nzpadded);
		c13temp = alloc2float(nx,nz);
		c15 = alloc2float(nxpadded,nzpadded);
		c15temp = alloc2float(nx,nz);
		c35 = alloc2float(nxpadded,nzpadded);
		c35temp = alloc2float(nx,nz);
		ngrids+=4;
	}

	/* calculate the number of bytes in all the grids */
	nbytes=ngrids*nxpadded*nzpadded*FSIZE;
	warn("%d grids allocated for model",ngrids);
	warn("%f MBytes allocated for model",0.000001*nbytes);

	/* set e_*/
	strcpy(tfile,hsfile);
	mfile=strcat(tfile,".mod");

	/* open mfile */	
	if((mfp=fopen(mfile,"w+"))==NULL)
		err("cannot open mfile=%s\n",mfile);

	/***** get density and stiffnesses from input files ***********/
	/* ... densities optional */
	getparstring("rhofile",&rhofile);

 	/*** begin density **/
	/* if specified, read density rho */
	if (*rhofile!='\0') {
		if((rhofilefp=fopen(rhofile,"r"))==NULL)
			err("cannot open rhofile=%s\n",rhofile);
		if (fread(rhotemp[0],sizeof(float),nx*nz,rhofilefp)!=nx*nz)
			err("error reading rhofile=%s\n",rhofile);
		fclose(rhofilefp);

		/* pad the density */
		pad_float_2D(nx,nz,wbc,rhotemp,rho); 
		free2float(rhotemp);



		/* we need to invert density */
		/* find min and max of 1/density */
		rhomin = rhomax = 1/rho[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {

					/* rho = 1/rho */
					rho[iz][ix] = 1/rho[iz][ix];

					rhomin= MIN(rhomin,rho[iz][ix]);
					rhomax = MAX(rhomax,rho[iz][ix]);
			}
		}
			
	 }

	 /* if densities not specified or constant, make densities = 10^{-3} */
	 if (*rhofile=='\0' || rhomin==rhomax ) {
		  for (iz=0; iz<nzpadded; ++iz)
			   for (ix=0; ix<nxpadded; ++ix)
				    rho[iz][ix] = .001;
		  rhomin = rhomax = .001 ;
	 }
 	/** end density **/


	/***** begin c11  */
	if((c11filefp=fopen(c11file,"r"))==NULL)
		err("cannot open c11file=%s\n",c11file);
	if (fread(c11temp[0],sizeof(float),nx*nz,c11filefp)!=nx*nz)
		err("error reading c11file=%s\n",c11file);
	fclose(c11filefp);

	/* pad the c11 */
	pad_float_2D(nx,nz,wbc,c11temp,c11); 
	free2float(c11temp);

	/* find min and max of c11 */
	c11min = c11max = c11[0][0];
	for (iz=0; iz<nzpadded; ++iz) {
		for (ix=0; ix<nxpadded; ++ix) {
			    c11min = MIN(c11min,c11[iz][ix]);
			    c11max = MAX(c11max,c11[iz][ix]);
		}
	}
			
 	/**** end ....... c11 **/

	/***** begin ...........c55  */
	if((c55filefp=fopen(c55file,"r"))==NULL)
		err("cannot open c55file=%s\n",c55file);
	if (fread(c55temp[0],sizeof(float),nx*nz,c55filefp)!=nx*nz)
		err("error reading c55file=%s\n",c55file);
	fclose(c55filefp);

	/* pad the c55 */
	pad_float_2D(nx,nz,wbc,c55temp,c55); 
	free2float(c55temp);

	/* find min and max of c55 */
	c55min = c55max = c55[0][0];
	for (iz=0; iz<nzpadded; ++iz) {
		for (ix=0; ix<nxpadded; ++ix) {
			    c55min = MIN(c55min,c55[iz][ix]);
			    c55max = MAX(c55max,c55[iz][ix]);
		}
	}
			
 	/**** end ....... c55 **/

	if (aniso==1) {  /* if anisotropy */

		/***** begin ...........c13  */
		if((c13filefp=fopen(c13file,"r"))==NULL)
			err("cannot open c13file=%s\n",c13file);
		if (fread(c13temp[0],sizeof(float),nx*nz,c13filefp)!=nx*nz)
			err("error reading c13file=%s\n",c13file);
		fclose(c13filefp);

		/* pad the c13 */
		pad_float_2D(nx,nz,wbc,c13temp,c13); 
		free2float(c13temp);

		/* find min and max of c13 */
		c13min = c13max = c13[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    c13min = MIN(c13min,c13[iz][ix]);
				    c13max = MAX(c13max,c13[iz][ix]);
			}
		}
			
 		/**** end ....... c13 **/

		/***** begin ...........c33  */
		if((c33filefp=fopen(c33file,"r"))==NULL)
			err("cannot open c33file=%s\n",c33file);
		if (fread(c33temp[0],sizeof(float),nx*nz,c33filefp)!=nx*nz)
			err("error reading c33file=%s\n",c33file);
		fclose(c33filefp);

		/* pad the c33 */
		pad_float_2D(nx,nz,wbc,c33temp,c33); 
		free2float(c33temp);

		/* find min and max of c33 */
		c33min = c33max = c33[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    c33min = MIN(c33min,c33[iz][ix]);
				    c33max = MAX(c33max,c33[iz][ix]);
			}
		}
			
 		/**** end ....... c33 **/

		/***** begin ...........c15  */
		if((c15filefp=fopen(c15file,"r"))==NULL)
			err("cannot open c15file=%s\n",c15file);
		if (fread(c15temp[0],sizeof(float),nx*nz,c15filefp)!=nx*nz)
			err("error reading c15file=%s\n",c15file);
		fclose(c15filefp);

		/* pad the c15 */
		pad_float_2D(nx,nz,wbc,c15temp,c15); 
		free2float(c15temp);

		/* find min and max of c15 */
		c15min = c15max = c15[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    c15min = MIN(c15min,c15[iz][ix]);
				    c15max = MAX(c15max,c15[iz][ix]);
			}
		}
 		/**** end ....... c15 **/

		/***** begin ...........c35  */
		if((c35filefp=fopen(c35file,"r"))==NULL)
			err("cannot open c35file=%s\n",c35file);
		if (fread(c35temp[0],sizeof(float),nx*nz,c35filefp)!=nx*nz)
			err("error reading c35file=%s\n",c35file);
		fclose(c35filefp);

		/* pad the c35 */
		pad_float_2D(nx,nz,wbc,c35temp,c35); 
		free2float(c35temp);

		/* find min and max of c35 */
		c35min = c35max = c35[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    c35min = MIN(c35min,c35[iz][ix]);
				    c35max = MAX(c35max,c35[iz][ix]);
			}
		}
			
 		/**** end ....... c35 **/
	}

	warn(" c11max=%f  c11min=%f",c11max,c11min);
	warn(" rhomax=%f  rhomin=%f",rhomax,rhomin);
	warn(" wbc[0]=%d  wbc[1]=%d wbc[2]=%d wbc[3]=%d",
			wbc[0],wbc[1],wbc[2],wbc[3]);

	vlim[0] = sqrt(c11min*rhomax);
	vlim[1] = sqrt(c11max*rhomin);


	/* close mfile */
	fclose(mfp);

	/* set arrays for plane boundary conditions */
	ttl = alloc1float(nzpadded); 
	ttb = alloc1float(nxpadded);
	ttr = alloc1float(nzpadded);
	vl = alloc1float(nzpadded);
	vb = alloc1float(nxpadded);
	vr = alloc1float(nzpadded);

	/* note that rho = 1/density now */
	 for (j=0; j<nzpadded; ++j) {
		vl[j]=sqrt(c11[j][wbc[1]]*rho[j][wbc[1]]);
		vr[j]=sqrt(c11[j][nxpadded-wbc[3]-1]*rho[j][nxpadded-wbc[3]-1]);
	}
	 	for (i=0; i<nxpadded; ++i) {
		vb[i]=sqrt(c11[nzpadded-wbc[2]-1][i]*rho[nzpadded-wbc[2]-1][i]);
	}

	/* ttb: traveltime from source to location on bottom */
	sa=sin(PI*sang/180.);ca=cos(PI*sang/180.);

	if(sang>=0&&sang<90){
		ttb[0]=(-wbc[1])*sa*dx/vb[wbc[1]];
	 		for (i=1; i<nxpadded; ++i) {
			ttb[i]=ttb[i-1]+dx*sa/vb[i];
		}
	}

	if(sang>-90&&sang<0){
		ttb[nxpadded-1]=-(-wbc[3])*sa*dx/vl[nxpadded-wbc[3]-1];
	 		for (i=nxpadded-1; i>0; --i) {
			ttb[i-1]=ttb[i]-dx*sa/vb[i];
		}
	}

	/* ttl, ttr: traveltimes from source to location on LHS, RHS */
	sa=sin(PI*sang/180.);
	ca=cos(PI*sang/180.);
	p=sin(sang*PI/180.)/vl[nzpadded-wbc[2]-1];

	ttl[nzpadded-1]=(-wbc[2])*ca*dz/vl[nzpadded-wbc[2]-1];
	 	for (j=nzpadded-1; j>0; --j) {
		sa=sin(asin(p*vl[j]));
		ca=cos(asin(p*vl[j]));
		ttl[j-1]=ttl[j]+dz*ca/vl[j];
	}

	sa=sin(PI*sang/180.);
	ca=cos(PI*sang/180.);
	p=sin(sang*PI/180.)/vr[nzpadded-wbc[2]-1];
	ttr[nzpadded-1]=(-wbc[2])*ca*dz/vr[nzpadded-wbc[2]-1];

	 	for (j=nzpadded-1; j>0; --j) {
		sa=sin(asin(p*vr[j]));
		ca=cos(asin(p*vr[j]));