inversion

su 的源代码库

#! /bin/sh
# Invert the data created by CSHOT

recordfile=record

echo "Run time information being stored in file:  $recordfile"
>$recordfile

################
cd inversion.dir	# where inversion files are
################
/bin/rm -f getcosines
echo "Please standby. The inversion will take about 5 min."
make

(
### Preprocessing
/bin/rm -f DATA1
cz1in <data1in

/bin/rm -f PTRACES
cz1

### Inversion
cxzcs

) > $recordfile 2>&1

sleep 20

/bin/mv -f $recordfile ..
echo "run time info now in file:  $recordfile"

### Strip the Fortran bytes
ftnstrip <reflect >reflect.strip;       mv -f reflect.strip reflect
ftnstrip <reflectcos >reflectcos.strip; mv -f reflectcos.strip reflectcos

### Plot the sections
xmovie <reflect n1=301 n2=50 d1=25 d2=80 f2=3000 perc=99 style=seismic \
	-geometry 560x430+50+570 \
	title=" Inversion (R)" \
	label1=" Depth(ft)"  label2=" Midpoint(ft)" &

xmovie <reflectcos n1=301 n2=50 d1=25 d2=80 f2=3000 perc=99 style=seismic \
	-geometry 560x430+650+570 \
	title=" Inversion (R * cos)" \
	label1=" Depth(ft)"  label2=" Midpoint(ft)" &


### Extract the cosines from known locations on the sections
subset <reflectcos n1=301 if1s=60 n1s=1 n2=50 if2s=0 n2s=50 >a1
subset <reflect    n1=301 if1s=60 n1s=1 n2=50 if2s=0 n2s=50 >b1

subset <reflectcos n1=301 if1s=100 n1s=1 n2=50 if2s=0 n2s=50 >a2
subset <reflect    n1=301 if1s=100 n1s=1 n2=50 if2s=0 n2s=50 >b2
 
farith op=div <a1 in2=b1>cos1
farith op=div <a2 in2=b2>cos2

# Run fortran program to compute the exact cosine of angle
#   produces files out1, out2 for the 2 reflectors
rm out1 out2
getcosines
ftnstrip <out1 >out1.strip; /bin/mv -f out1.strip out1
ftnstrip <out2 >out2.strip; /bin/mv -f out2.strip out2

# Plot computed against exact cosines
cat cos1 out1 |
xgraph nplot=2 n=50,50  d1=80,80 f1=3000 style=normal \
	-geometry 560x430+50+85 \
	x2beg=0. x2end=1.1  \
	linecolor=2,4 linewidth=3 \
	title="First interface (blue=exact)"  \
	label1="midpoint(ft)" label2="Cosine of angle" &

cat cos2 out2 |
xgraph nplot=2 n=50,50  d1=80,80 f1=3000 style=normal \
	-geometry 560x430+650+85 \
	x2beg=0. x2end=1.1  \
	linecolor=2,4 linewidth=3 \
	title="Second interface (blue=exact)" \
	label1="midpoint(ft)" label2="Cosine of angle" &
  
exit 0