vipick.c

seismic software,very useful

			XDrawLine(dpy,win,gcvip,xx1,yy1,xx2,yy2);

			yy1=(zp[ip-1]-x1begb)/(x1endb-x1begb)*height+y;
			xx1=(vp[ip-1]-x2begb)/(x2endb-x2begb)*width+x;
			yy2=(zp[ip]-x1begb)/(x1endb-x1begb)*height+y;
			xx2=(vp[ip]-x2begb)/(x2endb-x2begb)*width+x;
			XDrawLine(dpy,win,gcvp,xx1,yy1,xx2,yy2);
		}
		/* free resources before returning */
                XFreeGC(dpy,gcvp);
                XFreeGC(dpy,gcvip);

	}

	/* reference velocity curve */
	if(nr>1) {
		/* create graphics contexts */
                gcvr = XCreateGC(dpy,win,0,values);
                if (XAllocNamedColor(dpy,cmap,vrefcolor,&scolor,&ecolor))
                        XSetForeground(dpy,gcvr,ecolor.pixel);
                else
                        XSetForeground(dpy,gcvr,1L);
                XSetLineAttributes(dpy,gcvr,2,LineSolid,CapButt,JoinMiter);

		for(ip=1;ip<nr;ip++) {
			yy1=(zr[ip-1]-x1begb)/(x1endb-x1begb)*height+y;
			xx1=(vr[ip-1]-x2begb)/(x2endb-x2begb)*width+x;
			yy2=(zr[ip]-x1begb)/(x1endb-x1begb)*height+y;
			xx2=(vr[ip]-x2begb)/(x2endb-x2begb)*width+x;
			XDrawLine(dpy,win,gcvr,xx1,yy1,xx2,yy2);
		}
		/* free resources before returning */
                XFreeGC(dpy,gcvr);

	}

	/* previous picks read from mpicks */
	if(npcrd>1) {
		/* create graphics contexts */
                gcvpcrd = XCreateGC(dpy,win,0,values);
                if (XAllocNamedColor(dpy,cmap,vpcrdcolor,&scolor,&ecolor))
                        XSetForeground(dpy,gcvpcrd,ecolor.pixel);
                else
                        XSetForeground(dpy,gcvpcrd,1L);
                XSetLineAttributes(dpy,gcvpcrd,2,LineSolid,CapButt,JoinMiter);

		gcvip = XCreateGC(dpy,win,0,values);
                if (XAllocNamedColor(dpy,cmap,vpcrdcolor,&scolor,&ecolor))
                        XSetForeground(dpy,gcvip,ecolor.pixel);
                else
                        XSetForeground(dpy,gcvip,1L);
                XSetLineAttributes(dpy,gcvip,1,LineSolid,CapButt,JoinMiter);
                xw = width /25;
                xh = height /25;

                yy1=(zpcrd[0]-x1begb)/(x1endb-x1begb)*height+y-xh/2;
                xx1=(vpcrd[0]-x2begb)/(x2endb-x2begb)*width+x;
                yy2=yy1+xh;
                xx2=xx1;
                XDrawLine(dpy,win,gcvip,xx1,yy1,xx2,yy2);
                yy1=(zpcrd[0]-x1begb)/(x1endb-x1begb)*height+y;
                xx1=(vpcrd[0]-x2begb)/(x2endb-x2begb)*width+x-xw/2;
                yy2=yy1;
                xx2=xx1+xw;
                XDrawLine(dpy,win,gcvip,xx1,yy1,xx2,yy2);


		for(ip=1;ip<npcrd;ip++) {
			yy1=(zpcrd[ip]-x1begb)/(x1endb-x1begb)*height+y-xh/2;
                        xx1=(vpcrd[ip]-x2begb)/(x2endb-x2begb)*width+x;
                        yy2=yy1+xh;
                        xx2=xx1;
                        XDrawLine(dpy,win,gcvip,xx1,yy1,xx2,yy2);
                        yy1=(zpcrd[ip]-x1begb)/(x1endb-x1begb)*height+y;
                        xx1=(vpcrd[ip]-x2begb)/(x2endb-x2begb)*width+x-xw/2;
                        yy2=yy1;
                        xx2=xx1+xw;
                        XDrawLine(dpy,win,gcvip,xx1,yy1,xx2,yy2);


			yy1=(zpcrd[ip-1]-x1begb)/(x1endb-x1begb)*height+y;
			xx1=(vpcrd[ip-1]-x2begb)/(x2endb-x2begb)*width+x;
			yy2=(zpcrd[ip]-x1begb)/(x1endb-x1begb)*height+y;
			xx2=(vpcrd[ip]-x2begb)/(x2endb-x2begb)*width+x;
			XDrawLine(dpy,win,gcvpcrd,xx1,yy1,xx2,yy2);
		}
		/* free resources before returning */
                XFreeGC(dpy,gcvpcrd);
                XFreeGC(dpy,gcvip);

	}

	/* interval velocity curve of picked result */
	if(ppkey==4 && *np > 1 ) viplot(dpy, win, vpickcolor,  
					x, y, width, height,
					x1begb, x1endb, x2begb, x2endb,
					zp, vp, *np, zr, nr, ivmode, 
					ztop, zbot, v0, dvdz,
					vatop, vabot, *ia, ivof);

	/* interval velocity curve of previous picked result */
	if(ppkey==4 && npcrd > 1 ) viplot(dpy, win, vpcrdcolor,  
					x, y, width, height,
					x1begb, x1endb, x2begb, x2endb,
					zpcrd, vpcrd, npcrd, zr, nr, ivmode, 
					ztop, zbot, v0, dvdz,
					vatop, vabot, *ia, ivof);

	/* interval velocity curve of reference velocity input */
	if(ppkey==4 && nr > 1 ) {
		/* create graphics contexts */
                gcvir = XCreateGC(dpy,win,0,values);
                if (XAllocNamedColor(dpy,cmap,vrefcolor,&scolor,&ecolor))
                        XSetForeground(dpy,gcvir,ecolor.pixel);
                else
                        XSetForeground(dpy,gcvir,1L);
                XSetLineAttributes(dpy,gcvir,2,LineOnOffDash,CapButt,JoinMiter);

		vi = (float *) malloc(nr*sizeof(float));

		vi[0] = vr[0];
		for(ip=1;ip<nr;ip++) {
		   	vi[ip] = (zr[ip]*vr[ip]-zr[ip-1]*vr[ip-1])
					/(zr[ip]-zr[ip-1]);
	   	} 
		for(ip=1;ip<nr;ip++) {
			yy1=(zr[ip-1]-x1begb)/(x1endb-x1begb)*height+y;
			xx1=(vi[ip-1]-x2begb)/(x2endb-x2begb)*width+x;
			yy2=(zr[ip]-x1begb)/(x1endb-x1begb)*height+y;
			xx2=(vi[ip]-x2begb)/(x2endb-x2begb)*width+x;
			XDrawLine(dpy,win,gcvir,xx1,yy1,xx2,yy2);
		}
		/* free resources before returning */
                XFreeGC(dpy,gcvir);
		free(vi);
	}
	/* delete dvdz analysis point */
	if( ppkey==6 && *ia > 0 ) {
		ipmin = 0;
		for(ip=0;ip<*ia;ip++) {
			if( zptop[ip] <= x1 && x1 <= zpbot[ip] ) {
				ipmin = ip;
				break;
			}
		}
		*ia= *ia - 1;
		for(ip=ipmin;ip<*ia;ip++) {
			zptop[ip] = zptop[ip+1];
			zpbot[ip] = zpbot[ip+1];
			vatop[ip] = vatop[ip+1];
			vabot[ip] = vabot[ip+1];
			dvdz[ip] = dvdz[ip+1];
			v0[ip] = v0[ip+1];
		}
	}
	/* velocity gradient picking plot */
	if( (ppkey==5 || ppkey==3 || ppkey==6 ) && *ia > 0 ) {
		zz = (float*) malloc(maxpks*sizeof(float));
		vaa = (float*) malloc(maxpks*sizeof(float));
		for(ip=0;ip<*ia;ip++) {
			zandva(zptop[ip],zpbot[ip],zp,vp,*np,zz,vaa,&nn);
			if(nn>1) {
				ztop[ip] = zz[0];
				zbot[ip] = zz[nn-1];
				vatop[ip] = vaa[0];
				vabot[ip] = vaa[nn-1];
				estvi_(zz,vaa,&nn,&ibfit,&v00,&aa);
				dvdz[ip] = aa;
				v0[ip] = v00;
			} else {
				ztop[ip] = x1begb - 1000. ;
				zbot[ip] = x1begb - 1000.;
				vatop[ip] = vp[0];
				vabot[ip] = vp[0];
				dvdz[ip] = 0.;
				v0[ip] = vp[0];
			}
		}
		free(zz);
		free(vaa);
		/* create graphics contexts */
                gcvip = XCreateGC(dpy,win,0,values);
                if (XAllocNamedColor(dpy,cmap,vpickcolor,&scolor,&ecolor))
                        XSetForeground(dpy,gcvip,ecolor.pixel);
                else
                        XSetForeground(dpy,gcvip,1L);
                XSetLineAttributes(dpy,gcvip,2,LineOnOffDash,CapButt,JoinMiter);

                gcvir = XCreateGC(dpy,win,0,values);
                if (XAllocNamedColor(dpy,cmap,vpickcolor,&scolor,&ecolor))
                        XSetForeground(dpy,gcvir,ecolor.pixel);
                else
                        XSetForeground(dpy,gcvir,1L);
                XSetLineAttributes(dpy,gcvir,1,LineSolid,CapButt,JoinMiter);

		for(ip=0;ip<*ia;ip++) {
			if(ztop[ip] < zbot[ip]) {
				yy1=(ztop[ip]-x1begb)/(x1endb-x1begb)*height+y;
				xx1=(vatop[ip]-x2begb)/(x2endb-x2begb)*width+x;
				yy2=yy1;
				xx2=width+x;
				XDrawLine(dpy,win,gcvir,xx1,yy1,xx2,yy2);

				yy1=(zbot[ip]-x1begb)/(x1endb-x1begb)*height+y;
				xx1=(vabot[ip]-x2begb)/(x2endb-x2begb)*width+x;
				yy2=yy1;
				xx2=width+x;
				XDrawLine(dpy,win,gcvir,xx1,yy1,xx2,yy2);

				yy1=(ztop[ip]-x1begb)/(x1endb-x1begb)*height+y;
				xx1=(v0[ip]-x2begb)/(x2endb-x2begb)*width+x;
				yy2=(zbot[ip]-x1begb)/(x1endb-x1begb)*height+y;
				xx2=(v0[ip]+dvdz[ip]*(zbot[ip]-ztop[ip])-x2begb)
						/(x2endb-x2begb)*width+x;
				XDrawLine(dpy,win,gcvip,xx1,yy1,xx2,yy2);
			}
		}
		/* free resources before returning */
                XFreeGC(dpy,gcvip);
                XFreeGC(dpy,gcvir);
	}
}

