Viewing contents of file '../idllib/iuedac/iuelib/pro/extfits.pro'
;******************************************************************************
;+
;*NAME:
;
;  	EXTFITS         8/30/90
;  
;*CLASS:
;
;	File Conversion 
;  
;*CATEGORY:
;
;	IUESIPS
;  
;*PURPOSE:
;
;       Converts IUESIPS extracted files (i.e., MEHI & MELO files) to FITS 
;       files. MODE = 1 (default) creates a binary table file containing 
;	order number (M), number of good data points (NPT), wavelength vector
;       (W), gross (G), background (B), net (N), and abnet (A) flux vectors and 
;       epsilons (E) (written in that order) for each spectral order requested.
;       MODE = 2 is similiar except wavelengths are linearized so W (above) is 
;       replaced by the startig wavelength (W0) and wavelength increment 
;       (DELTAW) in the binary table. In MODE = 3, the wavelengths are 
;       linearized and stored as FITS keywords, and a single abnet vector is 
;       written as a primary array FITS file (this is the IRAF-compatible mode).
;  
;*CALLING SEQUENCE:
;
;  	EXTFITS,IMAGET,ORDER,mode,header
;  
;*PARAMETERS:
;
;	IMAGET	(REQ) (I) (0) (S)
;		Input IUESIPS file name (without .dat extension).
;
;       ORDER   (REQ) (I) (0,1) (B I)
;               Order number(s) to be extracted.  For low dispersion images,
;		ORDER should be set to 1. For high dispersion, ORDER can
;               be a vector with up to 60 entries (not recommended).
;   
;       MODE 	(OPT) (I) (0) (B I)
;               MODE = 1 (default) - creates a binary table file with
;		   M,NPT and vectors W,G,B,N,A,E (written in that order) for
;		   each spectral order requested.
;               MODE = 2 - wavelengths are linearized so W (above) is replaced
;		   by W0 and DELTAW in the 3D table.
;               MODE = 3 - wavelengths are linearized, and a single abnet
;		   vector is written as a primary array (this is the
;		   IRAF-compatible mode)
;
;       HEADER  (OPT) (O) (1) (S)
;               If specified, FITS header is output as a string array.
;
;*SYSTEM VARIABLES USED:
;
;	!iuer.dat
;
;*INTERACTIVE INPUT:
;
;	none
;  
;*SUBROUTINES CALLED:
;
;    	PARCHECK
;       PCHECK
;       DECOMPOSE
;       GETIUE
;       ADDPAR  - adds a keyword to FITS header
;       KEYGEN  - converts header to FITS keywords
;       FITSLAB - converts IUE label to FITS HISTORY keywords
;       TABINV
;	TRANS_BYTES
;	PLATFORM
;  
;*FILES USED:
;  
;        MAIN_HD.TXT - input file: basic main fits header
;        3DMEHI.TXT  - input file: basic FITS extension header for mode 1
;        3DMELO.TXT  - input file: basic FITS extension header for mode 1 (low)
;        3DMEHI2.TXT - input file: basic FITS extension header for mode 2
;        3DMELO2.TXT - input file: basic FITS extension header for mode 2 (low)
;        IMAGET.FIT  - output file: FITS data with header (currently written 
;                      to same directory as input file)
;  
;*SIDE EFFECTS:
;  
;*RESTRICTIONS:
;
;       For modes = 2 & 3, input file wavelengths must be monotonic.
;       For mode = 3 and if more than one order is specified, only the
;          first one listed is used.
;
;	Fits formatted files should not be used as input files.
;
;*PROCEDURE:
;
;       Calls KEYGEN to convert header record entries to FITS keywords 
;       and stores them in a string array containing 80-byte rows.
;       For mode 1, reads 3DMELO.TXT or 3DMEHI.TXT for table extension
;       header.  Additional keywords are created from H vector using
;       ADDPAR.  Vectors (i.e., W,G,B,N,A,E) are then extracted using GETIUE.
;       For mode 2, 3DMELO2.TXT or 3DMEHI2.TXT is used for extension
;       header.  GETIUE is used as above.  The wavelength array is checked 
;       for monotonicity, and, if monotonic, the remaining vectors are 
;       resampled to linear wavelengths (i.e. w0 = min(w) and deltaw =
;       max(w)-min(w)/npts).  Flux values are derived using simple linear 
;       interpolation.  Data vectors are then converted to proper format 
;       and written to disk.  For mode =3, the abnet (or rnet) vector is 
;       extracted and written to a primary array file after values are 
;       resampled to produce linear wavelengths.
;  
;*I_HELP nn:
;  
;*EXAMPLES:
;
;	extfits,'swp32199l',1,2            ; convert MELO file to FITS 3D table
;					    after linearizing wavelengths
;       to read above file:
;           ifitsrd,'swp32199l.fit',1,hd,ext,m,npt,w0,dw,g,b,n,a,e
;           w = w0 + findgen(npt)*dw       ; wavelength array
;           plot,w,a                       ; plot data
;     
;
;	extfits,'lwp3415h',[90,91,92],1,h  ; extract orders 90-92 and store as
;					     a 3D table with non-linear W
;					     keywords are also written as
;					     string array H.
;       to read order 91 (2nd row) data from above file:
;           ifitsrd,'lwp3415h.fit',2,hd,ext,m,npt,w,g,b,net,rnet,eps
;           w = w(0:npt-1)                 ; remove padded zeroes
;           f = rnet(0:npt-1)              ; remove padded zeroes
;           plot,w,f                       ; plot data (note f=rnet not abnet
;                                            flux for early MEHI file)
;
;
;       extfits,'lwr15248llg',1,3          ; write just abnet vector to primary
;                                            array FITS file for compatibility
;                                            with IRAF
;       to read data in above example:
;           ifitsrd,'lwr15248llg.fit',-1,hd,ext,flux
;           plot,flux
;
;*NOTES:
;        
;       In low dispersion, the number of non-zero points in each vector is
;	equal to the number of extracted data points.  In high dispersion,
;       all vectors are padded with zeroes to have 1000 elements. The number
;       of good points is stored in the 2nd field of each row. (See examples
;       above for how to remove padded zeroes.)
;
;	GETIUE and FITSLAB check the file format of the .dat and .lab input
;	files.  If they are in fits format, a retall is performed.
;
;       In low dispersion, the flux vector stored in the output FITS file 
;       will have units of ergs/cm^2/A representing time-integrated fluxes. 
;       This is identical to how data is stored in the IUESIPS melo files. 
;       IUEDAC routines however (such as IUESPEC) produce fluxes in units of 
;       ergs/cm^2/A/sec by dividing the IUESIPS fluxes by the user-specified 
;       exposure time. 
;
;       In high dispersion, the absolutely-calibrated data will NOT be
;       included IF the IUESIPS file was processed before ~ 1990. 
;       For these older files, EXTFITS writes out the ripple-corrected net flux
;       vector (RNET) instead of the the ABNET vector.
;
;       Additional differences between EXTFITS fluxes and IUESPEC fluxes 
;       include the following:
;         - IUESPEC may apply additional corrections such as sensitivity
;           degradation correction and THDA sensitivity corrections.
;         - EXTFITS wavelength coverage is slightly smaller (i.e. fluxes
;           at end points are set to 0).
;         - resampling in EXTFITS may cause slight changes in fluxes.
;         - The interpolation used by IUESPEC is slightly different than
;           that used by IUESIPS causing additional small differences in flux.
;
;*MODIFICATION HISTORY:
;
;      	Written by R. Thompson 8/30/90
;	 4-16-91 PJL removed vtos
;        5-17-91 GRA added vers to DECOMPOSE call statement
;        5-17-91 GRA added branch for VMS file name processing
;        7-24-91 rwt rename temporary variable A to ATMP to avoid
;                setting abnet vector to 0 when mode=2, and replace vtos
;                with dtrans for compatibility with ultrix, vms, & dos.
;        8-5-91 rwt allow version numbers (for vms), make all string
;               keywords at least 8 characters, correct numrec for multiple
;               orders, and correct starting index value (st).
;        8-7-91 rwt correct values for numrec and index parameters (again).
;       11/7/91 gra defined cpupar = !version.arch for trans_bytes
;               parameter; tested on sun, vax, and dec.
;      11/12/91 gra removed h(3) test for vector length and padded
;               all vectors to 1000 elements; cleaned up; tested.
;      11/14/91 gra padded output row for high dispersion orders to
;               a multiple of 2880 bytes; each order begins on a file
;               record boundry; tested on sun, vax, dec.
;       3/06/92 rwt add /none to openw statements to prevent vms carriage
;               control commands from being added to output data.
;       3/10/92 rwt remove change made on 11/14/91 to comply with proposed
;               binary table proposal.
;       4/10/92 rwt check for valid order numbers
;	2/24/93 PJL added version.os branch to handle Version 3 UNIX/Ultrix
;		problem with the /none keyword
;       9/14/93 rwt add mode=3 for IRAF-compatibility
;      10/07/93 rwt correct mode=3 when converting < 1 full data record 
;      11/29/93 PJL update prolog
;       2/17/94 rwt send output to fname + '.fit' not prefix + '.fit'.
;               This will write FITS file to working directory not the 
;               directory of the input file. Also updated prolog, and commented
;               out code to create the .hhh file.
;       6/07/94 rwt fix bug when multiple rows are written, write header
;               after data is read (for possibly modifying keywords), properly
;               write RNET rather than ABNET vector for early MEHI files and 
;               expand prolog information.
;       8 Jun 94 PJL replaced !version with PLATFORM
;       2 Sep 94 replace IUER_DAT with !iuer.dat
;       2 Apr 95 RWT subtract 1 in denominator of delta w calculation
;-
;******************************************************************************
 pro extfits,imaget,order,mode,header
