rayt3d.c

seismic software,very useful

			nxyzot = nxot*nyot*nzot;
			out_bytes = (long long)floor((double)(num_points[0]*num_points[1]*num_points[2]*sizeof(float))/(double)compress);
			in_data = (float *)malloc(num_points[0]*num_points[1]*num_points[2]*sizeof(float));
		}
	} 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.) {
				fwrite(&out_bytes, sizeof(long long), 1, tfp);
			}
		}
		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);
		filevyz[lenvyz]='\0';

		lenvzz = strlen(diskvyz); 
		filevzz = (char*)emalloc(lenvzz+1);
		bcopy(diskvzz,filevzz,lenvzz);
		filevzz[lenvzz]='\0';
	}

  
 	/* compute slowness squares  */
 	ov2 = ealloc1float(nzo*nyo*nxo);
	ov2int_(v,&nx,&ny,&nz,&fx,&fy,&fz,&dx,&dy,&dz,ov2,&nxo,&nyo,&nzo,
		&fxo,&fyo,&fzo,&dxo,&dyo,&dzo);

	if(nev){
    	    if (!getparfloat("v0",&v0)) v0 = 5000;
    	    if (!getparfloat("dvz",&dvz)) dvz = 0.;
	    dr = (dxo<dyo)?dxo/3.:dyo/3.;
	    rmax = sqrt((exs-fxo)*(exs-fxo)+(eys-fyo)*(eys-fyo));
	    if(rmax<sqrt((exo-fxs)*(exo-fxs)+(eyo-fys)*(eyo-fys)))
		rmax = sqrt((exo-fxs)*(exo-fxs)+(eyo-fys)*(eyo-fys));
	    nr = 1.5+rmax/dr;
	    tb = alloc1float(nr*nzo);
 	    timeb_(&nr,&nzo,&dr,&dzo,&fzo,&dvz,&v0,tb);
 	}

	/* allocate velocity derivative arrays for ray tracing */
	/* memory for this vt group of arrys: 7*nzyxt*ncpu*sizeof(float) bytes */
  	vt = ealloc1float(nzyxt*ncpu);
 	vxxt = ealloc1float(nzyxt*ncpu);
	vxyt = ealloc1float(nzyxt*ncpu);
	vxzt = ealloc1float(nzyxt*ncpu);
	vyyt = ealloc1float(nzyxt*ncpu);
	vyzt = ealloc1float(nzyxt*ncpu);
	vzzt = ealloc1float(nzyxt*ncpu);

/*
	fprintf(stderr," after alloc  vzzt nzyxt=%d \n",nzyxt);
*/
   
	/* allocate time arrays   */
	/* this t group of arrays size: 2*(nzo+1)*nyo*nxo*ncpu*sizeof(float) */
 	tt1 = ealloc1float(nyo*nxo*ncpu);
 	tt2 = ealloc1float(nyo*nxo*ncpu);
	t = ealloc1float(nzo*nyo*nxo*ncpu);
	s = ealloc1float(nzo*nyo*nxo*ncpu);
	ysp = ealloc1float(ncpu);
	xsp = ealloc1float(ncpu);
	azhnp = ealloc1float(ncpu);
	azhxp = ealloc1float(ncpu);
	fxtp = ealloc1float(ncpu);
	fytp = ealloc1float(ncpu);
	nxtp = ealloc1int(ncpu);
	nytp = ealloc1int(ncpu);

