
!---------------------------------------------------------------------------
! Build a synthetic velocity model.
!---------------------------------------------------------------------------
program Velmod
  use sep

  implicit none
  integer 				:: n1, n2, ns, is
  real, dimension (:,:), allocatable	:: slow
  real 					:: o1,d1, s0,ds, x0, dx
  real					:: xmin, xmax, vmin, vmax, tmax
  real 					:: v1,v2,v3, s1,s2,s3

  call sep_init( )
  call from_param ("nt", n1, 100)
  call from_param ("nx", n2, 30)
  call from_param ("ns", ns, 25)

  call to_history ("n1",n1) 
  call to_history ("n2",ns)
  call to_history ("esize",4)

  call from_param ("tmax", tmax, 4.)
  call from_param ("xmax", xmax, 4.5 ) ; call from_param ("xmin", xmin, -1.)
  call from_param ("vmax", vmax, 100.) ; call from_param ("vmin", vmin, 1.2)

  call from_param ("o1", o1, 0.)
  call from_param ("x0", x0, xmin)
  call from_param ("s0", s0, 1./vmax)

  call from_param ("d1", d1, tmax/n1)
  call from_param ("dx", dx, (xmax-xmin)/(n2-1.))
  call from_param ("ds", ds, (1./vmin-1./vmax)/(ns-1.))

  call to_history ("o2",s0)
  call to_history ("d2",ds)
  call to_history ("label2","slowness(km/sec)")
  call sep_close ()

  allocate (slow(n1,ns)) ; slow = 0.

  v1=0.30;	v2=2.2;		v3=1.5;
  s1=1/v1;	s2=1/v2;	s3=1/v3;
  is = 1.5 + (s2-s0)/ds;	
  s2 = s0 + ds*(is-1)	
  is= max(1., min(1.*ns, 1.5+(s2-s0)/ds)); slow(int(.25 * tmax/d1), is) = 1.
  is= max(1., min(1.*ns, 1.5+(s2-s0)/ds)); slow(int(.45 * tmax/d1), is) = 1.
  is= max(1., min(1.*ns, 1.5+(s2-s0)/ds)); slow(int(.75 * tmax/d1), is) = 1.
  is= max(1., min(1.*ns, 1.5+(s2-s0)/ds)); slow(int(.94 * tmax/d1), is) = 1.
  is= max(1., min(1.*ns, 1.5+(s3-s0)/ds)); slow(int(.15 * tmax/d1), is) = 2.

  call sep_write (slow)

  deallocate (slow)
  call exit (0)

end program Velmod


