function sunrot, image, from_time, to_time, radius, center, $
	parameters = parameters, solid = solid, B0 = B0, _extra = _extra


;+
;  NAME:
; 	SUNROT
;
;  PURPOSE:
;  Differential or solid rotation of the Sun.
;
;  CATEGORY:
;  Mapping, Image processing.
;
;  CALLING SEQUENCE:
;  RESULT = SUNROT(Image, From_time, To_time, Radius [, Center, $
;		PARAMETERS = parameters, MISSING = missing, CUBIC = cubic, SOLID = solid, B0 = b0])
;
;  INPUTS:
;  Image: Array to be wrapped. This array may be of any numerical type.
;
;  From_time: Time of the observation according to SolarSoft conventions.
;
;  To_time: Time to which rotation is desired (also according to SolarSoft conventions).
;
;  Radius: Solar radius in pixels.
;
;
;  OPTIONAL INPUT PARAMETERS:
;
;  Center:  the solar center relative to the lower left corner. By default the center
; 	is the center of the array.
;
;  KEYWORD PARAMETERS:
;
;  PARAMETERS: rotational paramaters.
; 	By default, the following expressions are used:
;
; 	delta_longitude = (13.39 - 2.7*(sin(phi))^2)*time_interval - for the differential rotation
; 	delta_longitude = 360/27.2753*time_interval - for the solid rotation
;
; 	Using the keyword PARAMETERS, these values can be modified.
;
;  SOLID:   Sets solid rotation with a period of 27.2753 days by default or uses the value
; 	specified by the keyword PARAMETERS, i.e., PARAMETERS = 360/period.
;
;  B0: 	Two-element array to explicitly specify the latitude of the solar disk center for
; 	the observation time (B0[0]) and for the desired time to which the image should be
; 	rotated (B0[1]). By default, the B0 is calculated using the GET_RB0P SolarSoft
;	function.
;
;  Keyword parameters of the INTERPOLATE function (see IDL documentation):
;
;  CUBIC: Cubic convolution parameter for the cubic convolution interpolation method.
;         Set this keyword to a value between -1 and 0 to use the cubic convolution interpolation method.
;
;  MISSING: The value to return for elements outside the bounds of the array.
;
;  GRID: Don't use it! Setting GRID = 1 causes an attempt to allocate a huge array in memory.
;
;
;  OUTPUTS:
; 	Wrapped image.
;
;  COMMON BLOCKS:
; 	None.
;
;  SIDE EFFECTS:
; 	The current graphics device is switched to Z buffer, where the MAP_SET procedure is called.
;	Then the initial graphics device is restored.
;
;  RESTRICTIONS:
; 	The SUNROT function uses ANYTIM and GET_RB0P functions; thus, basic SolarSoft routines
;	should be installed.
;
;	The keyword parameter GRID = 1 (/GRID) should not be set.
;
;
;  PROCEDURE:
;	The current graphics device is switched to Z buffer to enable operability even if
;	no graphics connection with the computer is available. The MAP_SET function
;	establishes the orthographic projection corresponing to the observation time.
;	Coordinates of all pixels inside the disk area are converted to longitudes and
;	latitudes. The longitudes are shifted according to the differential or solid
;	rotation. Then, the MAP_SET function re-establishes the orthographic projection
;	corresponing to the time, to which the rotation is performed, and longitudes
;	and latitudes are converted back to the rectangular coordinates of the pixels.
;	The interior of the disk is replaced with the initial image interpolated
;	to the new coordinates, while the outer part of the image is left unchanged.
;
;  EXAMPLE:
;
; Specify the center and radius of the sun in pixels
;
; 	Center = [200, 300.]
; 	Radius = 180.
;
; Specify times of observation and rotation
;
;	From_time = '12-sep-95 03:00:00'
;	To_time = '14-sep-95 17:00:00'
;
; Calculate B0 for the observation time
;
;	B0 = (get_rb0p(anytim(From_time, /yy)))[1]*!radeg
;
; Simulate a solar image in a graphics window:
;
;	window, 0, xsize = 512, ysize = 512
;
; Set the orthographic projection
;	map_set, B0, 0, 0, /ortho, /nobor
;
; Set proper conversion between DEVICE and NORMAL coordinates:
;
;	!x.s=[Center[0], Radius] / float(!d.x_size)
;	!y.s=[Center[1], Radius] / float(!d.y_size)
;
; Plot heliographic grid
;
;	map_grid, latdel = 15, londel = 15, glinestyle = 0, /horizon, /label
;
; Read the image from the graphics window
;
;	image = tvrd()
;
; Wrap image
;
;	rotated_image = sunrot(image, From_time, To_time, Radius, Center, cubic = -0.5)
;
; Display wrapped image in another window
;
;	window, 1, xsize = 512, ysize = 512
;	tvscl, rotated_image
;
;  ACKNOWLEDGMENTS: VG is grateful to D.Stern, V.Garaimov, M.Aschwanden, V.Slemzin,
;	E.Golubeva, and S.White for fruitful ideas and suggestions, which are used in
;	this code.
;
;  MODIFICATION HISTORY:
;
;	ISTP SD RAS, 1996.
;	Victor Grechnev (grechnev@iszf.irk.ru): Initially written.
;
;	ISTP SD RAS, Feb. 2006.
;	VG: Essentially rewritten to introduce interpolation and simplify the routine.
;-


