;IDL program to create polar plots of CHAMP density data mapped ;to 360km using MSIS. ;time interval for the main phase time_mp1 = 7.3;217.6;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 = 8.3;218.4;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 = 7.39375;217.796;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 = 2.e-14;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 = 2005; 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 ;;;;;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 ;;;;% of CHAMP orbit reaching polar latitudes >= 75 MLAT dayofyear = findgen(366)+1.5 ; [1.5:1:366.5] per_mlat = dayofyear*!values.f_nan for idoy = 0, n_elements(dayofyear)-1 do begin fp_doy = where((abs(timemp-dayofyear(idoy)) le 0.5), count_doy) fp_mlat = where((abs(timemp-dayofyear(idoy)) le 0.5) and (abs(mlatmp) ge 75.), count_mlat) if (count_doy ne 0.) then per_mlat(idoy) = count_mlat*100./count_doy ;print, count_mlat, count_doy, per_mlat(idoy) endfor ;set_plot, 'ps' ;device, filename='CHAMP_360km_AACGM_orbit_75MLAT'+strtrim(start_end_year(jj),2)+'.ps',/portrait ;device,/color plot, dayofyear, per_mlat, xtitle='DOY',ytitle='percent of CHAMP orbit>75MLAT', $ xstyle=1, yr = [10, 15], psym = 4 ;device,/close ;set_plot,'x' den_sm = smooth(den_champ_360mp,30,/nan);~90 MLAT window(30) 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,2] 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], thick = 2.5 loadct, 18 ;oplot, [0, 366],[abs_cut1, abs_cut1], color = 30 loadct, 39 oplot, [mp1, mp1],[0,1], color = 250, thick = 6 oplot, timemp,den_sm,color = 250, thick = 2.5 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, thick = 2.5 oplot, [mp1, mp1],[-1.,1.], color = 250, 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, thick = 2.5 loadct, 18 oplot, [0, 366], [percent_cut, percent_cut], color = 30 loadct, 39 oplot, [mp1, mp1],[-200,200], color = 250, 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) !p.multi = 0 device,/close set_plot,'x' stop maxden_n = fltarr(npeak)*!values.f_nan deltaden_n = maxden_n percent_n = maxden_n 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 deltaden_s = maxden_n percent_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))) deltaden_n(ii) = alog10(delta_den(fp_peak(ii))) percent_n(ii) = alog10(percent(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))) deltaden_s(ii) = alog10(delta_den(fp_peak(ii))) percent_s(ii) = alog10(percent(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 ;;;------------------- ;;;Write the data into a file ;;;-------------------- ; openw, 22,'CHAMP_360km_AACGM_peaks_30%_'+strtrim(start_end_year(jj),2)+'_180window.dat', width = 800 ; printf, 22, 'time maxden delta_den percent glat glon glt mlat mlon mlt by bz Press SYMH' ; for i=0L, n_elements(timen)-1 do begin ; if (finite(timen(i))) then printf, 22, timen(i), maxden_n(i),deltaden_n(i),percent_n(i),glatn(i),glonn(i),gltn(i), $ ; mlatn(i),mlonn(i),mltn(i),byn(i),bzn(i), Presn(i),symhn(i) ; endfor ; for i=0L, n_elements(times)-1 do begin ; if (finite(times(i))) then printf, 22, times(i), maxden_s(i),deltaden_s(i),percent_s(i),glats(i),glons(i),glts(i), $ ; mlats(i),mlons(i),mlts(i),bys(i),bzs(i), Press(i),symhs(i) ; endfor ; close,22 maxden_n = deltaden_n maxden_s = deltaden_s ;;;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) set_plot,'ps' device,filename='max_CHAMP_deltaden_360km_30%_'+strtrim(start_end_year(jj),2)+'.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,'CHAMP density',charsize=1.1 xyouts,-49.0,40.0,'Northern Hemisphere',charsize=1.1 ;;;;;trajectory-------------------------------------------------------------------------- ;indexn = where((time ge time_mp1) and (time le time_mp2) and (mlat ge 50.),nptsn) ;NH ;if (nptsn gt 0) then begin ; mlatn_orbit = mlat(indexn) ; mltn_orbit = mlt(indexn) ; ; for kk = 0L, nptsn-2 do begin ;Take care of MLT time singularity. ; if (abs(mltn_orbit(kk)-mltn_orbit(kk+1)) gt 12) then begin ; ikk = where(mltn_orbit ge 12) ; mltn_orbit(ikk) = mltn_orbit(ikk)-24. ; endif ; endfor ; ineg = where(mltn_orbit lt 0.,nptneg) ;Make MLT positive ; if(nptneg gt 0) then mltn_orbit(ineg) = mltn_orbit(ineg) + 24. ; ; loadct, 0 ; plots,(90.-mlatn_orbit)*cos((mltn_orbit-6.)*!pi/12.),(90.-mlatn_orbit)*sin((mltn_orbit-6.)*!pi/12.),color = 210,psym = 7, symsize = 0.5 ;endif 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) ;index_finite = where((finite(dcolorn) eq 1) and (byn lt -10.) and (bzn lt 0.0),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(abs_cut2);max(maxden_n,/nan);1.5e-15;(Aug) ;200.;max(binave_Ti); max(alog10(temps),/nan) minData = alog10(abs_cut1);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=1.1,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,'CHAMP density',charsize=1.1 xyouts,-49.0,40.0,'Southern Hemisphere',charsize=1.1 ;;;;;trajectory ;indexs = where((time ge time_mp1) and (time le time_mp2) and (mlat le -50.),nptss) ;SH ;if (nptss gt 0) then begin ; mlats_orbit = mlat(indexs) ; mlts_orbit = mlt(indexs) ; ; for kk = 0L, nptss-2 do begin ;Take care of MLT time singularity. ; if (abs(mlts_orbit(kk)-mlts_orbit(kk+1)) gt 12) then begin ; ikk = where(mlts_orbit ge 12) ; mlts_orbit(ikk) = mlts_orbit(ikk)-24. ; endif ; endfor ; ineg = where(mlts_orbit lt 0.,nptneg) ;Make MLT positive ; if(nptneg gt 0) then mlts_orbit(ineg) = mlts_orbit(ineg) + 24. ; ;loadct, 0 ;plots,(90.+mlats_orbit)*cos((mlts_orbit-6.)*!pi/12.),(90.+mlats_orbit)*sin((mlts_orbit-6.)*!pi/12.),color = 210,psym = 7, symsize = 0.5 ;endif 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) ;index_finite = where((finite(dcolors) eq 1) and (bys lt -10.) and (bzs lt 0.0),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 = alog10(abs_cut2);max(maxden_s,/nan);1.5e-15;(Aug) ;200.;max(binave_Ti); max(alog10(temps),/nan) minData = alog10(abs_cut1);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=1.1,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 device,/close set_plot,'x' endfor ;end major analysis loop. end