program geodip
!
!infile=grd.H outfile=grdflt.H angles=xhh.H steer_lags=15 niter=20 epspoly=20. tri_width=6350
!infile=grd.H outfile=grdflt.H angles=xhh.H steer_lags=40 15 niter=20 epspoly=40. tri_width=130000
!infile=grd.H outfile=grdflt.H angles=xhh.H steer_lags=30 niter=12 epspoly=80. tri_width=130000
!infile=grd.H outfile=grdflt.H angles=xhh.H mask=1.H steer_lags=40 20 niter=12 epspoly=160. tri_width=1000000
!
  use sep
  use inv_tomo_mod
  use filters_mod

  implicit none
  type (tomo_param_type)                        :: param
  type (model_param_type)                       :: model_param
  real, allocatable, dimension (:,:,:)          :: slow0,slow,tempa
  character(len=128)                            :: model_tag
  integer                                       :: ierr,nelem_ray,nelem_model,verb,frame
  integer                                       :: n_h,n_mid,auxout,i1,i2,tempi,temp2i,temp3i,type,test_oper

!FILTER PARAMETERS
  type(filter)                                  :: filt
  integer                                       :: lap_lags,steer_lags
  real                                          :: tri_width,lap_factor,steer_factor,eps,epspoly
  real,allocatable                              :: angles(:),tangle(:,:)

!MASK
  real,allocatable                              :: mask(:),mask_in(:,:)

 call sep_init()
 call sep_begin_prog()

 call from_param('new',param%new,0)
 call from_param('verb',verb,1)
 call from_param('frame',frame,50)
 call from_param("type",type,2)
 call from_param('verb_solver',param%verb_solver,verb)
 param%approx=0
 param%test_oper=0
 call from_param('dot_test',param%dot_test,0)
 call from_param('test_oper',param%test_oper,0)
 call from_param('niter',param%niter,10)
 call from_param('epspoly',epspoly,1.)
 call from_param('ndim_mod',model_param%ndim_mod,1)
 call from_aux("infile","n2",model_param%nx_mod)
 call from_aux("infile","n1",model_param%nz_mod)
 call from_aux("infile","n3",model_param%ny_mod ,1)
 call from_aux("infile","o2",model_param%x0_mod)
 call from_aux("infile","o1",model_param%z0_mod)
 call from_aux("infile","d2",model_param%dx_mod)
 call from_aux("infile","d1",model_param%dz_mod)

 call to_history("n2",model_param%nx_mod,"outfile")
 call to_history("n1",model_param%nz_mod,"outfile")
 call to_history("n3",model_param%ny_mod,"outfile")
 call to_history("o2",model_param%x0_mod,"outfile")
 call to_history("o1",model_param%z0_mod,"outfile")
 call to_history("d2",model_param%dx_mod,"outfile")
 call to_history("d1",model_param%dz_mod,"outfile")

 write(0,*) "Done with initial parameter reading"
!ALLOCATE MODEL
  allocate (slow0(model_param%nz_mod+frame*2,model_param%nx_mod,model_param%ny_mod))
  allocate (slow(model_param%nz_mod+frame*2,model_param%nx_mod,model_param%ny_mod));slow=0.
  call sreed('infile',slow0(frame+1:frame+model_param%nz_mod,:,:),model_param%ny_mod*&
            model_param%nz_mod*model_param%nx_mod*4)
  do i1=1,frame
    slow0(i1,:,:)=slow0(frame+1,:,:)
    slow0(frame+model_param%nz_mod+i1,:,:)=slow0(frame+model_param%nz_mod,:,:)
  end do
  model_param%nz_mod=model_param%nz_mod+frame*2
!##################
!CONSTRUCT FILTERS#
!##################
 call from_param('lap_lags',lap_lags,0)
 call from_param('steer_lags',steer_lags,0)
 call from_param('tri_width',tri_width,2.5)
 call from_param('lap_factor',lap_factor,0.)
 call from_param('steer_factor',steer_factor,.999)
 call from_param('eps',eps,1.)
 write(0,*) "width lap_factor lap_lags",tri_width,lap_factor,lap_lags
 write(0,*) "steerlags steer_factor eps",steer_lags,steer_factor,eps
 call init_filters(lap_scale=lap_factor,steer_scale=steer_factor&
    &,lap_lags=lap_lags,steer_lags=steer_lags,triangle_width=tri_width&
    &,run_triangle=.true.)