;
 npar = n_params(0)
 if (npar eq 0) then begin
    print,'EXTFITS,IMAGET,ORDER,mode,header'
    print,'   mode=1 - binary table containing M,NPT,W,G,B,N,A,E (default)'
    print,'   mode=2 - binary table as above with linear wavelengths'
    print,'            (i.e., M,NPT,W0,DW,G,B,N,A,E)'
    print,'   mode=3 - primary array file of abnet fluxes, w/linear wavelengths'
    print,'            stored as keywords'
    retall
 endif  ; npar eq 0
 parcheck,npar,[2,3,4],'EXTFITS'
 pcheck,imaget,1,100,1000
 pcheck,order,2,110,0111
;
 idfmt = ' '
 platform,dummy,idfmt=idfmt,sysos=sys
;
; process input parameters & initialize variables
;
 if (npar lt 3) then mode = 1         ; mode = 1 is default
 decompose,imaget,disk,account,fname,ext,vers   ; parse file name
 prefix = disk + account + fname
 s = size(order)
 if (s(0) eq 0) then begin            ; make order number a vector
    m = intarr(1)
    m(0) = fix(order)
 endif else m = fix(order)
 if (n_elements(m) gt 60) then begin
    print,'Number of requested orders truncated to 60'
    m = m(0:59)
 endif  ; n_elements(m) gt 60
 if (m(0) eq 1) then m = m(0) else $
 if ( (max(m) gt 125) or (min(m) lt 66) ) then begin
    print,'Invalid order number(s) specified, EXTFITS returning'
    retall
 endif  ; (max(m) gt 125) or (min(m) lt 66)
