function turn, dgram,x0,y0,r,tim, asec_per_pix = asec_per_pix,  $
               b = b, xmn = xmn, ymn = ymn
;+
; NAME:
;      turn
; PURPOSE:
;      This function makes the appropriate transformation for a solar
;      rotation of rot radians with the solar rotation axis inclined b
;      radians out of the plane-of-the-sky (at the north pole).
; CALLING SEQUENCE:
;      rotdg = turn(dgram, x0, y0, r, tim [,asec_per_pix=asec_per_pix,
;                   b = b, xmn = xmn, ymn = ymn])
; INPUTS:
;      dgram        = the original(unrotated) dopplergram array
;      x0,y0        = the e/w and n/s solar coordinates(in pixels) of
;                     the lower-left (south-ease) corner of the 
;                     original dopplergram
;      r            = the solar radius (in arcsec)
;      tim          = the projection time, pos. or neg. (in hours)
; OPTIONAL INPUTS:
;      asec_per_pix = arcsec per pixel, default value is 1.
;      b            = the solar b-angle (in degrees), default value is 0
;      xmn, ymn     = the x/y - coordinate of the south-east corner of
;                     the projected image (in pixels)
;                     if not present, the field-of-vies is not to be
;                     prescribed
; OUTPUTS:
;      rotdg        = the rotated dopplergram array
; COMMON BLOCKS:
;      none
; RESTRICTIONS:
;
; PROCEDURE:
;
; PROCEDURE CALLED:
;      none
; MODIFICATION HISTORY:
;      Convert Jim Klimchuk's FORTRAN subroutine TURN into IDL function
;      by Fei-Mei Lee Chou on NOV. 5, 1991
;-
;   
;   
;
if (n_params() eq 0) then return, $
 'turn( dgram,x0,y0,r,tim [,asec_per_pix=asec_per_pix,b=b,xmn=xmn,ymn=ymn])'
pi = 3.14159265
pihalf = 1.57079633
ttim = tim * .261799388    ;2*pi/24. = 0.261799388
if KEYWORD_SET(asec_per_pix) then rsun = r / asec_per_pix else rsun=r
;rsun=float(fix(rsun))
rsunp = rsun + 1
rsun2 = rsun * rsun
if KEYWORD_SET(b) then bb = b / 180. * pi else bb = 0.
;
sz = size(dgram)
if (sz(0) eq 2) then begin
   ix = sz(1)
   iy = sz(2)
endif else begin
   message,' dgram dimension error'
endelse
if(sz(sz(0)+1) ne 4) then dgram=float(dgram)
;
temp = fltarr(ix, iy)
if (tim eq 0.) then begin
   rotdg = dgram
   imin = 0
   jmin = 0
endif else begin
;
;  this section determines the transformed coordinates for each point
;  in the original image
;
;  rotates cartesian axes about x-axis
;
   rotdg = replicate(1.e+36, ix, iy)
   x_indx = indgen(ix)
   y_indx = transpose(indgen(iy))
   x_indx = rebin(x_indx, ix, iy, /sample)    
   y_indx = rebin(y_indx, ix, iy, /sample)    
   x  = x0 + x_indx
   y  = y0 + y_indx
   xy2 = x * x + y * y                     ;x**2 + y**2
   ss = where(xy2 le rsun2)
   ssout = where(xy2 gt rsun2)
   if (ss(0) ne -1) then begin
      z = fltarr(ix, iy)
      z(ss) = sqrt(rsun2 - xy2(ss))
      p = fltarr(ix, iy)
      p(ss) = sqrt(y(ss)*y(ss) + z(ss)*z(ss))
;     compute alpha
      ss0 = where((y eq 0.) and (xy2 le rsun2))
      ss1 = where((y ne 0.) and (xy2 le rsun2))
      alpha = fltarr(ix, iy)
      if(ss1(0) ne -1) then begin
         alpha(ss1) = atan(z(ss1)/y(ss1))
         ssm = where(alpha lt 0.)
         if(ssm(0) ne -1) then alpha(ssm) = alpha(ssm) + pi
      endif 
      if(ss0(0) ne -1) then alpha(ss0) = pihalf
      sxx = where(alpha gt pi)
