pro read_princeton, file, $
                    data, $
                    header=header, $
                    x_calibration=x_calibration, $
                    y_calibration=y_calibration
;+
; NAME:
;    READ_PRINCETON
;
; PURPOSE:
;    This procedure reads data files written by Princeton Instruments'
;    WinSPEC and WinVIEW software.
;
; CATEGORY:
;    File input.
;
; CALLING SEQUENCE:
;    READ_PRINCETON, File, Data, Header
;
; INPUTS:
;    File:    The name of the data file to read.
;
; OUTPUTS:
;    Data[nx, ny, nframes]:
;        The output data array. The array will be 1, 2 or 3 dimensions
;        (depending upon whether ny and nframes are 1) and can be integer,
;        long or float data type.
;
; KEYWORD OUTPUTS:
;    Header: The 4100 byte header from the file. This header can be used to
;            extract additional information about the file. See the
;            Princteon Instruments "PC Interface Library Programming Manual"
;            for the description of the header structure, and this procedure
;            for examples of how to extract information from the header.
;
;    X_calibration[nx]:
;            An array of calibrated values for each pixel in the X
;            direction.
;    Y_calibration[nx]:
;            An array of calibrated values for each pixel in the Y
;            direction.
;
; RESTRICTIONS:
;        This procedure currently only extracts the data size from the header.
;        It should be exhanced to extract calibration information, etc.
;
; EXAMPLE:
;    Read a data file:
;
;    IDL> READ_PRINCETON, 'test.spe', data, header=header, x_cal=x_cal
;    IDL> plot, x_cal, data
;    IDL> clock_speed = float(header, 1428)
;    IDL> print, 'Vertical clock speed (microseconds) = ', clock_speed
;
; MODIFICATION HISTORY:
;    Written by:    Mark Rivers, 11/4/97
;-

openr, lun, /get, file

header = bytarr(4100)
readu, lun, header

; Get the image size from the header

nx = fix(header, 42)
ny = fix(header, 656)
nframes = long(header, 1446)
type = fix(header, 108)
case type of
    0: data = fltarr(nx, ny, nframes)
    1: data = lonarr(nx, ny, nframes)
    2: data = lonarr(nx, ny, nframes)
    3: data = intarr(nx, ny, nframes)
    default: message, 'Unknown data type'
endcase

offset = 3000
xcal = { $
    offset: double(header, offset), $
    factor: double(header, offset+8), $
    current_unit: byte(header, offset+16), $
    reserved1: byte(header, offset+17), $
    string1: byte(header, offset+18, 40), $
    reserved2: byte(header, offset+58, 40), $
    calib_valid: byte(header, offset+98), $
    input_unit: byte(header, offset+99), $
    polynom_unit: byte(header, offset+100), $
    polynom_order: byte(header, offset+101), $
    calib_count: byte(header, offset+102), $
    pixel_pos: double(header, offset+103, 10), $
    calib_value: double(header, offset+183, 10), $
    polynom_coeff: double(header, offset+263, 6), $
    laser_position: double(header, offset+311), $
    reserved3: byte(header, offset+319), $
    new_calib_flag: byte(header, offset+320), $
    calib_label: byte(header, offset+321, 81), $
    expansion: byte(header, offset+402, 87) $
}

offset = 3489
ycal = { $
    offset: double(header, offset), $
    factor: double(header, offset+8), $
    current_unit: byte(header, offset+16), $
    reserved1: byte(header, offset+17), $
    string1: byte(header, offset+18, 40), $
    reserved2: byte(header, offset+58, 40), $
    calib_valid: byte(header, offset+98), $
    input_unit: byte(header, offset+99), $
    polynom_unit: byte(header, offset+100), $
    polynom_order: byte(header, offset+101), $
    calib_count: byte(header, offset+102), $
    pixel_pos: double(header, offset+103, 10), $
    calib_value: double(header, offset+183, 10), $
    polynom_coeff: double(header, offset+263, 6), $
    laser_position: double(header, offset+311), $
    reserved3: byte(header, offset+319), $
    new_calib_flag: byte(header, offset+320), $
    calib_label: byte(header, offset+321, 81), $
    expansion: byte(header, offset+402, 87) $
}

x_calibration = poly(findgen(nx), xcal.polynom_coeff(0:xcal.polynom_order))
y_calibration = poly(findgen(ny), ycal.polynom_coeff(0:ycal.polynom_order))

readu, lun, data
data = reform(data) ; Eliminate trailing dimensions if 1
free_lun, lun
end