;
; determine value for cpupar based on current cpu architecture
;
 cpupar = 0
;;; case !version.arch of
;;;   'sparc':  cpupar = 0
;;;   'vax':    cpupar = 3
;;;   'mipsel': cpupar = 5
;;;   '386':    cpupar = 5
;;;   else:  begin
;;;             print,"!version.arch not equal 'sparc', 'vax', 'mispel', " +  $
;;;                "nor '386'."
;;;             print,'ACTION:  retall
;;;             retall
;;;          end  ; !version.arch else
;;; endcase  ; !version.arch
 case (idfmt) of
   'ieeebig':  cpupar = 0
   'vaxformat':  cpupar = 3
   'ieeelittle':  cpupar = 5
 endcase  ; idfmt
;
; call getiue to get h vector and data if mode=3
; also determine number of extracted data points
;
 im = prefix + '.dat' + vers
 ord = m(0)
 abst = (ord eq 1)              ; distinguish low & high dispersion
 getiue,im,4,ord,h,w,f,e,unit
 type = 4 + (h(5) eq 7)         ; if h(5) = 7, mehi file has abscal data
 if (type eq 5) then getiue,im,type,ord,h,w,f,e,unit
 ;
 if (mode eq 3) then r = n_elements(w) else r = 1000  ; fixed vector length
 if (ord eq 1) then begin
   m = intarr(1) + 1
   r = n_elements(w) 
 endif 
