;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
;
; NAME:
;    EQ2SURVEY
;       
; PURPOSE:
;    Convert from ra, dec to lambda, eta (SDSS survey coordinates)
;
; CALLING SEQUENCE:
;    eq2survey, ra, dec, lambda, eta
;
; INPUTS: 
;    ra: Equatorial latitude in degrees 
;    dec: Equatorial longitude in degrees
;
; OPTIONAL INPUTS:
;    None
;
; KEYWORD PARAMETERS:
;    None
;       
; OUTPUTS: 
;    lambda: Survey longitude in degrees
;    eta: Survey latitude in degrees
;
; OPTIONAL OUTPUTS:
;    None
;
; CALLED ROUTINES:
;    ATBOUND
;    ATBOUND2
; 
; PROCEDURE: 
;    
;	
;
; REVISION HISTORY:
;    Written: 5/15/2000  Erin Scott Sheldon
;                        Taken from astrotools.
;       
;                                      
;-                                       
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

PRO atbound, angle, min, max

  w=where(angle LT min,nw)
  WHILE nw NE 0 DO BEGIN 
      angle[w] = angle[w] + 360.0
      w=where(angle LT min,nw)
  ENDWHILE 

  w=where(angle GE max,nw)
  WHILE nw NE 0 DO BEGIN 
      angle[w] = angle[w] - 360.0
      w=where(angle GT max,nw)
  ENDWHILE 

  return
END 

PRO atbound2, theta, phi

  atbound, theta, -180.0, 180.0
  w = where( abs(theta) GT 90.,nw)
  IF nw NE 0 THEN BEGIN
      theta[w] = 180. - theta[w]
      phi[w] = phi[w] + 180.
  ENDIF 
  atbound, theta, -180.0, 180.0
  atbound, phi, 0.0, 360.0

  w=where( abs(theta) EQ 90., nw)
  IF nw NE 0 THEN phi[w] = 0.0

  return
END 

PRO eq2survey, ra_in, dec_in, slong, slat

  IF n_params() LT 2 THEN BEGIN 
      print,'-Syntax: eq2survey, ra, dec, lambda, eta'
      print,' ra, dec in degrees'
      return
  ENDIF 
  
  ;; slong=lambda
  ;; slat=eta

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; Some parameters
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  at_surveyCenterRa  =  185.0
  at_surveyCenterDec =   32.5

  at_deg2Rad = !dpi/180.
  at_rad2Deg = 180./!dpi
 
  node = at_surveyCenterRa - 90.
  node = node * at_deg2Rad

  etaPole = at_surveyCenterDec * at_deg2Rad

  nra = n_elements(ra_in)
  ndec = n_elements(dec_in)
  IF nra NE ndec THEN BEGIN 
      print,'ra and dec must be same size'
      return
  ENDIF 

  ;; Convert to radians
  ra = double(ra_in*at_deg2Rad)
  dec = double(dec_in*at_deg2Rad)

  x1 = cos(ra-node)*cos(dec)
  y1 = sin(ra-node)*cos(dec)
  z1 = sin(dec)

  slong = -asin( temporary(x1) )             ;lambda
  slat = atan( temporary(z1), temporary(y1) ) - etaPole ;eta

  ;; convert to degrees
  slong = slong * at_rad2Deg
  slat = slat * at_rad2Deg

  atbound2, slong, slat
  atbound, slat, -180.0, 180.0
  
  w=where(slat GT (90. - at_surveyCenterDec) , nw)
  IF nw NE 0 THEN BEGIN 
      slat[w] = slat[w] - 180.
      slong[w] = 180. - slong[w]
  ENDIF 
  atbound, slong, -180.0, 180.0
  
  return
END