function sunrotate,sun,x0_,y0_,r0,b0,p0,rdeg,sld,b1=b1,noint=noint
dim=size(sun)
rsun=float(sun)
image=rsun
zmin=min(rsun)
nx=dim(1)
ny=dim(2)
if n_elements(noint) ne 0 then nint=noint(0) else nint=0
if n_elements(b1) ne 0 then b1=b1(0) else b1=b0
l0=double(rdeg)/24.0d*360.0d/27.2753d
x0	=x0_-1.   ;IDL convention
y0	=y0_-1.	;IDL convention 
ir	=r0 ;solar radius in units of EW pixels
pi	=acos(-1.)
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.0,p0,b1,0.0,x,yy,ir,hlong,hlat
      diffrot=sld*l0*(3.0/13.45)*(sin(hlat*pi/180.))^2
      hlong  =hlong+diffrot     ;differential rotation + rotation
      heliotrans2,x0,y0,0.0,0.0,b0,l0,hlong,hlat,ir,ix2,iy2
      if nint eq 1 then begin
        ix2=round(ix2+1) & iy2=round(iy2+1) 
        ix2=ix2>0<nx-1 & iy2=iy2>0<ny-1
        rsun(ind,iy)=image(ix2,iy2)
      endif else begin
       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)
      endelse
   endif 
endfor

return,rsun
end