To create value bars for raster plots
This is probably superceded by Hector's bar feature. This one
gives the author nice control over the bar values.
#
# Bar generates a array that is used to show the scale of raster image.
#
# Bar n1=2 n2=2 minval=0. maxval=1. type=h > out.h
#
# FILES
# OUTPUT:
# stdout: scale data
#
# PARAMETERS
# FROM PAR:
# n1=2,n2=2 array size
# minval=0. minimum value
# maxval=1. maximum value
# type=h(v) horizontal bar
# vertical bar
# LIMITATIONS:
#
#
#%end of self doc
# ----------
# Keyword: scale raster
# ----------
integer n1,n2
from par: integer n1=2,n2=2
allocate: real array(n1,n2)
subroutine submain(n1,n2,array)
integer n1,n2, i1,i2
integer outfd,output
character*10 type
character CHAR
real d1,d2,d3,o1,o2,o3
real minval,maxval
real array(n1,n2)
from par: real minval=0.,maxval=1.
from par: character type='h'//CHAR(0)
if((type(1:1) != 'h')&(type(1:1) != 'v'))
{
write(0,*)'type=',type
call erexit('type option not implemented')
}
if(type(1:1) == 'h')
{
d1=1.
o1=0.
d2=(maxval-minval)/(n2-1)
o2=minval
}
else
{
d1=(maxval-minval)/(n1-1)
o1=minval
d2=1.
o2=0.
}
to history: integer n1,n2,n3:1
to history: real d1,d2,d3:1.
to history: real o1,o2,o3:0.
outfd=output()
call hclose()
#assign values to the array
if(type(1:1) == 'h')
do i2 = 1,n2
do i1 = 1,n1
array(i1,i2) = minval+(i2-1)*(maxval-minval)/(n2-1)
else
do i2 = 1,n2
do i1 = 1,n1
array(i1,i2) = minval+(i1-1)*(maxval-minval)/(n1-1)
call rite(outfd,array,n1*n2*4)
return
end
Finally we add the bar to our figure in our makerule:
ttmlay.v& FIGDIR/ttmlay.v : DATDIR/Strvel/strvel_lay.h DATDIR/Ttm/ttmlay.h \
DATDIR/Ttm/ttmlay_max.h
Taplot < DATDIR/Strvel/strvel_lay.h bias=3. clip=3. | \
Ta2vplot title="Traveltime" label1="Depth (km)" label2="Surface (km)" \
wheretitle=t TLSZFAT color=G d1num=1. d2num=1. o1num=0 o2num=0 grid=n \
polarity=1 xinch=10.24 yinch=5.12 out=back.v > Junk
Bar n2=100 type=h minval=1.499 maxval=4.501 | Taplot bias=3. clip=3.| \
Ta2vplot title=" " label2="Velocity (km/s)" wherexlabel=b TLSZFAT \
color=G wantaxis1=n d2num=1. o2num=1.5 grid=n polarity=1 xinch=10.24 \
yinch=0.5 out=junk.v > Junk
vppen junk.v yshift=-3 > bar.v
Contour nc=24 c0=0 dc=0.1 < DATDIR/Ttm/ttmlay.h ALL0_10 d1num=0 \
d2num=0 title=" " label1=" " label2=" " wantframe=n wantaxis=n pad=n \
xinch=10.24 yinch=5.12 out=junka.v > Junk
Contour nc=24 c0=0 dc=0.1 < DATDIR/Ttm/ttmlay.h ALL5_2 d1num=0 \
d2num=0 title=" " label1=" " label2=" " wantframe=y wantaxis=n pad=n \
xinch=10.24 yinch=5.12 out=junkb.v > Junk
vp_Overlay back.v bar.v junka.v junkb.v > junk1.v
Contour nc=30 c0=0 dc=0.1 < DATDIR/Ttm/ttmlay_max.h ALL0_10 d1num=0 \
d2num=0 title=" " label1=" " label2=" " wantframe=n wantaxis=n pad=n \
xinch=10.24 yinch=5.12 out=junka.v > Junk
Contour nc=30 c0=0 dc=0.1 < DATDIR/Ttm/ttmlay_max.h ALL4_2 d1num=0 \
d2num=0 title=" " label1=" " label2=" " wantframe=y wantaxis=n pad=n \
xinch=10.24 yinch=5.12 out=junkb.v > Junk
vp_Overlay back.v bar.v junka.v junkb.v > junk2.v
vp_Movie junk1.v junk2.v > FIGDIR/ttmlay.v
rm back.v junka.v junkb.v junk.v junk1.v junk2.v bar.v Junk
To draw a line animation in vplot
# Purpose:Draws raypath geometry for an ellipse generated by
# an impulse in time on a common-offset section
# Reference: Deregowski and Rocca (Classic DMO paper)
# Author: Jim Black
#
# ellipse.x par=parafile.p out=file.v
#
# FILES:
# INPUT:
# OUTPUT:
# stdout: vplot file
#
# PARAMETERS:
# FROM INPUT:
# FROM PAR:
# xll,yll,xur,yur the true coordinates of the vplot
# plotcol=6 color of plot
# plotfat=1 fat of plot
# wantframe=0 =0 no frame; =1 frame
# narray =41 number of points along ellipse
# h=1.0 half-offset
# y=0.0 midpoint
# time=1.6 two-way time (S to R to G)
# v=2.0 velocity
# option=1 =1 Ellipse, all S-R-G paths
# =2 Ellipse, one theta, zero-offset geometry
# =3 Ellipse, one theta, ellipse parametrization
# animate=0 =1 Animates ellipse
#
# NOTE:
#% end of self doc
#
#-------------
# keyword: dmo geometry ellipse
#-------------
integer narray
from par: integer narray=41
allocate: real xval(narray),zval(narray)
subroutine submain(xval,zval,narray)
integer narray,ith
real xval(narray), zval(narray)
integer plotcol,plotfat,fat
integer wantframe,option,animate
integer outfd, output, iflag, emph
real xll,yll,xur,yur
real xmin,zmin,xmax,zmax,slop
real xcenter,zcenter,xscale,zscale
real x,z,x1,z1,xs,xg,xr,zr,xrp,zrp,xsp,y,h,surface,d,time,v
real a,b,dth,thtmp
real xstart,zstart,xend,zend,xz
real theta,tanth,sinth,costh
character*50 out
# PARAMETERS
# INPUT
from par: integer plotcol=6,plotfat=1
from par: integer wantframe=0,option=1,animate=0
from par: real xll=1.705,yll=1.37,xur=11.945,yur=8.87
from par: real h=1.0,y=0.0,time=1.6,v=2.0,theta=30
# OUTPUT
# CHECK
# FILES
# INPUT
# OUTPUT
outfd = output()
if(getch('out','s',out) == 0)
call vpfilep(outfd)
else
call vpfile(out)
call hclose()
# CALCULATION
# set parameters
#
a = v * time /2.
b = sqrt(a*a-h*h)
xs = y - h
xg = y + h
surface = 0.
slop = 0.05*a
xmin = -2.
zmin = -1.
xmax = 2.
zmax = 2.
xcenter=0.5*(xmin+xmax)
zcenter=0.5*(zmin+zmax)
# set vplot scales
call vporig(xll,yur)
call vpuorig(xmin,zmin)
xscale = (xur-xll)/(xmax-xmin)
zscale = -(yur-yll)/(zmax-zmin)
call vpscale(xscale,zscale)
# draw and annotate surface
fat = 2*plotfat
call vpfat(fat)
call vpcolor(7)
x = xmin
z = surface
call vpumove(x,z)
x = xmax
z = surface
call vpudraw(x,z)
call vpcolor(5)
call vpfat(plotfat)
call vptjust(4,6)
call vputext(xs-0.5*slop,surface-slop,8,0,'S')
call vputext(xg,surface-slop,8,0,'G')
call vputext(y,surface-slop,8,0,'M')
dth = 180./(narray-1)
call vpcolor(4)
# LOOP on theta (ith) to generate ellipse
# ---------------------------------------
do ith = 1,narray {
thtmp = 90. - (ith-1) * dth
thtmp = thtmp*4.*atan(1.)/180.
tanth = tan(thtmp)
sinth = sin(thtmp)
costh = cos(thtmp)
d = sqrt(a*a*sinth*sinth + b*b*costh*costh)
xr = y - d*sinth - h*h*costh*costh*sinth/d
zr = surface + d*costh - h*h*costh*sinth*sinth/d
xval(ith) = xr
zval(ith) = zr
} # End theta loop
# Draw ellipse
call vpcolor(5)
call vpupline(xval,zval,narray)
if( option != 1 ) goto 500
# Loop on theta to draw SRG raypaths
# ----------------------------------
iflag = 0
emph = 3
do ith = 1,narray,5 {
iflag = iflag + 1
thtmp = 90. - (ith-1) * dth
thtmp = thtmp*4.*atan(1.)/180.
tanth = tan(thtmp)
sinth = sin(thtmp)
costh = cos(thtmp)
d = sqrt(a*a*sinth*sinth + b*b*costh*costh)
xr = y - d*sinth - h*h*costh*costh*sinth/d
zr = surface + d*costh - h*h*costh*sinth*sinth/d
# Shot to reflection point
call vpcolor(4)
x = xs
z = surface
x1 = xr
z1 = zr
if(iflag == emph) call vpcolor(6)
if(ith != 1 && ith != narray) call vpuarrow(x,z,x1,z1,0.02*zmax)
if(iflag == emph) call vpcolor(4)
# Annotate beta
xend = (xs + xr)/2.
zend = (surface + zr)/2.
if(iflag == emph) {
call vpcolor(6)
call vputext(xend-slop, zend+0.8*slop,8,0,'\F9 b')
call vpcolor(4)
}
# reflection point to geophone
x = xr
z = zr
call vpumove(x,z)
x1 = xg
z1 = surface
if(iflag == emph) call vpcolor(6)
if(ith != 1 && ith != narray) call vpuarrow(x,z,x1,z1,0.02*zmax)
if(iflag == emph) call vpcolor(6)
# Annotate alpha
xend = (xg + xr)/2.
zend = (surface + zr)/2.
if(iflag == emph) {
call vpcolor(6)
call vputext(xend+slop, zend+0.5*slop,8,0,'\F9 a')
call vpcolor(4)
}
if(animate == 1) {
call vperase() # clear screen, then redraw basic stuff
# redraw and annotate surface
fat = 2*plotfat
call vpfat(fat)
call vpcolor(7)
x = xmin
z = surface
call vpumove(x,z)
x = xmax
z = surface
call vpudraw(x,z)
call vpcolor(5)
call vpfat(plotfat)
call vptjust(4,6)
call vputext(xs-0.5*slop,surface-slop,8,0,'S')
call vputext(xg,surface-slop,8,0,'G')
call vputext(y,surface-slop,8,0,'M')
# redraw ellipse
call vpcolor(5)
call vpupline(xval,zval,narray)
} # end animate=1 block
} # end loop on theta
#
goto 1000 # That's it for option 1
500 continue
# Begin work for drawing dipping layer
# for options 2 and 3
# ------------------------------------
theta = theta*4.*atan(1.)/180.
tanth = tan(theta)
sinth = sin(theta)
costh = cos(theta)
d = sqrt(a*a*sinth*sinth + b*b*costh*costh)
xr = y - d*sinth - h*h*costh*costh*sinth/d
zr = surface + d*costh - h*h*costh*sinth*sinth/d
xz = y - h*h/d*sinth
xsp = y - h - 2*(d-h*sinth)*sinth
xrp = y - d*sinth
zrp = surface + d*costh
xstart = xsp
zstart = surface + d/costh + tanth * (xstart-y)
xend = 0.5*(xg + y)
zend = surface + d/costh + tanth * (xend-y)
# Draw dipping reflector
call vpcolor(2)
x=xstart
z=zstart
call vpumove(x,z)
x=xend
z=zend
call vpudraw(x,z)
call vpdash(0.,0.,0.,0.)
# indicate angle theta
call vpcolor(7)
call vpdash(.1,.1,.1,.1)
x=xend
z=zend
call vpumove(x,z)
x=xend-0.5*h
z=zend
call vpudraw(x,z)
call vpdash(.0,.0,.0,.0)
# finite-offset raypaths
call vpcolor(4)
# Shot to reflection point
x = xs
z = surface
x1 = xr
z1 = zr
call vpuarrow(x,z,x1,z1,0.02*zmax)
# reflection point to geophone
x = xr
z = zr
x1 = xg
z1 = surface
call vpuarrow(x,z,x1,z1,0.02*zmax)
# text
call vpcolor(5)
call vpfat(plotfat)
call vptjust(4,6)
call vputext(xr,zr+2.0*slop,8,0,'R')
if(option == 3)
call vputext(xrp+0.5*slop,zrp+slop,8,0,'R\\F15 \\v131 \\-')
x=0.5*(xrp+y)
z=0.5*(zrp+surface)
if(option == 3) call vputext(x-slop,z,8,0,'d')
# theta
call vputext(xend-5*slop, zend-1.*slop,8,0,'\\F9 q')
if( option==2 ) {
call vputext(xz,surface-slop,8,0,'Z')
call vputext(xr,surface-slop,8,0,'C')
}
# draw dashed lines
call vpdash(.1,.1,.1,.1)
call vpfat(plotfat)
# R' to M
call vpcolor(6)
x=xrp
z=zrp
call vpumove(x,z)
x=y
z=surface
if(option == 3) call vpudraw(x,z)
# R to Z
call vpdash(0.,0.,0.,0.)
x=xr
z=zr
call vpumove(x,z)
x=xz
z=surface
if(option == 2) call vpudraw(x,z)
# Right triangle RCZ
if(option == 2) {
call vpdash(.1,.1,.1,.1)
call vpcolor(6)
call vpumove(xr,zr)
call vpudraw(xr,surface)
}
1000 continue
# draw a frame
if(wantframe != 0)
{
call vpdash(0.,0.,0.,0.)
call vpcolor(7)
call vpfat(1)
call vpumove(xmin,zmin)
call vpudraw(xmin,zmax)
call vpudraw(xmax,zmax)
call vpudraw(xmax,zmin)
call vpudraw(xmin,zmin)
}
stop
end