PRO Setup_astrom, telescope, sunephem, naxis1, naxis2, astrom, $
                  ocentre=ocentre, pscale=pscale, vangle=vangle

; Set double precision conversion from degrees to radians.

radeg = 180.0/!dpi

; Define angle of solar north to image vertical.

IF NOT keyword_set(vangle) THEN vangle = 0. $
ELSE IF n_elements(vangle) NE 1 THEN BEGIN
   message, 'Orientation angle must be a scalar!', /info
   message, 'Proceeding with default: solar north straight up!', /info
   vangle = 0.
ENDIF

; Set angle of ecliptic north to image vertical in radians.

crota = (vangle+sunephem.orient.ecliptic)/radeg

; Make CD matrix.

cd11 = cos(crota)
cd12 = sin(crota)
cd21 = -sin(crota)
cd22 = cos(crota)

; Insert fix to allow for different definition of CROTA in SOHO
; keywords than in standard FITS (Greisen & Calabretta 1995).

cd12 = -cd12
cd21 = -cd21

; End of fix.

cd = [[cd11, cd21], [cd12, cd22]]

; Define plate scale for the specified telescope in arcsec/pixel.

IF NOT keyword_set(pscale) THEN BEGIN
    CASE strlowcase(telescope) OF
        'c1': cdelt = [-5.6d0, 5.6d0]
        'c2': cdelt = [-11.9d0, 11.9d0]
        'c3': cdelt = [-56.0d0, 56.0d0]
        Else: message, 'Select Telescope: C1, C2 or C3 !!'
    ENDCASE
ENDIF ELSE IF n_elements(pscale) NE 1 THEN BEGIN
    message, 'Pixel Scale Size must be a scalar!', /info
    message, 'Proceeding with default for selected telescope!', /info
    CASE strlowcase(telescope) OF
        'c1': cdelt = [-5.6d0, 5.6d0]
        'c2': cdelt = [-11.9d0, 11.9d0]
        'c3': cdelt = [-56.0d0, 56.0d0]
        Else: message, 'Select Telescope: C1, C2 or C3 !!'
    ENDCASE
ENDIF ELSE cdelt = [-pscale, pscale]

; Convert to degrees/pixel.

cdelt = cdelt/3.6d3

; Define reference pixel location (tangent point for gnomonic
; projection). Assumed to be centre of sun.

IF NOT keyword_set(ocentre) THEN BEGIN
    crpix = [(float(naxis1)/2.)+0.5, (float(naxis2)/2.)+0.5]
ENDIF ELSE IF n_elements(ocentre) NE 2 THEN BEGIN
    message, 'Location of occulter must be a two-element vector!', $
      /info
    message, 'Proceeding with default: occulter at image centre!', $
      /info
    crpix = [(float(naxis1)/2.)+0.5, (float(naxis2)/2.)+0.5]
ENDIF ELSE crpix = [float(ocentre(0)), float(ocentre(1))]
;
; - CRPIX IS FITS COORDINATE, NOT IDL COORDINATE!!!      
;

; Define reference pixel value to be J2000 solar coordinates.

crval = [sunephem.long.J2000, sunephem.lat.J2000]

; Define projection to be tangent (gnomonic) projection.

ctype = ['RA---TAN', 'DEC--TAN']
longpole = 180.0
latpole  = 0.0
projp1 = -1.
projp2 = -2.

; Make astrometry structure.

astrom = {CD: double(cd), CDELT: double(cdelt), CRPIX: float(crpix), $
      CRVAL: double(crval), CTYPE: string(ctype), $
          LONGPOLE: longpole, LATPOLE: latpole, $
          PROJP1: projp1, PROJP2: projp2, PV2:[0.,0.] }

return
end