;$author: $ ;$Date: 2015/08/17 16:17:36 $ ;$Header: /home/cdaweb/dev/control/RCS/cdaweb_skymap.pro,v 1.9 2015/08/17 16:17:36 kovalick Exp kovalick $ ;$Locker: kovalick $ ;$Revision: 1.9 $ ;+------------------------------------------------------------------------ ; NAME: CDAWeb_Skymap ; PURPOSE: This Skymap code was given to us by the TWINS project in ; order to properly produce mapped image plots of their TWINS images. ; ; Modifications: TJK modified it slightly to allow a position to be ; defined, thus allowing it to be used from CDAWeb for making ; thumbnails and larger sized mapped images. I also changed the name ; of the routine from skymap to cdaweb_skymap, so there's no confusion ; between the two. ; ; CALLING SEQUENCE: (sample call from below, making large images) ; ; cdaweb_skymap, image, sc_pos = sc_posv_re_eci_dat[*,(frame-1)], $ ; spin_axis=spin_axis_eci_dat[*,(frame-1)], sun_pos=sun_posv_eci_dat[*,(frame-1)], $ ; prime_meridian=prime_meridian_eci_dat[*,(frame-1)], lonmin=-88, lonmax=268,latmin=4, $ ; latmax=88, colorbar=0, grid = 1, $ ; field=1, limb=1, sphere = 1, mag = mag, term=1, log=0, norotate=1, noerase=1,ctab=39, exobase=0 ; INPUTS: ; Required are: image, sc_pos, spin_axis, sun_pos, prime_meridian ; ; SEE DETAILS BELOW ; ; pro cdaweb_skymap, image, $ sc_pos = sc_pos, $ spin_axis = spin_axis, $ prime_meridian=prime_meridian, $ sun_pos = sun_pos, mag = mag, $ sphere = sphere, $ lonmin = lonmin, lonmax = lonmax, $ latmin = latmin, latmax = latmax, $ annotation = annotation, $ top_annotation = top_annotation, $ min = min, max = max, $ grid = grid, $ field = field, $ limb = limb, $ terminator = terminator, $ exobase = exobase, $ colorbar = colorbar, $ log = log, $ smooth = smooth, $ norotate = norotate, $ xmargin = xmargin, ymargin = ymargin, $ title = title, subtitle = subtitle, $ clip = clip, true = true, $ square = square, $ pathlogo = pathlogo, $ contour = contour, median = median, $ lee = lee, $ closeup = closeup, xcloseup = xcloseup, $ top = top, noerase = noerase, $ limit = limit, earth = earth, $ fullscreen = fullscreen, time = time, $ geogrid = geogrid, fill = fill, $ landcolor = landcolor, $ seacolor = seacolor, gei = gei, $ bartitle = bartitle, $ barlevels = barlevels, noimage = noimage, $ antiearthward = antiearthward, $ cbarposition = cbarposition, $ barcharsize = barcharsize, $ barformat = barformat,wim=wim, $ p0lat=p0lat,p0lon=p0lon, $ rot=rot, $ annosize=annosize, $ line_thick=line_thick, $ ps = ps, $ ctab=ctab, $ units = units, $ sun_spot = sun_spot, $ label_field_line = label_field_line, $ position = position, $; TJK add position for support of CDAWeb thumbnails thumb=thumb ; TJK add to indicate this is for a small thumbnail plot ;+ ; NAME: CDAWEB_SKYMAP.PRO ; ; ; PURPOSE: ; Project a lat-lon matrix onto the sphere of the sky and draw the ; earth with field lines. ; ; ; ; CATEGORY: ; ; ; CALLING SEQUENCE: ; CDAWEB_SKYMAP, IMAGE ; ; ; INPUTS: IMAGE[nlon,nlat] - An image with a column corresponding to ; one latitude and row corresponding to longtitude of the map. ; ; ; ; OPTIONAL INPUTS: ; ; ; ; KEYWORD PARAMETERS: ; ; SC_POS - A three element array of the S/C postition in some ; geocentric coordinate system. ; ; SPIN_AXIS - Spin axis direction of the satellite in the same ; geocentric coordinate system. ; ; PRIME_MERIDIAN - A 3 element vector that lies in the plane of the ; sphere pole and the zero reference angle of the ; sphere equator (the +x axis). This has nothing ; to do with the prime meridian of the earth; it is ; just the most convenient description. ; ; SUN_POS - The unit vector to the sun in the same system ; ; MAG - Magnetic dipole direction in the same geocentric system. ; ; All vectors and positions must be expressed in the same geocentric ; coordinate system. ; ; SPHERE - Set 0 gives us the full sphere. Set to 1 gives us the ; hemisphere lat and lon 0 deg in the center. ; ; LATMIN, LATMAX, LONMIN, LONMAX - Same as keywords to MAP_IMAGE.PRO ; ; ANNOTATION - An array of strings containing information such as time ; and position, to be output together with the plot in ; the assigned space. ; ; MIN - The minimum number on the colorbar. ; ; MAX - The max number of the color bar ; ; GRID - If set turns on lat-lon grid in S/C coordinates with the ; 'north pole' of the grid (+90 deg) being the spinaxis, and ; lon being the spin angle. ; ; FIELD - Plots, dipole field lines L=4 and 8 for noon, dusk, ; midnight, and dawn. ; ; LIMB - If set plots the Earth's limb. ; ; TERMINATOR - Turns on terminator line ; ; EXOBASE - If set plots the exobase line (dashed) at 300 km altitude. ; ; COLORBAR - If set plots a color bar. ; ; LOG - Turns on logarithmic scaling. ; ; SMOOTH - Applies a boxcar smoothing function to the mapped (output ; of MAP_IMAGE) image. The purpose is to smooth the visible ; image to aid the eye in interpreting the image as a 3D ; structure and not make the eye unconciously map it to the ; sphere. CAUTION!!: Localized features broadens. Care has to ; be taken when interpreting. ; ; NOROTATE - Turns off the fancy rotation that make some people ; seasick. If set, spin axis (lat = +90 deg) always points ; up (positive orbit normal down). ; ; XMARGIN, YMARGIN - Used in MAP_SET (Margin between map and image border) ; ; TITLE - Inserts title at top of SKYMAP Image ; ; SUBTITLE - Inserts subtitle under title at top of SKYMAP Image ; ; CLIP - All values GE than clip are set to 0. ; ; TRUE - Use TRUE Color representation ; ; SQUARE - Enable a square format window with minimum space for ; annotation. Default is landscape. The effect of this ; keyword is essentially to put the colorbar in the right ; place (vertically) and annotation under the map. XMARGIN ; and YMARGIN are used to set the position of the map. ; ; PATHLOGO - The absolut path to a GIF file containing an optional ; pathlogo to be placed. Does not work just yet. Seems that ; calling TVIMAGE the first time makes the position keyword ; not work. ; ; CONTOUR -'contourimage' is the image that is smoothed and median filtered. ; The CONTOUR procedure puts NLEVELS contours between min_value > ; min(data) and max_value < max(data). ; ; MEDIAN - Apply the MEDIAN filtering function with a ; nearest-neighbor window. This is applied to the raw ; pixelated image and not the mapped one. This can be used ; together with the SMOOTH keyword but not the LEE ; keyword. The purpose is to clean up statistical noise and ; outliers. CAUTION!!: This function alters the appearance ; of localized features such as the "low-altitude" ENA ; emissions in the IMAGE/HENA data. ; ; LEE - Apply the LEE filter to the raw image. This can be used ; together with the SMOOTH keyword but not the MEDIAN ; keyword. Should do a better job than MEDIAN but one needs to ; tweak the N and SIG arguments to the LEEFILT.PRO ; function. EXPERIMENTAL! ; ; CLOSEUP - A shorcut of setting the LIMIT keyword to MAP_SET. This ; gives you automatically 120 x 120 deg FOV. Equal to ; setting LIMIT=[60,-60,-60,60]. Overrides SPHERE. ; ; XCLOSEUP -A shorcut of setting the LIMIT keyword to MAP_SET. This ; gives you automatically 120 x 120 deg FOV. Equal to ; setting LIMIT= [45, -45, -45, 45]. Overrides SPHERE. ; ; TOP - Same as for BYTSCL.PRO. ; ; NOERASE - If set, skymap doesn't erase the previous plot. ; ; LIMIT - LAT, LON (in degrees) [LATMIN, LONMIN, LATMAX, LONMAX] ; ; EARTH - Adds the Earth with landcolor, seacolor, maps out continents ; ; FULLSCREEN - Eliminates colorbar, annotation. Entire screen filled ; by IMAGE ; ; TIME - Stucture with Year, Doy, Hour, Min, Sec Used with keyword EARTH ; ; GEOGRID - Used with keyword EARTH ; ; FILL - Used with keyword EARTH ; ; LANDCOLOR - Used with keyword EARTH, indicates color of continents ; ; SEACOLOR - Used with keyword EARTH, indicates color of water ; ; GEI - Indicates sc_pos (plus other vectors) is in GEI coordinates Used with keyword EARTH ; ; BARTITLE - Annotation to appear on colorbar, typically Units ; ; BARLEVELS - Number of divisions in the colorbar. (divisions + 1) annotations ; ; NOIMAGE - Indicates that no image should be made ; ; ANTIEARTHWARD - essentially sets ; p0lon=180 instead of 0 deg. ; ; CBARPOSITION - positions colorbar ; ; BARCHARSIZE - The character size of the color bar annotations. ; ; BARFORMAT - Format of the Numbers used in the colorbar labeling ; ; WIM - Passes back the map_image(image) to the calling program ; ; P0LAT - center of the plot being p0lat, p0lon in map coordinates ; ; P0LON - center of the plot being p0lat, p0lon in map coordinates ; ; ROT - If not set, then Rotate the map so that mag dipole axis always points ; upwards. If set use rot value to rotate the map. ; ; ANNOSIZE - Size of the text used in the annotation ; ; LINE_THICK - Thickness of lines used in images ; ; PS - Indicates whether we are making a postscript image or not ; **Might go away ; ; CTAB - Indicates which colortable is used to make the image ; ; UNITS - **Will most likely be replaced by bartitle ; ; SUN_SPOT - Adds the A, Dot annotation on the plots ; ; LABEL_FIELD_LINE - Adds 0, 6, 12, 18 labels to the field lines ; ; OUTPUTS: ; ; The input image projected onto the sphere of the sky, using the ; azimuthal projection with the 'north pole' of the sky sphere being ; the spin axis. The coordinats of the map are the spin based ; coordinates to represent what a spinning imager 'sees'. ; ; ; ; OPTIONAL OUTPUTS: ; ; ; ; COMMON BLOCKS: ; None! ; ; ; SIDE EFFECTS: ; ; KNOWN PROBLEMS/BUGS: - Field-line plotting crashes at sc_pos_x=0??? ; ; ; RESTRICTIONS: ; ; ; ; PROCEDURE: ; Set up the map by, setting 0 deg lat and lon as the ; subsatellite point of the skysphere, i.e. the center of ; the plot. The rotation (ROT argument to MAP_SET) is set ; such that MAG will always point upward on the plot ; window. Points a bit to the side for some positions but ; gives a continuous pointing position of the earth, so ; that it will work if one wants to do a movie. REFLECT_MAP ; is used to mirror turn the map since the map projections ; in IDL are used to put data on the surface of a ; sphere. We want to put data INSIDE a sphere. ; ; Field lines are plotted for fixed L-shell and MLT range. ; ; DEPENDENCIES (non-standard IDL routines): ; ; FIELD_LINE.PRO - To plot field_lines. ; EARTH_BOUNDARY.PRO - To plot Earth's limb. ; TERMINATOR_LINE.PRO - To plot terminator line. ; SPHERE_BASIS.PRO ; DOT.PRO ; SPHERE2SPHERE.PRO ; SPHERE_TO_XYZ.PRO ; XYZ_TO_SPHERE.PRO ; REFLECT_MAP.PRO ; XYTEXT.PRO ; COLORBAR.PRO ; ; EXAMPLE: ; Following examples lets you view the earth as the ; satellite passes it in a straight trajectory below the ; southpole. ; ;pro test ; mag = [0, 0, 1] ; sun_pos = [1, 0, 0] ; spin_axis = [-1, 1, 0] ; xs = transpose(interpol([6, -6], 50)) ; ys = xs ; zs = transpose(replicate(-2, 50)) ; sc_pos = [[xs, ys, zs]] ; for i = 0, 49 do begin ; SKYMAP, dist(10), sc_pos = reform(sc_pos(*, i)), spin_axis = spin_axis, $ ; sun = sun_pos, sphere = 0, mag = mag ; wait, .1 ; endfor ;end ; ; ; MODIFICATION HISTORY: ; ; Jun 19 2007, Jillian Redfern, SWRI ; - received new skymap.pro from Pontus Brandt ; incorporating his changes into our version of skymap ; Added in no_sun_spot, ps, line_thick, annosize, ctab keywords ; Added additional comments to the top of the program ; ; Mar 22, 2005. Michael Muller, SwRI ; - Commented out overplot of gridlines. ; ; Dec 06, 2004. Michael Muller, SwRI ; - originally created/tested Nov 07, 2002. ; - Added code to make map of pixels ; so can get values at any x,y grid ; point. Also, overlay/plot the grid. ; ; Thu Mar 4 09:32:56 2004, Pontus Brandt ; ; ; - Added keyword PRIME_MERIDIAN for plotting the limb. ; ; Thu Sep 25 15:18:24 2003, Pontus Brandt ; ; ; - Bug: LATMIN=0 (and some other keywords too) caused ; latmin to be set to -90, since I used KEYWORD_SET ; routine. Changed to check set keyowrds with ; N_ELEMENTS instead. ; ; Thu Sep 25 15:11:21 2003, Pontus Brandt ; ; ; - Added keywords P0LAT and P0LON to generalize type of display for ; other instrumental geometries. ; ; Tue Apr 17 14:21:16 2001, Pontus Brandt ; ; ; - Added MLTs to the tip of each field line ; plotted. Removed "SUN" annotation. ; ; Wed Mar 14 18:39:08 2001, Pontus Brandt ; ; ; - Moved call to COLORBAR last and instead saved system ; variables from the map coordinate system just after ; having created it. Restore them again just before ; exiting this procedure. ; ; Wed Mar 14 14:40:12 2001, Pontus Brandt ; ; ; - Added keyword CBARPOSITION to have an option to ; position the colorbar. ; ; Fri Jan 12 12:55:15 2001, Pontus Brandt ; ; ; - Added keyword ANTIEARTHWARD, which essentially sets ; p0lon=180 instead of 0 deg. Should carry with it the ; right LIMIT, but somehow it looks weird for CLOSEUP ; and XCLOSEUP. Hemisphere and full sphere seems to be ; fine. ; ; Fri Jan 12 10:23:02 2001, Pontus Brandt ; ; ; - Added keyword ; EARTH: gives outline of earth and continents. ; COLOREARTH: gives filled continents. ; Positioning and sizing of the earth is NOT solved ; yet, in other words, it doesnt work, but have to set ; manually in the program! It is EXPERIMENTAL. ; ; Wed Dec 6 12:48:36 2000, Pontus Brandt ; ; ; - Adding keyword LIMIT in order to bring the setting ; of CLOSEUP and XCLOSEUP external to SKYMAP. Keeping ; old keywords CLOSEUP, XCLOSEUP, SPHERE just in case ; someone needs to use them. ; ; Wed Oct 11 08:45:47 2000, Pontus Brandt ; ; ; - Changed the plotting order to be able to save the ; data coordinate system established by MAP_SET. Put ; all plotting not using the data system, first. Still ; problem with CONTOUR keyword. It is using the data ; system and shouldnt destroy it, but some how it gets ; meesed up. ; ; Thu Sep 28 13:42:27 2000, Pontus Brandt ; ; ; - Added NOERASE keyword. The erasing should be done in ; here, if any. ; ; Wed Sep 27 14:05:57 2000, Pontus Brandt ; ; ; - Added TOP keyword to be able to use only a specific ; region of the colortable. ; ; Mon Sep 25 10:01:29 2000, Pontus Brandt ; ; ; - Made filtering possible for 24-bit images. ; ; Wed Sep 20 07:19:20 2000, Pontus Brandt ; ; ; - Added check if field lines cannot be drawn. ; ; Tue Sep 19 10:10:30 2000, Pontus Brandt ; ; ; - Added MEDIAN and LEE keywords. Need to tweak their ; parameters a bit, especially LEE. ; - Also, set BILINEAR=1 if SMOOTH is set. This soften ; the edges in the MAP_IMAGE function. ; ; Wed Sep 6 12:01:07 2000, Pontus Brandt ; ; ; - Added clip keyword. All values GE than clip ; are set to 0. ; ; Wed Sep 6 10:00:31 2000, Pontus Brandt ; ; ; - Fixed format codes for color bar and logscale. ; ; Thu Jun 15 10:03:09 2000, Pontus Brandt ; ; ; - Added keywords XMARGIN, YMARGIN, and TITLE. ; ; Thu May 18 19:06:34 2000, Pontus Brandt ; ; ; - Added some keywords for logscale and general looks. ; ; Fri May 12 18:28:56 2000, Pontus Brandt ; ; ; - Added keyword MAX for the colorbar. ; ; Tue May 9 11:10:50 2000, Pontus Brandt ; ; ; - Plotting position of Sun and magnetic north pole. ; ; Mon May 8 18:07:00 2000, Pontus Brandt ; ; ; - Added keyword ANNOTATION. ; ; Sat Apr 15 16:58:21 2000, Pontus Brandt ; ; ; - Included keywords latmin, latmax, lonmin, lonmax. ; ; Fri Apr 14 13:14:47 2000, Pontus Brandt ; ; ; - First version done with spin based map ; ; coordinates. Next step is to have an option for ; having ; the map coordinates in a geocentric ; coordinate system ; such as Solar Magnetic ; coordinates. ; ;- DEBUG=0 if not keyword_set(sc_pos) then sc_pos = [0, 3, 0] ;TWINS Specific if not keyword_set(spin_axis) then spin_axis = [0, -1, 0] if not keyword_set(prime_meridian) then prime_meridian = -sc_pos if not keyword_set(sun_pos) then sun_pos = [1, 0, 0] if not keyword_set(mag) then mag = [0, 0, 1] if not keyword_set(sphere) then sphere = 0 if not keyword_set(lonmin) then lonmin = -90 ;TWINS Specific if not keyword_set(lonmax) then lonmax = 270 ;TWINS Specific if not n_elements(latmin) then latmin = 2 ;TWINS Specific if not n_elements(latmax) then latmax = 90 ;TWINS Specific if not keyword_set(min) then min = min(image) if not keyword_set(max) then max = max(image) if keyword_set(smooth) then bilinear = 1 else bilinear = 0 if not keyword_set(xmargin) then xmargin = 1 if not keyword_set(ymargin) then ymargin = [10,1] if not keyword_set(title) then title = '' if not keyword_set(subtitle) then subtitle = '' if not keyword_set(clip) then clip = max(image)+10. if not keyword_set(true) then true = 0 if not keyword_set(square) then square = 0 if not keyword_set(contour) then contour = 0 if n_elements(top) eq 0 then top = !d.table_size-4 if not keyword_set(noerase) then erase if keyword_set(fullscreen) then begin annotation = 0 top_annotation = 0 colorbar = 0 xmargin = 0 ymargin = 0 endif if (not(keyword_set(bartitle)) and not(keyword_set(units))) then bartitle = 'FLUX (cm!E2!N sr s)!E-1!N' if keyword_set(units) and not(keyword_set(bartitle)) then begin if keyword_set(log) then begin extra = 'LOG(' extra_end = ') ' endif else begin extra = '' extra_end = ' ' endelse case units of 1: begin bartitle = 'FLUX ' + $ extra + '(cm!E2!N sr s)!E-1!N' + $ extra_end ;integral number flux end 2: begin bartitle = 'E FLUX ' + $ extra + '(eV cm!E2!N sr s)!E-1!N' + $ extra_end ; differential number flux end 4: begin bartitle = extra +'DEGREES FROM SUN' + $ extra_end end 5: begin bartitle =extra + 'COUNTS' + $ extra_end end 10: begin bartitle = '!3Fractional Error !7D!3F/F' end else: begin bartitle = 'E FLUX ' + $ extra + '(eV cm!E2!N sr s)!E-1!N' + $ extra_end ; differential number flux end endcase endif if not keyword_set(barlevels) then barlevels = 2 if not keyword_set(p0lat) then p0lat = 90 ;TWINS Specific if not keyword_set(p0lon) then p0lon = 0 ;TWINS Specific if not keyword_set(line_thick) then line_thick = 2.0 if not keyword_set(ps) then ps = 0 if not keyword_set(ctab) then ctab = 0 ;grayscale if not keyword_set(annosize) then annosize = 1.0 ;TJK need to catch when sc_pos values aren't defined - plot w/o field ;lines,etc. as per Jillian's instructions if (sc_pos[0] eq -9.9999998e+30) then begin print, 'WARNING=CDAWeb_skymap sc_pos is not set, setting field, limb and term to 0' field = 0 limb = 0 terminator = 0 endif if n_elements(sc_pos) ne 3 then begin print, '%SKY.PRO: SC_POS must have three elements' return endif if sqrt(total(sc_pos^2)) le 1 then begin print, '%SKY.PRO: SC_POS is below the surface' return endif if n_params() lt 1 then begin print, 'Must have 1 argument' return endif s = size(image) if s[0] lt 2 or s[0] gt 3 then begin print, '%SKYMAP:Argument must have 2 or 3 dimensions' return endif ;; Check for mismatching keywords if keyword_set(lee) and keyword_set(median) then begin print, '%SKYMAP.PRO: Cannot do Lee and Median filtering together.' return endif if keyword_set(contour) and keyword_set(true) then begin print, '%SKYMAP.PRO: Keywords CONTOUR and TRUE not allowed...yet' contour = 0 endif ;; Annotation ;for CDAWeb, want background to be black which is 0 w/ C.T. 13 ;for CDAWeb, want foreground to be white which is 255 w/ C.T. 13 !p.color = 255 !p.background = 0 if keyword_set(annotation) then begin cs = PlotCharSize() if (ps EQ 1) then tvimage, (intarr(255,255)+1)*255 if (ps EQ 1) then !p.color = 0 if (ps EQ 1) then !p.background = 255 if keyword_set(square) then begin xyouts, .02, .95, strcompress(string("!5"+title)), $ align = 0, /normal, charsize = 1.5, color=!p.color xyouts, .02, .9, strcompress(string("!5"+subtitle)) $ , align = 0, $ /normal, charsize = 1, color=!p.color endif else begin xytext, .05, cs.y*annosize*6.5, annotation, /normal, spacing = -1.2, charsize = $ annosize, charthick = line_thick, color=!p.color ; xytext, .6, .7, annotation, /normal, spacing = -1.5, charsize = 1.2 xyouts, .28, .9, string("!5"+title), align = .5, /normal, charsize = annosize, color=!p.color xyouts, .28, .86, string("!5"+subtitle), align = .5, $ /normal, charsize = annosize*0.75, color=!p.color endelse endif ;Top Annotation if keyword_set(top_annotation) then begin xyouts, .985, .05, top_annotation, /normal, charsize = $ annosize*0.75, color=!p.color, orientation=90 endif ;; Set up the map ;; Rotate the map so that mag dipole axis always points ;; upwards. Sign of rotation is given by vector product between spin ;; and mag. If the vector product points towards the observer (the ;; same hemisphere) the longitudes go upward in the same direction ;; as the dipole, but if it is negative the longitudes go down and ;; so does the dipole. dd = sqrt(total(prime_meridian^2)*total(mag^2)) if n_elements(rot) eq 0 then begin if dd ne 0 then $ rot = acos(total(prime_meridian*mag)/dd)/!dtor $ else begin print, '%SKYMAP.PRO: WARNING: Spinaxis or Mag = 0!' rot = 0 endelse v = crossp(prime_meridian, mag) if total(v^2) ne 0 then $ rot = total(v*(-sc_pos))/ $ sqrt(total(v^2)*total(sc_pos^2))*rot $ else rot = rot endif ;If we don't want to rotate, set rot to be 0 if keyword_set(norotate) then rot = 0 if keyword_set(limit) then $ limit = limit else $ if keyword_set(xcloseup) then $ limit = [0, -45, 45, 45] else $ if keyword_set(closeup) then $ limit = [0, -60, 60, 60] else $ if keyword_set(sphere) then $ limit = [0, -180, 90, 180] $ ;; Hemisphere ;TWINS Specific else $ limit = [0, -180, 90, 180] ;; Hemisphere ;TWINS Specific ;; Preliminary for a cartwheel orbit with the center of the plot being ;; nadir, i.e. lat, lon = 0, 0 in map coordinates ;; (DEFINITION)...unless ANTIEARTHWARD keyword is set. if n_elements(p0lat) eq 0 and n_elements(p0lon) eq 0 then begin if keyword_set(antiearthward) then begin p0lat = 0 p0lon = 180 limit = limit*[1, -1, 1, -1] ;; reverse longitudes endif else begin p0lat = 0 p0lon = 0 endelse endif if n_elements(p0lat) ne 0 and n_elements(p0lon) eq 0 then begin p0lon=0 endif ;TJK 10/15/2009 - added this code to accept a position w/in a window ; so that we can make thumbnails... hopefully if keyword_set(thumb) then begin ;Have to switch to normalized coordinates - that's what this code uses ;vs. device x_loc = [position[0],position[2]] y_loc = [position[1],position[3]] n_loc = convert_coord(x_loc,y_loc,/device, /to_normal) ;convert_cord returns, x,y,z, x2,y2,z2 - we're ignoring the z's position=[n_loc[0],n_loc[1],n_loc[3],n_loc[4]] ; print, 'position values in normalized coord',n_loc ; print, 'the following values should be between 0 and 1',position map_set, p0lat, p0lon, rot, /azimuth, /iso, $ limit = limit, /noborder, /noerase,$ position=position ;TJK added used of position for thumbnails ;TJK removed xmargin = xmargin, ymargin = ymargin,$ endif else begin ;regular large plots ; print, 'thumb is NOT set in cdaweb_skymap' ;TJK have to clear out the !p.position and !p.region values, otherwise ;the position/size of the plot carries over from the previous ;requested plot, e.g. orbit plot... !p.position=[0,0,0,0] !p.region=[0,0,0,0] ;print, 'clear out !p.position = ',!p.position ;print, 'clear out !p.region = ',!p.region ;print, 'Above map_set call in cdaweb_skymap' ;print, 'limit = ',limit ;print, 'xmargin ',xmargin ;print, 'ymargin ',ymargin map_set, p0lat, p0lon, rot, /azimuth, /iso, $ limit = limit, xmargin = xmargin, ymargin = ymargin, $ /noborder, /noerase endelse reflect_map box = 10 ;; Filter the raw image if keyword_set(median) then begin if keyword_set(true) then begin if (n_elements(image[0, *, 0]) mod 2 eq 0) and $ (n_elements(image[0, 0, *]) mod 2 eq 0) then even = 1 image[0, *, *] = median(reform(image[0, *, *]), 3, even = even) image[1, *, *] = median(reform(image[1, *, *]), 3, even = even) image[2, *, *] = median(reform(image[2, *, *]), 3, even = even) endif else begin if (n_elements(image[*, 0]) mod 2 eq 0) and $ (n_elements(image[0, *]) mod 2 eq 0) then even = 1 image = median(image, 3, even = even) endelse endif if keyword_set(lee) then begin n = 2 sig = .4 if keyword_set(true) then begin image[0, *, *] = leefilt(reform(image[0, *, *]), n, sig) image[1, *, *] = leefilt(reform(image[1, *, *]), n, sig) image[2, *, *] = leefilt(reform(image[2, *, *]), n, sig) endif else $ image = leefilt(image, n, sig) endif ;; Warp and display image ;; For 24-bit images if keyword_set(true) then begin r = map_image(reform(image[0, *, *]), sx, sy, xsize, ysize, $ latmin = latmin, latmax = latmax, $ lonmin = lonmin, lonmax = lonmax, /compress) g = map_image(reform(image[1, *, *]), sx, sy, xsize, ysize, $ latmin = latmin, latmax = latmax, $ lonmin = lonmin, lonmax = lonmax, /compress) b = map_image(reform(image[2, *, *]), sx, sy, xsize, ysize, $ latmin = latmin, latmax = latmax, $ lonmin = lonmin, lonmax = lonmax, /compress) wim = bytarr(3, n_elements(r[*, 0]), n_elements(r[0, *]) ) wim[0, *, *] = r wim[1, *, *] = g wim[2, *, *] = b endif else begin ;; Save the image to contour for later if (n_elements(image[*, 0]) mod 2 eq 0) and $ (n_elements(image[0, *]) mod 2 eq 0) then even = 1 ;print, 'sx, sy and xsize and ysize are undefined before this call' ;print, 'above map_image calls latmin, latmax = ',latmin, latmax ;print, 'above map_image calls lonmin, lonmax = ',lonmin, lonmax contourimage = median(image, 3, even = even) contourimage = map_image(contourimage, sx, sy, xsize, ysize, $ latmin = latmin, latmax = latmax, $ lonmin = lonmin, lonmax = lonmax, /compress, $ /bilinear) ;TJK changed from edge to edge_truncate because IDL8.1 offers 3 edge keywords ; contourimage = smooth(contourimage, box, /edge, /nan) contourimage = smooth(contourimage, box, /edge_truncate, /nan) if keyword_set(log) then contourimage = alog10(contourimage+1e-6) ;print,'After return from 1st map_image, sx and sy are ', sx, sy ;print, 'above map_image calls xsize and ysize = ',xsize, ysize ;print, 'above 2nd map_image calls latmin, latmax = ',latmin, latmax ;print, 'above 2nd map_image calls lonmin, lonmax = ',lonmin, lonmax wim = map_image(image, sx, sy, xsize, ysize, $ latmin = latmin, latmax = latmax, $ lonmin = lonmin, lonmax = lonmax, /compress, $ bilinear = bilinear, missing = 0) ;!p.background) ;print,'After return from 2nd map_image, sx and sy are ', sx, sy ;print, 'and xsize and ysize = ',xsize, ysize endelse if min(image) lt 0 then print, '% SKYMAP: WARNING: Found pixels < 0.' ;; Top clip image ind = where(wim ge clip, cts) if cts gt 0 then begin ;print, 'cliping/setting ',cts,' values to 0' wim[ind] = 0 endif if keyword_set(log) then begin enamax = alog10(max > 1e-1) enamin = alog10(min > 1e-2) if keyword_set(units) then begin if (units eq 2) then begin enamax = alog10(max > 1e-3) enamin = alog10(min > 1e-4) endif endif wim = alog10(wim+0.000001) endif else begin enamin=min enamax=max endelse if enamin ge enamax then begin if (keyword_set(log)) then $ enamin = enamax-1.5 $ else begin enatmp = enamin enamin = enamax enamax = enatmp if (enamin eq enamax) then begin enamin=0 enamax = ((enamin+.001) > enamax) endif endelse endif ;bad = where(wim le 0b,cts) ;if cts gt 0 then wim(bad) = !p.background if keyword_set(smooth) then begin if smooth gt 1 then box = smooth $ else box = 10 ;TJK changed from edge to edge_truncate because IDL8.1 offers 3 edge keywords ; wim = smooth(TEMPORARY(wim), box, /edge, /nan) wim = smooth(TEMPORARY(wim), box, /edge_truncate, /nan) endif ;; We bytscale WIM here, so that the IMAGE matrix always contains the ;; data values byteimage = bytscl(TEMPORARY(wim), min = enamin, max = enamax, top = top) ind = where(byteimage eq 0B, cts) if cts gt 0 then byteimage[ind] = !p.background if not keyword_set(noimage) then $ tv, byteimage, sx, sy, $ xsize = xsize, ysize = ysize, true = true if keyword_set(grid) then map_grid, /label, glinethick = 1, $ charsize=annosize*0.75, glinestyle = 1 ;; Save map coordinate system here. Restore it before exiting this ;; procedure. Then we can plot whatever we want (colorbars, ;; contours, etc.) draw_st = {p:!p, x:!x, y:!y, z:!z, map:!map} if keyword_set(limb) then begin ; Earth boundary with spin axis in geocentric coordinates pointing to ; the +z direction in viewing sphere. eb = earth_boundary(sc_pos, spin_axis, $ prime_meridian = prime_meridian, /deg, $ earth = earth) plots, eb[2, *], eb[1, *], thick = line_thick, /data;, linestyle=1 endif if keyword_set(exobase) then begin ; Exosphere base at 300km eb = earth_boundary(sc_pos, spin_axis, $ prime_meridian = prime_meridian, /deg, $ radius = 1.+300./6378.) plots, eb[2, *], eb[1, *] ;, linestyle = 5 endif if keyword_set(terminator) then begin ; Terminator w = terminator_line({c:{axis:spin_axis, phase:prime_meridian, sun:sun_pos}, $ p:sc_pos*6378.}, num = 10) if n_elements(w) ge 2 then begin if annosize ge 4. then $ plots, w[2, *], w[1, *], $ noclip = 0, psym=2 $ else $ plots, w[2, *], w[1, *], $ noclip = 0, linestyle=1 endif endif if keyword_set(sun_spot) then begin if keyword_set(annotation) then begin ;; Get Sun and Mag Pole in instrumental coordinates ;; First, define the axis of the instrumental cooridnate system in ;; whatever geocentric system you've defined spin_axis, sun, sc_pos, ;; and mag in. z = spin_axis/norm(spin_axis) x = prime_meridian x = x/norm(x) y = crossp(z, x) ;; And, second, convert sun and mag to the instrumental coordinate system. p = sphere2sphere(xyz_to_sphere([[sun_pos], [mag]], /degree), $ diag([1, 1, 1]), [[x], [y], [z]], /degree) npts = 20 x = cos(interpol([0, 2*!pi], npts)) y = sin(interpol([0, 2*!pi], npts)) usersym, x, y, /fill plots, p[2, 0], p[1, 0], psym = 8, /data, symsize = 2 xyouts, (p[2, 0]+180.) mod 360., -p[1, 0], $ align = .5, 'A', charsize = 1 endif ;annotation endif ;sunspot ;NO WORKING, SWRI Doesn't have a copy of all needed programs ; if keyword_set(earth) then begin ; ;; Calculate angular extent of earth. Ofcourse the earth is NOT ; ;; warped onto the existing map. We assume the region where earth is ; ;; plotted is rectangular enough for it to look real. in other ; ;; words, if the earth is small enough it can be approximated by a ; ;; circle. ; alpha = atan(1., sqrt(total(sc_pos^2)))/!dtor ; x0 = alpha ; y0 = x0 ; x1 = -alpha ; y1 = x1 ; position = [x0, y0, x1, y1] ; ;; p0lat and p0lon is the footprint on the earth in GEOGRAPHIC ; ;; coordinates. Thus, take the SM lat and lon and convert to GEO. ; if not keyword_set(gei) then begin ; pos = gei2sm(sc_pos, time.year, time.doy, time.hour, $ ; time.min, time.sec, /inv) ; endif else pos = sc_pos ; date = doy2date(time.doy, time.year, /num) ; posgeo = gei2geo(pos, $ ; [time.hour+time.min/60., date[2], date[1], time.year]) ; posgeo_sphere = xyz2sphere(posgeo, /deg) ; p0lat = posgeo_sphere(1) ; p0lon = posgeo_sphere(2) ; ;; Calculate the rotation by first calculate the rotation to line ; ;; up with spin axis. The vgeo=Geonorth-(p0lat,p0lon)_geo vector is ; ;; pointing down on the screen for rot=0. ; ;; Rotation between spinaxis and vgeo ; if not keyword_set(gei) then begin ; spin = gei2sm(spin_axis, time.year, time.doy, time.hour, $ ; time.min, time.sec, /inv) ; endif else spin = spin_axis ; spingeo = gei2geo(spin, $ ; [time.hour+time.min/60., date[2], date[1], time.year]) ; spingeo_sphere = xyz2sphere(spingeo, /deg) ; rot = spingeo_sphere(2)-p0lon ; rot = 172. ; earth, rot = rot, p0lat = p0lat, p0lon = p0lon, position = position, $ ; data=1, geogrid = geogrid, fill = fill, landcolor = landcolor, $ ; seacolor = seacolor ;, /pole, /foot ; endif ;TJK add code to save off the colortable before these fieldlines and ;others are drawn and then restore them at the bottom. Using a ;different colortable than original, so purple is a different #. tvlct, orig_r,orig_g,orig_b, /get if keyword_set(field) then begin ;; Field lines ; Red(color=254)=noon=12LocalTime. Yellow(90)=dusk=18LT. st = {gei: sc_pos*6378., mag_gei: mag, sun_gei: sun_pos, $ axis_gei: spin_axis, prime_gei: prime_meridian} l = [4, 8] flon = [0, 90, 180, 270] ; purple = 253 purple = 203 ; in CDAWeb using yellow because purple clashes w/ data ; tvlct, 180, 120, 200, purple tvlct, 180, 120, 213, purple ; TJK guessing 200 needs to be raised for i = 0, n_elements(l)-1 do $ for j = 0, n_elements(flon)-1 do begin if flon[j] eq 0 then ccolor = 254 $ ;254 is red in CDAWeb else if flon[j] eq 90 then ccolor = purple else ccolor=!p.color line = field_line(st, l[i], flon[j]) if size(line, /n_dim) eq 2 then begin plots, line[2, *], line[1, *], thick = line_thick, color=ccolor ;if (flon[j] eq 90 ) then plots, line(2, *), line(1, *), linestyle=1, thick = line_thick, color= !p.color ;if (flon[j] eq 0 ) then plots, line(2, *), line(1, *), linestyle=1, thick = line_thick, color= !p.color endif else begin print, '%SKYMAP: Cannot draw field line. Function returns', line endelse endfor if keyword_set(label_field_line) then begin ;; Put MLT on the tips of each field line. mltname = ['12', '18', '00', '06'] for i = 0, n_elements(flon)-1 do begin line = field_line(st, l[n_elements(l)-1], flon[i]) if size(line, /n_dim) eq 2 then begin r = line[2, *]^2+line[1, *]^2 ind = where(r eq max(r)) phi = line[2, ind[0]] phi = 1.1*phi theta = line[1, ind[0]] theta = 1.1*theta xyouts, phi, theta, mltname[i], $ charsize = annosize endif endfor endif endif ; ;mam Jun21,2005. ; ;This lanl field lines/satellite ; ;section is not used (lanl=0 is set). ; ;* It also has a bug! * ; ;The part: ; ; lanl_mltlon = lanl_lon ; ; for kk = 0, n_elements(lanl_lon) - 1 do begin ; ; if (lanl_lon[kk] le 0) then $ ; ; lanl_mltlon[kk] = ut - lanl_lon[kk]/15. $ ; ; else $ ; ; lanl_mltlon[kk] = ut - lanl_lon[kk]/15. ; ; lanl_mltlon[kk] = lanl_mltlon[kk] * 15. ; ; endfor ; ;can be rewritten lanl_mltlon = ut*15 - lanl_lon ; ;It is unknown what the exact equation is or ; ;if the intent was to use the if() and then ; ;have one branch a + and the other a -. ; ; ; lanl = 0 ; if (lanl) then begin ; ;; Field lines ; ;if (!d.name eq 'PS') then loadct,0,/silent ; st = {gei: sc_pos*6378., mag_gei: mag, sun_gei: sun, $ ; axis_gei: spin_axis, prime_gei: -sc_pos} ; l = [6.6] ; lanl_lon = [-165.1,8.5,69.5,-38.1,103.4] ; lanl_mltlon = lanl_lon ; for kk = 0, n_elements(lanl_lon) - 1 do begin ; if (lanl_lon[kk] le 0) then $ ; lanl_mltlon[kk] = ut - lanl_lon[kk]/15. $ ; else $ ; lanl_mltlon[kk] = ut - lanl_lon[kk]/15. ; lanl_mltlon[kk] = lanl_mltlon[kk] * 15. ; endfor ; for i = 0, n_elements(l)-1 do $ ; for j = 0, n_elements(lanl_mltlon)-1 do begin ; line = field_line(st, l(i), lanl_mltlon(j)) ; plots, line(2, *), line(1, *), thick = line_thick, color=254 ; endfor ; endif if keyword_set(contour) then begin ;; Contour position = [sx, sy, sx+xsize, sy+ysize] ;; 'contourimage' is the image that is smoothed and median filtered. ;; The CONTOUR procedure puts NLEVELS contours between min_value > ;; min(data) and max_value < max(data). contour, contourimage, /follow, nlevels = 9, $ /noerase, xstyle = -1, ystyle = -1, position = position, dev=1, $ c_charsize = 1.8, min_value = enamin, max_value = enamax, c_thick = 2 endif if keyword_set(colorbar) then begin ;; Color bar if enamin ge enamax then begin print, '% SKYMAP: Keyword MIN must be less than MAX!' enamin = 0 endif if keyword_set(square) then begin position = [0.92, 0.25, 0.95, 0.75] vertical = 1 endif else begin ;; Default position = [0.92, 0.05, 0.95, 0.35] vertical = 1 endelse if keyword_set(cbarposition) then $ position = cbarposition if not keyword_set(barcharsize) then barcharsize = 1. if not keyword_set(barformat) then begin barformat='(f5.2)' if (enamax ge .001 and enamax lt 1) then barformat='(f6.3)' if (enamax ge 1 and enamax lt 10) then barformat='(f6.3)' if (enamax ge 10 and enamax lt 100) then barformat='(f6.2)' if (enamax ge 100 and enamax lt 1000) then barformat='(f7.2)' if (enamax ge 1000 and enamax lt 10000) then barformat='(f6.1)' if (enamax ge 10000 and enamax lt 100000) then barformat='(f7.1)' if (enamax ge 100000 and enamax lt 1000000) then barformat='(f8.1)' if (enamin ge .0001) then begin if (enamin ge .001 and enamin lt .01 ) then barformat='(f5.3)' if (enamin ge .0001 and enamin lt .001) then barformat='(f6.4)' endif if (enamax ge 1000000) then barformat='(e7.1)' endif if (DEBUB) then print,' ***### min,max',min,max twinscolorbar, position = position, vertical = vertical, $ minrange = enamin, maxrange = enamax, format = barformat, title = bartitle, $ ncolors = top, ps = (!d.name eq 'PS'), color = !p.color, $ divisions = barlevels > 4, charsize = annosize*.7, minor=0 endif ;; Restore coordinate system variables !p = draw_st.p !x = draw_st.x !y = draw_st.y !z = draw_st.z !map = draw_st.map ;TJK restore the color table prior to the use of purple tvlct, orig_r,orig_g,orig_b return end ;pro skymap