;+
; NAME:
;	OCCTIME
;
; PURPOSE:
;	This procedure computes the times of X-ray occultation of a source
;	cosmic or solar during a BATSE orbit.
;
; CATEGORY:
;	CGRO/BATSE
;
; CALLING SEQUENCE:
;	OCCTIME,SPACECRAFT=SC,UTIME=UT,DAYMARKS=DAY,NIGHTMARKS=NIGHT,$
;	S_ANGLE=S_ANGLE,O_ANGLE=O_ANGLE,NTRANS=NTRANS, XYZ_ORBIT=XYZ_ORBIT,$
;	INRADEC=INRADEC, D_PERIOD=D_PERIOD, BAD_DATA=BAD_DATA
;
; INPUTS:
;	SPACECRAFT:  spacecraft vector in 0.25 units, x,y,z in celestial J2000 coords.
;	UTIME:       time of spacecraft vector, anytim units.
;	XYZ_ORBIT:3x3 array describing orbit in inertial coordinates, normalized
;           first 3 elements, x vector giving s/c vector at ut(good(0))
;	    second 3 elements, y vector in orbital plane
;	    third 3 elements, z vector normal to orbital plane
;	D_PERIOD:	orbital period in seconds 
;
; OPTIONAL INPUTS:
;	INRADEC: Input source RA and DEC in decimal degrees, J2000
;	
; OUTPUTS:
;	DAYMARKS: times of start of spacecraft xray day in sec since 1-jan-1979
;	NIGHTMARKS: times of start of spacecraft xray night
;	S_ANGLE:angle formed by dotting earth-sun and earth-sc vectors, 
;	O_ANGLE: occultation angle as a function of time
;	NTRANS: dummy argument, no longer used, input will be set to 0
;	BAD_DATA: set if data is too bad to find occultations
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
; 	VERSION 2, ras, 29-mar-1996
;	Version 3, ras, 3-June-1996, changed minimum o_angle to 100 from 105 degrees
;-
pro occtime,spacecraft=sc,utime=ut,daymarks=day,nightmarks=night,$
	s_angle=s_angle,o_angle=o_angle,ntrans=ntrans, xyz_orbit=xyz_orbit,$
	inradec=inradec, d_period=d_period, bad_data=bad_data

;o_angle - sun_earth_sc angle of day or night transitions, about 110 deg
	o_angle=180.-horizonangle(sqrt(total((sc/4.)^2,2))) 
	a_oangle=avg(o_angle(where(o_angle gt 100.0))) 
;s_angle - angle formed by dotting earth-sun and earth-sc vectors, 
	s_angle=sun_earth_sc(sc,ut) ;(x_pos,y_pos,z_pos,ut) 
;inradec - optional
;	   2 element vector, ra and dec in degrees of a point source instead of the sun
;
;FIND TIMES OF DAY AND NIGHT BOUNDARIES FOR 30 ORBITS CENTERED AROUND THE
;INITIAL TIME PASSED.
;

	bad_data = 0
	ntrans=0
	nut = n_elements(ut)        
	if n_elements(inradec) eq 2 then sut = inradec else sut = solephut(avg(ut))
        sunxyz= sphcart( sut(0), sut(1))
	a1a2= (sunxyz#xyz_orbit)(0:1)
	a1 = a1a2(0)
	a2 = a1a2(1)

;THE EQUATION FOR THE DAY/NIGHT BOUNDARIES IS GIVEN BY
;	A1 * COS( WT ) + A2 * SIN( WT ) = COS( O_ANGLE/!RADEG )
;SET CWT = COS( WT), SUBSTITUTE FOR THE SIN( WT ) GIVES 
;AN EQUATION QUADRATIC IN CWT.
;SOLVE THE QUADRATIC EQUATION FOR THE TANGENT POINTS, CWT1 AND CWT2.

	a = a1^2 + a2^2
	b = (-2)*a1*cos(a_oangle/!radeg)
	c = cos(a_oangle/!radeg)^2 - a2^2
        temp = b^2 - 4. * a * c
        if temp le 0. then begin
           print,'Error calculating day/night boundaries in OCCTIME.'
           bad_data = 1
           goto, getout
        endif
	cwt1 = (sqrt(temp) - b) / (2.*a)
	cwt2 = (0.0 - sqrt(temp) - b) / (2.*a)
	if a2 ne 0.0 then begin 
		swt1 = (cos(a_oangle/!radeg) - a1*cwt1)/a2
		swt2 = (cos(a_oangle/!radeg) - a1*cwt2)/a2 
		wt1  = atan(swt1,cwt1)
		wt2  = atan(swt2,cwt2) 
	endif else begin
		wt1 = acos( cos(a_oangle/!radeg) / a1)
		wt2 = -wt1
	endelse
	
	t1 = ((wt1  / (!pi *2.0) + 1.0) mod 1.0 )* d_period  
	t2 = ((wt2  / (!pi *2.0) + 1.0) mod 1.0 )* d_period   
	
;IS T1 GOING INTO DAY OR NIGHT?

	if a1 * sin( wt1 ) - a2 * cos( wt1) lt 0 then day=1 else day=0
	
 ;GET DAY AND NIGHT TRANSITIONS FOR 31 ORBITS

	t1 = ut(0) + t1 + d_period*(findgen(31)-15)
	t2 = ut(0) + t2 + d_period*(findgen(31)-15)
	if day then begin 
		day  = t1
		night= t2
	endif else begin
		day  = t2
		night= t1
	endelse

;FIT A CURVE TO A SIX MINUTE STRETCH SPANNING ANY TRANSITIONS WITHIN THE
;TIME INTERVAL CHOSEN, WILL ONLY USE ONE DAY AND ONE NIGHT
;
;O_ANGLE - SUN_EARTH_SC ANGLE OF DAY OR NIGHT TRANSITIONS, ABOUT 110 DEG
;S_ANGLE - ANGLE FORMED BY DOTTING EARTH-SUN AND EARTH-SC VECTORS,
;

;FIND TIMES OF DAY AND NIGHT BOUNDARIES
        diff=s_angle-o_angle ;positive is night, negative is day

        wday = where( (day ge ut(0)-1.024) and (day le ut(nut-1)+1.024),nday)
	if nday ge 1 then begin
	        for i=0,nday-1 do begin
		  find_tangent, ut, diff, day(wday(i)), new_day, error=error
		  if ((not error) and (day(wday(i)) ne new_day)) then goto, NEW_DAY
		endfor
	NEW_DAY:	
		if not error then day= new_day + d_period*(findgen(31)-15)
	endif

	wnight = where( (night ge ut(0)-1.024) and $
			(night le ut(nut-1)+1.024),nnight)
	if nnight ge 1 then begin
	        for i=0,nnight-1 do begin
		  find_tangent, ut, diff, night(wnight(i)), new_night, error=error
		  if ((not error) and (night(wnight(i)) ne new_night)) $
			then goto, NEW_NIGHT
		endfor
	NEW_NIGHT:	
		if not error then night = new_night + d_period*(findgen(31)-15)
	endif

getout:
return
end