program Qdome
  ! Qdome -- Make quarter dome model.
  !
  ! Qdome [n1=400 n2=100 n3=50 o1=0 d1=.004 ... ] > out.H
  !%
  use sep
  use random_mod
  use triangle

  implicit none
  real                             :: ranget,o1,d1,o2,d2,o3,d3
  real                             :: gaussvel, throw, rand, t1, t2
  real                             :: frac1, frac2, f1, zrad, xnew
  real                             :: pi, thetac, dr, x, y, z, rmax
  real                             :: xp, yp, rcap, rad1, rad2, dz
  real                             :: r, t, frac, zcap, ynew, interval
  integer                          :: large, n1, n2, n3, endtaper, stat
  integer                          :: i1,  i2, i3, n1a, n1b, it1, it2, i
  real, dimension (:),     pointer :: imp, refl, refl1, refl2
  real, dimension (:,:,:), pointer :: earth
real :: rat

  call sep_init ()
  call from_par ("ranget", ranget, 3.)
  call from_either ("n1", n1, 400)
  call from_either ("n2", n2, 100)
  call from_either ("n3", n3, 50 )
  call from_either ("large", large, int(n1*ranget))
  call to_history ("esize",4)

  call from_par ("d1", d1, .004) ; call from_par ("o1", o1, 0.   )  
  call from_par ("d2", d2, .01 ) ; call from_par ("o2", o2, -3*d2)
  call from_par ("d3", d3, .02 ) ; call from_par ("o3", o3, -3*d3)

   rat=d1/d2
  call from_param("ratio",rat,rat)
write(0,*)  "CHECK RAT",rat

  call from_par ("gaussvel", gaussvel,2.5)
  call from_par ("throw", throw, .01)
  call from_par ("endtaper", endtaper, 20)

  call to_history ("label1","sec")
  call to_history ("label2","W-E,km")
  call to_history ("label3","N-S,km")
  call sep_close

  allocate (imp(large), earth(n1,n2,n3), refl (n1))

  call random_init (19921995)

  imp(1) = 0.
  do i1= 2, large 
     call random (rand)
     if (rand > .3) then
        imp(i1) = imp(i1-1)
     else 
        call random (rand)
        imp(i1) = .1 * rand
     end if
  end do

  do i1= 2, large  
     imp(i1) = imp(i1) + imp(i1-1)
  end do

  earth = 1.0

  pi     = 4. * atan(1.)
  thetac = 45 * 2*pi/360.
  dr = abs(d2) + abs(d1)
  x = o2+d2*(n2-1)
  y = o3+d3*(n3-1)
  rmax = sqrt( x**2 + y**2 ) + dr
  xp = 0.6 * x
  yp = 0.75 * y
  rcap = 0.2 * rmax
  rad1 = (xp**2 + yp**2) - rcap
  rad2 = (xp**2 + yp**2)
  n1a = .15 * n1
  n1b = .8  * n1
  dz = 1./(n1b-n1a)

  do i3= 1, n3   
     y = o3+d3*(i3-1)
     do i2= 1, n2
        x = o2+d2*(i2-1)
        r = sqrt( x**2 + y**2)  + dr/2.
        earth (1 : n1a-1,i2,i3) = 1. + (imp(2:n1a) - imp(:n1a-1))
        ! shallow beds
        do i1=    n1b+1, n1      ! deep beds
           t = i1 +rat * i2 +.2 * i3
           it1  = t
           it2  = t+1
           frac = t - it1
           interval = &
           ((imp( it2) * (1.-frac) + imp( it2+1) * frac) - &
           (imp( it1) * (1.-frac) + imp( it1+1) * frac))
           earth( i1, i2, i3) =  1.0   + interval
        end do

        do i1= n1a, n1b 
           z = float( i1 - n1a) / (n1b - n1a) ! 0 to 1
           zcap = z * gaussvel /4.0
           z = dz + (1-2*dz) * z   ! dz to 1-dz

           t1 = (large/gaussvel) * (r/rmax) / sqrt( -alog( z-dz/2))
           t2 = (large/gaussvel) * (r/rmax) / sqrt( -alog( z+dz/2))

           t1 = (large/gaussvel) * (r/rmax) * sqrt( -alog( z+dz/2))
           t2 = (large/gaussvel) * (r/rmax) * sqrt( -alog( z-dz/2))

           if( t1 > t2) call erexit('misconception')

           xnew = cos(thetac) * (x-xp) + sin(thetac) * (y-yp)
           xnew = xnew * 0.75
           ynew = -sin(thetac) * (x-xp) + cos(thetac) * (y-yp)
           if ( rcap**2 > ( xnew**2 + ynew**2 + zcap**2 )  ) then
              if (r < rad2) then
                 zrad = (rad2 - r)/rcap 
              else     
                 zrad = 0.0  ! no faulting after center 
              end if
              t1  =  t1  + throw * n1 * zrad 
              t2  =  t2  + throw * n1 * zrad 
           end if
           it1  = t1
           it2  = t2
           if( it1 >=1 .and. it1+1 < large) then
              frac1 = t1 - it1
              if( it2 >=1 .and. it2+1 < large) then
                 frac2 = t2 - it2
                 interval = &
                 ((imp( it2) * (1.-frac2) + imp( it2+1) * frac2) - &
                 (imp( it1) * (1.-frac1) + imp( it1+1) * frac1))
                 earth( i1, i2, i3) =  1.0   + interval / (t2 - t1)
              end if
           end if
        end do
     end do
  end do


  !   Form reflectivity from impedance.
  refl (1) = 0.
  do i3 = 1, n3
     do i2 = 1, n2
        refl2 => earth (2:    ,i2,i3)
        refl1 => earth ( :n1-1,i2,i3)
        refl (2:) = (refl1 - refl2) / (refl1 + refl2)
        earth (:,i2,i3) = refl
     end do
  end do


  ! temporal smoothing
  do i = 2, 4
     call triangle_init (i, n1)
     do i3 = 1, n3
        do i2 = 1, n2
           stat = triangle_lop (.true.,.false.,refl,earth (:,i2,i3))
           earth (:,i2,i3) = refl
        end do
     end do
     call triangle_close ()
  end do

  refl (1) = 0.
  do i3 = 1, n3
     do i2 = 1, n2
        refl2 => earth (2:    ,i2,i3)
        refl1 => earth ( :n1-1,i2,i3)
        refl (2:) = (refl2 - refl1)
        earth (:,i2,i3) = refl
     end do
  end do

  f1 = real(endtaper)
  do i1 = 1, endtaper
     earth (     i1   ,:,:) = earth (     i1    ,:,:) * i1  / f1
     earth (n1 - i1 +1,:,:) = earth (n1 - i1 + 1,:,:) * i1  / f1
  end do


  call sep_write (earth)

  deallocate (earth, imp, refl)

  call exit (0)
end program Qdome

