c geometries.f : subroutines which define the geometry of the
c                sdss 2df_sgp 2df_ngp survey regions.
c
c USAGE : each subroutine takes the arguments (in,xo,yo,zo,r)
c
c  INPUT:
c    xo,yo,zo  (REAL*4) Cartesian galaxy position
c    r         (REAL*4) related distance r=sqrt(xo**2+yo**2+zo**2).  
c
c OUTPUT:
c    in   (LOGICAL)    set to .true.  if the galaxy is in the survey
c                      set to .false. if the galaxy is outside the survey
c
c---------------------------------------------------------------------
      subroutine geometry_sdss(in,xo,yo,zo,r)
****************************variables*********************************
      implicit none
      integer ifirst
      real r,xo,yo,zo,fx1,fx2,fy1,fy2,fz1,fz2,dmajor,
     & smajor,degtorad,sminor,cosfocus,cos1,cos2,acos1,acos2
      real PI
      parameter(PI=3.1415926)
      logical in
      save fx1,fx2,fy1,fy2,fz1,fz2,ifirst,dmajor
      data ifirst/1/
**********************************************************************
      if (ifirst.eq.1) then
	 degtorad=PI/180. 
	 smajor=130*degtorad/2. 
	 dmajor=2.*smajor
	 sminor=110*degtorad/2.
	 cosfocus=cos(smajor)/cos(sminor)
	 fz1=cosfocus
	 fy1=0.0
	 fx1=sqrt(1.-fz1*fz1)
	 fz2=fz1
	 fy2=fy1
	 fx2=-fx1
	 ifirst=0
	 write(0,*) 'sdss boundary unit vectors:'
         write(0,*) fx1,fy1,fz1
	 write(0,*) fx2,fy2,fz2
      end if
c
c     rotate to system in which z-axis points to centre of the survey
      
      in=.false.
d     if (r.le.0.0) then
d        write(0,*) 'geometry() r=',r
d     end if
      cos1=(fx1*xo+fy1*yo+fz1*zo)/r
      cos2=(fx2*xo+fy2*yo+fz2*zo)/r
d     if (abs(cos1).gt.1.0.or.abs(cos2).gt.1.0) 
d    &     write(0,*) 'cos1=',cos1,' cos2=',cos2
      if (cos1.ge.1.0) then
         acos1=0.0   !in case of rounding error
      else if (cos1.le.-1.0) then
         acos1=PI   !in case of rounding error
      else
         acos1=acos(cos1)
      end if
      if (cos2.ge.1.0) then
         acos2=0.0   !in case of rounding error
      else if (cos2.le.-1.0) then
         acos2=PI   !in case of rounding error
      else
         acos2=acos(cos2)
      end if

      if (acos1+acos2.le.dmajor)  in=.true.
c
      return
      end 
c---------------------------------------------------------------
c------------------------------------------------------------------------------
c  2dF SGP strip
c
c subroutine that specifies the geometrical constraints of the
c survey, in = .true. is returned if the coordinate is in the survey
c The coordinates are also transformed to the specified basis and
c returned in xo,yo,zo, if in the survey.
c
c  2dF survey bounds from Steve Maddox 27 May 1997
c
c 
c 781,782,783,784,785,786,787,788,789,790,791,792,793,794,795,
c 853,854,855,856,857,858,859,860,861,862,863,864,865,866,867
c
c and 
c
c 349,350,351,352,353,354,355,356,357,404,405,406,407,408,
c 409,410,411,412,413,414,415,416,417,418,466,467,468,469,470,471,
c 472,473,474,475,476,477,478,479,480,481,532,533,534,535,536,537
c
c
c If I have done my sums right this should correspond to the following
c boundaries defined in 1950 coords. 
c
c In the NGP   dec between +2 30 0 and -7 30 0  
c             ra  between  9 50 0 to 14 50 0 
c
c In the SGP the three dec strips have different ra limits
c
c dec -22 30 0 to -27 30 0   ra between 21 48 00 to 3 24 00 
c dec -27 30 0 to -32 30 0   ra between 21 39 30 to 3 43 30 
c dec -32 30 0 to -37 30 0   ra between 21 49 00 to 3 29 00 
c

