FUNCTION destretch_spot, p, ang, scl, NOCLIP=noclip
;+
; NAME:
;       Destretch_Spot
; PURPOSE:
;       Remove geometrical distortions due to position on sun
; CALLING SEQUENCE:
;       Res = Destretch_Spot( Img, Angle, Scale )
; INPUTS:
;       Img:   Image to correct
;       Angle: Angle (degree) between Y-axis and direction towards disc center
;       Scale: Scale factor.  Equal to mu=cosine(theta)
; OUTPUTS:
;       Res:   resampled image.  The pixel size stays unchanged,
;              i.e. is still the same as before, but now in both
;              directions, so we get 'true kilometers' on the surface.
; PROCEDURE:
;       
; MODIFICATION HISTORY:
;       25-Apr-2001  P.Suetterlin, SIU
;       03-May-2001  xp, yp have to be float arrays
;       09-May-2001  New (much faster) way to compute xp,yp
;-

  ;;; convert angle from degree to radian
a = ang*!dtor
  ;;; rotation matrices for left/right rotation
r1 = [[cos(a), sin(a)], [-sin(a), cos(a)]]
r2 = [[cos(-a), sin(-a)], [-sin(-a), cos(-a)]]
  ;;; matrices for scaling
t = [[1., 0], [0, 1./scl]]
t1 = [[1., 0], [0, scl]]

  ;;; compute the destretched positions of the 3 image corners (0,0 stays)
s0 = (size(p))(1:2)
pp1 = [0, 0]
pp2 = fix(r2#t#r1#[s0(0), 0])
pp3 = fix(r2#t#r1#[0, s0(1)])
pp4 = fix(r2#t#r1#[s0(0), s0(1)])

  ;;; find extrema of size
xmi = min([pp1(0), pp2(0), pp3(0), pp4(0)], max=xma)
ymi = min([pp1(1), pp2(1), pp3(1), pp4(1)], max=yma)
  ;;; and compute new image dimensions
sx = xma-xmi+1
sy = yma-ymi+1

xp = fltarr(sx, sy)
yp = fltarr(sx, sy)

  ;;; must be easier, but I'm too stupid to see:
  ;;; fill the x,y coordinates of the new image from backtransforming
  ;;; their coordinates to 'observed' ones

;FOR i=xmi, xma DO BEGIN
;    FOR j=ymi, yma DO BEGIN 
;        tmp = r2#(t1#(r1#[i, j]))
;        xp(i-xmi, j-ymi) = tmp(0)
;        yp(i-xmi, j-ymi) = tmp(1)
;    ENDFOR 
;ENDFOR

np1 = r2#(t1#(r1#[xmi, ymi]))
np2 = r2#(t1#(r1#[xma, ymi]))
np3 = r2#(t1#(r1#[xmi, yma]))
np4 = r2#(t1#(r1#[xma, yma]))

xp(*, 0) = np1(0)+dindgen(sx)/(sx-1)*(np2(0)-np1(0))
dy = double(np3(0)-np1(0))/(sy-1)
FOR i=1, sy-1 DO xp(*, i) = xp(*, 0)+dy*i

yp(0, *) = np1(1)+dindgen(sy)/(sy-1)*(np3(1)-np1(1))
dx = double(np2(1)-np1(1))/(sx-1)
FOR i=1, sx-1 DO yp(i, *) = yp(0, *)+dx*i

  ;;; interpolate the image
res = interpolate(p, xp, yp, miss=0, cub=-.5)

  ;;; find the rectangular area that contains no missing data
  ;;; 1 pixel safety border to prevent zero values at the corners
x0 = max([pp1(0)+1, pp3(0)+1])-xmi
y0 = max([pp1(1)+1, pp2(1)+1])-ymi
x1 = min([pp2(0)-1, pp4(0)-1])-xmi
y1 = min([pp3(1)-1, pp4(1)-1])-ymi

IF keyword_set(noclip) THEN $
  return, res $
ELSE $
  return, res(x0:x1, y0:y1)

END