;IDL program to create polar plots of CHAMP density data mapped
;to 360km using MSIS.

;time interval for the main phase
time_mp1 = 0.;21.60;7.3
;0.0;217.6;215;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 = 366.;22.40;8.3
;366.;218.4;223;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
mp1 = 0.0;21.7243; 7.39375;plot vertical onset time
;make sure to change the year for read in data file!!!!!!

abs_cut1 = 1.e-15;1.5e-15(Oct);1.e-15(Sep);.6e-15(Aug);1.0e-15(Jan2012)
abs_cut2 = 5.0e-15;2.e-14;3.e-14;1.5e-14;2.e-14;polar plot limits
percent_cut = 30.;30.(Oct);35(Sep);40.(Aug);30.(Jan2012)
plot_cut = 0.0;0.0;2.0e-15;1.5e-15;1.0e-15
Vsc = 7415. ; m/sec 

start_end_year = 2008;indgen(10)+2001;2006;indgen(10)+2001

for jj = 0, n_elements(start_end_year)-1 do begin ;begin major analysis loop, 1 file/loop.
print, start_end_year(jj)
;;;;;read in omni2 minutely By and Bz----------------------------------------------
omni_minute = read_omni_minutely_ByBz('../omni_min_ByBz_'+strtrim(start_end_year(jj),2)+'.lst')
year_omni = omni_minute.year
doy_omni = omni_minute.doy
hour_omni = omni_minute.hour
minute_omni = omni_minute.minute

time_omni = double(doy_omni) + (hour_omni+minute_omni/60.)/24. ;unit:day

imf = omni_minute.imf
by = omni_minute.by
bz = omni_minute.bz
V_sw = omni_minute.V_sw
Pres_sw = omni_minute.Pres_sw
symh = omni_minute.sym_h

;remove bad data points
fp = where(imf ge 9999.99)
imf(fp) = !Values.F_NAN
fp = where(by ge 9999.99)
by(fp) = !Values.F_NAN
fp = where(bz ge 9999.99)
bz(fp) = !Values.F_NAN
fp = where(V_sw ge 99999.9)
V_sw(fp) = !Values.F_NAN
fp = where(Pres_sw ge 99.99)
Pres_sw(fp) = !Values.F_NAN
;fp = where(abs(symh) ge 9999.)
;symh(fp) = !Values.F_NAN

stop

;;;;;read in CHAMP density----------------------------------------------
filename = file_search('./CHAMP_MSIS_400_360_320km_'+strtrim(start_end_year(jj),2)+'.dat')

    data = read_champ_msis(filename)
    
    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_champ = data.den_champ
    den_msis = data.den_msis
    den_champ_400 = data.den_champ_400
    den_champ_360 = data.den_champ_360
    den_champ_320 = data.den_champ_320

;;;;;;;interpol omni miniutely data to CHAMP data time step(30~90sec)--------------
by_int = interpol(by, time_omni, time)
bz_int = interpol(bz, time_omni, time)
Vsw_int = interpol(V_sw, time_omni, time)
Pres_int = interpol(Pres_sw, time_omni, time)
symh_int = interpol(symh, time_omni, time)
        
;;;;;;;;;;;;;;;;;;;;;;;;;;;;CHAMP 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)
        glatmp = lat(indexmp)
        glonmp = lon(indexmp)
        gltmp = stl(indexmp)
        mlatmp = mlat(indexmp)
        mlonmp = mlon(indexmp)
        mltmp = mlt(indexmp)
        den_champ_360mp = den_champ_360(indexmp)
        by_intmp = by_int(indexmp)
        bz_intmp = bz_int(indexmp)
        Vsw_intmp = Vsw_int(indexmp)
        Pres_intmp = Pres_int(indexmp)
        symh_intmp = symh_int(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_champ_360mp,30,/nan)
    delta_den = den_champ_360mp - den_sm
    percent = (den_champ_360mp-den_sm)/den_sm*100
;;;;;;;;finding peaks;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    d0 = den_champ_360mp - shift(den_champ_360mp,1);shift to right
    d1 = den_champ_360mp - shift(den_champ_360mp,-1);shift to left
;    fp_peak = where((d0 gt 0) and (d1 gt 0) and (percent ge percent_cut) and (den_champ_360mp ge abs_cut1),npeak)
    fp_peak = where((d0 gt 0) and (d1 gt 0) and (percent ge percent_cut),npeak)

;;;;;;;;;;calculate hemisphric average without peaks(NH/SH)
    den_withpeak = den_champ_360mp
    den_nopeak = den_withpeak
    den_nopeak(fp_peak) = !values.f_nan

    set_plot,'ps'
    device,filename='CHAMP_den_360km_mp.ps',/landscape ;/portrait
    device,/color
    
    loadct, 39
    !p.multi = [0,1,3]
    
    plot, timemp, den_champ_360mp,ystyle = 8,xstyle = 1,ytitle='Den_360km', $
      xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5,yr=[0,abs_cut2]
    loadct, 18
    ;oplot, [0, 366],[abs_cut1, abs_cut1], color = 30
    loadct, 39
    oplot, [mp1, mp1],[0,1], color = 125, thick = 6
    oplot, timemp,den_sm,color = 250
    oplot, timemp(fp_peak),den_champ_360mp(fp_peak),color = 250, psym = 1, symsize = 1.
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 1.5, yticks = 4
    oplot,timemp,mlatmp,color = 50
    oplot,[0, 365],[0,0],linestyle = 2,color = 50
    ;xyouts, timemp(150), 45., 'threshold = '+string(abs_cut1, format = '(e7.1)')

