pro mkmulti, input_file, help=help, norefsub=norefsub, numsamrd=numsamrd, numutr=numutr, frmtime=frmtime, output_file=output_file, mimic_multiaccum=mimic_multiaccum, sca_type=sca_type

        ;; NAME:
        ;;       pro mkmulti.pro
        ;;
        ;; PURPOSE:
        ;;       This program calibrates MULTIACCUM exposures. Calibration
        ;;       includes fitting the slope and performing reference pixel correction.
        ;;
        ;; CALLING SEQUENCE:
        ;;       mkmulti, input_file
        ;;
        ;; INPUTS:
        ;;       input_file  - A FITS data cube. The file must use IDTL header keywords.
        ;;       numutr      - Number of up the ramp sample groups.
        ;;       frmtime     - Time required to sample one frame in seconds
        ;;       output_file - The name of the output file
        ;;                     sca_type    - Set this to specify the type of
        ;;                                   SCA on the command line
        ;;       
        ;;
        ;; KEYWORD PARAMETERS:
        ;;       help             - Print a help message
        ;;       norefsub         - Do not perform reference pixel correction
        ;;       mimic_multiaccum - Mimic MULTIACCUM readout. This updates
        ;;                          NUMSAMRD and NUMUTR to the values supplied
        ;;                          on the command line (if supplied) and
        ;;                          adjusts the exposure time.
        ;;
        ;; EXAMPLE
        ;;
        ;;      mkmulti, my_multiaccum_image, my_resultant_image
        ;;
        ;;
        ;; REFERENCE:
        ;;
        ;;
        ;; MODIFICATION HISTORY:
        ;;       Written by:
        ;;       Modified:    B.J. Rauscher, IDTL, May 6, 2002
        ;;                    Initial release
        ;;
        ;;                    B.J. Rauscher, 21 January 2003
        ;;                    --- Major revisions. Automated calibration of general images
        ;;
        ;;                    B.J. Rauscher, 17 February 2003
        ;;                    --- Read large fits files by slice. This avoids the problem
        ;;                        of running out of memory on large files.
        ;;                    --- Add ability to mimic multiaccum exposures
        ;;
        ;;                    B.J.Rauscher, 8 April 2003
        ;;                    --  Add on the fly reference pixel correction
        ;;
        ;;                    B.J. Rauscher, 15 April 2003
        ;;                    --  Added sca_type parameter to command line
        ;;-
        ;;


        ;; This version of mkmulti
        version = '2.2'

        ;; Print a help message if requested
        if (not keyword_set(help)) then help = 0
        if (help) then begin
            print, 'MKMULTI version '+strtrim(string(version),2)
            stop, 'USAGE: mkmulti, input_file [, /help, /norefsub, numsamrd=, numutr=, frmtime=, sca_type=, /mimic_multiaccum]'
        endif

        ;; Set keyword defauts
        if (not keyword_set(norefsub)) then norefsub = 0; Default is to perform reference pixel correction
        if (not keyword_set(mimic_multiaccum)) then mimic_multiaccum = 0; Default is not to mimic
                                                                        ; multiaccum
        ;; Open the FITS image and read the header. We will be reading
        ;; the FITS file slice-by-slice. Close it now so that the next
        ;; slice will work
        fits_open, input_file, fcb
        fits_read, fcb, data, header, /header_only
        fits_close, fcb

        ;; Load variables using header information
        print, 'Reading header information'
        xsize = sxpar(header, 'NAXIS1', 'Keyword NAXIS1 not found') ; x-dimension of image
        ysize = sxpar(header, 'NAXIS2', 'Keyword NAXIS2 not found') ; y-dimension of image
        zsize = sxpar(header, 'NAXIS3', 'Keyword NAXIS3 not found') ; z-dimension of image
        exp_time = sxpar(header,'EXP_TIME', 'Keyword EXP_TIME not found')
        
        ;; If user has not specified the sca_type, get it from the header
        if (not keyword_set(sca_type)) then begin
            sca_type = '' ; Create an empty string to hold the SCA type
            getscatype, header, sca_type
        endif


        ;; Try to get the following variables from the command
        ;; line. If not found there, try to get them from the header

        ;; NUMSAMRD is the number of samples per group
        if (not keyword_set(numsamrd)) then $
          numsamrd = sxpar(header, 'NUMSAMRD', 'Keyword NUMSAMRD not found')

        ;; FRMTIME is the time required to read one frame
        ;; This is used only when mimicing MULTIACCUM
        if (mimic_multiaccum) then begin
            if (not keyword_set(frmtime)) then begin
                frmtime = sxpar(header, 'FRMTIME', 'Keyword FRMTIME not found')
            endif
        endif

        ;; NUMUTR is the number of groups going up the ramp
        if (not keyword_set(numutr)) then $
            numutr = sxpar(header, 'NUMUTR', 'Keyword NUMUTR not found')

        ;; Set up to mimic MULTIACCUM readout if desired
        if (mimic_multiaccum) then begin

            ;; MULTIACCUM exposure time is the difference between the
            ;; first sample of the 0th group and the first sample of
            ;; the last group. This will in general be shorter than
            ;; the SUTR exposure time.
            exp_time = numsamrd * (numutr - 1) * frmtime

            ;; Update header fields
            sxdelpar, header, 'NUMUTR'
            sxdelpar, header, 'NUMSAMRD'
            sxdelpar, header, 'EXP_TIME'
            sxaddpar, header, 'NUMUTR', numutr
            sxaddpar, header, 'NUMSAMRD', numsamrd
            sxaddpar, header, 'EXP_TIME', exp_time

        endif


        ;; Calculate the time between groups of samples
        group_time = exp_time / (numutr - 1)         ; Time required to read a group

        ;; Create and initialize an array to hold the time of the read
        ;; We only need one of these if we update it on the fly.
        time_of_read = dblarr(xsize, ysize)


        ;; To keep it simple, use the standard formula from Numerical
        ;; Recipes.
        ;;
        ;; slope = N * S[x_i * y_i] - S[x_i] * S[y_i] / ( N * S[x_i^2] -
        ;;         S[x_i]^2 )

        ;; First initialize all sums that we will need. To keep it
        ;; simple, I'll keep some redundant sums (e.g., the
        ;; information already is available to the program)
        N = 0
        Sxy = dblarr(xsize, ysize)
        Sx = dblarr(xsize, ysize)
        Sy = dblarr(xsize, ysize)
        Sxx = dblarr(xsize, ysize)

        frame0 = fltarr(xsize, ysize); Need this to subtract residual bias during reference pixel correction.

        ;; Compute sums
        print, 'Fitting least squares'
        for i = 0, numutr - 1 do begin

            ;; Update the time of read for the current slice
            time_of_read[*,*] = i * group_time

            ;; Average current group of samples
            group_mean = dblarr(xsize, ysize); Initialize an array to
                                             ; hold running sum
            for j = 0, numsamrd -1 do begin
                ;; Print status information
                current_slice = i * numsamrd + j
                print, 'Reading in frame # '+strtrim(string(current_slice),2)
                ;; Get data for the current slice
                i1 = (xsize * ysize) * numsamrd * i + (xsize * ysize) * j
                i2 = i1 + xsize * ysize - 1

                ;; Unfortunately, it seems like you have to open the
                ;; file every read. I tried leaving it open, but this
                ;; caused errors.
                fits_open, input_file, fcb
                fits_read, fcb, data, header, first=i1, last=i2, /data_only
                fits_close, fcb

                data = reform(data, xsize, ysize, /overwrite)

                ;; Do reference pixel correction if requested frame-by-frame
                if (not norefsub) then begin
                    if (i eq 0) then begin

                        ;; Update FITS header and get necessary information
                        refsub_version = ' ' ; Create a string to hold refsub version
                        refsub, version=refsub_version
                        sxaddhist, 'Processed by refsub_'+refsub_version, header
                        sxaddhist, 'Refsub using Spatial Averaging', header

                        ;; Save the 0th frame away to allow for subtracting residual bias
                        frame0 = data

                    endif else begin

                        data = data - frame0 ;Subtract residual bias
                        refsub, data, sca=sca_type; Perform reference pixel correction
                        data = data + frame0 ; Add residual bias back in

                    endelse
                endif

                ;; Update group mean
                group_mean = group_mean + data

            endfor
            group_mean[*,*] = group_mean[*,*] / numsamrd ; Take the average

            ;; Compute sums for least square fit line
            N = N + 1
            Sxy = Sxy + time_of_read[*,*] * group_mean[*,*]
            Sx = Sx + time_of_read[*,*]
            Sy = Sy + group_mean[*,*]
            Sxx = Sxx + time_of_read[*,*] * time_of_read[*,*]

            ;;;print, time_of_read[0,0], group_mean[0,0]

        endfor

        ;; Output file should be in ADU, not ADU/s
        signal = fltarr(xsize, ysize)
        signal = exp_time * (N * Sxy - Sx * Sy) / ( N * Sxx - Sx * Sx )

        ;; Update FITS header for output file
        sxaddhist, 'Calibrated using mkmulti_'+version, header

        ;; Perform reference pixel correction, if requested, using
        ;; Spatial Averaging technique
        if (not norefsub) then begin
            refsub, signal, sca=sca_type; This over-writes signal
        endif

        ;; Create output filename if one is not supplied on command line
        if (not keyword_set(output_file)) then $
          output_file = strmid(input_file, 0, strpos(input_file, '.fits', /REVERSE_SEARCH))+'.cal.fits'

        ;; Write the FITS file
        print, 'Writing Fits file'
        fits_write, output_file, signal, header

end