C*********************************************************************
C
C Function: Three dimension seismic ray tracing
C
program main
C
USE MSFLIB
C
real slow[allocatable](:,:,
real v[allocatable](:,:,
real time[allocatable](:,:,
real ft[allocatable](:)
integer*2 nflag[allocatable](:,:,:)
integer*2 nfrom[allocatable](:,:,:,:)
real tl[allocatable](:)
integer*2 lx[allocatable](:)
integer*2 ly[allocatable](:)
integer*2 lz[allocatable](:)
integer*4 nsn[allocatable](:,:,:)
real rmodel[allocatable](:,:,:)
integer*2 modelf[allocatable](:,:,:)
integer*2 mray(0:2000,1:3)
integer*2 rcx(0:2100),rcy(0:2100),rcz(0:2100)
character file1*64,file2*64,file3*64,file4*64,file5*64
external mdflag
c
CHARACTER MESSAGE0*40,MESSAGE1*40,MESSAGE2*40,MESSAGE3*50
CHARACTER START1*20,START2*20,START3*20
INTEGER*2 DUMMY,XA1,XA2,YA1,YA2,DYA,XA,YA
INTEGER*4 ZPS,JSXS,JSDS,JSDZS,WGDSX,WGDSY,WGDSZ
REAL DJJ,ZXPJJ,ZDPJJ,JLCD,CYJG,WGDXX,WGDXY,WGDXZ,MXQDSD
LOGICAL statusmode
RECORD /xycoord/ xy
TYPE (windowconfig) wc
c
open(10,file=\'trace3d1.dat\',status=\'old\')
read(10,10) file1,file2,file3,file4,file5
read(10,*) isource,iorder
10 format(A/A/A/A/A)
close(10)
CALL CANSHU(ZPS,JSXS,JSDS,JSDZS,DJJ,ZXPJJ,ZDPJJ,JLCD,CYJG,
* WGDSX,WGDSY,WGDSZ,WGDXX,WGDXY,WGDXZ,MXQDSD)
c
iix=WGDSX
iiy=WGDSY
iiz=WGDSZ
ix=iix-1
iy=iiy-1
iz=iiz-1
dx=WGDXX
dy=WGDXY
dz=WGDXZ
mioffset=INT(ZXPJJ/DJJ)
maoffset=INT(ZDPJJ/DJJ)
c
allocate(slow(0:iix,0:iiy,0:iiz))
allocate(v(0:iix,0:iiy,0:iiz))
allocate(time(0:iix,0:iiy,0:iiz))
allocate(ft(JSDS))
allocate(nflag(0:iix,0:iiy,0:iiz))
allocate(nfrom(0:iix,0:iiy,0:iiz,1:3))
allocate(tl(0:iix*iiy*iiz))
allocate(lx(0:iix*iiy*iiz))
allocate(ly(0:iix*iiy*iiz))
allocate(lz(0:iix*iiy*iiz))
allocate(nsn(0:iix,0:iiy,0:iiz))
allocate(rmodel(-iix:iix,-iiy:iiy,-iiz:iiz))
allocate(modelf(-iix:iix,-iiy:iiy,-iiz:iiz))
c
ip=4
mr1=ip
mr2=ip
mr3=ip
c
C
C Set the vedio mode to a graphics mode
C
wc.numxpixels =640
wc.numypixels =480
wc.numtextcols = -1
wc.numtextrows = -1
wc.numcolors = -1
wc.title=\"Three Dimension Seimic Ray Tracing\" c
statusmode=setwindowconfig(wc)
if(.not. statusmode) THEN
statusmode=setwindowconfig(wc)
endif
dummy=INITIALIZEFONTS()
C
MESSAGE0=\'Three Dimension\'
MESSAGE1=\'Seismic Ray Tracing Processing\'
MESSAGE2=\'Processing Finished\'
MESSAGE3=\'Close this window to return the main menu\'
START1=\'0%\'
START2=\'rate of progress\'
START3=\'100%\'
DYA=30
C
C Write this processing name on screen
C
DUMMY=SETCOLOR(15)
DUMMY=RECTANGLE($GFILLINTERIOR,0,0,639,479)
DUMMY=SETCOLOR(4)
XA1=100
YA1=50
DUMMY=SETFONT(\'T\'\'ARIAL\'\'H20W12\')
CALL MOVETO(XA1,YA1,XY)
CALL OUTGTEXT(MESSAGE0)
YA1=80
DUMMY=SETFONT(\'T\'\'ARIAL\'\'H20W12\')
CALL MOVETO(XA1,YA1,XY)
CALL OUTGTEXT(MESSAGE1)
XA1=100
YA1=130
DUMMY=SETFONT(\'T\'\'ARIAL\'\'H16W10\')
CALL MOVETO(XA1,YA1,XY)
CALL OUTGTEXT(START1)
XA1=220
YA1=130
CALL MOVETO(XA1,YA1,XY)
CALL OUTGTEXT(START2)
XA1=460
YA1=130
CALL MOVETO(XA1,YA1,XY)
CALL OUTGTEXT(START3)
DUMMY=SETCOLOR(14)
XA1=99
YA1=149
XA2=501
YA2=YA1+DYA+2
CALL MOVETO(XA1,YA1,XY)
DUMMY=LINETO(XA1,YA2)
DUMMY=LINETO(XA2,YA2)
DUMMY=LINETO(XA2,YA1)
DUMMY=LINETO(XA1,YA1)
DUMMY=SETCOLOR(7)
XA1=XA1+1
YA1=YA1+1
XA2=XA2+1
YA2=YA2+1
CALL MOVETO(XA2,YA1,XY)
DUMMY=LINETO(XA2,YA2)
DUMMY=LINETO(XA1,YA2)
XA1=XA1+1
YA1=YA1+1
XA2=XA2+1
YA2=YA2+1
CALL MOVETO(XA2,YA1,XY)
DUMMY=LINETO(XA2,YA2)
DUMMY=LINETO(XA1,YA2)
XA1=XA1+1
YA1=YA1+1
XA2=XA2+1
YA2=YA2+1
CALL MOVETO(XA2,YA1,XY)
DUMMY=LINETO(XA2,YA2)
DUMMY=LINETO(XA1,YA2)
C
lbyte=(iz+1)*4
open(11,file=file1,status=\'old\',form=\'unformatted\',
& access=\'direct\',recl=lbyte)
c
c input slowness file
c
do 100 i=0,ix
do 110 j=0,iy
irec=1
read(11,rec=irec) (v(i,j,k),k=0,iz)
irec=irec+1
110 continue
XA1=100
YA1=150
XA=XA1+(i*10/ix+0.5)
YA=YA1+DYA
DUMMY=SETCOLOR(1)
DUMMY=RECTANGLE($GFILLINTERIOR,XA1,YA1,XA,YA)
100 continue
close(11)
c
c initialize
c
do 120 i=0,ix
do 130 j=0,iy
do 140 k=0,iz
slow(i,j,k)=1/v(i,j,k)
time(i,j,k)=10000.
140 continue
130 continue
120 continue
c
c prepare
c
do 150 i=-mr1,mr1
do 160 j=-mr2,mr2
do 170 k=-mr3,mr3
if((i.eq.0).and.(j.eq.0).and.(k.eq.0)) then
rmodel(i,j,k)=0.
modelf(i,j,k)=0
goto 170
endif
rmodel(i,j,k)=sqrt(i*i*dx*dx+j*j*dy*dy+k*k*dz*dz)
modelf(i,j,k)=mdflag(i,j,k)
170 continue
160 continue
150 continue
c
open(12,file=file2,status=\'old\')
open(14,file=file4,status=\'unknown\')
open(15,file=file5,status=\'unknown\',form=\'binary\')
c
c input source interface file
c start ray tracing
c
if(isource.ne.1) then
do 205 mmm=1,isource-1
read(12,*)msh1,msh2,nsy0,nsx0,nsz0
205 continue
else
endif
c
mmmm=0
do 200 mm=isource,ZPS
mmmm=mmmm+1
read(12,*)msh1,msh2,nsy0,nsx0,nsz0
my1=nsy0
my2=nsy0+maoffset
ns=0
nb=1
c
do 210 i=0,ix
do 220 j=my1,my2
do 230 k=0,iz
time(i,j,k)=100000.
nfrom(i,j,k,1)=0
nfrom(i,j,k,2)=0
nfrom(i,j,k,3)=0
nsn(i,j,k)=0
nflag(i,j,k)=0
230 continue
220 continue
210 continue
c
nsx=nsx0
nsy=nsy0
nsz=nsz0
time(nsx,nsy,nsz)=0.
nfrom(nsx,nsy,nsz,1)=nsx
nfrom(nsx,nsy,nsz,2)=nsy
nfrom(nsx,nsy,nsz,3)=nsz
nflag(nsx,nsy,nsz)=1
111 call trace(slow,time,nfrom,nflag,ix,my1,my2,iz,
$ rmodel,modelf,mr1,mr2,mr3,
$ nsx,nsy,nsz,iix,iiy,iiz,ns,tl,lx,ly,lz,nsn,nb)
if(nb.gt.WGDSX*maoffset*WGDSZ)goto 222
nsx=lx(nb)
nsy=ly(nb)
nsz=lz(nb)
nflag(nsx,nsy,nsz)=1
nb=nb+1
goto 111
222 open(13,file=file3,status=\'old\')
c
do 240 nn=1,JSXS
write(14,45) mm,nn
45 format(2I7)
jj=1
do 250 i=1,JSDZS
read(13,*)dd,rcy(i),rcx(i),rcz(i)
if(rcy(i).ge.nsy0+mioffset.and.rcy(i).le.nsy0+maoffset) then
k=0
mray(k,1)=rcx(i)
mray(k,2)=rcy(i)
mray(k,3)=rcz(i)
333 kx=mray(k,1)
ky=mray(k,2)
kz=mray(k,3)
k=k+1
mray(k,1)=nfrom(kx,ky,kz,1)
mray(k,2)=nfrom(kx,ky,kz,2)
mray(k,3)=nfrom(kx,ky,kz,3)
if(mray(k,1).ne.mray(k-1,1).or.mray(k,2).ne.mray(k-1,2)
$ .or.mray(k,3).ne.mray(k-1,3))then
goto 333
else
goto 444
endif
444 write(15)k
do 260 j=0,k-1
write(15)mray(j,2),mray(j,1),mray(j,3)
260 continue
ft(jj)=time(rcx(i),rcy(i),rcz(i))
jj=jj+1
endif
250 continue
if(iorder.eq.0) then
write(14,55)(ft(i),i=1,JSDS)
else
write(14,55)(ft(i),i=JSDS,1,-1)
endif
55 format(10f7.3)
240 continue
close(13)
C
XA1=110
YA1=150
