!
!Mon Jan 13 08:01:39 PST 2003
!
module tomo_op_mod
  use sep3d_struct_mod
  use calc_sigma_mod
  use basic_util_mod
  use sep
  implicit none
  logical, private,         save :: tau,first_sreed
  integer, private, pointer,save :: nlens(:),nsect(:),nsize(:),lens(:)
  real,    private, pointer,save :: scale1(:),scale2(:),scale3(:)
  real,    pointer, private,save :: tempm(:),tempm2(:)
  integer, private, pointer,save :: iloc(:)
  integer, private,         save :: in_mem,ndims,verb

  contains
  logical function load_tomo_op(slow_s,err0)
    type(sep3d) :: slow_s,err_s
    real, pointer :: err0(:)
    integer,external :: sreed

     load_tomo_op=.false.
     call from_param("tau",tau,.false.)
     call from_param("verb",verb,0) 
     ndims=sep3d_ndims(slow_s)
     first_sreed=.true.

     if(.not. calc_load_op()) then
       write(0,*) "trouble in reading operator"
       return
     end if

     call init_sep3d("err",err_s,"INPUT")
     allocate(err0(product(err_s%n(1:sep3d_ndims(err_s)))))
     if(4*size(err0) /= sreed("err",err0,4*size(err0))) then
       write(0,*) "trouble reading err0"
       return
     end if

     
     if(tau) allocate(tempm(product(slow_s%n(1:ndims))))
     if(tau .and. ndims==3) allocate(tempm(product(slow_s%n(1:ndims))))
     in_mem=0
  

     load_tomo_op=.true.
  end function

  logical function calc_load_op()
    type(sep3d) :: nlen_s
    integer     :: mem,mem_left,isect,mem_sect,i,fact
    integer, pointer :: temp(:)
    integer,external :: sreed
    
   
     calc_load_op=.false.

     call from_param("maxsize",mem,900)
     mem=mem*250000
     call init_sep3d("nlen",nlen_s,"INPUT")
     allocate(nlens(product(nlen_s%n(1:sep3d_ndims(nlen_s)))))
     if(4*size(nlens) /= sreed("nlen",nlens,4*size(nlens))) then
       write(0,*) "trouble reading nlens"
       return
     end if


     fact=2
     if(tau) then
       fact=3
       if(ndims==3) fact=4
     end if
     i=1
     allocate(nsize(1))
     allocate(nsect(1))
     do while  (i <= size(nlens)) 
       nsize(size(nsize))=0
       nsect(size(nsect))=0
       do while  (i <= size(nlens) .and. nsize(size(nsize))*fact<mem ) 
          nsize(size(nsize))=nsize(size(nsize))+nlens(i)
          nsect(size(nsect))=nsect(size(nsect))+1
          i=i+1
       end do
       if(i<size(nlens)) then
         allocate(temp(size(nsize)))
         temp=nsize
         deallocate(nsize)
         allocate(nsize(size(temp)+1))
         nsize(1:size(temp))=temp
         deallocate(temp)
         allocate(temp(size(nsect)))
         temp=nsect
         deallocate(nsect)
         allocate(nsect(size(temp)+1))
         nsect(1:size(temp))=temp
         deallocate(temp)
       end if 
     end do
     if(verb>0) write(0,*) "nsections =",size(nsect)

     allocate(scale1(maxval(nsize)))
     allocate(lens(maxval(nsect)))
     allocate(iloc(maxval(nsize)))
     call auxinout("scale")
     call auxinout("iloc")
     if(tau) then
       call auxinout("scale_sigx")
       allocate(scale2(maxval(nsize)))
       if(ndims==3) then
         call auxinout("scale_sigy")
         allocate(scale3(maxval(nsize)))
       end if
     end if
     calc_load_op=.true.
  end function

  subroutine load_section(isect)
    integer, external :: sreed
    integer ::  i,ibeg,iend,isect
    if(isect==0) then
      call  sseek("iloc",0,0)
      call  sseek("scale",0,0)
      if(tau) call  sseek("scale_sigx",0,0)
      if(tau .and. ndims==3) call  sseek("scale_sigy",0,0)
    end if
    if(isect==1 .and.  .not. first_sreed ) then
      call sseek("iloc",0,0)
      call sseek("scale",0,0)
      if(tau) then 
        call sseek("scale_sigx",0,0)
        if(ndims==3) then
          call sseek("scale_sigy",0,0)
        end if
      end if
    end if
    first_sreed=.false.
    if(4*nsize(isect)/=sreed("iloc",iloc,4*nsize(isect))) &
      call seperr("trouble loading section(iloc)")
    if(4*nsize(isect)/=sreed("scale",scale1,4*nsize(isect))) &
      call seperr("trouble loading section(scale)")
    call sect_lim(isect,ibeg,iend)
    lens=nlens(ibeg:iend)
    if(tau) then
      if(4*nsize(isect)/=sreed("scale_sigx",scale2,4*nsize(isect))) &
        call seperr("trouble loading section(scale_sigx)")
      if(ndims==3) then
        if(4*nsize(isect)/=sreed("scale_sigy",scale2,4*nsize(isect))) &
          call seperr("trouble loading section(scale_sigy)")
      end if 
    end if
    in_mem=isect
  end subroutine

  subroutine sect_lim(isect,ibeg,iend)
    integer :: isect,i,ibeg,iend
    ibeg=1
    ibeg=1
    do i=1,isect-1
      ibeg=ibeg+nsect(i)
    end do
    iend=ibeg+nsect(isect)-1
  end subroutine
  integer function tomo_op(adj,add,model,data)
    logical,intent(in) :: adj,add
    real :: model(:),data(:)
    integer :: isect,ierr
    tomo_op=0 