;
 totrec = 0                         ; parameter for total # of records
 mnum = n_elements(m)               ; array of order numbers
 npts = intarr(mnum)                ; array for # of extracted pts. per order
 ehead = strarr(72)                 ; = two records of keywords
 a = ''
;
; read appropriate file for generating extension header keywords, and
; determine number of bytes per extracted order & other record parameters
; 
 da = !iuer.dat 
 if (mode eq 1) then begin               ; mode = 1
    ntotr = 2 + 2 + r*6*4                ; bytes per extracted order
    if (ord eq 1) then begin             ; low dispersion
       fil = da + '3dmelo.txt'           ; text file of default keywords
    endif else begin                     ; high dispersion
       fil = da + '3dmehi.txt'          
    endelse  ; ord=1
 endif ; mode=1
if (mode eq 2) then begin                        ; mode = 2
    ntotr = 2 + 2 + 2*4 + r*5*4          ; bytes per extracted order
    if (ord eq 1) then begin             ; low dipsersion
       fil = da + '3dmelo2.txt'          ; text file of default keywords
    endif else begin                     ; high dispersion
       fil = da + '3dmehi2.txt'
    endelse  ; ord=1
 endif ; mode=2 
if (mode eq 3) then ntotr = r*4
;
 left = ntotr mod 2880                   ; leftover bytes
 nrec = fix(ntotr/2880)                  ; # of full fits records per order
 record = bytarr(ntotr)                  ; byte array for 1 spectral order
 buff = bytarr(2880)                     ; buffer for storing leftovers
;
; convert h vector entries to fits keywords
;
 keygen,h,mhead                          ; create keywords
 addpar,mhead,'FILENAME',fname           ; add file name
;
; create basic extension header 
;
if (mode lt 3) then begin
   openr,lu1,fil,/get_lun
   for n=0,2 do readf,lu1,a                ; skip 1st 3 lines of text file
   i = 0                                   
   while not eof(lu1) do begin
     tempa = string(bytarr(80) + 32b)
     readf,lu1,a
     strput,tempa,a
     ehead(i) = tempa                     ; store in string array
     i = i + 1
   endwhile  ; not (eof)
   free_lun,lu1
   if (i le 71) then ehead(i:71) = string(bytarr(80) + 32b)
;
; add keywords to main & or extension headers
;
   tmp = strtrim(string(r),2)+'E'
   form = string(bytarr(8) + 32b)
   strput,form,tmp,0                       ; embed data in 8 character string
   addpar,ehead,'NAXIS1',ntotr
   addpar,ehead,'NAXIS2',mnum
endif else ehead = strarr(1)               ; no ext head for mode=3
;
 if (mode eq 1) then begin 
    addpar,ehead,'TFIELDS',8
    addpar,ehead,'TFORM1','1I      '
    addpar,ehead,'TFORM2','1I      '
    addpar,ehead,'TFORM3',form
    addpar,ehead,'TFORM4',form
    addpar,ehead,'TFORM5',form
    addpar,ehead,'TFORM6',form
    addpar,ehead,'TFORM7',form
    addpar,ehead,'TFORM8',form
    if (type eq 4) and (abst eq 0) then begin
        addpar,ehead,'TTYPE7','RNET    ',' ripple corrected net flux'
        addpar,ehead,'TUNIT7','FN      ',' flux numbers'
    endif
 endif  ; mode eq 1
