iipick.c

seismic software,very useful

	sprintf(string,"(%0.4g,%0.4g)",x1,x2);
	XTextExtents(fs,string,strlen(string),&dummy,&dummy,&dummy,&overall);
	XDrawImageString(dpy,win,gc,xoffset,yoffset+overall.ascent,
		string,strlen(string));
}


void xMousePrint(Display *dpy, Window win, XEvent event, int style,
	FILE *mpicksfp, int x, int y, int width, int height,
	float x1begb, float x1endb, float x2begb, float x2endb,
	float *x1picks, float *x2picks, float *veltop, float *velbot, int *dif,
	float *velavg, float *dvdz,
	int *npicks, int nhs, int ism, 
	float xmin, float xmax, float zmin, float zmax,
	int cdpxmin, int cdpxmax, float dcdp, 
	int *hnums, char hnames[maxhs][80],
	char *dunits, char *zattrib, int idif, int ihfile) {


	int ic,ip,np,ih,jp,nhout,h0; 
	int p1,p2,i,ismx,ismxh, izt;
	float vmax, vmin, s;
	float *w, *f;
	float vt, vb, g;



	/* lateral smooth velocities */
	if(ism==0) {
	/* take average */ 
		for(ih=0;ih<nhs;ih++) {
			vt = 0.;
			vb = 0.;
			g = 0.;
			np = npicks[ih];
			for(ip=0;ip<np;ip++) {
				vt += veltop[ip+ih*maxpks];
				vb += velbot[ip+ih*maxpks];
				g += dvdz[ip+ih*maxpks];
			}
			if(np>0) {
				vt = vt/np;
				vb = vb/np;
				g = g/np;
			}
			for(ip=0;ip<np;ip++) {
				veltop[ip+ih*maxpks] = vt;
				velbot[ip+ih*maxpks] = vb;
				dvdz[ip+ih*maxpks] = g;
			}
		}
	} else if (ism>1) {
	/* smoothing */
		for(ih=0;ih<nhs;ih++) {
			h0 = ih * maxpks;
			np = npicks[ih];
			if(np<2) continue;
			ismx = ism;
			ismx = (ismx/2)*2+1;
			if(ismx>np) ismx=(np/2)*2+1;
			f = (float *) malloc(ismx*sizeof(float));
			w = (float *) malloc(np*sizeof(float));
			/* triangular smoothing coefficients */
			ismxh = (ismx+1)/2;
			s = 0.;
			for(i=0;i<ismx;i++) {
				f[i] = ismxh - abs(i+1-ismxh);
				s += f[i];
			}
			s = 1./s;
			for(i=0;i<ismx;i++) {
				f[i]=f[i]*s;
			}
			/* apply smoothing to veltop, velbot and dvdz */
			ic = 1;
			smth1d_(veltop+h0,f,w,&np,&ismx,&ic);
			for(ip=0;ip<np;ip++) veltop[ip+h0] = w[ip];
			smth1d_(velbot+h0,f,w,&np,&ismx,&ic);
			for(ip=0;ip<np;ip++) velbot[ip+h0] = w[ip];
			smth1d_(dvdz+h0,f,w,&np,&ismx,&ic);
			for(ip=0;ip<np;ip++) dvdz[ip+h0] = w[ip];

			free(w);
			free(f);
		}
	}

	/* find maximum and minimum velocities */
	vmin = veltop[0];
	vmax = veltop[0];
	for (ih=0;ih<nhs;ih++) {
		for(ip=0;ip<npicks[ih];ip++) {
			jp = ip + ih * maxpks;
			if(vmin>veltop[jp]) vmin = veltop[jp];
			if(vmin>velbot[jp]) vmin = velbot[jp];
			if(vmax<veltop[jp]) vmax = veltop[jp];
			if(vmax<velbot[jp]) vmax = velbot[jp];
		}
	}

	nhout = 0;
	for(ih=0;ih<nhs;ih++) {
		if(npicks[ih]>1) nhout = nhout + 1;
	}


	/* output picks */

	if(ihfile==0) {
		fprintf(mpicksfp, 
		"IFILE   <--- FILE FORMAT DESCRIPTOR (INTERFACE FILE) \n");
	} else {
		fprintf(mpicksfp, "HFILE  \n");
	}

	fprintf(mpicksfp,"\n");
	fprintf(mpicksfp, "MODEL.X-MIN         %8.1f\n",xmin);
	fprintf(mpicksfp, "MODEL.X-MAX         %8.1f\n",xmax);
	fprintf(mpicksfp, "MODEL.Z-MIN         %8.1f \n",zmin);
	fprintf(mpicksfp, "MODEL.Z-MAX         %8.1f \n",zmax);
	fprintf(mpicksfp, "MODEL.CDP#-AT-X-MIN %6d\n",cdpxmin);
	fprintf(mpicksfp, "MODEL.CDP#-AT-X-MAX %6d\n",cdpxmax);
	fprintf(mpicksfp, "MODEL.MIN-VELOCITY  %8.1f\n",vmin); 
	fprintf(mpicksfp, "MODEL.MAX-VELOCITY  %8.1f\n",vmax); 
	fprintf(mpicksfp, "MODEL.#HORIZONS     %6d\n",nhout); 
	fprintf(mpicksfp, "MODEL.VELOCITY-TYPE  INTERVAL\n"); 
	fprintf(mpicksfp, "MODEL.DIST-UNITS     %s\n",dunits); 
	fprintf(mpicksfp, "MODEL.Z-ATTRIBUTE    %s\n",zattrib); 
	fprintf(mpicksfp, "MODEL.CDP-SPACING   %8.1f\n",dcdp); 
	fprintf(mpicksfp, "\n"); 

	if( zattrib=="DEPTH SECTION" ) {
		izt = 1; 
	} else {
		izt = 0; 
	}

	
	for(ih=0;ih<nhs;ih++) {

		np = npicks[ih];
		if(np<=1) continue;

		fprintf(mpicksfp, "\n"); 
		fprintf(mpicksfp, "\n"); 
		fprintf(mpicksfp, "HORIZON-NUMBER %3d\n", hnums[ih]); 
		fprintf(mpicksfp, "HORIZON-NAME     %s\n", hnames[ih]); 
		fprintf(mpicksfp, "* \n"); 
		fprintf(mpicksfp, 
"* X-VALUE   Z-VALUE   VELTOP   VELBOT  DIF  DENS  QVAL  POIS   VELAVG   DVDZ \n");
		fprintf(mpicksfp, 
"* =======  ========  =======  =======  ===  ====  ====  ====  =======  ======= \n");

		h0 = ih*maxpks;
		for(ip=1;ip<np-1;ip++) dif[h0+ip]=idif;
		dif[h0] = 1;
		dif[h0+np-1] = 1;
		
		for (ip=0;ip<npicks[ih];ip++) {

			if(izt==1) {
   			fprintf(mpicksfp, 
"%9.1f %9.1f %8.0f %8.0f   %1d   0.00    0.  .000 %8.0f %8.0f \n",
				x2picks[ip+h0],x1picks[ip+h0],
				veltop[ip+h0],velbot[ip+h0],dif[ip+h0],
				velavg[ip+h0],dvdz[ip+h0]);
			} else {
   			fprintf(mpicksfp, 
"%9.1f %9.3f %8.0f %8.0f   %1d   0.00    0.  .000 %8.0f %8.0f \n",
				x2picks[ip+h0],x1picks[ip+h0],
				veltop[ip+h0],velbot[ip+h0],dif[ip+h0],
				velavg[ip+h0],dvdz[ip+h0]);
			}
		}
	}
	fflush(mpicksfp);
}

