rsta2d.c

seismic software,very useful

	a3 = (double*) malloc(nint*maxspl*sizeof(double));
	sign = (float*)malloc(nint*sizeof(float));
	norder = (int*)malloc(nint*sizeof(int));

 	t = (float*)malloc(nzo*nxo*sizeof(float));
	work = (float*)malloc(nzo*nxo*sizeof(float));

	/* amplitude array and angle array */
	if(iamp==1) { 
 		a = (float*)malloc(nzo*nxo*sizeof(float));
    		p = (float*)malloc(nzo*nxo*sizeof(float));
	}
	  
	flag = 1;
	one = 1;
	dx = dxo / 2.;
	amin = angmin * 3.141592654/180.;
	amax = angmax * 3.141592654/180.;

	xmin = xint[0];
	xmax = xint[npicks[0]-1];

	/* loop over sources */
	for (is=is0;is<ns;is++) {
		xsi = xs[is];
		zsi = 0.;

		if(fabs(xsi-xmin<0.001) &&  xsi<=xmin) xsi = xsi + 0.001;
		if(fabs(xsi-xmax<0.001) &&  xsi>=xmax) xsi = xsi - 0.001;

		/* velocities at this shot */
		for(i=0;i<nint;i++) {
			nn = npicks[i];
			bisear_(&nn,&one,xpicks+i*maxpks,&xsi,&j);
			v[i] = vint[i+(j-1)*nint];
		} 

	   	/* compute times */
		cshoot(xint, zint, vint, maxspl, nint, npicks,
			a0, a1, a2, a3,
			sign, norder, flag, 
			xsi, angmin, dang, nang, v,
			xstart, xend, dx,
			xray, zray, tray, npts);
		flag = 0;
	
		/* ray checking */
		for(i=0;i<nang;i++) {
			k = 0;
			nc = i * nint;
			for(j=0;j<npts[i];j++) {
				if(xray[nc+j]>=xstart && xray[nc+j]<=xend &&
				   zray[nc+j]<=zmax) {
					xray[nc+k] = xray[nc+j];
					zray[nc+k] = zray[nc+j];
					tray[nc+k] = tray[nc+j];
					k = k + 1;
				}
			}
			npts[i] = k;
		}

		/* plot ray paths */
		if(rayplot!=0) {
			j = 0;
			ncurve = 0;
			for(i=0;i<nhs;i++) {
				j = j + npicks[i];
				ncurve = ncurve + 1;
			} 
			for(i=0;i<nang;i++) {
				j = j + npts[i];
				ncurve = ncurve + 1;
			}
			plotdata = (float*) malloc(j*2*sizeof(float));
			nplots = (int*) malloc(ncurve*sizeof(int));
			k = 0;
			nc = 0; 
			for(i=0;i<nhs;i++) {
				for(j=0;j<npicks[i];j++) {
					plotdata[k] = zpicks[j+i*maxpks];
					plotdata[k+1] = xpicks[j+i*maxpks];
					k = k + 2;
				}
				if(npicks[i]>0) {
					nplots[nc] = npicks[i];
					nc = nc + 1;
				}
			}
			for(i=0;i<nang;i++) {
				for(j=0;j<npts[i];j++) {
					plotdata[k] = zray[j+i*nint];
					plotdata[k+1] = xray[j+i*nint];
					k = k + 2;
				}
				if(npts[i]>0) {
					nplots[nc] = npts[i];
					nc = nc + 1;
				}
			}

			bzero(title,80);
			sprintf(title,"ray paths at source=%d xs=%g \n",
				is+1,xs[is]);
			if(rayplot==1) {
				dump2xgraph(plotdata,nplots,nc,title,
					label1,label2,style);
			}else if(rayplot==-1) {
				dump2psgraph(plotdata,nplots,nc,title,
					label1,label2,style);
			}
			
			free(nplots);
			free(plotdata);
			
		}
		/*	
		for(i=0;i<nang;i++) {
			for(j=0;j<npts[i];j++) {
	fprintf(stderr,"xray=%g zray=%g tray=%g point=%d ray=%d \n",
					xray[i*nint+j],
					zray[i*nint+j],
					tray[i*nint+j],
					j, i);
			}
		}
		*/

	   	/* interpolation time to grid */
		vs = v[0];
		ray2grd(xray, zray, tray, npts,
			nang, nint, 
			xsi, 0, treplace, 
			fxo, fzo, dxo, dzo, nxo, nzo, 
			t, maxfac); 
	
		/* extrapolation travel time if needed in shadow zones */
		if(fillsd==1) {
			for(i=0;i<nxo;i++) {
				for(j=0;j<nzo;j++) {
					zz = fzo + j*dzo;
					if(t[i*nzo+j]==treplace &&
						zz>=zminfl && zz<=zmaxfl) {
						disz = nxo; 
						disx = nzo;
						ix = nxo + 1;
						iz = nzo + 1;
						for(jj=0;jj<nzo;jj++) {
							tmp = jj - j;
							if(t[i*nzo+jj]!=
								treplace &&
								fabs(tmp)<disz){
								iz = jj;
								disz=fabs(tmp);
							}
						}
						for(ii=0;ii<nxo;ii++) {
							tmp = ii - i;
							if(t[ii*nzo+j]!=
								treplace &&
						   		fabs(tmp)<disx){
								ix = ii;
								disx=fabs(tmp);
							}
						}
						if(disx<nxo && disz<nzo) {
							tmp = disx*dxo+disz*dzo;
							work[i*nzo+j] =  
						   	(t[i*nzo+iz]*disx*dxo + 
						   	t[ix*nzo+j]*disz*dzo)
							/tmp;
						} else {
							work[i*nzo+j] =
								t[i*nzo+j];
						}
					} else {
						work[i*nzo+j] = t[i*nzo+j];
					}		
				}
			}
			bcopy(work,t,nxo*nzo*sizeof(float));
		} 
		/* output grid */
	        fwrite(t,sizeof(float),nxo*nzo,tfp);
		fminmax(t,nxo*nzo,&gmin,&gmax);
	   	if(itg==0) {
			tgmin = gmin;
			tgmax = gmax;
			itg = 1;
	   	} else {
			if(tgmin>gmin) tgmin = gmin;
			if(tgmax<gmax) tgmax = gmax;
	   	}

	   	/* compute ampitudes and angles */
	   	if ( iamp == 1 ) {

			/* normalize distance */
			xsi1 = xsi / dzo;
			zsi1 = zsi / dzo;
			fxo1 = fxo / dzo;
			fzo1 = fzo / dzo;
			dzo1 = dzo / dzo;
			dxo1 = dxo / dzo;

	        	amp2d_ (t,a,p,&nzo,&nxo,&amin,&amax,work,
				&amptype,&xsi1,&zsi1,&fxo1,&fzo1,&dxo1,&dzo1);

	        	/* output amplitude */
	        	if (afile[0]!='\0') { 
	         		fwrite(a,sizeof(float),nxo*nzo,afp);
				fminmax(a,nxo*nzo,&gmin,&gmax);
	             		if(iag==0) {
					agmin = gmin;
					agmax = gmax;
					iag = 1;
	   			} else {
					if(agmin>gmin) agmin = gmin;
					if(agmax<gmax) agmax = gmax;
				}
	   		}

	        	/* output propagation angle */
	        	if (pfile[0]!='\0') { 
	         		fwrite(p,sizeof(float),nxo*nzo,pfp);
				fminmax(p,nxo*nzo,&gmin,&gmax);
	             		if(ipg==0) {
					pgmin = gmin;
					pgmax = gmax;
					ipg = 1;
	   			} else {
					if(pgmin>gmin) pgmin = gmin;
					if(pgmax<gmax) pgmax = gmax;
				}
	   		}
	   	}

   fprintf(stderr,"travel time computed at source=%d xs=%g \n",is+1,xs[is]);

	   	if( fabs(smax0-xs[is]) < 0.1*dxs ) {
			scale = 1.e-6;
			dtype = 4;
			n1 = nzo;
			n2 = nxo;
			n3 = ns0;
			n4 = 1;
			n5 = 1;
			d1 = dzo;
			d2 = dxo;
			d3 = dxs;
			d4 = 0.;
			d5 = 0.;
			o1 = fzo; 
			o2 = fxo; 
			o3 = os0; 
			o4 = 0.;
			o5 = 0.;
			ocdp2 = 0.;
			oline3 = 0.;
			dcdp2 = 0.;
			dline3 = 0.;
			gmin = 0.;
			gmax = 0.;
		
			fflush(tfp);
			toghdr(&gh,&scale,&dtype,&n1,&n2,&n3,&n4,&n5,
                		&d1,&d2,&d3,&d4,&d5,&o1,&o2,&o3,&o4,&o5,
                		&dcdp2,&dline3,&ocdp2,&oline3,&tgmin,&tgmax);
			fputghdr(tfp,&gh);

			if (afile[0]!='\0') { 
				fflush(afp);
				putgval(&gh,"gmin",agmin);
				putgval(&gh,"gmax",agmax);
				fputghdr(afp,&gh);
			}
			if (pfile[0]!='\0') {
				fflush(pfp);
				putgval(&gh,"gmin",pgmin);
				putgval(&gh,"gmax",pgmax);
				fputghdr(pfp,&gh);
			}

			/* output hisfile */
			fprintf(hfp,
			"time/amplitude table parameter for KIRMIG \n"); 
			fprintf(hfp,
			"xmint=%f zmint=%f smint=%f \n",fxo,fzo,os0); 
			fprintf(hfp,"dxt=%f dzt=%f dst=%f \n",dxo,dzo,dxs); 
			fprintf(hfp,"nxt=%d nzt=%d nst=%d \n",nxo,nzo,ns0);
			if(amptype==3) {
				fprintf(hfp,"amptype=%d \n",0);
			} else {
				fprintf(hfp,"amptype=%d \n",amptype);
			}

	   	}
	}

	fclose(hfp);
	fclose(tfp);
	if (afile[0]!='\0') fclose(afp);
	if (pfile[0]!='\0') fclose(pfp);
	/* free space */
	free(xpicks);
	if(rayplot!=0) free(zpicks);
	free(npicks);
        free(t);
	free(work);
	if(iamp==1) {
		free(a);
		free(p);
	}
	free(a0);
	free(a1);
	free(a2);
	free(a3);
	free(sign);
	free(norder);
	free(xray);
	free(zray);
	free(tray);
	free(npts);
	
	return 0;
}