; Main: rampcontour.pro
;
; Hawkeye ramp contours for ions and electrons
;
; to run
;
; IDL>.run rampcontour.pro
;
; and follow the instruction
;
; SHC 11/1/96 adapted from S. Boardsen's binpolar.pro
; Last revision 11/14/96
;
;--------------------------------------------------------------------------
function pr_phase_space,rad,n
si=size(rad)
rad_v=fltarr(si(1))
for i=0,n-1 do begin
if (rad(i) ge 0) and (rad(i) lt 1) then rad_v(i)=148.6 else $
if (rad(i) ge 1) and (rad(i) lt 2) then rad_v(i)=159.7 else $
if (rad(i) ge 2) and (rad(i) lt 3) then rad_v(i)=174.0 else $
if (rad(i) ge 3) and (rad(i) lt 4) then rad_v(i)=189.2 else $
if (rad(i) ge 4) and (rad(i) lt 5) then rad_v(i)=206.5 else $
if (rad(i) ge 5) and (rad(i) lt 6) then rad_v(i)=222.5 else $
if (rad(i) ge 6) and (rad(i) lt 7) then rad_v(i)=242.6 else $
if (rad(i) ge 7) and (rad(i) lt 8) then rad_v(i)=261.1 else $
if (rad(i) ge 8) and (rad(i) lt 9) then rad_v(i)=292.7 else $
if (rad(i) ge 9) and (rad(i) lt 10) then rad_v(i)=357.2 else $
if (rad(i) ge 10) and (rad(i) lt 11) then rad_v(i)=461.8 else $
if (rad(i) ge 11) and (rad(i) lt 12) then rad_v(i)=647.3 else $
if (rad(i) ge 12) and (rad(i) lt 13) then rad_v(i)=966.0 else $
if (rad(i) ge 13) and (rad(i) lt 14) then rad_v(i)=1414.1 else $
if (rad(i) ge 14) and (rad(i) lt 15) then rad_v(i)=2065.4 else $
if (rad(i) ge 15) and (rad(i) lt 16) then rad_v(i)=2920.9
endfor
return,rad_v
end
;--------------------------------------------------------------------------
;
; convert 2-d distribution function to energy-angle
;
function pr_energy_angle,rad,n
si=size(rad)
rad_ev=fltarr(si(1))
for i=0,n-1 do begin
if (rad(i) ge 0) and (rad(i) lt 1) then rad_ev(i)=116. else $
if (rad(i) ge 1) and (rad(i) lt 2) then rad_ev(i)=134. else $
if (rad(i) ge 2) and (rad(i) lt 3) then rad_ev(i)=159. else $
if (rad(i) ge 3) and (rad(i) lt 4) then rad_ev(i)=188. else $
if (rad(i) ge 4) and (rad(i) lt 5) then rad_ev(i)=224. else $
if (rad(i) ge 5) and (rad(i) lt 6) then rad_ev(i)=260. else $
if (rad(i) ge 6) and (rad(i) lt 7) then rad_ev(i)=309. else $
if (rad(i) ge 7) and (rad(i) lt 8) then rad_ev(i)=358. else $
if (rad(i) ge 8) and (rad(i) lt 9) then rad_ev(i)=450. else $
if (rad(i) ge 9) and (rad(i) lt 10) then rad_ev(i)=670. else $
if (rad(i) ge 10) and (rad(i) lt 11) then rad_ev(i)=1120. else $
if (rad(i) ge 11) and (rad(i) lt 12) then rad_ev(i)=2200. else $
if (rad(i) ge 12) and (rad(i) lt 13) then rad_ev(i)=4900. else $
if (rad(i) ge 13) and (rad(i) lt 14) then rad_ev(i)=10500. else $
if (rad(i) ge 14) and (rad(i) lt 15) then rad_ev(i)=22400. else $
if (rad(i) ge 15) and (rad(i) lt 16) then rad_ev(i)=44800.
rad_ev(i)=alog10(rad_ev(i))
endfor
return,rad_ev
end
;--------------------------------------------------------------------------
;
; convert 2-d distribution function to energy-angle
;
function el_phase_space,rad,n
si=size(rad)
rad_v=fltarr(si(1))
for i=0,n-1 do begin
if (rad(i) ge 0) and (rad(i) lt 1) then rad_v(i)=1. else $
if (rad(i) ge 1) and (rad(i) lt 2) then rad_v(i)=148. else $
if (rad(i) ge 2) and (rad(i) lt 3) then rad_v(i)=161. else $
if (rad(i) ge 3) and (rad(i) lt 4) then rad_v(i)=175. else $
if (rad(i) ge 4) and (rad(i) lt 5) then rad_v(i)=191. else $
if (rad(i) ge 5) and (rad(i) lt 6) then rad_v(i)=206. else $
if (rad(i) ge 6) and (rad(i) lt 7) then rad_v(i)=224. else $
if (rad(i) ge 7) and (rad(i) lt 8) then rad_v(i)=241. else $
if (rad(i) ge 8) and (rad(i) lt 9) then rad_v(i)=269. else $
if (rad(i) ge 9) and (rad(i) lt 10) then rad_v(i)=329. else $
if (rad(i) ge 10) and (rad(i) lt 11) then rad_v(i)=428. else $
if (rad(i) ge 11) and (rad(i) lt 12) then rad_v(i)=601. else $
if (rad(i) ge 12) and (rad(i) lt 13) then rad_v(i)=894. else $
if (rad(i) ge 13) and (rad(i) lt 14) then rad_v(i)=1302. else $
if (rad(i) ge 14) and (rad(i) lt 15) then rad_v(i)=1907. else $
if (rad(i) ge 15) and (rad(i) lt 16) then rad_v(i)=2697.
endfor
return,rad_v
end
;--------------------------------------------------------------------------
;
; convert 2-d distribution function to energy-angle
;
function el_energy_angle,rad,n
si=size(rad)
rad_ev=fltarr(si(1))
for i=0,n-1 do begin
if (rad(i) ge 0) and (rad(i) lt 1) then rad_ev(i)=1. else $
if (rad(i) ge 1) and (rad(i) lt 2) then rad_ev(i)=115. else $
if (rad(i) ge 2) and (rad(i) lt 3) then rad_ev(i)=136. else $
if (rad(i) ge 3) and (rad(i) lt 4) then rad_ev(i)=160. else $
if (rad(i) ge 4) and (rad(i) lt 5) then rad_ev(i)=191. else $
if (rad(i) ge 5) and (rad(i) lt 6) then rad_ev(i)=222. else $
if (rad(i) ge 6) and (rad(i) lt 7) then rad_ev(i)=264. else $
if (rad(i) ge 7) and (rad(i) lt 8) then rad_ev(i)=306. else $
if (rad(i) ge 8) and (rad(i) lt 9) then rad_ev(i)=380. else $
if (rad(i) ge 9) and (rad(i) lt 10) then rad_ev(i)=570. else $
if (rad(i) ge 10) and (rad(i) lt 11) then rad_ev(i)=960. else $
if (rad(i) ge 11) and (rad(i) lt 12) then rad_ev(i)=1900. else $
if (rad(i) ge 12) and (rad(i) lt 13) then rad_ev(i)=4200. else $
if (rad(i) ge 13) and (rad(i) lt 14) then rad_ev(i)=8900. else $
if (rad(i) ge 14) and (rad(i) lt 15) then rad_ev(i)=19100. else $
if (rad(i) ge 15) and (rad(i) lt 16) then rad_ev(i)=38200.
rad_ev(i)=alog10(rad_ev(i))
endfor
return,rad_ev
end
;--------------------------------------------------------------------------
;
;
function bin1,y0,eprange,dang,interp=interp
y=y0
if n_elements(interp) eq 0 then interp = 0
if eprange(1) le eprange(0) then begin
dummy = min(abs(y.epoch-eprange(0)),istart)
period = y(istart).period
dt = (1./(11.52/period-1)+0)*11.52*1000 ;time for one complete sweep
eprange(1)=eprange(0) + dt
endif
ii=where( (y.epoch ge eprange(0)) and (y.epoch lt eprange(1)) )
if ii(0) ne -1 then begin
y = y(ii)
n=n_elements(ii)
period = y(0).period
dang= (11.52/period - 1)*360
REC = {EPOCH:0D0,flux:fltarr(16),xsunang:fltarr(16) $
,bave:0.,blong:0.,bcolat:0.}
if interp eq 1 then begin
m = (n/2)*2
g = REPLICATE(REC,M)
; bininterp,y,epoch,xsunang,flux
g.flux = flux
g.xsunang = xsunang
g.epoch = epoch
endif else begin
g = REPLICATE(REC,n)
g.flux=y.flux
g.xsunang=y.xsunang
g.epoch=y.epoch
endelse
g.epoch = y.epoch
ii = where( y.bave ge 0)
n= n_elements(ii)
g.bave = total( y(ii).bave )/(n)
g.blong = total( y(ii).blong )/(n)
g.bcolat = total( y(ii).bcolat )/(n)
endif
return,g
end
;--------------------------------------------------------------------------
;
;
pro clip,z,min_vel,min_value
si=size(z)
for i=0,si(1)-1 do $
for j=0,si(2)-1 do $
if(sqrt(i*i+j*j) lt min_vel) and (z(i,j) lt 8.) then z(i,j)=min_value
return
end
;--------------------------------------------------------------------------
;
; calibration for low energy channel
;
pro calibration,z,m,n
factor=fltarr(6)
factor(0)=22.52
factor(1)=11.06
factor(2)=5.645
factor(3)=2.538
factor(4)=1.7
factor(5)=0.8761
for i=0,5 do begin
for j=0,n-1 do begin
z(i,j)=z(i,j)+alog10(factor(i))
endfor
endfor
return
end
;--------------------------------------------------------------------------
pro pregion,nx,ny,jplot
iplot= jplot mod (nx*ny)
common deviceTypeC, deviceType
;
;if (deviceType eq 0) or (deviceType eq 1) then ratio=8.5/11.0*22.85/21. $
;if (deviceType eq 0) or (deviceType eq 1) then ratio=8.5/11.0*0.991 $
;
if (deviceType ge 1) and (deviceType le 3) then ratio=8.5/11.0*0.993 $
else if (deviceType eq 0 or deviceType eq 4) then ratio=8.5/11.0*1.05 $
else ratio=1.0
;
case devicetype of
0:begin
; xm= [.05,.95*ratio]
; ym= [.05,.95]
xm= [0.05,1.*ratio]
ym= [0.05,1.]
end
4:begin
; xm= [.1,.9*ratio]
; ym= [.1,.9]
xm= [0.05,1.*ratio]
ym= [0.05,1.]
end
else: begin
; xm= [.05,.9*ratio]
; ym= [.05,.9]
xm= [0.05,1.*ratio]
ym= [0.05,1.]
end
endcase
dy = (ym(1)-ym(0))/ny*0.8
;dx = (xm(1)-xm(0))/nx
dx = (xm(1)-xm(0))/ny*0.8
iy = fix(iplot/nx)
ix = iplot-iy*nx
y0 = ym(0) + (ny-iy-1)*dy+0.05
y1 = y0+dy*0.95
x0 = xm(0) + ix*dx
x1 = x0+dx*0.95
!p.region=[x0,y0,x1,y1]
!p.position=[x0,y0,x1,y1]
print,'plot position:',x0,y0,x1,y1
end
;--------------------------------------------------------------------------
pro arrow,xv,yv,color=color
if n_elements(color) eq 0 then color = 0
p0=[xv(0),yv(0)]
p1=[xv(1),yv(1)]
;draw arrow
dx = xv(1)-xv(0)
dy = yv(1)-yv(0)
c = [ [dx,dy],[-dy,dx] ]
tip = [1,0]
r = tip + [-1,1]/8./sqrt(2)
l = tip + [-1,-1]/8./sqrt(2)
tip = p0 + c#tip
r = p0 + c#r
l = p0 + c#l
oplot,[xv(0),tip(0)],[yv(0),tip(1)],color=color ;draw shaft
oplot,[tip(0),r(0)],[tip(1),r(1)],color=color ;draw arrow head
oplot,[tip(0),l(0)],[tip(1),l(1)],color=color ;draw arrow head
end
;--------------------------------------------------------------------------
pro grid,chan,rotate,reflect,black,white
;draw polar grid
if n_elements(rotate) eq 0 then rotate = 0
if n_elements(reflect) eq 0 then reflect = 0
n1=n_elements(chan)
a= indgen(73)*2*!pi/72
x = cos(a) & y = sin(a)
for i=0,n1-1,4 do begin
if(i eq n1-1) then begin
linestyle=0
color=black ;black
endif else begin
linestyle=1
color=white ;white
endelse
oplot,chan(i)*x,chan(i)*y,color=color,linestyle=linestyle,thick=1
endfor
;ang=[0,90,180,270]
;lab=['X!DGSM!N',' ',' ',' ']
;ali=[0,.5,1,.5]
;for i=0,3 do begin
; x = [0,cos(ang(i)/!radeg+rotate)]*chan(n1-1)
; y = [0,sin(ang(i)/!radeg+rotate)]*chan(n1-1)
; if reflect eq 1 then y=-y
; iali = fix( ((atan(y(1),x(1))*!radeg +405) mod 360)/90 )
; oplot,x,y,color=255,linestyle=1
; xyouts,x(1),y(1),lab(i), align=ali(iali)
;endfor
end
;--------------------------------------------------------------------------
function smooth2d,z
si=size(z)
zbuf=fltarr(si(1),si(2))
s1m=si(1)-1
s2m=si(2)-1
for i=1,si(1)-2 do begin
for j=1,si(2)-2 do begin
zbuf(i,j)=0.5*z(i,j)+$
0.5*(z(i-1,j)+z(i+1,j)+z(i,j-1)+z(i,j+1))/4.
endfor
endfor
for i=1,si(1)-2 do begin
zbuf(i,0)=0.5*z(i,0)+0.5*(z(i-1,0)+z(i+1,0)+z(i,1))/3.
zbuf(i,s2m)=0.5*z(i,s2m)+0.5*(z(i-1,s2m)+z(i+1,s2m)+z(i,s2m-1))/3.
endfor
for j=1,si(2)-2 do begin
zbuf(0,j)=0.5*z(0,j)+0.5*(z(0,j-1)+z(0,j+1)+z(1,j))/3.0
zbuf(s1m,j)=0.5*z(s1m,j)+0.5*(z(s1m,j-1)+z(s1m,j+1)+z(s1m-1,j))/3.
endfor
zbuf(0,0)=0.5*z(0,0)+0.5*(z(1,0)+z(0,1))/2.
zbuf(0,s2m)=0.5*z(0,s2m)+0.5*(z(1,s2m)+z(0,s2m-1))/2.
zbuf(s1m,0)=0.5*z(s1m,0)+0.5*(z(s1m,1)+z(s1m-1,0))/2.
zbuf(s1m,s2m)=0.5*z(s1m,s2m)+0.5*(z(s1m-1,s2m)+z(s1m,s2m-1))/2.
return,zbuf
end
;------------------------------------------------------------------
; average
;
pro rampave,y,ybuf
si=size(y)
for i=0,si(1)-1 do begin
ybuf(i)=0.
for j=0,si(2)-1 do begin
ybuf(i)=ybuf(i)+10.^y(i,j)
endfor
if(ybuf(i) ne 0) then ybuf(i)=ybuf(i)/si(2) $
else ybuf(i)=ymin
endfor
return
end
;------------------------------------------------------------------
;
; line plot of flux density vs. energy
;
pro lineplot,x,y,xmin,xmax,ymin,ymax,title
si=size(y)
plot,x(1:si(2)-1),y(0,1:si(2)-1),$
xrange=[xmin,xmax],yrange=[ymin,ymax],/ylog,/xlog,$
xtitle='Log Energy (eV)',$
ytitle='Log differential Flux (/cm!U2!N/str/eV)',$
position=[0.15,0.15,0.9,0.9],title=title
xyouts,x(11),y(0,11),string(1,format='(i1.1)')
for i=1,si(1)-1 do begin
oplot,x(1:si(2)-1),y(i,1:si(2)-1)
xyouts,x(11),y(i,11),string(i+1,format='(i1.1)')
endfor
end
;
;---------------------------------------------------------------
;main
;---------------------------------------------------------------
close,/all
on_ioerror, badinput
rotate=!pi ;rotate 180 so tw flow is to the right
nxy=3
common deviceTypeC, deviceType
;!p.multi=[0,nxy,2]
!x.omargin = 0
!y.omargin = 0
deviceopen
if (devicetype ge 1 and devicetype le 3) then begin
print,' '
print,'The output file will be idl.ps'
print,' '
endif
;if (devicetype eq 1) then begin
; !x.omargin=[1,1]
; !y.omargin=[5,5]
;endif
path=' '
if n_elements(yload) lt 10 then begin
;
print,' '
print,'Enter the path (directory) to the ramp data file?'
print,' (Hit return if under current directory) '
read,path
print,' '
print,'input 0 to plot ion or '
read, ' 1 to plot electrons:',ielect
if(ielect lt 0) then ielect=0 $
else if(ielect gt 1) then ielect=1
; ielect=0
if ielect eq 1 then rampload,path=path,yload,/electrons,/xdr $
else rampload,yload,path=path,/xdr
;
print,' '
print,'input 1 to calibrate data to [s^3/km^6] or'
read, ' 2 to calibrate data to [/cm^2/str/eV]:',icalibrate
if(icalibrate lt 1) then icalibrate=1 $
else if(icalibrate gt 2) then icalibrate=2
if (icalibrate eq 1) then begin
if(ielect eq 0) then lepcal,yload,/ions,/psd,/ramp $
else lepcal,yload,/electrons,/psd,/ramp
endif else begin
if(icalibrate eq 2) then begin
if(ielect eq 0) then lepcal,yload,/ions,/ramp $
else lepcal,yload,/electrons,/ramp
endif
endelse
print,' '
print,'input 1 to plot flux vs. velocity space'
read, ' 2 to plot flux vs. Log energy ?',phase_energy
lplot_yon=1
if(ielect eq 1 and icalibrate eq 2 and phase_energy eq 2) then begin
print,' '
read,'Do you want to do xyplot for flux vs. energy?(0- yes, 1- no)? ',$
lplot_yon
endif
; read,'input 1 to despin data:',ispin
; if ispin eq 1 then despin,yload
nel = n_elements(yload.epoch)
endif
if ielect eq 1 then species='ELECTRONS' else species='IONS'
p0 = yload(10).sunang(0)/!radeg
xp0 = yload(10).xsunang(0)/!radeg
ysunlat = 90-asin(cos(p0)/cos(xp0))*!radeg
yr=0&doy=0&hr=0&min=0
iset = 1
i0 = 0
if iset eq 0 then begin
; print,' time interval:',pdate(yload(0).epoch),' ',pdate(yload(nel-1).epoch)
read,' type in start time:',yr,doy,hr,min
ep0 = epoch0(yr,doy,hr*3600d0+min*60d0)
dummy=min(abs(ep0-g.epoch),i0)
endif
erase
ztop = 1
;read,'input 1 to place z-axis at top:',ztop
reflect=1
;read,'input 1 to reflect image about x-axis:',reflect
ismooth=0
read,'do smoothing? (0 - yes, 1 - no)',ismooth
;read,'input # of plots:',np
xsquare=fltarr(4)
ysquare=fltarr(4)
xsquare(0)=-2.
xsquare(1)=2.
xsquare(2)=2.
xsquare(3)=-2.
ysquare(0)=-1.
ysquare(1)=-1.
ysquare(2)=1.
ysquare(3)=1.
labels=' '
aa=strarr(6)
aa(0)='a'
aa(1)='b'
aa(2)='c'
aa(3)='d'
aa(4)='e'
aa(5)='f'
np=6
zb_buf1=fltarr(6,16)
for ij=0,np-1 do begin
; !p.noerase = 1
ik=i0+ij
; if devicetype eq 0 then if (ij mod nxy^2) eq 0 then begin
; window,ij/nxy^2
; erase
; endif else if (ij mod nxy^2) eq 0 then erase
;
; if devicetype eq 1 then if (ij mod nxy^2) eq 0 then erase
if iset eq 1 then begin
; print,' time interval:',pdate(yload(0).epoch),' ',pdate(yload(nel-1).epoch)
read,' type in start time yr,doy,hr,min,sec:',yr,doy,hr,min,sec
; read,' type in start time hr, min, sec:',hr,min,sec
ep0 = epoch0(yr,doy,hr*3600d0+min*60d0+sec)
; read,' type in stop time yr,doy,hr,min,sec:',yr,doy,hr,min,sec
; ep1 = epoch0(yr,doy,hr*3600d0+min*60d0+sec)
yr=0
if yr eq 0 then begin
dummy = min(abs(ep0-yload.epoch),ik)
nt = fix(1./( (11.52/yload(ik).period) mod 1) +1)
ep1 = ep0 + nt*11.52d3
endif
dummy=min(abs(ep0-yload.epoch),ik)
endif
;*********************************************
; read,' label (one character)? ',labels
eprange=[ep0,ep1]
print,pdate(eprange(0))
g=bin1(yload,eprange,dang)
dt=dang/!radeg
ng=n_elements(g.epoch)
epoch=g.epoch
;*********************************************
; pregion,nxy,nxy,ij
;
pregion,3,2,ij
print,'plot #',ij
tsum = string(fix( (g(ng-1).epoch-g(0).epoch)/1000. ), $
format='("!7D!xT:",i4," s")' )
title = pdate(g(0).epoch,icase=2) +' '
sunlat = lep_ang_1(epoch=g(0).epoch,northgsm=north,yp_gsm=dusk)
if ztop eq 1 then begin ;place + z axis at top of plot
az = north(1) + rotate*!radeg
aa = 90
if reflect eq 1 then aa = -90
rotate1 = aa - az
rotate = (rotate + rotate1/!radeg)
endif
; plot flux part of code
if(icalibrate eq 1) then z=alog10(g.flux)+30. $
else z=alog10(g.flux)
; print,z
;***********************************************
; calibration,z,si(1),si(2)
;**********************************************
; zrange=[min(alog10(g.flux),max=max),max]
; zrange=[-.2,4.4]
;
zrange=[0.,6.]
if(icalibrate eq 1) or (icalibrate eq 2) then zb=z $
else zb = byteconvert(z,zrange=zrange)
si = size(zb)
chan = indgen(17)
ysc = (!d.y_vsize/!d.y_px_cm)/(!d.x_vsize/!d.x_px_cm)
yrange = [-1,1]*20
xrange = [-1,1]*20/ysc
;
;********************************************************************
sunlat = lep_ang_1(epoch=g(0).epoch,northgsm=north,yp_gsm=dusk)
print,' sunlat:',sunlat,ysunlat
anorth= north(1)/!radeg+rotate
adusk= dusk(1)/!radeg+rotate
print,'Zgsm= ',north(0),north(1),' Ygsm= ',dusk(0),dusk(1),rotate
ijv=0
zb_buf2=fltarr(16)
zbuf=fltarr(16*ng)
zbuf0=fltarr(16*ng)
rad=fltarr(16*ng)
theta=fltarr(16*ng)
for iv=0,15 do begin
t=g.xsunang(iv)/!radeg
t=t+rotate
for jv=0,ng-1 do begin
; ta=fix(t(jv)/!pi/2.*ng)
; if (ta lt 0) then ta=ta+ng $
; else if(ta gt ng) then ta=ta-ng
rad(ijv)=iv*1.0+0.5
theta(ijv)=t(jv)
zbuf(ijv)=zb(iv,jv)
ijv=ijv+1
endfor
endfor
; for iv=0,15 do begin
; for jv=1,ng-2 do begin
; if(zbuf(iv,jv) le 1.) then $
; zbuf(iv,jv)=(zbuf(iv,jv-1)+zbuf(iv,jv+1))*0.5
; endfor
; endfor
; rad=findgen(16)+0.5
; theta=findgen(ng)*2.*!pi/ng
; zbuf1=polar_surface(zbuf,rad,theta,/grid)
;
if(phase_energy eq 1) then begin
; convert to phase space
if(ielect eq 0) then rad_v=pr_phase_space(rad,16*ng) $
else rad_v=el_phase_space(rad,16*ng)
endif else begin
; convert to log energy
if(ielect eq 0) then rad_logeV=pr_energy_angle(rad,16*ng) $
else rad_logeV=el_energy_angle(rad,16*ng)
endelse
;
; do the averaging of ramp contours
;
rampave,zb,zb_buf2
for i=0,15 do zb_buf1(ij,i)=zb_buf2(i)
;
;
; output for Lun
;openw,50,'lun.asc'
;
; ijv=0
; for iv=0,15 do begin
; for jv=1,ng-2 do begin
; printf,50,rad_v(ijv),theta(ijv),zbuf(ijv)
; ijv=ijv+1
; endfor
; endfor
; printf,50,99999.,99999.,99999.
;
if(phase_energy eq 1) then begin
;
; for velocity-angle-flux.
;
if(ielect eq 0) then begin ; for protons
zbuf1=polar_surface(zbuf,rad_v,theta,$
bounds=[-1000.,-1000.,1000.,1000.],spacing=25.)
zbuf3=reverse(congrid(zbuf1,41,41),2);reflect along x (Horizontal)-axis
endif else begin ; for electrons
zbuf1=polar_surface(zbuf,rad_v,theta,$
bounds=[-2000.,-2000.,2000.,2000.],spacing=50.)
zbuf3=reverse(congrid(zbuf1,41,41),2);reflect along x (Horizontal)-axis
endelse
endif else begin
;
; for energy-angle-flux.
;
if(ielect eq 0) then begin ; for protons
zbuf1=polar_surface(zbuf,rad_logeV,theta,$
bounds=[-5.,-5.,5.,5.],spacing=0.125)
zbuf3=reverse(congrid(zbuf1,41,41),2);reflect along x (Horizontal)-axis
endif else begin ; for electrons
zbuf1=polar_surface(zbuf,rad_logeV,theta,$
bounds=[-5.,-5.,5.,5.],spacing=0.125)
zbuf3=reverse(congrid(zbuf1,41,41),2);reflect along x (Horizontal)-axis
endelse
endelse
;
; output data for future use
;
; if(ij eq 0) then begin
; zbuf2=polar_surface(zbuf,rad_v,theta,$
; bounds=[-2000.,-2000.,2000.,2000.],spacing=25.)
; zbuf5=reverse(congrid(zbuf2,81,81),2)
;
; openw,30,'dist_fun.asc'
; for i=0,80 do $
; for j=0,80 do $
; printf,30,i*50.-2000.0,j*50.-2000.0,zbuf5(i,j)
; close,30
; endif
;
; cut out the contours at lower energy channels
;
; minvalue=0.0
; min_vel=8
; clip,zbuf3,min_vel,minvalue
;
clabels=indgen(12)*0
cthick=indgen(12)*0+1
clabels(11)=1
cthick(11)=2
clabels(9)=1
cthick(9)=2
clabels(7)=1
cthick(7)=2
clabels(5)=1
cthick(5)=2
clabels(3)=1
cthick(3)=2
clabels(1)=1
cthick(1)=2
;
; zbuf=min_curve_surf(zrec)
; zbuf2=min_curve_surf(zbuf)
;
if(ismooth eq 0) then zbuf4=smooth2d(zbuf3) $
else zbuf4=zbuf3
; print,"**********",zbuf4
;
if(icalibrate eq 1) then begin
; nlev=16
; exp=1.2
; levels=findgen(nlev+1)*0.25+2.
; color=(findgen(nlev))^exp*254./((nlev-1)^exp)
; levels(0)=0
nlev=3
levels=findgen(nlev+1)*0.5+7
color=fltarr(nlev)
color(0)=180
color(1)=150
color(2)=100
endif else begin
nlev=18
exp=0.8
levels=findgen(nlev+1)*10.+50
; color=findgen(19)*255/18
color=(findgen(nlev))^exp*255./((nlev-1)^exp)
levels(0)=0
endelse
;
; for color shading (disabled 11/5/96 for tek users)
; contour,zbuf4,levels=levels,c_color=color,xstyle=4,ystyle=4,/fill,/noerase
;
; overlay contour curves
;
; if(icalibrate eq 1) then levels=findgen(11)*0.25+6.25 $
; else if(icalibrate eq 2) then levels=findgen(11)*0.25+2.75
;
if(ielect eq 0) then begin ; for protons
if(icalibrate eq 1) then levels=findgen(11)*0.5+5.5 $
else if(icalibrate eq 2) then levels=findgen(11)*0.5+2.5
endif else begin ; for electrons
if(icalibrate eq 1) then levels=findgen(11)*1. $
else if(icalibrate eq 2) then levels=findgen(11)*1.-2.
endelse
contour,zbuf4,levels=levels,c_thick=cthick,c_labels=clabels,$
xstyle=4,ystyle=4,/noerase
;
; contour,zrec,levels=levels,max_value=500,xstyle=4,ystyle=4,/fill,/noerase
; levels=findgen(10)*10.+150.
; contour,zrec,levels=levels,max_value=500,xstyle=4,ystyle=4,/noerase
;
fullwidth=40.
half=fullwidth/2.
quar=fullwidth/4.
xcir0=cos(findgen(73)*!pi/36.)
ycir0=sin(findgen(73)*!pi/36.)
xcir1=half*xcir0+half
ycir1=half*ycir0+half
xcir2=quar*1.5*xcir0+half
ycir2=quar*1.5*ycir0+half
xcir3=quar*xcir0+half
ycir3=quar*ycir0+half
xcir4=0.5*quar*xcir0+half
ycir4=0.5*quar*ycir0+half
if(deviceType eq 4) then begin
white=0
black=15
color=black
endif else begin
white=255
black=0
color=black
endelse
a=indgen(36)*2.*!pi/35
;
; block out the unmeasured low energy portion of pie plot
;
; xcir=0.45*quar*cos(a)+half
; ycir=0.45*quar*sin(a)+half
; polyfill,xcir,ycir,color=white
; xcir=3.*cos(a)+half
; ycir=3.*sin(a)+half
; oplot,xcir,ycir,color=color
a=indgen(91)*0.5*!pi/89
xconer=half*cos(a)+half
yconer=half*sin(a)+half
xconer(90)=fullwidth
yconer(90)=fullwidth
polyfill,xconer,yconer,color=white
a=(indgen(91)+90.)*0.5*!pi/89
xconer=half*cos(a)+half
yconer=half*sin(a)+half
xconer(90)=0.
yconer(90)=fullwidth
polyfill,xconer,yconer,color=white
a=(indgen(91)-90.)*0.5*!pi/89
xconer=half*cos(a)+half
yconer=half*sin(a)+half
xconer(90)=fullwidth+0.5
yconer(90)=-0.5
polyfill,xconer,yconer,color=white
a=(indgen(91)+180.)*0.5*!pi/89
xconer=half*cos(a)+half
yconer=half*sin(a)+half
xconer(90)=-0.5
yconer(90)=-0.5
polyfill,xconer,yconer,color=white
; oplot,xcir1,ycir1,color=color,thick=2
oplot,xcir1,ycir1,thick=2
oplot,xcir2,ycir2,linestyle=1
oplot,xcir3,ycir3,linestyle=1
oplot,xcir4,ycir4,linestyle=1
;
; if(ij eq 0) then begin
; oplot,[half+quar+3,half+1.5*quar+3],[2.,2.],$
; color=black,linestyle=0
; oplot,[half+quar+3,half+quar+3],[4.,0.],color=black,linestyle=0
; oplot,[half+1.5*quar+3,half+1.5*quar+3.],[4.,0.],$
; color=black,linestyle=0
; endif else begin
; if(ij eq 3) then begin
; xyouts,half+quar+3,fullwidth-2.5,'250'
; endif else begin
; if(ij eq 4) then xyouts,-1.,fullwidth-2.5,'km/s'
; endelse
; endelse
;
; grid,chan,rotate
;
; draw magnetic field line projection
;
longb= g.blong(0)/!radeg+rotate
aveb=g.bave(0)
colatb=g.bcolat(0)/!radeg
; longb0=0.0
; aveb0=0.0
; colatb0=0.0
; nelem=n_elements(longb)
; sort_index=reverse(sort(aveb))
; for i=0,nelem/2-1 do begin
; longb0=longb0+longb(sort_index(i))
; colatb0=colatb0+colatb(sort_index(i))
; aveb0=aveb0+aveb(sort_index(i))
; endfor
; longb0=longb0/nelem*2
; colatb0=colatb0/nelem*2
; aveb0=aveb0/nelem*2
; aveb0=aveb0*sin(colatb0) ;the mag. of the projection in the spin plane
print,"B ave=",aveb(0)," long=",longb(0)," colat=",colatb(0)
bvecx = [0,cos(longb(0))*aveb(0)*half/50.]
bvecy = [0,sin(longb(0))*aveb(0)*half/50.]
; while (longb0 lt 0.) do longb0 =longb0+!pi*2.0
;
; if(longb0 ge 0.) and (longb0 lt !pi/2.) then begin
; longb0=longb0*(1.-adusk/!pi)+adusk
; bvecx = [0,cos(longb0)*aveb0*half/50.]
; bvecy = [0,sin(longb0)*aveb0*half/50.]
; endif else begin
; if(longb0 ge !pi/2.) and (longb0 lt !pi) then begin
; longb0=!pi-longb0
; longb0=longb0*(1.-adusk/!pi)-adusk
; bvecx = [0,-cos(longb0)*aveb0*half/50.]
; bvecy = [0, sin(longb0)*aveb0*half/50.]
; endif else begin
; if(longb0 ge !pi) and (longb0 lt !pi*1.5) then begin
; longb0=longb0-!pi
; longb0=longb0*(1.+adusk/!pi)+adusk
; bvecx = [0,-cos(longb0)*aveb0*half/50.]
; bvecy = [0,-sin(longb0)*aveb0*half/50.]
; endif else begin
; if(longb0 ge !pi*1.5) and (longb0 lt !pi*2.) then begin
; longb0=!pi*2.-longb0
; longb0=longb0*(1.+adusk/!pi)-adusk
; bvecx = [0, cos(longb0)*aveb0*half/50.]
; bvecy = [0,-sin(longb0)*aveb0*half/50.]
; endif
; endelse
; endelse
; endelse
;
if reflect eq 1 then bvecy=-bvecy
print,'adusk=',adusk,' bvecx=',bvecx(1),' bvecy=',bvecy(1)
bvecx=bvecx+half
bvecy=bvecy+half
; oplot,bvecx,bvecy,color=color
; arrow,bvecx,bvecy,color=color
oplot,bvecx,bvecy
arrow,bvecx,bvecy
x0 = -24 ;!x.crange(0) -3
charsize=.9
if (ij eq 0) then begin
if(icalibrate eq 1) then begin
if(phase_energy eq 1)then begin
xyouts,10,fullwidth+10.,$
species+ ' Log Phase Space Density f(v) (s!U3!N/km!U6!N)',charsize=1.2
endif else begin
xyouts,5,fullwidth+10.,$
species+' Log Density (s!U3!N/km!U6!N) vs. Log energy (eV)',charsize=1.2
endelse
endif else begin
if(phase_energy eq 1)then begin
xyouts,5,fullwidth+10.,$
species+' Log differential Flux (/cm!U2!N/str/eV) vs Velocity',$
charsize=1.2
endif else begin
xyouts,5,fullwidth+10.,$
species+' Log differential Flux (/cm!U2!N/str/eV) vs. Log energy (eV)',$
charsize=1.2
endelse
endelse
endif
if (ij lt 3) then xyouts,quar,fullwidth+5.,title,charsize=charsize $
else xyouts,quar,-5.,title,charsize=charsize
;
; xyouts,x0,16,title,charsize=charsize
; xyouts,x0,14,'RAMP ' + species,charsize=charsize
; xyouts,x0,12,tsum,charsize=charsize
; xyouts,x0,-16,string(north,format='("Z_GSM:(",i3,",",i3,")")'),charsize=charsize
; xyouts,x0,-14,string(dusk,format='("Y_GSM:(",i3,",",i3,")")'),charsize=charsize
; xyouts,x0,-12,string(sunlat,format='("Slat:",i3)'),charsize=charsize
; xyouts,x0,-10,string(90-g(ik).bcolat,format='("Blat:",i3)'),charsize=charsize
a = chan(n_elements(chan)-1)+.5
;
;projected z gsm axis onto s/c spin plane
;
lpos = cpos( (anorth)*!radeg,a)
x0 = lpos(0)
y0 = lpos(1)
if reflect eq 1 then y0=-y0
ali=[0,.5,1,.5]
iali = fix( ((atan(y0,x0)*!radeg +405) mod 360)/90 )
; xyouts,x0,y0,'Z!dGSM!n',align=ali(iali)
; xp =[-1,1]*half*cos(anorth)
; yp =[-1,1]*half*sin(anorth)
xp =[half,half]
yp =[0,fullwidth]
; if reflect eq 1 then yp = -yp
; oplot,xp,yp,color=color,linestyle=1
oplot,xp,yp,linestyle=1
;
;projected y gsm axis onto s/c spin plane
;
lpos = cpos( (adusk)*!radeg,a)
x0 = lpos(0)
y0 = lpos(1)
; if reflect eq 1 then y0=-y0
ali=[0,.5,1,.5]
iali = fix( ((atan(y0,x0)*!radeg +405) mod 360)/90 )
; xyouts,x0,y0,'Y!dGSM!n',align=ali(iali)
xp =[-1,1]*half*cos(adusk)
yp =[-1,1]*half*sin(adusk)
; if reflect eq 1 then yp = -yp
; oplot,xp+half,yp+half,color=color,linestyle=1
oplot,xp+half,yp+half,linestyle=1
;
; sun-earth line
;
; xp=[-1,1]*half*cos(rot_gsm)
; yp=[-1,1]*half*sin(rot_gsm)
; if reflect eq 1 then yp=-yp
; oplot,xp+half,yp+half,color=white,linestyle=1,thick=2
; oplot,[25,25],[0,50],color=white,linestyle=1,thick=2
;
; colorbar,zrange,'log10( cnts/s )',xloc=18.5,ycharsize = .5
;
if(ij eq 0) then begin
gridlabel=strarr(3)
if(ielect eq 0) then begin ; protons
if(phase_energy eq 1) then begin
gridlabel(2)='750'
gridlabel(1)='500'
gridlabel(0)='250'
endif else begin
gridlabel(2)='1000'
gridlabel(1)='100'
gridlabel(0)='10'
endelse
endif else begin ; electrons
if(phase_energy eq 1) then begin
gridlabel(2)='1500'
gridlabel(1)='1000'
gridlabel(0)='500'
endif else begin
gridlabel(2)='1000'
gridlabel(1)='100'
gridlabel(0)='10'
endelse
endelse
if(phase_energy eq 1) then begin
xyouts,half-1,fullwidth+1.,'Vz'
xyouts,fullwidth,half+2.5,'Vy'
endif else begin
xyouts,half-1,fullwidth+1.,'z'
xyouts,fullwidth,half+2.5,'y'
endelse
polyfill,xsquare+half,ysquare+0.5*quar,color=white
xyouts,half-1.5,0.5*quar-0.5,gridlabel(2),charsize=0.8
polyfill,xsquare+half,ysquare+quar,color=white
xyouts,half-1.5,quar-0.5,gridlabel(1),charsize=0.8
polyfill,xsquare+half,ysquare+1.5*quar,color=white
xyouts,half-1.5,1.5*quar-0.5,gridlabel(0),charsize=0.8
endif
; xyouts,2.,38.,string(labels,format='("(",a1,")")')
xyouts,2.,38.,string(ij+1,format='("(",i1.1,")")')
;**********************************************************************
nxy2=nxy^2
if (ij mod nxy2) eq (nxy2-1) then $
if devicetype eq 1 then begin
epa= epoch(max([i0,ik-nxy^2+1]))
epb= epoch(ik+1)
fps = 'hkpolar_'+ rdate(epa) + '_' + rdate(epb) + '.ps'
print,'****creating postscript file:',fps
deviceclose
; spawn,'rename idl.ps '+fps
endif
if (ij mod nxy2) eq (nxy2-1) then $
if( (devicetype eq 0) or (devicetype eq 5) )then begin ;store image in byte map
bytemap = tvrd()
epa= epoch(max([i0,ik-nxy2+1]))
epb= epoch(ik+1)
epr= [epa,epb]
; fgif = 'hkpolar_'+ rdate(epa) + '_' + rdate(epb) + '.gif'
; print,'****creating gif file:',fgif
; write_gif,'scratch:[schen.gif]'+fgif,bytemap
if devicetype eq 5 then erase
endif
endfor ;----------------------- end of main loop -------------------
;
badinput: print,'*** Incorrent input parameters - program terminated ***'
;
; line plot
;
si=size(zb_buf1)
if(si(1) gt 0 and lplot_yon eq 0) then begin
window,1
rad_integer=indgen(16)
rad_eV=el_energy_angle(rad_integer,16)
for i=0,15 do rad_eV(i)=10.^rad_eV(i)
lineplot,rad_eV,zb_buf1,10.,1.e5,0.01,1.e10,species
endif
;close,50
deviceclose
end