rayt3d.c

seismic software,very useful

	    nout = alloc1int(nev);
	    dxout = alloc1float(nev);
	    dyout = alloc1float(nev);
	    xbeg = alloc1float(nev);
	    xend = alloc1float(nev);
	    ybeg = alloc1float(nev);
	    yend = alloc1float(nev);
	    tofile = (string *) malloc(nev*sizeof(string));
	    tofp = (FILE**) malloc(sizeof(FILE *)*nev);
  	    for(iev=0; iev<nev; ++iev) {
		if(!getnparint(iev+1,"nout",&nout[iev])) 
			err("cannot get nout for event %i. \n",1+iev);
		if(maxno<nout[iev]) maxno = nout[iev];
		if(!getnparfloat(iev+1,"xbeg",&ybeg[iev])) 
			err("cannot get xbeg for event %i. \n",1+iev);
		if(!getnparfloat(iev+1,"xend",&yend[iev])) 
			err("cannot get xend for event %i. \n",1+iev);
		if(!getnparfloat(iev+1,"ybeg",&xbeg[iev])) 
			err("cannot get ybeg for event %i. \n",1+iev);
		if(!getnparfloat(iev+1,"yend",&xend[iev]))
			err("cannot get yend for event %i. \n",1+iev);
		if(xbeg[iev]<fx || xbeg[iev]>ex ||
		   xend[iev]<fx || xend[iev]>ex ||
		   ybeg[iev]<fy || ybeg[iev]>ey ||
 		   yend[iev]<fy || yend[iev]>ey)
		    	err("event %i is out of output range!\n",1+iev);
		dxout[iev] = (xend[iev]-xbeg[iev])/(nout[iev]-1);
		dyout[iev] = (yend[iev]-ybeg[iev])/(nout[iev]-1);
		if(!getnparstring(iev+1,"tofile",&tofile[iev]))
			err("cannot get output file for event %i. \n",1+iev);
 	    }
	    tout = alloc1float(maxno*nzo);
	    xout = alloc1float(maxno);
	    yout = alloc1float(maxno);
	}
/* compute x-y geometry of extraploation/eikonal time table */
	if(nev==0){
		temp = dxot/dx;
		ixot = NINT(temp);
		ixot = MAX(ixot,1);
		dxo = dxot/ixot;
		temp = dyot/dy;
		iyot = NINT(temp);
		iyot = MAX(iyot,1);
		dyo = dyot/iyot;
		nxo = (nxot-1)*ixot+1;
		nyo = (nyot-1)*iyot+1;
		fxo = fxot;
		fyo = fyot;
	} else {
		fminmax(xbeg,nev,&xbmin,&xbmax);		
		fminmax(xend,nev,&xemin,&xemax);		
		fminmax(ybeg,nev,&ybmin,&ybmax);		
		fminmax(yend,nev,&yemin,&yemax);		
		xomin = MIN(xbmin,xemin);
		xomax = MAX(xbmax,xemax);
		yomin = MIN(ybmin,yemin);
		yomax = MAX(ybmax,yemax);
		dxo = dx;
		dyo = dy;
		nxo = (xomax-xomin)/dxo+1.5;
		nyo = (yomax-yomin)/dyo+1.5;
		fxo = xomin;
		fyo = yomin;
/* added for nev>0  z.li 10/2/96 */
		fxot = xomin;
		fyot = yomin;
	} 

	if (nxo<3) {
		nxo = 3;	
		fxo = fxot - dxo;
	}  
	if (nyo<3) {
		nyo = 3; 
		fyo = fyot - dyo;
	}

	exo = fxo+(nxo-1)*dxo;
	eyo = fyo+(nyo-1)*dyo;
	

 	if (!getparint("nys",&nxs)) nxs = 1;
	if (!getparint("nxs",&nys)) nys = 1;
	if (!getparfloat("fys",&fxs)) fxs = fx;
	if (!getparfloat("fxs",&fys)) fys = fy;
	if (!getparfloat("dys",&dxs)) dxs = dxot * 2.;
	if (!getparfloat("dxs",&dys)) dys = dyot * 2.;
	exs = fxs+(nxs-1)*dxs;
	eys = fys+(nys-1)*dys;
	if (!getparfloat("apery",&aperx)) aperx = 0.5*nx*dx;
	if(nxs>1) aperx += dxs;
	if (!getparfloat("aperx",&apery)) apery = 0.5*ny*dy;
	if(nys>1) apery += dys;
	if (!getparfloat("fa",&fa)) fa = 0;
	if (fa<0) err("fa must be not negative!\n");
	if (!getparfloat("da",&da)) da = 1.;
	if (da<0.0001) err("da must be positive!\n");
	if (!getparint("na",&na)) na = 61;
	if (!getparfloat("amin",&amin)) amin = 0;
	if (!getparfloat("amax",&amax)) amax = 90;
	if (amax>180 || amin<0 ) 
		err("polar angle must be within 0 to 180 degrees!\n");	
	if (!getparfloat("azhmin",&azhmin)) azhmin = 0;
	if (!getparfloat("azhmax",&azhmax)) azhmax = 360;
	if (azhmax>360 || azhmin<0 ) 
		err("azimuth angle must be within 0 and 360 degrees!\n");	
 	fa *= PI/180.;
	da *= PI/180.;
	amin *= PI/180.;
	amax *= PI/180.;
	azhmin *= PI/180.;
	azhmax *= PI/180.;
  	if (!getparfloat("fac",&fac)) fac = 0.01;
	fac *= 1.4141;
	if (!getparint("ek",&ek)) ek = 1;
 	if (!getparstring("restart",&restart)) restart = "n"; 
 	if (!getparint("isres",&isres)) isres = -999;
	if (!getparint("ms",&ms)) ms = 1;
	if (!getparfloat("compress",&compress)) compress = 0.;

