c string.f
C- 01 NOV 1997 Seth Stein
C- Modified to be more compatible with GMT
C- 02 NOV 1997 JED
C- f77 -o STRING string.f -lm -lmath
c version of istring.f: find “snapshots” of displacement u(x)
c at various times t for waves propagating on
c inhomogeneous string. normal modes solution
c following exact solution in Geller & Stein 1978
c assumes junction at alngth/2. between string of
c density rho1, stiffness c1 and density rho2, stiffness c2
c outputs successive frames to file fort.N
c needs math library subroutine zreal, so compile with -lmath
character lab(70)
dimension para(22)
dimension omega(300),asym(300)
dimension aa(300),bb(300),qq(300)
dimension u(200),x(200),xs(2),ys(2)
dimension uu(200),xx(200)
common /consts/ rho1,rho2,c1,c2,alngth
external root
pi = 3.1415927
write (6,*) ‘ synthetic seismogram for string’
c length, number of length increments
read(5,*) alngth, nlen
write (6,*) ‘ string length ‘,alngth,’number of steps’, nlen
c space step
dx=alngth/float(nlen)
write (6,*) ‘space step ‘, dx
c time step (sec), number of time steps
read(5,*) dt, ntstep, tzero
write (6,*) ‘time step ‘,dt,’number of steps’, ntstep
write (6,*) ‘time start ‘,tzero
c source position
read(5,*) xsrc
write (6,*) ‘ source at ‘,xsrc
c number of modes
read(5,*) nmode
write (6,*) ‘ number of modes ‘,nmode
c source shape term
read(5,*) tau
write (6,*) ‘ source shape term’,tau
c string constants
read(5,*) c1, rho1
read(5,*) c2, rho2
c compute acoustic impedences
aip1=c1*rho1
aip2=c2*rho2
c compute reflection & transmission coefficients
refl12=(aip1-aip2) /(aip1+aip2)
trans12=2.*aip1 /(aip1+aip2)
refl21=(aip2-aip1) /(aip1+aip2)
trans21=2.*aip2 /(aip1+aip2)
c compute tensions
ten1=c1*c1*rho1
ten2=c2*c2*rho2
c list parameters
write (6,*) ‘lhs: c, rho,impedance, tension : ‘,c1,rho1,aip1,ten1
write (6,*) ‘rhs: c, rho,impedance, tension : ‘,c2,rho2,aip2,ten2
write (6,*) ‘reflection coefficient l-r’, refl12
write (6,*) ‘transmission coefficient l-r’, trans12
write (6,*) ‘reflection coefficient r-l’, refl21
write (6,*) ‘transmission coefficient r-l’, trans21
c parameters for root searching
eps=.1
eps2=.1
eta=0.2
nsig=5
itmax=6
pi=3.141596
c set up initial guesses for eigenfrequencies
caverg=2./(1./c1 + 1./c2)
do 10 i=1,nmode
asym(i)=i*pi*caverg/alngth
omega(i)=asym(i)
10 continue
c find eigenfrequencies
write(6,*) ‘calling zreal2’
call zreal2(root,eps,eps2,eta,nsig,nmode,omega,itmax,ier)
write(6,*) ‘roots found: ier=’,ier
write(6,*) ‘mode asymptote root’
c find constants for eigenfunctions
do 20 i=1,nmode
c write(6,*) i, asym(i), omega(i)
call eigcon (omega(i),aa(i),bb(i),qq(i))
20 continue
c
nlen1=nlen+1
c loop over times
do 50 j = 1, ntstep
t=dt*j + tzero
c initialize displacement
do 15 k = 1, nlen1
15 u(k) = 0.0
c
c outer loop over modes
do 40 i = 1, nmode
c mode frequency
wn = omega(i)
dmp = (tau*wn)**2
c time dependant term
scale = exp(-dmp)*cos(wn*t)
c source position term
phis= phi (xsrc,i,wn,aa(i),bb(i),qq(i))
c compute displacement at each point
do 40 k = 1, nlen1
x(k)=(k-1)*dx
phix= phi (x(k),i,wn,aa(i),bb(i),qq(i))
40 u(k) = u(k)-phis*phix*scale
c find value to normalize u(x) values from first time step
if (j.eq.1) then
call maxmin(nlen,u,amax,amin,imax,imin)
anorm= max(abs(amax),abs(amin))*2.
write (6,*) ‘wave normalization’,anorm
endif
c output this frame to fort.ifile
C- First, output the left-hand side data
ifile=9+j
write(6,*) ‘frame number ‘,j,ifile
write(ifile,*) j,t
do 42 k = 1, (nlen/2)
write(ifile,*) x(k),u(k)/anorm,0.5
42 continue
close(ifile)
C- Now output the right-hand side data
ifile=109+j
write(6,*) ‘frame number ‘,j,ifile
write(ifile,*) j,t
do 43 k = (nlen/2), (nlen1-1)
write(ifile,*) x(k),u(k)/anorm,1.5
43 continue
close(ifile)
50 continue
stop
end
real function root (freq)
c secular equation for eigenfrequencies assuming junction at alngth/2.
common /consts/ rho1,rho2,c1,c2,alngth
oml=freq*alngth
oml1=oml/(2.*c1)
oml2=oml/(2.*c2)
ratio=rho1*c1/(rho2*c2)
root=sin(oml1)*cos(oml2) + ratio*cos(oml1)*sin(oml2)
return
end
subroutine eigcon (freq,a,b,q)
c given eigenfrequency compute constants for eigenfunctions
c assuming junction at alngth/2.
common /consts/ rho1,rho2,c1,c2,alngth
oml=freq*alngth/2.
oml1=oml/c1
oml2=oml/c2
ratio=rho1*c1/(rho2*c2)
a=sin(oml1)*sin(oml2) + ratio*cos(oml1)*cos(oml2)
b=sin(oml1)*cos(oml2) – ratio*cos(oml1)*sin(oml2)
c normalization
omll=freq*alngth
omll1=omll/c1
omll2=omll/c2
x1=alngth/4. – c1*sin(omll1)/(4.*freq)
x2=alngth *(a**2.0 + b**2.0)/4.
x3=(c2 / (4.*freq))*(b**2.0 – a**2.0)
x4=sin(2.*omll2)- sin(omll2)
x5=c2*a*b / (2.*freq)
x6=cos(2.*omll2)- cos(omll2)
q= rho1*x1 + rho2 * (x2 + x3 * (x4 – x5*x6))
c check by computing left and right displacements
ald=sin(oml1)
ard=a*sin(oml2) + b*cos(oml2)
diff=ald-ard
write(6,*) ‘a,b,left,right,difference,normalization’
write(6,*) a,b,ald,ard,diff,q
return
end
real function phi (x,i,freq,a,b,q)
c evaluate ith eigenfunction with eigenfrequency freq at point x
common /consts/ rho1,rho2,c1,c2,alngth
sq=sqrt(q)
if (x.le.alngth/2.) then
phi=sin(freq*x/c1)/sq
else
phi=(a*sin(freq*x/c2) +b*cos(freq*x/c2))/sq
endif
return
end
subroutine maxmin(ndat,dat,amax,amin,imax,imin)
c find max and min and their indicies
dimension dat(1)
amax=dat(1)
imax=1
amin=dat(1)
imin=1
do 10 i=2,ndat
if(dat(i).le.amax) go to 15
amax=dat(i)
imax=i
go to 10
15 if(dat(i).ge.amin) go to 10
amin=dat(i)
imin=i
10 continue
return
end