c     release
c
c     Example code for using the 2dFGRS survey masks.
c
c     Given ra and dec in radians (B1950), it returns the completeness, the 
c     magnitude limit, the mu value and the APM plate number of that 
c     particular area.
c     N.B.:
c      - when compl = -1., the given position is outside the main 2dF boundary
c      - when compl = -2., the given position is inside a drill hole 
c      
c     Usage: release ra dec  (where ra and dec are in radians)

      program release
*************************************variables********************************
      implicit none

      integer iplate,ix,iy
      real ra,dec,compl,maglim,muval
      real*8 decdble,radble
      character reg*3
******************************************************************************

      call getdble(1,radble)
      call getdble(2,decdble)

      ra=real(radble)
      dec=real(decdble)

      call which_reg(ra,dec,reg)
      write(0,*) 'Region: ',reg
      call in_2df_mask_reg(radble,decdble,reg,compl,maglim,muval,ix,iy) 
      call fnumber(ra,dec,iplate)

      write(0,*) 'Position: ra = ',ra,' ; dec = ',dec
      write(0,*) 'Compl:',compl,' maglim:',maglim,' muval:',muval
     &,          ' plate:',iplate

      end

c-----------------------------------------------------------------------------
c     Reads an input as a real

      subroutine getdble(i,d)
*************************************variables********************************
      implicit none
      integer i
      real*8 d
      character a*100
******************************************************************************

      call getarg(i,a)
      read(a,*) d

      return
      end

c----------------------------------------------------------------------------
c
c                       PUBLIC MASK SOFTWARE SUBROUTINES
c
c----------------------------------------------------------------------------
c     Given ra & dec determines which region to look at
      subroutine which_reg(ra,dec,reg)
*************************************variables********************************
      implicit none

      real ra,dec
      character reg*3

      real EPS
      parameter (EPS=1.e-5)
      integer i,ifirst,nb_strip,nb_stripsgp,nb_stripngp
      real ra_min(4),ra_max(4),dec_min(4),dec_max(4)
      data ifirst/1/
      save ra_min,ra_max,dec_min,dec_max,nb_strip,nb_stripsgp
     &,    nb_stripngp,ifirst
******************************************************************************

      if (ifirst.eq.1) then
         ifirst=0
         nb_stripsgp=3
         nb_stripngp=4
         nb_strip=5
c        SGP strips: Thanks Andrew H. for finding this bug.
         ra_min(3)= 5.707227    ! 21h48m
         ra_max(3)= 0.890118    !  3h24m
         dec_min(1)= -0.479966  ! -27.5 deg
         dec_max(1)= -0.392699  ! -22.5 deg
         ra_min(2)= 5.670138    ! 21h39.5m
         ra_max(2)= 0.975203    !  3h43.5m
         dec_min(2)= -0.567232  ! -32.5 deg
         dec_max(2)= -0.479966  ! -27.5 deg
         ra_min(1)= 5.711590    ! 21h49m
         ra_max(1)= 0.911935    !  3h29m
         dec_min(3)= -0.654498  ! -37.5 deg
         dec_max(3)= -0.567232  ! -32.5 deg
c        NGP strips
         ra_min(4)= 2.574361    !  9h50m
         ra_max(4)= 3.883358    ! 14h50m
         dec_min(4)= -0.130900  !  -7.5 deg
         dec_max(4)= +0.043633  !  +2.5 deg
c        Make boundaries more secure
         do i=1,4
            ra_min(i)=ra_min(i)-EPS
            ra_max(i)=ra_max(i)+EPS
            dec_min(i)=dec_min(i)-EPS
            dec_max(i)=dec_max(i)+EPS
         enddo
         call set_pivals
      endif

      i=0
      do while (i.lt.nb_stripsgp)
         i=i+1
         if (((ra_min(i).le.ra).or.(ra_max(i).ge.ra)).and.
     &        (dec_min(i).le.dec).and.(dec_max(i).ge.dec)) then
            i=nb_strip
            reg='sgp'
         endif
      enddo

      do while (i.lt.nb_stripngp)
         i=i+1
         if ((ra_min(i).le.ra).and.(ra_max(i).ge.ra).and.
     &        (dec_min(i).le.dec).and.(dec_max(i).ge.dec)) then
            i=nb_strip
            reg='ngp'
         endif
      enddo
      
      if (i.ne.nb_strip) reg='ran'

      return
      end