/* Set up pthreads */
     	pthread_attr_init(&attr);
	pthread_attr_setscope(&attr, PTHREAD_SCOPE_SYSTEM);
 
	if(nev==0) {
if(compress<=1. && compress!=0.) err(" compress=%f must be => 1. \n",compress);
		if(compress>1.) {
			num_points[0] = nxot;
			num_points[1] = nyot;
			num_points[2] = nzot;
			nxyzot = nxot*nyot*nzot;
			out_bytes = (long long)floor((double)(num_points[0]*num_points[1]*num_points[2]*sizeof(float))/(double)compress);
		work_bytes = compress_ttc_bytes(num_points,1);

	/*
		fprintf(stderr," compression used \n");
		fprintf(stderr,"num_points[0] = %d\n",num_points[0]);
		fprintf(stderr,"num_points[1] = %d\n",num_points[1]);
		fprintf(stderr,"num_points[2] = %d\n",num_points[2]);
		fprintf(stderr,"compress = %f\n",compress);
		fprintf(stderr,"null = %f\n",null);
		fprintf(stderr,"out_bytes = %lld\n",out_bytes);
	*/

		}
	} else {
		compress = 0.;
	}

	if(!nev){
 	    ixs0 = 0;
	    iys0 = 0;
	    isize = 0;

	    nsize = nzot*nyot*nxot;
	    if( !getparstring("tfile",&tfile) ) {
		tfp = stdout;
		if(restart[0]=='y') { 
			err(" tfile is needed for restart job \n");
		}
	    } else {
		if((tfp = fopen(tfile,"r"))!=NULL) {
			fclose(tfp);
			tfp = fopen(tfile,"r+");	
		} else {
	   		tfp = fopen(tfile,"w");
		}
		if(restart[0]=='y') { 
			fseek2g(tfp,0,SEEK_END);
			isize = ftell64(tfp);

			if(compress<=1.) {
				isize=isize/sizeof(float)/nsize;
			} else {
				isize=(isize-8)/out_bytes;
			}
			if(isize<0) isize = 0;
			if(isize>isres && isres>0) isize=isres-1;
			isres = isize + 1;
			ixs0 = isize/nys;
			iys0 = isize-ixs0*nys;
			/* if(isize==ns) return 0; */
		} else {
			fclose(tfp);
			tfp = fopen(tfile,"w");	
			isres = 1;
		}
		if(isize==0) {
			if(compress>1.) {
#ifndef WORDS_BIGENDIAN
  				/* Swap bytes before writing. */
    				ttc_byteswap64(out_bytes);
#endif

      				fwrite(&out_bytes, sizeof(long long), 1, tfp);

#ifndef WORDS_BIGENDIAN
  				/* Swap bytes back after writing. */
	  			ttc_byteswap64(out_bytes);
#endif

			}
		}
		if(compress<1.0) {
			isize = isize*nsize*sizeof(float);
		} else {
			isize = isize*out_bytes+8;
		}
		fseek2g(tfp,isize,0);
	    }
	} else {
 	    ixs0 = nxs;
	    iys0 = nys;
	    for(iev=0; iev<nev; ++iev){
			ix = 0;
			iy = 0;
			isize = 0;
			nsize = nzo*nout[iev];

 			if((tofp[iev] = fopen(tofile[iev],"r"))!=NULL) {
				fclose(tofp[iev]);
				tofp[iev] = fopen(tofile[iev],"r+");	
			} else {
	   			tofp[iev] = fopen(tofile[iev],"w");
			}
			if(restart[0]=='y') { 
				fseek2g(tofp[iev],0,SEEK_END);
				isize = ftell64(tofp[iev]);
				isize = isize/sizeof(float)/nsize;
				ix = isize/nys;
				/* iy = isize-ixs0*nys; */
				iy = isize-ix*nys;
				/*
				fprintf(stderr,"isize=%f ix=%d iy=%d \n",(float)isize,ix,iy);
				*/
				/* if(isize==ns) return 0; */
			} else {
				fclose(tofp[iev]);
				tofp[iev] = fopen(tofile[iev],"w");	
				isres = 1;
			}
	    	if(ixs0>ix) {
			ixs0=ix;
	    		iys0=iy;
		}
	    }
	    is0 = ixs0*nys + iys0;
	    if(is0>isres && isres>0) is0 = isres - 1;
	    isres = is0 + 1;
	    ixs0 = is0/nys;
	    iys0 = is0-ixs0*nys;
	    for(iev=0; iev<nev; ++iev){
			nsize = nzo*nout[iev];
			isize = ixs0*nys+iys0;
			fseek2g(tofp[iev],isize*nsize*sizeof(float),0);
	    }
	}

 
	/* ensure sources is in grid */

