rayt3d.c

seismic software,very useful

		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 */
	if(compress<=1.) { 
  		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);
   	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");
	fprintf(jpfp," nxs=%d nys=%d fxs=%g fys=%g dxs=%g dys=%g \n",
		nys,nxs,fys,fxs,dys,dxs);			  
	fprintf(jpfp," aperx=%g apery=%g \n",apery,aperx);
	fprintf(jpfp," fa=%g da=%g na=%d \n",fa*180./PI,da*180./PI,na);
	fprintf(jpfp," amin=%g amax=%g azhmin=%g azhmax=%g \n",
		amin*180./PI,amax*180./PI,azhmin*180./PI,azhmax*180./PI);
	fprintf(jpfp," ek=%d ms=%d ncpu=%d jpfile=%s restart=%s isres=%d\n",
		ek,ms,ncpu,jpfile,restart,isres);
	fprintf(jpfp,"compress = %f\n",compress);
	
	if(compress>1) {
		fprintf(jpfp,"num_points[0] = %d\n",num_points[0]);
		fprintf(jpfp,"num_points[1] = %d\n",num_points[1]);
		fprintf(jpfp,"num_points[2] = %d\n",num_points[2]);
		fprintf(jpfp,"null = %f\n",null);
		fprintf(jpfp,"out_bytes = %lld\n",out_bytes);
	}
	

   	fprintf(jpfp," \n");
   	fprintf(jpfp," finish velocity input ... \n");
   	fprintf(jpfp," begin traveltime calculation ... \n");
	fflush(jpfp);

	is0 = ixs0*nys + iys0;
	nzyx = nz * ny * nx;
	nzyxo = nzo * nyo * nxo;