c-----------------------------------------------------------------------------
c     For a given position in ra & dec, this subroutine tells if this
c     position is inside the 2df_mask or not, by returning the actual
c     completeness. It returns in the same time the corresponding 
c     magnitude limit, and mu mask value.

      subroutine in_2df_mask_reg(ra,dec,reg,compl,maglim,muval,ix,iy)
*************************************variables********************************
      implicit none

      real compl,maglim,muval
      real*8 ra,dec

      integer np_x,np_y
      parameter (np_x=2*1800,np_y=2*600)
      integer np_xx,np_yy,ifirst,ix,iy,np,sgn
      real rpix(np_x,np_y),magpix(np_x,np_y),mupix(np_x,np_y)
      real*8 x,y,z,rx,ry
      
      character name_mask*20,last_reg*3,reg*3
      save ifirst,np_xx,np_yy,np,rpix,magpix,mupix,last_reg
      data ifirst/1/
******************************************************************************

      if (reg.eq.'ran') then
         call in_2df_ran(ra,dec,compl,maglim,muval,ix,iy)
         return
      endif

c     On the first call, the program reads the mask
      if ((ifirst.eq.1).or.(last_reg.ne.reg)) then
         ifirst=0
         last_reg=reg
         name_mask='mask.'//reg//'.dat'
         call read_mask(rpix,np_xx,np_yy,name_mask,np_x,np_y,np)
         name_mask='maglim.'//reg//'.dat'
         call read_mask(magpix,np_xx,np_yy,name_mask,np_x,np_y,np)
         name_mask='mumask.'//reg//'.dat'
         call read_mask(mupix,np_xx,np_yy,name_mask,np_x,np_y,np)
         write(0,*) 'Read masks for ',reg,' region!'
      endif

      call radec_xyz_dbl(ra,dec,x,y,z)
      call eq_2dfrx(np,x,y,z,rx,ry,sgn)
      call rx_ix_map(rx,ry,ix,iy)

      if ((ix.ge.1).and.(ix.le.np_xx).and.
     :    (iy.ge.1).and.(iy.le.np_yy)) then 
         compl=rpix(ix,iy)
         maglim=magpix(ix,iy)
         muval=mupix(ix,iy)
      else
         compl= -1.0  !default value (i.e. outside 2dF boundary)    
         maglim=-2.0  !default value (i.e. outside 2dF boundary)
         muval=19.0   !default value (i.e. outside 2dF boundary)
      endif
                
      return
      end

c-----------------------------------------------------------------------------
c     For a given position in ra & dec, this subroutine tells if this 
c     position is inside one of the used random fields or not, by returning 
c     the actual completeness of the field.

      subroutine in_2df_ran(ra,dec,compl,maglim,muval,ix,iy)
*************************************variables********************************
      implicit none

      integer ix,iy
      real*8 ra,dec
      real compl,maglim,muval

      integer nmax,NP,NPX,NPY
      real theta_min
      parameter (nmax=98,NP=2*2400,NPX=NP/48,NPY=NP/48,theta_min=1.014)
      logical in
      integer i,ifirst,nb_field
      real*8 costheta_r,r,x,y,z,costheta_min,costheta
     &,      xc(nmax),yc(nmax),zc(nmax)
      real comp(nmax),magpix(nmax,NPX,NPY),mupix(nmax,NPX,NPY)
      character name_mask*20
      data ifirst/1/
      save ifirst,nb_field,xc,yc,zc,comp,costheta_min,magpix,mupix
******************************************************************************

