	PROGRAM SIR2BYTE
C
C	THIS PROGRAM ILLUSTRATES THE READING OF A FILE IN THE BYU MERS SIR
C       FILE FORMAT.  IT READS THE FILE AND
C       WRITES OUT A BYTE-ARRAY COPY OF THE IMAGE
C
C	WRITTEN BY D.LONG: MARCH 1997
C
        CHARACTER*70 FNAME,ONAME
C
C	IMAGE STORAGE VARIABLES
C
	PARAMETER (MAXSIZE=4500000)       ! MAXIMUM PIXELS
	REAL 	STVAL(MAXSIZE)            ! READ ARRAY
	BYTE	B(MAXSIZE)                ! OUTPUT BYTE ARRAY
C
C	IMAGE FILE HEADER INFORMATION
C
	CHARACTER SENSOR*40,TITLE*80
        CHARACTER TAG*40,TYPE*138,CRPROC*100,CRTIME*28
        INTEGER   NSX,NSY,IOPT,IDATATYPE,IOFF,ISCALE
        INTEGER   NHEAD,NDES,NHTYPE,LDES,NIA,IPOL,IFREQHM,ISPARE1
        INTEGER   IREGION,ITYPE,IYEAR,ISDAY,ISMIN,IEDAY,IEMIN
        REAL      XDEG,YDEG,ASCALE,BSCALE,A0,B0,ANODATA,VMIN,VMAX
C
C     OPTIONAL HEADER INFO
C
        PARAMETER (MAXI=256)
        CHARACTER DESCRIP*1024
	DIMENSION IAOPT(MAXI)
C
C	PROMPT USER FOR FILE NAME INPUT
C
	WRITE(*,*)
	WRITE(*,15)
15	FORMAT(' Enter Input SIR File Name: ',$)
	READ(*,20) FNAME
20	FORMAT(A70)
C
40	CONTINUE
	WRITE(*,*) 'SIR In= "',FNAME(1:LENGTH(FNAME)),'"'
C
C	READ SIR FILE HEADER
C
	MAXDES=LEN(DESCRIP)
	CALL READSIRHEAD(FNAME,IU,IERR,NHEAD,NDES,NHTYPE,IDATATYPE,
     *     NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE,
     *     IYEAR,ISDAY,ISMIN,IEDAY,IEMIN,IREGION,ITYPE,
     *     IPOL,IFREQHM,ISPARE1,ANODATA,VMIN,VMAX,
     *     SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,
     *     MAXDES,DESCRIP,LDES,MAXI,IAOPT,NIA)
C
C     READ IMAGE CONTENTS
C
	IF (IERR.GE.0) THEN
	   CALL READSIRF(IU,IERR,NHEAD,NHTYPE,IDATATYPE,STVAL,NSX,NSY,
     *	         IOFF,ISCALE,SMIN,SMAX,NCNT,ANODATA,VMIN,VMAX)
	ELSE
	   WRITE (*,*) '*** ERROR OPENING FILE',IERR
	ENDIF
	IF (IERR.LT.0) STOP
C
C	WRITE OUT IMAGE HEADER INFORMATION
C
	WRITE (*,*)
	WRITE (*,*) 'SIR File Header Information:'
	CALL PRINTHEAD(NHEAD,NDES,NHTYPE,IDATATYPE,
     *     NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE,
     *     IYEAR,ISDAY,ISMIN,IEDAY,IEMIN,IREGION,ITYPE,
     *     IPOL,IFREQHM,ISPARE1,ANODATA,VMIN,VMAX,
     *     SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,
     *     DESCRIP,LDES,IAOPT,NIA)
	WRITE (*,*)
C
C	DETERMINE MAX/MIN OF DATA ARRAY
C
	SMIN=1.E25
	SMAX=-1.E25
	DO NN1=0,NSY-1
		DO NN2=1,NSX
			INDEX = NN1*NSX + NN2
			IF (STVAL(INDEX).NE.ANODATA) THEN
				SMIN=MIN(SMIN,STVAL(INDEX))
				SMAX=MAX(SMAX,STVAL(INDEX))
			ENDIF
		END DO
	END DO