/*
	fprintf(stderr,"is0=%d ixs0=%d nys=%d iys0=%d \n",is0,ixs0,nys,iys0);
*/


	/* loop over sources */
	for (is=is0;is<nxs*nys;is=is+ncpu) {
        np = ncpu;
	    if(is+np>nxs*nys) np = nxs*nys - is;
	    for(ip=0;ip<np;ip++) {
			ixs = (is+ip)/nys;
			iys = (is+ip) - ixs*nys; 
			ysp[ip] = fys + iys*dys;		 
			xsp[ip] = fxs + ixs*dxs;		 
		/*
		fprintf(stderr,"ysp=%f xsp=%f ip=%d fxs=%f ixs=%d dxs=%f iys=%d dys=%f\n",
			ysp[ip],xsp[ip],ip,fxs,ixs,dxs,iys,dys);
		*/
	    }
		
/*
	fprintf(stderr,"nzyx=%d nzyxo=%d nz=%d ny=%d nx=%d \n",
		nzyx,nzyxo,nz,ny,nx);
	fprintf(stderr,"nzo=%d nyo=%d nxo=%d ek=%d nt=%d dt=%g tmax=%g\n",
		nzo,nyo,nxo,ek,nt,dt,tmax);
fprintf(stderr,"na=%d fa=%g da=%g amin=%g amax=%g azhmin=%g azhmax=%g \n",
		na,fa,da,amin,amax,azhmin,azhmax);
fprintf(stderr,"fxo=%g fyo=%g fzo=%g dxo=%g dyo=%g dzo=%g fac=%g \n",
		fxo,fyo,fzo,dxo,dyo,dzo,fac);
fprintf(stderr,"fx=%g fy=%g fz=%g dx=%g dy=%g dz=%g aperx=%g apery=%g \n",
		fx,fy,fz,dx,dy,dz,aperx,apery);
		for(ip=0;ip<np;ip++) 
			fprintf(stderr,"ysp=%g xsp=%g \n",xsp[ip],ysp[ip]);
*/


/*
	fprintf(stderr,"before raytracing \n");
*/

/* allocate velocity derivative arrays for ray tracing */
/* memory for this vt group of arrys: 7*nzyxt*ncpu*sizeof(float) bytes */
	if(compress>1.0) {
  		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);
	}

	    /* parallel computation */
    rt3dp_(&np,&nzyx,&nzyxo,&nz,&ny,&nx,&nzo,&nyo,&nxo,&ek,&nt,
       &na,&fa,&da,&amin,&amax,&azhmin,&azhmax,&dt,&tmax,&fxo,&fyo,&fzo,
       &dxo,&dyo,&dzo,&fac,&fx,&fy,&fz,&dx,&dy,&dz,&aperx,&apery,v,vxx,vxy,vxz,
	   vyy,vyz,vzz,ov2,vt,vxxt,vxyt,vxzt,vyyt,vyzt,vzzt,tt1,tt2,t,s,
       	   xsp,ysp,azhnp,azhxp,fxtp,fytp,nxtp,nytp,&nzyxt,
	   xp,yp,zp,pxp,pyp,pzp,e1xp,e1yp,e1zp,e2xp,e2yp,e2zp,
	   q111p,q112p,q121p,q122p,p211p,p212p,p221p,p222p,
	   q211p,q212p,q221p,q222p,vp,dvdxp,dvdyp,dvdzp,nrsp,
	   a0p,azh0p,&n1,&n2,p2p,q2p,hp,gradvp,d2tp,&i2,&i3,&i6,
	   map,vs,dvdxs,dvdys,dvdzs,uxx,uxy,uxz,uyy,uyz,uzz,tzt,
	   xx,yy,zz,pxs,pys,pzs,e1xs,e1ys,e1zs,e2xs,e2ys,e2zs,
	   p111s,p112s,p121s,p122s,q111s,q112s,q121s,q122s,
   	   p211s,p212s,p221s,p222s,q211s,q212s,q221s,q222s,
	   dxx,dyy,dzz,dpxs,dpys,dpzs,de1x,de1y,de1z,de2x,
	   de2y,de2z,dp111,dp112,dp121,dp122,dq111,dq112,dq121,dq122,
	   dp211,dp212,dp221,dp222,dq211,dq212,dq221,dq222,
	   xt,yt,zt,pxt,pyt,pzt,e1xt,e1yt,e1zt,e2xt,e2yt,e2zt,
	   p111t,p112t,p121t,p122t,q111t,q112t,q121t,q122t,p211t,
	   p212t,p221t,p222t,q211t,q212t,q221t,q222t,dxt,dyt,
	   dzt,dpxt,dpyt,dpzt,de1xt,de1yt,de1zt,de2xt,de2yt,de2zt,
	   dp111t,dp112t,dp121t,dp122t,dq111t,dq112t,dq121t,dq122t,
	   dp211t,dp212t,dp221t,dp222t,dq211t,dq212t,dq221t,dq222t,
	   vv,&idisk,filevxx,filevxy,filevxz,filevyy,filevyz,filevzz,
	   &lenvxx,&lenvxy,&lenvxz,&lenvyy,&lenvyz,&lenvzz);
				   
	/* free memory */
	if(compress>1.0) {
  		free1float(vt);
 		free1float(vxxt);
 		free1float(vxyt);
 		free1float(vxzt);
 		free1float(vyyt);
 		free1float(vyzt);
 		free1float(vzzt);
	}
	   /* check the travel time calculation */
	   for(ip=0;ip<np;ip++) {
		isnow = is + ip;
		itmp = nxo * nyo * nzo;
	   	fminmax(t+ip*itmp,itmp, &gmin, &gmax);
	   	/* print warning message if needed */
		if(gmin > 9999.) {
 		fprintf(jpfp,
		" NULL travel time table at source is=%d xs=%g ys=%g\n",
		isnow+1,ysp[ip],xsp[ip]);
   	      	fflush(jpfp);
 		warn(
		" NULL travel time table at source is=%d xs=%g ys=%g\n",
		isnow+1,ysp[ip],xsp[ip]);
		}
	    }


	    /* output travel time table */
	    if(!nev){