c     On the first call, reads the file containing the random fields
      if (ifirst.eq.1) then
         ifirst=0
         name_mask='mask.ran.dat'
         call ran2df_cen_used_dbl(nb_field,xc,yc,zc,comp,name_mask)
         name_mask='maglim.ran.dat'
         call readranmask(name_mask,NPX,NPY,nmax,magpix)
         name_mask='mumask.ran.dat'
         call readranmask(name_mask,NPX,NPY,nmax,mupix)
         costheta_min=dble(cos(theta_min*atan(1.)/45.))
         write(0,*) 'Read mask for ran region!'
      endif

      i=0
      in=.false.
      call radec_xyz_dbl(ra,dec,x,y,z)
      r=sqrt(x**2+y**2+z**2)
      costheta_r=costheta_min*r
      do while ((.not.in).and.(i.lt.nb_field))
         i=i+1
         costheta=(x*xc(i)+y*yc(i)+z*zc(i))
         if (costheta.gt.costheta_r) then   ! We use here the fact the random
            in=.true.                       ! fields doesn't overlapp with 
         endif                              ! each other...
      enddo
      
      if (in) then
         call get_pixelran(NP,NPX,NPY,x,y,z,xc(i),yc(i),zc(i),maglim
     &,                    muval,nmax,magpix,mupix,i,ix,iy)

         if (maglim.le.-1) then 
            compl=-2.           ! default value in holes
         else
            compl=comp(i)
         endif
      else
         compl=-1.              ! default value (ie. outside 2dF boundary)
         maglim=-2.             ! default value (when not specified)
         muval=19.0             ! default value (when not specified)
      endif

      return
      end

c----------------------------------------------------------------------------
c     Given x,y,z position and centre of the corresponding random field,
c     returns the magnitude limit at that position
      subroutine get_pixelran(NP,NPX,NPY,x,y,z,xc,yc,zc,maglim
     &,                       muval,nmax,magpix,mupix,ifield,ix,iy)
      
      implicit none

      integer NP,NPX,NPY,nmax,ifield
      real*8 x,y,z,xc,yc,zc
      real maglim,magpix(nmax,NPX,NPY),muval,mupix(nmax,NPX,NPY)
      
      integer sgn,ixc_min,iyc_min,ix,iy
      real rx,ry


      call eq_2dfrx(NP,xc,yc,zc,rx,ry,sgn)
      call bound_map_ran(NPX,NPY,rx,ry,ixc_min,iyc_min)

      call eq_2dfrx(NP,x,y,z,rx,ry,sgn)
      call rx_ix_map(rx,ry,ix,iy)
      
      ix=min(max(1,ix-ixc_min),NPX)
      iy=min(max(1,iy-iyc_min),NPY)

      if ((ix.ge.1).and.(ix.le.NPX).and.(iy.ge.1).and.(iy.le.NPY)) then  
         maglim=magpix(ifield,ix,iy)
         muval=mupix(ifield,ix,iy)
      else
         maglim=-2.0
         muval=19.0
      endif

      return
      end

c-----------------------------------------------------------------------------
c This subroutine transforms ra & dec into x,y & z cartesian coordinates.
      subroutine radec_xyz_dbl(ra,dec,x,y,z)

      implicit none

      real*8 x,y,z,ra,dec
      
      x=cos(dec)*cos(ra)
      y=cos(dec)*sin(ra)
      z=sin(dec)
      
      return
      end

c-----------------------------------------------------------------------------
c     Convert an equatorial cartesian coordinate to the rx,ry pixel
c     coordinate used for the 2df mask. This subroutine is valid for both 
c     sgp and ngp elements.
c
c
      subroutine  eq_2dfrx(NP,x,y,z,rx,ry,sgn)
*************************************variables********************************
      implicit none
      real*8 x,y,z,xx,yy,zz,r,phi,costheta,sintheta_2
      real PI,EPS
      parameter (EPS=3.e-13)
      integer ifirst,sgn,NP,ix_min,iy_min
      real racen,sc,deccen,rx,ry
      real xc,yc,zc,dec_gc,ra_gc,xgc,ygc,zgc
     & ,xgp,ygp,zgp,phi0
      save ifirst,sc,xc,yc,zc,xgc,ygc,zgc,xgp,ygp,zgp,phi0,PI !ix_min,iy_min
      data ifirst/1/