;;;;;delta density
    plot, timemp, delta_den,ystyle = 8,xstyle = 1,ytitle='delta_Den_360km', $
      yr=[-5.e-15,5.e-15], xr = [time_mp1, time_mp2], charsize = 1.5,xticks = 5
    oplot, [mp1, mp1],[-1.,1.], color = 125, thick = 6
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 1.5, yticks = 4
    oplot,timemp,mlatmp,color = 50
    oplot,[0, 365],[0,0],linestyle = 2,color = 50
    
;;;;;percentage
    plot, timemp, percent, ystyle = 9,xstyle = 1,ytitle='%Den_360km',xtitle = 'DOY', $
      yr = [-150,150], xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5
    loadct, 18
    oplot, [0, 366], [percent_cut, percent_cut], color = 30
    loadct, 39
    oplot, [mp1, mp1],[-200,200], color = 125, thick = 6
    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 = 1.5, 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)
  
;IMF/SW
    plot, timemp, bz_intmp, ystyle = 1,xstyle = 1,ytitle='Bz',xtitle = 'DOY', $
      yr = [-5,5], xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5
    oplot, timemp(fp_peak),bz_intmp(fp_peak),color = 250, psym = 1, symsize = 1.5
    plot, timemp, Vsw_intmp, ystyle = 1,xstyle = 1,ytitle='Vsw',xtitle = 'DOY', $
      yr = [250,350], xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5
    oplot, timemp(fp_peak),Vsw_intmp(fp_peak),color = 250, psym = 1, symsize = 1.5
    plot, timemp, symh_intmp, ystyle = 1,xstyle = 1,ytitle='SYMH',xtitle = 'DOY', $
      yr = [-10,10], xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5
    oplot, timemp(fp_peak),symh_intmp(fp_peak),color = 250, psym = 1, symsize = 1.5

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

    maxden_n = fltarr(npeak)*!values.f_nan
    timen = maxden_n
    glatn  = maxden_n
    glonn  = maxden_n
    gltn  = maxden_n
    mlatn  = maxden_n
    mlonn  = maxden_n
    mltn  = maxden_n
    byn = maxden_n
    bzn = maxden_n
    Presn = maxden_n
    symhn = maxden_n

    maxden_s = maxden_n
    times = maxden_n
    glats  = maxden_n
    glons  = maxden_n
    glts  = maxden_n
    mlats  = maxden_n
    mlons  = maxden_n
    mlts  = maxden_n
    bys = maxden_n
    bzs = maxden_n
    Press = maxden_n
    symhs = maxden_n

    for ii = 0, npeak-1 do begin
        if (mlatmp(fp_peak(ii)) gt 0.0) then begin
            maxden_n(ii) = alog10(den_champ_360mp(fp_peak(ii)))
            timen(ii) = timemp(fp_peak(ii))
            glatn(ii) = glatmp(fp_peak(ii))
            glonn(ii) = glonmp(fp_peak(ii))
            gltn(ii) = gltmp(fp_peak(ii))
            mlatn(ii) = mlatmp(fp_peak(ii))
            mlonn(ii) = mlonmp(fp_peak(ii))
            mltn(ii) = mltmp(fp_peak(ii))
            byn(ii) = by_intmp(fp_peak(ii))
            bzn(ii) = bz_intmp(fp_peak(ii))
            Presn(ii) = Pres_intmp(fp_peak(ii))
            symhn(ii) = symh_intmp(fp_peak(ii))
        endif else begin
            maxden_s(ii) = alog10(den_champ_360mp(fp_peak(ii)))
            times(ii) = timemp(fp_peak(ii))
            glats(ii) = glatmp(fp_peak(ii))
            glons(ii) = glonmp(fp_peak(ii))
            glts(ii) = gltmp(fp_peak(ii))
            mlats(ii) = mlatmp(fp_peak(ii))
            mlons(ii) = mlonmp(fp_peak(ii))
            mlts(ii) = mltmp(fp_peak(ii))
            bys(ii) = by_intmp(fp_peak(ii))
            bzs(ii) = bz_intmp(fp_peak(ii))
            Press(ii) = Pres_intmp(fp_peak(ii))
            symhs(ii) = symh_intmp(fp_peak(ii))
        endelse
    endfor

;;;set the cut for plot in polar contour
    fp_cut = where(maxden_n lt alog10(plot_cut))
    maxden_n(fp_cut) = !Values.F_NaN
    fp_cut = where(maxden_s lt alog10(plot_cut))
    maxden_s(fp_cut) = !Values.F_NaN

;;;;;;;;change coordinate and sort by magnitude;;---------------------
    radiusn = 90.-mlatn
    anglen = (mltn-6.)*!pi/12.
    xposn = radiusn*cos(anglen)
    yposn = radiusn*sin(anglen)
;    dcolorn = (maxden_n-alog10(abs_cut1))/(max(maxden_n,/nan)-alog10(abs_cut1))*256
    dcolorn = (maxden_n-alog10(abs_cut1))/(alog10(abs_cut2)-alog10(abs_cut1))*256
    fp_sortn = sort(dcolorn)
    xposn = xposn(fp_sortn)
    yposn = yposn(fp_sortn)
    dcolorn = dcolorn(fp_sortn)
    byn = byn(fp_sortn)
    bzn = bzn(fp_sortn)
    Presn = Presn(fp_sortn)

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

endfor                          ;end major analysis loop.

end