/* allow some round-off errors */
/*	if(fx>fxs || ex<exs || fy>fys || ey<eys || fz>0)  */
	if(fx>fxs+0.1 || ex<exs-0.1 || fy>fys+0.1 || ey<eys-0.1 || fz>0) {
		if(fx>fxs+0.1) warn("fy=%f fys=%f \n", fx,fxs);
		if(ex<exs-0.1) warn("ey=%f eys=%f \n", ex,exs);
		if(fy>fys+0.1) warn("fx=%f fxs=%f \n", fy,fys);
		if(ey<eys-0.1) warn("ex=%f exs=%f \n", ey,eys);
		if(fz>0) warn("fz=%f \n", fz);

		err("source lies outside of specified v(x,y,z) grid\n");
	}
  
	if(fx>fxo || ex<exo  || fy>fyo || ey<eyo
		|| fz>fzo || ez<fzo+(nzo-1)*dzo) { 
	warn("fxo=%g fx=%g exo=%g ex=%g \n",fyo,fy,eyo,ey);
	warn("fyo=%g fy=%g eyo=%g ey=%g \n",fxo,fx,exo,ex);
	warn("fzo=%g fz=%g ezo=%g ez=%g \n",fzo,fz,fzo+(nzo-1)*dzo,ez);
	err("extrapolation lies outside of specified v(x,y,z) grid \n");
	}
 
	/* compute nxtmax and nytmax for ***t arrays */
	temp = 2.*aperx/dx + 2.5;
	nxtmax = temp;
	if(nxtmax>nx) nxtmax = nx;
	temp = 2.*apery/dy + 2.5;
	nytmax = temp;
	if(nytmax>ny) nytmax = ny;
	nzyxt = nxtmax * nytmax * nz;

	if(!( getparstring("diskvxx",&diskvxx) &&
	      getparstring("diskvxy",&diskvxy) &&
	      getparstring("diskvxz",&diskvxz) &&
	      getparstring("diskvyy",&diskvyy) &&
	      getparstring("diskvyz",&diskvyz) &&
	      getparstring("diskvzz",&diskvzz) )) {
		idisk = 0;
	} else {
		idisk = 1;
	    if(!getparint("idiskc",&idiskc)) idiskc = 0;
	}

