program stacking_nmaps
  
  implicit none
  integer::i,ix1,ix2,iargc,npixel,fpix,nbuf,readwrite,blocksize,status,na,group,nfound
  integer::naxis,bitpix,morekeys,keysexist,keysadd
  integer,allocatable::unit(:),naxes(:,:)
  real,allocatable::array(:,:,:),buf(:)
  real::nullval
  character(len=80)::arg,comment,ci,keyval
  character(len=80),allocatable,dimension(:)::file
  logical anynull,simple,extend
  
  na = iargc()
  
  allocate(unit(1:na))
  readwrite = 0
  allocate(file(na-1))
  do i = 1,na-1
     status = 0
     call getarg(i,arg)
     file(i)=arg
     call FTGIOU(unit(i),status)
     call FTOPEN(unit(i),arg,readwrite,blocksize,status)
  end do
  readwrite = 1
  call getarg(na,arg)
  call FTGIOU(unit(na),status)
  call deletefile(arg,status)
  call FTINIT(unit(na),arg,blocksize,status)


  allocate(naxes(1:2,1:na))
  do i = 1,na-1
     call FTGKNJ(unit(i),'NAXIS',1,2,naxes(1:2,i),nfound,status)
  end do

  ! check size of fits image
  do i = 1,na-2
     if(naxes(1,i)/=naxes(1,i+1))goto 999
     if(naxes(2,i)/=naxes(2,i+1))goto 999
  end do

  npixel = naxes(1,1)*naxes(2,1)
  group = 1
  nullval = -999
  
  allocate(buf(1:naxes(1,1)))
  allocate(array(1:na,1:naxes(1,1),1:naxes(2,1)))
  do i = 1,na-1
     fpix = 1
     nbuf = naxes(2,1)
     do ix1 = 1,naxes(1,1)
        call FTGPVE(unit(i),group,fpix,nbuf,nullval,buf,anynull,status)
        fpix = fpix + nbuf
        array(i,ix1,1:naxes(2,1))=buf(1:naxes(2,1))
     end do
  end do
  
  do ix1 = 1,naxes(1,1)
     do ix2 = 1,naxes(2,1)
        do i = 1,na-1
           array(na,ix1,ix2)=array(na,ix1,ix2)+array(i,ix1,ix2)/dble(na-1)
        end do
     end do
  end do


  !!! copy header keywords
  call FTCPHD(unit(1) , unit(na), status)

  !!! add history
  comment = "Created by cygnus:work/ISWANA/src/stackn.f90"
  call  FTPHIS(unit(na), comment, status)
  write(comment,"(I3)")na-1
  comment = trim(comment)//" files are stacked--"
  call  FTPHIS(unit(na), comment, status)
  do i = 1,na-1
  write(ci,"(I3.3)")i
  comment = "FILE"//trim(ci)//" : "//trim(file(i))
  call  FTPHIS(unit(na), comment, status)
  end do


  !!! write image
  fpix = 1
  do ix1 = 1,naxes(1,1)
     fpix = (ix1-1)*naxes(2,1)+1
     call FTPPRE(unit(na),group,fpix,naxes(2,1),array(na,ix1,1:naxes(2,1)),status)
  end do
  
  do i = 1,na
     call FTCLOS(unit(i),status)
     call FTFIOU(unit(i),status)
  end do
  deallocate(unit,array,naxes,file)

  stop
999 stop "image sizes must be same for all images"

end program stacking_nmaps

! *************************************************************************
subroutine deletefile(filename,status)
  integer status,unit,blocksize
  character*(*) filename
  if (status > 0)return
  call ftgiou(unit,status)
  call ftopen(unit,filename,1,blocksize,status)
  if (status == 0)then
     call ftdelt(unit,status)
  else if (status .eq. 103)then
     status=0
     call ftcmsg
  else
     status=0
     call ftcmsg
     call ftdelt(unit,status)
  end if  
  call ftfiou(unit, status)
  return
end subroutine deletefile