iipick.c

seismic software,very useful

		z2 = (float *) malloc(nhs*sizeof(float));
		/* spline interpolation of velocities --- not a good idea 
		vt2 = (float *) malloc(maxpks*sizeof(float));
		vb2 = (float *) malloc(maxpks*sizeof(float));
		va2 = (float *) malloc(maxpks*sizeof(float));
		dvdz2 = (float *) malloc(maxpks*sizeof(float));
		*/

	

		dz = zmax/(nzgrid-2);
		for(iz=0;iz<nzgrid;iz++) zgrid[iz] = iz * dz;
		/* compute average velocity and gradient at velo location and
		   at depth zgrid */


		for(i2=0;i2<ncdp_velo;i2++) {
			n2 = i2*nzgrid;
			/* find gradient */
			x = x_velo[i2];

			if(ncdp_dvdz==0) {
				for(iz=0;iz<nzgrid;iz++) 
					dvdzgrid[iz+n2] = 0.;
			} else if (ncdp_dvdz==1) {
				n = nps_dvdz[0];
				dvdzint(zt_dvdz,zb_dvdz,dvdz_dvdz,n,
					zgrid,dvdzgrid+n2,nzgrid);
			} else {
				bisear_(&ncdp_dvdz,&one,x_dvdz,&x,&j2);
				if(x<=x_dvdz[0]) {
					n = nps_dvdz[0];
					dvdzint(zt_dvdz,zb_dvdz,dvdz_dvdz,n,
					     zgrid,dvdzgrid+n2,nzgrid);
				} else if (x>=x_dvdz[ncdp_dvdz-1]) {
					n = nps_dvdz[ncdp_dvdz-1];
					i1 = (ncdp_dvdz-1)*maxp;
					dvdzint(zt_dvdz+i1,zb_dvdz+i1,
						dvdz_dvdz+i1,n,zgrid,
						dvdzgrid+n2,nzgrid);
				} else {
					j1 = j2 - 1;
					n = nps_dvdz[j1];
					i1 = j1*maxp;
					dvdzint(zt_dvdz+i1,zb_dvdz+i1,
						dvdz_dvdz+i1,n,
						zgrid,g1,nzgrid);
					n = nps_dvdz[j2];
					i1 = j2*maxp;
					dvdzint(zt_dvdz+i1,zb_dvdz+i1,
						dvdz_dvdz+i1,n,
						zgrid,g2,nzgrid);
					tmp = (x-x_dvdz[j1]) /
					      (x_dvdz[j2]-x_dvdz[j1]);
				
					for(iz=0;iz<nzgrid;iz++) {
						dvdzgrid[iz+n2] = g1[iz]+
						tmp * (g2[iz]-g1[iz]);
					}
				}
			}
			/* find average velocity at zgrid */
			n = nps_velo[i2];
			n1 = i2*maxp;
			j1 = 0;
			j2 = 0;
			vconint(z_velo+n1,v_velo+n1,n,zgrid,vagrid+n2,nzgrid,
				j1,j2,dvdzgrid+n2);

		}
			


		/* compute vt and vb along horizons at velan locations */

		one = 1;
		for(ih=0;ih<nhs;ih++) nvh[ih] = 0;
		for(i2=0;i2<ncdp_velo;i2++) {
			x = x_velo[i2];	
			jh = 0;
			/* find cross points of horizons at velan locations */

			for(ih=0;ih<nhs;ih++) {
				xmin = xpicks[ih*maxpks]-xtol;
				xmax = xpicks[ih*maxpks+npicks[ih]-1]+xtol;
				n1 = npicks[ih]; 
				n2 = ih * maxpks;
				if(x>=xmin && x<=xmax) {
					ivs[jh] = ih;
					/* find z */
					bisear_(&n1,&one,xpicks+n2,&x,&i1);
					if(x<=xpicks[n2]) {
						z = zpicks[n2];
					} else if(x>=xpicks[n1+n2-1]) {
						z = zpicks[n1+n2-1];
					} else {
						/*
						if(n1==2) {
						*/
						/* linear interpolation */
						   z = zpicks[n2+i1-1] +
						    (x-xpicks[n2+i1-1])*
						    (zpicks[n2+i1]-
						     zpicks[n2+i1-1])/
						    (xpicks[n2+i1]-
						     xpicks[n2+i1-1]);
						/*
						} else {
						*/
						/* spline interpolation */
						/* not a good idea */
						/*
						   z21 = 0.;
                                 		   z2n = 0.;
                                                   spline_(xpicks+n2,zpicks+n2,
                                         		   &n1,&z21,&z2n,z2);
                                 		   splint_(xpicks+n2,zpicks+n2,
							   z2,&n1,&x,&z);
						}
						*/
					}
					zvs[jh] = z; 
					jh = jh + 1;
				}
			}

			/* sort zvs into ascending order */
			for (iz=0;iz<jh;iz++) {
				isort[iz] = iz;
				zsort[iz] = zvs[iz];
			}

			qkisort(jh,zsort,isort);
			for(iz=0;iz<jh;iz++) zvs[iz] = zsort[isort[iz]];

			/* find vt and vb */
			for(iz=0;iz<jh;iz++) {

				z = zvs[iz];

				if(iz==0) {
					zt = z-0.1*dz;
				} else {
					zt = zvs[iz-1];
				}	
				if(iz==jh-1) {
					zb = z+dz*0.1;
				} else {
					zb = zvs[iz+1];
				}

				tmp = zt/dz;
				i1 = tmp;
				if (i1<0) {
					i1 = 0;
				} else if (i1>=nzgrid) {
					i1 = nzgrid-1;
				}
				vat = vagrid[i2*nzgrid+i1];
				gt = dvdzgrid[i2*nzgrid+i1];

				tmp = z/dz;
				i1 = tmp;
				if (i1<0) {
					i1 = 0;
				} else if (i1>=nzgrid) {
					i1 = nzgrid-1;
				}
				va = vagrid[i2*nzgrid+i1];
				g = dvdzgrid[i2*nzgrid+i1];
				
				tmp = zb/dz;
				i1 = tmp;
				if (i1<0) {
					i1 = 0;
				} else if (i1>=nzgrid) {
					i1 = nzgrid-1;
				}
				vab = vagrid[i2*nzgrid+i1];
				gb = dvdzgrid[i2*nzgrid+i1];

				/*
				fprintf(stderr,
					"zt=%f zb=%f x=%f z=%f \n",
					zt,zb,x,z);

				fprintf(stderr,
					"va=%f vab=%f vat=%f \n",
					va,vab,vat);
				*/
				
				i1 = ivs[isort[iz]]*ncdp_velo
					+nvh[ivs[isort[iz]]];

				vth[i1] = (z*va-zt*vat)/(z-zt) + 0.5*gt*(z-zt);
				vbh[i1] = (zb*vab-z*va)/(zb-z) - 0.5*gb*(zb-z);
				vah[i1] = va;
				dvdzh[i1] = g;
				xvh[i1] = x;
				zvh[i1] = z;
				nvh[ivs[isort[iz]]] += 1;
				if(iz==0) {
					vth[i1] = va;
				} else if (iz==jh-1) {
					vbh[i1] = vth[i1];
				}
				/*
				fprintf(stderr,
			"vt=%6.0f vb=%6.0f va=%6.0f g=%6.0f x=%6.0f z=%6.0f \n",
					vth[i1],vbh[i1],vah[i1],dvdzh[i1],x,z);
				*/
			}
		}


		/* fing vt, bt, dvdz and va at pick positions */
		for(ih=lhnvelo;ih<nhs;ih++) {
			tty = efopen("/dev/tty","r");
                        fprintf(stderr,
      "Replace Top velocities in horizon %d with those in VELO? (y/n) ->",ih+1);
                        vtopyes = getc(tty);
                        efclose(tty);
			tty = efopen("/dev/tty","r");
                        fprintf(stderr,
  "Replace Bottom velocities in horizon %d with those in VELO? (y/n) ->",ih+1);
                        vbotyes = getc(tty);
                        efclose(tty);

                        if(vtopyes=='n' && vbotyes=='n') continue;

			n2 = ih*maxpks;
			n1 = ih*ncdp_velo;
			nvs = nvh[ih];
			if(nvs==0) {
				fprintf(stderr, 
				"warning: no velan at horizon %d\n",ih+1); 
				fprintf(stderr, 
				"warning: editing of ifile needed !!!\n"); 
				for(ip=0;ip<npicks[ih];ip++) {
					if(vtopyes=='y') veltops[ip+n2] = 0.;
                                	if(vbotyes=='y') velbots[ip+n2] = 0.;
                                	velavgs[ip+n2] = 0.;
                                	dvdzs[ip+n2] = 0.;
				}
			} else if (nvs==1) {
				for(ip=0;ip<npicks[ih];ip++) {
					if(vtopyes=='y') veltops[ip+n2]=vth[n1];
                                	if(vbotyes=='y') velbots[ip+n2]=vbh[n1];
                                	velavgs[ip+n2] = vah[n1];
                                	dvdzs[ip+n2] = dvdzh[n1];
				}
			} else {
				/* not a good idea to use spline interpolation
					here */
				/*
				if(nvs>2) {
					z21 = 0.;
                                        z2n = 0.;
                                        spline_(xvh+n1,vth+n1,
						&nvs,&z21,&z2n,vt2,vb2);
                                        spline_(xvh+n1,vbh+n1,
						&nvs,&z21,&z2n,vb2);
                                        spline_(xvh+n1,vah+n1,
						&nvs,&z21,&z2n,va2);
                                        spline_(xvh+n1,dvdzh+n1,
						&nvs,&z21,&z2n,dvdz2);
				} 
				*/
				for(ip=0;ip<npicks[ih];ip++) {
					x = xpicks[ip+ih*maxpks];
					bisear_(&nvs,&one,xvh+n1,&x,&i2);
					if(x<=xvh[n1]) {
						if(vtopyes=='y') 
						veltops[ip+n2] = vth[n1];
						if(vbotyes=='y')
                                		velbots[ip+n2] = vbh[n1];
                                		velavgs[ip+n2] = vah[n1];
                                		dvdzs[ip+n2] = dvdzh[n1];
					} else if(x>=xvh[n1+nvs-1]) {
						if(vtopyes=='y') 
						veltops[ip+n2] = vth[n1+nvs-1];
						if(vbotyes=='y')
                                		velbots[ip+n2] = vbh[n1+nvs-1];
                                		velavgs[ip+n2] = vah[n1+nvs-1];
                                		dvdzs[ip+n2] = dvdzh[n1+nvs-1];
					} else { 
					/*} else if(nvs==2) {  */
						tmp = (x-xvh[n1+i2-1])/
						      (xvh[n1+i2]-xvh[n1+i2-1]);
						if(vtopyes=='y') 
						veltops[ip+n2] = vth[n1+i2-1] +
						  tmp*(vth[n1+i2]-vth[n1+i2-1]);
						if(vbotyes=='y')
						velbots[ip+n2] = vbh[n1+i2-1] +
						  tmp*(vbh[n1+i2]-vbh[n1+i2-1]);
                                		velavgs[ip+n2] = vah[n1+i2-1] +
						  tmp*(vah[n1+i2]-vah[n1+i2-1]);
                                		dvdzs[ip+n2] = dvdzh[n1+i2-1] +
						  tmp*
						  (dvdzh[n1+i2]-dvdzh[n1+i2-1]);
					/*
					} else {
						if(vtopyes=='y') 
                                        	splint_(xvh+n1,vth+n1,vt2,&nvs,
							&x,&veltops[ip+n2]);
						if(vbotyes=='y')
                                        	splint_(xvh+n1,vbh+n1,vb2,&nvs,
							&x,&velbots[ip+n2]);
                                        	splint_(xvh+n1,vah+n1,va2,&nvs,
							&x,&velavgs[ip+n2]);
                                        	splint_(xvh+n1,dvdzh+n1,dvdz2,
							&nvs,&x,&dvdzs[ip+n2]);
					*/
					}
				}
			}
			/* add velan locations to horizon picks */
			j2 = nvs + npicks[ih];
			for(i2=0;i2<npicks[ih];i2++) { 
				xnew[i2] = xpicks[i2+ih*maxpks];
				znew[i2] = zpicks[i2+ih*maxpks];
				vtnew[i2] = veltops[i2+ih*maxpks];
				vbnew[i2] = velbots[i2+ih*maxpks];
				vanew[i2] = velavgs[i2+ih*maxpks];
				gnew[i2] = dvdzs[i2+ih*maxpks];
			}
			for(i2=0;i2<nvs;i2++) { 
				xnew[i2+npicks[ih]] = xvh[n1+i2];
				znew[i2+npicks[ih]] = zvh[n1+i2];
				vtnew[i2+npicks[ih]] = vth[n1+i2];
				vbnew[i2+npicks[ih]] = vbh[n1+i2];
				vanew[i2+npicks[ih]] = vah[n1+i2];
				gnew[i2+npicks[ih]] = dvdzh[n1+i2];
			}

			if(vtopyes=='n') {
				i2 = npicks[ih];
				lin1d_(xpicks+ih*maxpks,veltops+ih*maxpks,&i2,
					xnew,vtnew,&j2,indx);
			}
			if(vbotyes=='n') {
				i2 = npicks[ih];
				lin1d_(xpicks+ih*maxpks,velbots+ih*maxpks,&i2,
					xnew,vbnew,&j2,indx);
			}
			for(i2=0;i2<j2;i2++) indx[i2] =  i2;
			qkisort(j2,xnew,indx);
			xpre = xnew[indx[0]] - 999.;
			j1 = 0;
			for(i2=0;i2<j2;i2++) {
				inow = indx[i2];
			   	xnow = xnew[inow];
				i1 = j1 + ih*maxpks;
			   	if(abs(xnow-xpre)>=1.) {
					xpicks[i1] = xnow;
					zpicks[i1] = znew[inow];
					veltops[i1] = vtnew[inow];
					velbots[i1] = vbnew[inow];
					velavgs[i1] = vanew[inow];
					dvdzs[i1] = gnew[inow];
					xpre = xnow;
					j1 = j1 + 1;
				}
			}
			if(j1>maxpks) 
			err("too many picks!! npick=%d at horizon %d",j1,ih+1); 
			npicks[ih] = j1;
		}

		free(indx);
		free(xnew);
		free(znew);
		free(vtnew);
		free(vbnew);
		free(vanew);
		free(gnew);

		free(zvs);
		free(ivs);
		free(isort);
		free(zsort);
		free(nvh);
		free(xvh);
		free(zvh);
		free(vbh);
		free(vth);
		free(vah);
		free(dvdzh);
		free(zgrid);
		free(vagrid);
		free(dvdzgrid);
		free(g1);
		free(g2);

		free(z2);
		/*
		free(vt2);
		free(vb2);
		free(va2);
		free(dvdz2);
		*/
	}
}

/* interface velocity edit */ 
void ivedit(float *xpick, float *vpick, char *hcolor, int nvi, char *title,
		float xmin, float dcdp, int cdpxmin, int *isave) {

	Display *dpy1;
	Window win1;
	XEvent event;
	KeySym keysym;

	GC gci, gci2;
        XGCValues *values;
        XColor scolor,ecolor;
        XWindowAttributes wa;
        Colormap cmap;
	int scr;
	XComposeStatus keystat;
	char keybuf[256];

	unsigned long black, white;

	int xbox, ybox, wbox, hbox, showloc=0;
	int xx1, yy1, xx2, yy2, i, n1tic, n2tic;
	float d1num, d2num, f1num, f2num, x1begb, x2begb, x1endb, x2endb;
	float vmin, vmax;
	int x, y, width, height, style ;
	int grid1=SOLID, grid2=SOLID,ip;
	float x1,x2,dis,dismin;
	float *vps;
	float cdpbegb, cdpendb;


	char *label1="cdp",*label2="velocity",
		*labelfont="Erg14",*titlefont="Rom22",
		*labelcolor="blue",*titlecolor="red",
		*g