pro interpolate_j,j_call,j_smootha ;routine to fill data gaps by interpolation, can be call iteratively ;can be used for e-,ions,energy flux, number flux nmlt = 96 nmlat = 160 ;since we are interpolating in MLT at fixed latitude, latitude is the ;outer loop (most other precip model code has lat as inner loop) ;function mlt_plus(mlt_current,delta_t) ;(delta_t can be positive or negative) j_smootha = fltarr(nmlt,nmlat) jea_new = j_call for j = 0, nmlat -1 do begin for i = 0, nmlt-1 do begin if( jea_new(i,j) gt 0. )then begin j_smootha(i,j) = jea_new(i,j) goto, next_mlt endif ;find closest neighbor to the right t_plus = 0 t_plus_bin = i while( (t_plus lt 8) and (jea_new(t_plus_bin,j) eq 0) )do begin t_plus = t_plus + 1 t_plus_bin = mlt_plus(i,t_plus) endwhile ;find earlier mlt neighbor (to the left) t_minus = 0 t_minus_bin = i while( (t_minus gt -8) and (jea_new(t_minus_bin,j) eq 0.) )do begin t_minus = t_minus - 1 t_minus_bin = mlt_plus(i,t_minus) endwhile if( (t_plus ge 8) or (t_minus le -8) )then goto, next_mlt if( t_plus eq 0 )then begin print,'Error: t_plus = 0 ',t_plus,i,j,jea_new(i,j),t_plus_bin endif if( (t_plus lt 8) and (t_minus gt -8) )then begin p1 = t_plus_bin p2 = mlt_plus(t_plus_bin,1) p3 = mlt_plus(t_plus_bin,2) jep1 = jea_new(p1,j) jep2 = jea_new(p2,j) jep3 = jea_new(p3,j) n1 = t_minus_bin n2 = mlt_plus(t_minus_bin,-1) n3 = mlt_plus(t_minus_bin,-2) jen1 = jea_new(n1,j) jen2 = jea_new(n2,j) jen3 = jea_new(n3,j) wp1 = 1./(t_plus*t_plus) wp2 = 1./((t_plus+1)*(t_plus+1)) wp3 = 1./((t_plus+2)*(t_plus+2)) wn1 = 1./(t_minus*t_minus) wn2 = 1./((t_minus-1)*(t_minus-1)) wn3 = 1./((t_minus-2)*(t_minus-2)) wsum = wp1 + wp2 + wp3 + wn1 + wn2 + wn3 j_smooth = jep1*wp1 + jep2*wp2 + jep3*wp3 + wn1*jen1 + wn2*jen2 + wn3*jen3 ;norm = sqrt(wp1*wp1 + wp2*wp2 + wp3*wp3 + wn1*wn1 + wn2*wn2 + wn3*wn3) ;norm = norm*sqrt(6.0) norm = wsum j_smooth = j_smooth/norm j_smootha(i,j) = j_smooth endif next_mlt: dummy = 0 endfor ;i loop (MLT) endfor ;j loop (MLAT) return end pro smooth_j,j_call,j_smootha,n3,n3_new ;routine to do a nearest neighbor smoothing, can be called iteratively ;can be used for e-,ions,energy flux, number flux ;common total_obs,n3 nmlt = 96 nmlat = 160 n3_new = n3 j_smootha = j_call for i = 0, nmlt-1 do begin ;for j = 1, nmlat -2 do begin ;for j=2, nmlat - 3 do begin for j = 0,nmlat-1 do begin ;if( j eq ((nmlat/2) - 2) )then j=j+4 ;if (j eq ((nmlat/2) - 3) )then j=j+6 jplus1 = mlat_plus(j,1) jplus2 = mlat_plus(j,2) jminus1 = mlat_plus(j,-1) jminus2 = mlat_plus(j,-2) p1 = mlt_plus(i,1) p2 = mlt_plus(i,2) je0 = j_call(i,j) jep1 = j_call(p1,j) jep2 = j_call(p2,j) n30 = n3(i,j) n3p1 = n3(p1,j) n3p2 = n3(p2,j) m1 = mlt_plus(i,-1) m2 = mlt_plus(i,-2) jem1 = j_call(m1,j) jem2 = j_call(m2,j) n3m1 = n3(m1,j) n3m2 = n3(m2,j) ja = j_call(i,jplus1) jaa = j_call(i,jplus2) jb = j_call(i,jminus1) jbb = j_call(i,jminus2) n3a = n3(i,jplus1) n3b = n3(i,jminus1) n3aa = n3(i,jplus2) n3bb = n3(i,jminus2) w0 = 3. wp1 = 1./2. wp2 = 1./8. wm1 = 1./2. wm2 = 1./8. wa = 1./3. wb = 1./3. waa = 1./27. wbb = 1./27. if( n30 lt 3 )then w0 = 0. if( n3p1 lt 3 )then wp1 = 0. if( n3p2 lt 3 )then wp2 = 0. if( n3m1 lt 3 )then wm1 = 0. if( n3m2 lt 3 )then wm2 = 0. if( n3a lt 3 )then wa = 0. if( n3b lt 3 )then wb = 0. if( n3aa lt 3)then waa = 0. if( n3bb lt 3)then wbb = 0. wsum = w0 + wp1 + wp2 + wm1 + wm2 + wa + wb + waa + wbb if( wsum gt 0. )then begin j_smooth = w0*je0 + wp1*jep1 + wp2*jep2 + wm1*jem1 + wm2*jem2 + $ wa*ja + wb*jb + waa*jaa + wbb*jbb n_smooth = w0*n30 + wp1*n3p1 + wp2*n3p2 + wm1*n3m1 + wm2*n3m2 + $ wa*n3a + wb*n3b + waa*n3aa + wbb*n3bb norm = wsum j_smooth = j_smooth/norm n_smooth = fix(0.5 + n_smooth/norm) j_smootha(i,j) = j_smooth n3_new(i,j) = n_smooth endif next_mlt: dummy = 0 endfor ;i loop (MLT) endfor ;j loop (MLAT) return end pro draw_je,out_file,jea_erg,n_or_s,je_minc,je_maxc,je_title,je_type,dS,sf, $ NO_SHOW_STATS_YEAR_RANGE=NO_SHOW_STATS_YEAR_RANGE, FORECAST_TIME=FORECAST_TIME, FORECAST_STD=FORECAST_STD, HEMISPHERIC_POWER=power_hemi ;Output goes directly to IDL file out_file. Input from calling routine ;Patrick Newell, September 2008 ;updated variously ;updated July 2009 to add difference between hemispheres ;n_or_s = 1 northern hemisphere ;n_or_s = 2 southern ;n_or_s = 3 both north and south ;;n_or_s= 4 difference between hemispheres ;jea_erg = fltarr(nmlt,nmlat) (no auroral type dependence) ;je_title = title used for plot ;je_type: 1 if actually je-electron,2 for je-ions, 3-je number flux, ;4-ji number flux, 5-eave e- average enegy, 6-iave ion average energy ;7-correlation coefficient, electrons, 8-correlation coefficient ions ;year0,doy0 = start, year1,doy1=end (for labeling plot only) ;dS = Earth surface integral of jea_erg ;sf = seasonal flag.0=winter,1=spring,2=summer,3=fall,4=all,none=all common data_interval,y0,d0,yend,dend,files_done nmlt = 96 nmlat = 160 je_min = je_minc je_max = je_maxc year0 = y0 doy0 = d0 year1 = yend doy1 = dend jea_new = fltarr(nmlt,nmlat/2) if( n_elements(sf) le 0 )then sf = 4 ;default is all seasons power_mlt = fltarr(24) ;E or n flux integ ovr lat set_plot,'ps' device,/color,bits=8 device, filename = out_file dot_spot = strpos(out_file,'.') text_file = strmid(out_file,0,dot_spot+1) text_file = text_file + 'txt' openw,12,text_file printf,12,text_file je_sum = fltarr(nmlt) out1 = '(f6.2,1x,f6.1,1x,f6.2)' if( (je_type eq 3) or (je_type eq 4) )then begin out1 = '(f6.2,1x,f6.1,1x,e)' endif out2 = '(i4,2x,i3)' r = 40. ANGLE = FINDGEN(361) THETA = (2.*!PI*ANGLE)/360. x = r*cos(theta) y = r*sin(theta) ; open color bar, read, and load ;(establishes the link between colors and numbers) red = intarr(256) green = intarr(256) blue = intarr(256) openr,2,'data/wondr2.col' idummy = 1 ry = 1 gy = 1 by = 1 for i=0,255 do begin readf,2,idummy,ry,by,gy red(i)=ry green(i)=gy blue(i)=by endfor close,2 tvlct,red,green,blue erase,0 mlat = 0. mlt = 0. t = 0. r = 0. ;delta_r = .16667 ;delta_t = !pi/24. if( n_or_s eq 4 )then begin ; je_min = -je_maxc/3. ; je_max = je_maxc/3. je_min = -je_maxc/2. je_max = je_maxc/2. endif efact = 250./(je_max-je_min) power_hemi = 0. ;power if je_type=1-2, 3-4=nflux device,/inches,xsize=8.0,ysize=8.0,yoffset=2.0,xoffset=0.5 plot,x,y,xsty=4,ysty=4,thick=2.0,xrange=[-1.0,1.0],yrange=[-1.0,1.0], $ position=[0.2, 0.2, 0.8, 0.8],color=254 for i=0,nmlt-1 do begin mlt = 0.001 + 0.25*float(i) mlt_real = mlt mlt = mlt - 6.0 ; rotate noon to top if( mlt lt 0.0 ) then mlt=mlt + 24.0 if( mlt ge 24.0 )then begin print,'odd mlt ',mlt,' i,j = ',i,j endif for j=0,(nmlat/2)-1 do begin mlat = 50.001 + 0.5*float(j) mlat_real = mlat if( n_or_s eq 2 )then mlat_real = -mlat mlat1 = mlat mlat2 = mlat + 0.5 ibin = mlt_bin(mlt_real) jbin = mlat_bin(mlat_real) jeave = jea_erg(ibin,jbin) if( n_or_s eq 3 )then begin jbin2 = mlat_bin(-mlat_real) je1 = jeave je2 = jea_erg(ibin,jbin2) jeave = (je1 + je2)/2. if( je1*je2 eq 0. )then jeave = je1 + je2 endif if( n_or_s eq 4 )then begin jbin2 = mlat_bin(-mlat_real) je1 = jeave je2 = jea_erg(ibin,jbin2) jeave = 0. if( je1*je2 gt 0.)then jeave = (je1 - je2) ;; if( je2 gt (2.*je1) )then jeave = 0. ;; if( je1 gt (2.*je2) )then jeave = 0. endif ang1 = ((mlt)/24.0) * (2.0*!pi) ang2 = ((mlt+0.25)/24.0) * (2.0*!pi) r1=abs(mlat1-90.)/40. r2=abs(mlat2-90.)/40. x1 = cos(ang1)*r1 y1 = sin(ang1)*r1 x2 = cos(ang1)*r2 y2 = sin(ang1)*r2 x3 = cos(ang2)*r2 y3 = sin(ang2)*r2 x4 = cos(ang2)*r1 y4 = sin(ang2)*r1 fill_col = efact*(jeave - je_min) if( jeave eq 0. )then fill_col = 0 if( fill_col gt 250.) then fill_col = 250 if( fill_col lt 0.) then fill_col = 0 polyfill,[x1,x2,x3,x4],[y1,y2,y3,y4],col=fill_col printf,12,format=out1,mlt_real,mlat_real,jeave je_sum(i) = je_sum(i) + jeave if( (mlt_real ge 7.) and (mlt_real le 19.) and $ (abs(mlat_real) le 63.) )then begin jeave = 0. endif mlt1 = mlt_real mlt2 = mlt1 + 0.25 A_earth = area_sphere(mlt1,mlt2,mlat1,mlat2) ;1 erg = 10**-7 joule but converting cm2 to km2 gives 10**10 ;converting watts to GW is a factor of 10-9 if( n_or_s ne 4 )then begin if( (je_type eq 1) or (je_type eq 2) )then begin power_hemi = power_hemi + (1.0e-6)*A_earth*jeave ;in GW imlt = fix(mlt_real) power_mlt(imlt) = power_mlt(imlt) + (1.0e-6)*A_earth*jeave endif if( (je_type eq 3) or (je_type eq 4) )then begin power_hemi = power_hemi + (1.0e10)*A_earth*jeave imlt = fix(mlt_real) power_mlt(imlt) = power_mlt(imlt) + (1.0e10)*A_earth*jeave endif endif endfor endfor dS = power_hemi ; ; draw mlat circles for r=0.125, 1.0, 0.125 do begin x = r*COS(THETA) y = r*SIN(THETA) oplot, X,Y,thick=2.0,color=255 endfor mlt = 1.2 - 6.0 ; rotate noon to top ang = mlt/24.0 * 2.0*!pi for mlat=50,80,5 do begin ; label mlat ticks r=abs(mlat-89.0)/40.0 x = cos(ang)*r y = sin(ang)*r xyouts,x,y,string(mlat,format='(i2.2)'),color=255 endfor ; draw mlt lines and label delta_theta = 2*!pi/(24.) t = 0. for i = 0,23 do begin t = i*delta_theta - !pi/2.0 if ((i mod 3) eq 0) then begin x= cos(t) y= sin(t) oplot,[0.0,x],[0.0,y],thick=2.0,color=255 x= 1.05*cos(t) y= 1.05*sin(t) xyouts,x,y,string(i,format='(i2.2)'),align=.5, size=1.5 endif else begin x = cos(t)*1.02 y = sin(t)*1.02 oplot,[0.0,x],[0.0,y],thick=2.0,color=255 endelse endfor ;... title map_name = je_title ;xyouts,!d.x_ch_size*2,!d.y_size-4*!d.y_ch_size,size=2,/dev, $ ; map_name black = 252B white = 251B colmax = 250B colmin = 0B maxv = je_max minv = je_min ;if( (je_type eq 3) or (je_type eq 4) )then begin ; maxv = alog10(maxv) ; if( je_min gt 0. )then minv = alog10(je_min) ;endif cb_h = 0.25*!d.y_vsize cb_w = float(!d.x_ch_size) cb = bytarr(cb_w,cb_h) tmp = colmin + findgen(colmax-colmin+1) cbtmp = interpol(tmp,cb_h) for i=0,cb_w-1 do cb(i,*) = cbtmp !P.POSITION = [0.05,0.4,0.05+float(cb_w-1)/!d.x_vsize,0.65] ;oplot,[minv,maxv], /nodata,color=black,xstyle=5, ystyle=5 tv, cb,0.05,0.4, /normal, xsize=cb_w/!d.x_vsize,ysize=cb_h/!d.y_vsize if( je_type le 2 )then begin xyouts,0.07,0.4, /normal, string(minv,format='(f5.2)') xyouts,0.07,0.65, /normal, string(maxv,format='(f5.2)') xyouts,0.07,0.525, /normal, string((maxv+minv)/2.,format='(f5.2)') if( n_or_s ne 4)then begin xyouts,0.05,0.70, /normal, size=2, 'ergs/cm2s' endif endif if( (je_type eq 3) or (je_type eq 4) )then begin xyouts,0.07,0.4, /normal, string(minv,format='(e7.1)') xyouts,0.07,0.65, /normal, string(maxv,format='(e7.1)') xyouts,0.07,0.525, /normal, string((maxv+minv)/2.,format='(e7.1)') if( n_or_s ne 4 )then begin xyouts,0.07,0.70, /normal, '1/cm2 s' endif endif if( (je_type eq 5) or (je_type eq 6) )then begin xyouts,0.07,0.4, /normal, string(minv,format='(f6.0)') xyouts,0.07,0.65, /normal, string(maxv,format='(f6.0)') xyouts,0.07,0.525, /normal, string((maxv+minv)/2.,format='(f6.0)') if( n_or_s ne 4 )then begin xyouts,0.05,0.70, /normal, 'eV' endif endif if( (je_type le 4) and (n_or_s ne 4) )then begin printf,12,'Totals:' grand_total = 0. for i=0,nmlt-1 do begin mlt = float(i)*0.25 printf,12,mlt,je_sum(i) grand_total = grand_total + je_sum(i) endfor printf,12,'Grand total: ',grand_total if(je_type le 2 )then printf,12,'Hemispherical power (GW): ',power_hemi ;endfor printf,12,'Flux integrated over latitude (fixed 1-h MLT):' for i=0,23 do begin printf,12,i,power_mlt(i) endfor power_day = total(power_mlt[6:17]) power_night = total(power_mlt[0:5]) + total(power_mlt[18:23]) printf,12,'Total day = ',power_day printf,12,'Total night= ',power_night endif close,12 ;xyouts,!d.x_ch_size*2,!d.y_size-4*!d.y_ch_size,size=2,/dev, $ ; map_name xyouts,0.07,0.92,/normal, size=2,map_name if ( n_elements( FORECAST_TIME ) eq 1 ) then begin xyouts, 0.07, 0.88, /normal, size=2, 'Forecast '+string( forecast_time, format='(I03)' ) + $ ' minutes (+/- ' + string( forecast_std, format='(F4.1)' ) + '), ' endif if( n_or_s ne 4 )then begin if( je_type le 2 )then begin xyouts,0.70,0.88, /normal, size=2, string(power_hemi,format='(f6.1)') xyouts,0.85,0.88, /normal, size=2, 'GW' endif if( (je_type eq 3) or (je_type eq 4) )then begin xyouts,0.63,0.88, /normal, size=2, string(power_hemi,format='(e8.1)') xyouts,0.80,0.88, /normal, size=2, '/s' endif endif season_string = ['winter','spring','summer','fall',' ','n. wint',$ 'n. spring','n. summer','n. fall'] ; FIXME: correct this next line. season_epoch.pro calls with an incorrect value of always 3, i.e. 'fall' ;xyouts,0.16,0.82, /normal, size=2, season_string(sf) if ( ~ n_elements( NO_SHOW_STATS_YEAR_RANGE ) ) then begin xyouts,0.65,0.82, /normal, size=2, string(year0,format='(i4)') xyouts,0.80,0.82, /normal, size=2, string(doy0,format='(i3)') xyouts,0.65,0.77, /normal, size=2, string(year1,format='(i4)') xyouts,0.80,0.77, /normal, size=2, string(doy1,format='(i4)') endif xyouts,0.05,0.08, /normal, size=2, string(files_done,format='(i5)') xyouts,0.17,0.08, /normal, size=2, 'Satellite-days' device,/close return END pro draw_je_geo,out_file,jea_erg,n_or_s,je_minc,je_maxc,je_title,je_type,dS,sf, $ NO_SHOW_STATS_YEAR_RANGE=NO_SHOW_STATS_YEAR_RANGE, FORECAST_TIME=FORECAST_TIME, FORECAST_STD=FORECAST_STD, HEMISPHERIC_POWER=power_hemi ;Output goes directly to IDL file out_file. Input from calling routine ;Patrick Newell, September 2008 ;updated variously ;updated July 2009 to add difference between hemispheres ;INPUTS ; out_file name of output plot file ; n_or_s = 1 northern hemisphere ; = 2 southern ; = 3 both north and south ; = 4 difference between hemispheres ; jea_erg = fltarr(nmlt,nmlat) (no auroral type dependence) ; je_minc The minimum flux for a given auroral type- used for the plotting color table ; je_maxc The maximum flux for a given auroral type- used for the plotting color table ; je_title = title used for plot ; je_type: 1 if actually je-electron,2 for je-ions, 3-je number flux, ; 4-ji number flux, 5-eave e- average enegy, 6-iave ion average energy ; 7-correlation coefficient, electrons, 8-correlation coefficient ions ; year0,doy0 = start, year1,doy1=end (for labeling plot only) ; OUTPUTS ;dS = Earth surface integral of jea_erg ;sf = seasonal flag.0=winter,1=spring,2=summer,3=fall,4=all,none=all common data_interval,y0,d0,yend,dend,files_done ;nmlt is the number of mltbins nmlt = 96 ;nmlat is the number of magnetic latitude bins nmlat = 160 ;These are to set the min and max for the colorbar je_min = je_minc je_max = je_maxc year0 = y0 doy0 = d0 year1 = yend doy1 = dend ;I'm not sure what jea_new is jea_new = fltarr(nmlt,nmlat/2) if( n_elements(sf) le 0 )then sf = 4 ;default is all seasons power_mlt = fltarr(24) ;E or n flux integ ovr lat ;Now set the plot to be a postscript file set_plot,'PS' ;This sets the plot to color device,/color,bits_per_pixel=8 ;This sets the filename of the output postscript plot file device, filename = out_file ;This is just creating a text file to output the data dot_spot = strpos(out_file,'.') text_file = strmid(out_file,0,dot_spot+1) text_file = text_file + 'txt' openw,12,text_file printf,12,text_file ;I'm not sure what je_sum is je_sum = fltarr(nmlt) ;Set the format for the output text file out1 = '(f6.2,1x,f6.1,1x,f6.2)' if( (je_type eq 3) or (je_type eq 4) )then begin out1 = '(f6.2,1x,f6.1,1x,e)' endif out2 = '(i4,2x,i3)' ;This is the setup to eventually draws an outer circle of radius 40 degrees r = 40. ANGLE = FINDGEN(361) THETA = (2.*!PI*ANGLE)/360. x = r*cos(theta) y = r*sin(theta) ; open color bar, read, and load ;(establishes the link between colors and numbers) red = intarr(256) green = intarr(256) blue = intarr(256) openr,2,'data/wondr2.col' idummy = 1 ry = 1 gy = 1 by = 1 for i=0,255 do begin readf,2,idummy,ry,by,gy red(i)=ry green(i)=gy blue(i)=by endfor close,2 ;This loads the new color table tvlct,red,green,blue erase,0 mlat = 0. mlt = 0. t = 0. r = 0. ;delta_r = .16667 ;delta_t = !pi/24. if( n_or_s eq 4 )then begin ; je_min = -je_maxc/3. ; je_max = je_maxc/3. je_min = -je_maxc/2. je_max = je_maxc/2. endif efact = 250./(je_max-je_min) power_hemi = 0. ;power if je_type=1-2, 3-4=nflux device,/inches,xsize=8.0,ysize=8.0,yoffset=2.0,xoffset=0.5 ;This plots a circle with radius equal to 40 degrees plot,x,y,xsty=4,ysty=4,thick=2.0,xrange=[-1.0,1.0],yrange=[-1.0,1.0], $ position=[0.2, 0.2, 0.8, 0.8],color=254 ;Now loop through all the mltbins for i=0,nmlt-1 do begin ;It's looping through .25 mlt bins mlt = 0.001 + 0.25*float(i) ;This creates a vector mlt_real with the actual mlt of the data ;because in the next few lines he subtracts 6 so that the plot puts ;12 mlt at the top of the page mlt_real = mlt ;mlt=0 would plot to the right at x=0,y=0 ;Subtracting 6 puts mlt=0 at the bottom of the page ;JGREEN: I commented out the mlt-6 for the geographic plot. I add the rotation later so ;that the noon is at the top of the page ;mlt = mlt - 6.0 ; rotate noon to top if( mlt lt 0.0 ) then mlt=mlt + 24.0 if( mlt ge 24.0 )then begin print,'odd mlt ',mlt,' i,j = ',i,j endif for j=0,(nmlat/2)-1 do begin ;Now loop through the mlat bins. ;The mlat starts at 50 ; And it increaases in .5 degree bins mlat = 50.001 + 0.5*float(j) mlat_real = mlat if( n_or_s eq 2 )then mlat_real = -mlat mlat1 = mlat mlat2 = mlat + 0.5 ;Not sure what this is for ibin = mlt_bin(mlt_real) jbin = mlat_bin(mlat_real) jeave = jea_erg(ibin,jbin) if( n_or_s eq 3 )then begin jbin2 = mlat_bin(-mlat_real) je1 = jeave je2 = jea_erg(ibin,jbin2) jeave = (je1 + je2)/2. if( je1*je2 eq 0. )then jeave = je1 + je2 endif if( n_or_s eq 4 )then begin jbin2 = mlat_bin(-mlat_real) je1 = jeave je2 = jea_erg(ibin,jbin2) jeave = 0. if( je1*je2 gt 0.)then jeave = (je1 - je2) ;; if( je2 gt (2.*je1) )then jeave = 0. ;; if( je1 gt (2.*je2) )then jeave = 0. endif ;================================================================= ;JGREEN: This is where the polygons are made to do the plot ;and where the change to geographic coordinates is done ;First define all thes variables ;This is the location of the sun in GSM coordinates sunGSM=[1.0D,0.0D,0.0D] ;This is the location of the sun in cartesion GEO coordinates sunGEO=[0.0D,0.0D,0.0D] ;This is the location of the sun in MAG coordinates sunMAG=[0.0D,0.0D,0.0D] ;This is the position of the sun in geographic spherical coords rsun=0.0D sunlati=0.0D sunlongi=0.0D ;To do the rotation, I change the vertices of each polygon from ;MLT and mlat coords to MAG coords and then geographic coords ;This is the position of the first polygon vertice in MAG xMAG1=[0.0D,0.0D,0.0D] ;This is the position of the first polygon vertice in cartesian GEO xGEO1=[0.0D,0.0D,0.0D] xMAG2=[0.0D,0.0D,0.0D] xGEO2=[0.0D,0.0D,0.0D] xMAG3=[0.0D,0.0D,0.0D] xGEO3=[0.0D,0.0D,0.0D] xMAG4=[0.0D,0.0D,0.0D] xGEO4=[0.0D,0.0D,0.0D] ;This is the 2 magnetic longitudes that define the polygon maglon1=0.D0 maglon2=0.D0 ;This is the 2 magnetic latitudes that define the polygon maglat1=0.D0 maglat2=0.D0 ;This is the coordinates of each of the polygon vertices in sperical geographic coords rgeo=0.0D rlat1=0.0D rlong1=0.0D rlat2=0.0D rlong2=0.0D rlat3=0.0D rlong3=0.0D rlat4=0.0D rlong4=0.0D ;psi is the angle of the dipole that is output as part of the coordiante transformation psi=0.0D ; FIXME: I put in a fixed time for debuggin purposes. This should be replaced by the forecast time iyr=2009 idoy=19 secs=18.0*3600.0D ;To change the polygon vertices from MLT and magnetic latitude we first need the magnetic longitude of the sun ;This changes the GSM coords of the sun to geographic result = call_external('c:\onera\onera_desp_lib_Win32_x86.dll', 'gsm2geo_', iyr,idoy,secs,psi,sunGSM,sunGEO, /f_value) ;Then change the geographic coordinates of the sun to MAG coords result = call_external('c:\onera\onera_desp_lib_Win32_x86.dll', 'geo2mag_', iyr,sunGEO,sunMAG, /f_value) ;Then change the MAG coordinates of the sun from cartesiam to spherical result = call_external('c:\onera\onera_desp_lib_Win32_x86.dll', 'car2sph_', sunMAG,rsun,sunlati,sunlongi, /f_value) ;Now change the MLT coords if the polygon to magnetic longitude maglon1=(mlt-12.0)*15.0+sunlongi ;This is the magnetic longitude for mlt+.25 maglon2=(mlt+0.250-12.0)*15.0+sunlongi r=1.0D ;change the two magnetic latitudes to a double maglat1=mlat1*1.0D maglat2=mlat2*1.0D ;Now I find the GEO cooords for maglat1 at maglon1 ;First change spherical MAG coords to cartesian result = call_external('c:\onera\onera_desp_lib.dll', 'sph2car_', r,maglat1,maglon1,xMAG1, /f_value) ;Now change from MAG to GEO result = call_external('c:\onera\onera_desp_lib.dll', 'mag2geo_', iyr,xMAG1,xGEO1, /f_value) ;Now change from GEO cartesian to GEO spherical result = call_external('c:\onera\onera_desp_lib.dll', 'car2sph_', xGEO1,rgeo,rlat1,rlong1, /f_value) ;Now I find the GEO cooords for maglat1 at maglon2 result = call_external('c:\onera\onera_desp_lib.dll', 'sph2car_', r,maglat1,maglon2,xMAG2, /f_value) result = call_external('c:\onera\onera_desp_lib.dll', 'mag2geo_', iyr,xMAG2,xGEO2, /f_value) result = call_external('c:\onera\onera_desp_lib.dll', 'car2sph_', xGEO2,rgeo,rlat2,rlong2, /f_value) ;Now I find the GeO coords for maglat2 at maglon1 result = call_external('c:\onera\onera_desp_lib.dll', 'sph2car_', r,maglat2,maglon1,xMAG3, /f_value) result = call_external('c:\onera\onera_desp_lib.dll', 'mag2geo_', iyr,xMAG3,xGEO3, /f_value) result = call_external('c:\onera\onera_desp_lib.dll', 'car2sph_', xGEO3,rgeo,rlat3,rlong3, /f_value) ;Now I find the GeO coords for maglat2 at maglon2 result = call_external('c:\onera\onera_desp_lib.dll', 'sph2car_', r,maglat2,maglon2,xMAG4, /f_value) result = call_external('c:\onera\onera_desp_lib.dll', 'mag2geo_', iyr,xMAG4,xGEO4, /f_value) result = call_external('c:\onera\onera_desp_lib.dll', 'car2sph_', xGEO4,rgeo,rlat4,rlong4, /f_value) ;ang1 is just mlt in radians ;ang1 = ((mlt)/24.0) * (2.0*!pi) ;For the plot it makes teh most sens to keep the sun in ;a fixed location and rotate the Earth. I will keep the sun at ;the top of the plot ;ang1-ang4 are the angles of the polygon vertices. I subtract 90degrees ; so that the sun is at the top. I add teh secs of day so that the coordiantes ; rotate appropriately ang1=(rlong1-90.0+secs*180.0/(3600*12.0))*(!pi/180.0); ;ang2 is the current mlt +.25 ;ang2 = ((mlt+0.25)/24.0) * (2.0*!pi) ang2=(rlong2-90.0+secs*180.0/(3600.0*12.0))*(!pi/180.0) ang3=(rlong3-90.0+secs*180.0/(3600.0*12.0))*(!pi/180.0) ang4=(rlong4-90.0+secs*180.0/(3600.0*12.0))*(!pi/180.0) ;r1=abs(mlat1-90.)/40. ;r2=abs(mlat2-90.)/40. ; ;r1-r4 is the radial distance in degrees r1=abs(rlat1-90.)/40. r2=abs(rlat2-90.)/40. r3=abs(rlat3-90.)/40. r4=abs(rlat4-90.)/40. ;Now this is the x and y values of each polygon x1 = float(cos(ang1)*r1) y1 = float(sin(ang1)*r1) x2 = float(cos(ang2)*r2) y2 = float(sin(ang2)*r2) x3 = float(cos(ang3)*r3) y3 = float(sin(ang3)*r3) x4 = float(cos(ang4)*r4) y4 = float(sin(ang4)*r4) fill_col = efact*(jeave - je_min) if( jeave eq 0. )then fill_col = 0 if( fill_col gt 250.) then fill_col = 250 if( fill_col lt 0.) then fill_col = 0 ;I chose fill_col gt 10 so that all the black is not plotted if (fill_col gt 10.) then begin polyfill,[x1,x2,x4,x3],[y1,y2,y4,y3],col=fill_col endif ;printf,12,format=out1,mlt_real,mlat_real,jeave ;JGREEN============================================= printf,12,format=out1,mlt_real,mlat_real,jeave je_sum(i) = je_sum(i) + jeave if( (mlt_real ge 7.) and (mlt_real le 19.) and $ (abs(mlat_real) le 63.) )then begin jeave = 0. endif mlt1 = mlt_real mlt2 = mlt1 + 0.25 A_earth = area_sphere(mlt1,mlt2,mlat1,mlat2) ;1 erg = 10**-7 joule but converting cm2 to km2 gives 10**10 ;converting watts to GW is a factor of 10-9 if( n_or_s ne 4 )then begin if( (je_type eq 1) or (je_type eq 2) )then begin power_hemi = power_hemi + (1.0e-6)*A_earth*jeave ;in GW imlt = fix(mlt_real) power_mlt(imlt) = power_mlt(imlt) + (1.0e-6)*A_earth*jeave endif if( (je_type eq 3) or (je_type eq 4) )then begin power_hemi = power_hemi + (1.0e10)*A_earth*jeave imlt = fix(mlt_real) power_mlt(imlt) = power_mlt(imlt) + (1.0e10)*A_earth*jeave endif endif endfor endfor dS = power_hemi ; ; draw mlat circles ;JGREEN for r=0.125, 1.0, 0.125 do begin x = r*COS(THETA) y = r*SIN(THETA) oplot, X,Y,thick=2.0,color=255 endfor mlt = 1.2 - 6.0 ; rotate noon to top ang = mlt/24.0 * 2.0*!pi for mlat=50,80,5 do begin ; label mlat ticks r=abs(mlat-89.0)/40.0 x = cos(ang)*r y = sin(ang)*r xyouts,x,y,string(mlat,format='(i2.2)'),color=255 endfor ; draw mlt lines and label delta_theta = 2*!pi/(24.) t = 0. for i = 0,23 do begin ;t = i*delta_theta - !pi/2.0 t = i*delta_theta-!pi/2.0 if ((i mod 3) eq 0) then begin ;x= cos(t) x= cos(t+secs*!pi/(3600.0*12.0)) y= sin(t+secs*!pi/(3600.0*12.0)) ;y= sin(t) oplot,[0.0,x],[0.0,y],thick=2.0,color=255 x= 1.05*cos(t+secs*!pi/(3600.0*12.0)) y= 1.05*sin(t+secs*!pi/(3600.0*12.0)) xyouts,x,y,string(i*180/12,format='(i3.2)'),align=.5 endif else begin x = cos(t+secs*!pi/(3600.0*12.0))*1.02 y = sin(t+secs*!pi/(3600.0*12.0))*1.02 oplot,[0.0,x],[0.0,y],thick=2.0,color=255 endelse endfor ;... title map_name = je_title ;xyouts,!d.x_ch_size*2,!d.y_size-4*!d.y_ch_size,size=2,/dev, $ ; map_name black = 252B white = 251B colmax = 250B colmin = 0B maxv = je_max minv = je_min ;if( (je_type eq 3) or (je_type eq 4) )then begin ; maxv = alog10(maxv) ; if( je_min gt 0. )then minv = alog10(je_min) ;endif cb_h = 0.25*!d.y_vsize cb_w = float(!d.x_ch_size) cb = bytarr(cb_w,cb_h) tmp = colmin + findgen(colmax-colmin+1) cbtmp = interpol(tmp,cb_h) for i=0,cb_w-1 do cb(i,*) = cbtmp !P.POSITION = [0.05,0.4,0.05+float(cb_w-1)/!d.x_vsize,0.65] ;oplot,[minv,maxv], /nodata,color=black,xstyle=5, ystyle=5 tv, cb,0.05,0.4, /normal, xsize=cb_w/!d.x_vsize,ysize=cb_h/!d.y_vsize if( je_type le 2 )then begin xyouts,0.07,0.4, /normal, string(minv,format='(f5.2)') xyouts,0.07,0.65, /normal, string(maxv,format='(f5.2)') xyouts,0.07,0.525, /normal, string((maxv+minv)/2.,format='(f5.2)') if( n_or_s ne 4)then begin xyouts,0.05,0.70, /normal, size=1.5, 'ergs/cm2s' endif endif if( (je_type eq 3) or (je_type eq 4) )then begin xyouts,0.07,0.4, /normal, string(minv,format='(e7.1)') xyouts,0.07,0.65, /normal, string(maxv,format='(e7.1)') xyouts,0.07,0.525, /normal, string((maxv+minv)/2.,format='(e7.1)') if( n_or_s ne 4 )then begin xyouts,0.07,0.70, /normal, '1/cm2 s' endif endif if( (je_type eq 5) or (je_type eq 6) )then begin xyouts,0.07,0.4, /normal, string(minv,format='(f6.0)') xyouts,0.07,0.65, /normal, string(maxv,format='(f6.0)') xyouts,0.07,0.525, /normal, string((maxv+minv)/2.,format='(f6.0)') if( n_or_s ne 4 )then begin xyouts,0.05,0.70, /normal, 'eV' endif endif if( (je_type le 4) and (n_or_s ne 4) )then begin printf,12,'Totals:' grand_total = 0. for i=0,nmlt-1 do begin mlt = float(i)*0.25 printf,12,mlt,je_sum(i) grand_total = grand_total + je_sum(i) endfor printf,12,'Grand total: ',grand_total if(je_type le 2 )then printf,12,'Hemispherical power (GW): ',power_hemi ;endfor printf,12,'Flux integrated over latitude (fixed 1-h MLT):' for i=0,23 do begin printf,12,i,power_mlt(i) endfor power_day = total(power_mlt[6:17]) power_night = total(power_mlt[0:5]) + total(power_mlt[18:23]) printf,12,'Total day = ',power_day printf,12,'Total night= ',power_night endif close,12 ;xyouts,!d.x_ch_size*2,!d.y_size-4*!d.y_ch_size,size=2,/dev, $ ; map_name xyouts,0.07,0.92,/normal, size=2, map_name if ( n_elements( FORECAST_TIME ) eq 1 ) then begin xyouts, 0.07, 0.88, /normal, size=2, 'Forecast '+string( forecast_time, format='(I03)' ) + $ ' minutes (+/- ' + string( forecast_std, format='(F4.1)' ) + '), ' endif if( n_or_s ne 4 )then begin if( je_type le 2 )then begin xyouts,0.70,0.88, /normal, size=2, string(power_hemi,format='(f6.1)') xyouts,0.85,0.88, /normal, size=2, 'GW' endif if( (je_type eq 3) or (je_type eq 4) )then begin xyouts,0.63,0.88, /normal, size=2, string(power_hemi,format='(e8.1)') xyouts,0.80,0.88, /normal, size=2, '/s' endif endif season_string = ['winter','spring','summer','fall',' ','n. wint',$ 'n. spring','n. summer','n. fall'] ; FIXME: correct this next line. season_epoch.pro calls with an incorrect value of always 3, i.e. 'fall' ;xyouts,0.16,0.82, /normal, size=2, season_string(sf) if ( ~ n_elements( NO_SHOW_STATS_YEAR_RANGE ) ) then begin xyouts,0.65,0.82, /normal, size=2, string(year0,format='(i4)') xyouts,0.80,0.82, /normal, size=2, string(doy0,format='(i3)') xyouts,0.65,0.77, /normal, size=2, string(year1,format='(i4)') xyouts,0.80,0.77, /normal, size=2, string(doy1,format='(i4)') endif xyouts,0.05,0.08, /normal, size=2, string(files_done,format='(i5)') xyouts,0.17,0.08, /normal, size=2, 'Satellite-days' Map_Set, /Stereographic, 90, 255, /Continents, /USA, /Isotropic, $ con_color=50, /Horizon, limit=[50,-180,90,180] device,/close return END pro remove_salt_pepper,jin,jout,n3 ;routine to reduce salt-and-pepper noise ;jin is nmlt by nmlat matrix in, jout is returned ;n3 is observations in each bin ;Patrick Newell, January 2009 nmlt = 96 nmlat = 160 jout = jin for i = 0, nmlt-1 do begin for j = 0,nmlat-1 do begin if( n3(i,j) eq 0 )then goto,next_lat jplus1 = mlat_plus(j,1) jplus2 = mlat_plus(j,2) jminus1 = mlat_plus(j,-1) jminus2 = mlat_plus(j,-2) p1 = mlt_plus(i,1) p2 = mlt_plus(i,2) j0 = jin(i,j) jp1 = jin(p1,j) jp2 = jin(p2,j) n30 = n3(i,j) n3p1 = n3(p1,j) n3p2 = n3(p2,j) m1 = mlt_plus(i,-1) m2 = mlt_plus(i,-2) jm1 = jin(m1,j) jm2 = jin(m2,j) n3m1 = n3(m1,j) n3m2 = n3(m2,j) ja = jin(i,jplus1) jaa = jin(i,jplus2) jb = jin(i,jminus1) jbb = jin(i,jminus2) n3a = n3(i,jplus1) n3b = n3(i,jminus1) n3aa = n3(i,jplus2) n3bb = n3(i,jminus2) n10 = 0 if( (j0 ge 10.*ja) and (n3a gt 0) )then n10 = n10 + 1 if( (j0 ge 10.*jaa) and (n3aa gt 0) )then n10 = n10 + 1 if( (j0 ge 10.*jb) and (n3b gt 0) )then n10 = n10 + 1 if( (j0 ge 10.*jbb) and (n3bb gt 0) )then n10 = n10 + 1 if( (j0 ge 10.*jp1) and (n3p1 gt 0) )then n10 = n10 + 1 if( (j0 ge 10.*jp2) and (n3p2 gt 0) )then n10 = n10 + 1 if( (j0 ge 10.*jm1) and (n3m1 gt 0) )then n10 = n10 + 1 if( (j0 ge 10.*jm2) and (n3m2 gt 0) )then n10 = n10 + 1 if( n10 ge 7 )then begin jout(i,j) = 0. w0 = 3. wp1 = 1./2. wp2 = 1./8. wm1 = 1./2. wm2 = 1./8. wa = 1./3. wb = 1./3. waa = 1./27. wbb = 1./27. if( n30 lt 3 )then w0 = 0. if( n3p1 lt 3 )then wp1 = 0. if( n3p2 lt 3 )then wp2 = 0. if( n3m1 lt 3 )then wm1 = 0. if( n3m2 lt 3 )then wm2 = 0. if( n3a lt 3 )then wa = 0. if( n3b lt 3 )then wb = 0. if( n3aa lt 3)then waa = 0. if( n3bb lt 3)then wbb = 0. wsum = w0 + wp1 + wp2 + wm1 + wm2 + wa + wb + waa + wbb if( wsum gt 0. )then begin j_smooth = w0*j0 + wp1*jp1 + wp2*jp2 + wm1*jm1 + wm2*jm2 + $ wa*ja + wb*jb + waa*jaa + wbb*jbb n_smooth = w0*n30 + wp1*n3p1 + wp2*n3p2 + wm1*n3m1 + wm2*n3m2 + $ wa*n3a + wb*n3b + waa*n3aa + wbb*n3bb norm = wsum j_smooth = j_smooth/norm n_smooth = fix(0.5 + n_smooth/norm) jout(i,j) = j_smooth n3(i,j) = n_smooth next_lat: dummy = 0 endif endif endfor ;j loop (MLT) endfor ;i loop (MLAT) return end ;pro make_movie,