rayt3d.c

seismic software,very useful

	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(!nev){
 	    ixs0 = 0;
	    iys0 = 0;
	    isize = 0;

	    nsize = nzot*nyot*nxot;
	    if( !getparstring("tfile",&tfile) ) {
		tfp = stdout;
	    } 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);
			isize=isize/sizeof(float)/nsize;
			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; 
		}
		fseek2g(tfp,isize*nsize*sizeof(float),SEEK_SET);
	    }
	} 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-ix*nys;
			/* 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 */

	if(fx>fxs || ex<exs || fy>fys || ey<eys || fz>0) 
	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;
/* 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);
	/*
 	v = ealloc1float(nx*ny*nz);
 	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);
	*/
 	v = emalloc(nx*ny*nz*sizeof(float));
 	vxx = emalloc(nx*ny*nz*sizeof(float));
 	vxy = emalloc(nx*ny*nz*sizeof(float));
 	vxz = emalloc(nx*ny*nz*sizeof(float));
 	vyy = emalloc(nx*ny*nz*sizeof(float));
 	vyz = emalloc(nx*ny*nz*sizeof(float));
	fprintf(stderr,"before ealloc vzz nzyxt=%d ncpu=%d\n",nzyxt,ncpu);
 	vzz = emalloc(nx*ny*nz*sizeof(float));
	/*
	vzz = ealloc1float(nx*ny*nz);
	*/
  	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);
 	ov2 = ealloc1float(nzo*nyo*nxo);

/*
	fprintf(stderr,"before read vfp nx=%d ny=%d nz=%d \n",nx,ny,nz);
*/
 
	/* read velocities */
	fseek(vfp,0,0);
	fread(v,sizeof(float),nx*ny*nz,vfp);
  
	/* compute second derivatives of velocity  */
	dv2_(&nx,&ny,&nz,&dx,&dy,&dz,v,vxx,vxy,vxz,vyy,vyz,vzz);
  
 	/* compute slowness squares  */
	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 time arrays   */
 	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);

    i2 = 2;
	i3 = 3;
	i6 = 6;
	n1 = 128;
	n2 = 1001*n1; 
	/* allocate other working arrays */
	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(jpfp," \n");
   	fprintf(jpfp," RAYT3D parameters \n");
   	fprintf(jpfp,"===================\n");
   	fprintf(jpfp," vfile=%s \n",vfile);
   	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);
   	fprintf(jpfp," dz=%g dx=%g dy=%g \n",dz,dy,dx);
   	fprintf(jpfp," nzo=%d nxo=%d nyo=%d \n",nzot,nyot,nxot);
   	fprintf(jpfp," fzo=%g fxo=%g fyo=%g \n",fzot,fyot,fxot);
   	fprintf(jpfp," dzo=%g dxo=%g dyo=%g \n",dzot,dyot,dxot);
   	fprintf(jpfp," nzoe=%d nxoe=%d nyoe=%d \n",nzo,nyo,nxo);
   	fprintf(jpfp," fzoe=%g fxoe=%g fyoe=%g \n",fzo,fyo,fxo);
   	fprintf(jpfp," dzoe=%g dxoe=%g dyoe=%g \n",dzo,dyo,dxo);
	if(nev){
   	fprintf(jpfp," \n");
   	fprintf(jpfp," ntline=%d v0=%g dvz=%g \n",nev,v0,dvz);
	for(iev=0; iev<nev; ++iev) {
   	fprintf(jpfp," \n");
   	fprintf(jpfp," itline=%d tofile=%s \n",iev+1,tofile[iev]);
   	fprintf(jpfp," xbeg=%g xend=%g ybeg=%g yend=%g \n",
		ybeg[iev],yend[iev],xbeg[iev],xend[iev]);
   	fprintf(jpfp," nout=%d dxout=%g dyout=%g \n",
		nout[iev],dyout[iev],dxout[iev]);
	}
	}
   	fprintf(jpfp," \n");