subr rot_shift, x, y, tai1, tai2, dx, dy
 ;returns expected rotation (rigid) on solar disk for a point at (x,y)
 ;12/13/2004 - upgrade for differential using global flag, if this
 ;flag is not defined, default to original behavior
 ;from tai1 to tai2 where (x,y) is in arc seconds from center of disk
 b = sun_b_tai(tai1)*#d.r
 r = sun_r_tai(tai1)
 if sqrt(x*x + y*y) ge r then {
  ;oops, we have a point off the sun, just return 0
  dx = 0	dy = 0
 } else {
  lat=(cos(B)*y/r + sin(B)*sqrt(1.-(x*x+y*y)/(r*r)))
  lat=asin(lat)
  longt=asin(x/cos(lat)/r)
  ;now rotate the longitude
  zeroifundefined, $solar_rotation_type
  ;$solar_rotation_type = 0 is rigid rotation (old behavior), 1 is
  ;differential and 2 is none
  if $solar_rotation_type lt 0 or $solar_rotation_type gt 2 then $solar_rotation_type=0
  ncase $solar_rotation_type
    rate = $solar_rigid_rate	;usually 13.39
    rate = !solar_rotate_a - lat*lat*!solar_rotate_b - (lat^4)*!solar_rotate_c
    rate = 0.0
  endcase
  dlon = rate * (tai2 - tai1)/86400.
  ;;ty,'longitude shift =', dlon
  longt = longt + dlon*#d.r
  ;use updated r and b to go back
  b = sun_b_tai(tai2)*#d.r
  r = sun_r_tai(tai2)
  x2 = r*cos(lat)*sin(longt)
  y2 = 0.5*r*(sin(lat-b)*(1.+cos(longt))+sin(lat+b)*(1.-cos(longt)))
  dx = x2 - x	dy = y2 - y
 }
 ty,'rotation shift =',dx, dy
 endsubr
 ;=======================================================================