rayt3d.c

seismic software,very useful

		if(compress<=1.0) {
	    		/* 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);
				}
			}
			}
		} else {

		/* Set up some thread structures. */
		in_data = (float *)emalloc(nxyzot*np*sizeof(float));
		for(ip=0; ip<np; ip++)
		{
		call_struct[ip].data = in_data 
			+ num_points[0]*num_points[1]*num_points[2]*ip;
		call_struct[ip].num_points = num_points;
		call_struct[ip].ratio = compress;
		call_struct[ip].null_value = null;
		call_struct[ip].thread_number = ip;
		call_struct[ip].buffer = (unsigned char *)emalloc(out_bytes);
		call_struct[ip].work_bytes=work_bytes;
		call_struct[ip].work_space=(unsigned char *)emalloc(work_bytes);
		}
			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;
      				bcopy(t+itemp+ip*nzyxo,
				      in_data+(ixt*nyot+iyt)*nzo+ip*nxyzot,
						nzo*sizeof(float));
				}
			}
			}
			/* compress the travel time */
			for (ip=0;ip<np;ip++) {
	pthread_create(&thread_id[ip], &attr, compress_ttc_pt, (void *)(&call_struct[ip]));
			}
			/* join the threads, write results to outfile */
			for (ip=0;ip<np;ip++) {
			pthread_join(thread_id[ip], NULL);
			/* Write the compressed data to the output stream. */
			fwrite(call_struct[ip].buffer,1,out_bytes,tfp);
			}

		  	/* Flush the output stream. */
		    	fflush(tfp);

       			/* Free the compressed data structure's memory. */
			for (ip=0;ip<np;ip++) {
 				free(call_struct[ip].buffer);
	 			free(call_struct[ip].work_space);
			}
	 		free(in_data);
		}

 	    } 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],isize);
	    	tofp[iev] = fopen(tofile[iev],"r+");	
	    	fseek2g(tofp[iev],0,SEEK_END);
	    }
 	    fputusghdr(tofp[iev],&ugh); 
 	}
 		 

 	if(!nev) {
		if(compress==0) {
	    		if(isres>nys*nxs) {
				fclose(tfp);
	    			if( getparstring("tfile",&tfile) ) {
	    				truncate64(tfile,isize);
	    				tfp = fopen(tfile,"r+");	
    					fseek2g(tfp,0,SEEK_END);
				}
	     		}
			fputusghdr(tfp,&ugh); 
		} else {
			usghed2ghed(&ugh, &gh);
			fwrite(&gh,1,100,tfp);
		}
 		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);
	}
	if(compress<1.0) {
 		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;
}