C
	SMIN1=SMIN
	SMAX1=SMAX
C
	WRITE (*,*)
	WRITE (*,*) 'Array min, max values:',SMIN1,SMAX1
C
    	WRITE (*,34) 'Saturation Max?'
34	FORMAT(X,'New ',A,X,$)
        READ (*,*) SMAX
	WRITE (*,34) 'Saturation Min?'
        READ (*,*) SMIN
C
C	SCALE ARRAY AND PUT INTO BYTE ARRAY
C
C     NOTE THAT SIR IMAGE DATA HAS THE ORIGIN AT THE LOWER-LEFT OF THE
C     IMAGE.  THE BYTE ARRAY HAS THE ORIGIN AT THE UPPER-LEFT.
C
	SS=(SMAX-SMIN)/255.0
	IF (SS.GT.0.0) THEN
		SS=1.0/SS
	ELSE
		SS=1.0
	ENDIF
	DO NN1=0,NSY-1
		DO NN2=1,NSX
			INDEX = NN1*NSX + NN2
			INDEX2= (NSY-NN1)*NSX + NN2
			VAL=STVAL(INDEX)
			IVAL=(VAL-SMIN)*SS
			IF (IVAL.LT.0) IVAL=0
			IF (IVAL.GT.255) IVAL=255
			IF (IVAL.GT.127) IVAL=IVAL-256
			B(INDEX2)=IVAL
		END DO
	END DO
C
C     GET OUTPUT FILE NAME FOR BYTE FILE
C
	WRITE(*,25)
25	FORMAT(' Enter BYTE Output File Name: ',$)
	READ(*,20) ONAME
	WRITE(*,*) 'BYTE Out="',ONAME(1:LENGTH(ONAME)),'"'
C
C	NOW WRITE OUT BYTE FILE
C
	WRITE(*,*) 'Output Image Size: (X x Y)',NSX,NSY
	NS=NSX*NSY
C
	OPEN(UNIT=1,FILE=ONAME,STATUS='NEW',FORM='UNFORMATTED')
C
	WRITE (1) (B(I),I=1,NS)
        CLOSE (1)
C
360	STOP
	END
C
C **************************************************************************
C
	SUBROUTINE PRINTHEAD(NHEAD,NDES,NHTYPE,IDATATYPE,
     *     NSX,NSY,XDEG,YDEG,ASCALE,BSCALE,A0,B0,IOPT,IOFF,ISCALE,
     *     IYEAR,ISDAY,ISMIN,IEDAY,IEMIN,IREGION,ITYPE,
     *     IPOL,IFREQHM,ISPARE1,ANODATA,VMIN,VMAX,
     *     SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,
     *     DESCRIP,LDES,IAOPT,NIA)
C
C	PRINTS OUT SIR FILE HEADER
C
	CHARACTER*(*) SENSOR,TITLE,TYPE,TAG,CRPROC,CRTIME,DESCRIP
	INTEGER IAOPT(1)
C
	IF (NHTYPE.LT.16) WRITE (*,*) 'Old style header',nhtype