c

c
      subroutine geometry_2df_sgp(in,xo,yo,zo,r)
****************************variables*********************************
      implicit none
      integer ifirst
      real r,xo,yo,zo,sindecr,ra,PI
      real dec1,dec2,dec3,dec4,ra1a,ra2a,ra1b,ra2b,ra1c,ra2c,off
      parameter(PI=3.1415926)
      logical in
      real sin1,sin2,sin3,sin4
      save ifirst,sin1,sin2,sin3,sin4,ra1a,ra2a,ra1b,ra2b,ra1c,ra2c,off
      data ifirst/1/
**********************************************************************
c   Set angular boundary limits and evaluate 
c   sin(dec) values on the first call only
      if (ifirst.eq.1) then
         ifirst=0
c
         ra1a=21.8      * PI/12.0
         ra2a= 3.4      * PI/12.0
         ra1b=21.658333 * PI/12.0
         ra2b= 3.725    * PI/12.0
         ra1c=21.816667 * PI/12.0
         ra2c= 3.48333  * PI/12.0
         ra1a=mod(ra1a+pi,2*pi)-pi !switch to interval [-pi,pi]
         ra2a=mod(ra2a+pi,2*pi)-pi 
         ra1b=mod(ra1b+pi,2*pi)-pi !switch to interval [-pi,pi]
         ra2b=mod(ra2b+pi,2*pi)-pi 
         ra1c=mod(ra1c+pi,2*pi)-pi !switch to interval [-pi,pi]
         ra2c=mod(ra2c+pi,2*pi)-pi 
         if (ra1a.gt.ra2a) then
            off=2.0*pi
         else
            off=0.0
         end if
c
         dec1 = -37.5 *PI/180.0
         dec2 = -32.5 *PI/180.0
         dec3 = -27.5 *PI/180.0
         dec4 = -22.5 *PI/180.0
         sin1=sin(dec1)
         sin2=sin(dec2)
         sin3=sin(dec3)
         sin4=sin(dec4)
         write(0,*) 'sin1=',sin1
         write(0,*) 'sin2=',sin2
         write(0,*) 'sin3=',sin3
         write(0,*) 'sin4=',sin4
      end if
c
c  Check galaxy has a declination that overlaps with the region
      
      sindecr=zo
       if ((sindecr-sin1*r)*(sin4*r-sindecr).le.0.0) then
         in=.false.
       else
         in=.true.
       end if
c 
c  If so check that it has an RA within the range of the corresponding
c  declination strip
      if (in) then
         in=.false.
         
         ra=atan2(yo,xo)
         if      ((sindecr-sin1*r)*(sin2*r-sindecr).ge.0.0) then
            if((ra-ra1a)*(ra2a-ra)*(pi-off).ge.0.0) in=.true.
         else if ((sindecr-sin2*r)*(sin3*r-sindecr).ge.0.0) then
            if((ra-ra1b)*(ra2b-ra)*(pi-off).ge.0.0) in=.true.
         else if ((sindecr-sin3*r)*(sin4*r-sindecr).ge.0.0) then
            if((ra-ra1c)*(ra2c-ra)*(pi-off).ge.0.0) in=.true.         
         end if
      end if
c
      return
      end 
c------------------------------------------------------------------------------
c  2dF SGP strip
c
c subroutine that specifies the geometrical constraints of the
c survey, in = .true. is returned if the coordinate is in the survey
c The coordinates are also transformed to the specified basis and
c returned in xo,yo,zo, if in the survey.
c
c  2dF survey bounds from Steve Maddox 27 May 1997
c
c 
c 781,782,783,784,785,786,787,788,789,790,791,792,793,794,795,
c 853,854,855,856,857,858,859,860,861,862,863,864,865,866,867
c
c and 
c
c 349,350,351,352,353,354,355,356,357,404,405,406,407,408,
c 409,410,411,412,413,414,415,416,417,418,466,467,468,469,470,471,
c 472,473,474,475,476,477,478,479,480,481,532,533,534,535,536,537
c
c
c If I have done my sums right this should correspond to the following
c boundaries defined in 1950 coords. 
c
c In the NGP   dec between +2 30 0 and -7 30 0  
c              ra  between  9 50 0 to 14 50 0 
c
c In the SGP the three dec strips have different ra limits
c
c dec -22 30 0 to -27 30 0   ra between 21 48 00 to 3 24 00 
c dec -27 30 0 to -32 30 0   ra between 21 39 30 to 3 43 30 
c dec -32 30 0 to -37 30 0   ra between 21 49 00 to 3 29 00 
c