void xMousePicks(Display *dpy, Window win, XEvent event, int style,
	char hcolors[maxhs][10],
	FILE *mpicksfp, int x, int y, int width, int height,
	float x1begb, float x1endb, float x2begb, float x2endb,
	float *x1picks, float *x2picks, float *veltop, float *velbot, int *dif,
	int *npicks, int nhs, int ih, GC gc, int pkey, int *savebg)
{
	float x1,x2;
	int ipicks;
	int xx1,yy1,xx2,yy2,temp,ip,jp;
	int dismin,dis;
	int ipmin, ipins, offset, jh, itmp;
	char hnum[3];

	GC gchi;
        XGCValues *values;
        XColor scolor,ecolor;
        XWindowAttributes wa;
        Colormap cmap;
        int scr;

	/* save bitmap of window to background for retrival */ 
	if(*savebg) {
		fg2bg(dpy,win);
		*savebg = 0;
	}

	/* convert mouse location to (x1,x2) coordinates */
	if (style==NORMAL) {
		x1 = x1begb+(x1endb-x1begb)*(event.xmotion.x-x)/width;
		x2 = x2endb+(x2begb-x2endb)*(event.xmotion.y-y)/height;
	} else {
		x1 = x1begb+(x1endb-x1begb)*(event.xmotion.y-y)/height;
		x2 = x2begb+(x2endb-x2begb)*(event.xmotion.x-x)/width;
	}

	offset = (ih-1)*maxpks;
	jh = ih - 1;
	/* save to x1picks and x2picks */
	ipicks = npicks[jh];
	if( pkey == 0 ) {
		jp = ipicks + offset;
		x1picks[jp] = x1;
		x2picks[jp] = x2;
		npicks[jh] += 1;
	/* delete nearest picks from x1picks and x2picks*/
	} else if ( pkey==1 && ipicks > 0 ) {
		dismin=(x1-x1picks[offset])*(x1-x1picks[offset])+
		       (x2-x2picks[offset])*(x2-x2picks[offset]);
		ipmin = 0;
		for(ip=1;ip<ipicks;ip++) {
			jp = ip + offset;
		   	dis=(x1-x1picks[jp])*(x1-x1picks[jp])+
		       		(x2-x2picks[jp])*(x2-x2picks[jp]);
			if(dis<dismin) {
				dismin = dis;
				ipmin = ip;
			}
		}
		ipicks = ipicks - 1;
		for(ip=ipmin;ip<ipicks;ip++) {
			jp = ip + offset;
			x1picks[jp] = x1picks[jp+1];
			x2picks[jp] = x2picks[jp+1];
		}
		npicks[jh] = ipicks;
	/* insert the pick */
	} else if ( pkey==2 ) {
		dismin=(x2-x2picks[offset])*(x2-x2picks[offset]);
		ipmin = 0;
		for(ip=1;ip<ipicks;ip++) {
			jp = ip + offset;
			dis=(x2-x2picks[jp])*(x2-x2picks[jp]);
			if(dis<dismin) {
				dismin = dis;
				ipmin = ip;
			}
		}
		if(x2>x2picks[ipmin+offset]) { 
			ipins = ipmin+1;
			ipicks +=1;
			for(ip=ipicks;ip>ipins;ip--) {
				jp = ip+offset;
				x1picks[jp] = x1picks[jp-1];
				x2picks[jp] = x2picks[jp-1];
			}
			x1picks[ipins+offset] = x1;
			x2picks[ipins+offset] = x2;
		} else if(x2<x2picks[ipmin+offset]) {
			ipins = ipmin;
			ipicks +=1;
			for(ip=ipicks;ip>ipins;ip--) {
				jp = ip + offset;
				x1picks[jp] = x1picks[jp-1];
				x2picks[jp] = x2picks[jp-1];
			}
			x1picks[ipins+offset] = x1;
			x2picks[ipins+offset] = x2;
		}
		npicks[jh] = ipicks;
	}
	if ( pkey==5 ) {
		for(ip=0;ip<npicks[jh];ip++) {
			jp = ip + offset;
			x1picks[jp] = 0.;
			x2picks[jp] = 0.;
		}
		npicks[jh] = 0;
		ipicks = 0;
	}

       	XClearWindow(dpy, win);

	/* draw lines between picks */
	for(jh=0;jh<nhs;jh++) {
	
		/* get screen */
                scr = DefaultScreen(dpy);
                /* determine window's current colormap */
                XGetWindowAttributes(dpy,win,&wa);
                cmap = wa.colormap;
		/* create graphics contexts */
                gchi = XCreateGC(dpy,win,0,values);
                if (XAllocNamedColor(dpy,cmap,hcolors[jh],&scolor,&ecolor))
                        XSetForeground(dpy,gchi,ecolor.pixel);
                else
                        XSetForeground(dpy,gchi,1L);
                XSetLineAttributes(dpy,gchi,2,LineSolid,CapButt,JoinMiter);
		ipicks = npicks[jh];
		offset = jh*maxpks;
		itmp = ipicks / 5;
		if(itmp<2) itmp=ipicks/2;
		if(ipicks>1) {
			for(ip=1;ip<ipicks;ip++) {
				jp = ip + offset;
				/*
				if(x1picks[jp-1]<x1begb ||
				   x1picks[jp-1]>x1endb ||
				   x1picks[jp]<x1begb ||
				   x1picks[jp]>x1endb ||
				   x2picks[jp-1]<x2begb ||
				   x2picks[jp-1]>x2endb ||
				   x2picks[jp]<x2begb ||
				   x2picks[jp]>x2endb) continue;
				*/

				if (style==NORMAL) {
			xx1=(x1picks[jp-1]-x1begb)/(x1endb-x1begb)*width+x;
			yy1=(x2picks[jp-1]-x2endb)/(x2begb-x2endb)*height+y;
			xx2=(x1picks[jp]-x1begb)/(x1endb-x1begb)*width+x;
			yy2=(x2picks[jp]-x2endb)/(x2begb-x2endb)*height+y;
				} else {
			yy1=(x1picks[jp-1]-x1begb)/(x1endb-x1begb)*height+y;
			xx1=(x2picks[jp-1]-x2begb)/(x2endb-x2begb)*width+x;
			yy2=(x1picks[jp]-x1begb)/(x1endb-x1begb)*height+y;
			xx2=(x2picks[jp]-x2begb)/(x2endb-x2begb)*width+x;
				}
				XDrawLine(dpy,win,gchi,xx1,yy1,xx2,yy2);
				if(ip==1) {
					XDrawLine(dpy,win,gchi,
						xx1-4,yy1,xx1+4,yy1);
					XDrawLine(dpy,win,gchi,
						xx1,yy1-4,xx1,yy1+4);
				}
				XDrawLine(dpy,win,gchi,xx2-4,yy2,xx2+8,yy2);
				XDrawLine(dpy,win,gchi,xx2,yy2-4,xx2,yy2+4);
				if( (ip%itmp)==0 ) {
					sprintf(hnum,"%-3d",jh+1);
					XDrawString(dpy,win,gchi,
                        		xx2,yy2,hnum,strlen(hnum));
				} 
			}
			sprintf(hnum,"%-3d",jh+1);
			XDrawString(dpy,win,gchi,
                        xx2,yy2,hnum,strlen(hnum));
			
			XFreeGC(dpy,gchi);
		}
	}
}


