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