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

;time interval for the main phase
time_mp1 = 217.6;22.1;269.4;297.6;217.6;69.0479;67.1840;(68.4625;67.000);22.24;297.800;269.526;217.796  ; unit:day
time_mp2 = 218.4;23.6;270.2;298.5;218.4;69.3424;67.3861;(69.3438;67.625);23.22;298.055;269.889;218.144
;make sure to change the year for read in data file!!!!!!

abs_cut = 0.6e-15;1.5e-15(Oct);1.e-15(Sep);.6e-15(Aug);1.0e-15(Jan2012)
percent_cut = 30.;30.(Oct);35(Sep);40.(Aug);30.(Jan2012)

;total_den = fltarr(20000L)*!Values.F_NaN
;total_mlat = fltarr(20000L)*!Values.F_NaN
;total_mlt = fltarr(20000L)*!Values.F_NaN

;;;;;read in GRACE density----------------------------------------------
;filename = file_search('./nonlinfit/GRACE_A_MSIS_470km_2011.dat')
filename = file_search('./GRACE_A_MSIS_470km_2011.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
    
    den_sm = smooth(den_grace_470mp,30,/nan)
    delta_den = den_grace_470mp - den_sm
    percent = (den_grace_470mp-den_sm)/den_sm*100
;;;;;;;;finding peaks;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    d0 = den_grace_470mp - shift(den_grace_470mp,1);shift to right
    d1 = den_grace_470mp - shift(den_grace_470mp,-1);shift to left
    fp_peak = where((d0 gt 0) and (d1 gt 0) and (percent ge percent_cut) and (den_grace_470mp ge abs_cut),npeak)

      

    set_plot,'ps'
    device,filename='GRACE_den_470km_mp.ps',/landscape ;/portrait
    device,/color
    
    loadct, 39
    !p.multi = [0,1,3]
    
    plot, timemp, den_grace_470mp,ystyle = 8,xstyle = 1,ytitle='Den_470km', $
      xr = [time_mp1, time_mp2],charsize = 2.,xticks = 4 ;, psym = -2,symsize = 0.2
    loadct, 18
    oplot, [0, 366],[abs_cut, abs_cut], color = 30
    loadct, 39
    oplot, timemp,den_sm,color = 250
    oplot, timemp(fp_peak),den_grace_470mp(fp_peak),color = 250, psym = 1, symsize = 1.
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 2., yticks = 4
    oplot,timemp,mlatmp,color = 50
    oplot,[0, 365],[0,0],linestyle = 2,color = 50
    xyouts, timemp(150), 45., 'threshold = '+string(abs_cut, format = '(e7.1)')
    
    plot, timemp, delta_den,ystyle = 8,xstyle = 1,ytitle='delta_Den_470km', $
      yr=[-2.e-15,2.e-15], xr = [time_mp1, time_mp2], charsize = 2.,xticks = 4
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 2., yticks = 4
    oplot,timemp,mlatmp,color = 50
    oplot,[0, 365],[0,0],linestyle = 2,color = 50
    
    plot, timemp, percent, ystyle = 8,xstyle = 1,ytitle='%Den_470km',xtitle = 'DOY', $
      yr = [-150,150], xr = [time_mp1, time_mp2],charsize = 2.,xticks = 4
    loadct, 18
    oplot, [0, 366], [percent_cut, percent_cut], color = 30
    loadct, 39
    oplot, timemp(fp_peak),percent(fp_peak),color = 250, psym = 1, symsize = 1.
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 2., yticks = 4
    oplot,timemp,mlatmp,color = 50
    oplot,[0, 365],[0,0],linestyle = 2,color = 50
    xyouts, timemp(150), 45., 'threshold = ' + strmid(percent_cut,6,2)
    
    !p.multi = 0
    device,/close
    set_plot,'x'
;stop    

    maxden_n = fltarr(npeak)*!values.f_nan
    mlatn  = maxden_n
    mltn  = maxden_n
    maxden_s = maxden_n
    mlats  = maxden_n
    mlts  = maxden_n

    for ii = 0, npeak-1 do begin
        if (mlatmp(fp_peak(ii)) gt 0.0) then begin
            maxden_n(ii) = alog10(den_grace_470mp(fp_peak(ii)))
            mlatn(ii) = mlatmp(fp_peak(ii))
            mltn(ii) = mltmp(fp_peak(ii))
        endif else begin
            maxden_s(ii) = alog10(den_grace_470mp(fp_peak(ii)))
            mlats(ii) = mlatmp(fp_peak(ii))
            mlts(ii) = mltmp(fp_peak(ii))
        endelse
    endfor


    radiusn = 90.-mlatn
    anglen = (mltn-6.)*!pi/12.
    xposn = radiusn*cos(anglen)
    yposn = radiusn*sin(anglen)
    dcolorn = (maxden_n-alog10(abs_cut))/(max(maxden_n,/nan)-alog10(abs_cut))*256
    fp_sortn = sort(dcolorn)
    xposn = xposn(fp_sortn)
    yposn = yposn(fp_sortn)
    dcolorn = dcolorn(fp_sortn)

    radiuss = 90.+mlats
    angles = (mlts-6.)*!pi/12.
    xposs = radiuss*cos(angles)
    yposs = radiuss*sin(angles)
    dcolors = (maxden_s-alog10(abs_cut))/(max(maxden_s,/nan)-alog10(abs_cut))*256
    fp_sorts = sort(dcolors)
    xposs = xposs(fp_sorts)
    yposs = yposs(fp_sorts)
    dcolors = dcolors(fp_sorts)

endfor                          ;end major analysis loop.

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

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

indexn = where((time ge time_mp1) and (time le time_mp2) and (mlat ge 50.),nptsn) ;NH
if (nptsn gt 0) then begin    
    mlatn = mlat(indexn)
    mltn = mlt(indexn)
    
    for kk = 0L, nptsn-2 do begin ;Take care of MLT time singularity.
        if (abs(mltn(kk)-mltn(kk+1)) gt 12) then begin
            ikk = where(mltn ge 12)
            mltn(ikk) = mltn(ikk)-24.
        endif
    endfor
    ineg = where(mltn lt 0.,nptneg) ;Make MLT positive
    if(nptneg gt 0) then mltn(ineg) = mltn(ineg) + 24.
endif
;loadct, 33
;plots,(90.-mlatn)*cos((mltn-6.)*!pi/12.),(90.-mlatn)*sin((mltn-6.)*!pi/12.),color = 0,psym = 7, symsize = 0.3
loadct, 0
plots,(90.-mlatn)*cos((mltn-6.)*!pi/12.),(90.-mlatn)*sin((mltn-6.)*!pi/12.),color = 210,psym = 7, symsize = 0.5



loadct, 39
;additional parts
pi = 3.141592654
angle = 2.*pi*indgen(721)/720.
for i=10.,40.,10. do begin
    radius = replicate(i,721)
    oplot,/polar,radius,angle;,linestyle =2
endfor
oplot,[0.,0.],[-45.,45.]
oplot,[-45.,45.],[0.,0.]

for i=10,40,10 do begin
    tang = -45.*pi/180.
    
    txt = ' '+strtrim(string(90-i),2) + '!9%!5' ;northern hemisphere
;   txt = ' '+strtrim(string(i-90),2) + '!9%!5'
    
    if (i eq 40) then txt = txt + ' MLAT'
    xyouts,float(i)*cos(tang),float(i)*sin(tang),txt,charsize=csize
endfor
xyouts,0.,-47.5,'0!CMLT',charsize=csize,alignment=0.5
xyouts,46.,0.,'6',charsize=csize,alignment=0.
xyouts,0.,46.,'12',charsize=csize,alignment=0.5
xyouts,-46.,0.,'18',charsize=csize,alignment=1.

loadct, 33
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 (sqrt(xpos(in)*xpos(in)+ypos(in)*ypos(in)) le 40.) then begin
        plots, xpos(in), ypos(in), psym = 8, symsize = 0.8
    endif
endfor

;;colorbar
number_levels = 33
levels_contour = dblarr(number_levels)
colors_contour = dblarr(number_levels)
color_table_offset =.0
maxData = alog10(5.e-15);max(maxden_n,/nan);1.5e-15;(Aug) ;200.;max(binave_Ti); max(alog10(temps),/nan)
minData = alog10(0.5e-15);abs_cut;1.5e-16;(Aug) ; 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, [[10^levels_contour],[10^levels_contour]],10^levels_contour,indgen(2),charsize=0.8,charthick=1, $
levels=10^levels_contour,c_color=colors_contour,/fill,position=[0.55,0.05,0.85,0.08],/noerase,ystyle=4, $
xr=[10^minData,10^maxData],xstyle=1,xticks = 2,/xlog


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

indexs = where((time ge time_mp1) and (time le time_mp2) and (mlat le -50.),nptss) ;SH
if (nptss gt 0) then begin    
    mlats = mlat(indexs)
    mlts = mlt(indexs)
    
    for kk = 0L, nptss-2 do begin ;Take care of MLT time singularity.
        if (abs(mlts(kk)-mlts(kk+1)) gt 12) then begin
            ikk = where(mlts ge 12)
            mlts(ikk) = mlts(ikk)-24.
        endif
    endfor
    ineg = where(mlts lt 0.,nptneg) ;Make MLT positive
    if(nptneg gt 0) then mlts(ineg) = mlts(ineg) + 24.
endif
;loadct, 33
;plots,(90.+mlats)*cos((mlts-6.)*!pi/12.),(90.+mlats)*sin((mlts-6.)*!pi/12.),color = 0,psym = 7, symsize = 0.3
loadct, 0
plots,(90.+mlats)*cos((mlts-6.)*!pi/12.),(90.+mlats)*sin((mlts-6.)*!pi/12.),color = 210,psym = 7, symsize = 0.5



loadct, 39
;additional parts
pi = 3.141592654
angle = 2.*pi*indgen(721)/720.
for i=10.,40.,10. do begin
    radius = replicate(i,721)
    oplot,/polar,radius,angle;,linestyle =2
endfor
oplot,[0.,0.],[-45.,45.]
oplot,[-45.,45.],[0.,0.]

for i=10,40,10 do begin
    tang = -45.*pi/180.
    
;    txt = ' '+strtrim(string(90-i),2) + '!9%!5' ;northern hemisphere
   txt = ' '+strtrim(string(i-90),2) + '!9%!5'
    
    if (i eq 40) then txt = txt + ' MLAT'
    xyouts,float(i)*cos(tang),float(i)*sin(tang),txt,charsize=csize
endfor
xyouts,0.,-47.5,'0!CMLT',charsize=csize,alignment=0.5
xyouts,46.,0.,'6',charsize=csize,alignment=0.
xyouts,0.,46.,'12',charsize=csize,alignment=0.5
xyouts,-46.,0.,'18',charsize=csize,alignment=1.

loadct, 33
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 (sqrt(xpos(in)*xpos(in)+ypos(in)*ypos(in)) le 40.) then begin
        plots, xpos(in), ypos(in), psym = 8, symsize = 0.8
    endif
endfor

;;colorbar
number_levels = 33
levels_contour = dblarr(number_levels)
colors_contour = dblarr(number_levels)
color_table_offset =.0
maxData = 5.e-15;max(maxden_s,/nan);1.5e-15;(Aug) ;200.;max(binave_Ti); max(alog10(temps),/nan)
minData = 0.5e-15;abs_cut;1.5e-16;(Aug) ; 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 = 2

device,/close
set_plot,'x'

end