!Compilar con /opt/bin/intel_fc80/bin/ifort -o reshuffle reshuffle.f90 -lcfitsio -static-libcxa

program reshuffle
  use IFPOSIX
  implicit none
  integer status,blocksize,unit,unitout,readwrite,naxes(3),hdunum
  integer hdutype,nfound,donde,iscan,idx
  integer i,istk,ikin,dist_off,j,k,count,nburst,nstk,nkin
  character filename*80,fileout*80,nstks*4,nkins*4,fileout2*80
  integer, dimension(:,:,:), allocatable :: perfout,profiles
  integer, dimension(:), allocatable:: array
  logical anyf

  nstk=6
  nkin=5
  
  status=0
  readwrite=0  ! Read only


  donde=ipxfargc()
  if (donde.eq.0) stop 'No file supplied'
  if (donde.ne.3) stop 'No info on number of stokes and number of acqs'
  call pxfgetarg(1,filename,donde,status)
!  open (1,file='bin2sp.cfg',status='old')
!  read (1,*) filename
!  close(1)
 
  call pxfgetarg(2,nstks,donde,status)
  nstk=int(rnum(nstks))
  print*,'Number of Stokes ',nstk
  call pxfgetarg(3,nkins,donde,status)
  nkin=int(rnum(nkins))
  print*,'Number of exps in kinetic series',nkin
  !#########################################################################
  ! All the FITS initialization stuff
 
  call ftgiou(unit,status)  !get free I/O unit for data file
  if (status.ne.0)  then
     call ftrprt('STDOUT',status)
     stop
  endif

  call ftopen(unit,filename,readwrite,blocksize,status) !open FITS file
  if (status.ne.0)  then
     call ftrprt('STDOUT',status)
     stop
  endif
  call ftgknj(unit,'NAXIS',1,3,naxes,nfound,status) !get naxes values from


  if (status.ne.0)  then
     call ftrprt('STDOUT',status)
     stop
  endif
  call ftthdu(unit,hdunum,status) !get number of FITS extensions,
                                  ! i.e. numbers of scans
  if (status.ne.0)  then
     call ftrprt('STDOUT',status)
     stop
  endif

!#################################################################33
! Initialization for the output FITS file

!!! New name for the output FITS file
  
  donde=index(filename,'/',BACK=.TRUE.)
  fileout2=trim(filename(donde+1:len(filename)))
  fileout='test.fits'

   call ftgiou(unitout,status)  !get free I/O unit for data file
  if (status.ne.0)  then
     call ftrprt('STDOUT',status)
     stop
  endif

  
  call ftopen(unitout,fileout,1,blocksize,status)

  if (status .eq. 0)then
     !         file was opened;  so now delete it
     call ftdelt(unitout,status)
  else if (status .eq. 103)then
     !         file doesn't exist, so just reset status to zero and clear errors
     status=0
     call ftcmsg
  else
     !         there was some other error opening the file; delete the file anyway
     status=0
     call ftcmsg
     call ftdelt(unitout,status)
  end if
  call ftfiou(unitout, status)

  call ftgiou(unitout,status)  !get free I/O unit for output file


  call ftinit(unitout,fileout,blocksize,status) !create  FITS output file
  if (status.ne.0)  then
     call ftrprt('STDOUT',status)
     stop
  endif

 
  call ftcphd(unit,unitout,status) ! Write a header for the primary extension
  ! from the header of the input data file
  if (status.ne.0)  then
     print*,'?'
     call ftrprt('STDOUT',status)
     stop
  endif

  call ftgkyj(unit,'NBURST',nburst,'',status)
  call ftmkyj(unitout,'NBURST',(nburst-1)*nkin,'',status)

!##################################################################333

  allocate(profiles(naxes(1),naxes(2),nstk*nkin))
  allocate(perfout(naxes(1),naxes(2),nstk*nkin))
  allocate(array(naxes(1)*naxes(2)*nstk*nkin))


  print*,naxes(1:3)
  do iscan=1,hdunum
     print*,iscan
    
     call ftmahd(unit,iscan,hdutype,status) !move to the iscan extension

!     call ftg3dd(unit,0,0,naxes(1),naxes(2),naxes(1),naxes(2),naxes(3),profiles,anyf,status)     ! read the stokes image
  

     do idx=1,naxes(3),nstk*nkin

 
        call ftgsvj(unit,0,3,naxes,[1,1,idx],[512,512,idx+nstk*nkin-1],[1,1,1],0,array,anyf,status)
        
        count=1
        do i=1,nstk*nkin
           do j=1,512
              do k=1,512
                 profiles(k,j,i)=array(count)
                 count=count+1
              enddo
           enddo
        enddo

       do istk=0,nstk-1

           do ikin=0,nkin-1
              
              perfout(1:naxes(1),1:naxes(2),ikin*nstk+istk+1)=profiles(1:naxes(1),1:naxes(2),istk*nkin+ikin+1)
           enddo
        enddo
           
        count=1
        do i=1,nstk*nkin
            do j=1,512
              do k=1,512
                 array(count)=perfout(k,j,i)
                 count=count+1
              enddo
           enddo
        enddo
        print*,'writing ',idx
        call ftpssj(unitout,0,3,naxes,[1,1,idx],[512,512,idx+nstk*nkin-1],array,status)
     enddo

      if (iscan.gt.1) call ftiimg(unitout,-32,3,(/ naxes(1),naxes(2),naxes(3)/),status) 
 
     
print*,'Scan #',iscan,' of',hdunum     
!!     call ftp3dd(unitout,0,naxes(1),naxes(2),naxes(1),naxes(2),naxes(3),perfout,status) 
                                                         !write intensity
     



  enddo
  print*,'Done: ',fileout


  
! Close the FITS data files
  call ftclos(unit,status)
  call ftclos(unitout,status)
  call ftfiou(-1,status) !-1 = all units freed

  print*,len(trim(fileout2))
  call pxfrename('test.fits',9,fileout2,len(trim(fileout2)),status)
  print*,'Output file ',trim(fileout2),' created'
end program reshuffle