C
 	WRITE (*,*) 'Title:  "',TITLE(1:LENGTH(TITLE)),'"'
	WRITE (*,*) 'Sensor: "',SENSOR(1:LENGTH(SENSOR)),'"'
	IF (NHTYPE.GT.16) THEN
	   WRITE (*,*) 'Type:   "',TYPE(1:LENGTH(TYPE)),'"'
	   WRITE (*,*) 'Tag:    "',TAG(1:LENGTH(TAG)),'"'
	   WRITE (*,*) 'Creator:"',CRPROC(1:LENGTH(CRPROC)),'"'
	   WRITE (*,*) 'Created:"',CRTIME(1:LENGTH(CRTIME)),'"'
	ENDIF
	IF (IOPT.EQ.-1) THEN		! IMAGE ONLY
	 WRITE (*,*) 'Rectangular image-only form: ',ioff,iscale,iopt
	 WRITE (*,*) ' x,y scale   ',ascale,bscale
	 WRITE (*,*) ' x,y offsets:',a0,b0
	 WRITE (*,*) ' x,y span:   ',xdeg,ydeg
	ELSE IF (IOPT.EQ.0) THEN	! LAT/LONG
	 WRITE (*,*) 'Lat/Lon Rectangular form: ',ioff,iscale,iopt
	 WRITE (*,*) ' X and Y scales (in pixels/deg)',ascale,bscale
	 WRITE (*,*) ' Offset values:',a0,b0
	 WRITE (*,*) ' Degrees:      ',xdeg,ydeg
	ELSE IF (IOPT.EQ.1.OR.IOPT.EQ.2) THEN	! LAMBERT
	 WRITE (*,*) 'Lambert form: ',ioff,iscale,iopt
	 WRITE (*,*) ' X and Y scales (in km/pixel)',ascale,bscale
	 WRITE (*,*) ' Lower Left Corner:',a0,b0
	 WRITE (*,*) ' Center Point:     ',xdeg,ydeg
	 IF (IOPT.EQ.2) WRITE (*,*) ' Local radius used'
	ELSE IF (IOPT.EQ.5) THEN	! POLAR STEREOGRAPHIC
	 WRITE (*,*) 'Polar Sterographic form:     ',ioff,iscale,iopt
	 WRITE (*,*) ' X and Y scales (in km/pixel)',ascale,bscale
	 WRITE (*,*) ' Lower Left Corner: (km)     ',a0,b0
	 WRITE (*,*) ' Center Lon, Lat:            ',xdeg,ydeg
	ELSE IF (IOPT.EQ.11.OR.IOPT.EQ.12) THEN	! EASE NORTH OR SOUTH
	 WRITE (*,*) 'EASE polar azimuthal form: ',ioff,iscale,iopt
	 WRITE (*,*) ' A, B scales:',ascale,bscale
	 WRITE (*,*) ' Map origin (col/row):',a0,b0
	 WRITE (*,*) ' Map center (col/row):',xdeg,ydeg
	ELSE IF (IOPT.EQ.13) THEN	! EASE CYLINDRICAL
	 WRITE (*,*) 'EASE cylindrical form: ',ioff,iscale,iopt
	 WRITE (*,*) ' A, B scales:',ascale,bscale
	 WRITE (*,*) ' Map origin (col/row):',a0,b0
	 WRITE (*,*) ' Map center (col/row):',xdeg,ydeg
	ELSE				! UNKNOWN TRANSFORMATION
	 WRITE (*,*) '*** Unrecognized SIR file form ***',ioff,iscale,iopt
	 WRITE (*,*) ' a,b scale   ',ascale,bscale
	 WRITE (*,*) ' a,b offsets:',a0,b0
	 WRITE (*,*) ' x,y span:   ',xdeg,ydeg
	ENDIF
	WRITE (*,*) ' Pixels:       ',nsx,' x ',nsy
	WRITE (*,*) ' Year:  ',iyear,'  Region: ',iregion
	WRITE (*,*) ' Start: ',isday,ismin
	WRITE (*,*) ' End:   ',ieday,iemin
	WRITE (*,*) ' Type:  ',itype,'  Transform option: ',iopt
	IF (NHTYPE.GT.16) THEN
	   WRITE (*,*) ' Frequency:  ',ifreqhm*0.1,' MHz  Pol: ',ipol
	   WRITE (*,*) ' Data:  ',idatatype,'  Header: ',nhead,nhtype,ndes
	   WRITE (*,*) ' Offset:',ioff,'  Scale:',iscale
	   WRITE (*,*) ' Nodata:',anodata,'  Vn,Vx:',vmin,vmax
	   IF (LDES.GT.0) WRITE (*,*) 'Description: "',DESCRIP(1:LDES),'"'
	   IF (NIA.GT.0) THEN
	      DO I=1,NIA
		 WRITE (*,*) ' Aopt: ',i,' ',iaopt(i)
	      END DO
	   ENDIF
	ENDIF
	RETURN
	END
C
C **************************************************************************
C
	INTEGER FUNCTION LENGTH(LINE)
