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

;time interval for the main phase
time_mp1 = 1.;217.796;69.0479;67.1840;(68.4625;67.000);22.24;297.800;269.526;217.796  ; unit:day
time_mp2 = 360.;218.144;69.3424;67.3861;(69.3438;67.625);23.22;298.055;269.889;218.144

;;;;;read in GRACE density----------------------------------------------
filename = file_search('./GRACE_A_MSIS_470km_2012.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

;;;quaility control
        
;;;;;;;;;;;;;;;;;;;;;;;;;;;;GRACE density for mp;;;;;;;;;;;;;;;;;;;;;;;;;
    indexmp = where((time ge time_mp1) and (time le time_mp2),nptsmp) ;mp
    if (nptsmp gt 0) then begin
        timemp = time(indexmp)
        mlatmp = mlat(indexmp)
        mltmp = mlt(indexmp)
        den_grace_470mp = den_grace_470(indexmp)

        for kk = 0L, nptsmp-2 do begin ;Take care of MLT time singularity.
            if (abs(mltmp(kk)-mltmp(kk+1)) gt 12) then begin
                ikk = where(mltmp ge 12)
                mltmp(ikk) = mltmp(ikk)-24.
            endif
        endfor
        ineg = where(mltmp lt 0.,nptneg) ;Make MLT positive
        if(nptneg gt 0) then mltmp(ineg) = mltmp(ineg) + 24.
    endif

set_plot,'ps'
device,filename='GRACE_den_470km_mp_Aug2011.ps',/landscape;/portrait
device,/color

loadct, 39
!p.multi = [0,1,3]

den_sm = smooth(den_grace_470mp,40,/nan)
percent = (den_grace_470mp-den_sm)/den_sm*100

plot, timemp, den_grace_470mp,ystyle = 8,xstyle = 1,ytitle='Den_470km'
oplot, timemp,den_sm,color = 250
axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT'
oplot,timemp,mlatmp,linestyle = 2,color = 50
oplot,[0, 365],[0,0],linestyle = 2,color = 50

plot, timemp, den_grace_470mp-den_sm,ystyle = 8,xstyle = 1,ytitle='delta_Den_470km',yr=[-2.e-15,2.e-15]
axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT'
oplot,timemp,mlatmp,linestyle = 2,color = 50
oplot,[0, 365],[0,0],linestyle = 2,color = 50

plot, timemp, percent, ystyle = 8,xstyle = 1,ytitle='%Den_470km',xtitle = 'DOY',yr = [-100,100]
axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT'
oplot,timemp,mlatmp,linestyle = 2,color = 50
oplot,[0, 365],[0,0],linestyle = 2,color = 50

!p.multi = 0
device,/close
set_plot,'x'
    
    mlat_next = timemp*!values.f_nan
    mlat_next(0:nptsmp-2) = mlatmp(1:nptsmp-1)
    multip = mlatmp*mlat_next
    fp_traj = where((multip le 0.0),ntraj)
    maxden_n = fltarr(ntraj)*!values.f_nan
    mlatn  = maxden_n
    mltn  = maxden_n
    maxden_s = maxden_n
    mlats  = maxden_n
    mlts  = maxden_n
    for ii = 0, ntraj-2 do begin
        temp_mlat = mlatmp(fp_traj(ii):fp_traj(ii+1))
        temp_den = den_grace_470mp(fp_traj(ii):fp_traj(ii+1))
        temp_percent = percent(fp_traj(ii):fp_traj(ii+1))
        temp_mlt = mltmp(fp_traj(ii):fp_traj(ii+1))
        if (mlatmp(fp_traj(ii)+5) gt 0.0) then begin
            fp_n = where(temp_mlat ge 1.)
            temp_mlat_n = temp_mlat(fp_n)
            temp_den_n = temp_den(fp_n)
            temp_percent_n = temp_percent(fp_n)
            temp_mlt_n = temp_mlt(fp_n)
            maxden_n(ii) = max(temp_den_n,/nan,indexn)
            mlatn(ii) = temp_mlat_n(indexn)
            mltn(ii) = temp_mlt_n(indexn)
            if (temp_percent_n(indexn) lt 50.) then maxden_n(ii) = !values.f_nan
        endif else begin
            fp_s = where(temp_mlat le -1.)
            temp_mlat_s = temp_mlat(fp_s)
            temp_den_s = temp_den(fp_s)
            temp_percent_s = temp_percent(fp_s)
            temp_mlt_s = temp_mlt(fp_s)
            maxden_s(ii) = max(temp_den_s,/nan,indexs)
            mlats(ii) = temp_mlat_s(indexs)
            mlts(ii) = temp_mlt_s(indexs)
            if (temp_percent_s(indexs) lt 50.) then maxden_s(ii) = !values.f_nan
        endelse
    endfor

    radiusn = 90.-mlatn
    anglen = (mltn-6.)*!pi/12.
    xposn = radiusn*cos(anglen)
    yposn = radiusn*sin(anglen)
    dcolorn = (maxden_n-1.5e-16)/(1.5e-15-1.5e-16)*256
    radiuss = 90.+mlats
    angles = (mlts-6.)*!pi/12.
    xposs = radiuss*cos(angles)
    yposs = radiuss*sin(angles)
    dcolors = (maxden_s-1.5e-16)/(1.5e-15-1.5e-16)*256

endfor                          ;end major analysis loop.

set_plot,'ps'
device,filename='max_GRACE_den_470km_mp.ps',/portrait
device,/color

loadct, 33
polar_dial_n,1,1
;xyouts,.0,54.0,title1,alignment=0.5,charsize=1.2
xyouts,-49.0,44.0,'density',charsize=1.1
xyouts,-49.0,40.0,'Northern Hemisphere',charsize=1.1

index_finite = where(finite(dcolorn) eq 1,nfinite)
xpos = xposn(index_finite)
ypos = yposn(index_finite)
dcolor = dcolorn(index_finite)
for in = 0, nfinite-1 do begin
    usersym,[-1,1,1,-1,-1],[1,1,-1,-1,1],/fill,color =  dcolor(in);filled square
    if (abs(xpos(in)) le 40.) then begin
        plots, xpos(in), ypos(in), psym = 8, symsize = 1.
    endif
endfor

;;colorbar
number_levels = 33
levels_contour = dblarr(number_levels)
colors_contour = dblarr(number_levels)
color_table_offset =.0
maxData = 1.5e-15 ;200.;max(binave_Ti); max(alog10(temps),/nan)
minData = 1.5e-16; min(alog10(temps),/nan)
for i = 0, number_levels-1 do begin
   levels_contour(i) = (maxData-minData)*i/(number_levels-1.0)+minData
   colors_contour(i) = (255.0-color_table_offset)*i/(number_levels-1.0)+color_table_offset
endfor
contour, [[levels_contour],[levels_contour]],levels_contour,indgen(2),charsize=0.8,charthick=1, $
levels=levels_contour,c_color=colors_contour,/fill,position=[0.55,0.05,0.85,0.08],/noerase,ystyle=4, $
xr=[minData,maxData],xstyle=1,xticks = 3


polar_dial_s,1,2
;xyouts,.0,54.0,title1,alignment=0.5,charsize=1.2
xyouts,-49.0,44.0,'density',charsize=1.1
xyouts,-49.0,40.0,'Southern Hemisphere',charsize=1.1

index_finite = where(finite(dcolors) eq 1,nfinite)
xpos = xposs(index_finite)
ypos = yposs(index_finite)
dcolor = dcolors(index_finite)
for in = 0, nfinite-1 do begin
    usersym,[-1,1,1,-1,-1],[1,1,-1,-1,1],/fill,color =  dcolor(in);filled square
    if (abs(xpos(in)) le 40.) then begin
        plots, xpos(in), ypos(in), psym = 8, symsize = 1.
    endif
endfor

;;colorbar
number_levels = 33
levels_contour = dblarr(number_levels)
colors_contour = dblarr(number_levels)
color_table_offset =.0
maxData = 1.5e-15 ;200.;max(binave_Ti); max(alog10(temps),/nan)
minData = 1.5e-16; min(alog10(temps),/nan)
for i = 0, number_levels-1 do begin
   levels_contour(i) = (maxData-minData)*i/(number_levels-1.0)+minData
   colors_contour(i) = (255.0-color_table_offset)*i/(number_levels-1.0)+color_table_offset
endfor
contour, [[levels_contour],[levels_contour]],levels_contour,indgen(2),charsize=0.8,charthick=1, $
levels=levels_contour,c_color=colors_contour,/fill,position=[0.55,0.05,0.85,0.08],/noerase,ystyle=4, $
xr=[minData,maxData],xstyle=1,xticks = 3

device,/close
set_plot,'x'

end