PRO HELIOTRANS,X0,Y0,CROTA2,POS,BLAT,BLONG,IX,IY,IR,HLONG,HLAT
;transforms spherical coordinates [IX-X0,IY-Y0,IR] into cartesian coordinates
;[HLONG,HLAT] of heliografic longitude/latitude.
;
;X0, Y0 are pixel coords of disk center
;POS is position angle to rotate, CROTA2 is rotation angle of image: both zero
;BLAT,BLONG is heliografic longitude and latitude of disk center.
;IX, IY are coords to be rotated: one may be an array (pixels)
;IR is solar radius in units of pixels: height of rotating surface
;HLONG, HLAT are helio lat and long - the output

PI	=ACOS(-1.)		&DPOS	=POS+CROTA2
X	=FLOAT(IX-X0)		&Y	=FLOAT(IY-Y0)
POSRAD	=-DPOS*PI/180.		&BLATRAD=BLAT*PI/180.
SINPOS	=SIN(POSRAD)		&COSPOS	=COS(POSRAD)
SINBLAT =SIN(BLATRAD)		&COSBLAT=COS(BLATRAD)
RXY	=SQRT(X^2+Y^2)
RR	=FLOAT(IR) > RXY
Z2	=RR^2-Y^2-X^2		;z-coordinate squared
Z	=SQRT(Z2 > 0)		;z-coordinate
XX	=X*COSPOS-Y*SINPOS	;rotation position angle
Y1	=X*SINPOS+Y*COSPOS
YY	=Z*SINBLAT+Y1*COSBLAT	;rotation by BLAT
V2	=(RR^2-YY^2)
V	=SQRT(V2 > 0)		 ;radius proj in equator-plane
SINPHI	=IX*0.+1.
ind	=where(v gt 0)
SINPHI(ind)=(XX(ind)/V(ind))	;longitude difference from center
SINPHI	=SINPHI > (IX*0.-1.)
SINPHI	=SINPHI < (IX*0.+1.)
DLON	=XX*0.
ind	=where(sinphi ne 0)
DLON(ind)=(180./PI)*ASIN(SINPHI(ind)) ;longitude difference in degree
HLONG	=BLONG+DLON		 ;heliographic longitude
SINLAT	=(YY/RR)
HLAT	=(180./PI)*ASIN(SINLAT)  ;heliographic latitude
END

; *******************************************

PRO HELIOTRANS2,X0,Y0,CROTA2,POS,BLAT,BLONG,HLONG,HLAT,IR,IX,IY
;transforms cartesian coordinates [HLONG,HLAT] of heliografic
;longitude/latitude into spherical coordinates [X,Y,R]=[IX-X0,IY-Y0,IR]
;POS is position angle, CROTA2 = rotation angle of image
;BLAT,BLONG is heliografic longitude and latitude of disk center.

PI	=ACOS(-1.)		&EPS	=1.E-8
DPOS	=POS+CROTA2		&DLON	=HLONG-BLONG
POSRAD	=+DPOS*PI/180.		&BLATRAD=+BLAT*PI/180.
SINPOS	=SIN(POSRAD)		&COSPOS	=COS(POSRAD)
SINBLAT =SIN(BLATRAD)		&COSBLAT=COS(BLATRAD)
SINPHI	=SIN(DLON*PI/180.)	&SINLAT	=SIN(HLAT*PI/180.)
Y1	=IR*SINLAT		;HLONG-SIN equatorial coord
X1	=SQRT(IR^2-Y1^2)*SINPHI	;HLAT-SIN  equatorial coord.
Z1	=SQRT(IR^2-Y1^2-X1^2 > EPS);z-coordinate
X2	=X1			;x-coordinate
Y2	=-Z1*SINBLAT+Y1*COSBLAT	;disk center at BLAT,BLONG
X3	=X2*COSPOS-Y2*SINPOS	;position angle rotation
Y3	=X2*SINPOS+Y2*COSPOS	;position angle rotation
IX	=X3+X0			;RA-SIN with image center at X0
IY	=Y3+Y0			;DEC-SIN with image center at Y0
END

; *******************************************



function sunrotas, data, from_time, to_time, Radius

dim = size(data)

x0_ = (dim[1]-1)*0.5
y0_ = (dim[2]-1)*0.5
r0  = float(Radius)

rb0p = get_rb0p(from_time)

b0 = rb0p[1]*!radeg

p0 = 0.

rdeg = (anytim(to_time) - anytim(from_time))/3600.

rsun = float(data)
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

   yy	=iy
   xx	=ix
   ind	=where((xx-x0)^2+(yy-y0)^2 lt (ir^2), count)

   if count gt 0   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 = 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