******************************************************************************

c     On the first call define the values which define the
c     centre, orientation and scale of the transformation
      if (ifirst.eq.1) then
         ifirst=0
         PI= atan(1.)*4.
c        This is the direction of the Galactic Centre as adopted
c        by Steve Maddox
         racen  = 12.3*PI/180.0  
         deccen = -27.5*PI/180.0  
         sc = 120.0*PI/180.0
         xc=cos(deccen)*cos(racen)
         yc=cos(deccen)*sin(racen)
         zc=sin(deccen)
c        This is a direction perpendicular to above but otherwise
c        arbitrary provided that the phi0 has then been correctly
c        set to be the offset between this arbitrary direction and
c        that used for the projection of the APM SGP region 
         ra_gc   = -94.40593*PI/180.0 
         dec_gc  = -28.90771*PI/180.0 
         phi0=236.97694*PI/180.0
         xgc=cos(dec_gc)*cos(ra_gc)
         ygc=cos(dec_gc)*sin(ra_gc)
         zgc=sin(dec_gc)
c        Generate the mutually perpendicular unit vector via
c        the cross product
         xgp = yc*zgc-zc*ygc
         ygp = zc*xgc-xc*zgc
         zgp = xc*ygc-yc*xgc
      endif
     
c     Section utilized on each call
c     Find components in the new 3D Cartesian system      
      r=sqrt(x**2+y**2+z**2)
      zz=(x*xc+y*yc+z*zc)
      xx=(x*xgc+y*ygc+z*zgc)
      yy=(x*xgp+y*ygp+z*zgp)

c     Compute corresponding spherical polar angles
      costheta= abs(zz/r)
      sgn=int(sign(1.,real(zz/r))) ! tells if x,y,z is in the ngp or sgp region
      phi=phi0-atan2(yy,xx)
      if (phi.gt.2.0*PI) phi=phi-2.0*PI

c     Apply the Zenithal Equal Area Projection
      if (real(dble(1.0)-costheta).ge.EPS) then
         sintheta_2=dsqrt(dble(0.5)*(dble(1.0)-costheta))
      else
         sintheta_2=0.
      endif
      r=2.*sintheta_2 !If sc changes, change this like:r=sintheta_2/sin(sc/4.)

c     convert to x,y coordinate
      xx=-r*sin(phi)*dble(NP)*0.5
      yy= r*cos(phi)*dble(NP)*0.5
      
      call bound_map_2df(NP,sgn,ix_min,iy_min)

      rx=real(xx+0.5*dble(NP)-1.*dble(ix_min)) 
      ry=real(-yy+0.5*dble(NP)-1.*dble(iy_min))

      return
      end

c-----------------------------------------------------------------------------
c     This subroutine contains the boundaries of the sgp and ngp pixel map.
c     It returns the ix_min and iy_min which are the offsets of the two maps.

      subroutine bound_map_2df(NP,sgn,ix_min,iy_min)
******************************************************************************
      implicit none

      integer NP,sgn,ix_min,iy_min
******************************************************************************

      if (sgn.eq.1) then        ! use the sgp map
         ix_min=NP*15/100
         iy_min=NP*44/100
      else if (sgn.eq.-1) then  ! use the ngp map
         ix_min=NP*7/100
         iy_min=NP*60/100
      endif

      return
      end

c----------------------------------------------------------------------------
c This subroutine gives the right boundaries for each random field mask, ie.
c calculates the offset needed such that the pixel center of each random 
c field, given by (rxc,ryc), is located at (NPX_ran/2, NPY_ran/2) .
c
      subroutine bound_map_ran(NPX_ran,NPY_ran,rxc,ryc,ixc_min,iyc_min)

      implicit none

      integer NPX_ran,NPY_ran,ixc_min,iyc_min
      real rxc,ryc

      integer ixc,iyc

      
      call rx_ix_map(rxc,ryc,ixc,iyc)
      ixc_min=ixc-int(NPX_ran/2)
      iyc_min=iyc-int(NPY_ran/2)

      return
      end