C
C	COMPUTES THE LENGTH (NON-BLANK) OF A STRING
C
	CHARACTER*(*) LINE
	N=LEN(LINE)
	DO I=N,1,-1
		IF (LINE(I:I).NE.' '.AND.LINE(I:I).NE.CHAR(0)) THEN
			LENGTH=I
			GOTO 10
		ENDIF
	END DO
	LENGTH=1
10	RETURN
	END
C
C **************************************************************************
C
C STANDARD SIR IMAGE FORMAT READ/WRITE ROUTINES
C
C **************************************************************************
C
      subroutine readsirhead(fname,iu,ierr,nhead,ndes,nhtype,idatatype,
     *     nsx,nsy,xdeg,ydeg,ascale,bscale,a0,b0,iopt,ioff,iscale,
     *     iyear,isday,ismin,ieday,iemin,iregion,itype,
     *     ipol,ifreqhm,ispare1,anodata,vmin,vmax,
     *     sensor,title,type,tag,crproc,crtime,
     *     maxdes,descrip,ldes,maxi,iaopt,nia)
c
c     Version 2 standard SIRF image format header read routine
c
c     written by D.G. Long  14 Mar 1997
c
c inputs:
c     fname    file name to read from
c     iu       fortran logical unit to use
c     maxdes   maximum number of characters in descrip
c     maxi     maximum number of integers in iaopt
c
c outputs:
c     ierr   file read error code
c             set to 0 for successful read of header
c             set to -1 for open error
c             set to -2 for header read error
c     nhead  number of 512 byte blocks in header
c     nhtype file header type code
c     idatatype file data type code (0,2=i*2,1=i*1,4=i*4) def=i*2
c     nsx,nsy image dimensions
c     xdeg..b0 projection factors
c     iopt    projection type
c     ioff,iscale image storage scale factors
c     iyear,isday,ismin  year,day,minute of start of data used in image
c     ieday,iemin day,minute of end of data used in image
c     iregion study region id number
c     itype   image type code (=0 none or unknown)
c             typically: 1's digit is polarization (0=H, 1=V)
c             10's digit and higher is frequency in hundreds of MHz
c     ispare* spare integers
c     sensor  20 character string: describes sensor
c     title   80 character string: image title
c     type    80 character string: description of image type
c     tag     creator byline
c     crproc  100 character string: name and version of creating program
c     crtime  40 character string: date and time of creation
c     descrip variable length character string: description and annontation
c     ldes    number of characters in descrip
c     iaopt   optional integers
c       by convention, first value iaopt is a code telling how to interpret
c       the rest of the array if nia>0
c     nia     number of integers
c
c     SIR files have headers which are multiples of 512 bytes followed by
c     image data stored as 2 byte integers zero padded to be a multiple
c     of 512 bytes.
c
	character*(*) 	fname
	integer		iu,nhead,nhtype,idatatype
	integer 	nsx,nsy,iopt
	real		xdeg,ydeg,ascale,bscale,a0,b0
	integer		ioff,iscale,iyear,isday,ismin,ieday,iemin
        integer		iregion,itype,ipol,ifreqhm,ispare1
	real            anodata,vmin,vmax
	character*(*)	sensor,title,type,tag,crproc,crtime,descrip
	integer		iaopt(*)
c
c     working variables
c
	integer*2       i2(2)
	real            f2
	equivalence (i2,f2)      ! used in i*2 to real mapping
	integer*2	temp(256)
	character	sens1*41,tit1*81,type1*139
        character       tag1*41,crp1*101,crt1*29
c
c     note: the bswap and vms flags should be set in the block data
c
	logical bswap,vms
	common /vaxorunix/bswap,vms
	save /vaxorunix/