;
 if (mode eq 2) then begin
    addpar,ehead,'TFIELDS',9
    addpar,ehead,'TFORM1','1I      '
    addpar,ehead,'TFORM2','1I      '
    addpar,ehead,'TFORM3','1E      '
    addpar,ehead,'TFORM4','1E      '
    addpar,ehead,'TFORM5',form
    addpar,ehead,'TFORM6',form
    addpar,ehead,'TFORM7',form
    addpar,ehead,'TFORM8',form
    addpar,ehead,'TFORM9',form
    if (type eq 4) and (abst eq 0) then begin
        addpar,ehead,'TTYPE8','RNET    ',' ripple corrected net flux'
        addpar,ehead,'TUNIT8','FN      ',' flux numbers'
    endif
 endif  ; mode eq 2
;
if (mode eq 3) then begin
    ; 
    ; in this mode, data must be processed before header is completed
    ;
    ; test for wavelength monotonicity
    ;
    monoinc = 0
    is = w - shift(w,1)
    atmp = where(is gt 0)
    if (n_elements(atmp) eq (n_elements(w)-1)) then monoinc = 1
    atmp = where(is lt 0)     
    if (n_elements(atmp) eq (n_elements(w)-1)) then monoinc = 2
    if (monoinc eq 0) then begin
       print,'input wavelength array is not monotonic'
       print,'either modify input array or rerun extfits using mode = 1'
       retall
    endif  ; monoinc
    ;
    ; linearize wavelengths
    ;
    npt = n_elements(w)
    w0 = min(w)
    delta = (max(w) - w0) / (npt-1)       ; same # of pts equally-spaced
    neww = w0 + delta * indgen(npt)   ; new wavelengths
    newf = f
    tabinv,w,neww,ri                      ; index of interpolated pts.
    index = fix(ri)
    ri = ri - index
    newf = f(index+1) * ri + f(index) * (1-ri)    
    ind = where(neww,npt)
    ;
    ; convert output data to appropriate format byte arrays
    ;
    notr = npt*4
    bflux = byte(newf,0,notr)
    left = ntotr mod 2880                ; leftover bytes
    nrec = fix(ntotr/2880)               ; # of full fits records per order
    record = bytarr(ntotr)               ; byte array for 1 spectral order
    if (cpupar gt 0) then trans_bytes,bflux,4,cpupar
    record(0) = bflux
    ;
    addpar,mhead,'BITPIX',-32,' Floating point data'
    addpar,mhead,'NAXIS',1,' One-dimensional array'
    addpar,mhead,'NAXIS1',npt,' Number of elements','EXTEND'
    addpar,mhead,'ORDNUM',m(0),' Spectral order number','EXTEND'
    addpar,mhead,'EXTEND','F',' No extensions'
    addpar,mhead,'CRVAL1',neww(0),' Starting wavelength','EXTEND'
    addpar,mhead,'CDELT1',neww(1)-neww(0),' Wavelength increment','EXTEND'
    addpar,mhead,'CTYPE1','WAVELENGTH',' X-axis','EXTEND'
    if ( (type eq 5) or (abst eq 1) ) then begin
       addpar,mhead,'COMMENT','File contains absolutely-calibrated net Fluxes'
       addpar,mhead,'BUNIT','ERGS/CM2/A',' Calibrated flux','EXTEND'
    endif  ; (type eq 5) or (abst eq 1)
    if ( (type eq 4) and (abst eq 0) ) then begin
       addpar,mhead,'COMMENT','File contains ripple-corrected net Fluxes' 
       addpar,mhead,'BUNIT','FN      ',' Flux Numbers','EXTEND'
    endif  ; (type eq 4) and (abst eq 0)
 endif   ; mode eq 3
