rayt3d.c.6-19-2001

seismic software,very useful

   	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," \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");
*/
	    /* 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);
				   
	   /* 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){ 
	    	/* write traveltime  	*/
	        for(ip=0;ip<np;ip++) {
		for(ixt=0;ixt<nxot;ixt++) {
			ixo = ixt*ixot;
			if(nxot<3) ixo +=1;
			for(iyt=0;iyt<nyot;iyt++) {
				iyo = iyt*iyot;
				if(nyot<3) iyo +=1;
				itemp = (ixo*nyo+iyo)*nzo;
  	    		fwrite(t+itemp+ip*nzyxo,sizeof(float),nzo,tfp);
			}
		}
		}
   	      	fflush(tfp);
 	    } else { 
	        for(ip=0;ip<np;ip++) {
		xs = xsp[ip];
		ys = ysp[ip]; 
 		resit_(t+ip*nzyxo,&nyo,&nxo,&nzo,&nr,&dyo,&dxo,&dr,&fyo,&fxo,
			&ys,&xs,tb);
 	    	for(iev=0; iev<nev; ++iev){

		    for(iout=0; iout<nout[iev]; ++iout){
			xout[iout] = xbeg[iev]+iout*dxout[iev];
			yout[iout] = ybeg[iev]+iout*dyout[iev];
		    }
 
		    interp_(&nyo,&nxo,&nzo,&nout[iev],&fyo,&fxo,&dyo,&dxo,
		    	yout,xout,t+ip*nzyxo,tout);

  		    recot_(tout,&nout[iev],&nzo,&nr,&dr,yout,xout,&ys,&xs,tb);
 
   	    	    fwrite(tout,sizeof(float),nout[iev]*nzo,tofp[iev]);
   	      	    fflush(tofp[iev]);

/*
		itmp = nout[iev]*nzo;
		fminmax(tout,itmp, &gmin, &gmax);
		fprintf(stderr," ixs=%d gmin=%g gmax=%g ip=%d \n",
			ixs,gmin,gmax,ip);
*/

	        }
	        }
	    }

	    for(ip=0;ip<np;ip++) {
		isnow = is + ip;
	    if((isnow+1)%ms==0 && (isnow+1)>=ms ) {
fprintf(jpfp," travel time computed at source is=%d xs=%g ys=%g\n",
		isnow+1,ysp[ip],xsp[ip]);
   	      fflush(jpfp);
	    }
	    }
	}

	fprintf(jpfp," RAYT3D job done\n");

	ugh.scale = 1.;
	ugh.dtype = 4;
 	ugh.n1 = nzot;
	ugh.n2 = nyot;
 	ugh.n3 = nxot;
 	ugh.n4 = nys;
	ugh.n5 = nxs;
 	ugh.d1 = dzot;
	ugh.d2 = dyot;
 	ugh.d3 = dxot;
 	ugh.d4 = dys;
	ugh.d5 = dxs;
	ugh.o1 = fzot; 
	ugh.o2 = fyot; 
	ugh.o3 = fxot; 
	ugh.o4 = fys;
	ugh.o5 = fxs;
	ugh.ocdp2 = 0.;
	ugh.oline3 = 0.;
	ugh.dcdp2 = 0.;
	ugh.dline3 = 0.;
	ugh.gmin = 0.;
	ugh.gmax = tmax;
	ugh.orient = 1;
	ugh.gtype = 0;

	isize = nzot*nyot*nxot;
	isize = isize*nys*nxs*sizeof(float);

	for(iev=0; iev<nev; ++iev){
	    ugh.o2 = ybeg[iev];
	    ugh.o3 = xbeg[iev];

	    if(dxout[iev]==0){
		ugh.n3 = 1;
		ugh.n2 = nout[iev];
		ugh.d2 = dyout[iev];
		ugh.d3 = dxo;
	    } else if(dyout[iev]==0){
		ugh.n2 = 1;
		ugh.n3 = nout[iev];
		ugh.d3 = dxout[iev];
		ugh.d2 = dyo;
 	    } else {
		ugh.n2 = nout[iev];
		ugh.n3 = 1;
		ugh.d2 = sqrt(dxout[iev]*dxout[iev]+dyout[iev]*dyout[iev]);
		ugh.d3 = 999999.;
 	    }
	    if(isres>nys*nxs) {
	    	fclose(tofp[iev]);
	    	isize = ugh.n1*ugh.n2*ugh.n3;
	    	isize = isize*ugh.n4*ugh.n5*sizeof(float);
	    	truncate64(tofile[iev],(off64_t)isize);
	    	tofp[iev] = fopen(tofile[iev],"r+");	
	    	fseek2g(tofp[iev],0,SEEK_END);
	    }
 	    fputusghdr(tofp[iev],&ugh); 
 	}
 		 

 	if(!nev) {
	    	if(isres>nys*nxs) {
			fclose(tfp);
	    		if( getparstring("tfile",&tfile) ) {
	    			truncate64(tfile,(off64_t)isize);
	    			tfp = fopen(tfile,"r+");	
    				fseek2g(tfp,0,SEEK_END);
			}
	     	}
		fputusghdr(tfp,&ugh); 
 		fclose(tfp);
	} else {
		for(iev=0; iev<nev; ++iev) fclose(tofp[iev]);
	}
 
 	/* free space */
 	free1float(v);
	if(idisk==0) {
  		free1float(vxx);
 		free1float(vxy);
 		free1float(vxz);
 		free1float(vyy);
 		free1float(vyz);
 		free1float(vzz);
	} else {
 		free1float(vv);
	}
 	free1float(vt);
   	free1float(vxxt);
 	free1float(vxyt);
 	free1float(vxzt);
 	free1float(vyyt);
 	free1float(vyzt);
 	free1float(vzzt);
  	free1float(t);
	free1float(s);
	free1float(xsp);
	free1float(ysp);
	free1float(azhnp);
	free1float(azhxp);
	free1float(fxtp);
	free1float(fytp);
	free1int(nxtp);
	free1int(nytp);
	if(nev){
		free1int(nout);
		free1float(dxout);
		free1float(dyout);
		free1float(xbeg);
		free1float(ybeg);
		free1float(xend);
		free1float(xend);
		free1float(tout);
	}

	free1float(xp);
	free1float(yp);
	free1float(zp);
	free1float(pxp);
	free1float(pyp);
	free1float(pzp);
	free1float(e1xp);
	free1float(e1yp);
	free1float(e1zp);
	free1float(e2xp);
	free1float(e2yp);
	free1float(e2zp);
	free1float(q111p);
	free1float(q112p);
	free1float(q121p);
	free1float(q122p);
	free1float(p211p);
	free1float(p212p);
	free1float(p221p);
	free1float(p222p);
	free1float(q211p);
	free1float(q212p);
	free1float(q221p);
	free1float(q222p);
	free1float(vp);
	free1float(dvdxp);
	free1float(dvdyp);
	free1float(dvdzp);

	free1int(nrsp);
	free1float(a0p);
	free1float(azh0p);
	free1float(p2p);
	free1float(q2p);
	free1float(hp);
	free1float(gradvp);
	free1float(d2tp);

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

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

	return 0;
}