;IDL program to create polar plots of GRACE density data mapped
;to 470km using MSIS.

;;;;;read in GRACE density----------------------------------------------
filename = file_search('./GRACE_A_MSIS_470km_2005.dat')
nfiles = n_elements(filename)
print, 'number of files = ',nfiles

for jj = 0, nfiles-1 do begin ;begin major analysis loop, 1 file/loop.

    data = read_grace_msis(filename(jj))
    
    time = data.doy
    alt = data.alt
    lon = data.lon
    lat = data.lat
    stl = data.stl
    mlon = data.mlon
    mlat = data.mlat
    mlt = data.mlt
    den_grace = data.den_grace
    den_msis = data.den_msis
    den_msis_470 = data.den_msis_470
    den_grace_470 = data.den_grace_470
endfor

    set_plot,'ps'
    device,filename='GRACE_den_470km_200501.ps';,/landscape ;/portrait
    device,/color
    
    time_mp11 = 7.39375
    time_mp12 = 8.00347
    time_mp21 = 21.7243
    time_mp22 = 22.2528

    loadct, 39
    ;!p.multi =[0,2,1]
    plot, time, den_grace_470,ystyle = 1,xstyle = 1,ytitle='Den_470km(GRACE)', $
      xr = [6.,9.],yr = [0., 3.5e-15], charsize = 1.5,xticks = 5,thick = 2
    oplot, [time_mp11, time_mp11],[0,10],linestyle = 2

    plot, time, den_grace_470,ystyle = 1,xstyle = 1,ytitle='Den_470km(GRACE)', $
      xr = [20,23],yr = [0., 3.5e-15],charsize = 1.5,xticks = 5,thick=2
    oplot, [time_mp21, time_mp21],[0,10],linestyle = 2
    !p.multi = 0
    device,/close
    set_plot,'x'

stop
    set_plot,'ps'
    device,filename='GRACE_den_470km.ps',/landscape ;/portrait
    device,/color
    
    loadct, 39
    !p.multi =[0,1,2]
    plot, time, den_grace_470,ystyle = 8,xstyle = 1,ytitle='Den_470km(GRACE)', $
      xr = [7.25,8.5],yr = [0., 2e-15], charsize = 1.5,xticks = 5,thick = 2
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 1.5, yticks = 4
    oplot,time,mlat,color = 50

    plot, time, den_grace_470,ystyle = 8,xstyle = 1,ytitle='Den_470km(GRACE)', $
      xr = [21.5,22.75],yr = [0., 4e-15],charsize = 1.5,xticks = 5,thick=2
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 1.5, yticks = 4
    oplot,time,mlat,color = 50

    !p.multi = 0
    device,/close
    set_plot,'x'


end