c-----------------------------------------------------------------------------
c     This subroutine converts real pixel coordinates rx,ry to integer 
c     ix,iy pixel coordinates. 

      subroutine rx_ix_map(rx,ry,ix,iy)
******************************************************************************
      implicit none
      
      integer ix,iy
      real rx,ry
******************************************************************************

      ix=int(rx+0.5)
      iy=int(ry+0.5)
         
      return
      end
      
c-----------------------------------------------------------------------------
c     This subroutine sets and checks that NP and NPX_NPY are correct and that
c     they correspond to the value used for NPX and NPY, the dimensions of the
c     pixel grid.
c     N.B.: NPX_NPY is also accepted if equal to 0 (default value if not 
c     used!)

      subroutine NP_set(NPX,NPY,NP,NPX_NPY)
******************************************************************************
      implicit none
      
      integer NPX,NPY,NP,NPX_NPY

      integer num
******************************************************************************

      NP=NPX*4/3 
      if (NP.ne.NPY*4) then
         write(0,*) 'Distortion of the pixel grid: NPX= ',NPX,' NPY='
     &,              NPY,' do not obey NPX=3*NPY!' 
         stop 'NP_set'
      endif

c     We check also that NPX_NPY is NPX*NPY or NPX_NPY = 0
      num=NPX*NPY
      if ((num.ne.NPX_NPY).and.(NPX_NPY.ne.0)) then
         write(0,*) 'Wrong dimension on NPX_NPY; should be ',num
         stop 'NP_set'
      endif

      return
      end

c----------------------------------------------------------------------------
c  Subroutine to read the 2dF mask
c                           
      subroutine  read_mask(mask,np_xx,np_yy,name_mask,np_x,np_y,np)
*************************************variables********************************
      implicit none

      integer np_x,np_y,np_xx,np_yy,np
      real mask(np_x,np_y)
      character name_mask*(*)

      integer ix,iy
******************************************************************************

      open(11,file=name_mask,status='old',form='unformatted')
      read(11) np_xx,np_yy
      call NP_set(np_xx,np_yy,np,0)
      if ((np_xx.le.np_x).and.(np_yy.le.np_y)) then
         do ix=1,np_xx
            read(11) (mask(ix,iy),iy=1,np_yy)
         enddo
      else
         write(0,*) 'Dimension of the mask too big!'
         write(0,*) ' Mask: NPX= ',np_xx,' NPY= ',np_yy
         stop 'read_mask!'
      endif
      close(11)

      return
      end

c-----------------------------------------------------------------------------
      subroutine readranmask(namemask,NPX,NPY,nbf,rpix)

c     This subroutine reads the mask stored in the file name_mask.
c     N.B.: This subroutine contains the same stuff as mask_2df, but has the 
c     advantage of giving back also sgn!

      implicit none

      integer nbf,NPX,NPY
      real rpix(nbf,NPX,NPY)
      character namemask*(*)
      
      integer i,j,k,npxx,npyy,nbff


      open(unit=14,file=namemask,status='old',form='unformatted')
      read(14) nbff,npxx,npyy
      if ((npxx.ne.NPX).or.(npyy.ne.NPY).or.(nbff.ne.nbf)) then
        write(*,*) 'Wrong dimensions for ran mask. Should have:'
     &,             npxx,npyy,nbff
        stop 'readranmask'
      endif
      do k=1,nbff
         do i=1,npxx
            read(14) (rpix(k,i,j),j=1,npyy)
         enddo
      enddo
      close(14)

      return
      end

c----------------------------------------------------------------------------
c     This subroutine reads the used 2dF random field centres from the file 
c     name, which contains also the overall completeness of the field (with 
c     respect to the underlying density field). It is a real*8 version 
c     of ran2df_cen_used.
      subroutine ran2df_cen_used_dbl(nb_field,xc,yc,zc,comp,name)
*************************************variables********************************
      implicit none

      integer nb_field
      real comp(*)
      real*8 xc(*),yc(*),zc(*)
      character name*(*)

      integer i