;
      x1 = fltarr(ix, iy)
      y1 = x1
      z1 = x1
;
      x1(ss) = x(ss)
      y1(ss) = p(ss) * cos(alpha(ss) - bb)
      z1(ss) = p(ss) * sin(alpha(ss) - bb)
; transforms to spherical coordinates
      theta = fltarr(ix, iy)
      phi   = fltarr(ix, iy)
      alat  = fltarr(ix, iy)
      w     = fltarr(ix, iy)
      rot   = fltarr(ix, iy)
      theta1= fltarr(ix, iy)
      phi1  = fltarr(ix, iy)
;
      theta(ss) = acos(y1(ss) / rsun)
      phi(ss) = acos(z1(ss) / (rsun * sin(theta(ss))))
      ssm = where(x1 lt 0.)
      if(ssm(0) ne -1) then phi(ssm) = -phi(ssm)
; rotates about z-axis (solar rotation)
; "w" is the synodic rotation rate (rotations/day)
      alat(ss) = pihalf - theta(ss)
      w(ss) = (1.-125.*(sin(alat(ss)))^2 -.16667*(sin(alat(ss)))^4) / 27.
      rot(ss) = w(ss) * ttim
      theta1(ss) = theta(ss)
      phi1(ss) = phi(ss) + rot(ss)
; transforms back to 'tilted' cartesian coordinates
      x1(ss) = rsun * sin(theta1(ss)) * sin(phi1(ss))
      z1(ss) = rsun * sin(theta1(ss)) * cos(phi1(ss))
; rotates back to original (plane-of-sky) cartesian coordinates
      p(ss)  = sqrt(y1(ss) * y1(ss) + z1(ss) * z1(ss))
      ss0 = where(( y1 eq 0.) and ( z1 ge 0.) and (xy2 le rsun2))
      ss1 = where(( y1 eq 0.) and ( z1 lt 0.) and (xy2 le rsun2))
      ss2 = where(( y1 ne 0.) and (xy2 le rsun2))
      gama = fltarr(ix, iy)
      if(ss2(0) ne -1) then begin
         gama(ss2) = atan(z1(ss2) / y1(ss2))
         ss3 = where(gama lt 0.)
         if(ss3(0) ne -1) then gama(ss3) = gama(ss3) + pi
         ss3 = where(z1 lt 0.)
         if(ss3(0) ne -1) then gama(ss3) = gama(ss3) + pi
      endif 
      if(ss0(0) ne -1) then gama(ss0) = pihalf
      if(ss1(0) ne -1) then gama(ss1) = -pihalf
;
      alpha(ss) = gama(ss) + bb
      sxx = where(alpha gt pi)
      y(ss) = p(ss) * cos(alpha(ss))
      z(ss) = p(ss) * sin(alpha(ss))
      x(ss) = x1(ss)
      ss0 = where((z lt 0.) and (xy2 le rsun2))
      if(ss0(0) ne -1) then begin
         print,'size of ss0=', size(ss0)
;         x00=x(0,0)
         x(ss0) = 9999.
;         x(0,0) = x00
      endif
   endif    ;end of if(ss0(0) ne -1), (xy2 le rsun2)
   if(ssout(0) ne -1) then rotdg(ssout) = 1.e+36
;
   inew = intarr(ix, iy)
   jnew = inew
   imin = fix(x(0,0) - x0 + 0.5)
   jmin = fix(y(0,0) - y0 + 0.5)