c
c	header mapping (see code for scaling of <= variables)
c       note that character strings do not have to be null terminated
c
c	temp(1)	= nsx			! pixels in x direction
c	temp(2) = nsy			! pixels in y direction
c       temp(3) <= xdeg			! span of deg in x
c	temp(4) <= ydeg			! span of deg in y
c	temp(5) = nhtype                ! header type (old<15)
c	temp(6) <= ascale               ! x scaling
c	temp(7) <= bscale               ! y scaling
c	temp(8) <= a0                   ! x (lon) origin
c         note: longitudes should be in the range -180 to 180
c	temp(9) <= b0                   ! y (lat) origin
c	temp(10) = ioff			! offset to be added to scaled val
c	temp(11) = iscale		! scale factor ival=(val-ioff)/iscale
c	temp(12) = iyear		! year for data used
c	temp(13) = isday		! starting JD
c	temp(14) = ismin		! time of day for first data (in min)
c	temp(15) = ieday		! ending JD
c	temp(16) = iemin		! time of day for last data (in min)
c	temp(17) = iopt			! projection type
c					!  -1 = no projection, image only
c					!   0 = rectalinear lat/lon
c					!   1 = lambert equal area
c					!   2 = lambert equal area (local rad)
c					!   5 = polar stereographic
c					!  11 = EASE north equal area grid
c					!  12 = EASE south equal area grid
c					!  13 = EASE cylindrical grid
c	temp(18) = iregion		! region id code
c	temp(19) = itype		! image type code
c                                       ! standard values: 0=unknown or n/a
c                                       ! 1 = scatterometer A (dB)
c                                       ! 2 = scatterometer B (dB/deg)
c                                       ! 3 = radiometer Tb (K)
c                                       ! 9 = topography (m)
c	temp(20)-temp(39) 40 chars of sensor
c       temp(40) = 0 
c       temp(41) = nhead                ! number of 512 byte header blocks
c       temp(42) = ndes                 ! number of 512 byte blocks description
c       temp(43) = ldes                 ! number of bytes of description
c       temp(44) = nia                  ! number of optional integers
c       temp(45) = ipol                 ! polarization (0=n/a,1=H,2=V)
c       temp(46) = ifreqhm              ! frequency in 100's MHz (0 if n/a)
c       temp(47) = ispare1              ! spare
c       temp(48) = idatatype            ! data type code 0,2=i*2,1=i*1,4=f
c
c       the value of idata type determines how data is stored and how
c       anodata, vmin, and vmax are stored.
c
c       if idatatype = 1 data is stored as bytes and minv=128
c       if idatatype = 2 data is stored as 2 byte integers and minv=32766
c       if idatatype = 4 data is stored as IEEE floating point
c
c       if idatatype = 1,2 anodata,vmin,vmax are stored as 2 byte integers
c         in temp(49)..temp(51)  minv, ioff and iscal used to convert
c         integers or bytes into floating point values
c         nodata, vmin, and vmax must be representable with ioff and iscale
c            temp(*) = (value-ioff)*iscale-minv
c            value = float(temp(*)+minv)/float(iscale)+ioff
c       idatatype=2 is considered the SIR standard format
c
c       if idatatype = f anodata,vmin,vmax are stored as floating points
c         in temp(52)..temp(57) and minv, ioff and iscale are ignored here
c         and when reading the file.
c         floating point numbers are NOT standard across platforms and
c         is not recommended
c
c       temp(49) <= anodata            ! value representing no data
c       temp(50) <= vmin               ! minimum useful value from creator prg
c       temp(51) <= vmax               ! maximum useful value from creator prg
c       temp(52,53) = anodata          ! IEEE floating value of no data
c       temp(54,55) = vmin             ! IEEE floating minimum useful value
c       temp(56,57) = vmax             ! IEEE floating maximum useful value
c
c	temp(58)-temp(126) 138 chars of type
c       temp(127) = 0 
c       temp(128) = 0 
c	temp(129)-temp(168) 80 chars of title
c       temp(169) = 0 
c	temp(170)-temp(189) 40 chars of tag
c       temp(190) = 0 
c	temp(191)-temp(240) 100 chars of crproc
c       temp(241) = 0 
c	temp(242)-temp(255) 28 chars of crtime
c       temp(256) = 0 
c
c     optional header blocks:
c
c	ndes header blocks of 512 bytes: chars of description
c	nhead-ndes-1 header blocks of 512 bytes: values of iaopt
c       by convention, first value iaopt is a code telling how to interpret
c       the rest of the array if nia>0
c
c       remainder of file is image data in a multiple of 512 byte blocks
c       two byte integer and byte scaling is
c          intval = (fvalue-ioff)*iscale-minv
c          fvalue = float(intval+minv)/float(iscale)+ioff
c
c	open input file
c
	if (vms) then
	  open(unit=iu,name=fname,status='old',form='formatted',
     *		recl=512,err=99) !,