!write(0,*) "IF TOMO",adj
    call adjnull(adj,add,model,data)

    if(.not. adj .and. tau) then
      ierr=calc_delta_sigmax(.false.,.false.,model,tempm)
      if(ndims==3) ierr=calc_delta_sigmay(.false.,.false.,model,tempm2)
    end if

   ! if(adj)  call srite("sss.H",model,size(model)*4)

    call adjnull(adj,add,model,data)
    do isect=1,size(nsect)
      if(isect/=in_mem) call load_section(isect)
      ierr=apply_section(adj,add,model,data,scale1)
      if(.not. adj  .and. tau) then
        ierr=apply_section(adj,.true.,tempm,data,scale2)
        if(ndims==3) ierr=apply_section(adj,.true.,tempm2,data,scale3)
      else if(tau) then
        ierr=apply_section(adj,add,tempm,data,scale2)
        if(ndims==3) ierr=apply_section(adj,add,tempm2,data,scale3)
      end if
    end do
          
    if(adj .and. tau) then
      ierr=calc_delta_sigmax(.true.,.true.,model,tempm)
      if(ndims==3) ierr=calc_delta_sigmay(.true.,.true.,model,tempm2)
    end if
    if(adj) call srite("sss.H",model,size(model)*4)
  end function

  integer function apply_section(adj,add,model,data,scale)
    logical, intent(in) :: adj,add
    real :: model(:),data(:),scale(:)
    real :: tots,totm,totd
    integer :: i,ibeg,iend,ibig,ismall,itrace,iray

    call sect_lim(in_mem,ibeg,iend)
    i=0
    ibig=0
    tots=0
    totm=0
    totd=0
    do iray=ibeg,iend
      ibig=ibig+1
!write(0,*) "BABB",ibeg,iend,lens(ibig),i,iray
      do itrace=1,lens(ibig)
        i=i+1
        if(.not. adj) then 
           data(ibig)=data(ibig)+model(iloc(i))*scale(i)
        else
           model(iloc(i))=model(iloc(i))+data(ibig)*scale(i)
           tots= tots+scale(i)
        end if
      end do
!write(0,*) "ABB",iray,adj,data(ibig)
    end do
    apply_section=0  
  end function

end module