;   if((KEYWORD_SET(xmn)) and (KEYWORD_SET(ymn))) then begin
   if((n_elements(xmn) ne 0) and (n_elements(ymn) ne 0)) then begin
      inew_o = fix(x(0,0) - xmn + 0.5)
      jnew_o = fix(y(0,0) - ymn + 0.5)
      if((inew_o lt -(ix-1)) or (inew_o gt (ix-1))) then $
         message,'rotation entirely out of field of view'   
      if((jnew_o lt -(iy-1)) or (jnew_o gt (iy-1))) then $
         message,'rotation entirely out of field of view'   
      inew = fix(x - xmn + 0.5)
      jnew = fix(y - ymn + 0.5)
   endif else begin
      inew = fix((x - x0 + 0.5)) - imin 
      jnew = fix((y - y0 + 0.5)) - jmin 
   endelse
;
   ss0 = where((inew ge 0) and (inew lt ix) and (jnew ge 0) and   $
               (jnew lt iy) and (dgram ge 1.e+35))
   ss1 = where((inew ge 0) and (inew lt ix) and (jnew ge 0) and   $
               (jnew lt iy) and (dgram lt 1.e+35))
   if(ss0(0) ne -1) then begin
      rotdg(inew(ss0), jnew(ss0))  = 1.e+36
      temp(inew(ss0), jnew(ss0)) = 1.
   endif
   if(ss1(0) ne -1) then begin
      knew=intarr(ix,iy)
      knew(ss1) = 1
      numss = n_elements(ss1)
      for s = 0, numss-1 do begin
        if(rotdg(inew(ss1(s)),jnew(ss1(s))) eq 1.e+36) then begin
          rotdg(inew(ss1(s)), jnew(ss1(s))) = dgram(ss1(s))
          temp(inew(ss1(s)), jnew(ss1(s))) = 1.
        endif else begin
          rotdg(inew(ss1(s)),jnew(ss1(s))) =       $
             rotdg(inew(ss1(s)),jnew(ss1(s))) + dgram(ss1(s))
          temp(inew(ss1(s)),jnew(ss1(s)))=temp(inew(ss1(s)),jnew(ss1(s)))+1.
        endelse
      endfor

   endif
;
   iminn = imin + x0
   jminn = jmin + y0
   print,' new x0, y0 = ', iminn, jminn
   ss0 = where(temp gt 1.)
   if(ss0(0) ne -1) then rotdg(ss0) = rotdg(ss0) / temp(ss0)
   x = x0 + imin + x_indx
   y = y0 + jmin + y_indx
   rho = fltarr(ix, iy)
   rho = sqrt(x * x + y * y)
   ss0 = where((rho lt rsunp) and (rho ge rsun))
   if(ss0(0) ne -1) then begin
      rotdg(ss0) = 1.e+36
      temp(ss0) = 1.
   endif
;
;  this section determines average values for 'blank' pixels in the
;  new image (rotdg)
;
   temp2 = temp
   temp2(0,*) = 2.
   temp2(ix-1,*) = 2.
   nbl = intarr(ix, iy)
   nbu = intarr(ix, iy)
   ss0 = where(temp2 eq 0.)
   if (ss0(0) ne -1) then begin
      numss0 = n_elements(ss0)
      for k = 0, numss0 - 1 do begin
         ll = (ss0(k)/ix) * ix
         uu = ll + ix - 1
         ssl = where(temp(ll:ss0(k)-1) ge 1.0)
         ssu = where(temp(ss0(k)+1:uu) ge 1.0)
         if(ssl(0) ne -1) then begin
            numl = n_elements(ssl)
            nbl(ss0(k)) = ssl(numl-1)+ll
         endif
         if(ssu(0) ne -1) then nbu(ss0(k))=ssu(0)+ss0(k)+1
      endfor
   endif
;
   temp2 = temp
   ss1 = where((nbl ne 0) and (nbu ne 0))
   if(ss1(0) ne -1) then begin
      rotdg(ss1) = (rotdg(nbl(ss1)) + rotdg(nbu(ss1))) / 2.
      temp2(ss1) = 2.
   endif
   ss2 = where(rotdg gt 1.e+30)
   if(ss2(0) ne -1) then rotdg(ss2) = 1.e+36
   ss3 = where(temp2 eq 0.)
   if(ss3(0) ne -1) then rotdg(ss3) = 1.e+36
endelse
return, rotdg
end