/* compute velocities from VELO and DVDZ cards */
void vtopbot(int maxp,int maxcdp, 
     int ncdp_velo,float *x_velo,float *z_velo,float *v_velo,int *nps_velo,
     int ncdp_dvdz,float *x_dvdz,float *zt_dvdz,float *zb_dvdz, 
     float *dvdz_dvdz, int *nps_dvdz,
     int nhs,float *xpicks,float *zpicks,int *npicks,float xtol, 
     float *veltops,float *velbots,float *velavgs,float *dvdzs,int lhnvelo) {

	float x, z, va, g, vat, gt, vab, gb;
	int ih, ip, jh;
	float zmax, xmin, xmax, zt, zb, tmp;
	float *vth, *vbh, *vah, *dvdzh, *xvh, *zvs, *zvh;
	float *zgrid, *vagrid, *dvdzgrid, dz, *g1, *g2, *z2, z21, z2n;
	/* float *vt2, *vb2, *va2, *dvdz2; */
	int i1,i2,n1,n2,one,nzgrid;
	int *ivs, *nvh;
	int iz, n, j2, j1, nvs;
	float *xnew, *znew, *vanew, *vtnew, *vbnew, *gnew, xpre, xnow;
	int *indx, inow;

	float *zsort;
	int *isort;

	FILE *tty;
        char vtopyes,vbotyes;
	
	
	if(ncdp_velo==0) {
		/* if no velo card input, simply return with warning message */
		warn("no VELO card input, ifile editing may be needed ! \n");
	} else {
		/* find maximum depth */
		zmax = 0.;
		for(ih=0;ih<nhs;ih++) {
			for(iz=0;iz<npicks[ih];iz++) {
				if(zmax<zpicks[ih*maxpks+iz])
				  	zmax=zpicks[ih*maxpks+iz];
			}
		}


		for(i2=0;i2<ncdp_velo;i2++) {
			for(iz=0;iz<nps_velo[i2];iz++) {
				if(zmax<z_velo[i2*maxp+iz])
				  	zmax=z_velo[i2*maxp+iz];
			}
		}




		nzgrid = 1024;
		zgrid = (float *) malloc(nzgrid*sizeof(float));
		vagrid = (float *) malloc(nzgrid*ncdp_velo*sizeof(float));
		dvdzgrid = (float *) malloc(nzgrid*ncdp_velo*sizeof(float));
		g1 = (float *) malloc(nzgrid*sizeof(float));
		g2 = (float *) malloc(nzgrid*sizeof(float));

		vth = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		vbh = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		vah = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		dvdzh = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		xvh = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		zvh = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		nvh = (int *) malloc(nhs*sizeof(int));

		zvs = (float *) malloc(nhs*sizeof(float));
		ivs = (int *) malloc(nhs*sizeof(int));
		zsort = (float *) malloc(nhs*sizeof(float));
		isort = (int *) malloc(nhs*sizeof(int));

		vtnew = (float *) malloc(maxpks*sizeof(float));
		vbnew = (float *) malloc(maxpks*sizeof(float));
		indx = (int *) malloc(maxpks*sizeof(int));
		vanew = (float *) malloc(maxpks*sizeof(float));
		gnew = (float *) malloc(maxpks*sizeof(float));
		xnew = (float *) malloc(maxpks*sizeof(float));
		znew = (float *) malloc(maxpks*sizeof(float));