;-----------------------------------------------------------------------------
; analyze_fig4.pro:
; Reproduces the results shown in Figure 4 of Cranmer et al (2014).
;
; Inputs:  ASCII file 'cranmer_granule13_kepler.dat' (in current directory)
; Outputs: Simplified plot to current device (no color)
;-----------------------------------------------------------------------------

; read data file

    FILE1 = 'cranmer_granule13_kepler.dat'

    get_lun,unit & openr,unit,FILE1 & readf,unit,npts
    KIC         = dblarr(npts)
    Teff        = dblarr(npts)
    logg        = dblarr(npts)
    Mstar       = dblarr(npts)
    flag_mass   = intarr(npts)
    tau_eff     = dblarr(npts)
    Mach        = dblarr(npts)
    sigma       = dblarr(npts)
    F8ppt_model = dblarr(npts)
    F8ppt_obs   = dblarr(npts)
    Rppt_obs    = dblarr(npts)
    ZC_obs      = dblarr(npts)
    for i=0,npts-1 do begin
      readf,unit,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12
      KIC(i)         = x1
      Teff(i)        = x2
      logg(i)        = x3
      Mstar(i)       = x4
      flag_mass(i)   = x5
      tau_eff(i)     = x6
      Mach(i)        = x7
      sigma(i)       = x8
      F8ppt_model(i) = x9
      F8ppt_obs(i)   = x10
      Rppt_obs(i)    = x11
      ZC_obs(i)      = x12
    endfor
    close,unit & free_lun,unit

; reproduce a simplified version of Figure 4 (on default plotting device)

    !p.multi=[0,1,2]

    xr = [0., 0.355]
    yr = [4.8, 2.25]

;(a)

    plot,xr,yr, /nodata, $
     xsty=1,ysty=1,xrange=xr,yrange=yr

    npoly    = 200
    fmin     = 0.014
    fmax     = 0.34
    eee      = 0.2
    f8poly   = dindgen(npoly)/double(npoly-1)*(fmax-fmin)+fmin
    apoly    = alog10(f8poly)
    loggpoly = 1.15136 - 3.59637*apoly - 1.40002*apoly^2 - 0.22993*apoly^3
    xpoly    = [f8poly,reverse(f8poly),f8poly(0)]
    ypoly    = [loggpoly-eee,reverse(loggpoly)+eee,loggpoly(0)-eee]
    oplot,xpoly,ypoly

    oplot,F8ppt_obs,logg,psym=4,symsize=0.5

;(c)

    plot,xr,yr, /nodata, $
     xsty=1,ysty=1,xrange=xr,yrange=yr

    oplot,xpoly,ypoly

    oplot,F8ppt_model,logg,psym=4,symsize=0.5

    end