/* plot interval velocity curve of picked result */
void viplot(Display *dpy, Window win, char *vicolor,  
	int x, int y, int width, int height,
	float x1begb, float x1endb, float x2begb, float x2endb,
	float *zp, float *vp, int np, float *zr, int nr, int ivmode,
	float *ztop, float *zbot, float *v0, float *dvdz,
	float *vatop, float *vabot, int ia, int ivof) {
	
	int ni, *indx, xx1, xx2, yy1, yy2;
	float *va, *vi, *zi;
	int ip; 
	float z1,z2,v1,v2,aa;

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

        /* get screen */
        scr = DefaultScreen(dpy);
        /* determine window's current colormap */
        XGetWindowAttributes(dpy,win,&wa);
        cmap = wa.colormap;

	/* create graphics contexts */
       	gci = XCreateGC(dpy,win,0,values);
       	if (XAllocNamedColor(dpy,cmap,vicolor,&scolor,&ecolor))
       		XSetForeground(dpy,gci,ecolor.pixel);
       	else
        	XSetForeground(dpy,gci,1L);
        XSetLineAttributes(dpy,gci,2,LineOnOffDash,CapButt,JoinMiter);


	if (ivmode==1 && ivof==0) {
		if( nr > 1 && ia == 0 ) {
			ni = nr;
			va = (float *) malloc(ni*sizeof(float));
			vi = (float *) malloc(ni*sizeof(float));
			zi = (float *) malloc(ni*sizeof(float));
			indx = (int *) malloc(ni*sizeof(int));

			for(ip=0;ip<ni;ip++) zi[ip]=zr[ip];

			lin1d_(zp,vp,&np,zi,va,&ni,indx);

			vi[0] = va[0];
			for(ip=1;ip<ni;ip++) {
				if(zi[ip]<=zp[np-1]) {  
			   		vi[ip] = (zi[ip]*va[ip]
						-zi[ip-1]*va[ip-1])
							/(zi[ip]-zi[ip-1]);
		   		} else {
					vi[ip] = vi[ip-1];
		   		}
			}
			for(ip=0;ip<ni;ip++) {
				va[ip] = vi[ip];
			}
		} else {
			ni = np;
			va = (float *) malloc(ni*sizeof(float));
			vi = (float *) malloc(ni*sizeof(float));
			zi = (float *) malloc(ni*sizeof(float));
			for(ip=0;ip<ni;ip++) {
				zi[ip] = zp[ip];
				xx1 = -1;
				aa = 0.;
				for(yy1=0;yy1<ia;yy1++) {
					if(zi[ip]>ztop[yy1] &&	
					   zi[ip]<=zbot[yy1]) {
						xx1 = yy1;
						break;
					} 	
				}
				if (xx1 != -1) aa = dvdz[xx1];
				if(v0[xx1]>0.) {
					vi[ip] = v0[xx1]+(zp[ip]-ztop[xx1])*aa;
					xx2 = ip - 1;
					if(xx2<0) xx2 = 0;
					if(zi[xx2]>ztop[xx1] &&	
					   zi[xx2]<=zbot[xx1]) {
						va[ip] = v0[xx1]+
							(zp[xx2]-ztop[xx1])*aa;
					} else {
						va[ip] = v0[xx1];
					}
				} else {
					z2 = zp[ip];
					v2 = vp[ip];
					if(ip==0) {
						v1 = v2;
						z1 = 0.;
					} else {
						v1 = vp[ip-1];	
						z1 = zp[ip-1];
					}
					if(z2>z1) {
						va[ip]=(z2*v2-z1*v1)/(z2-z1) - 
							0.5*aa*(z2-z1);
					} else {
						va[ip] = v2;
					}
					vi[ip] = va[ip] + aa * (z2-z1); 
				}
			}
			
		}
	
	} else if(ivof==0 && ivmode!=1) {
		ni = np;
		va = (float *) malloc(ni*sizeof(float));
		vi = (float *) malloc(ni*sizeof(float));
		zi = (float *) malloc(ni*sizeof(float));

		vi[0] = vp[0];
		va[0] = vi[0];
		zi[0] = zp[0];
		for(ip=1;ip<ni;ip++) {
			zi[ip] = zp[ip];	
		   	vi[ip] = (zp[ip]*vp[ip]-zp[ip-1]*vp[ip-1])
					/(zp[ip]-zp[ip-1]);
			va[ip] = vi[ip];
		}
	} else if(ivof==1) {
		ni = np;
		va = (float *) malloc(ni*sizeof(float));
		vi = (float *) malloc(ni*sizeof(float));
		zi = (float *) malloc(ni*sizeof(float));

		vi[0] = vp[0];
		va[0] = vi[0];
		zi[0] = zp[0];
		for(ip=1;ip<ni;ip++) {
			zi[ip] = zp[ip];	
		   	vi[ip] = (zp[ip]*vp[ip]*vp[ip]
					-zp[ip-1]*vp[ip-1]*vp[ip-1])
					/(zp[ip]-zp[ip-1]);
			if(vi[ip]>0.) {
				vi[ip] = sqrt(vi[ip]);
			} else {
				vi[ip] = 0.;
			}
			va[ip] = vi[ip];
		}
	}

	
	for(ip=1;ip<ni-1;ip++) {
	/*
		fprintf(stderr,"draw point ip=%d \n",ip);
	*/

		yy1=(zi[ip-1]-x1begb)/(x1endb-x1begb)*height+y;
		xx1=(va[ip]-x2begb)/(x2endb-x2begb)*width+x;
		yy2=(zi[ip]-x1begb)/(x1endb-x1begb)*height+y;
		xx2=(vi[ip]-x2begb)/(x2endb-x2begb)*width+x;
		XDrawLine(dpy,win,gci,xx1,yy1,xx2,yy2);

		yy1=(zi[ip]-x1begb)/(x1endb-x1begb)*height+y;
		xx1=(vi[ip]-x2begb)/(x2endb-x2begb)*width+x;
		yy2=yy1;
		xx2=(va[ip+1]-x2begb)/(x2endb-x2begb)*width+x;
		XDrawLine(dpy,win,gci,xx1,yy1,xx2,yy2);
	}

	yy1=(zi[ni-2]-x1begb)/(x1endb-x1begb)*height+y;
	xx1=(va[ni-1]-x2begb)/(x2endb-x2begb)*width+x;
	yy2=(zi[ni-1]-x1begb)/(x1endb-x1begb)*height+y;
	xx2=(vi[ni-1]-x2begb)/(x2endb-x2begb)*width+x;
	XDrawLine(dpy,win,gci,xx1,yy1,xx2,yy2);

	if(zi[0]>x1begb) {
		yy1=y;
		xx1=(va[0]-x2begb)/(x2endb-x2begb)*width+x;
		yy2=(zi[0]-x1begb)/(x1endb-x1begb)*height+y;
		xx2=xx1;
		XDrawLine(dpy,win,gci,xx1,yy1,xx2,yy2);

		yy1=(zi[0]-x1begb)/(x1endb-x1begb)*height+y;
		xx1=(vi[0]-x2begb)/(x2endb-x2begb)*width+x;
		yy2=yy1;
		xx2=(va[1]-x2begb)/(x2endb-x2begb)*width+x;
		XDrawLine(dpy,win,gci,xx1,yy1,xx2,yy2);
	}
	if(zi[ni-1]<x1endb) {
		yy1=(zi[ni-1]-x1begb)/(x1endb-x1begb)*height+y;
		xx1=(vi[ni-1]-x2begb)/(x2endb-x2begb)*width+x;
		yy2=height+y;
		xx2=xx1;
		XDrawLine(dpy,win,gci,xx1,yy1,xx2,yy2);
	}
		

	/* free resources before returning */
        XFreeGC(dpy,gci);
	free(indx);
	free(vi);
	free(va);
	free(zi);

}