c     *          readonly,recordtype='fixed',carriagecontrol='none')
	else
	  open(unit=iu,name=fname,status='old',form='unformatted',
     *          recl=512,access='direct',err=99)
	endif
c
c	read first header block
c
	if (vms) then
		read (iu,414) temp
 414		format(256a2)
		call swapbuf(temp,256)
	else
		read (iu,rec=1) temp
                if (bswap) call swapbuf(temp,256)
	endif
c
c	decode first header block
c
	nsx = temp(1)			! dimension
	nsy = temp(2)
	iopt = temp(17)			! transformation option
        nhtype = temp(5)                ! header type
        if (nhtype.lt.20) nhtype=1
c
c	get default transformation decoding
c
	xdeg = 0.01 * temp(3)		! center point or span
	ydeg = 0.01 * temp(4)
	ascale = temp(6) * 0.001	! scale factor
	bscale = temp(7) * 0.001
	a0 = temp(8) * 0.01		! origin
	b0 = temp(9) * 0.01
c
c	get special cases which depend on transformation option
c
	if (iopt.eq.-1) then				! image only
		xdeg = 0.1 * temp(3)
		ydeg = 0.1 * temp(4)
	else if (iopt.eq.0) then			! rectalinear
		xdeg = 0.01 * temp(3)+100.0
	else if (iopt.eq.1.or.iopt.eq.2) then		! lambert
		ascale=1000.0/float(temp(6))
		bscale=1000.0/float(temp(7))
		a0=temp(8)
		b0=temp(9)
	else if (iopt.eq.5) then			! polar stereographic
		ascale = temp(6) * 0.01	
		bscale = temp(7) * 0.01
		a0 = temp(8)
		b0 = temp(9)
		xdeg = 0.01 * temp(3)+100.0
	else if (iopt.eq.11.or.iopt.eq.12.or.iopt.eq.13) then ! EASE grid
		ascale = 2.0d0*(0.001*float(temp(6)))*6371.228d0/
     *                  25.067525d0
		bscale = 2.0d0*(0.001*float(temp(7)))*25.067525d0
		xdeg = 0.1 * temp(3)
		ydeg = 0.1 * temp(4)
		a0 = 0.1 * temp(8)
		b0 = 0.1 * temp(9)
	else
		write (*,*)'*** Unrecognized SIR option in READSIRHEAD ***'
	endif
c
c	decode scaling, data time, region, and type information
c
	ioff = temp(10)
	iscale = temp(11)
	if (iscale.eq.0) iscale=1
	iyear = temp(12)
	isday = temp(13)
	ismin = temp(14)
	ieday = temp(15)
	iemin = temp(16)
	iregion = temp(18)
	itype = temp(19)
        nhead = temp(41)        ! total number of header blocks
        if (nhead.eq.0) nhead=1
        ndes = temp(42)         ! number of description blocks
        ldes = temp(43)
        nia = temp(44)
        ipol = temp(45)
        ifreqhm = temp(46)
        ispare1 = temp(47)
        idatatype = temp(48)    ! data type code
	if (idatatype.eq.0) idatatype=2
	soff=32767.0/float(iscale)
	if (idatatype.eq.1) soff=128.0/float(iscale)
        anodata = temp(49)/float(iscale)+ioff+soff
        vmin = temp(50)/float(iscale)+ioff+soff
        vmax = temp(51)/float(iscale)+ioff+soff
c
	if (idatatype.eq.4) then
	   i2(1)=temp(52)
	   i2(2)=temp(53)
	   anodata=f2
	   i2(1)=temp(54)
	   i2(2)=temp(55)
	   vmin=f2
	   i2(1)=temp(56)
	   i2(2)=temp(57)
	   vmax=f2
	endif
