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

; First read in theoretical 2D grid of ZAMS models.

    FILE1 = 'zams_girardi_grid.dat'

    get_lun,unit & openr,unit,FILE1 & readf,unit,nteff,nrot
    Mstar    = dblarr(nteff)
    Rstar    = dblarr(nteff)
    Lstar    = dblarr(nteff)
    Protday  = dblarr(nrot)
    for i=0,nteff-1 do begin
     for j=0,nrot-1 do begin
       readf,unit,x1,x2,x3,x4
       Mstar(i)    = x1
       Rstar(i)    = x2
       Lstar(i)    = x3
       Protday(j)  = x4
     endfor
    endfor
    close,unit & free_lun,unit

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

    Mdot_theory = dblarr(nteff,nrot)
    FeHnow      = 0.0

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

; Compute effective temperature for use as x-axis coordinate.

    xRsun  = 6.96d+10
    xLsun  = 3.826d+33
    stefan = 5.66961d-5
    xRcgs  = Rstar * xRsun
    xLcgs  = Lstar * xLsun
    Teff   = (xLcgs/4./!PI/xRcgs/xRcgs/stefan)^0.25

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

    xr = [3200., 9400.]
    yr = [-21.5, -10.]
    plot,xr,yr, /nodata, $
      xstyle=1,ystyle=1,xrange=xr,yrange=yr

    for j=0,nrot-1 do begin
      oplot,Teff,alog10(Mdot_theory(*,j))
    endfor

    aa = findgen(16) * (!PI * 2 / 15.)
    usersym,cos(aa),sin(aa)
    oplot,[5770.],[alog10(2.5d-14)], psym=8,symsize=1.3

    end