program aa
  use sep
  integer :: n(5)
  integer :: i,npts,ipt,ia,ib,i1,i2
  real :: o1,o2,d1(5),d2(5)
  real :: mymax(2)
  character(len=128) :: label
  real, allocatable :: map1(:,:),map2(:,:),count(:,:)
  real, allocatable :: pts(:,:)
  real, allocatable :: err(:)

  call sep_init()
  call from_param("max1",mymax(1))
  call from_param("max2",mymax(2))
  call from_param("min1",o1)
  call from_param("min2",o2)

  call from_history("n2",npts)
    
  n=(/10,20,40,80,160/)
  
  do i=1,5
    d2(i)=(mymax(2)-o2)/(n(i)-1)
    d1(i)=(mymax(1)-o1)/(n(i)-1)
  end do
  call sep_put_data_axis_par("out",1,n(5),o1,d1(5),"z")
  call sep_put_data_axis_par("out",2,n(5),o2,d2(5),"x")
  allocate(pts(2,npts),err(npts))

write(0,*) "CHECK ME",npts

  call sreed("in",pts,npts*8)
  call sreed("err",err,4*npts)

  do i=1,5
     allocate(map1(n(i),n(i)),count(n(i),n(i)))
     count=0; map1=0;

     do ipt=1,npts
       i1=max(1,min(n(i),nint((pts(1,ipt)-o1)/d1(i))))
       i2=max(1,min(n(i),nint((pts(2,ipt)-o2)/d2(i))))
       count(i1,i2)=count(i1,i2)+1
       map1(i1,i2)=map1(i1,i2)+err(ipt)
     end do
     if(i>1) then
       ib=0
       do i2=1,n(i)
         if(mod(i2,2)==1) ib=ib+1
         ia=0
         do i1=1,n(i)
            if(mod(i1,2)==1) ia=ia+1
            if(count(i1,i2) >0) then 
              map1(i1,i2)=map1(i1,i2)/count(i1,i2)
            else
              map1(i1,i2)=map2(ia,ib)
            end if
         end do 
      end do
      deallocate(map2)
      else 
        where(count>0) map1=map1/count
    end if 
call srite("xx.H",map1,size(map1)*4)
    allocate(map2(n(i),n(i)))
    map2=map1
    deallocate(map1,count)
  end do
  call srite("out",map2,size(map2)*4)


end program
