pro euvi_heliorot,ix,iy,ir,pos,blat,blong,pos_,blat_,blong_,ix_,iy_
;transforms pixel coordinates [IX,IY,IR] with center [IX,IY,IR]=[0,0,0]
;into heliografic coordinates, performs differential rotation from
;BLONG to BLONG_, and transform back into pixel coordinates [IR,IY,IR]
;with POS_,BLAT_,BLONG_ (angles given in degrees)
;POS is "effective" position angle of evt. rotated image (POS=POS0+CROTA2)
;BLAT,BLONG is heliografic longitude and latitude of disk center.
;Vectorized version for IX and IY

PI	=ACOS(-1.)		&EPS	 =1.E-8		
POSRAD	=-POS*PI/180.		&POSRAD_ =+POS_*PI/180.		
BLATRAD =BLAT*PI/180.		&BLATRAD_=BLAT_*PI/180.
SINPOS	=SIN(POSRAD)		&SINPOS_ =SIN(POSRAD_)		
COSPOS	=COS(POSRAD)		&COSPOS_ =COS(POSRAD_)
SINBLAT =SIN(BLATRAD)		&SINBLAT_=SIN(BLATRAD_)		
COSBLAT =COS(BLATRAD)		&COSBLAT_=COS(BLATRAD_)

R	=SQRT(IX^2+IY^2) < IR
R	=R > EPS
Z	=SQRT(IR^2-R^2)				;z-coordinate 
X	=IX*COSPOS-IY*SINPOS			;rotation position angle
YY	=IX*SINPOS+IY*COSPOS	
Y	=Z*SINBLAT+YY*COSBLAT			;rotation by BLAT
V	=SQRT((IR^2-Y^2) > EPS)
SINPHI  =(X/V)  <   1.				;longitude difference from cent
SINPHI  =SINPHI > (-1.)				;longitude difference from cent
DLON  	=(180./PI)*ASIN(SINPHI) 		;longitude difference in deg
HLONG	=BLONG+DLON		 		;heliographic longitude
SINLAT	=(Y/IR) <   1.
SINLAT 	=SINLAT > (-1.)
HLAT	=(180./PI)*ASIN(SINLAT)  		;heliographic latitude
SIN2	=(SIN(HLAT*PI/180.))^2
DIFFROT	=(BLONG_-BLONG)*(3.0/13.45)*SIN2
HLONG 	=HLONG+DIFFROT				;differential rotation 
DLON	=HLONG-BLONG_
SINPHI	=SIN(DLON*PI/180.)	
SINLAT	=SIN(HLAT*PI/180.)
Y2	=IR*SINLAT				;HLONG-SIN equatorial coord
X1	=SQRT((IR^2-Y2^2) > EPS)*SINPHI		;HLAT-SIN  equatorial coord.
Z1	=SQRT((IR^2-Y2^2-X1^2) > EPS)		;z-coordinate
Y1	=-Z1*SINBLAT_+Y2*COSBLAT_		;disk center at BLAT,BLONG
IX_	=X1*COSPOS_-Y1*SINPOS_			;position angle rotation
IY_	=X1*SINPOS_+Y1*COSPOS_			;position angle rotation
END