;
; append label information (if .lab file exists)
;
 nkm = n_elements(mhead)                 ; find end of main header
 last = 0
 while ( ( strpos(mhead(last),'END') ne 0 ) and (last le nkm) ) do last=last+1
 imagelab = prefix + '.lab' + vers
 on_ioerror,contin                       ; skip to contin if no label file
 file = findfile(imagelab)
 on_ioerror,null
 fitslab,imagelab,lab,flab               ; add label info
 nfl = n_elements(flab)
 ntot = last + nfl + 1
 hrec = fix(ntot / 36) + ((ntot mod 36) ne 0)  ; # of header records
 tmp = mhead(0:last)
 mhead = strarr(hrec*36)
 mhead(0:last) = tmp 
 mhead(last:last+nfl-1) = flab           ; want to write over END
 if (ntot-1 lt hrec*36-1) then mhead(ntot-1:hrec*36-1) =  $
    string(bytarr(80) + 32b)
 tmp = string(bytarr(80) + 32b)          ; add end to end of main header
 strput,tmp,'END',0
 mhead(ntot-1) = tmp
 contin:
;
; append extension header (if available) to main header 
;
 nkm = n_elements(mhead)
 nke = n_elements(ehead) 
 ntot = nkm + nke * (mode ne 3)          ; mode = 3 has no extension
 hrec = fix(ntot / 36) + ((ntot mod 36) ne 0)
 header = strarr(hrec*36)
 header(0:nkm-1) = mhead
 if (mode ne 3) then header(nkm:ntot-1) = ehead
 if (ntot lt hrec*36-1) then header(ntot:hrec*36-1) = string(bytarr(80) + 32b)
;
; use mode to determine data record format 
;
;     binary table with original wavelength array
;     convert floating point values to ansi standard
;     store in order m,npt,w,g,b,n,a,e
;     use fixed length vectors of 'r' elements
;
 offset = 0
 if (mode eq 1) then begin
    for i=0,mnum-1 do begin              ; go through each order
       wpad = fltarr(r)
       epad = wpad
       flux = fltarr(4*r)
       for nr=1,4 do begin               ; extract g,b,n,a flux vectors
          if (nr eq 4) then getiue,im,type,m(i),h,w,f,e,unit $
          else getiue,im,nr,m(i),h,w,f,e,unit
          flux(r*(nr-1)) = f
       endfor  ; nr
       epad(0) = float(e)
       wpad(0) = w
       morass = [wpad,flux,epad]         ; all floating pt. vectors
       ind = where(w ne 0)
       npts(i) = n_elements(ind)
       npt = npts(i)
       ord = m(i)
      ;
      ; convert output data to appropriate format byte arrays
      ;
       bord = byte(ord,0,2)
       bnpt = byte(npt,0,2)
       bmorass = byte(morass,0,r*6*4)
       if (cpupar gt 0) then begin
          trans_bytes,bord,2,cpupar
          trans_bytes,bnpt,2,cpupar
          trans_bytes,bmorass,4,cpupar
       endif  ; cpupar gt 0
       record(0) = bord 
       record(2) = bnpt
       record(4) = bmorass
     ;
     ; write all keywords to disk file with .fit extension
     ;
       if (i eq 0) then begin
         get_lun,lun2
         rech = assoc(lun2,bytarr(2880))
;;;         if (!version.os eq 'vms') then openw,lun2,fname+'.fit',2880,/none  $
;;;           else openw,lun2,fname+'.fit',2880
         if (sys eq 'vms') then openw,lun2,fname+'.fit',2880,/none  $
           else openw,lun2,fname+'.fit',2880
         for j=0,hrec-1 do begin                 ; write out header
           st = (j * 36)                        ; 36 keywords per record
           rech(j) = byte(header(st:st+35))     ; convert to byte array
         endfor  ; j 
         numrec = hrec                           ; cumulative record total
       endif
      ;
      ; write data records to disk file
      ;
       nrec = fix((ntotr+offset)/2880)         ; # of full records
       for ir=0,nrec-1 do begin                ; write out whole records
          off = offset * (ir eq 0)             ; keep leftover bytes
          st = (ir * 2880) - (offset * (ir ne 0)) ;start record index
          eod = st + 2879 - off                ; end record index
          buff(off) = record(st:eod)           ; store at index offset
          rech(ir+numrec) = buff               ; write 2880 byte record
       endfor  ; ir
       offset = (ntotr+offset) mod 2880        ; where to start next row
       print,'extracting order',byte(m(i)),' containing',fix(npts(i)),' points'
       numrec = numrec + nrec                  ; total record count
   ;
   ; write any remaining bytes to last data record
   ;
       if (offset ne 0) then begin         ; leftover bytes
           buff = bytarr(2880)             ; create new buffer
           buff(0) = record(eod+1:*)       ; store remaining bytes in buff
           if (i eq (mnum-1)) then begin   ; if last row of file
              rech(numrec) = buff
              numrec = numrec + 1
           endif                           ; write out leftover if last row
       endif                               ; save leftover data in buff
    endfor  ; selected orders
 endif  ; mode = 1
