function sza,log,mtime

; Test a given test log and mission time, calculate the current rotation period
; and solar zenith angle (using Marty Tomasko's formulas)
; For DISR#2.

common already_done,doneit,t,t2,in_lock_time,out_of_lock_time,lock_group,oldlog

if n_elements(oldlog) eq 0 then goto,initialize
if log ne oldlog then goto,initialize
if n_elements(doneit) ne 0 then goto,process

initialize:
; Does the sun pulse file already exist?
fn=findfile(log+'/post/anal_results/sun_pulses')
if fn(0) eq '' then list_sun_pulses,log

str=''
t=dblarr(3,2000)
ttemp=dblarr(3)
openr,unit,log+'/post/anal_results/sun_pulses',/get_lun
readf,unit,str
i=0
while not eof(unit) do begin
  readf,unit,ttemp
  t(*,i)=ttemp
  i=i+1
endwhile
free_lun,unit

t=t(*,0:i-1)
t2=reform(t(1,*))

; Read the message_log to find in-lock and out-of-lock times
in_lock_time=dblarr(100)
out_of_lock_time=dblarr(100)
cmd=string(log+'/message_log',format='("grep using_sun_pulses ",a)')
spawn,cmd,results
n_usings=n_elements(results)
for i=0,n_elements(results)-1 do begin
  reads,results(i),hr,min,sec,format='(56x,f2.0,1x,f2.0,1x,f7.4)'
  in_lock_time(i)=3600.D0*double(hr)+60.D0*double(min)+double(sec)
endfor
in_lock_time=in_lock_time(0:i-1)

cmd=string(log,format='("tail -n1 ",a,"/post/anal_results/sun_pulses")')
spawn,cmd,results
reads,results,o_time,format='(13x,f12.0)'
last_time=o_time+0.01
cmd=string(log+'/message_log',format='("grep SunSensor_sig_lost ",a)')
spawn,cmd,results
n_lost=n_elements(results)
if results(0) ne '' then begin
  k=0
  for j=0,(n_usings<n_lost)-1 do begin
    reads,results(j),hr,min,sec,format='(56x,f2.0,1x,f2.0,1x,f7.4)'
    olt=3600.D0*double(hr)+60.D0*double(min)+double(sec)
    if j eq 0 then begin 
      out_of_lock_time(k)=olt
      k=k+1
    endif
    out_of_lock_time(k)=last_time
    if j gt 0 and olt gt in_lock_time(k) then begin
      out_of_lock_time(k)=olt
      k=k+1
    endif
  endfor
  out_of_lock_time=out_of_lock_time(0:i-1)
endif else begin ; if never lost lock!
  cmd=string(log,format='("tail -n1 ",a,"/post/anal_results/sun_pulses")')
  spawn,cmd,results
  reads,results,o_time,format='(13x,f12.0)'
  out_of_lock_time(0)=o_time+0.01D0
endelse
if out_of_lock_time(i-1) eq 0.0 then begin
  cmd=string(log,format='("tail -n1 ",a,"/post/anal_results/sun_pulses")')
  spawn,cmd,results
  reads,results,o_time,format='(13x,f12.0)'
  out_of_lock_time(i-1)=o_time+0.01D0
endif

; Assign a "lock group" number to each triplet
lock_group=intarr(n_elements(t))-1
for j=0,i-1 do begin
  index=where(t2 ge in_lock_time(j) and t2 le out_of_lock_time(j))
  if index(0) ne -1 then lock_group(index)=j
endfor
doneit=1
oldlog=log

process:
for l=0,n_elements(in_lock_time)-1 do begin
  if mtime ge in_lock_time(l) and mtime le out_of_lock_time(l) then goto,locked
endfor
results=[-1.0,-1.0]
return,results

locked:
index=where(t2 ge mtime)
index=index(0)
if index eq -1 then return,[-1.0,-1.0]
; Case 1 -- three consecutive pulses in lock
if (lock_group(index-1) eq lock_group(index)) and (lock_group(index) eq $
  lock_group(index+1)) then begin
  p1=t2(index)-t2(index-1)
  p2=t2(index+1)-t2(index)
  tt=dblarr(3)
  p=dblarr(3)
  tt(0)=t2(index-1)
  tt(1)=t2(index)
  tt(2)=t2(index+1)
  p(0)=(t2(index-1)-0.5D0*(tt(0)+tt(1)))/(0.5D0*(tt(2)-tt(0)))*(p2-p1)+p1
  p(1)=(t2(index)-0.5D0*(tt(0)+tt(1)))/(0.5D0*(tt(2)-tt(0)))*(p2-p1)+p1
  p(2)=(t2(index+1)-0.5D0*(tt(0)+tt(1)))/(0.5D0*(tt(2)-tt(0)))*(p2-p1)+p1
  deltaz=dblarr(3)
  deltaz(0)=(t(0,index-1)-t(2,index-1))/p(0)*360.D0
  deltaz(1)=(t(0,index)-t(2,index))/p(1)*360.D0
  deltaz(2)=(t(0,index+1)-t(2,index+1))/p(2)*360.D0
  sza=100.91D0+4.2575D0*deltaz+0.11029D0*deltaz^2+0.0014921D0*deltaz^3+ $
    8.0614D-06*deltaz^4
  coefs=poly_fit(tt-tt(0),sza,2)
  x=coefs(0)+coefs(1)*(mtime-tt(0))+coefs(2)*(mtime-tt(0))^2
  p=(mtime-0.5D0*(tt(0)+tt(1)))/(0.5D0*(tt(2)-tt(0)))*(p2-p1)+p1
  results=[x,p]
  return,results
endif

; Case 2 only first two pulses in lock
if (lock_group(index-1) eq lock_group(index)) and (lock_group(index) lt $
  lock_group(index+1)) then begin
  if mtime gt t2(index) then begin
    results=[-1.0,-1.0] & STOP
    return,results
  endif
  p=t2(index)-t2(index-1)
  deltaz=0.5D0*(t(0,index-1)-t(2,index-1)+t(0,index)-t(2,index))
  ; deltaz should now be converted to degrees
  sza=100.91D0+4.2575D0*deltaz+0.11029D0*deltaz^2+0.0014921D0*deltaz^3+ $
    8.0614D-06*deltaz^4
  results=[sza,p]
  return,results
endif

; Case 3 only last two pulses in lock
if (lock_group(index-1) lt lock_group(index)) and (lock_group(index) eq $
  lock_group(index+1)) then begin
  if mtime lt t2(index) then begin
    results=[-1.0,-1.0]
    return,results
  endif
  p=t2(index+1)-t2(index)
  deltaz=0.5D0*(t(0,index)-t(2,index)+t(0,index+1)-t(2,index+1))
  ; deltaz should now be converted to degrees
  sza=100.91D0+4.2575D0*deltaz+0.11029D0*deltaz^2+0.0014921D0*deltaz^3+ $
    8.0614D-06*deltaz^4
  results=[sza,p]
  return,results
endif
return,[-1.0,-1.0]
end