!!
!! The routine calculates a velocity model with sinusoidal
!! structure with varying period
!!          
!!
!!

  program velsin

  use sep

  implicit none
  
  real		  		    		:: pi,x,vi
  real						:: den,z0,z1,z2,v0,v1,v2,g0,alpha
  integer                         		:: i1,i2,ia,ib,ic
  integer					:: n1,n2,n3
  real                            		:: o1,o2,o3,d1,d2,d3
  real,dimension(:),allocatable 		:: vel

  call sep_init()
  call sep_begin_prog()

  if(exist_file("indata")) then
  !! get info from history file
    call from_aux("indata","n1",n1);	call from_aux("indata","n2",n2);	call from_aux("indata","n3",n3,1)
    call from_aux("indata","d1",d1);	call from_aux("indata","d2",d2);	call from_aux("indata","d3",d3,1.)
    call from_aux("indata","o1",o1);	call from_aux("indata","o2",o2);	call from_aux("indata","o3",o3,1.)
  else
  !! get info from parameters
    call from_param("n1",n1,500);   call from_param("n2",n2,500);   call from_param("n3",n3,1)
    call from_param("d1",d1,10.);   call from_param("d2",d2,10.);   call from_param("d3",d3,1.)
    call from_param("o1",o1,0.);   call from_param("o2",o2,0.);   call from_param("o3",o3,1.)
  end if

  !! get info from parameter
  call from_param("den",den,50000.);	call from_param("depth",z0,500.)
  call from_param("v0",v0,2000.);	call from_param("g0",g0,0.5)
  call from_param("v1",v1,4000.);	call from_param("v2",v2,3500.)
  call from_param("z1",z1,2000.);	call from_param("z2",z2,3000.)
  call from_param("alpha",alpha,0.1)


  !! write info out to history file
  call to_history("n1",n1);call to_history("n2",n2);call to_history("n3",1)
  call to_history("d1",d1);call to_history("d2",d2);call to_history("d3",1.)
  call to_history("o1",o1);call to_history("o2",o2);call to_history("o3",1.)
  call to_history ("label1","z");call to_history ("label2","cmpx");call to_history ("label3"," ")

!! allocate working space
  allocate(vel(n1))

pi = 4*atan(1.)
ia = z0/d1 + 1.5
ib = z1/d1 + 1.5
ic = z2/d1 + 1.5
if (ia .gt. n1) call seperr("Reflector 1 is deeper than the maximum depth")
if (ib .gt. n1) call seperr("Reflector 2 is deeper than the maximum depth")
if (ic .gt. n1) call seperr("Reflector 3 is deeper than the maximum depth")
if (ia .gt. ib) call seperr("Reflector 1 is deeper than Reflector 2")
if (ia .gt. ic) call seperr("Reflector 1 is deeper than Reflector 3")
if (ib .gt. ic) call seperr("Reflector 2 is deeper than Reflector 3")
do i2=1,n2
  x=(i2-1)*d2+o2
  vel=0.;vel(1:ia)=v0;vel(ib+1:ic)=v1;vel(ic+1:n1)=v2
  if (mod(i2,10).eq.0) write(0,*) "CMP  ", i2, " of ", n2

  do i1=ia+1,ib
	vi=v0*(1+alpha*sin((x*x*pi/180)/den))
	vel(i1)=vi+g0*((i1-1)*d1+o1)
  end do
  call sep_write(vel)
end do
!! clean up
deallocate(vel)
call sep_end_prog()
end program velsin
