program OW_reflect

use sep

implicit none
real			:: o1,o2,d1,d2,r3,f,fz,gz
integer			:: i1,i2,n1,n2,i,im
real,allocatable		:: refl(:,:)

call sep_init()
call sep_begin_prog()

call from_param("n1",n1);call from_param("d1",d1);call from_param("o1",o1)
call from_param("n2",n2);call from_param("d2",d2);call from_param("o2",o2)

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("eout",4);call to_history("label3","")

allocate(refl(n1,n2))

do i2=1,n2
   r3=1000-((((i2-1)*d2+o2)/4000)**4*1000)
   f=r3/d1
   i=f
   im=1+i
   if ( 1<im .and. im < n1) then
      fz=f-i
      gz= 1.-fz
      refl(im,i2)=refl(im,i2)+gz
      refl(im+1,i2)=refl(im+1,i2)+fz
   end if
   r3=1400-((((i2-1)*d2+o2)/5000)**4*1000)
   f=r3/d1
   i=f
   im=1+i
   if ( 1<im .and. im < n1) then
      fz=f-i
      gz= 1.-fz
      refl(im,i2)=refl(im,i2)-gz
      refl(im+1,i2)=refl(im+1,i2)-fz
   end if
!
!      i1=int((1000.-(r3/4000.)**4*1000.)/d1)+1.5
!   if (i1 > 1 .and. i1 <= n1) refl(i1,i2)=1
!   i1=int((1400.-(r3/5000.)**4*1000.)/d1)+1.5
!   if (i1 > 1 .and. i1 <= n1) refl(i1,i2)=-1
end do

call srite("out",refl,4*n1*n2)
deallocate(refl)
call sep_end_prog()
end program