;
;     binary table with linear wavelengths
;     convert e to floating point vector
;     convert floating point values to ansi standard
;     store in order m,npt,w0,delta,[g0,g1,...],...[e0,e1,...]
;     
 if (mode eq 2) then begin
    for i=0,mnum-1 do begin                    ; go through each order
       flux = fltarr(5*r)
       getiue,im,1,m(i),h,w,g,e,unit
       getiue,im,2,m(i),h,w,b,e,unit
       getiue,im,3,m(i),h,w,n,e,unit
       getiue,im,type,m(i),h,w,a,e,unit
      ;
      ; test for wavelength monotonicity
      ;
       monoinc = 0
       is = w - shift(w,1)
       atmp = where(is gt 0)
       if (n_elements(atmp) eq (n_elements(w)-1)) then monoinc = 1
       atmp = where(is lt 0)     
       if (n_elements(atmp) eq (n_elements(w)-1)) then monoinc = 2
       if (monoinc eq 0) then begin
          print,'Input wavelength array is not monotonic'
          print,'either modify input array or rerun extfits using mode = 1'
          retall
       endif  ; monoinc
       e = float(e)
       npts(i) = n_elements(w)
       w0 = min(w)
       delta = (max(w) - w0) / (npts(i)-1)   ; same # of pts equally-spaced
       neww = w0 + delta * indgen(npts(i))   ; new wavelengths
       newg = g
       newb = b
       newn = n
       newa = a
       newe = e
      ;
      ; linearize wavelengths
      ;
       tabinv,w,neww,ri                      ; index of interpolated pts.
       index = fix(ri)
       ri = ri - index
       newg = g(index+1) * ri + g(index) * (1-ri)    ; interpolate fluxes
       newb = b(index+1) * ri + b(index) * (1-ri)    
       newn = n(index+1) * ri + n(index) * (1-ri)    
       newa = a(index+1) * ri + a(index) * (1-ri)    
       newe = e(index+1) < e(index)                  ; take lowest neighbor
       flux(0) = newg
       flux(r) = newb
       flux(r*2) = newn
       flux(r*3) = newa
       flux(r*4) = newe
       ind = where(neww ne 0)
       npts(i) = n_elements(ind)
       ord = m(i)
       npt = npts(i)
      ;
      ; convert output data to appropriate format byte arrays
      ;
       bord = byte(ord,0,2)
       bnpt = byte(npt,0,2)
       bw0 = byte(w0,0,4)
       bdelta = byte(delta,0,4)
       bflux = byte(flux,0,r*5*4)
       if (cpupar gt 0) then begin
          trans_bytes,bord,2,cpupar
          trans_bytes,bnpt,2,cpupar
          trans_bytes,bw0,4,cpupar
          trans_bytes,bdelta,4,cpupar
          trans_bytes,bflux,4,cpupar
       endif  ; cpupar gt 0
       record(0) = bord 
       record(2) = bnpt
       record(4) = bw0
       record(8) = bdelta
       record(12) = bflux 
     ;
     ; write all keywords to disk file with .fit extension
     ;
       if (i eq 0) then begin
         get_lun,lun2
         rech = assoc(lun2,bytarr(2880))