c

c
      subroutine geometry_2df_ngp(in,xo,yo,zo,r)
****************************variables*********************************
      implicit none
      integer ifirst
      real r,xo,yo,zo,sindecr,ra,PI
      real dec1,dec2,ra1,ra2,off
      parameter(PI=3.1415926)
      logical in
      real sin1,sin2
      save ifirst,sin1,sin2,ra1,ra2,off
      data ifirst/1/
**********************************************************************
c   Set angular boundary limits and evaluate 
c   sin(dec) values on the first call only
      if (ifirst.eq.1) then
         ra1= 9.83333 *     PI/12.0
         ra2=14.83333 *     PI/12.0
         ra1=mod(ra1+pi,2*pi)-pi !switch to interval [-pi,pi]
         ra2=mod(ra2+pi,2*pi)-pi 
         if (ra1.gt.ra2) then
            off=2.0*pi
         else
            off=0.0
         end if
c
         dec1 = -5.0 *PI/180.0
         dec2 =  2.5 *PI/180.0
         ifirst=0
         sin1=sin(dec1)
         sin2=sin(dec2)
         write(0,*) 'sin1=',sin1
         write(0,*) 'sin2=',sin2
      end if
c
c  Check galaxy has a declination that overlaps with the region
      in=.false.
      sindecr=zo
       if ((sindecr-sin1*r)*(sin2*r-sindecr).ge.0.0) then
c        if in dec bounds check RA bounds
          ra=atan2(yo,xo)
         if((ra-ra1)*(ra2-ra)*(pi-off).ge.0.0) in=.true.
       end if
c
      return
      end 
c--------------------------------------------------------------------------
c  Converts galactic coordinate version of x,y,z,
c  to Eclipltic coordinates and then to the corresponding mask
c  bin number and then returns in=.true. if the mask is .true. in 
c  this bin.
c
       subroutine geometry_pscz(in,x,y,z,r)
******************************************************************
       implicit none
       real elon,elat,size,x,y,z,xe,ye,ze,re,re2,r
       integer i,j,n,nbin,maxbin(182),ifirst,istat
       real PI,R2D
       parameter(R2D=57.29577951,PI=3.141592654)
       logical lmask(41167),in
       save ifirst,maxbin,lmask
       data ifirst/1/
       data maxbin
     #  /0,1,7,19,37,62,93,130,173,223,279,341,409,483,563,
     #  650,743,842,947,1058,1175,1298,1427,1561,1701,1847,1999,
     #  2156,2319,2488,2662,2841,3026,3216,3412,3613,3819,4030,
     #  4246,4467,4693,4924,5160,5400,5645,5895,6149,6407,6670,
     #  6937,7208,7483,7762,8045,8332,8623,8917,9215,9516,9821,
     #  10129,10440,10754,11071,11391,11714,12040,12368,12699,
     #  13032,13368,13706,14046,14388,14732,15078,15425,15774,
     #  16124,16476,16829,17183,17538,17894,18251,18609,18967,
     #  19326,19685,20044,20403,20763,21122,21481,21840,22199,
     #  22557,22915,23272,23628,23983,24337,24690,25042,25392,
     #  25741,26088,26434,26778,27120,27460,27798,28134,28467,
     #  28798,29126,29452,29775,30095,30412,30726,31037,31345,
     #  31650,31951,32249,32543,32834,33121,33404,33683,33958,
     #  34229,34496,34759,35017,35271,35521,35766,36006,36242,
     #  36473,36699,36920,37136,37347,37553,37754,37950,38140,
     #  38325,38505,38679,38848,39011,39168,39320,39466,39606,
     #  39740,39869,39992,40109,40220,40325,40424,40517,40604,
     #  40684,40758,40826,40888,40944,40994,41037,41074,41105,
     #  41130,41148,41160,41166,41167/