; Save coordinate system:

X_save = !X
Y_save = !Y
Z_save = !Z
Map_save = !map

Sz = size(image)

if n_elements(center) lt 2 then center = [(Sz[1]-1)*0.5, (Sz[2]-1)*0.5]

	if n_elements(B0) ne 2 then begin

B00 = (get_rb0p(anytim(from_time, /yy)))[1]*!radeg
B01 = (get_rb0p(anytim(to_time, /yy)))[1]*!radeg

	endif else begin
	B00 = B0[0]
	B01 = B0[1]
	endelse

solid_rotation = keyword_set(solid)

	CASE solid_rotation OF

0: if n_elements(parameters) lt 2 then parameters = [13.39, -2.7]
1: if n_elements(parameters) lt 1 then parameters = 360/27.2753

	ENDCASE

; Save the name of the initial graphics device:
init_dev = !d.name

; Set Z buffer graphics device:
set_plot, 'Z'

device, SET_RESOLUTION = Sz[1:2]

; Set the orthographic projection for the observation time:
map_set, B01, 0, 0, /ortho, /noerase, pos=[0,0,1,1], /nobor

; Set the proper conversion between DEVICE and NORMAL coordinates:
!x.s=[Center[0], Radius] / float(!d.x_size)
!y.s=[Center[1], Radius] / float(!d.y_size)

; Prepare X and Y to specify a circle parametrically:

N_arg = mean(Sz[1:2])*4

arg = findgen(N_arg)/(N_arg-1)*2*!pi

x_arg = Center[0] + Radius*cos(arg)
y_arg = Center[1] + Radius*sin(arg)

; Interior of the solar disk, 1d index
disk_index = polyfillv(x_arg, y_arg, Sz[1], Sz[2])

; Convert 1d index to 2d:
disk_x = fix(disk_index mod Sz[1])
disk_y = fix(disk_index / Sz[1])

; Convert coordinates from pixels to longitudes & latitudes:
lonlat = convert_coord(disk_x, disk_y, /dev, /to_data)

; Time difference in seconds:
delta_t = float(anytim(from_time)-anytim(to_time))/86400.

; Compute the longitudinal shift: scalar if solid, array if differential:
if solid_rotation then diffrot = delta_t*parameters[0] else $
	diffrot = (parameters[0]+parameters[1]*(sin(lonlat[1,*]*!dtor))^2)*delta_t

; shift longitudes and reset multiples of 360 degree:
lonlat[0,*] = (lonlat[0,*] + diffrot) mod 360

; Set the orthographic projection for the desired time to account for the B0 difference:
map_set, B00, 0, 0, /ortho, /noerase, pos=[0,0,1,1], /nobor

; Set proper conversion between DEVICE and NORMAL coordinates:
!x.s=[Center[0], Radius] / float(!d.x_size)
!y.s=[Center[1], Radius] / float(!d.y_size)

; Convert NEW coordinates from longitudes & latitudes to pixels:
New_plane_coord = convert_coord(lonlat, /data, /to_dev)

; Fix possible hits beyond the spherical surface:
bad_ind = where(finite(New_plane_coord) ne 1, count)

if count ne 0 then New_plane_coord[bad_ind] = -1

; New image is the same as initial (outside of the disk):
new_image = image

; Replace  the interior of the disk with the old image interpolated to the new coordinates:
new_image[disk_index] = reform(interpolate(image, New_plane_coord[0,*], New_plane_coord[1,*], _extra = _extra))

; Return back to the initial graphics device:
set_plot, init_dev

; Restore initial coordinate system:

!X = X_save
!Y = Y_save
!Z = Z_save
!map = Map_save


return, new_image

end