;If there is an asc node within [t_start,t_start+1] then that is the first asc node,
;Assumes orbit period is between 5000 and 10000 seconds (valid for AIM)
;otherwise first asc node is first asc node before t_start
pro find_nodes,t_start,t_stop,a_node=a_node,d_node=d_node,n_term=n_term,s_term=s_term,sunrise=sunrise,sunset=sunset,rise150=rise150,set150=set150
  if t_start lt 1d9 then message,'Called with GPS seconds. Please change calling routine to pass GPS usec.'
  if t_stop  lt 1d9 then message,'Called with GPS seconds. Please change calling routine to pass GPS usec.'
  a_node=dblarr(16)
  a_node_count=0
  d_node=dblarr(16)
  d_node_count=0
  time_ofs=t_start
  t_lo=0d6
  t_hi=1d6
  t_max=t_stop-t_start
  f_lo=find_node_f(t_lo,time_ofs=time_ofs)
  f_hi=find_node_f(t_hi,time_ofs=time_ofs)
  while(f_hi*f_lo gt 0) do begin
    t_hi=t_lo+(t_hi-t_lo)*2
    f_hi=find_node_f(t_hi,time_ofs=time_ofs)
  end
  while (t_hi lt t_max) do find_node,time_ofs,t_lo,t_hi,a_node,a_node_count,d_node,d_node_count
  ;Three times more to make sure we include the orbit with the end time in it.
  find_single_node,time_ofs,t_lo,t_hi,a_node,a_node_count,d_node,d_node_count
  find_single_node,time_ofs,t_lo,t_hi,a_node,a_node_count,d_node,d_node_count
  find_single_node,time_ofs,t_lo,t_hi,a_node,a_node_count,d_node,d_node_count
  a_node=a_node[0:a_node_count-1]
  d_node=d_node[0:d_node_count-1]
  ;Orbits start on an ascending node
  if(d_node[0] lt a_node[0]) then d_node=d_node[1:n_elements(d_node)-1]

  ;Now, find all the terminator crossings
  n_term=dblarr(n_elements(a_node<d_node))
  s_term=dblarr(n_elements(a_node<d_node))
  for i=0,n_elements(n_term)-1 do begin
    n_term[i]=rtbis(0,d_node[i]-a_node[i],!dpi/2,"term_f",0.01,state=state,time_ofs=a_node[i])+a_node[i]
    s_term[i]=rtbis(0,-(d_node[i]-a_node[i]),!dpi/2,"term_f",0.01,state=state,time_ofs=d_node[i])+d_node[i]
  end
  ;And now for the sunrises
  sunrise=dblarr(n_elements(n_term))
  sunset= dblarr(n_elements(s_term))
  rise150=dblarr(n_elements(sunrise))
  set150= dblarr(n_elements(sunset))
  for i=0,n_elements(sunrise)-1 do begin
    sunrise[i]=sofie_time_for_alt(n_term[i]-20d*60d6,n_term[i],  0)
    rise150[i]=sofie_time_for_alt(n_term[i]-20d*60d6,n_term[i],150)
    sunset[i] =sofie_time_for_alt(s_term[i],s_term[i]+20d*60d6,  0)
    set150[i] =sofie_time_for_alt(s_term[i],s_term[i]+20d*60d6,150)
  end

end