**********************************************************************
       if (ifirst.eq.1) then !read mask on first call
         do i = 1, 41167
            lmask(i) = .true.
         end do
         open(10,file='mask_pscz.data',status='old',form='formatted')
         istat=0
         do while (istat.eq.0)
            read (10,*,iostat=istat) j
            if (istat.eq.0) then
               if (j.le.1 .or. j.gt.41167) stop 'mask bin out of range'
               lmask(j) = .false.
            end if
         end do
         close(10)
         write(0,*) 'mask has been read from mask_pscz.data'
         ifirst=0
       end if
c
c  Convert x,y,z to ecliptic coordinates
       xe=x !copy
       ye=y
       ze=z
       call gal_to_eclip(1,xe,ye,ze)
       re2=xe**2+ye**2+ze**2
       if (re2.gt.0.0) then
          re=sqrt(re2)
          elon=atan2(ye,xe)
          if (elon.le.0.0) elon=elon+2.0*PI
          elat=asin(ze/re)
c  and then convert to mask lune bin number nbin
          I=int(90.0-ELAT*R2D+2.5)
          N=MAXBIN(I)-MAXBIN(I-1)
          SIZE=360.0/REAL(N)
          J=int(ELON*R2D/SIZE+1.0)
          NBIN=MAXBIN(I-1)+J
c         Look up mask in this bin
          in=lmask(nbin)
       else
          in=.false.
       end if
       RETURN
       END
c--------------------------------------------------------------------
c  These 4 subroutines have been tested against the the
c  STARLINK utility COCO for barycentric 1950 coordinates
c

c------------------------------------------------------------------------------
c Given a cartesian x,y,z coordinate in which the axes are
c aligned with the equatorial coordinate system
c rotate to align with the galactic coordinate system.
c
c  (RA,dec)==>(l,b) 
c     
      subroutine  eq_to_gal(n,x,y,z)
*************************************variables*********************************
      implicit none
      integer n,i,ifirst
      real x(n),y(n),z(n),PI,rot(3,3),xt,yt,zt,ra_ngp,dec_ngp,
     & ra_gc,dec_gc,ra_gp,dec_gp
      parameter(PI=3.141592654)
      save ifirst,rot
      data ifirst/1/
*******************************************************************************
      if (ifirst.eq.1) then
         ifirst=0
         ra_ngp  = 192.24999787*PI/180.0 !RA  of NGP
         dec_ngp =  27.39994929*PI/180.0 !dec of NGP
         ra_gc   = 265.61073954*PI/180.0 !RA  of Galactic Centre
         dec_gc  = -28.91678503*PI/180.0 !dec of Galactic Centre
         ra_gp   = 317.57300839*PI/180.0
         dec_gp  =  48.12347118*PI/180.0
         rot(1,1)=cos(dec_gc)*cos(ra_gc)  
         rot(1,2)=cos(dec_gc)*sin(ra_gc)  
         rot(1,3)=sin(dec_gc)             
         rot(2,1)=cos(dec_gp)*cos(ra_gp)  
         rot(2,2)=cos(dec_gp)*sin(ra_gp)  
         rot(2,3)=sin(dec_gp)             
         rot(3,1)=cos(dec_ngp)*cos(ra_ngp)
         rot(3,2)=cos(dec_ngp)*sin(ra_ngp)
         rot(3,3)=sin(dec_ngp)            
      end if
      do i=1,n
         xt= rot(1,1)*x(i) + rot(1,2)*y(i) + rot(1,3)*z(i)
         yt= rot(2,1)*x(i) + rot(2,2)*y(i) + rot(2,3)*z(i)
         zt= rot(3,1)*x(i) + rot(3,2)*y(i) + rot(3,3)*z(i)
         x(i)=xt
         y(i)=yt
         z(i)=zt
      end do
      return
      end
