subr trace_rot_shifts, xc, yc, tai, dx, dy, xref, yref, tairef
 ;does a vector of shifts given a single ref
 ;1/4/2001  allow arbitrary reference in (x,y,t) by adding tairef
 ;compute the shifts to correct for solar rotation and re-pointings
 ;often assume that the reference lat/long is the center of the first
 ;image but these are not entered as separate arguments
 ;The pixel scale doesn't come into play here since xc, yc are assumed to
 ;be in arc sec already. The computed (dx,dy) will be in arc sec also.
 ;the first dimensions of xc and yc must agree with tai
 ;need B angles and radii
 n = dimen(tai, 0)
 if n ne dimen(xc, 0) then {
 	ty,'dimension error in trace_rot_shifts'
 	return }
 b = sun_b_tai(tai)*#d.r
 r = sun_r_tai(tai)
 bref = b(0)	rref = r(0)
 ;10/19/2000 - redim b and r to allow for 2-D xc, yc, tai
 redim, b, dimen(tai)
 redim, r, dimen(tai)
 ;check if on the sun
 if sqrt(xref*xref + yref*yref) ge rref then {
   ;oops, we have a ref point off the sun, just apply pointing changes
   x = xref	y = yref
 } else {
   latref = (cos(bref)*yref/rref + sin(bref)*sqrt(1.-(xref*xref+yref*yref)/(rref*rref)))
   ;rate includes correction for differential rotation
   ;;rate = 13.27 - latref*latref*1.68 - (latref^4)*2.40
   ;added 8/16/2009 to be consistent with rot_shift in sunlib except that
   ;type is 1 by default
   if defined($solar_rotation_type) eq 0 then $solar_rotation_type = 1
   if notdefined($solar_rigid_rate) then $solar_rigid_rate = 13.39
   ;$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 - latref*latref*!solar_rotate_b - (latref^4)*!solar_rotate_c
     rate = 0.0
   endcase
   ;;rate = !solar_rotate_a - latref*latref*!solar_rotate_b - (latref^4)*!solar_rotate_c
   ty, 'trace_rot_shifts solar rotation rate =', rate
   latref = asin(latref)
   longref = asin(xref/cos(latref)/rref)
   ;use rigid rotation to update the long for each time
   ;;dlon = 13.39 * (tai - tai(0))/86400.
   ;now upgraded to include differential rotation
   dlon = rate * (tai - tairef)/86400.
   longt = longref + dlon*#d.r
   ;1/4/2001 - these are presently NOT clipped to a +- 90 range so a ref that
   ;would roll over the limb during the time span will have interesting effects
   ;now convert the lat/long back to (x,y) for all times
   x = r*cos(latref)*sin(longt)
   y = 0.5*r*(sin(latref-b)*(1.+cos(longt))+sin(latref+b)*(1.-cos(longt)))
 }
 ;xc and yc could be 2-D arrays, if so we have to adjust x and y
 if (num_dim(xc) eq 2) and (num_dim(x) eq 1) then { redim, x, n, 1   redim, y, n, 1 }
 dx = x - xc	dy = y - yc
 endsubr
 ;=====================================================================