!############
!Angles     #
!############
 call from_aux("angles","n2",temp2i)
 call from_aux("angles","n1",tempi)
 if (tempi+frame*2.ne.model_param%nz_mod .or. temp2i.ne.model_param%nx_mod) then
      write(0,*) "Angles",tempi,temp2i
      write(0,*) "model",model_param%nz_mod,model_param%nx_mod
      call erexit("angles not the same size as model")
 end if
 allocate(angles((tempi+frame*2)*temp2i),tangle(tempi,temp2i))
 call sreed('angles',tangle,tempi*temp2i*4)
 do i2=1,temp2i
   do i1=1,tempi
     angles(i1+frame+(i2-1)*(2*frame+tempi))=tangle(i1,i2)
   end do
   do i1=1,frame
     angles(i1+(i2-1)*(2*frame+tempi))=angles(frame+1+(i2-1)*(2*frame+tempi))
   end do
   do i1=1,frame
      angles(i1+frame+tempi+(i2-1)*(2*frame+tempi))= angles(frame+tempi+(i2-1)*(2*frame+tempi))
   end do
 end do
 filt=new_steering((/model_param%nz_mod,model_param%nx_mod/),angles,1)
 deallocate(angles,tangle)
!############
!ARRAY SIZES#
!############
  nelem_model=model_param%ndim_mod*model_param%nz_mod*model_param%nx_mod*model_param%ny_mod
!###############
!INIT OPERATORS#
!###############
  write(0,*) "initialize operators"
  call helix_chain_init(nelem_model)
  call polydiv_init(filt)
!#################
!RUN THE OPERATOR#
!#################
  if (exist_file("mask")) then
     call from_aux("mask","n2",temp2i)
     call from_aux("mask","n1",tempi)
     if (tempi+frame*2.ne.model_param%nz_mod .or. temp2i.ne.model_param%nx_mod) then
        write(0,*) "mask",tempi,temp2i
        write(0,*) "model",model_param%nz_mod,model_param%nx_mod
        call erexit("angles not the same size as model")
     end if
     allocate(mask((tempi+frame*2)*temp2i),mask_in(tempi,temp2i))
     call sreed('mask',mask_in,tempi*temp2i*4)
     do i2=1,temp2i
       do i1=1,tempi
         mask(i1+frame+(i2-1)*(2*frame+tempi))=mask_in(i1,i2)
       end do
       do i1=1,frame
         mask(i1+(i2-1)*(2*frame+tempi))=mask(frame+1+(i2-1)*(2*frame+tempi))
       end do
       do i1=1,frame
          mask(i1+frame+tempi+(i2-1)*(2*frame+tempi))= mask(frame+tempi+(i2-1)*(2*frame+tempi))
       end do
     end do
     deallocate(mask_in)
     call mask_init(mask)
  end if
  if (type==1) then
     call inv_tomo (param,model_param,nelem_model,slow,slow0,epspoly,polydiv_lop,i_op,mask_op)
  else if(type==2) then
     call inv_tomo (param,model_param,nelem_model,slow,slow0,epspoly,i_op,polydiv_lop,mask_op)
  end if
  if(allocated(mask)) then
    deallocate(mask)
    call mask_close()
  end if
!#################
!WRITE RESULTS   #
!#################
  model_param%nz_mod=model_param%nz_mod-frame*2
  call srite ("outfile",slow(frame+1:model_param%nz_mod+frame,:,:),4*model_param%nz_mod*model_param%nx_mod*model_param%ny_mod)
end program

!  allocate(tempb(nelem_model),tempc(nelem_model))
!  tempi=0
!  do i2=1,model_param%nx_mod
!     do i1=1,model_param%nz_mod
!        tempi= tempi+1
!        tempb(tempi)=d_slow(i1,i2,1)
!     end do
!  end do
!  ierr=polydiv_lop(.false.,.false.,tempb,tempc)
!  tempi=0
!  do i2=1,model_param%nx_mod
!     do i1=1,model_param%nz_mod
!        tempi= tempi+1
!        tempa(i1,i2,1)=tempc(tempi)
!     end do
!  end do
!  slow=slow0+tempa
!  slow=1./slow
!  deallocate(tempa,stat=ierr)
!  if (0.ne. putch('n1','d',model_param%nz_mod )) call erexit('Trouble writing  n1 to history  ')
!  if (0.ne. putch('n2','d',model_param%nx_mod )) call erexit('Trouble writing  n2 to history  ')
!  if (0.ne. putch('n3','d',3 )) call erexit('Trouble writing  n3 to history  ')
!   allocate (tempa(2*frame+model_param%nz_mod,model_param%nx_mod,model_param%ny_mod))