c-----------------------------------------------------------------------------
c Given a cartesian x,y,z coordinate in which the axis are
c aligned with the galactic coordinate system
c rotate to align with the equatorial coordinate system.
c
c  (l,b) ==> (RA,dec)
c     
      subroutine  gal_to_eq(n,x,y,z)
*************************************variables*********************************
      implicit none
      integer n,i,ifirst
      real x(n),y(n),z(n),PI,rot(3,3),xt,yt,zt,ra_ngp,dec_ngp,
     & ra_gc,dec_gc,ra_gp,dec_gp
      parameter(PI=3.141592654)
      save ifirst,rot
      data ifirst/1/
*******************************************************************************
      if (ifirst.eq.1) then
         ifirst=0
         ra_ngp  = 192.24999787*PI/180.0 !RA  of NGP
         dec_ngp =  27.39994929*PI/180.0 !dec of NGP
         ra_gc   = 265.61073954*PI/180.0 !RA  of Galactic Centre
         dec_gc  = -28.91678503*PI/180.0 !dec of Galactic Centre
         ra_gp   = 317.57300839*PI/180.0
         dec_gp  =  48.12347118*PI/180.0
         rot(1,1)=cos(dec_gc)*cos(ra_gc)  
         rot(1,2)=cos(dec_gc)*sin(ra_gc)  
         rot(1,3)=sin(dec_gc)             
         rot(2,1)=cos(dec_gp)*cos(ra_gp)  
         rot(2,2)=cos(dec_gp)*sin(ra_gp)  
         rot(2,3)=sin(dec_gp)             
         rot(3,1)=cos(dec_ngp)*cos(ra_ngp)
         rot(3,2)=cos(dec_ngp)*sin(ra_ngp)
         rot(3,3)=sin(dec_ngp)            
      end if
      do i=1,n
         xt= rot(1,1)*x(i) + rot(2,1)*y(i) + rot(3,1)*z(i)
         yt= rot(1,2)*x(i) + rot(2,2)*y(i) + rot(3,2)*z(i)
         zt= rot(1,3)*x(i) + rot(2,3)*y(i) + rot(3,3)*z(i)
         x(i)=xt
         y(i)=yt
         z(i)=zt
      end do
      return
      end
c-----------------------------------------------------------------------------
c aligned with the equatorial coordinate system
c rotate to align with the ecliptic coordinate system
c
c  (RA,dec) ==> (elon,elat)
c     
      subroutine  eq_to_eclip(n,x,y,z)
*************************************variables*********************************
      implicit none
      integer n,i,ifirst
      real x(n),y(n),z(n),PI,rot(3,3),xt,yt,zt,ra_nep,dec_nep,
     & ra_ec,dec_ec,ra_ep,dec_ep
      parameter(PI=3.141592654)
      save ifirst,rot
      data ifirst/1/
*******************************************************************************
      if (ifirst.eq.1) then
         ifirst=0
         ra_nep  = 269.99962*PI/180.0  
         dec_nep = 66.554187*PI/180.0  
         ra_ec   = 359.9998359*PI/180.0  
         dec_ec  = -0.00000794*PI/180.0  
         ra_ep   = 89.99995571*PI/180.0
         dec_ep  = 23.44579351*PI/180.0
         rot(1,1)=cos(dec_ec)*cos(ra_ec)  
         rot(2,1)=cos(dec_ec)*sin(ra_ec)  
         rot(3,1)=sin(dec_ec)             
         rot(1,2)=cos(dec_ep)*cos(ra_ep)  
         rot(2,2)=cos(dec_ep)*sin(ra_ep)  
         rot(3,2)=sin(dec_ep)             
         rot(1,3)=cos(dec_nep)*cos(ra_nep)
         rot(2,3)=cos(dec_nep)*sin(ra_nep)
         rot(3,3)=sin(dec_nep)
      end if
      do i=1,n
         xt= rot(1,1)*x(i) + rot(2,1)*y(i) + rot(3,1)*z(i)
         yt= rot(1,2)*x(i) + rot(2,2)*y(i) + rot(3,2)*z(i)
         zt= rot(1,3)*x(i) + rot(2,3)*y(i) + rot(3,3)*z(i)
         x(i)=xt
         y(i)=yt
         z(i)=zt
      end do
      return
      end
