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