c
c	get character variables
c
	do i=1,20
		j=(i-1)*2+1
		sens1(j:j)=char(mod(temp(i+19),256))
		sens1(j+1:j+1)=char(temp(i+19)/256)
	end do
	sens1(41:41)=char(0)
	sensor=sens1
	do i=1,40
		j=(i-1)*2+1
		tit1(j:j)=char(mod(temp(i+128),256))
		tit1(j+1:j+1)=char(temp(i+128)/256)
	end do
	tit1(81:81)=char(0)
	title=tit1
        if (nhtype.gt.1) then
           do i=1,69
              j=(i-1)*2+1
              type1(j:j)=char(mod(temp(i+57),256))
              type1(j+1:j+1)=char(temp(i+57)/256)
           end do
	   type1(139:139)=char(0)
           do i=1,20
              j=(i-1)*2+1
              tag1(j:j)=char(mod(temp(i+169),256))
              tag1(j+1:j+1)=char(temp(i+169)/256)
           end do
	   tag1(41:41)=char(0)
           do i=1,50
              j=(i-1)*2+1
              crp1(j:j)=char(mod(temp(i+190),256))
              crp1(j+1:j+1)=char(temp(i+190)/256)
           end do
	   crp1(101:101)=char(0)
           do i=1,14
              j=(i-1)*2+1
              crt1(j:j)=char(mod(temp(i+241),256))
              crt1(j+1:j+1)=char(temp(i+241)/256)
           end do
	   crt1(29:29)=char(0)
	else
	   type1=' '
	   tag1=' '
	   crp1=' '
	   crt1=' '
        endif
        type=type1
        tag=tag1
        crproc=crp1
        crtime=crt1
        descrip=' '
c
c	get extra storage information
c
        if (nhead.gt.1) then
           if (ndes.gt.0) then
c
c	read next header blocks
c
              icnt=0
              do irec=2,ndes+1
                 if (vms) then
                    read (iu,414) temp
                    call swapbuf(temp,256)
                 else
                    read (iu,rec=irec) temp
                    if (bswap) call swapbuf(temp,256)
                 endif
                 do j=1,512
                    icnt=icnt+1
                    if (icnt.le.ldes.and.icnt.le.maxdes) then
                       i=(j-1)/2+1
                       iv=temp(i)
                       if (mod(j,2).eq.0) then
                          iv=iv/256
                       else
                          iv=mod(iv,256)
                       endif
                       descrip(icnt:icnt)=char(iv)
                    endif
                 end do                
              end do
              if (maxdes.lt.ldes) ldes=maxdes
           endif
           if (nhead-ndes-1.gt.0) then
              icnt=0
              do irec=ndes+2,nhead
                 if (vms) then
                    read (iu,414) temp
                    call swapbuf(temp,256)
                 else
                    read (iu,rec=irec) temp
                    if (bswap) call swapbuf(temp,256)
                 endif
                 do i=1,256
                    icnt=icnt+1
                    if (icnt.le.nia.and.icnt.le.maxi) 
     *                 iaopt(icnt) = temp(i)
                 end do
              end do
              if (nia.gt.maxi) nia=maxi
           endif
        endif
	return
c     
 300    continue
	write (*,398)
398	format(' *** ERROR reading header of input SIR file ***')
c
c	close file on header read error
c
	close(unit=iu)
c
	ierr=-2
	return

99	continue
	write (*,198)
198	format(' *** ERROR opening input SIR file ***')
	ierr=-1
	return
	end
c
c
	subroutine readsirf(iu,ierr,nhead,nhtype,idatatype,stval,nsx,nsy,
     *       ioff,iscale,smin,smax,ncnt,anodata,vmin,vmax)