c-----------------------------------------------------------------------------
c-----------------------------------------------------------------------------
c Given a cartesian x,y,z coordinate in which the axis are
c aligned with the galactic coordinate system
c rotate to align with the ecliptic coordinate system.
c
c  (l,b) ==> (elon,elat)
c     
      subroutine  gal_to_eclip(n,x,y,z)
*************************************variables*********************************
      implicit none
      integer n,i,ifirst,j,k
      real x(n),y(n),z(n),PI,rot(3,3),xt,yt,zt,ra_ngp,dec_ngp,
     & ra_gc,dec_gc,ra_gp,dec_gp,rot_1(3,3),rot_2(3,3),ra_nep,
     & dec_nep,ra_ec,dec_ec,ra_ep,dec_ep
      parameter(PI=3.141592654)
      save ifirst,rot
      data ifirst/1/
*******************************************************************************
      if (ifirst.eq.1) then
         ifirst=0
c        Generate the two rotation matrices used in the preceeding
c        subroutines to rotate from:
c        i) galactic --> equatorial coords
c        ii) equatorial --> ecliptic coords
c        and form the matric product to yield the matrix that rotates
c        directly from galactic --> ecliptic coordinates
         ra_ngp  = 192.24999787*PI/180.0 !RA  of NGP
         dec_ngp =  27.39994929*PI/180.0 !dec of NGP
         ra_gc   = 265.61073954*PI/180.0 !RA  of Galactic Centre
         dec_gc  = -28.91678503*PI/180.0 !dec of Galactic Centre
         ra_gp   = 317.57300839*PI/180.0
         dec_gp  =  48.12347118*PI/180.0
         rot_1(1,1)=cos(dec_gc)*cos(ra_gc)  
         rot_1(1,2)=cos(dec_gc)*sin(ra_gc)  
         rot_1(1,3)=sin(dec_gc)             
         rot_1(2,1)=cos(dec_gp)*cos(ra_gp)  
         rot_1(2,2)=cos(dec_gp)*sin(ra_gp)  
         rot_1(2,3)=sin(dec_gp)             
         rot_1(3,1)=cos(dec_ngp)*cos(ra_ngp)
         rot_1(3,2)=cos(dec_ngp)*sin(ra_ngp)
         rot_1(3,3)=sin(dec_ngp)            
c
         ra_nep  = 269.99962*PI/180.0  
         dec_nep = 66.554187*PI/180.0  
         ra_ec   = 359.9998359*PI/180.0  
         dec_ec  = -0.00000794*PI/180.0  
         ra_ep   = 89.99995571*PI/180.0
         dec_ep  = 23.44579351*PI/180.0
         rot_2(1,1)=cos(dec_ec)*cos(ra_ec)  
         rot_2(2,1)=cos(dec_ec)*sin(ra_ec)  
         rot_2(3,1)=sin(dec_ec)             
         rot_2(1,2)=cos(dec_ep)*cos(ra_ep)  
         rot_2(2,2)=cos(dec_ep)*sin(ra_ep)  
         rot_2(3,2)=sin(dec_ep)             
         rot_2(1,3)=cos(dec_nep)*cos(ra_nep)
         rot_2(2,3)=cos(dec_nep)*sin(ra_nep)
         rot_2(3,3)=sin(dec_nep)
c
c       Matrix multiplication  rot = rot_2 * rot_1
         do i=1,3
            do j=1,3
               rot(i,j)=0.0
               do k=1,3
                  rot(i,j)=rot(i,j)+rot_2(k,j)*rot_1(i,k)
               end do
            end do
         end do
      end if
c
c   On subsequent calls we eexcute only the following code segment
      do i=1,n
         xt= rot(1,1)*x(i) + rot(2,1)*y(i) + rot(3,1)*z(i)
         yt= rot(1,2)*x(i) + rot(2,2)*y(i) + rot(3,2)*z(i)
         zt= rot(1,3)*x(i) + rot(2,3)*y(i) + rot(3,3)*z(i)
         x(i)=xt
         y(i)=yt
         z(i)=zt
      end do
      return
      end