
function sunrota, sun, date0, time0, date1, time1, Radius

Sz = size(sun)

x0_ = float(Sz(1)/2)
y0_ = float(Sz(2)/2)
r0  = float(Radius)

if strlen(date0) eq 8 then date6_0 = date0 else $
	date6_0 = strmid(date0, 8, 2) + '/' + strmid(date0, 5, 2) + '/' + strmid(date0, 2,2)

	day0 = strmid(date6_0, 0,2)
	month0 = monthnames(fix(strmid(date6_0, 3,2)))
	year0 = strmid(date6_0, 6,2)

un_time0 = day0 + '-' + month0 + '-' + year0 + ' ' + time0

if strlen(date1) eq 8 then date6_1 = date1 else $
	date6_1 = strmid(date1, 8, 2) + '/' + strmid(date1, 5, 2) + '/' + strmid(date1, 2,2)

	day1 = strmid(date6_1, 0,2)
	month1 = monthnames(fix(strmid(date6_1, 3,2)))
	year1 = strmid(date6_1, 6,2)

un_time1 = day1 + '-' + month1 + '-' + year1 + ' ' + time1

rb0p = get_rb0p(un_time0)

b0 = rb0p(1)*!radeg

p0 = 0.

rdeg = float(time_difference(date6_0, time0, date6_1, time1)/3600.)

dim=size(sun)
rsun=float(sun)
image=rsun
zmin=min(rsun)
nx=dim(1)
ny=dim(2)

l0=rdeg*360.0/27.2753/24.0
x0	=x0_-1.   ;IDL convention
y0	=y0_-1.	;IDL convention
ir	=r0 ;solar radius in units of EW pixels
ix	=findgen(nx)

for iy=0,ny-1 do begin
;   if (long(iy/50) eq float(iy)/50.) then PRINT,'processing line =',iy
   yy	=iy
   xx	=ix
   ind	=where((xx-x0)^2+(yy-y0)^2 lt (ir^2))

   if (ind(0) ne -1)   then begin
      x =xx(ind)
      heliotrans, x0, y0, 0., p0, b0, 0., x, yy, ir, hlong, hlat

;	  diffrot=(2.7*(sin(hlat*!pi/180.))^2)*rdeg/24.
	  diffrot=-(2.7*abs(sin(hlat*!pi/180.)))*rdeg/24.*0.


      hlong  = hlong + diffrot     ;differential rotation + rotation
      heliotrans2, x0, y0, 0., 0., b0, l0, hlong, hlat, ir, ix2, iy2

	i1=long(ix2-0.5) > 0
	i2=i1+1 < (nx-1)
	j1=long(iy2-0.5) > 0
	j2=j1+1 < (ny-1)

	z1=image(i1,j1)
	z2=image(i2,j1)
	z3=image(i2,j2)
	z4=image(i1,j2)

	t=ix2-0.5-float(i1)
	u=iy2-0.5-float(j1)

      zz=(1-t)*(1-u)*z1+t*(1-u)*z2+t*u*z3+(1-t)*u*z4  ;bilinear interpol.

      rsun(ind,iy)=float(zz)

   endif
endfor

return,rsun

end
