;Upsample an image. This is intended to be the inverse operation of downsample
; (1) To make a digital audio signal higher quality by increasing the sample rate,
;  for example, from 44.1 kHz to 88.2 kHz, and interjecting new samples in
;  between existing samples. The sample size (word size) is also increased for
;  finer granularity. The objective is to have a smoother digital wave going into
;  the digital-to-analog converter. (2) To increase the color depth of a
;  digital image, for example, from 16 bits per pixel to 24 bits per pixel.
;   
; Inputs
;   image_downsample -  two-dimensional array of data to be upsampled. This
;                       is usually the output of downsample. It will have a size [xd+1,yd+1]
;   amount - integer factor by which to upsample. This should be the same number
;            that was passed to downsample when it was called. Special case: if amount is 1,
;            the input image is immediately returned unchanged.
;   /longitude -   set this flag for longitude-type variables only!
; Returns
;   An upsampled version of the image. This will have a size [xd*amount,yd*amount]
function upsample,image_downsample,amount,longitude=longitude
  ;special case: nothing to do if amount is 1
  if amount eq 1 then return,image_downsample
  if keyword_set(longitude) then begin
    ;compute the sines and cosines of the input, then interpolate each of them
    ;and then take the atan of the results. This avoids all singlarities and branch
    ;cuts related to longitude.
    slon_downsample=sin(image_downsample*!const.dtor)
    clon_downsample=cos(image_downsample*!const.dtor)
    slon_upsample=upsample(slon_downsample,amount) ;We specifically do NOT set /longitude here
    clon_upsample=upsample(clon_downsample,amount) ;otherwise this is an infinite loop.
    image_upsample=atan(slon_upsample,clon_upsample)*!radeg
    return,image_upsample
  end else begin
    xd=(size(image_downsample,/dimensions))[0]-1
    yd=(size(image_downsample,/dimensions))[1]-1
    x=xd*amount
    y=yd*amount

    image_upsample=interpolate(image_downsample,findgen(x+1)/amount,findgen(y+1)/amount,/grid,cubic=-0.5)
    image=interpolate(image_upsample,findgen(x)*x/(x-1),findgen(y)*y/(y-1),/grid)
    return,image
  end
end