  program refcoef

  use sep
  use tcai1_mod

  implicit none

  real, dimension(:,:),allocatable	:: dado, sism
  real,dimension(:),allocatable		:: w, rc
  integer                      		:: n1, n2, n3, i2, i3, it, stat
  real                        		:: o1, o2, o3, d1, d2, d3 
  real                        		:: gama, offset
  real                        		:: zr, v1, v2
  real, parameter             		:: PI = 3.14159265359
  
  allocate(w(31))
  w(1)=105;w(2)=544;w(3)=2156;w(4)=6888;w(5)=18292;w(6)=41120;w(7)=79002;w(8)=129948;w(9)=181506;w(10)=209664;w(11)=185640
  w(12)=90624;w(13)=-68184;w(14)=-251328;w(15)=-398517;w(16)=-454950;w(17)=-398517;w(18)=-251328;w(19)=-68184;w(20)=90624;w(21)=185640
  w(22)=209664;w(23)=181506;w(24)=129948;w(25)=79002;w(26)=41120;w(27)=18292;w(28)=6888;w(29)=2156;w(30)=544;w(31)=105

  call sep_init()
  call sep_begin_prog()

  call initpar()

  if(exist_file("in"))  then
	call from_history("n1", n1); call from_history("d1", d1); call from_history("o1", o1 )
	call from_history("n2", n2); call from_history("d2", d2); call from_history("o2", o2 )  
	call from_history("n3", n3); call from_history("d3", d3); call from_history("o3", o3 )
  else
	call from_param("n1", n1, 501);	call from_param("d1",d1,0.004);	call from_param("o1",  o1,  0.)
	call from_param("n2", n2, 401);	call from_param("d2", d2, 20.);	call from_param("o2",o2,-4000.)  
	call from_param("n3", n3, 500);	call from_param("d3", d3, 20.);	call from_param("o3",  o3,  0.)
  end if

  call from_param("zr", zr, 1000.)
  call from_param("v1", v1, 4000.)
  call from_param("v2", v2, 3000.)

  call to_history("n1", n1); call to_history("d1", d1); call to_history("o1", o1 )
  call to_history("n2", n2); call to_history("d2", d2); call to_history("o2", o2 )  
  call to_history("n3", n3); call to_history("d3", d3); call to_history("o3", o3 )

  call tcai1_init (size(w),w)

  allocate(dado(n1,n2),sism(size(w)+n1,n2),rc(n1))
  sism=0.;dado=0.
  do i2=1,n2
	rc=0.
	offset = (i2-1)*d2 + o2
	gama=atan((offset)/(2*zr))
	it=int(sqrt((2*zr/v1)**2+(offset/v1)**2)/d1+1.5)
        if ((it <= n1).and.(it >= 1)) then
	   rc(it)=(v2*cos(gama)-sqrt(v1**2-(v2*sin(gama))**2))/(v2*cos(gama)+sqrt(v1**2-(v2*sin(gama))**2))
	   rc(it)=rc(it)/sqrt(2*zr/cos(gama))
	end if
	stat = tcai1_lop ( .false., .true., rc, sism(:,i2))
  end do

  dado=sism(int(size(w)/2+2.5):int(size(w)/2+2.5)+n1,:)

  do i3=1,n3
	call srite("out", dado, 4*n1*n2)
  end do
  deallocate(sism,dado,rc,w)
  call sep_end_prog()
  end program