******************************************************************************

      open(unit=11,file=name,status='old',form='unformatted')
      read(11) nb_field
      read(11) (comp(i),i=1,nb_field)
      do i=1,nb_field
         xc(i)=dble(comp(i))
      enddo
      read(11) (comp(i),i=1,nb_field)
      do i=1,nb_field
         yc(i)=dble(comp(i))
      enddo
      read(11) (comp(i),i=1,nb_field)
      do i=1,nb_field
         zc(i)=dble(comp(i))
      enddo
      read(11) (comp(i),i=1,nb_field)
      close(11)

      return
      end

c-----------------------------------------------------------------------------
c
c     DO NOT ADD THE FOLLOWING ROUTINES IF THIS FILE IS COMBINED WITH HOLES.F
c
c-----------------------------------------------------------------------------
      subroutine fnumber(rar,decr,ifield)
c
c     given ra,dec in radians returns corresponding field no.
c     
      implicit none

      integer ifield
      real decr,rar

      integer iygrid,ixgrid,i
      real ra,dec,raminus,yspace,ygrid,xspace,xgrid,xoff,yoff
     &,    xspace348,ratest348,edge348

      integer*4 field0h(22)
      real*4 ragaps(22), decbands(22)

      data (ragaps(i),i=1,22)/ 0.,144.,90.,66.,52.,44.,38.,33.,30.,28.,
     +                 26.,24.,23.,22.,21.,20.,20.,20.,20.,20.,20.,20./
      data (decbands(i),i=1,22)/-90.,-85.,-80.,-75.,-70.,-65.,-60.,-55.,
     +   -50.,-45.,-40.,-35.,-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15./
      data (field0h(i),i=1,22)/1,2,12,28,50,78,111,149,193,241,293,349,
     +   409,472,538,607,679,751,823,895,967,1039 / 

      real pival,degs_rads,ramins_rads,radfac,tpi
      common/pivals/ pival,degs_rads,ramins_rads,radfac,tpi


c     find the nearest grid centre

      ra = rar*180./pival
      dec = decr*180./pival
      
      if (dec.gt.15.) then
         ifield = 0
         return
      end if
      
      if (ra.lt.0.) ra = ra + 360.
      if (ra.gt.180.) then
         raminus = ra - 360.
      else
         raminus = ra
      end if
      
      yspace = 5.0
      iygrid = nint(dec/5.)+19
      ygrid = decbands(iygrid)
      
      xspace = ragaps(iygrid)
      if (abs(raminus).lt.xspace/4.0/2.0) then ! 4 for degs, 2 for half field
c     it's a 0hr field
         xgrid = 0.
         ixgrid = 0
      else 
         xspace = ragaps(iygrid) / 4.0 ! in degs
         ixgrid = int(ra/xspace + 0.5)
         xgrid = ixgrid * xspace
      end if
      
      ifield = field0h(iygrid) + ixgrid

c     Do we really need this part???? Apparently for the SGP region...
      xspace348   = 26
      ratest348   = 0.5*xspace348*ramins_rads
      edge348     = -10.*ramins_rads +ratest348
      if(ifield.eq.293 .and. raminus.le.edge348) then
         ifield=348
         xgrid = -10*xspace348
      end if
c     ?????

      xoff = ra - xgrid
      if (xoff.gt.180.) xoff = xoff - 360.
      yoff = dec - ygrid
c     type*,'Field no. = ',ifield
c     type*,'RA , Dec offset = ',xoff,yoff
      
      return
      end

c---------------------------------------------------------------------------
      subroutine set_pivals 
c     Initialisation of some constants for these routines
      implicit none

      real pival,degs_rads,ramins_rads,radfac,tpi
      common/pivals/ pival,degs_rads,ramins_rads,radfac,tpi

      pival       = atan(1.)*4.
      degs_rads   = pival/180.0
      ramins_rads = pival/(12.*60.) 
      radfac      = 2.*pival/360./60.
      tpi         = 2.*pival 
      return
      end