  program  read_magdat

  implicit none
  integer(4) :: nside
  integer(8) :: i, npix
  real(4), allocatable :: wl_data(:,:)
  real(8) :: theta, phi
  character(len=200) :: fname
 
  fname="allskymap_nres12r000.zs16.mag.dat";
  open(12,file=fname,status="old",form='unformatted')

  read(12) nside,npix

  allocate(wl_data(0:npix-1,1:4))

  read(12) wl_data(0:npix-1,1)  ! convergence
  read(12) wl_data(0:npix-1,2)  ! shear1
  read(12) wl_data(0:npix-1,3)  ! shear2
  read(12) wl_data(0:npix-1,4)  ! rotation
  close(12)

  do i=0,npix-1   ! Healpix ring scheme
    !call pix2ang(i,nside,theta,phi)
    write(*,*) i,wl_data(i,1),wl_data(i,2),wl_data(i,3),wl_data(i,4)
  end do

  deallocate(wl_data)


  contains
  subroutine pix2ang(pix,nside,theta,phi) ! convert a pixel index (pix) to an angular position (theta,phi)[rad] in spherical coordinates
  integer(8) :: pix, npix, ncap, i, j, pp, s
  integer(4) :: nside
  real(8) :: theta, phi, ph, z, pi

  pi=4.d0*datan(1.d0)
  npix=12*nside*int8(nside); ncap=2*nside*int8(nside-1)

  if(pix.lt.ncap) then  ! North polar cap
    ph=0.5d0*(pix+1.d0)
    i=int8(sqrt(ph-sqrt(dble(int8(ph)))))+1
    j=pix+1-2*i*(i-1)
    z=1.d0-i*i/(3.d0*nside*nside)
    phi=0.5d0*pi/i*(j-0.5d0)
  elseif(pix.lt.(npix-ncap)) then ! Equatorial belt 
    pp=pix-ncap
    i=int(0.25d0*pp/nside)+nside
    j=mod(pp,4*nside)+1
    s=mod(i+nside,2)+1
    z=4.d0/3.d0-2.d0*i/(3.d0*nside)
    phi=0.5d0*pi/nside*(j-0.5d0*s)
  else    ! South polar cap 
    ph=0.5d0*(npix-pix)  
    i=int8(sqrt(ph-sqrt(dble(int8(ph)))))+1
    j=4*i+1-(npix-pix-2*i*(i-1))
    z=-1.d0+(i*i)/(3.d0*nside*nside)
    phi=0.5d0*pi/i*(j-0.5d0)
  endif
  theta=dacos(z)

  end subroutine

  end program
