;+
; NAME: TMAGZEN
;	
; PURPOSE:
;	To produce an array that masks the pixels that would be found
;       on an image within N degrees of the magnetic zenith at Eureka.
;
; CATEGORY:
;	
; CALLING SEQUENCE:
;       TMAGZEN,YEAR,MONTH,DAY,CAM,FILT,DEG,MASK,CCP,CRP [,SMALL=SMALL] [,HEIGHT=HEIGHT]	
;
; INPUTS:
;       YEAR,MONTH,DAY: The year,month and day the image that the mask will be used on 
;                       was created.
;
;       CAM,FILT: The camera and filter number of the image that the mask will be used on.
;
;       DEG: The number of degrees from the magnetic zenith at Eureka.
;
; KEYWORD PARAMETERS:
;       SMALL: If set, then the mask will be a 256 by 256 BYTE array. Otherwise, a
;              512 by 512 BYTE array will be returned.
;
;       HEIGHT: If set, then the calculations will be performed using the paramters
;               from the magnetic zenith vector at 300 km. Otherwise the magnetic vector
;               at a height of 1 km will be used.
;               
; OUTPUTS:
;       MASK: A BYTE array, that contains 1's were the pixels of an image would be withing
;             DEG degrees of the magnetic zenith and 0's elsewhere.
;
;       CCP: Returns the column number of the pixel that corresponds to the magnetic zenith.
;       
;       CRP: Returns the row number of the pixel that corresponds to the magnetic zenith.
;
; COMMON BLOCKS:
;       None.
;
; SIDE EFFECTS:
;       None.
;
; RESTRICTIONS:
;       None.
;
; PROCEDURE:
;       ....
;
; MODIFICATION HISTORY:
;       Written August 1994, by T A Oliynyk.	
;-

PRO tmagzen,year,month,day,cam,filt,deg,mask,ccp,crp,small=small,height=height

IF KEYWORD_SET(height) THEN BEGIN
  hmag=1515.      ; horizontal magnitude of the zenith magnetic vector at 300 km.
  vmag=49775.     ; vertical magnitude of the zenith magnetic vector at 300 km.
  rot=-83.476     ; rotation of the zenith magnetic vector from north. Positive quantity
		  ; reflects an eastward rotation.
ENDIF ELSE BEGIN
  hmag=2015.      ; same as above except for the magnetic vector at 1 km
  vmag=56351.
  rot=-94.446
ENDELSE

IF KEYWORD_SET(small) THEN edge=256 ELSE edge=512

col=INDGEN(edge)#REPLICATE(1,1,edge)
row=TRANSPOSE(col)

IF NOT KEYWORD_SET(small) THEN small=0
; transform row and column numbers to zenith and azimuth angles
cr2az,year,month,day,cam,filt,col,row,az,za,small=small,site=3 $
     ,/degrees

lambda=(90*!DTOR-ATAN(hmag/vmag))  ; latitude of zenith magnetic vector
phi=(rot+180)*!DTOR                ; east longtitude of zenith magnetic vector

beta=(90-za)*!DTOR  ; zenith angle converted to latitude 
alpha=az*!DTOR      ; azmuth angle is already equal to east longtitude

; transformation produced angles, in degrees, that the points labeled in the 
; zenith array, za, and the azimuth array, az, lie , along great circles, from
; the magnetic zenith vector.
del=ACOS(SIN(lambda)*SIN(beta)+COS(phi-alpha)*COS(lambda)*COS(beta))*!RADEG
mask=BYTARR(edge,edge)
; produce array to located pixels in an image that lie DEG degrees from the magnetic
; zenith.
mask(WHERE(del LE deg))=1B

cpix=MIN(WHERE(del EQ MIN(del)))
crp=cpix/edge
ccp=cpix MOD edge

RETURN
END