;-----------------------------------------------------------------------------
; figure_10a.pro:
; Reproduces the results shown in Figure 10a of Cranmer & Saar (2011).
;
; Calls:   IDL procedure boreas.pro
; Inputs:  ASCII file titled 'mdot_obs_table2.dat' (in current directory)
; Outputs: Plot to current device
;-----------------------------------------------------------------------------

; First read in empirical data for database of 47 stars.

    FILE1 = 'mdot_obs_table2.dat'

    get_lun,unit & openr,unit,FILE1 & readf,unit,npts
    Mstar    = dblarr(npts)
    Rstar    = dblarr(npts)
    Lstar    = dblarr(npts)
    Mdot_obs = dblarr(npts)
    Protday  = dblarr(npts)
    FeH      = dblarr(npts)
    for i=0,npts-1 do begin
      readf,unit,x1,x2,x3,x4,x5,x6
      Mstar(i)    = x1
      Rstar(i)    = x2
      Lstar(i)    = x3
      Mdot_obs(i) = x4
      Protday(i)  = x5
      FeH(i)      = x6
    endfor
    close,unit & free_lun,unit

; For each star, call BOREAS to compute the theoretical mass loss rate.

    Mdot_theory = dblarr(npts)

    for i=0,npts-1 do begin
      BOREAS,Mstar(i),Rstar(i),Lstar(i),Protday(i),FeH(i), $
        Rossby,fstar,Mdot,Mdot_hot,Mdot_cold,MATR
      Mdot_theory(i) = Mdot
    endfor

; Make a simple plot similar to Figure 10a of Cranmer & Saar (2011).

    plot,alog10(Lstar),alog10(Mdot_obs),psym=7, symsize=1.3, $
      xstyle=1,ystyle=1,xrange=[-2.2,4.5],yrange=[-19.,-3.]

    aa = findgen(16) * (!PI * 2 / 15.)
    usersym,cos(aa),sin(aa)
    oplot,alog10(Lstar),alog10(Mdot_theory),psym=8, symsize=1.4

    end