/*
	fprintf(stderr," after alloc  nytp \n");
*/

    i2 = 2;
	i3 = 3;
	i6 = 6;
	n1 = 128;
	if(n1<na) n1 = na;
	n2 = 1001*n1;
	if(nt>1001) n2 = nt * n1;
	/* allocate other working arrays */
	/* this p group of arrays size: 28*128*1001*ncpu*sizeof(float) */
	xp = ealloc1float(n2*ncpu);
	yp = ealloc1float(n2*ncpu);
	zp = ealloc1float(n2*ncpu);
	pxp = ealloc1float(n2*ncpu);
	pyp = ealloc1float(n2*ncpu);
	pzp = ealloc1float(n2*ncpu);
	e1xp = ealloc1float(n2*ncpu);
	e1yp = ealloc1float(n2*ncpu);
	e1zp = ealloc1float(n2*ncpu);
	e2xp = ealloc1float(n2*ncpu);
	e2yp = ealloc1float(n2*ncpu);
	e2zp = ealloc1float(n2*ncpu);
	q111p = ealloc1float(n2*ncpu);
	q112p = ealloc1float(n2*ncpu);
	q121p = ealloc1float(n2*ncpu);
	q122p = ealloc1float(n2*ncpu);
	p211p = ealloc1float(n2*ncpu);
	p212p = ealloc1float(n2*ncpu);
	p221p = ealloc1float(n2*ncpu);
	p222p = ealloc1float(n2*ncpu);
	q211p = ealloc1float(n2*ncpu);
	q212p = ealloc1float(n2*ncpu);
	q221p = ealloc1float(n2*ncpu);
	q222p = ealloc1float(n2*ncpu);
	vp = ealloc1float(n2*ncpu);
	dvdxp = ealloc1float(n2*ncpu);
	dvdyp = ealloc1float(n2*ncpu);
	dvdzp = ealloc1float(n2*ncpu);

	nrsp = ealloc1int(n1*ncpu);
	a0p = ealloc1float(n1*ncpu);
	azh0p = ealloc1float(n1*ncpu);
	p2p = ealloc1float(i2*i2*ncpu);
	q2p = ealloc1float(i2*i2*ncpu);
	hp = ealloc1float(i3*i3*ncpu);
	gradvp = ealloc1float(i3*ncpu);
	d2tp = ealloc1float(i6*ncpu);

	map = ealloc1int(n1*ncpu);
	vs = ealloc1float(n1*ncpu);
	dvdxs = ealloc1float(n1*ncpu);
	dvdys = ealloc1float(n1*ncpu);
	dvdzs = ealloc1float(n1*ncpu);
	uxx = ealloc1float(n1*ncpu);
	uxy = ealloc1float(n1*ncpu);
	uxz = ealloc1float(n1*ncpu);
	uyy = ealloc1float(n1*ncpu);
	uyz = ealloc1float(n1*ncpu);
	uzz = ealloc1float(n1*ncpu);
	tzt = ealloc1float(n1*ncpu);
	xx = ealloc1float(n1*ncpu);
	yy = ealloc1float(n1*ncpu);
	zz = ealloc1float(n1*ncpu);
	pxs = ealloc1float(n1*ncpu);
	pys = ealloc1float(n1*ncpu);
	pzs = ealloc1float(n1*ncpu);
	e1xs = ealloc1float(n1*ncpu);
	e1ys = ealloc1float(n1*ncpu);
	e1zs = ealloc1float(n1*ncpu);
	e2xs = ealloc1float(n1*ncpu);
	e2ys = ealloc1float(n1*ncpu);
	e2zs = ealloc1float(n1*ncpu);
	p111s = ealloc1float(n1*ncpu);
	p112s = ealloc1float(n1*ncpu);
	p121s = ealloc1float(n1*ncpu);
	p122s = ealloc1float(n1*ncpu);
	q111s = ealloc1float(n1*ncpu);
	q112s = ealloc1float(n1*ncpu);
	q121s = ealloc1float(n1*ncpu);
	q122s = ealloc1float(n1*ncpu);
	p211s = ealloc1float(n1*ncpu);
	p212s = ealloc1float(n1*ncpu);
	p221s = ealloc1float(n1*ncpu);
	p222s = ealloc1float(n1*ncpu);
	q211s = ealloc1float(n1*ncpu);
	q212s = ealloc1float(n1*ncpu);
	q221s = ealloc1float(n1*ncpu);
	q222s = ealloc1float(n1*ncpu);
	dxx = ealloc1float(n1*ncpu);
	dyy = ealloc1float(n1*ncpu);
	dzz = ealloc1float(n1*ncpu);
	dpxs = ealloc1float(n1*ncpu);
	dpys = ealloc1float(n1*ncpu);
	dpzs = ealloc1float(n1*ncpu);
	de1x = ealloc1float(n1*ncpu);
	de1y = ealloc1float(n1*ncpu);
	de1z = ealloc1float(n1*ncpu);
	de2x = ealloc1float(n1*ncpu);
	de2y = ealloc1float(n1*ncpu);
	de2z = ealloc1float(n1*ncpu);
	dp111 = ealloc1float(n1*ncpu);
	dp112 = ealloc1float(n1*ncpu);
	dp121 = ealloc1float(n1*ncpu);
	dp122 = ealloc1float(n1*ncpu);
	dq111 = ealloc1float(n1*ncpu);
	dq112 = ealloc1float(n1*ncpu);
	dq121 = ealloc1float(n1*ncpu);
	dq122 = ealloc1float(n1*ncpu);
	dp211 = ealloc1float(n1*ncpu);
	dp212 = ealloc1float(n1*ncpu);
	dp221 = ealloc1float(n1*ncpu);
	dp222 = ealloc1float(n1*ncpu);
	dq211 = ealloc1float(n1*ncpu);
	dq212 = ealloc1float(n1*ncpu);
	dq221 = ealloc1float(n1*ncpu);
	dq222 = ealloc1float(n1*ncpu);

	xt = ealloc1float(n1*ncpu);
	yt = ealloc1float(n1*ncpu);
	zt = ealloc1float(n1*ncpu);
	pxt = ealloc1float(n1*ncpu);
	pyt = ealloc1float(n1*ncpu);
	pzt = ealloc1float(n1*ncpu);
	e1xt = ealloc1float(n1*ncpu);
	e1yt = ealloc1float(n1*ncpu);
	e1zt = ealloc1float(n1*ncpu);
	e2xt = ealloc1float(n1*ncpu);
	e2yt = ealloc1float(n1*ncpu);
	e2zt = ealloc1float(n1*ncpu);
	p111t = ealloc1float(n1*ncpu);
	p112t = ealloc1float(n1*ncpu);
	p121t = ealloc1float(n1*ncpu);
	p122t = ealloc1float(n1*ncpu);
	q111t = ealloc1float(n1*ncpu);
	q112t = ealloc1float(n1*ncpu);
	q121t = ealloc1float(n1*ncpu);
	q122t = ealloc1float(n1*ncpu);
	p211t = ealloc1float(n1*ncpu);
	p212t = ealloc1float(n1*ncpu);
	p221t = ealloc1float(n1*ncpu);
	p222t = ealloc1float(n1*ncpu);
	q211t = ealloc1float(n1*ncpu);
	q212t = ealloc1float(n1*ncpu);
	q221t = ealloc1float(n1*ncpu);
	q222t = ealloc1float(n1*ncpu);
	dxt = ealloc1float(n1*ncpu);
	dyt = ealloc1float(n1*ncpu);
	dzt = ealloc1float(n1*ncpu);
	dpxt = ealloc1float(n1*ncpu);
	dpyt = ealloc1float(n1*ncpu);
	dpzt = ealloc1float(n1*ncpu);
	de1xt = ealloc1float(n1*ncpu);
	de1yt = ealloc1float(n1*ncpu);
	de1zt = ealloc1float(n1*ncpu);
	de2xt = ealloc1float(n1*ncpu);
	de2yt = ealloc1float(n1*ncpu);
	de2zt = ealloc1float(n1*ncpu);
	dp111t = ealloc1float(n1*ncpu);
	dp112t = ealloc1float(n1*ncpu);
	dp121t = ealloc1float(n1*ncpu);
	dp122t = ealloc1float(n1*ncpu);
	dq111t = ealloc1float(n1*ncpu);
	dq112t = ealloc1float(n1*ncpu);
	dq121t = ealloc1float(n1*ncpu);
	dq122t = ealloc1float(n1*ncpu);
	dp211t = ealloc1float(n1*ncpu);
	dp212t = ealloc1float(n1*ncpu);
	dp221t = ealloc1float(n1*ncpu);
	dp222t = ealloc1float(n1*ncpu);
	dq211t = ealloc1float(n1*ncpu);
	dq212t = ealloc1float(n1*ncpu);
	dq221t = ealloc1float(n1*ncpu);
	dq222t = ealloc1float(n1*ncpu);

/*
	fprintf(stderr," after alloc  all arrays \n");
*/
 
   	fprintf(jpfp," \n");
   	fprintf(jpfp," RAYT3D parameters \n");
   	fprintf(jpfp,"===================\n");
   	fprintf(jpfp," vfile=%s \n",vfile);
   	if(nev==0) fprintf(jpfp," tfile=%s \n",tfile);
   	fprintf(jpfp," one-way time: dt=%g nt=%d \n",dt,nt);
   	fprintf(jpfp," nz=%d nx=%d ny=%d \n",nz,ny,nx);
   	fprintf(jpfp," fz=%g fx=%g fy=%g \n",fz,fy,fx);