PRO rotflat,img,c0,r0,az,za,xang,zang,yang

; procedure to find the center azimuth and zenith angle of an LBS image 
;through an iterative process
; program finds a centre of the LBS, rotates the coordinate system so the z-axis
; points through the centre; finds a new centre, and iterates until the 
; difference between iterations is negligible.
A=FLTARR(30) & B=A
A(0)=0.*!DTOR & B(0)=0.*!DTOR
n=0

IF min(xang) LT 0 THEN xang(where(xang lt 0))=xang(where(xang lt 0))+360.

; start iterations
REPEAT BEGIN
   ;change to polar coords and rotate over angles A & B
   polrec3d,1.,za,-az,x0,y0,z0,/DEGREES
   rot_3d,3,x0,y0,z0,A(n),x1,y1,z1
   rot_3d,2,x1,y1,z1,B(n),x2,y2,z2

   ;display image and azimuth/zenith angle grid 
   ;image_cont,img,/asp
   ;contour,x2,/follow,level=findgen(5)*.25-0.5,c_label=[1,1,1,1,1],/overplot,$
	;color=!d.n_colors/2.
   ;contour,y2,/follow,level=findgen(5)*.25-0.5,c_label=[1,1,1,1,1],/overplot,$
	;color=!d.n_colors/2.

   n=n+1

   lbs=WHERE(img GT MAX(img)/2.)

   ;centre of image located at mean of x, y and z coords of LBS part of image
   i=FLTARR(256,256) & i(lbs)=img(lbs)
   mx=mean(x2(lbs))
   my=mean(y2(lbs))
   mz=SQRT(1-mx^2.-my^2.)

   ;PRINT,mx,my,mz,$
   ;      FORM='("Approximate centre of image located at: ",2(f7.4,","),f7.4)'

   ;recalculation of rotation angles - z axis should be orthogonal to a plane
   ; tangent to the centre of the LBS image
   delta=ACOS(mz)
   zz=(ATAN(my,mx))
   ;PRINT,STRING(delta*!RADEG,zz*!RADEG $
   ;     ,FORMAT='("delta: ",F7.2,"  zz: ",F7.2)')
   B(n)=ACOS(COS(B(n-1))*(COS(delta))-SIN(B(n-1))*SIN(delta)*(COS(zz)))
   d=ASIN(sin(delta)*(SIN(zz)/sin(B(n))))
   A(n)=A(n-1)+d
   ;PRINT,A(n)*!RADEG,B(n)*!RADEG,FORM='("A=",f7.2," B=",f7.2)'
   ;stop iterations when difference between last two iterations is small
ENDREP UNTIL((ABS(A(n)-A(n-1)) LT .001) AND (ABS(B(n)-B(n-1)) LT .001))
zang=A(n)*!RADEG & yang=B(n)*!RADEG
PRINT,zang,yang $
     ,FORMAT='("Angle about z-axis =",F7.2,"; angle about y-axis = ",f7.2)'
RETURN
END