;+
; NAME: 
;	MIGRF_MAP 
;
; PURPOSE:
;	To display maps of the X, Y, Z, H, I, and D components of the
;	earth's magnetic field as given by the International Geomagnetic
;	Reference Field.
;
; CATEGORY:
;	Geophysical models.
;
; CALLING SEQUENCE:
;	MIGRF_MAP, year
;
; INPUTS:
;	YEAR: The year for which the field is to be evaluated.
;             (Last two digits).
; KEYWORDS:
;	LIMIT:	A four-element vector = [latmin, lonmin, latmax, lonmax],
;		specifying the boundaries of the region to be mapped.
;		(See also the MAP_SET procedure).
;
; OUTPUTS:
;	A map of the desired component, chosen from a menu.
;
; COMMON BLOCKS:
;	None.	
;
; SIDE EFFECTS:
;	None.
;
; RESTRICTIONS:
;	None.
;
; MODIFICATION HISTORY:
;	September 14, 1993. Written by A. Louro.
;-

pro map_plot, var, x, y, levels, title, limit = limit, equ = equ	; Map plotting routine.
tek_color
erase
if keyword_set(limit) then begin
  p0lat = (limit(2) + limit(0)) / 2.
  p0lon = (limit(3) + limit(1)) / 2.
  map_set, p0lat, p0lon, /grid, /continents, /orthographic, color = 6, title = title, $
    limit = limit, /label, sat_p = [1., 0., 0.], mlinethick = 4.
endif else begin
  map_set, /grid, /continents, /mercator, color = 6, title = title, $
    /label, mlinethick = 4.
endelse
contour, var, x, y, /overplot, levels = levels, /follow, $
  color = 7, charsize = 1.5
if keyword_set(equ) then contour, var, x, y, /overplot, levels = [0.], $
  c_thick = 2., c_annotation = 'dip equator', color = 1, charsize = 1.5
return
end

pro migrf_map, year, limit = limit

if not(keyword_set(limit)) then begin
  lat = findgen(91) * 2. - 86		; Construct latitude and
  lon = findgen(91) * 4. - 180.		; and longitude arrays.
  n_lat = 91
  n_lon = 91
endif else begin
  read, 'Minimum latitude: ', lat_min
  read, 'Maximum latitude: ', lat_max
  read, 'Minimym longitude: ', lon_min
  read, 'Maximum longitude: ', lon_max
  limit = [lat_min, lon_min, lat_max, lon_max]
  lat_step = min([(lat_max - lat_min) / 10., 2.])
  lon_step = min([(lon_max - lon_min) / 10., 2.])
  lat = gen(lat_min, lat_max, lat_step)
  lon = gen(lon_min, lon_max, lon_step)
  n_lat = n_elements(lat)
  n_lon = n_elements(lon)
endelse

x = fltarr(n_lon, n_lat)			; Calculate field components at
y = fltarr(n_lon, n_lat)			; grid points.
z = fltarr(n_lon, n_lat)
h = fltarr(n_lon, n_lat)
dip = fltarr(n_lon, n_lat)
dec = fltarr(n_lon, n_lat)
print, 'Calculating...'
for i = 0, n_lat - 1 do begin
  for j = 0, n_lon - 1 do begin
    print, 'Done ', 100. * (i * n_lat + j) / (n_lat * n_lon), '%...'
    igrf, lat(i), lon(j), year, aux_x, aux_y, aux_z, aux_hor, aux_dip, aux_dec
    x(j,i) = aux_x
    y(j,i) = aux_y
    z(j,i) = aux_z
    h(j,i) = aux_hor
    dip(j,i) = aux_dip
    dec(j,i) = aux_dec
  endfor
endfor

title = 'IGRF ' + strcompress(1900 + year,/remove_all)
again:
print, 'Plot options:'			; Inquire plot option, or exit.
print
print, '1. X component'
print, '2. Y component'
print, '3. Z component'
print, '4. H component'
print, '5. Dip angle I'
print, '6. Declination D'
print, '7. Exit'
print
read, 'Enter option: ',index
if keyword_set(limit) then begin

case index of
1: map_plot, x, lon, lat, genn(min(x), max(x), 20), title + ' - X ', limit = limit
2: map_plot, y, lon, lat, genn(min(y), max(y), 20), title + ' - Y ', limit = limit
3: map_plot, z, lon, lat, genn(min(z), max(z), 20), title + ' - Z', limit = limit, /equ 
4: map_plot, h, lon, lat, genn(min(h), max(h), 20), title + ' - H', limit = limit
5: map_plot, dip, lon, lat, genn(min(dip), max(dip), 20), title + ' - I', limit = limit, /equ
6: map_plot, dec, lon, lat, genn(min(dec), max(dec), 20), title + ' - D', limit = limit
7: goto, exit
endcase
endif else begin
case index of
1: map_plot, x, lon, lat, gen(-.1,.4,.02), title + ' - X '
2: map_plot, y, lon, lat, gen(-.16,.4, .02), title + ' - Y '
3: map_plot, z, lon, lat, gen(-.56, .56, .04), title + ' - Z'
4: map_plot, h, lon, lat, gen(0., .4, .02), title + ' - H'
5: map_plot, dip, lon, lat, findgen(19) * 10. - 90., title + ' - I'
6: map_plot, dec, lon, lat, findgen(19) * 10. - 90., title + ' - D'
7: goto, exit
endcase
endelse
if index ne 7 then goto, again
exit:
return
end