;;;         if (!version.os eq 'vms') then openw,lun2,fname+'.fit',2880,/none  $
;;;           else openw,lun2,fname+'.fit',2880
         if (sys eq 'vms') then openw,lun2,fname+'.fit',2880,/none  $
           else openw,lun2,fname+'.fit',2880
         for j=0,hrec-1 do begin                 ; write out header
           st = (j * 36)                        ; 36 keywords per record
           rech(j) = byte(header(st:st+35))     ; convert to byte array
         endfor  ; j 
         numrec = hrec                           ; cumulative record total
       endif
      ;
      ; write data records to disk file
      ;
       nrec = fix((ntotr+offset)/2880)         ; # of full records
       for ir=0,nrec-1 do begin                ; write out whole records
          off = offset * (ir eq 0)             ; keep leftover bytes
          st = (ir * 2880) - (offset * (ir ne 0)) ;start record index
          eod = st + 2879 - off                ; end record index
          buff(off) = record(st:eod)           ; store at index offset
          rech(ir+numrec) = buff               ; write 2880 byte record
       endfor  ; ir
       offset = (ntotr+offset) mod 2880        ; where to start next row
       print,'extracting order',byte(m(i)),' containing',fix(npts(i)),' points'
       numrec = numrec + nrec                  ; total record count
;    endfor  ; selected orders
   ;
   ; write any remaining bytes to last data record
   ; (increase record count if last row written)
   ;
       if (offset ne 0) then begin         ; leftover bytes
           buff = bytarr(2880)             ; create new buffer
           buff(0) = record(eod+1:*)       ; store remaining bytes in buff
           if (i eq (mnum-1)) then begin   ; if last row of file
              rech(numrec) = buff
              numrec = numrec + 1
           endif                           ; write out leftover if last row
       endif                               ; save leftover data in buff
     endfor  ; selected orders
endif  ; mode = 2
;
; primary array with linear wavelengths specified as keywords
; data already converted to IEEE format 
;
 if (mode eq 3) then begin
     ;
     ; write all keywords to disk file with .fit extension
     ;
       get_lun,lun2
       rech = assoc(lun2,bytarr(2880))
;;;       if (!version.os eq 'vms') then openw,lun2,fname+'.fit',2880,/none  $
;;;         else openw,lun2,fname+'.fit',2880
       if (sys eq 'vms') then openw,lun2,fname+'.fit',2880,/none  $
         else openw,lun2,fname+'.fit',2880
       for j=0,hrec-1 do begin                 ; write out header
         st = (j * 36)                        ; 36 keywords per record
         rech(j) = byte(header(st:st+35))     ; convert to byte array
       endfor  ; j 
       numrec = hrec                           ; cumulative record total
      ;
      ; write data records to disk file
      ;
       nrec = fix((ntotr+offset)/2880)         ; # of full records
       if (nrec eq 0) then begin               ; less than one record
          buff(0) = record
          rech(numrec) = buff
       endif else begin
          for ir=0,nrec-1 do begin                ; write out whole records
             off = offset * (ir eq 0)             ; keep leftover bytes
             st = (ir * 2880) - (offset * (ir ne 0)) ;start record index
             eod = st + 2879 - off                ; end record index
             buff(off) = record(st:eod)           ; store at index offset
            rech(ir+numrec) = buff               ; write 2880 byte record
          endfor  ; ir
       endelse  ; nrec eq 0
       offset = ((ntotr+offset) mod 2880) * (nrec ne 0) ; where start next row
       print,'extracting order',byte(m(0)),' containing',fix(npt),' points'
       numrec = numrec + nrec + (nrec eq 0)      ; total record count
   ;
   ; write any remaining bytes to last data record
   ;
    if (offset ne 0) then begin
       buff = bytarr(2880)
       buff(0) = record(eod+1:*)
       rech(numrec) = buff
       numrec = numrec + 1
    endif  ; offset ne 0
 endif  ; mode = 3
;
 free_lun,lun2
 free_lun,unit
 print,'Output fits file has',byte(hrec),' header records and', $
    byte(numrec-hrec), ' data records.'
;
 return  
 end  ; extfits