;$Author: kenb $
;$Date: 2006-10-11 13:32:51 -0700 (Wed, 11 Oct 2006) $
;$Header: /home/rumba/cdaweb/dev/control/RCS/vectplt.pro,v 1.12 2001/11/26 18:32:30 johnson Exp johnson $
;$Locker: johnson $
;$Revision: 8 $
;+
; NAME: VECTPLT.PRO
;
; PURPOSE: Multi-purpose vector plotting routine for CDF data file
; representation
;
; CALLING SEQUENCE:
;
; vectplt,mlats,mlons,malts,vest,vnrt,mytimes,qflgs,Stitle=station,$
; mcors=mltin,Qmin=qmin,Qmax=qmax,nopolar=nop,Scale=scale,nobin=nob,$
; p0lat=p0lat,p0lon=p0lon,rot=rot,binlat=bin_lat,binlon=bin_lon,$
; limit=limit,latdel=latdel,londel=londel,Alt=alt,Ttitle=thetitle,$
; lthik=lthik,symsiz=symsiz,symcol=symcol,_extra=extras
;
; VARIABLES:
;
; Input:
;
; mlats(*) - an N element array of geographic latitudes
; mlons(*) - an N element array of geographic longitudes
; malts(*) _ an N element array of geographic altitudes
; vest(*) - an N element array of the eastward component velocity vector
; vnrt(*) - an N element array of the northward component velocity vector
; mytimes - an IDL structure comprised of the following: year, day of year,
; month, day of month, and an N element array of seconds of day
; qflgs - an N element array of the quality flag
;
; Keyword Parameters:
;
; Stitle=station - Observing Station Name
; mcors=mltin - Coordinate Transformation Flag
; 0 - Eccentric Dipole Magnetic Local Time
; 1 - Altitude Adjusted Corrected Geomagnetic
; Coordinates (AACGC) MLT
; 2 - Geographic Coordinates
; Qmin=qmin - Minimum exceptable quality flag
; Qmax=qmax - Maximum exceptable quality flag
; nopolar=nop - Disables clockdial display; enables user defined
; projections
; Scale=scale - Set the size of vectors plotted (1000 default)
; nobin=nob - turns off binning and averaging of vectors
; binlat=bin_lat - latitude bin interval
; binlon=bin_lon - longitude bin interval
; p0lat=p0lat - map_set argument; latitude center of map
; p0lon=p0lon - map_set argument; longitude center of map
; rot=rot - map_set argument; rotation of map
; limit=limit - map_set limits
; latdel=latdel - latitude interval
; londel=londel - longitude interval
; Alt=alt - Altitude of coordinate transformation
; Ttitle=thetitle - Title for plot
; lthik=lthik - vector line thickness
; symsiz=symsiz - vector position marker
; symcol=symcol - vector position marker color
;
;
; REQUIRED PROCEDURES:
; This procedure requires the routine GENLIB.PRO
; and a shared object module, LIB_PGM.so, of C and fortran source.
;
;-------------------------------------------------------------------
; History
;
; Revision 1.1 94/10/31 09:49:24 09:49:24 baker (Kile Baker S1G)
;
; Initial version: Based on dn_cdf.pro
;
; 1.0 R. Baldwin HSTX 9/28/95
; Initial version: Based on dn_cdf.pro 10/31/94 (Kile Baker)
; 1.1 R. Baldwin HSTX 10/23/95
; Options included for eccentric & AACGC MLT coordinate
; transformation; geographic coordinates; polar and non-polar
; representations; vector magnitude; quality flag bounds;
; option for binning;
; 1.2 R. Baldwin HSTX 11/15/95
; Title, margin options added; londel & latdel variables
; revised
; 1.3 R. Baldwin HSTX 12/18/95
; Rotation of vectors is a magnetic system from
; geographic to geomagnetic
;
;-------------------------------------------------------------------
function vectplt,mlats,mlons,malts,vest,vnrt,mytimes,qflgs,Stitle=station,$
mcors=mltin,Qmin=qmin,Qmax=qmax,nopolar=nop,Scale=scale,nobin=nob,$
p0lat=p0lat,p0lon=p0lon,rot=rot,binlat=bin_lat,binlon=bin_lon,$
limit=limit,latdel=latdel,londel=londel,Alt=alt,xmargin=xmargin,$
ymargin=ymargin,Ttitle=thetitle,lthik=lthik,symsiz=symsiz,symcol=symcol,$
_extra=extras
; Establish error handler
catch, error_status
if(error_status ne 0) then begin
if((error_status eq -108) or (error_status eq -133)) then begin
print, 'STATUS= No plottable data.'
print, 'ERROR=Error number: ',error_status,' in vectplt.'
print, 'ERROR=Error Message: ', !ERR_STRING
endif else begin
print, 'STATUS= Data cannot be plotted.'
print, 'ERROR=Error number: ',error_status,' in vectplt.'
print, 'ERROR=Error Message: ', !ERR_STRING
endelse
return,-1
endif
icnt=n_elements(mlats)
if (n_elements(mlons) ne icnt or n_elements(malts) ne icnt or $
n_elements(vest) ne icnt or n_elements(vnrt) ne icnt or $
n_elements(mytimes.times) ne icnt or n_elements(qflgs) ne icnt) then $
message, 'Arrays have incorrect size.'
if(n_elements(nop) eq 0) then nop=0
if(n_elements(lthik) eq 0)then lthik = 1.7
if(n_elements(symcol) eq 0)then symcol = 100
if(n_elements(symsiz) eq 0)then symsiz = 0.5
; Set Scale
if(n_elements(scale) ne 0) then Scale=scale else scale=1000
; Set Title
yrst=strmid(strcompress(string(mytimes.year),/remove_all),2,2)
monst=string(mytimes.mon)
if mytimes.mon lt 10 then monst = '0'+monst
dayst=string(mytimes.day)
if mytimes.day lt 10 then dayst = '0'+dayst
; patch
if(!d.name eq 'X') then bkgrd=255 & dtxt=0.5
if(!d.name eq 'Z') then bkgrd=0 & dtxt=1.0
; new code
if n_elements(mltin) eq 0 then mltin=0
if mltin eq 2 then nop=1
if(keyword_set(qmin)) then Qmin=qmin else qmin=0
if(keyword_set(qmax)) then Qmax=qmax else qmax=4
if(n_elements(limit) ne 0) then begin
ymin=limit(0) & xmin=limit(1) & ylim=limit(2) & xlim=limit(3)
endif else begin
ymin=60. & xmin=-180. & ylim=90. & xlim=180.
endelse
; Check for S polar plot
spole=0
if(ymin eq -90. and nop eq 0) then spole=1
if(n_elements(latdel) eq 0) then latdel=10.0
if(n_elements(londel) eq 0) then londel=45.0
; Set Binning intervals
if keyword_set(nobin) then begin
if(n_elements(bin_lat) ne 0) then bin_lat=bin_lat else bin_lat=1.
if(n_elements(bin_lon) ne 0) then bin_lon=bin_lon else bin_lon=2.
endif
; Set Altitude
if(n_elements(alt) ne 0) then Alt=alt else alt=400.
;if(mltin eq 1) then begin
;for i=0,icnt-1 do begin
; malts(i)=alt
;opos=cnvcoord(mlats(i),mlons(i),malts(i))
;mlats(i) = opos(0)
;mlons(i) = opos(1)
;endfor
;endif
mglons=mlons
vcount=0
if mltin eq 1 then begin
print, 'AACGC ', mltin
isoy=0L
FOR I=0,ICNT-1 DO BEGIN
malts(i)=alt
opos=cnvcoord(mlats(i),mlons(i),malts(i))
mlats(i) = opos(0)
mlons(i) = opos(1)
mglons(i) = opos(1)
temp = mytimes.times(i)+(mytimes.doy-1)*24*3600
isoy = long(temp)
if( spole eq 0) then begin
mlons(i) = 180.0 + 15.*mlt(mytimes.year,isoy,mlons(i))
endif else begin
mlons(i) = 15.*mlt(mytimes.year,isoy,mlons(i))
endelse
if mlons(i) gt 180. then mlons(i)=mlons(i)-360.
; printf,4, mytimes.year,isoy,mytimes.times(i),mytimes.doy,mlons(i)
endfor
utalon = (mytimes.times(0)/3600.0)*15.0
zlon = mlons(0) - utalon
print, zlon,'0 UT',utalon
noonstr='0 UT'
endif
if mltin eq 0 then begin
; print, 'ECC ', mltin
sod=0L
for i=0,icnt-1 do begin
malts(i)=malts(i)+6371.2
sod = long(mytimes.times(i))
;printf, 4, mytimes.year,mytimes.doy,sod,malts(i),mlats(i),mlons(i),vest(i),$
; vnrt(i),qflgs(i)
opos = eccmlt(mytimes.year,mytimes.doy,sod,malts(i),mlats(i),mlons(i))
mglons(i) = opos(0)
mlats(i) = opos(1)
mlons(i) = opos(2)
if( spole eq 0) then begin
mlons(i)= 180.0 + mlons(i)*15.0
endif else begin
mlons(i)= mlons(i)*15.0
endelse
iF mlons(i) gt 180. then mlons(i)=mlons(i)-360.
endfor
utalon = (mytimes.times(0)/3600.0)*15.0
zlon = mlons(0) - utalon
; print, zlon,'0 UT',utalon
noonstr='0 UT'
endif
if(keyword_set(nob)) then begin
b_vest=vest
b_vnrt=vnrt
angs = b_vest
lats = mlats
lons = mlons
b_cont=vnrt
b_cont(*)=1
b_data=b_cont
b_data(*)=0
nvects = icnt
b_mlat=mlats
b_mlon=mglons
b_malt=malts
endif else begin
bin_lat = 1.
bin_lon = 2.
nbin_lat = abs(ylim-ymin)/bin_lat
nbin_lon = abs(xlim-xmin)/bin_lon
lats = fltarr(nbin_lat*nbin_lon)
lons = fltarr(nbin_lat*nbin_lon)
; print, ylim,ymin,xlim,xmin,nbin_lat,nbin_lon
; print, nop,mltin,qmin,qmax
for j = 0,nbin_lat-1 do begin
lat = float(j)*bin_lat + ymin + bin_lat/2.
for i=0,nbin_lon-1 do begin
lon = float(i)*bin_lon + bin_lon/2.
lats(i+j*nbin_lon) = lat
lons(i+j*nbin_lon) = lon - 180.
endfor
endfor
nvects = nbin_lat*nbin_lon
b_vest = fltarr(nvects)
b_vnrt = fltarr(nvects)
b_cont = lonarr(nvects)
b_data = lonarr(nvects)
b_mlat = fltarr(nvects)
b_mlon = fltarr(nvects)
b_malt = fltarr(nvects)
b_rot = fltarr(nvects)
;patch
; print, ymin, bin_lat, bin_lon, nbin_lon
for i = 0,icnt-1 do begin
lat_bin = abs(fix((mlats(i)-ymin)/bin_lat))
lon_bin = fix((mlons(i)+180.)/bin_lon)
indx = fix(lon_bin + lat_bin*nbin_lon)
if qflgs(i) ge qmin and qflgs(i) le qmax then begin
; printf,4, indx,i,b_vest(indx),vest(i)
b_vest(indx) = b_vest(indx) + vest(i)
b_vnrt(indx) = b_vnrt(indx) + vnrt(i)
b_cont(indx) = b_cont(indx) + 1
b_mlat(indx) = mlats(i)
b_mlon(indx) = mglons(i)
b_malt(indx) = malts(i)
endif
if qflgs(i) gt qmax then b_data(indx)=1
endfor
angs = fltarr(nvects)
w = where(b_cont ne 0,wc)
if( wc gt 0) then begin
b_vest(w) = b_vest(w)/b_cont(w)
b_vnrt(w) = b_vnrt(w)/b_cont(w)
endif
endelse
; Compute angles
w = where(b_vnrt ne 0,wc)
angs(w) = atan(b_vest(w),b_vnrt(w))
; Compute vector rotation angles for desired coordinate system
; print, 'Computing angle adjustment'
for i=0, nvects-1 do begin
if(b_cont(i) ne 0) then begin
b_rot(i)=angadj(mltin,mytimes.year,b_mlat(i),b_mlon(i),b_malt(i))
b_rot(i)=b_rot(i)*(3.141592/180.)
endif
endfor
w=where(b_rot ne 0, wc)
if(wc gt 0) then angs(w)=angs(w)-b_rot(w)
; Compute magnitudes
mags = sqrt(b_vest^2 + b_vnrt^2)
if keyword_set(nop) then begin
map_set,p0lat,p0lon,rot,limit=[ymin,xmin,ylim,xlim],/GRID,$
xmargin=xmargin,ymargin=ymargin,glinethick=0.5,/cyl,$
color=bkgrd,$
/noborder,latdel=latdel,londel=londel,_extra=extras
if(mltin eq 2) then begin
axis,xmin,ymin,xax=0,xrange=[xmin,xlim],xsty=1
axis,xmin,ymin,yax=0,yrange=[ymin,ylim],ysty=1
endif else begin
axis,xmin,ymin,xax=0,xrange=[0,24],xsty=1
axis,xmin,ymin,yax=0,yrange=[ymin,ylim],ysty=1
endelse
endif else begin
map_set,p0lat,p0lon,rot,limit=[ymin,xmin,ylim,xlim],/GRID,$
xmargin=xmargin,ymargin=ymargin,glinethick=0.5,/stereo,$
color=bkgrd,$
/noborder,latdel=latdel,londel=londel,_extra=extras
; MLT scale
for j=0,3 do begin
tstr = strtrim(string(6*j),2)
if( spole eq 0) then begin
lon=-180+j*90
xyouts,lon,ymin-2,tstr,alignment=0.5
endif else begin
lon=0-j*90
xyouts,lon,ylim+2,tstr,alignment=0.5
endelse
endfor
; Invariant Latitude scale
for j=0,3 do begin
lat= fix(ymin) + 10*j
tstr = strtrim(string(lat),2)
lon=zlon
xyouts,lon,lat,tstr,alignment=0.5
endfor
; UT scale
for j=0,7 do begin
utstr=strtrim(string(3*j),2)
lon=zlon+j*45.0
if( spole eq 0) then lat=78. else lat=-78.
xyouts,lon,lat,utstr,alignment=0.5
endfor
; Make tick marks on both UT and MLT scales
for j=0,23 do begin
lon=zlon+j*15.
olon=j*15
if( spole eq 0) then begin
plots,[lon,lon],[79.5,80.5]
plots,[olon,olon],[ymin,(ymin+0.5)]
endif else begin
plots,[lon,lon],[-79.5,-80.5]
plots,[olon,olon],[ylim,(ylim-0.5)]
endelse
endfor
; Axis
if( spole eq 0) then begin
utlat=replicate(80.0,360)
endif else begin
utlat=replicate(-80.0,360)
endelse
utlon=zlon+findgen(360.0)
oplot,utlon,utlat,line=0
; Axis
plots,[zlon,zlon],[ymin,ylim],thick=1.7
; Labels
xyouts,0.48,0.07,'MLT',/normal
xyouts,0.50,0.53,'UT',/normal
endelse
; vectors determined
S = findgen(16)*(!PI*2/16.)
usersym,cos(S),sin(S),/FILL
Re = 6362.
clon = (xlim+xmin)/2.
clat = (ylim+ymin)/2.
side_scale = 120./scale
cside = (90. - clat)*3.141592/180.
bside = scale/Re*side_scale
ang = 90.*3.141592/180.
arg = cos(bside)*cos(cside) + sin(bside)*sin(cside)*cos(ang)
if( arg lt -1)then arg=-1.
if(arg gt 1)then arg=1.
aside = acos(arg)
tlat = 90. - aside*180./3.141592
arg = (cos(bside)-cos(aside)*cos(cside))/(sin(aside)*sin(cside))
if( arg lt -1)then arg=-1.
if(arg gt 1)then arg=1.
bang = acos(arg)*180./3.141592
tlon = clon + bang
nx = convert_coord(clon,clat,/data,/to_normal)
ndx = convert_coord(tlon,tlat,/data,/to_normal)
delta = sqrt((ndx(0)-nx(0))^2+(ndx(1)-nx(1))^2)
side_scale = side_scale*.08/delta
cside = (90. - lats)*3.141592/180.
bside = mags/Re*side_scale
arg = cos(bside)*cos(cside) + sin(bside)*sin(cside)*cos(angs)
arg = arg < 1. > (-1.)
aside = acos(arg)
tlat = 90. - aside*180./3.141592
arg = (cos(bside)-cos(aside)*cos(cside))/(sin(aside)*sin(cside))
arg = arg < 1. > (-1.)
bang = acos(arg)*180./3.141592
q = where (angs lt 0,icnt)
if icnt ne 0 then bang(q) = -1.*bang(q)
tlon = lons + bang
for i = 0,nvects-1 do begin
if lons(i) gt xmin and lons(i) lt xlim and $
lats(i) gt ymin and lats(i) lt ylim and $
b_cont(i) ne 0 then begin
; ind = max(where(LVL le mags(i))) + 1
plots,lons(i),lats(i),psym=8,symsize=symsiz,color=symcol
plots,[lons(i),tlon(i)],[lats(i),tlat(i)]$
,thick=lthik
endif
if lons(i) gt xmin and lons(i) lt xlim and $
lats(i) gt ymin and lats(i) lt ylim and $
b_data(i) ne 0 then begin
plots,lons(i),lats(i),psym=3
endif
endfor
plots,.84,.92,psym=8,symsize=symsiz,/normal,color=symcol
plots,[.84,.92],[.92,.92],thick=lthik,/normal
scstring = strtrim(string(scale),2)+'m/s'
xyouts,.84,.88,scstring,/normal
if(thetitle eq ' ') then begin
if mltin eq 0 then begin
thetitle='Velocity vs Eccentric MLT/UT and Magnetic Latitude'
endif
if mltin eq 1 then begin
thetitle='Velocity vs AACGC MLT/UT and Magnetic Latitude'
endif
if mltin eq 2 then begin
thetitle='Geographic Velocity '
endif
endif
xyouts,0.15,0.95,thetitle,/normal
if(spole eq 1) then station=strmid(station,0,14)
xyouts,0.05,0.85,station,/normal
date = strcompress(string(fix(mytimes.mon))+'/'$'
+string(fix(mytimes.day))+$
'/'+string(fix(mytimes.year)),/remove_all)
xyouts,0.05, 0.9,date,/normal
time_string=systime()
disclaimer="Key Parameter and Survey data (labels K0,K1,K2) are preliminary data."
disclaimer1="Generated by CDAWeb on: "+time_string
xyouts,0.01,0.035,disclaimer,charsize=dtxt,/normal
xyouts,0.01,0.01,disclaimer1,charsize=dtxt,/normal
end