/* velocity input */
 	/* allocate space */
	/*
	fprintf(stderr,"nxtmax=%d nytmax=%d nz=%d \n",nxtmax,nytmax,nz);
	fprintf(stderr,"before ealloc v nx=%d ny=%d nz=%d \n",nx,ny,nz);
	fprintf(stderr,"nzo=%d nxo=%d nyo=%d \n",nzo,nxo,nyo);
	*/


	/* the v group of arrays size: 7*nx*ny*nz*sizeof(float) */
 	v = ealloc1float(nx*ny*nz);

	if(idisk==0) { 
 		vxx = ealloc1float(nx*ny*nz);
		vxy = ealloc1float(nx*ny*nz);
		vxz = ealloc1float(nx*ny*nz);
		vyy = ealloc1float(nx*ny*nz);
		vyz = ealloc1float(nx*ny*nz);
		vzz = ealloc1float(nx*ny*nz);
	} else {
 		vv = ealloc1float(nx*ny*nz);
	}


	/* read velocities */
	fseek64(vfp,0,0);
	fread(v,sizeof(float),nx*ny*nz,vfp);
  
	/* compute second derivatives of velocity  */
	if(idisk==0) {
		dv2_(&nx,&ny,&nz,&dx,&dy,&dz,v,vxx,vxy,vxz,vyy,vyz,vzz);
	} else if(idisk==1 && idiskc==0)  {
		diskfp = fopen(diskvxx,"w");
		file2g(diskfp);
		dv2xx_(&nx,&ny,&nz,&dx,&dy,&dz,v,vv);
		fwrite(vv,sizeof(float),nx*ny*nz,diskfp);
		fclose(diskfp);

		diskfp = fopen(diskvxy,"w");
		file2g(diskfp);
		dv2xy_(&nx,&ny,&nz,&dx,&dy,&dz,v,vv);
		fwrite(vv,sizeof(float),nx*ny*nz,diskfp);
		fclose(diskfp);

		diskfp = fopen(diskvxz,"w");
		file2g(diskfp);
		dv2xz_(&nx,&ny,&nz,&dx,&dy,&dz,v,vv);
		fwrite(vv,sizeof(float),nx*ny*nz,diskfp);
		fclose(diskfp);

		diskfp = fopen(diskvyy,"w");
		file2g(diskfp);
		dv2yy_(&nx,&ny,&nz,&dx,&dy,&dz,v,vv);
		fwrite(vv,sizeof(float),nx*ny*nz,diskfp);
		fclose(diskfp);

		diskfp = fopen(diskvyz,"w");
		file2g(diskfp);
		dv2yz_(&nx,&ny,&nz,&dx,&dy,&dz,v,vv);
		fwrite(vv,sizeof(float),nx*ny*nz,diskfp);
		fclose(diskfp);

		diskfp = fopen(diskvzz,"w");
		file2g(diskfp);
		dv2zz_(&nx,&ny,&nz,&dx,&dy,&dz,v,vv);
		fwrite(vv,sizeof(float),nx*ny*nz,diskfp);
		fclose(diskfp);

	}

	if(idisk==1) {
		free(vv);
 		vv = ealloc1float(nz);
	}

	if(idisk==1) {
		lenvxx = strlen(diskvxx); 
		filevxx = (char*)emalloc(lenvxx+1);
		bcopy(diskvxx,filevxx,lenvxx);
		filevxx[lenvxx]='\0';

		lenvxy = strlen(diskvxy); 
		filevxy = (char*)emalloc(lenvxy+1);
		bcopy(diskvxy,filevxy,lenvxy);
		filevxy[lenvxy]='\0';

		lenvxz = strlen(diskvxy); 
		filevxz = (char*)emalloc(lenvxz+1);
		bcopy(diskvxz,filevxz,lenvxz);
		filevxz[lenvxz]='\0';

		lenvyy = strlen(diskvyy); 
		filevyy = (char*)emalloc(lenvyy+1);
		bcopy(diskvyy,filevyy,lenvyy);
		filevyy[lenvyy]='\0';

		lenvyz = strlen(diskvyz); 
		filevyz = (char*)emalloc(lenvyz+1);
		bcopy(diskvyz,filevyz,lenvyz);