c
c	read image data
c
c inputs:
c     iu       fortran logical unit of previously opened file
c     ierr     error code from header read
c     nhead    number of 512 byte blocks in header
c     nhtype   header type
c     nsx,nsy  image dimensions
c     ioff,iscale image storage scale factors
c     anodata   specify nodata value
c     vmin,vmax specify min, max valid range
c
c outputs:
c     ierr   file read error code
c             set to 0 for successful read of file
c             set to -3 for data read error
c     stval  output image array
c     smin,max min and max of image within valid range 
c              valid range = [vmin..vmax] and .ne. anodata
c     ncnt   count of pixels in valid range
c
	integer		iu,nhead,nhtype
	real		stval(*)
	integer 	nsx,nsy,ncnt
	integer		ioff,iscale
	real            anodata,vmin,vmax,smin,smax
c
c     working variables
c
	byte            tempb(512)
	integer*2	temp(256)
	real            tempf(128)
	integer         ntot,jj,irec
	real            s,soff
c
c     note: the bswap and vms flags should be set in the block data
c
	logical bswap,vms
	common /vaxorunix/bswap,vms
	save /vaxorunix/
c
	if (ierr.lt.0) return   ! return if file open previously failed
c
	ntot = nsx*nsy
	if (iscale.eq.0) iscale=1
	s=1./float(iscale)
	soff=32767.0/float(iscale)
	nrec=256
	if (idatatype.eq.1) then
	   soff=128.0/float(iscale)
	   nrec=512
	endif
	if (idatatype.eq.4) nrec=128
c
	jj=0
	irec=nhead
	smin=1.e25
	smax=-1.e25
	ncnt=0
200	continue
	        irec=irec+1
		if (vms) then
		   if (idatatype.eq.1) then
		      read (iu,413,end=299) tempb
 413		      format(512a1)
		   else if (idatatype.eq.4) then
		      read (iu,412,end=299) tempf
 412		      format(128a4)
		   else
		      read (iu,414,end=299) temp
 414		      format(256a2)
		      call swapbuf(temp,256)
		   endif
		else
		   if (idatatype.eq.1) then
		      read(iu,rec=irec,err=298) tempb
		   else if (idatatype.eq.4) then
		      read(iu,rec=irec,err=298) tempf
		   else
		      read(iu,rec=irec,err=298) temp
		      if (bswap) call swapbuf(temp,256)
		   endif
		endif
		do i=1,nrec
			jj=jj+1
			if (idatatype.eq.1) then
			   k=tempb(i)
			   stval(jj)=s*k+soff+ioff
			else if (idatatype.eq.4) then
			   stval(jj)=tempf(i)
			else
			   stval(jj)=s*temp(i)+soff+ioff
			endif
			if (stval(jj).ne.anodata.and.
     *			     stval(jj).ge.vmin.and.
     *			     stval(jj).le.vmax) then
			   smin=min(smin,stval(jj))
			   smax=max(smax,stval(jj))
			   ncnt=ncnt+1
			endif
			if (jj.ge.ntot) goto 299
		end do
		goto 200
298     continue
	write (*,98)
98	format(' *** ERROR reading input SIR file ***')
c
c	close file on read error
c
	close(unit=iu)
	ierr=-3
	return
299	continue
c
c	close file on successful completion
c
	close(unit=iu)
	ierr=0
	return
	end
c
	subroutine swapbuf(t,n)
c
c     swap byte order of 2 byte integers
c
	logical*1 t(1),s
	do i=1,n
		j=(i-1)*2+1
		s=t(j)
		t(j)=t(j+1)
		t(j+1)=s
	end do
	return
	end
c
c *************************************************************************
c
c       this block data defines the file access and byte order logicals
c       in their common block.  set the parameters according to the machine
c       code is execute.
c
	BLOCK DATA
	LOGICAL BSWAP,VMS
	COMMON /VAXORUNIX/BSWAP,VMS
C	DATA VMS/.TRUE./		! VAXVMS FILE ACCESS & BYTE ORDER
	DATA VMS/.FALSE./		! UNIX FILE ACCESS
	DATA BSWAP/.FALSE./		! HP, SUN (big-endian)
C	DATA BSWAP/.TRUE./		! VAX, PC (little-endian)
	SAVE /VAXORUNIX/
	END
C
C ***************************************************************************
