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,