function orbit_beta,time,alpha=alpha
  if n_elements(time) gt 1 then begin
    result=dblarr(n_elements(time))
    alpha=result
    for i=0,n_elements(time)-1 do begin
      result[i]=orbit_beta(time[i],alpha=this_alpha)
      alpha[i]=this_alpha
    end
    return,result
  end
  rv=scpv(time)
  r=rv[0:2]
  v=rv[3:5]
  ;print,r
  ;print,v
  n=crossp_grid(r,v)
  ;print,n
  sun=sunvec(time)
  ;print,sun
  zhat=n/norm(n)
  sun_z_comp=dotp(zhat,sun)
  beta_angle=asin(sun_z_comp)*!radeg
  sun_z_proj=zhat*sun_z_comp
  sun_z_projperp=sun-sun_z_proj
  xhat=sun_z_projperp/norm(sun_z_projperp)
  yhat=crossp_grid(xhat,zhat)
  xcomp=dotp(r,xhat)
  ycomp=dotp(r,yhat)
  alpha=atan(ycomp,xcomp)*!radeg
  return,beta_angle
end